In this webinar, we will see how to create all sorts of GIS maps with the packages sf, tmap, raster, leaflet, ggplot2, grid (part of Base R), and mapview:
Examples: temperature, air quality, elevation, water depth.
Vector data
Vector data come in multiple types:
point: single set of coordinates,
multi-point: multiple sets of coordinates,
polyline: multiple sets for which the order matters,
multi-polyline: multiple of the above,
polygon: same as polyline but first and last sets are the same,
multi-polygon: multiple of the above.
Raster data
Grid of equally sized rectangular cells containing values for some variables.
Size of cells = resolution.
For computing efficiency, rasters do not have coordinates of each cell, but the bounding box and the number of rows and columns.
Coordinate Reference Systems (CRS)
A location on Earth’s surface can be identified by its coordinates and some reference system called CRS.
The coordinates (x, y) are called longitude and latitude.
There can be a 3rd coordinate (z) for elevation or other measurement—usually a vertical one.
And a 4th (m) for some other data attribute—usually a horizontal measurement.
In 3D, longitude and latitude are expressed in angular units (e.g. degrees) and the reference system needed is an angular CRS or geographic coordinate system (GCS).
In 2D, they are expressed in linear units (e.g. meters) and the reference system needed is a planar CRS or projected coordinate system (PCS).
Datums
Since the Earth is not a perfect sphere, we use spheroidal models to represent its surface. Those are called geodetic datums.
Some datums are global, others local (more accurate in a particular area of the globe, but only useful there).
Examples of commonly used global datums:
WGS84 (World Geodesic System 1984),
NAD83 (North American Datum of 1983).
Angular CRS
An angular CRS contains a datum, an angular unit and references such as a prime meridian (e.g. the Royal Observatory, Greenwich, England).
In an angular CRS or GCS:
Longitude (\(\lambda\)) represents the angle between the prime meridian and the meridian that passes through that location.
Latitude (\(\phi\)) represents the angle between the line that passes through the center of the Earth and that location and its projection on the equatorial plane.
Longitude and latitude are thus angular coordinates.
Projections
To create a two-dimensional map, you need to project this 3D angular CRS into a 2D one.
Various projections offer different characteristics. For instance:
some respect areas (equal-area),
some respect the shape of geographic features (conformal),
some almost respect both for small areas.
It is important to choose one with sensible properties for your goals.
Examples of projections:
Mercator,
UTM,
Robinson.
Planar CRS
A planar CRS is defined by a datum, a projection and a set of parameters such as a linear unit and the origins.
Common planar CRS have been assigned a unique ID called EPSG code which is much more convenient to use.
In a planar CRS, coordinates will not be in degrees anymore but in meters (or other length unit).
Projecting into a new CRS
You can change the projection of your data.
Vector data won’t suffer any loss of precision, but raster data will.
→ best to try to avoid reprojecting rasters: if you want to combine various datasets which have different projections, reproject vector data instead.
1 Bivand, R.S. Progress in the R ecosystem for representing and handling spatial data. J Geogr Syst (2020). https://doi.org/10.1007/s10109-020-00336-0.
Well-known text (WKT) is a markup language for representing vector geometry objects according to those standards.
A compact computer version also exists—well-known binary (WKB)—used by spatial databases.
The package sp predates Simple Features.
sf—launched in 2016—implements these standards in R in the form of sf objects: data.frames (or tibbles) containing the attributes, extended by sfc objects or simple feature geometries list-columns.
2 RGI Consortium (2017). Randolph Glacier Inventory – A Dataset of Global Glacier Outlines: Version 6.0: Technical Report, Global Land Ice Measurements from Space, Colorado, USA. Digital Media. DOI: https://doi.org/10.7265/N5-RGI-60.
3 Fagre, D.B., McKeon, L.A., Dick, K.A. and Fountain, A.G., 2017, Glacier margin time series (1966, 1998, 2005, 2015) of the named glaciers of Glacier National Park, MT, USA: U.S. Geological Survey data release. DOI: https://doi.org/10.5066/F7P26WB1.
4 Farinotti, Daniel, 2019, A consensus estimate for the ice thickness distribution of all glaciers on Earth - dataset, Zurich. ETH Zurich. DOI: https://doi.org/10.3929/ethz-b-000315707.
The datasets can be downloaded as zip files from these websites
Packages
Packages need to be installed before they can be loaded in a session.
Packages on CRAN can be installed with:
install.packages("<package-name>")
basemaps is not on CRAN and needs to be installed from GitHub thanks to devtools:
We load all the packages that we will need at the top of the script:
library(sf) # spatial vector data manipulation
Linking to GEOS 3.12.0, GDAL 3.7.2, PROJ 9.2.1; sf_use_s2() is TRUE
library(tmap) # map production and tiled web map
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
library(dplyr) # non GIS specific (tabular data manipulation)library(magrittr) # non GIS specific (pipes)library(purrr) # non GIS specific (functional programming)library(rnaturalearth) # basemap data access functions
Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(rnaturalearthdata) # basemap data
Attaching package: 'rnaturalearthdata'
The following object is masked from 'package:rnaturalearth':
countries110
library(mapview) # tiled web maplibrary(grid) # (part of base R) used to create inset maplibrary(ggplot2) # alternative to tmap for map productionlibrary(ggspatial) # spatial framework for ggplot2library(terra) # gridded spatial data manipulation
terra 1.7.46
Attaching package: 'terra'
The following object is masked from 'package:grid':
depth
The following objects are masked from 'package:magrittr':
extract, inset
library(ggmap) # download basemap data
Error: package or namespace load failed for 'ggmap' in dyn.load(file, DLLpath = DLLpath, ...):
unable to load shared object '/home/marie/R/lib/stringi/libs/stringi.so':
libicui18n.so.72: cannot open shared object file: No such file or directory
library(basemaps) # download basemap datalibrary(magick) # wrapper around ImageMagick STL
This dataset contains the contour of all glaciers on Earth.
We will focus on glaciers in Western North America.
You can download and unzip 02_rgi60_WesternCanadaUS and 01_rgi60_Alaska from the Randolph Glacier Inventory version 6.0.
Reading in data
Data get imported and turned into sf objects with the function sf::st_read:
ak <-st_read("data/01_rgi60_Alaska")
Make sure to use the absolute paths or the paths relative to your working directory (which can be obtained with getwd).
ak <-st_read("data/01_rgi60_Alaska")
[Out]
Reading layer `01_rgi60_Alaska' from data source `./data/01_rgi60_Alaska' using driver `ESRI Shapefile'Simple feature collection with 27108 features and 22 fieldsGeometry type: POLYGONDimension: XYBounding box: xmin:-176.1425 ymin:52.05727 xmax:-126.8545 ymax:69.35167Geodetic CRS: WGS 84
Your turn:
Read in the data for the rest of north western America (from 02_rgi60_WesternCanadaUS) and create an sf object called wes.
Classes ‘sf’ and 'data.frame':27108 obs. of 23 variables:$ RGIId : chr "RGI60-01.00001""RGI60-01.00002""RGI60-01.00003" ...$ GLIMSId : chr "G213177E63689N""G213332E63404N""G213920E63376N" ...$ BgnDate : chr "20090703""20090703""20090703""20090703" ...$ EndDate : chr "-9999999""-9999999""-9999999""-9999999" ...$ CenLon : num -147-147-146-146-147 ...$ CenLat : num 63.763.463.463.463.6 ...$ O1Region: chr "1""1""1""1" ...$ O2Region: chr "2""2""2""2" ...$ Area : num 0.360.5581.6853.6812.573 ...$ Zmin : int 1936171316091273149412011918282617501280 ...$ Zmax : int 2725214421822317231735472811355525141998 ...$ Zmed : int 2385200518681944191417402194319519771666 ...$ Slope : num 42161819162223451817 ...$ Aspect : int 346162175195181331518027435 ...$ Lmax : int 839119721064175298110518181861322553332 ...$ Status : int 0000000000 ...$ Connect : int 0000000000 ...$ Form : int 0000000000 ...$ TermType: int 0000000000 ...$ Surging : int 9999999999 ...$ Linkages: int 9999999999 ...$ Name : chr NANANANA ...$ geometry:sfc_POLYGON of length 27108; first list element: List of 1..$: num [1:65, 1:2] -147-147-147-147-147 .....-attr(*, "class")= chr [1:3] "XY""POLYGON""sfg"-attr(*, "sf_column")= chr "geometry"-attr(*, "agr")= Factor w/3 levels "constant","aggregate",..:NANANA .....-attr(*, "names")= chr [1:22] "RGIId""GLIMSId""BgnDate""EndDate" ...
Inspect your data
Your turn:
Inspect the wes object you created.
Glacier National Park dataset
This dataset contains a time series of the retreat of 39 glaciers of Glacier National Park, MT, USA for the years 1966, 1998, 2005 and 2015.
You can download and unzip the 4 sets of files from the USGS website.
Read in and clean datasets
Create a function that reads and cleans the data:
prep <-function(dir) { g <-st_read(dir) g %<>%rename_with(~tolower(gsub("Area....", "area", .x))) g %<>% dplyr::select( year, objectid, glacname, area, shape_leng, x_coord, y_coord, source_sca, source )}
Create a vector of dataset names:
dirs <-grep("data/GNPglaciers_.*", list.dirs(), value = T)
Pass each element of that vector through prep() thanks to map():
gnp <-map(dirs, prep)
We use dplyr::select because terra also has a select function.
Combine datasets into one sf object
Check that the CRS are all the same:
all(sapply(list(st_crs(gnp[[1]]),st_crs(gnp[[2]]),st_crs(gnp[[3]]),st_crs(gnp[[4]])),function(x) x ==st_crs(gnp[[1]])))
[Out]
[1] TRUE
We can rbind the elements of our list:
gnp <-do.call("rbind", gnp)
You can inspect your new sf object by calling it or with str.
Estimate for ice thickness
This dataset contains an estimate for the ice thickness of all glaciers on Earth.
The nomenclature follows the Randolph Glacier Inventory.
Ice thickness being a spatial field, this is raster data.
We will use data in RGI60-02.16664_thickness.tif from the ETH Zürich Research Collection which corresponds to one of the glaciers (Agassiz) of Glacier National Park.
We need to reproject ag (remember that it is best to avoid reprojecting raster data):
ag %<>%st_transform(st_crs(ras))
The retreat and ice thickness layers will hide each other (the order matters!). One option is to use tm_borders for one of them, but we can also use transparency (alpha). We also adjust the legend:
