[R-sig-Geo] PBSmapping importShapefile question

[Julie Lee-Yaw]

Hi again,

Thanks to those who helped me with my earlier post. I've installed rgdal and have made some progress but I'm now stuck on an error message that I don't know how to interpret.

I have imported one shapefile that is projected in the CRS that I want. I then import another shapefile that is not projected using the PBSmapping function. To project it to the same CRS, I use use the spTranform function from rgdal. Plotting them together, I get a map that generally looks correct. My code is:

shape1=readShapePoly("file1", proj4string=CRS("+proj=laea +lat_0=90 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))   
shape1PS=SpatialPolygons2PolySet(shape1)
shape2<-importShapefile("file2.shp",readDBF=FALSE)
shape2Proj=spTransform(shape2,CRS("+proj=laea +lat_0=90 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
shape2PS=SpatialPolygons2PolySet(shape2Proj)
plotPolys(shape1PS)
addPolys(shape2PS)

This plots everything but returns the following warning message: 

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

My first question is whether this message is just saying the plot is a little off (e.g. one shapefile ends up outside the plot's boundaries) or whether my two shapefiles are still not in a common projection and thus will give me problems in my downstream polygon overlay analysis.

I'm also wondering if someone might be able to clarify what "projection=1" means when you use the PBSmapping function "importShapefile" and the .prj file does not specify 'GEO' nor 'PROJ' with 'UTM' elsewhere? What projection gets assigned in this case?

Thanks!

Julie

--- On Thu, 12/31/09, alexandre villers <alexandre.villers at cebc.cnrs.fr> wrote:

> From: alexandre villers <alexandre.villers at cebc.cnrs.fr>
> Subject: Re: [R-sig-Geo] stuck on importing and projecting shapefiles
> To: julleeyaw at yahoo.ca
> Cc: r-sig-geo at stat.math.ethz.ch
> Received: Thursday, December 31, 2009, 2:26 AM
> Hey Julie,
> 
> The main problem, for what I understood of working with R
> and spatial objects, is that you are using the wrong tools
> to import your data. Maptools is mostly deprecated for all
> the import export function with shapefile/rasters (but still
> holds interesting ones for other jobs).
> You'll have to struggle a bit and get rgdal to work on your
> Mac ! http://cran.r-project.org/bin/macosx/
> Once this is done, it will be much easier to do the job.
> GDAL is a very powerful, essential tool to import export
> data to/from R. So get to the help list for MAC, the forum
> and check if
> - rgdal is compatible for your version of R MAC
> - look for all the packages necessary to compile this
> packages (I couldn't download and compile rgdal with linux
> ubuntu because a small package was missing...). Just check
> through the documentation which packages are necessary to
> compile rgdal with MAC.
> 
> Once this is done, your work will be much more simple
> 
> You will find all the help you need on the rgdal/sp
> helpfile and on the Internet. But do not hesitate to ask
> questions on the list.
> 
> 
> Best regards and good luck.
> 
> Alex
> 
> PS: I'm pretty sure a R / rgdal /MAC user will help you
> with this !
> 
> -----Original Message-----
> From: Julie Lee-Yaw <julleeyaw at yahoo.ca>
> To: r-sig-geo at stat.math.ethz.ch
> Date: Wed, 30 Dec 2009 15:32:46 -0800 (PST)
> Subject: [R-sig-Geo] stuck on importing and projecting
> shapefiles
> 
> 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!
> 
> 
> 
> Alexandre Villers
> PhD. Candidate
> Team Agripop
> CEBC CNRS UPR 1934
> 79360 Beauvoir sur Niort
> 
> Phone: +33 (0)549 099 613
> 
> 
> __________ Information from ESET Mail Security, version of
> virus signature database 4730 (20091230) __________
> 
> The message was checked by ESET Mail Security.
> http://www.eset.com
> 
> 
> 


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