I have read in a shapefile of lodgepole pine seed planning units in R with readShapePoly
, like so
spu <- readShapePoly("spus.shp")
When I map the spus alone, they show up.
plot(spu, border=TRUE)
When I try to add the shapefile to a map of Canada, where I know the spus should go, they do not show up.
map("worldHires","Canada", xlim=c(-141,-110), ylim=c(45,65), col="gray90", fill=TRUE)
plot(spu, add=TRUE, border=TRUE)
I think this is because the shapefile is missing the projection information contained in the .prj file:
PROJCS["NAD_1983_Albers",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",1000000.0],
PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-126.0],PARAMETER["Standard_Parallel_1",50.0],
PARAMETER["Standard_Parallel_2",58.5],PARAMETER["Latitude_Of_Origin",45.0],UNIT["Meter",1.0]]
The partial output of
summary(spu)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 642162.8 1870556
y 457057.2 1421478
Is projected: NA
proj4string : [NA]
Data attributes:
shows that no projection information is associated with the spu object in R. I think I can fix this by using the correct proj4string
with readShapePoly,
but could not find how to translate the information in the .prj file into a proj4 string. How can I do this/is this the correct approach?