Thematic Maps with R

R User Group at the Harvard Data Science Initiative

Timothée Giraud

Apr 20, 2023

Thematic maps with R?

Reproducible cartography?

Maps, as other graphics or statistical outputs, are part of scientific studies.

The problem for cartographers

Many maps produced in an academic context are made with a set of software products that slices the cartographic process.

The problem for cartographers

Fully reproducible maps should be associated with code and data.

Giraud and Lambert (2017)

The problem for R users

A solution

Simplify the map making process in an unified workflow thanks to R and its spatial ecosystem !

The R spatial ecosystem

A bit of history

  • before 2003 : spatial, sgeostat, splancs, akima, geoR, spatstat, spdep, maptools.
  • 2003 : rgdal (Roger Bivand, Tim Keitt & Barry Rowlingson), interface between R and GDAL and PROJ4
  • 2005 : sp (Edzer Pebesma & Roger Bivand), classes methods for spatial objects, quickly adopted.
  • 2008 : sp surpport in ggplot2 (Hadley Wickham)
  • 2010 : rgeos (Roger Bivand & Colin Rundel), interface between R & GEOS.
  • 2010 : raster (Robert J. Hijmans), support for raster data
  • 2016 : sf (Edzer Pebesma), supersedes sp, rgdal & rgeos
  • 2018 : stars (Edzer Pebesma), supersedes raster
  • 2020 : terra (Robert J. Hijmans) also supersedes raster

Thematic cartography

  • 2014 : tmap (Martijn Tennekes)
  • 2017 : ggplot2 + ggspatial (Dewey Dunnington)
  • 2021 : mapsf (T. Giraud) supersedes cartography

Interactive cartography

  • 2015 : leaflet (Joe Cheng et al.), relies on the leaflet javascript library.
  • 2015 : mapview (Tim Appelhans et al.), relies on the leaflet package.
  • 2018 : mapdeck (David Cooley), relies on Mapbox GL & Deck.gl libraries.

Spatial statistics

  • spatstat : Point pattern analysis
  • gstat : Variogram & Krigeage
  • rgeoda : Geoda with R
  • GWmodel, spgwr : Geographically Weighted Models
  • …

sf, corner stone of the R spatial ecosystem

Interface between R and extensively used libraries in GIS:

  • GDAL - Geospatial Data Abstraction Library
  • PROJ - Coordinate Transformation Software
  • GEOS - Geometry Engine - Open Source

GitHub repo GitHub repo GitHub repo

Pebesma and Bivand (2023)

sf

library(sf)
mtq <- st_read(dsn = "data/mtq.gpkg",layer = "mtq")
Reading layer `mtq' from data source 
  `/home/tim/Documents/prz/RUG_HDSI/data/mtq.gpkg' using driver `GPKG'
Simple feature collection with 34 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 690574 ymin: 1592536 xmax: 735940.2 ymax: 1645660
Projected CRS: WGS 84 / UTM zone 20N
mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), cex = 1.2, col = "red", pch = 20, add = TRUE)

mat <- st_distance(x = mtq_c, y = mtq_c)
mat[1:5, 1:5]
Units: [m]
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35316.29  3019.561 12247.132 17181.321
[2,] 35316.293     0.00 38303.361 25478.399 18597.357
[3,]  3019.561 38303.36     0.000 15099.121 20200.611
[4,] 12247.132 25478.40 15099.121     0.000  7149.189
[5,] 17181.321 18597.36 20200.611  7149.189     0.000
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), lwd = 2, border = "red", add = TRUE)

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq_b), col = "grey", lwd = 2, border = "red")
plot(st_geometry(mtq), col = "lightblue", add = TRUE)
plot(st_geometry(mtq_u), lwd = 2, add = TRUE)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = TRUE)
plot(st_geometry(mtq_z), lwd = 2)

mtq_c |> 
  st_union() |> 
  st_voronoi() |> 
  st_collection_extract("POLYGON") |> 
  st_intersection(mtq_u) |> 
  st_sf() |> 
  st_join(mtq_c, st_intersects) |>
  st_cast("MULTIPOLYGON") |>
  st_geometry() |>
  plot(col = "ivory4")

plot(st_geometry(mtq))

plot(mtq)

mapsf

mapsf logo

Typical workflow

mf_map()

mf_map() is the main function of the package.

mf_map(x = objet_sf,
var = "variable",
type = "map type",
...)

Various map types

Map layout

Example

Base map

mtq <- mf_get_mtq()
# Plot the base map
mf_map(x = mtq)

Choropleth map

mtq <- mf_get_mtq()
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro")

Parameters customisation

mtq <- mf_get_mtq()
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright")

Map layout

mtq <- mf_get_mtq()
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright")
# Plot a layout elements
mf_title("Wealth in Martinique, 2015")
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
mf_scale(size = 5)
mf_arrow('topleft')

Theme and layout

mtq <- mf_get_mtq()
# Set a theme 
mf_theme(x = "green")
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright")
# Plot a layout elements
mf_title("Wealth in Martinique, 2015")
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
mf_scale(size = 5)
mf_arrow('topleft')

Inset

mtq <- mf_get_mtq()
# Set a theme 
mf_theme(x = "green")
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright")
# Add an inset world map
mf_inset_on(x = "worldmap", pos = "right")
mf_worldmap(mtq, col = "#0E3F5C")
mf_inset_off()
# Plot a layout elements
mf_title("Wealth in Martinique, 2015")
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
mf_scale(size = 5)
mf_arrow('topleft')

Margins

mtq <- mf_get_mtq()
# Set a theme 
mf_theme(x = "green")
# Start an empty map with extra margins
mf_init(x = mtq, expandBB = c(0,0,0,.3))
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright",
       add = TRUE)
# Add an inset world map
mf_inset_on(x = "worldmap", pos = "right")
mf_worldmap(mtq, col = "#0E3F5C")
mf_inset_off()
# Plot a layout elements
mf_title("Wealth in Martinique, 2015")
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
mf_scale(size = 5)
mf_arrow('topleft')

Export

mtq <- mf_get_mtq()
# Set a theme 
mf_theme(x = "green")
# Start the export with extra margins
mf_export(x = mtq, filename = "img/map.svg", 
          width = 5, expandBB = c(0,0,0,.3))
# Plot a choropleth map
mf_map(x = mtq, var = "MED", type = "choro",
       pal = "Dark Mint",
       breaks = "quantile",
       nbreaks = 6,
       leg_title = "Median Income\n(euros)",
       leg_val_rnd = -2,
       leg_pos = "topright",
       add = TRUE)
# Add an inset world map
mf_inset_on(x = "worldmap", pos = "right")
mf_worldmap(mtq, col = "#0E3F5C")
mf_inset_off()
# Plot a layout elements
mf_title("Wealth in Martinique, 2015")
mf_credits("T. Giraud\nSources: INSEE & IGN, 2018")
mf_scale(size = 5)
mf_arrow('topleft')
# close export
dev.off()

Themes

Vignette : How to Use Themes

Raster

library(terra)
osm <- rast("data/osm.tiff")
mf_theme("default", mar = c(0,0,0,0), 
         inner = TRUE, pos = "right")
mf_raster(osm)
mf_title("mf_raster()")
mf_credits(
  paste0("Map tiles by Stamen Design ", 
         "CC BY 3.0 — Map data ", 
         "© OpenStreetMap contributors"))

Documentation

A website

Documentation

Documentation

A cheat sheet

Development

GitHub repo GitHub repo GitHub repo

Other packages

tanaka

Also called “relief contours”, “illuminated contours” or “shaded contour lines”, the Tanaka method enhances the representation of topography on a map by using shaded contour lines. The result is a 3D-like map.

GitHub repo GitHub repo

tanaka

library(tanaka)
library(terra)
ras <- rast(system.file("tif/elev.tif", package = "tanaka"))
plot(ras)
tanaka(ras, breaks = seq(80,400,20), 
       legend.pos = "topright", legend.title = "Elevation\n(meters)")

tanaka

tanaka can be used with non-topographical data:

linemap

linemap displays a map made of lines using a data frame of gridded data.


GitHub repo GitHub repo

maptiles

To create maps from tiles, maptiles downloads, composes and displays tiles from a large number of providers (e.g. OpenStreetMap, Stamen, Esri, CARTO, or Thunderforest).


GitHub repo GitHub repo

maptiles

library(sf)
library(maptiles)
# import North Carolina counties
nc_raw <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
# Project to EPSG:3857
nc <- st_transform(nc_raw, "EPSG:3857")
# dowload tiles and compose raster (SpatRaster)
nc_osm <- get_tiles(nc, crop = TRUE)
# display map
mf_theme(mar = c(0,0,0,0))
mf_raster(nc_osm)
mf_map(nc, col = NA, add = TRUE)
# add credit
mf_credits(txt = get_credit("OpenStreetMap"))

maptiles

Mini maps of most of the tiles providers

cartogram

GitHub repo GitHub repo

cartogram

fisheye

Transform base maps using log-azimuthal projection

GitHub repo GitHub repo

fisheye

library(fisheye)
library(mapsf)
# Import dataset
ncraw <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
nc <- st_transform(ncraw, 3857)
mf_map(nc, col ="grey90")
mf_map(nc[51, ], add = TRUE, col = "grey40")
mf_title("Original Map")
# transform the basemap
nc_fe  <- fisheye(nc, centre = nc[51, ])
mf_map(nc_fe, col ="grey90")
mf_map(nc_fe[51, ], add = TRUE, col = "grey40")
mf_title("Log-Azimuthal Projection")

fisheye

popcircle

This one-function package computes circles with areas scaled to a variable and displays them using a compact layout. Original polygons are (roughly) scaled to fit inside these circles.

GitHub repo

remotes::install_github("rCarto/popcircle")

popcircle

library(rnaturalearth)
library(sf)
library(wbstats)
library(popcircle)
library(mapsf)

# Get countries
ctry <- ne_countries(scale = 50, returnclass = "sf")
ctry <- st_transform(ctry, "ESRI:54030")

# Only keep the largest polygons of multipart polygons for a few countries
# (e.g. display only continental US)
frag_ctry <- c("US", "RU", "FR", "IN", "ES", "NL", "CL", "NZ", "ZA")
largest_ring = function(x) {
  x$ids <- 1:nrow(x)
  pols = st_cast(x, "POLYGON", warn = FALSE)
  spl = split(x = pols, f = pols$ids)
  do.call(rbind, (lapply(spl, function(y)
    y[which.max(st_area(y)), ])))
}
st_geometry(ctry[ctry$iso_a2 %in% frag_ctry, ]) <-
  st_geometry(largest_ring(ctry[ctry$iso_a2 %in% frag_ctry, ]))


# Get and merge data
data_pop <-
  wb_data(indicator = "SP.POP.TOTL",
          start_date = 2021,
          end_date = 2021)
ctry_pop <-
  merge(ctry[, "iso_a2"], data_pop, by.x = "iso_a2", by.y = "iso2c")

# Computes circles and polygons
res_pop <- popcircle(x = ctry_pop, var = "SP.POP.TOTL")
circles_pop <- res_pop$circles
shapes_pop <- res_pop$shapes


# Create the figure
mf_theme(mar = c(0, 0, 0, 0),
         bg = "#e6ebe0",
         fg = "grey50")
mf_export(
  x = circles_pop,
  filename = "pop.png",
  width = 800,
  res = 100
)

# display circles and polygons
mf_map(
  circles_pop,
  col = "#9bc1bc",
  border = "white",
  add = T
)
mf_map(
  shapes_pop,
  col = "#ed6a5a95",
  border = "#ed6a5a",
  add = TRUE,
  lwd = .3
)
# labels
circles_pop$lab <- paste0(circles_pop$country,
                          '\n',
                          round(circles_pop$SP.POP.TOTL / 1000000))
mf_label(
  x = circles_pop[1:36, ],
  var = "lab",
  halo = TRUE,
  overlap = T,
  pos = 3,
  cex = seq(1, 0.4, length.out = 36),
  r = .15
)
# title
mtext(
  "Population",
  side = 3,
  adj = 0.01,
  padj = -1.5,
  col = "grey50",
  cex = 2
)
mtext(
  "Millions of inhabitants\nin 2021",
  side = 3,
  adj = 0.01,
  padj = 0.4,
  col = "grey50",
  cex = 1.2
)
# cerdits
mf_credits(txt = "T. Giraud, 2023 - World Development Indicators, 2023",
           pos = "bottomright")
dev.off()

spikemap

Lazaro Gamio, Karen Yourish and Bill Marsh, 2020-04-08

spikemap

It is possible to map quantities with circles, squares or other simple geometric symbols, spikemap uses spikes.

GitHub repo

remotes::install_github("rCarto/spikemap")

spikemap

library(sf)
library(spikemap)
library(mapsf)
# import the dataset from the package
com <- st_read(system.file("gpkg/com.gpkg", package="spikemap"))
# theme 
mf_theme(mar = c(0,0,0,0), bg = "#e1e5eb", fg =  "grey30")
# save figure as spiky.png in img folder
mf_export(x = com, filename = "img/spiky.png", width = 1000, res = 100)
# plot the base map
mf_map(com, col="#99aed1", border = "#e1e5eb", lwd = 0.2, add = T)
# display spikes for municipalities under 1000 inhabitants.
# use fixmax arg to allow multiple spike plots with the same scale.
spikemap(x = com[com$pop<=1000, ], var = "pop",
         inches = 2.3, fixmax = 500000,
         col = "#ffffff90", border = "#94000090",  lwd = .5,
         legend.pos = "x")
# display spikes for other municipalities
# use locator() to pick a place for the legend or use "bottomleft".
spikemap(x = com[com$pop>1000, ], var = "pop",
         inches = 2.3, fixmax = 500000,
         col = "#ffffff", border = "#940000", lwd = 1.1,
         legend.pos = c(799307.2, 6128000),
         legend.title.txt = "Population",
         legend.values.rnd = -3)

# get the tips of the spikes
lbl <- spikelabel(x = com, var = "pop",
                  inches = 2.3, fixmax = 500000)
lbl <- lbl[order(lbl$pop, decreasing = T),]
# display only the 12 first, use various cex and halo
mf_label(lbl[1:12,], var = "name",
           pos = 3, offset = .5,
           halo = T, bg = "#99aed150",
           cex = c(1.3, 1.1, 1, rep(.8,12)),
           col = "grey30")

# add scale bar, north arrow, title, sources...
mf_scale(size = 20, pos= c(629638.7 ,6136862.3 ), lwd = 1)
mf_arrow(pos = "topright", col = "grey60")
mf_credits(paste0("ADMIN EXPRESS COG édition 2019, IGN\n",
                  "T. Giraud, 2020 | spikemap 0.1.0"))
mf_title(txt = "Population \nin Occitanie", 
         inner = T, line = 8,  
         cex = 2.5, font = 3, fg = "grey30", bg = NA)
dev.off()

spikemap

spikemap

more…

Lambert (2021)

more…

Giraud and Lambert (2019)

more…

Le Goix et al. (2021)

Resources

Geocomputation with R

Lovelace, Nowosad, and Muenchow (2019)

Spatial Data Science with applications in R

Pebesma and Bivand (2023)

Spatial Data Science with R and “terra”

Robert J. Hijmans

Progress in the R ecosystem for representing and handling spatial data

Bivand (2021)

Cartographie avec R(fr)

Giraud and Pecout (2023a)

Géomatique avec R(fr)

Giraud and Pecout (2023b)

Misc

Thank you

rcarto.github.io/RUG_HDSI

@rcarto@fosstodon.org

rcarto.github.io

@rcarto

References

Bivand, Roger S. 2021. “Progress in the R Ecosystem for Representing and Handling Spatial Data.” Journal of Geographical Systems 23 (4): 515–46. https://doi.org/10.1007/s10109-020-00336-0.
Giraud, Timothée, and Nicolas Lambert. 2017. “Reproducible Cartography.” In Advances in Cartography and GIScience, 173–83. Springer International Publishing. https://doi.org/10.1007/978-3-319-57336-6_13.
———. 2019. “Reproducible Workflow for Cartography - Migrants Deaths in the Mediterranean.” In Proceedings of the ICA, 2:1–7. Copernicus GmbH. https://doi.org/10.5194/ica-proc-2-38-2019.
Giraud, Timothée, and Hugues Pecout. 2023a. “Cartographie avec R.” https://doi.org/10.5281/zenodo.7528161.
———. 2023b. “Géomatique avec R.” https://doi.org/10.5281/zenodo.7528145.
Lambert, Nicolas. 2021. “Le nouveau rideau de fer: Un Exemple de Carte En 2.5D.” FR2007 CIST. https://doi.org/10.48645/a4ra-yr11.
Le Goix, Renaud, Ronan Ysebaert, Timothée Giraud, Marc Lieury, Guilhem Boulay, Mathieu Coulon, Sébastien Rey-Coyrehourcq, et al. 2021. “Unequal Housing Affordability Across European Cities. The ESPON Housing Database, Insights on Affordability in Selected Cities in Europe.” Cybergeo, April. https://doi.org/10.4000/cybergeo.36478.
Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2019. Geocomputation with R. CRC Press. https://r.geocompx.org/.
Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With Applications in R. CRC Press. https://r-spatial.org/book/.