6.6 Reprojecting raster geometries

The projection concepts described in the previous section apply equally to rasters.However, there are important differences in reprojection of vectors and rasters:transforming a vector object involves changing the coordinates of every vertex but this does not apply to raster data.Rasters are composed of rectangular cells of the same size (expressed by map units, such as degrees or meters), so it is impossible to transform coordinates of pixels separately.Raster reprojection involves creating a new raster object, often with a different number of columns and rows than the original.The attributes must subsequently be re-estimated, allowing the new pixels to be ‘filled’ with appropriate values.In other words, raster reprojection can be thought of as two separate spatial operations: a vector reprojection of cell centroids to another CRS (Section 6.4), and computation of new pixel values through resampling (Section 5.3.3).Thus in most cases when both raster and vector data are used, it is better to avoid reprojecting rasters and reproject vectors instead.

The raster reprojection process is done with projectRaster() from the raster package.Like the st_transform() function demonstrated in the previous section, projectRaster() takes a geographic object (a raster dataset in this case) and a crs argument.However, projectRaster() only accepts the lengthy proj4string definitions of a CRS rather than concise EPSG codes.

It is possible to use a EPSG code in a proj4string definition with "+init=epsg:MY_NUMBER".For example, one can use the "+init=epsg:4326" definition to set CRS to WGS84 (EPSG code of 4326).The PROJ library automatically adds the rest of the parameters and converts them into "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0".

Let’s take a look at two examples of raster transformation: using categorical and continuous data.Land cover data are usually represented by categorical maps.The nlcd2011.tif file provides information for a small area in Utah, USA obtained from National Land Cover Database 2011 in the NAD83 / UTM zone 12N CRS.

cat_raster = raster(system.file("raster/nlcd2011.tif", package = "spDataLarge"))
crs(cat_raster)
#> CRS arguments:
#>  +proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
#> +no_defs

In this region, 14 land cover classes were distinguished (a full list of NLCD2011 land cover classes can be found at mrlc.gov):

unique(cat_raster)
#>  [1] 11 21 22 23 31 41 42 43 52 71 81 82 90 95

When reprojecting categorical rasters, the estimated values must be the same as those of the original.This could be done using the nearest neighbor method (ngb).This method assigns new cell values to the nearest cell center of the input raster.An example is reprojecting cat_raster to WGS84, a geographic CRS well suited for web mapping.The first step is to obtain the PROJ definition of this CRS, which can be done using the http://spatialreference.org webpage.The final step is to reproject the raster with the projectRaster() function which, in the case of categorical data, uses the nearest neighbor method (ngb):

wgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
cat_raster_wgs84 = projectRaster(cat_raster, crs = wgs84, method = "ngb")

Many properties of the new object differ from the previous one, including the number of columns and rows (and therefore number of cells), resolution (transformed from meters into degrees), and extent, as illustrated in Table 6.1 (note that the number of categories increases from 14 to 15 because of the addition of NA values, not because a new category has been created — the land cover classes are preserved).

Reprojecting numeric rasters (with numeric or in this case integer values) follows an almost identical procedure.This is demonstrated below with srtm.tif in spDataLarge from the Shuttle Radar Topography Mission (SRTM), which represents height in meters above sea level (elevation) with the WGS84 CRS:

con_raster = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
crs(con_raster)
#> CRS arguments:
#>  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

We will reproject this dataset into a projected CRS, but not with the nearest neighbor method which is appropriate for categorical data.Instead, we will use the bilinear method which computes the output cell value based on the four nearest cells in the original raster.The values in the projected dataset are the distance-weighted average of the values from these four cells:the closer the input cell is to the center of the output cell, the greater its weight.The following commands create a text string representing the Oblique Lambert azimuthal equal-area projection, and reproject the raster into this CRS, using the bilinear method:

equalarea = "+proj=laea +lat_0=37.32 +lon_0=-113.04"
con_raster_ea = projectRaster(con_raster, crs = equalarea, method = "bilinear")
crs(con_raster_ea)
#> CRS arguments:
#>  +proj=laea +lat_0=37.32 +lon_0=-113.04 +ellps=WGS84

Raster reprojection on numeric variables also leads to small changes to values and spatial properties, such as the number of cells, resolution, and extent.These changes are demonstrated in Table 6.2³⁰:

#> Warning: `data_frame()` is deprecated, use `tibble()`.#> This warning is displayed once per session.

Of course, the limitations of 2D Earth projections apply as much to vector as to raster data.At best we can comply with two out of three spatial properties (distance, area, direction).Therefore, the task at hand determines which projection to choose.For instance, if we are interested in a density (points per grid cell or inhabitants per grid cell) we should use an equal-area projection (see also Chapter 13).

There is more to learn about CRSs.An excellent resource in this area, also implemented in R, is the website R Spatial.Chapter 6 for this free online book is recommended reading — see: rspatial.org/spatial/rst/6-crs.html