Projecting sp objects in R

Question

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

Answer 1

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)

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

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:

• Choosing the correct value for proj4string for shapefile reading in R?

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

• Assigning CRS to shapefile when it doesn't have one, in R?

Answer 2

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.