CODEWITHSUNDEEP
rgis
digitalmaps
digitalmaps

Reputation: 2905

Combining geographic layers with different projections in R

EDIT: I have reworded the title question slightly, and adjusted the text to respond to the comment by @DWin.

Combining geographic layers that are projected and not projected can be challenging. Often, it seems, some transformation is necessary, as geographic layers come from different products and publishers.

I am aware that R has several tools to perform geographic transformations. For example:

  1. For objects of class Spatial* in the sp package, the spTransform() function in the rgdal package can be used; and,
  2. For objects of class Raster* in the raster package, the projectRaster() function can be used.

Here is a specific task that I would like to accomplish in R: Transform to UTM grid Zone 15N (Datum: NAD83) a polygons layer describing lakes in a UTM grid Zone 15N (Datum: NAD27) projection (this is in an ESRI shapefile format).

Upvotes: 3

Views: 1244

Answers (1)

Spacedman
Spacedman

Reputation: 94237

The useful thing here is the epsg database included in rgdal.

epsgs = make_EPSG()
subset(epsgs,grepl("15N",epsgs$note))

[etc]
      code
2703 26715                         # NAD27 / UTM zone 15N  [etc]
2851 26915                         # NAD83 / UTM zone 15N  [etc]
[etc]

Those codes are what you need in spTransform. If your lakes are in a shapefile with that NAD27 projection, then:

require(maptools)
lakes = readShapeSpatial("lakes.shp")
proj4string(lakes)=CRS("+init=epsg:26715")

should give you the lakes as supplied (note I dont think readShapeSpatial will read a .prj file with a shapefile set, so I've set it here explicitly)

Now to convert to NAD83 datum version of UTM zone 15N:

lakes83 = spTransform(lakes,CRS("+init=epsg:26915"))

Rasters are a bit trickier since they generally involve a warp so that you end up with a regular grid in your projected coordinate system - you can't just transform the coordinates of the corners...

Upvotes: 5

Related Questions