6 Spatial Data and Maps

In this chapter, we’ll explore the basics of

  • building and using feature-based (vector) data using the sf (simple features) package
  • raster data and methods, using the terra and raster packages
  • some common mapping methods in ggplot2, tmap, mapview, and leaflet

Especially for the feature data which is built on a database model (with geometry as one variable), the tools we’ve been learning in dplyr and tidyr will also be useful to working with attributes.

Caveat: Spatial data and analysis is very useful for environmental data analysis, and we’ll only scratch the surface. Readers are encouraged to explore other resources on this important topic, such as Geocomputation with R at https://geocompr.robinlovelace.net/ or https://r-spatial.github.io/sf/ .

6.1 Spatial Data

To work with spatial data requires extending R to deal with it using packages. Many have been developed, but the field is starting to mature using international open GIS standards.

sp (until recently, the dominant library of spatial tools)

  • Includes functions for working with spatial data
  • Includes spplot to create maps
  • Also needs rgdal package for readOGR – reads spatial data frames
  • Earliest archived version on CRAN is 0.7-3, 2005-04-28, with 1.0 released 2012-09-30. Authors: Edzer Pebesma & Roger Bivand

sf (Simple Features)

  • Feature-based (vector) methods
  • ISO 19125 standard for GIS geometries
  • Also has functions for working with spatial data, but clearer to use.
  • Doesn’t need many additional packages, though you may still need rgdal installed for some tools you want to use.
  • Replacing sp and spplot though you’ll still find them in code. We’ll give it a try…
  • Works with ggplot2 and tmap for nice looking maps.
  • Cheat sheet: https://github.com/rstudio/cheatsheets/raw/master/sf.pdf
  • Earliest CRAN archive is 0.2 on 2016-10-26, 1.0 not until 2021-06-29.
  • Author: Edzer Pebesma


  • The established raster package for R, with version 1.0 from 2010-03-20, Author: Robert J. Hijmans (Professor of Environmental Science and Policy at UC Davis)
  • Depends on sp
  • CRAN: “Reading, writing, manipulating, analyzing and modeling of spatial data. The package implements basic and high-level functions for raster data and for vector data operations such as intersections. See the manual and tutorials on https://rspatial.org/ to get started.”


  • Intended to eventually replace raster, which it imports, by the same author (Hijmans), with version 0.5-2 2020-03-20, and 1.0 2021-01-04. Does not depend on sp, but includes Bivand and Pebesma as contributors.
  • CRAN: “Methods for spatial data analysis with raster and vector data. Raster methods allow for low-level data manipulation as well as high-level global, local, zonal, and focal computation. The predict and interpolate methods facilitate the use of regression type (interpolation, machine learning) models for spatial prediction, including with satellite remote sensing data. Processing of very large files is supported. See the manual and tutorials on https://rspatial.org/terra/ to get started. ‘terra’ is very similar to the ‘raster’ package; but ‘terra’ can do more, is easier to use, and it is faster.”

6.1.1 Simple geometry building in sf

Simple Features (sf) includes a set of simple geometry building tools that use matrices or lists of matrices of XYZ (or just XY) coordinates (or lists of lists for multipolygon) as input, and each basic geometry type having single and multiple versions:

Table 6.1: Some common sf functions for building geometries from coordinates
sf function input
st_point() numeric vector of length 2, 3, or 4, represent a single point
st_multipoint() numeric matrix with points in rows
st_linestring() numeric matrix with points in rows
st_multilinestring() list of numeric matrices with points in rows
st_polygon() list of numeric matrices with points in rows
st_multipolygon() list of lists with numeric matrices

sf functions have the pattern st_* (st means “space and time”)

See Geocomputation with R at https://geocompr.robinlovelace.net/ or https://r-spatial.github.io/sf/ for more details, but here’s an example of manual feature creation of sf geometries (sfg) followed by building a simple feature geometry list column from the collection with st_sfc():

library(sf); library(tidyverse)
eyes <- st_multipoint(rbind(c(1,5), c(3,5)))
nose <- st_point(c(2,4))
mouth <- st_linestring(rbind(c(1,3),c(3, 3)))
border <- st_polygon(list(rbind(c(0,5), c(1,2), c(2,1), c(3,2), 
                              c(4,5), c(3,7), c(1,7), c(0,5))))
face <- st_sfc(eyes, nose, mouth, border)  # sfc = sf column 
Building simple geometries in sf

Figure 6.1: Building simple geometries in sf

The face was a simple feature column (sfc) built from the list of sfgs. An sfc just has the one column, so is not quite like a shapefile.

  • But it can have a coordinate referencing system CRS, and so can be mapped.
  • Kind of like a shapefile with no other attributes than shape

6.1.2 Building a mappable sfc from scratch

The face-building example above is really intended for helping to understand the concept and isn’t really how we’d typically build data. For one thing, it doesn’t have a real coordinate system that can be reprojected into locations anywhere on the planet. And most commonly we’ll be accessing data that have been created by government agencies such as those from the USGS in The National Map, or NGO, university or corporate sources who have an interest in collecting and providing data.

However, let’s take the sfc building a bit further to better understand it, and build a mappable sfc from scratch. We’ll build a generalized map of California and Nevada boundaries (and then later in the exercises, you’ll expand this to a map of the west.) We’ll use a shortcut to the coordinate referencing system (crs) by using a short numeric EPSG code (EPSG Geodetic Parameter Dataset (n.d.)), 4326, for a common geographic coordinate system (GCS) of latitude and longitude.

CA_matrix <- rbind(c(-124,42),c(-120,42),c(-120,39),c(-114.5,35),
NV_matrix <- rbind(c(-120,42),c(-114,42),c(-114,36),c(-114.5,36),
CA_list <- list(CA_matrix);       NV_list <- list(NV_matrix)
CA_poly <- st_polygon(CA_list);   NV_poly <- st_polygon(NV_list)
sfc_2states <- st_sfc(CA_poly,NV_poly,crs=4326)  # crs=4326 specifies GCS

Then we’ll use the geom_sf() function in ggplot2 to create a map:

ggplot() + geom_sf(data = sfc_2states)
A simple map built from scratch with hard-coded data as simple feature columns

Figure 6.2: A simple map built from scratch with hard-coded data as simple feature columns

sf class

So far, we’ve built an sfg (sf geometry) and then a collection of these in an sfc (sf column), and we were even able to make a map of it. But to be truly useful, we’ll want other attributes other than just where the feature is, so we’ll build an sf that has attributes, and that can also be used like a data frame – similar to a shapefile, since a data frame is like a database table .dbf. We’ll build one with the st_sf() function from attributes and geometry.

attributes <- bind_rows(c(abb="CA", area=423970, pop=39.56e6),
                        c(abb="NV", area=286382, pop=3.03e6))
twostates <- st_sf(attributes, geometry = sfc_2states)
ggplot(twostates) + geom_sf() + geom_sf_text(aes(label = abb))
Using an sf class to build a map, displaying an attribute

Figure 6.3: Using an sf class to build a map, displaying an attribute

Or we could let the plot system try to figure out what to do with it. To help it out a bit, we can select just one variable, for which we need to use the [] accessor. Try the $ accessor instead to see what you get, and consider why this might be.


6.1.3 Creating features from shapefiles

sf’s st_read reads shapefiles

  • shapefile is an open GIS format for points, polylines, polygons

You would normally have shapefiles (and all the files that go with them – .shx, etc.) stored on your computer, but we’ll access one from the iGIScData external data folder, and use that ex() function we used earlier with CSVs. Remember that we could just include that and the library calls just once at the top of our code like this…

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}

… then here’s the code that we’re talking about now:

CA_counties <- st_read(ex("CA_counties.shp"))
## Reading layer `CA_counties' from data source `C:\Users\900008452\Documents\R\win-library\4.1\iGIScData\extdata\CA_counties.shp' using driver `ESRI Shapefile'
## Simple feature collection with 58 features and 60 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.4152 ymin: 32.53427 xmax: -114.1312 ymax: 42.00952
## Geodetic CRS:  WGS 84

And with just one variable, and for the subset of Bay Area Counties:


There’s a lot more we could do with the base plot system, but we’ll mostly focus on some better options in ggplot2 and tmap.

6.1.4 Creating sf data from data frames

If your data frame has geospatial coordinates like LONGITUDE and LATITUDE …

## # A tibble: 62 x 9
##    <chr>                             <fct>         <dbl>    <dbl>     <dbl>         <dbl>       <dbl>  <dbl>   <dbl>
##  1 GROVELAND 2, CA US                Tuolumne      853.      37.8     -120.         176.          6.1 -0.574   6.67 
##  2 CANYON DAM, CA US                 Plumas       1390.      40.2     -121.         164.          1.4 -2.00    3.40 
##  3 KERN RIVER PH 3, CA US            Kern          824.      35.8     -118.          67.1         8.9  2.05    6.85 
##  4 DONNER MEMORIAL ST PARK, CA US    Nevada       1810.      39.3     -120.         167.         -0.9 -1.74    0.840
##  5 BOWMAN DAM, CA US                 Nevada       1641.      39.5     -121.         277.          2.9  1.03    1.87 
##  6 GRANT GROVE, CA US                Tulare       2012.      36.7     -119.         186.          1.7  2.09   -0.394
##  7 LEE VINING, CA US                 Mono         2072.      38.0     -119.          71.9         0.4  1.16   -0.760
##  8 OROVILLE MUNICIPAL AIRPORT, CA US Butte          57.9     39.5     -122.         138.         10.3 -1.23   11.5  
##  9 LEMON COVE, CA US                 Tulare        156.      36.4     -119.          62.7        11.3  0.373  10.9  
## 10 CALAVERAS BIG TREES, CA US        Calaveras    1431       38.3     -120.         254           2.7 -0.450   3.15 
## # ... with 52 more rows

… we have what we need to create spatial data from it. We’ll use st_as_sf() for that, but we’ll need to specify the coordinate referencing system (CRS), in this case GCS.We’ll only briefly explore how to specify the CRS here. For a thorough coverage, please see Robin Lovelace (n.d.) .

6.1.5 Coordinate Referencing System

There are many different coordinate referencing systems (CRS), some spherical like the geographic coordinate system (GCS) of longitude and latitude, others planar like UTM or state plane. Even for GCS, there are many varieties since geodetic datums can be chosen, and for very fine resolution work where centimetres or even millimetres matter, this decision can be important (and tectonic plate movement can play havoc with tying it down.)

There are also multiple ways of expressing the CRS, either to read it or to set it. The full specification of a CRS can be displayed for data already in a CRS:

## 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["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

… though the terra::crs tool displays it with the PROJ-string method, which in this case is considerably shorter:

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

An even shorter way to represent a CRS is with the EPSG numeric code. For instance, to convert the sierra data into geospatial data with st_as_sf, we might either do it with the reasonably short PROJ format …

GCS <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
wsta = st_as_sf(sierraFeb, coords = c("LONGITUDE","LATITUDE"), crs=GCS)

… or with the even shorter EPSG code (which we can find by Googling), which can either be provided as text “epsg:4326” or often just the number:

wsta <- st_as_sf(sierraFeb, coords = c("LONGITUDE","LATITUDE"), crs=4326)

6.1.6 Removing Geometry

There are many instances where you want to remove geometry from a sf data frame

  • Some R functions run into problems with geometry and produce confusing error messages, like “non-numeric argument”

  • You’re wanting to work with an sf data frame in a non-spatial way

One way to remove geometry:

myNonSFdf <- mySFdf %>% st_set_geometry(NULL)

6.2 Base plot extended a bit with terra

As we’ve seen in the examples above, the venerable plot system in base R can often do a reasonable job of displaying a map. The terra package extends this a bit by providing tools for overlaying features on other features (or rasters using plotRGB).

lines(vect(CA_counties), col="gray")

Alternatively, we could use plot’s add=TRUE parameter to add onto a plot that we’ve previously scaled with a plot that we’ll cover over:

plot(vect(wsta)) # simple way of establishing the coordinates, will be covered by the following:
plot(vect(CA_counties), col="lightgray", add=TRUE)
plot(vect(wsta), col="red", add=TRUE)

6.2.1 Using maptiles to create a basemap

For vector data, it’s often nice to display over a basemap by accessing raster tiles that are served on the internet by various providers. We’ll use the maptiles package to try displaying the CA and NV boundaries we created earlier. The maptiles package supports a variety of basemap providers, and I’ve gotten the following to work: “OpenStreetMap,” “Stamen.Terrain,” “Stamen.TerrainBackground,” “Esri.WorldShadedRelief,” “Esri.NatGeoWorldMap,” “Esri.WorldGrayCanvas,” “CartoDB.Positron,” “CartoDB.Voyager,” “CartoDB.DarkMatter,” “OpenTopoMap,” “Wikimedia,” however they don’t work at all scales and locations – you’ll often see an Error in grDevices, if so then try another provider – the default “OpenStreetMap” seems to work the most reliably.

We’ll also use the terra::vect() function to create a SpatVector that works with terra::lines to display the boundary lines in plot(). (More on terra later, in the raster section)

library(terra); library(maptiles)
calnevaBase <- get_tiles(twostates, provider="OpenTopoMap",)  # gets the raster that covers the extent of CANV
plotRGB(calnevaBase)  # starts plot with a raster basemap
lines(vect(twostates), col="black", lwd=2)

And for the sierraFeb data, we’ll start with st_as_sf and the coordinate system (4326 for GCS), then use a maptiles basemap again, and the terra::points method to add the points:

library(terra); library(maptiles)
sierraFebpts <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs=4326)
sierraBase <- get_tiles(sierraFebpts)

6.3 Geospatial abstraction and analysis tools

See Robin Lovelace (n.d.) and Hijmans (n.d.) for more information on geospatial methods, but we’ll look at a selection of useful geospatial abstraction and analysis methods.

But before we get into those specifically spatial operations, it’s important to remember that feature data at least are also data frames, so we can use the methods we already know on them. For instance, we can look at properties of variables and then select features that meet a criterion, like all climate station points at greater than 2000 m elevation.

library(terra); library(maptiles); library(dplyr); library(stringr)
sierraFebpts <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs=4326)
names(sierraFebpts) # to quickly see what its variables are
## [1] "STATION_NAME"  "COUNTY"        "ELEVATION"     "PRECIPITATION" "TEMPERATURE"   "resid"         "predict"       "geometry"

sierraFebpts2000 <- sierraFebpts %>%
  filter(ELEVATION >= 2000)
sierraBase <- get_tiles(sierraFebpts2000)
pts <- vect(sierraFebpts2000)
points(pts, col="red")
text(pts, labels=str_sub(pts$STATION_NAME, 1, str_locate(pts$STATION_NAME, ",")-1), cex = 0.7, pos=3)

The iGIScData::BayAreaTracts data has census data from 2010 by census tracts in the SF Bay area. The code below includes a few useful sf functions applied to that data:

  • st_centroid() : to derive a single point for each tract that should be approximately in its middle
  • st_make_valid() : to make an invalid geometry valid (or just check for it). This function that has just become essential now that sf supports spherical instead of just planar data, which is the case BayAreaTracts which ends up containing “slivers” where boundaries slightly under- or over-lap since they were originally built from a planar projection. The st_make_valid() appears to make minor (and typically invisible) corrections to these slivers.
  • st_bbox : reads the bounding box, or spatial extent, of any dataset, which we can then use to set the scale to focus on that area. In this case, we’ll focus on the census centroids instead of the statewide CAfreeways for instance.
BayAreaTracts <- st_make_valid(BayAreaTracts)
censusCentroids <- st_centroid(BayAreaTracts)
TRI_sp <- st_as_sf(TRI_2017_CA, coords = c("LONGITUDE", "LATITUDE"), 
       crs=4326) # simple way to specify coordinate reference
bnd <- st_bbox(censusCentroids)
ggplot() +
 geom_sf(data = BayAreaCounties, aes(fill = NAME)) +
 geom_sf(data = censusCentroids) +
 geom_sf(data = CAfreeways, color = "grey", alpha=0.5) +
 geom_sf(data = TRI_sp, color = "yellow") +
 coord_sf(xlim = c(bnd[1], bnd[3]), ylim = c(bnd[2], bnd[4])) +
   labs(title="Bay Area Counties, Freeways and Census Tract Centroids")
ggplot map of Bay Area TRI sites, census centroids, freeways

Figure 6.4: ggplot map of Bay Area TRI sites, census centroids, freeways

6.3.1 Creating a buffer

Creating buffers, or polygons defining the area within some distance of a feature, is commonly used in GIS. Since you need to specify that distance (or read it from a variable for each feature), you need to know what the horizontal distance units are for your data. If GCS, these will be decimal degrees and 1 degree is a long way, about 111 km (or 69 miles), though that differs for longitude where 1 degree of longitude is 111 km * cos(latitude). If in UTM, the horizontal distance will be in metres, but in the US, state plane coordinates are typically in feet. So let’s read the trails shapefile from the Marble Mountains and look for its units:

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
trails <- st_read(ex("trails.shp"))
## Reading layer `trails' from data source `C:\Users\900008452\Documents\R\win-library\4.1\iGIScData\extdata\trails.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 8 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 481903.8 ymin: 4599196 xmax: 486901.9 ymax: 4603200
## Projected CRS: NAD83 / UTM zone 10N

Then we know we’re in metres, so we’ll create a 100-m buffer this way:

trailbuff <- st_buffer(trails,100)
plot(trailbuff[0]) # for just the boundaries

6.3.2 Spatial join st_join

A spatial join can do many things, but we’ll just look at its simplest application – connecting a point with the polygon that it’s in. In this case, we want to join the attribute data in BayAreaTracts with EPA Toxic Release Inventory (TRI) point data (at factories and other point sources) and then display a few ethnicity variables from the census, associated spatially with the TRI point locations. We’ll be making better maps later; this is a quick plot() display to show that the spatial join worked. [air_quality]

TRI_sp <- st_as_sf(TRI_2017_CA, coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
  st_join(st_make_valid(BayAreaTracts)) %>%
  filter(!is.na(TRACT)) %>%  # removes points not in the Bay Area -- not in BayAreaTracts

6.4 Raster

Simple Features are feature-based, of course, so it’s not surprising that sf doesn’t have support for rasters. So we’ll want to use the raster package or its imminent replacement, terra:

6.4.1 Building rasters

We can start by building one from scratch. The following raster::raster() function…

new_ras <- raster(nrows = 6, ncols = 12, 
                  xmn = -180, xmx = 180, ymn = -90, ymx = 90,
                  vals = 1:72)
## class      : RasterLayer 
## dimensions : 6, 12, 72  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 72  (min, max)

… can also be done in the newer terra package with the terra::rast() function:

new_ras <- rast(ncol = 12, nrow = 6,
                xmin = -180, xmax = 180, ymin = -90, ymax = 90,
                vals = 1:72)
## class       : SpatRaster 
## dimensions  : 6, 12, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source      : memory 
## name        : lyr.1 
## min value   :     1 
## max value   :    72

Look at the coordinate reference. How does it know that the data are in GCS (or “+proj=longlat +datum=WGS84 +no_defs”)? The author of this tool let it assume GCS if it doesn’t exceed the limits of the planet: -180 to +180 longitude and -90 to +90 latitude.

See https://rspatial.org/terra/spatial/4-rasterdata.html for an example that builds a map of the world using similar code.

The terra package builds geometries called SpatVectors which appear to work like simple features. Read more about them at rspatial.org. In any case, simple feature data or SpatVectors can be plotted along with rasters using terra’s lines, polys, or points functions. Here we’ll use the twostates sf we created earlier (make sure to run that code again if you need it). Look closely at the map for our two states:

# add polygon and points
CANV <- vect(twostates) #, type="polygons", crs="+proj=longlat +datum=WGS84")

6.4.2 Map algebra in terra

Map algebra was originally developed by Dana Tomlin in the 1970s and 1980s (Tomlin (1990)), and was the basis for his Map Analysis Package. It works by assigning raster outputs from an algebraic expression of raster inputs. Map algebra was later incorporated in Esri’s Grid and Spatial Analyst subsystems of Arc/Info and ArcGIS. Its simple and elegant syntax makes it still one of the best ways to manipulate raster data.

We’ll start by reading in some elevation data from the Marble Mountains of California and use terra’s terrain function to derive slope and aspect rasters…

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
elev <- rast(ex("elev.tif"))

slope <- terrain(elev, v="slope")

aspect <- terrain(elev, v="aspect")

… then we’ll use terra::classify to make six discrete categorical slope classes (though the legend suggests it’s continuous)…

slopeclasses <-matrix(c(00,10,1, 10,20,2, 20,30,3,
                        30,40,4, 40,50,5, 50,90,6), ncol=3, byrow=TRUE)
slopeclass <- classify(slope, rcl = slopeclasses)

… then finally a couple of simple map algebra statements to derive some new rasters …

elevft <- elev / 0.3048

… including some that create and use Boolean (true-false) values where 1 is true and 0 is false, so answers the question “Is it steep?” (as long as we understand 1 means Yes or true) …

steep <- slope > 20

… or this map that shows all areas that are steeper than 20 degrees and above 2000 m elevation:

plot(steep * (elev > 2000))

Sinking Cove, Tennessee is a karst valley system carved into the Cumberland Plateau, a nice place to see the use of a hillshade raster created from a digital elevation model using raster or terra functions for slope, aspect, and hillshade:

The raster package has a rich array of raster GIS operations, but these should all be incorporated in its replacement, terra, so we’ll use that. Note that the shade() function in the terra package wants radians for slope and aspect:

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
library(sf); library(tidyverse); library(tmap)
DEM <- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))
slope <- terrain(DEM, v='slope')
aspect <- terrain(DEM, v='aspect')
hillsh <- shade(slope/180*pi, aspect/180*pi, angle=40, direction=330)
# Need to crop a bit since grid north != true north 
bbox0 <- st_bbox(DEM)
xrange <- bbox0$xmax - bbox0$xmin
yrange <- bbox0$ymax - bbox0$ymin
bbox1 <- bbox0
crop <- 0.05
bbox1[1] <- bbox0[1] + crop * xrange # xmin
bbox1[3] <- bbox0[3] - crop * xrange # xmax
bbox1[2] <- bbox0[2] + crop * yrange # ymin
bbox1[4] <- bbox0[4] - crop * yrange # ymax
bboxPoly <- bbox1 %>% st_as_sfc() # makes a polygon
tm_shape(hillsh, bbox=bboxPoly) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) +
  tm_raster(palette=terrain.colors(24), alpha=0.5) +

6.5 ggplot2 for maps

The Grammar of Graphics is the gg of ggplot.

  • Key concept is separating aesthetics from data
  • Aesthetics can come from variables (using aes()setting) or be constant for the graph

Mapping tools that follow this lead

  • ggplot, as we have seen, and it continues to be enhanced
  • tmap (Thematic Maps) https://github.com/mtennekes/tmap Tennekes, M., 2018, tmap: Thematic Maps in R, Journal of Statistical Software 84(6), 1-39
ggplot(CA_counties) + geom_sf()

Try ?geom_sf and you’ll find that its first parameters is mapping with aes() by default. The data property is inherited from the ggplot call, but commonly you’ll want to specify data=something in your geom_sf call.

Another simple ggplot, with labels

ggplot(CA_counties) + geom_sf() +
  geom_sf_text(aes(label = NAME), size = 1.5)

and now with fill color

ggplot(CA_counties) + geom_sf(aes(fill = MED_AGE)) +
  geom_sf_text(aes(label = NAME), col="white", size=1.5)

Repositioned legend, no “x” or “y” labels

ggplot(CA_counties) + geom_sf(aes(fill=MED_AGE)) +
  geom_sf_text(aes(label = NAME), col="white", size=1.5) +
  theme(legend.position = c(0.8, 0.8)) +

Map in ggplot2, zoomed into two counties [air_quality]: (Toxic Release Inventory (TRI) Program, n.d.)

library(tidyverse); library(sf); library(iGIScData)
census <- st_make_valid(BayAreaTracts) %>%
   filter(CNTY_FIPS %in% c("013", "095"))
TRI <- TRI_2017_CA %>%
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
  st_join(census) %>%
  filter(CNTY_FIPS %in% c("013", "095"),
         (`5.1_FUGITIVE_AIR` + `5.2_STACK_AIR`) > 0)
bnd = st_bbox(census)
ggplot() +
  geom_sf(data = BayAreaCounties, aes(fill = NAME)) +
  geom_sf(data = census, color="grey40", fill = NA) +
  geom_sf(data = TRI) +
  coord_sf(xlim = c(bnd[1], bnd[3]), ylim = c(bnd[2], bnd[4])) +
  labs(title="Census Tracts and TRI air-release sites") +
  theme(legend.position = "none")

6.5.1 Rasters in ggplot2

Raster display in ggplot2 is currently a little awkward, as are rasters in general in the feature-dominated GIS world.

We can use a trick: converting rasters to a grid of points [marbles]:

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
elev <- raster(ex("elev.tif"))
trails <- st_read(ex("trails.shp"))
## Reading layer `trails' from data source `C:\Users\900008452\Documents\R\win-library\4.1\iGIScData\extdata\trails.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 8 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 481903.8 ymin: 4599196 xmax: 486901.9 ymax: 4603200
## Projected CRS: NAD83 / UTM zone 10N
elevpts = as.data.frame(rasterToPoints(elev))
ggplot() +
  geom_raster(data = elevpts, aes(x = x, y = y, fill = elev)) +
  geom_sf(data = trails)

terra : You can use the rasterVis::gplot function to plot a raster in ggplot, since that gplot is a wrapper around ggplot, but I haven’t figured out yet how to also plot vectors on top. So this works just for rasters …

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
elev <- rast(ex("elev.tif"))
gplot(elev) + geom_tile(aes(fill=value))

6.6 tmap

The tmap package provides some nice cartographic capabilities. We’ll use some of its capabilities, but for a more thorough coverave, see Geocomputation with R Chapter 8 “Making Maps with R”: https://geocompr.robinlovelace.net/adv-map.html

The basic building block is tm_shape(data) followed by various layer elements such as tm_fill(). The tm_shape function can work with either features or rasters.

tm_shape(world) + tm_fill() + tm_borders()

Color by variable [air_quality]

tm_shape(st_make_valid(BayAreaTracts)) + tm_fill(col = "MED_AGE")

tmap of sierraFeb with hillshade and point symbols [sierra]

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
library(terra) # alt: library(raster)
tmap_options(max.categories = 8)
sierra <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
hillsh <- rast(ex("ca_hillsh_WGS84.tif"))   # alt: hillsh <- raster(...)
bounds <- st_bbox(sierra)
  tm_raster(palette="-Greys",legend.show=FALSE,n=10) + tm_shape(sierra) + tm_symbols(col="TEMPERATURE",
     palette=c("blue","red"), style="cont",n=8) +
  tm_legend() + 

Note: “-Greys” needed to avoid negative image, since “Greys” go from light to dark, and to match reflectance as with b&w photography, they need to go from dark to light.


From a hydrologic and geochemical study of a fluviokarstic valley system in Tennessee (Jerry D. Davis and Brook 1993):

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
library(sf); library(tidyverse); library(readxl); library(tmap)
wChemData <- read_excel(ex("SinkingCove/SinkingCoveWaterChem.xlsx")) %>%
  mutate(siteLoc = str_sub(Site,start=1L, end=1L))
wChemTrunk <- wChemData %>% filter(siteLoc == "T") %>% mutate(siteType = "trunk")
wChemDrip <- wChemData %>% filter(siteLoc %in% c("D","S")) %>% mutate(siteType = "dripwater")
wChemTrib <- wChemData %>% filter(siteLoc %in% c("B", "F", "K", "W", "P")) %>% mutate(siteType = "tributary")
wChemData <- bind_rows(wChemTrunk, wChemDrip, wChemTrib)
sites <- read_csv(ex("SinkingCove/SinkingCoveSites.csv"))
wChem <- wChemData %>%
  left_join(sites, by = c("Site" = "site")) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
library(terra) # alt: library(raster)
DEM <- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))   # alt: raster(...)
slope <- terrain(DEM, v='slope')   # alt w/raster opt="slope"
aspect <- terrain(DEM, v='aspect')
hillsh <- shade(slope/180*pi, aspect/180*pi, 40, 330)  # alt (raster): hillShade
bounds <- st_bbox(wChem)
xrange <- bounds$xmax - bounds$xmin
yrange <- bounds$ymax - bounds$ymin
xMIN <- as.numeric(bounds$xmin - xrange/10)
xMAX <- as.numeric(bounds$xmax + xrange/10)
yMIN <- as.numeric(bounds$ymin - yrange/10)
yMAX <- as.numeric(bounds$ymax + yrange/10)
newbounds <- st_bbox(c(xmin=xMIN, xmax=xMAX, ymin=yMIN, ymax=yMAX), crs= st_crs(4326))
tm_shape(hillsh,bbox=newbounds) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) + tm_raster(palette=terrain.colors(24), alpha=0.5,legend.show=F) +
  tm_shape(wChem) + tm_symbols(size="TH", col="Lithology", scale=2, shape="siteType") +
  tm_layout(legend.position = c("left", "bottom")) +
  tm_graticules(lines=F) Comparable methods using terra and plot

[under construction]

The terra package promises to be a major improvement in the r spatial world, so I thought we might see if we could do something similar to what we did in tmap. Learn more about it at https://rspatial.org but we’ll briefly try some of what what we just did in tmap. We’ll start again with a spreadsheet of data and use some dplyr methods to put a data frame together, which will also include longitude & latitude variables:

ex <- function(dta){system.file("extdata",dta,package="iGIScData")}
library(tidyverse); library(readxl); library(tmap)
wChemData <- read_excel(ex("SinkingCove/SinkingCoveWaterChem.xlsx")) %>%
  mutate(siteLoc = str_sub(Site,start=1L, end=1L))
wChemTrunk <- wChemData %>% filter(siteLoc == "T") %>% mutate(siteType = "trunk")
wChemDrip <- wChemData %>% filter(siteLoc %in% c("D","S")) %>% mutate(siteType = "dripwater")
wChemTrib <- wChemData %>% filter(siteLoc %in% c("B", "F", "K", "W", "P")) %>% mutate(siteType = "tributary")
wChemData <- bind_rows(wChemTrunk, wChemDrip, wChemTrib)
sites <- read_csv(ex("SinkingCove/SinkingCoveSites.csv"), show_col_types = FALSE)
wChemJoin <- wChemData %>%
  left_join(sites, by = c("Site" = "site"))

Then we’ll use terra::vect() to create vector (spatVector) data (note that terra likes crs specified in text):

lonlat <- cbind(wChemJoin$longitude, wChemJoin$latitude)
wChemV <- vect(lonlat, atts=wChemJoin, crs="+proj=longlat +datum=WGS84")
plot(wChemV, col=factor(wChemV$Lithology))

… and raster (spatRaster) data, same as we did previously:

DEM <- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))
slope <- terrain(DEM, v='slope')
aspect <- terrain(DEM, v='aspect')
hillshT <- shade(slope/180*pi, aspect/180*pi, 40, 330)
wChemUTM <- terra::project(wChemV, crs(DEM))  # since the crs for DEM is UTM
xrange <- xmax(wChemUTM)-xmin(wChemUTM)
yrange <- ymax(wChemUTM)-ymin(wChemUTM)
xMIN <- xmin(wChemUTM) - xrange/10
xMAX <- xmax(wChemUTM) + xrange/10
yMIN <- ymin(wChemUTM) - yrange/10
yMAX <- ymax(wChemUTM) + yrange/10
plot(hillshT, col = gray.colors(12), xlim=c(xMIN, xMAX), ylim=c(yMIN, yMAX))
points(wChemUTM, col=factor(wChemV$Lithology), cex=wChemV$TH/100)

[Not finished – more to figure out to create what we did in tmap – maybe have to convert some data – but we’ve created some of the elements]

6.7 Interactive Maps

The word “static” in “static maps” isn’t something you would have heard in a cartography class 30 years ago, since essentially all maps then were static. Very important in designing maps is considering your audience and their perception of symbology on a static medium.

  • Figure-to-ground relationships assume “ground” is a white piece of paper (or possibly a standard white background in a pdf), so good cartographic color schemes tend to range from light for low values to dark for high values.
  • Scale is fixed and there are no “tools” for changing scale, so a lot of attention must be paid to providing scale information.
  • Similarly, without the ability to see the map at different scales, inset maps are often needed to provide context.

Interactive maps change the game in having tools for changing scale, and always being “printed” on a computer or device where the color of the background isn’t necessarily white. We are increasingly used to using interactive maps on our phones or other devices, and often get frustrated not being able to zoom into a static map.

A widely used interactive mapping system is Leaflet, but we’re going to use tmap to access Leaflet behind the scenes and allow us to create maps with one set of commands. The key parameter needed is tmap_mode which must be set to “view” to create an interactive map.


With an interactive map, we do have the advantage of a good choice of base maps and the ability to resize and explore the map, but symbology is more limited, mostly just color and size, with only one variable in a legend.

bounds <- st_bbox()
wChem2map <- filter(wChem, Month == 8)
minVal <- min(wChem2map$TH); maxVal <- max(wChem2map$TH)
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="siteType", size="TH", scale=2) +
  tm_layout(title=paste("Total Hardness ",as.character(minVal),"-",as.character(maxVal)," mg/L", sep=""))
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="Lithology", size="TH", scale=2)


tm_shape(st_make_valid(BayAreaTracts)) + tm_fill(col = "MED_AGE", alpha = 0.5)