[R-sig-Geo] stuck on importing and projecting shapefiles

Julie Lee-Yaw julleeyaw at yahoo.ca
Thu Dec 31 00:32:46 CET 2009


Hi,
I am new to both R and spatial data and am very stuck on what I think should be a simple analysis. I am trying to do a polygon overlay analysis of two shapefiles from different sources. Before proceeding to my overlay functions, I am trying to plot my shapefiles to double check that they line up properly. In R, I am using PBSmapping and maptools and the following script:

## read in first shapefile using maptools so I can pull out only those polygons that are labeled "blue"

s1=readShapePoly("shapefile1")
s1blue=s1[s1$SYMB=="BLUE",]
s1PS=SpatialPolygons2PolySet(sp1blue)
shape1=combinePolys(sp1PS)
attr(shape1,"projection")<-"LL"

## read in second shapefile using PBSmapping routine (in this case, I need to use importShapefile because I have previously modified this function to import only select polygons from this file

shape2<-importShapefile("shapefile2", projection="LL")

## trying to plot both together:

plotPolys(shape1) ##produces reasonable looking plot
addPolys(shape2)

## I get the following error message

In .checkProjection(projectionPlot, projectionPoly) :
  The data's 'projection' attribute (LL) differs from the
projection of the plot region (NA).

This leads me to think that the projections are different and thus my downstream calculations of intersection areas will be off. The first shapefile is projected and has the following prj file attached:

PROJCS["Clarke_1866_Lambert_Conformal_Conic",GEOGCS["GCS_Clarke_1866",DATUM["D_Clarke_1866",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-95.0],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",77.0],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]

The second shapefile is unprojected and the metadata indicates:
projection parameters: units: decimal degrees
Datum: North America 1983

Since I'm using a MAC and can't figure out how to load rgdal, I have tried importing both shapefiles into QGIS (a separate GIS program), reprojecting them and exporting them from there and re-applying my R code but I get the same error message. Does anyone have any suggestions as to what I might be doing wrong and how I can get these shapefiles into a common projection so that I can calculate the area of their intersection? Thanks in advance for the help!



      __________________________________________________________________
Be smarter than spam. See how smart SpamGuard is at giving junk email the boot


More information about the R-sig-Geo mailing list