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

[alexandre villers]

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 4731 (20091231) __________

The message was checked by ESET Mail Security.
http://www.eset.com