34

I have a number of shapefiles in different CRSs (mostly WGS84 lat/lon) that I'd like to transform into a common projection (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.).

Edit

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

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
43

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")

If one is not certain which CRS to project from, refer to the following post:

And if one wants to define/assign a CRS when data doesn't have one, refer to:

  • 10
    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
6

Since the introduction of the sf-package (have a look at the vignettes sf1, sf2, sf3, sf4 and a migration guide here) you can use st_transform() for re-reprojecting your vector data:

require(sf)

hrr_sf = st_read('HRR_Bdry.shp', stringsAsFactors = FALSE,
    crs = 4326) # has +proj=longlat +datum=WGS84
plot(hrr_sf)

hrr_sf2 = st_transform(hrr_sf, "+init=epsg:26978") # 1st option sp::CRS() not working/ needed
hrr_sf2 = st_transform(hrr_sf, 26978) # 2nd option - EPSG code as an integer
plot(hrr_sf2)

# don't think about doing this:
hrr_sf3 = st_read('HRR_Bdry.shp', stringsAsFactors = FALSE,
    crs = 26978)

# Output layer
st_write(hrr_sf2, dsn = getwd(), layer = "HRR_Bdry_NAD83", driver = "ESRI Shapefile")

sf will replace sp in the future and has due to its simplicity and speed imho several advantages compared to sp.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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