Week 2: How Do Maps Work?

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2024

Jeff Jacobs

jj1088@georgetown.edu

Wednesday, September 4, 2024

Vector \(\rightarrow\) Raster Data

End of Last Week: Vector Data

library(sf)
library(tidyverse)
# Load DC tracts data
dc_sf_fpath <- "data/DC_Census_2020/Census_Tracts_in_2020.shp"
dc_sf <- st_read(dc_sf_fpath, quiet = TRUE)
cols_to_keep <- c("OBJECTID", "TRACT", "GEOID", "ALAND", "AWATER", "STUSAB", "SUMLEV", "GEOCODE", "STATE", "NAME", "POP100", "HU100", "geometry")
dc_sf <- dc_sf |> select(cols_to_keep)
# We can extract the geometry with the st_geometry function
dc_geo <- st_geometry(dc_sf)
# And plot the geometry with base R's plot() function
plot(dc_geo, mar = c(0,0,0,0))

Now with ggplot!

ggplot2::theme_set(ggplot2::theme_classic())
dc_sf |>
  ggplot() +
  geom_sf()

Raster Data

• Each DC Census Tract has its own (odd) shape, which can be described by discrete coordinates forming a POLYGON
• For geospatial analysis, however, we often need to compute over evenly-spaced grids rather than this odd collection of shapes
• POLYGONs may make sense for demographers, but how about someone studying air pollution in DC? (Smog, for example, does not confine itself to census tracts!)

Step 1: Union of All Tracts

dc_union_sf <- sf::st_union(dc_sf)
dc_union_sf |>
  ggplot() +
  geom_sf()

Step 2: Rasterize (terra)

library(terra)
dc_SpatVector <- terra::vect(dc_union_sf)
rast_template <- rast(ext(dc_SpatVector), resolution = 1000, crs = crs(dc_SpatVector))
dc_SpatRaster <- terra::rasterize(dc_SpatVector, rast_template)
dim(dc_SpatRaster)

[1] 29 23  1

plot(dc_SpatRaster)

And with tidyterra

library(tidyterra)
ggplot(data=dc_SpatRaster, aes(fill=layer)) +
  geom_spatraster(data=dc_SpatRaster)

Rasters From Scratch

Welcome to Gridtown!

set.seed(6805)
library(terra)
gridtown <- terra::rast(
  nrows = 4, ncols = 4,
  xmin = 0, xmax = 4, ymin = 0, ymax = 4,
  vals = sample(1:16)
)
plot(gridtown)
text(
  gridtown,
  labels=1:16,
  halo=TRUE, hc="black", col="white", hw=0.2
)

Figure 1: Gridtown Indices

• Raster indices vs. values: The above plot displays indices for each cell: since a raster is a regular grid, can achieve memory-efficient representation with a single index (rather than, e.g., \((x, y)\) coords). But what we really care about are…

Raster Layer Values

plot(gridtown)
text(gridtown, halo=TRUE, hc="black", col="white", hw=0.2)

Figure 2: Gridtown Values

First Things First: Coordinates

• We need a way to unambiguously specify where things are on the earth

Latitude and Longitude are Angles!

From Krygier and Wood (2016)

Angular Distance vs. Travel Distance

The Earth’s “width” is slightly greater than its “length” 😰

From Wikimedia Commons

Projections

Smooshing 3D into 2D

From Monmonier (2018)

Avoid Getting Lost in the Sauce

How To Avoid Getting Lost in the Sauce

Tissot Circles: Imagine infinitely small ellipses placed at regular intervals on the curved surface of the earth. Imagine these ellipses being projected along with the earth’s surface. When scaled up, changes in the ellipses show the location and quality of distortions on the projected map.

The Things We Put On Maps

• Scales
• Points
• Lines
• Areas

Scales (w/Rules of Thumb)

From Krygier and Wood (2016)

The Three Key Elements

From Monmonier (2018)

Thought Experiment: Zooming In and Out

• What happens to points, lines, and areas?

What Can Happen to Points?

From Monmonier (2018)

What Can Happen to Lines?

From Monmonier (2018)

What Can Happen to Lines?

Selection

In any map, most features from the human or natural environment are eliminated! (Krygier and Wood 2016)

Simplification

Within the features that do remain, eliminate details that are unnecessary with respect to audience/context (Krygier and Wood 2016)

What Can Happen to Areas?

From Monmonier (2018)

Areal Dimension Change

From Krygier and Wood (2016)

What Is Gained?

From Monmonier (2018)

What Is Lost?

From Monmonier (2018)

Who Is Affected?

From Monmonier (2018)

Choropleths for Good and Evil

Having Domain Knowledge of a Region

From XKCD

Crime in Mongolia

From Reddit

Population of Mongolia

From Wikimedia Commons

Exhibit A

From Monmonier (2018)

Exhibit B

From Monmonier (2018)

Finding Geospatial Data

• Draws heavily on Chapter 6: “R packages to download open spatial data”, in Moraga (2018), Spatial Statistics for Data Science: Theory and Practice with R

Country Boundaries

Key package: rnaturalearth

library(rnaturalearth)
library(sf)
library(ggplot2)
library(viridis)
library(patchwork)
map1 <- ne_countries(type = "countries", country = "Germany", scale = "medium", returnclass = "sf")
map2 <- rnaturalearth::ne_states("Germany", returnclass = "sf")
(ggplot(map1) + geom_sf()) + (ggplot(map2) + geom_sf())

Climate Data

Key package: geodata

library(geodata)
d <- worldclim_country(country = "Jamaica", var = "tmin",
                       path = tempdir())
terra::plot(mean(d), plg = list(title = "Min. temperature (C)"))

Elevation

Key packages: rnaturalearth + elevatr

library(rnaturalearth)
library(elevatr)
library(terra)
map <- ne_countries(type = "countries", country = "Switzerland",
                    scale = "medium", returnclass = "sf")
d <- get_elev_raster(locations = map, z = 9, clip = "locations")
terra::plot(rast(d), plg = list(title = "Elevation (m)"))

Street Maps

library(osmdata)
placebb <- getbb("Barcelona")
hospitals <- placebb %>% opq() %>%
  add_osm_feature(key = "amenity", value = "hospital") %>%
  osmdata_sf()
motorways <- placebb %>% opq() %>%
  add_osm_feature(key = "highway", value = "motorway") %>%
  osmdata_sf()
library(leaflet)
leaflet() |>  addTiles() |>
  addPolylines(data = motorways$osm_lines, color = "black") |>
  addPolygons(data = hospitals$osm_polygons,
              label = hospitals$osm_polygons$name)

World Bank Dataverse

Key package: wbstats

library(wbstats)
d <- wb_data(indicator = "MO.INDEX.HDEV.XQ",
             start_date = 2011, end_date = 2011)
library(rnaturalearth)
library(mapview)
map <- ne_countries(continent = "Africa", returnclass = "sf")
map <- dplyr::left_join(map, d, by = c("iso_a3" = "iso3c"))
mapview(map, zcol = "MO.INDEX.HDEV.XQ")

References

Krygier, John, and Denis Wood. 2016. Making Maps, Third Edition: A Visual Guide to Map Design for GIS. Guilford Publications.

Monmonier, Mark. 2018. How to Lie with Maps. University of Chicago Press.