Coordinate Reference Systems
A fundamental aspect of spatial data is defining its Coordinate Reference System (CRS). Briefly, a CRS defines how to represent the Earth’s surface on the plane. Until we define a projection, coordinates are just numbers.
A CRS is a combination of an ellipsoid, a datum, and a projection.
Ellipsoid
In reality, Earth has a irregular shape. An ellipsoid is a 3D simplified model that approximates the shape of the Earth.
Some ellipsoids are global, i.e. they try to model the entire globe. Among these, WGS84 is the most recent and more widely used global reference system. GRS80 is also widely used.
But global ellipsoids have larger errors than ellipsoids fitted locally. Local ellipsoids are used when higher precision is needed. Examples: ED50 (Europe), OSGB36 (United Kingdom).
Datum
The datum provides the information needed to “anchor” coordinates to the surface of the Earth. Defines the origin point and the direction of the coordinate axes.
The datum always specifies the ellipsoid, and sometimes have the same name.
Projection
Finally, the projection indicates how we represent locations on the surface of the earth. Basically we have two options:
- Use latitude and longitude, which are spherical coordinates that represent locations on a 3D globe defined by the datum and ellipsoid used. But remember:
- A set of latitude and longitude values can indicate different places when using different datums/ellipsoids. Differences can be important depending on the application.
- When we plot latitude/longitude values into a flat surface (like a computer screen or a piece of paper) there is an implicit transformation where we take them as Cartesian coordinates (x=longitude, y=latitude). This is the Equidistant Cylindrical projection.
- When we do not know the datum/ellipsoid, it is usually safe to assume it is WGS84. In particular for large scale mapping, when the differences among datums are not noticeable.
- Use an explicit transformation to project spherical coordinates into planar coordinates. These transformations are known as spatial projections.
Spatial projections
Different combinations of ellipsoid, datum and transformation result in many different spatial projections. They can be classified as planar, conical or cylindrical.
All map projections
introduce distortions in area, shape or distance. These properties are
mutually exclusive: a map projection that preserves one will distort the
other two.
For example, the Mercator projection is well-know for preserving shapes and exaggerating the area of objects near the poles. South America in, in reality, much larger than Greenland. In the Mercator projection, the loxodromes (or rhumb lines) are straight lines, so it is good for navigation.
On the other hands, the Azimuthal Equidistant projection does not preserve shapes (can you see Australia?). But it preserves distances to and from the center of the projection (it is an equidistant projection).
The Cylindrical Equal Area is an example of a projection that does not distort the areas (i.e is an equal area projection). Here the proportion between the areas of Greenland and South America is right… although the shapes are distorted near the poles.
Conformal projections preserve local shapes.
Which projection to use?
- Do not simply ignore CRSs. Select an appropriate projection for your map/project.
- Projections are not wrong nor right. They simply have different properties and are used for different purposes.
- For densities (points per grid cell), use equal-area projections. For example, use the Lambert azimuthal equal-area (LAEA) with center in your study area (select the lon_0 and lat_0 parameteres). LAEA maintains the area of the objects, but distorts the shape away from the center.
- Azimuthal equidistant (AEQD) projections if you need accurate straight-line distance between a point and the center point of the local projection.
- Lambert conformal conic (LCC) projections are useful for large areas (1000’s km). They preserve the shape of objects, so they are good for general mapping.
- Stereographic (STERE) projections for polar regions, but area and distances get distorted away from the centre.
- Many countries have “official” coordinate system. For example, Iceland’s ISN2004 uses a LAEA projection, centered in Iceland, and the GRS1980 ellipsoid.
- Use a projected CRS when doing geometric operations (polygon overlays, etc.), as the algorithms assume planar coordinates.
- http://spatialreference.org/ lists many, many projections
- http://projectionwizard.org/ and http://www.radicalcartography.net/index.html?projectionref may be useful to select a projection.
… or select one based on your personality: https://imgs.xkcd.com/comics/map_projections.png
CRS in R
The sf package links directly to the external library PROJ, a very complete cartographic projection library (http://trac.osgeo.org/proj/). In the PROJ library, projections are defined using projstring, a text string defining the projection, datum and associated parameters. Here are some examples:
For unprojected (geographical) data: “+proj=longlat +datum=WGS84”
For data in the UTM projection, zone 28: “+proj=utm +zone=28 +datum=WGS84”
For a Lamberts Equal Area projection centred in the Nordic Seas: “+proj=laea +lat_0=69 +lon_0=-4 +datum=WGS84”
Many common projections are included in the European Petroleum Survey Group (EPSG) database (http://www.epsg.org/), and can be specified only by using a reference number. For example, the EPSG code 4326 indicates urprojected data, so it is equivalent to the projstring “+proj=longlat +datum=WGS84”.
When we read a spatial object into, if the projection is included in the metadata then the projection will be included in the object. We can use the st_crs() function to extract it.
library(rnaturalearth)
# need to install below if not already on your computer
# install.packages("rnaturalearthhires", repos = "http://packages.ropensci.org", type = "source")
library(tidyverse)
library(sf)
iceland <- ne_countries(returnclass = "sf", country = "Iceland", scale = 10) %>%
st_geometry()
st_crs(iceland) # Find out the CRS## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]If the metadata does not includes the CRS, or if we create our own sf object, the CRS will not be defined.
mysf <- data.frame(lat = c(66,65.5), lon = c(-20, -19)) %>%
st_as_sf(coords = c("lon", "lat"))
st_crs(mysf)## Coordinate Reference System: NAWe can define it with st_crs() (or with st_set_crs() if we use it in a pipe)
mysf <- data.frame(lat = c(66,65.5), lon = c(-20, -19)) %>%
st_as_sf(coords = c("lon", "lat")) %>%
st_set_crs(4326)
st_crs(mysf)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]If I change the CRS of an object that already has it, I get a warning:
st_set_crs(mysf, 5325)## Simple feature collection with 2 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -20 ymin: 65.5 xmax: -19 ymax: 66
## Projected CRS: ISN2004 / Lambert 2004
## geometry
## 1 POINT (-20 66)
## 2 POINT (-19 65.5)By doing this we are changing the metadata of the object, but we are not reprojecting the data.
Changing the projection of vector data
Vector data is comprised of points and combination of points into lines and polygons. Reprojecting vector data consists on transforming the coordinates of these points.
In most cases we use the st_transform() function to project sf objects. This function can receive the proj4 string or the EPSG code of the target projection.
europe <- ne_countries (returnclass = "sf", continent = "Europe") %>%
st_geometry()
# Transform to the Mollweide projection
europe_moll <- europe %>% st_transform("+proj=moll")
plot(europe, axes = TRUE)
plot(europe_moll, axes = TRUE)
st_crs(europe)## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]st_crs(europe_moll)## Coordinate Reference System:
## User input: +proj=moll
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Mollweide"],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]Additional projections are available using the st_transform_proj() in the lwgeom package. The lwgeom package provides bindings to the liblwgeom library (part of the PostGis).
europe_win <- lwgeom::st_transform_proj(europe, "+proj=wintri")
plot(europe_win, axes = TRUE)
We can modify the projections by passing the required parameters. For example:
# Centred in Paris
europe_laea1 <- st_transform(europe,
crs = "+proj=laea +lon_0=-2.35 +lat_48.8")
# Centred in Moscow
europe_laea2 <- st_transform(europe,
crs = "+proj=laea +lon_0=37.6 +lat_0=55.7")
plot(europe_laea1, axes = TRUE)
plot(europe_laea2, axes = TRUE)
Projections in ggplot
One of the advantages of using the sf objects is that we can utilise the geom_sf() function when plotting with ggplot. geom_sf() uses by default the coord_sf() function to control the coordinate system of the plot. coord_sf() ensures that all layers use a common CRS (otherwise it is not possible to place them together in the same graph). You can either specify it using the CRS argument, or coord_sf() will take it from the first layer that defines a CRS.
Let’s load two sf objects with different CRS’s.:
coast <- read_sf("ftp://ftp.hafro.is/pub/data/shapes/iceland_coastline.gpkg")
st_crs(coast)## Coordinate Reference System:
## User input: unnamed
## wkt:
## PROJCRS["unnamed",
## BASEGEOGCRS["GRS 1980(IUGG, 1980)",
## DATUM["unknown",
## ELLIPSOID["GRS80",6378137,298.257222101,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of 1st standard parallel",64.25,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",65.75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Latitude of false origin",65,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-19,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Easting at false origin",500000,
## LENGTHUNIT["Meter",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",500000,
## LENGTHUNIT["Meter",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["Meter",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["Meter",1]]]depth <- read_sf("ftp://ftp.hafro.is/pub/data/shapes/iceland_contours.gpkg")
st_crs(depth)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]Note the differences:
ggplot() +
geom_sf(data = coast, fill = "darkgray") +
geom_sf(data = depth)
ggplot() +
geom_sf(data = depth) +
geom_sf(data = coast, fill = "darkgray")
ggplot() +
geom_sf(data = depth) +
geom_sf(data = coast, fill = "darkgray") +
coord_sf(crs = 4087) # World Equidistant Cylindrical projection
Here 1 degree north is the same as 1 degrees longitude.
The different Europe projections above via ggplot
p1 <-
ggplot() +
geom_sf(data = europe) +
ggtitle("longlat")
p2 <-
ggplot() +
geom_sf(data = europe_moll) +
ggtitle("moll")
p3 <-
ggplot() +
geom_sf(data = europe_laea1) +
ggtitle("Center Paris")
p4 <-
ggplot() +
geom_sf(data = europe_laea2) +
ggtitle("Center Moscow")