Two centuries of American settlement

HISDAC-US records, for every 250 m cell of the conterminous United States, the year its first structure appears in the property record. One raster, 1810–2020. Because a cell becomes built and stays built, the series is cumulative - which is what lets the whole animation come from a single layer and a moving threshold, rather than from 43 separate rasters.

Load the raster

Opens the GeoTIFF as a SpatRaster - a lazy handle, not 214 million values in memory - and silences terra’s progress bars, which otherwise print a row of pipes under every chunk.

rast() - open the GeoTIFF; terraOptions(progress = 0) - silence progress bars

library(terra)
library(here)

terraOptions(progress = 0)

FBUY <- here("dataverse-files", "v2", "FBUY.tif")
stopifnot(file.exists(FBUY))

Constrain, downsample, and build a land mask

clamp() sends anything outside 1810–2020 to NA. aggregate(fun = "min") reduces 250 m cells to 2.5 km blocks, keeping the earliest year in each - the minimum, never the mean, because averaging build years produces dates that are in no property record. rasterize() burns the state polygons onto that same grid so the two can never drift apart.

clamp() → aggregate(fun = ‘min’) → rasterize() - constrain, downsample, mask

r <- rast(FBUY)
r <- clamp(r, 1810, 2020, values = FALSE)

a <- aggregate(r, fact = params$agg_factor, fun = "min", na.rm = TRUE)

states <- rnaturalearth::ne_states(country = "United States of America",
                                   returnclass = "sf")
states <- states[!states$postal %in% c("AK", "HI"), ]
sv     <- vect(sf::st_transform(states, crs(a)))
land   <- rasterize(sv, a)

The animation

Draw one frame

ifel() is the whole trick: a vectorised if-else evaluated on the raster itself, so each frame is a masked to cells built at or before year y - no loop over cells, no melting to a data frame. breaks pins values to colour bins explicitly; range fails in the 1810 frame, where a single distinct value gives terra nothing to build a scale from. The HUD is drawn on a fresh 0–1 canvas because par("usr") describes the map’s coordinates, not the visible frame.

ifel() + breaks - threshold the cumulative raster and pin the colour bins

pal <- colorRampPalette(c("#C2410C", "#DB6A16", "#E9A23B",
                          "#F2C36F", "#F9DFAB", "#FFF3DC"))(100)
brk <- seq(1810, 2020, length.out = length(pal) + 1)

draw <- function(y) {
  op <- par(bg = "#10131A", xpd = FALSE, new = FALSE)
  on.exit(par(op), add = TRUE)

  plot(land, col = "#1B2130", legend = FALSE, axes = FALSE, mar = c(0, 0, 4, 0))
  plot(ifel(a <= y, a, NA), col = pal, breaks = brk, add = TRUE, legend = FALSE)

  par(new = TRUE, mar = c(0, 0, 0, 0), xpd = NA)
  plot.new()
  plot.window(xlim = c(0, 1), ylim = c(0, 1), xaxs = "i", yaxs = "i")

  text(0.02, 0.960, "Two centuries of American settlement",
       adj = c(0, 1), cex = 1.9, font = 2, col = "#E8E6E1")
  text(0.02, 0.920, "Year each cell was first built on, where records exist - HISDAC-US, 1810–2020",
       adj = c(0, 1), cex = 1.05, col = "#8A8F9A")

  text(0.02, 0.155, y, adj = c(0, 0), cex = 3.4, font = 2, col = "#E8E6E1")
  text(0.02, 0.115, "First built", adj = c(0, 0), col = "#A8AEB8", cex = 1.25)

  bx <- seq(0.02, 0.26, length.out = length(pal) + 1)
  rect(head(bx, -1), 0.075, bx[-1], 0.100, col = pal, border = NA)
  text(0.02, 0.040, "1810", adj = c(0, 0), col = "#A8AEB8", cex = 1.3)
  text(0.26, 0.040, "2020", adj = c(1, 0), col = "#A8AEB8", cex = 1.3)

  text(0.98, 0.040, "Ahn, Leyk, Uhl & McShane (2023) · CC BY 4.0",
       adj = c(1, 0), col = "#5C6472", cex = 0.95)

  segments(0.02, 0.018, 0.98, 0.018, col = "#2A3140", lwd = 4)
  segments(0.02, 0.018, 0.02 + 0.96 * (y - 1810) / (2020 - 1810), 0.018,
           col = "#E9A23B", lwd = 4)
}

Encode the frames to video

Steps annually rather than every five years. The data is unchanged - still 5-year buckets, each simply held for five frames - but 43 frames at 10 fps was over in 4.3 seconds, before a scrolling viewer had registered it was moving. av_capture_graphics() opens a graphics device, runs the loop, and encodes straight to H.264 without ever writing PNGs to disk.

av_capture_graphics() - render 241 frames straight to H.264

frames <- c(seq(1810, 2020, by = 1), rep(2020, 30))

invisible(av::av_capture_graphics(
  for (y in frames) draw(y),
  output    = here("settlement.mp4"),
  width     = 1200,
  height    = 860,
  framerate = params$fps
))

The finished map

draw(2020) - render the last frame as a static figure

draw(2020)

Final frame, 2020. The full animation is written to settlement.mp4 in the repository root.

Credit

Ahn, Leyk, Uhl & McShane (2023), Gridded Settlement Data for the Conterminous United States Over 200 Years, Harvard Dataverse, doi:10.7910/DVN/HHFM5E. Derived from Zillow ZTRAX, Parcel Atlas and BuildZero Open City Model. Dates are when a property record says a structure was built, not when a place was settled.