6.3 Which CRS to use?

The question of which CRS is tricky, and there is rarely a ‘right’ answer:“There exist no all-purpose projections, all involve distortion when far from the center of the specified frame” (Bivand, Pebesma, and Gómez-Rubio 2013).For geographic CRSs, the answer is often WGS84, not only for web mapping (covered in the previous paragraph) but also because GPS datasets and thousands of raster and vector datasets are provided in this CRS by default.WGS84 is the most common CRS in the world, so it is worth knowing its EPSG code: 4326.This ‘magic number’ can be used to convert objects with unusual projected CRSs into something that is widely understood.

What about when a projected CRS is required?In some cases, it is not something that we are free to decide:“often the choice of projection is made by a public mapping agency” (Bivand, Pebesma, and Gómez-Rubio 2013).This means that when working with local data sources, it is likely preferable to work with the CRS in which the data was provided, to ensure compatibility, even if the official CRS is not the most accurate.The example of London was easy to answer because (a) the British National Grid (with its associated EPSG code 27700) is well known and (b) the original dataset (london) already had that CRS.

In cases where an appropriate CRS is not immediately clear, the choice of CRS should depend on the properties that are most important to preserve in the subsequent maps and analysis.All CRSs are either equal-area, equidistant, conformal (with shapes remaining unchanged), or some combination of compromises of those.Custom CRSs with local parameters can be created for a region of interest and multiple CRSs can be used in projects when no single CRS suits all tasks.‘Geodesic calculations’ can provide a fall-back if no CRSs are appropriate (see proj4.org/geodesic.html).For any projected CRS the results may not be accurate when used on geometries covering hundreds of kilometers.

When deciding a custom CRS, we recommend the following:²⁷

• A Lambert azimuthal equal-area (LAEA) projection for a custom local projection (set lon_0 and lat_0 to the center of the study area), which is an equal-area projection at all locations but distorts shapes beyond thousands of kilometres.
• Azimuthal equidistant (AEQD) projections for a specifically accurate straight-line distance between a point and the centre point of the local projection.
• Lambert conformal conic (LCC) projections for regions covering thousands of kilometres, with the cone set to keep distance and area properties reasonable between the secant lines.
• Stereographic (STERE) projections for polar regions, but taking care not to rely on area and distance calculations thousands of kilometres from the center.
  A commonly used default is Universal Transverse Mercator (UTM), a set of CRSs that divides the Earth into 60 longitudinal wedges and 20 latitudinal segments.The transverse Mercator projection used by UTM CRSs is conformal but distorts areas and distances with increasing severity with distance from the center of the UTM zone.Documentation from the GIS software Manifold therefore suggests restricting the longitudinal extent of projects using UTM zones to 6 degrees from the central meridian (source: manifold.net).

Almost every place on Earth has a UTM code, such as “60H” which refers to northern New Zealand where R was invented.All UTM projections have the same datum (WGS84) and their EPSG codes run sequentially from 32601 to 32660 (for northern hemisphere locations) and 32701 to 32760 (southern hemisphere locations).

To show how the system works, let’s create a function, lonlat2UTM() to calculate the EPSG code associated with any point on the planet as follows:

lonlat2UTM = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if(lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

The following command uses this function to identify the UTM zone and associated EPSG code for Auckland and London:

epsg_utm_auk = lonlat2UTM(c(174.7, -36.9))
epsg_utm_lnd = lonlat2UTM(st_coordinates(london))
st_crs(epsg_utm_auk)$proj4string
#> [1] "+proj=utm +zone=60 +south +datum=WGS84 +units=m +no_defs"
st_crs(epsg_utm_lnd)$proj4string
#> [1] "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs"

Maps of UTM zones such as that provided by dmap.co.uk confirm that London is in UTM zone 30U.

Another approach to automatically selecting a projected CRS specific to a local dataset is to create an azimuthal equidistant (AEQD) projection for the center-point of the study area.This involves creating a custom CRS (with no EPSG code) with units of meters based on the centerpoint of a dataset.This approach should be used with caution: no other datasets will be compatible with the custom CRS created and results may not be accurate when used on extensive datasets covering hundreds of kilometers.

Although we used vector datasets to illustrate the points outlined in this section, the principles apply equally to raster datasets.The subsequent sections explain features of CRS transformation that are unique to each geographic data model, continuing with vector data in the next section (Section 6.4) and moving on to explain how raster transformation is different, in Section 6.6.