Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a number of shapefiles in different CRSs (mostly WGS84 lat/lon) that I'd like to transform into a common projectly (likely Albers Equal Area Conic, but I may ask for help on choosing in another question once my problem gets better-defined).

I spent a few months doing spatial stats stuff in R, but it was 5 years ago. For the life of me, I cannot remember how to transform an sp object (e.g. SpatialPolygonsDataFrame) from one projection to another.

Example code:

P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
hrr.shp <- readShapePoly("HRR_Bdry"), verbose=TRUE, proj4string=P4S.latlon) 
# Shapefile available at http://www.dartmouthatlas.org/downloads/geography/hrr_bdry.zip but you must rename all the filenames to have the same capitalization for it to work in R

Now I have a SpatialPolygonsDataFrame with appropriate projection information, but I'd like to transform it to the desired projection. I recall there being a somewhat unintuitively-named function for this, but I can't remember what it is.

Note that I do not want just to change the CRS but to change the coordinates to match ("reproject", "transform", etc.).

Thanks!

Edit

Excluding AK/HI which are annoyingly placed in Mexico for this shapfile:

library(taRifx.geo)
hrr.shp <- subset(hrr.shp, !(grepl( "AK-" , hrr.shp@data$HRRCITY ) | grepl( "HI-" , hrr.shp@data$HRRCITY )) )
proj4string(hrr.shp) <- P4S.latlon
share|improve this question
    
Previous answer on projecting using the proj4 package here. Haven't tried this with SpatialPolygonsDataFrame though. –  Simbamangu Aug 19 '12 at 7:06
    
Actually looks like proj4 doesn't work with Spatial objects - but see answer below. –  Simbamangu Aug 19 '12 at 7:50
2  
There's always the Spatial Task View: cran.r-project.org/web/views/Spatial.html and my notes on Spatial Data [shameless plug]: maths.lancs.ac.uk/~rowlings/Teaching/UseR2012 –  Spacedman Aug 19 '12 at 9:01
add comment

1 Answer 1

up vote 15 down vote accepted

You can use the spTransform() methods in rgdal - using your example, you can transform the object to NAD83 for Kansas (26978):

library(rgdal)
library(maptools)

P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
hrr.shp <- readShapePoly("HRR_Bdry", verbose=TRUE, proj4string=P4S.latlon)
plot(hrr.shp)

unprojected

hrr.shp.2 <- spTransform(hrr.shp, CRS("+init=epsg:26978"))
plot(hrr.shp.2)

projected

To save it in the new projection:

writePolyShape(hrr.shp.2, "HRR_Bdry_NAD83")

EDIT: Or, as per @Spacedman's suggestion (which writes a .prj file with the CRS info):

writeOGR(hrr.shp.2, dsn = getwd(), layer = "HRR_Bdry_NAD83", driver="ESRI Shapefile")
share|improve this answer
5  
note that writePolyShape DOES NOT write the .prj file! You should use writeOGR from rgdal (and use readOGR to read shapefiles) if you want to write and read the .prj file to set the CRS of your spatial objects in R! –  Spacedman Aug 19 '12 at 8:59
    
Much better (edited accordingly) - thanks; hadn't realised it creates the .prj file! Awesome cheatsheet on your page, by the way. –  Simbamangu Aug 19 '12 at 10:47
1  
It's strange how the projection in Mexico affects the appearance of the Alaska and Hawaii insets :-). –  whuber Aug 19 '12 at 18:03
    
@whuber - hmm, yes ... someone's edited my posting which didn't have the actual maps in it showing those rather inappropriate insets ... never plotted 'em myself to see that they were there. –  Simbamangu Aug 19 '12 at 18:30
    
@Simbamangu Sorry, forgot that this .shp file rather inappropriately included the insets when I attempted to be helpful in adding the graphs! –  Ari B. Friedman Aug 19 '12 at 21:35
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.