How to import projection in readShapePoly() function? - R

0
votes

I am importing a polygon shapefile from Arcmap that has already a projection set and comes with all its files (sbn,sbx,prj,etc).

However after I use the readShapePoly function, when I do a summary it looks like the projection information is blank. Is the projection included already or not recognized?

Object of class SpatialPolygonsDataFrame
Coordinates:
        min     max
x   35551.4 1585917
y 6318047.3 9408727
Is projected: NA 
proj4string : [NA]
Data attributes:

I know there is a proj4string attribute but it's not clear how to use it and the prj file it's already attached to the .shp file. I can consider another function that does a better job if there is one. Not sure if rgdal with readOGR does what I want.

edit follow up: I tried with readOGR thanks for the reply. I am using this code test<-readOGR(dsn=getwd(), layer="grid") and the shapefile is here speedy.sh/Ry8rU/grid.zip and it still doesn't read the projection.

1
rshapefilemap-projectionsmaptools
Yes. Use readOGR(), as it'll preserve/capture the projection info. – Josh O'Brien
@JoshO'Brien I just tried readOGR but it's not loading the projection. In the help page it states "PROJ4 string defining CRS, if default NULL, the value is read from the OGR data set" – Herman Toothrot
@JoshO'Brien I also tried with OGRSpatialRef but it doesn't pick up the .prj file. – Herman Toothrot
Hmm. Can't help you then. If there's not an error in the syntax you're using to call those functions, and if your shapefiles are correctly formed, this should work. Try, for example, map <- readOGR(dsn = system.file("vectors", package="rgdal"), layer="scot_BNG"); proj4string(map) or OGRSpatialRef(dsn = system.file("vectors", package="rgdal"), layer="scot_BNG") to see the expected behavior of those functions. The fact that it retains projection info is one big reason to prefer rgdal over maptools for data import. – Josh O'Brien

1 Answers

1
votes

If you must use RShapePoly, try something like:

to read the projection on your shapefile:

proj4string(yourshapefile.pr)

[1] "+init=epsg:2163 +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

now project it

yourshapefile.pr <- spTransform(yourshapefile, CRS( "+init=epsg:2163" )) #US National Atlas Equal Area