Chapter 3 Coordinate reference systems
3.1 Lesson Objectives
Describe and understand a CRS in R and components of a CRS
3.2 What makes up a CRS?
A CRS is made up of several components:
- Coordinate system: The x,y grid that defines where your data lies in space
- Horizontal and vertical units: The units describing grid along the x,y and possibly z axes
- Datum: The modeled version of the shape of the earth
- Projection details: If projected, the mathematical equation used to flatten objects from round surface (earth) to flat surface (paper or screen)
Figure 3.1: Source: https://mgimond.github.io/Spatial/coordinate-systems.html
3.3 The ellipsoid and geoid
The earth is a sphere, but more precisely, an ellipsoid, which is defined by two radii: - semi-major axis (equatorial radius) - semi-minor axis (polar radius)
The terms speroid and ellipsoid are used interchangeably. One particular spheroid is distinguished from another by the lengths of the semimajor and semiminor axes. Both GRS80 and WGS84 spheroids use semimajor and semiminor axes of 6,378,137 meters and 6,356,752 meters respectively (the semiminor axis differs by thousandths of meters between the two). You’ll encounter older spheroid / ellipsoids out there such as Clark 1866.
More precisely than an ellipsoid, though, we know that earth is a geoid - it is not perfectly smooth - and modelling the the undulations due to changes in gravitational pull for different locations is crucial to accuracy in a GIS. This is where a datum comes in - a datum is built on top of the selected spheroid and can incorporate local variations in elevation. We have many datums to choose from based on location of interest - in the US we would typically choose NAD83
3.4 Why you need to know about CRS working with spatial data in R:
library(Rspatialworkshop)
library(readr)
library(sf)
data(pnw)
gages = read_csv(system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop"),show_col_types = FALSE)
gages_sf <- gages %>%
st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE) %>%
dplyr::select(STATION_NM,LON_SITE, LAT_SITE)
# Awesome, let's plot our gage data and state boundaries!
plot(pnw$geometry, axes=TRUE)
plot(gages_sf$geometry, col='red', add=TRUE)
# um, what?There is no ‘on-the-fly’ projection in R - you need to make sure you specify the CRS of your objects, and CRS needs to match for any spatial operations or you’ll get an error
spatialreference.org is your friend in R - chances are you will use it frequently working with spatial data in R.
projection Wizard is also really useful, as is epsg.io
Useful
rgdalpackage functions:- projInfo(type=‘datum’)
- projInfo(type=‘ellps’)
- projInfo(type=‘proj’)
3.5 Changes to CRS recently in R in sf
It’s important to understand recent changes in handling of crs in sf
Prior to sf 0.9, crs were represented as lists with two components:
- epsg (European Petroleum Survey Group ) code, which could be NA
- proj4string - a projection string
library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
st_crs(4326)
# Coordinate Reference System:
# EPSG: 4326
# proj4string: "+proj=longlat +datum=WGS84 +no_defs"Now, however, sf represents crs as lists with two different components:
- input
- wkt
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
(x = st_crs(4326))
## 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]]You can still get the proj4string or epsg value by:
x$epsg
## [1] 4326
x$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"3.6 Projected coordinate systems
Typically we want to work with data that is projected. Projected coordinate systems (which are based on Cartesian coordinates) have: an origin, an x axis, a y axis, and a linear unit of measure. Going from geographic coordinates to a projected coordinate reference systems requires mathematical transformations.
Four spatial properties of projected coordinate systems that are subject to distortion are: shape, area, distance and direction. A map that preserves shape is called conformal; one that preserves area is called equal-area; one that preserves distance is called equidistant; and one that preserves direction is called azimuthal (from https://mgimond.github.io/Spatial/coordinate-systems.html.
The takeaway from all this is you need to be aware of the crs for your objects in R, make sure they are projected if appropriate and in a projection that optimizes properties you are interested in, and objects you are analyzing or mapping together need to be in same crs.
Going back to our original example, we can transform crs of objects to work with them together:
library(ggplot2)
# Check our coordinate reference systems
st_crs(gages_sf)## Coordinate Reference System:
## User input: EPSG:4269
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## 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["Geodesy."],
## AREA["North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands."],
## BBOX[14.92,167.65,86.46,-47.74]],
## ID["EPSG",4269]]st_crs(pnw)## Coordinate Reference System:
## User input: +proj=aea +lat_1=41 +lat_2=47 +lat_0=44 +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6269]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",44,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-120,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",41,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",47,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## 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]]]]# Are they equal?
st_crs(gages_sf)==st_crs(pnw)## [1] FALSE# transform one to the other
gages_sf <- st_transform(gages_sf, st_crs(pnw))
ggplot() +
geom_sf(data=gages_sf, color="blue") +
geom_sf(data=pnw, color="black", fill=NA) +
labs(title="USGS Stream Gages in the Pacific Northwest") +
theme_bw() 
3.7 Exercise
Re-project both the PNW state polygons and the stream gages to a different projected CRS - I suggest UTM zone 11, or perhaps Oregon Lambert, but choose your own. Use resources I list above to find a projection and get it’s specification to use in R to re-project both sf objects, then plot together either in base R or ggplot.
3.8 Solution
My solution:
library(cowplot)
# first I save the previous plot so we can view it alongside our update:
p1 <-ggplot() +
geom_sf(data=gages_sf, color="blue") +
geom_sf(data=pnw, color="black", fill=NA) +
labs(title="Albers Equal Area") +
theme_bw()
# You can fully specify the WKT:
utmz11 <- 'PROJCS["NAD83(CSRS98) / UTM zone 11N (deprecated)",GEOGCS["NAD83(CSRS98)",DATUM["NAD83_Canadian_Spatial_Reference_System",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6140"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4140"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","2153"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
# or you can simply provide the EPSG code:
utmz11 <- 2153
gages_sf <- st_transform(gages_sf, utmz11)
pnw <- st_transform(pnw, utmz11)
p2 <- ggplot() +
geom_sf(data=gages_sf, color="blue") +
geom_sf(data=pnw, color="black", fill=NA) +
labs(title="UTM Zone 11") +
theme_bw()
plot_grid(p1, p2)