Data Science II

Session 5: Rasters — Continuous Surfaces and Zonal Statistics

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

What You’ll Learn Today

  1. What a raster is, and how it differs from everything we have used so far
  2. Reading, inspecting, and mapping rasters with terra
  3. Raster algebra — arithmetic across layers
  4. extract — the operation that turns a raster into a column in your table
  5. Turning points into a raster of counts

Today’s data: temperature, elevation, and rainfall across Mexico.

Before We Start

#Step1: Libraries for today
library(terra)         # the modern raster package
library(sf)            # vectors, as before
library(dplyr)
library(ggplot2)
library(tidyterra)     # geom_spatraster, for ggplot
library(exactextractr) # area-weighted zonal statistics

sf_use_s2(FALSE)

New data in ../data/raster/ — three small GeoTIFFs, 320 KB in total.

Where We Left Off

Four sessions of vectors: points, lines, polygons. Discrete objects with sharp edges.

That works for a hospital, a road, a municipio.

It works badly for temperature, rainfall, elevation, population density, night-time lights — things that exist everywhere and change gradually.

For those you need the other data model.

Part A: What a Raster Is

Vector and Raster

The same river as a vector (left) and as a raster (right).

A Raster Is a Matrix with a Location

Strip away the jargon and a raster is three things:

  1. A matrix of numbers — the cell values
  2. An extent — the bounding box those cells cover
  3. A CRS — what those coordinates mean

From extent and matrix dimensions you get resolution: how much ground one cell covers.

Note

A raster has no “features” and no attribute table. There is nothing to join to by name — only by location.

Resolution Is a Choice

The same extent can be divided into few big cells or many small ones.

Resolution One cell is Good for
10 arc-min (today) ~18 km national climate
1 km 1 km² municipio analysis
30 m (Landsat) a city block land cover
10 cm (drone) a paving stone infrastructure

Finer is not better — it is bigger. Halving the cell size quadruples the file.

Raster File Formats

Format Extension Notes
GeoTIFF .tif the default; universally readable
Cloud-Optimised GeoTIFF .tif a GeoTIFF you can read partially, over HTTP
NetCDF .nc many dimensions — climate time series
ASCII grid .asc plain text, huge, avoid

Use GeoTIFF unless something forces you not to.

Which R Package?

Package Status
raster the old standard — retired, still everywhere online
terra its replacement, by the same author. Use this.
stars for data cubes: many dimensions, time series
exactextractr one job, done fast: zonal statistics

Warning

Just like sp in Session 2, most raster answers you find online are written for raster. terra uses rast() where raster used raster(), and objects are SpatRaster, not RasterLayer.

Today’s Data: WorldClim

WorldClim 2.1 — global climate surfaces interpolated from weather stations, free, at several resolutions.

  • mexico_tavg_monthly.tif — mean temperature, 12 layers, one per month
  • mexico_elevation.tif — metres above sea level
  • mexico_precip_annual.tif — annual rainfall in mm

Already cropped to Mexico. The global download is 37 MB; cropped, all three come to 320 KB.

https://worldclim.org

Part B: Reading and Inspecting

Reading a Raster

#Step1: rast() reads any raster format terra understands
tavg <- rast("../data/raster/mexico_tavg_monthly.tif")
tavg
class       : SpatRaster
size        : 112, 192, 12  (nrow, ncol, nlyr)
resolution  : 0.1666667, 0.1666667  (x, y)
extent      : -118.5, -86.5, 14.33333, 33  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source      : mexico_tavg_monthly.tif
names       :       Jan,       Feb,       Mar,       Apr,       May,       Jun, ...
min values  :     3.154,   4.25175,   5.23275,   7.44325,    9.8045,   9.08475, ...
max values  : 27.126429, 27.564501, 29.026001, 31.155251, 32.101749, 31.525249, ...

What That Told You

  • size 112, 192, 12 — rows, columns, and 12 layers
  • resolution 0.1667 degrees — about 18 km at this latitude
  • extent — the bounding box, in CRS units
  • coord. ref. — EPSG:4326, as usual
  • names — Jan through Dec; a layer is like a column

Tip

print() on a SpatRaster does not print the data — it prints the description. With 21,504 cells per layer that is a mercy.

Properties, One at a Time

dim(tavg)          # rows, columns, layers
[1] 112 192  12
res(tavg)          # cell size, in CRS units
[1] 0.1666667 0.1666667
nlyr(tavg)         # how many layers
[1] 12
names(tavg)[1:4]
[1] "Jan" "Feb" "Mar" "Apr"
ext(tavg)
SpatExtent : -118.5, -86.5, 14.333333333333334, 33 (xmin, xmax, ymin, ymax)

Getting at the Values

#Step1: One layer behaves like a matrix
summary(values(tavg[["Jan"]]))
      Jan        
 Min.   : 3.154  
 1st Qu.:11.342  
 Median :14.200  
 Mean   :15.044  
 3rd Qu.:18.929  
 Max.   :27.126  
 NAs    :14715

NA cells are the sea — we masked everything outside Mexico when the file was made. Most raster cells in a real analysis are NA, and every function needs na.rm = TRUE.

Mapping: Base plot

plot(tavg[["Jan"]],
     main = "Mean temperature, January")

terra’s plot is fast and needs no arguments. Use it constantly while you work.

Mapping: ggplot + tidyterra

ggplot() +
  geom_spatraster(data = tavg[["Jan"]]) +
  scale_fill_viridis_c(
    option = "inferno",
    na.value = NA,
    name = "°C") +
  theme_void()

geom_spatraster() is the raster equivalent of geom_sf(). Everything you know about scales and themes still applies.

All Twelve Months

ggplot() +
  geom_spatraster(data = tavg) +
  facet_wrap(~lyr, ncol = 6) +
  scale_fill_viridis_c(option = "inferno", na.value = NA, name = "°C") +
  theme_void(base_size = 7)

One facet_wrap over twelve layers. The north swings hard between winter and summer; the Yucatán barely moves.

Part C: Raster Algebra

Arithmetic Just Works

Rasters on the same grid can be added, subtracted, multiplied, compared — cell by cell.

#Step1: The annual mean across the 12 monthly layers
annual <- mean(tavg)
names(annual) <- "tavg"

#Step2: Seasonal amplitude -- how much the year swings
amplitude <- max(tavg) - min(tavg)
names(amplitude) <- "amplitude"

round(c(min = min(values(annual), na.rm = TRUE),
        max = max(values(annual), na.rm = TRUE)), 1)
 min  max 
 8.0 28.2

8.0 °C to 28.2 °C across the country.

Two Derived Surfaces

p1 <- ggplot() + geom_spatraster(data = annual) +
  scale_fill_viridis_c(option = "inferno", na.value = NA) + ggtitle("Annual mean")
p2 <- ggplot() + geom_spatraster(data = amplitude) +
  scale_fill_viridis_c(option = "mako", na.value = NA) + ggtitle("Seasonal swing")
p1 + p2

Two completely different maps from the same twelve layers.

Logical Operations Give You Masks

#Step1: A comparison returns TRUE/FALSE per cell
hot <- annual > 25

#Step2: TRUE counts as 1, so summing counts cells
sum(values(hot), na.rm = TRUE)
[1] 1130
sum(!is.na(values(annual)))
[1] 6789

1,130 of 6,789 land cells average above 25 °C — about 17% of the country.

Warning

Cell counts are not area unless every cell is the same size. In EPSG:4326 they are not — a degree of longitude shrinks toward the poles. Project first if you need real areas.

Coarsening: aggregate

#Step1: Combine blocks of 3x3 cells by taking their mean
coarse <- aggregate(annual, fact = 3, fun = mean)

c(cells_before = ncell(annual), cells_after = ncell(coarse))
cells_before  cells_after 
       21504         2432 
round(c(res_before = res(annual)[1], res_after = res(coarse)[1]), 3)
res_before  res_after 
     0.167      0.500

21,504 cells become 2,432. Use this when a raster is finer than your question — it makes everything downstream faster.

The reverse is disagg(), which invents detail. It does not add information.

Reprojecting a Raster

projected <- project(annual, "EPSG:6372")
round(res(projected))
[1] 17254 17254

Now each cell is about 17 km × 17 km in real metres.

Important

Reprojecting a vector moves exact coordinates. Reprojecting a raster must build a new grid and interpolate values into it — so it changes your data. Reproject once, as late as possible, and use method = "near" for categorical rasters like land cover.

Part D: Rasters Meet Vectors

Crop and Mask

Two operations that sound alike and are not:

  • crop cuts the raster down to a rectangle — the bounding box of your vector
  • mask sets cells outside the actual shape to NA
sta <- st_read("../data/mexico_homicides.gpkg", layer = "states", quiet = TRUE) %>%
  st_transform(4326)
yuc <- sta %>% filter(state_name == "Yucatán")

cropped <- crop(annual, vect(yuc))     # rectangle
masked  <- mask(cropped, vect(yuc))    # shape

You almost always want crop then mask, in that order — cropping first makes the mask far cheaper.

The Difference, Drawn

extract: The Operation That Matters

This is why rasters are worth learning. extract turns a surface into a column in your table.

#Step1: One value per polygon -- the mean of all cells inside it
e <- terra::extract(annual, vect(sta), fun = mean, na.rm = TRUE)
head(e, 3)
  ID     tavg
1  1 16.69675
2  2 18.36829
3  3 21.48210
#Step2: Attach it to the sf object
sta$tavg <- e$tavg

Hottest and Coolest States

sta %>% st_drop_geometry() %>% arrange(desc(tavg)) %>%
  transmute(state_name, tavg = round(tavg, 1)) %>% head(4)
    state_name tavg
1      Tabasco 26.5
2     Campeche 26.0
3      Yucatán 25.8
4 Quintana Roo 25.7
sta %>% st_drop_geometry() %>% arrange(tavg) %>%
  transmute(state_name, tavg = round(tavg, 1)) %>% head(4)
        state_name tavg
1         Tlaxcala 13.4
2 Ciudad de México 14.4
3           México 15.4
4        Chihuahua 16.5

Tabasco at 26.5 °C; Tlaxcala at 13.4 °C. Both are at similar latitudes — the difference is altitude.

A Choropleth of Extracted Values

ggplot(sta) +
  geom_sf(aes(fill = tavg),
          colour = "white",
          linewidth = 0.15) +
  scale_fill_viridis_c(option = "inferno",
                       name = "°C") +
  coord_sf(crs = 6372) +
  theme_void()

We are back in Session 3. Once extracted, a raster variable is just another column to map, classify, and join.

Cells Do Not Respect Borders

An 18 km cell straddling a state line belongs partly to each. terra::extract includes a cell if its centre falls inside.

exactextractr weights each cell by the fraction of it that the polygon covers:

sta$tavg_exact <- exact_extract(annual, sta, "mean", progress = FALSE)

round(max(abs(sta$tavg - sta$tavg_exact), na.rm = TRUE), 3)
[1] 0.169

A maximum difference of 0.17 °C here — small, because our states are far larger than the cells. For municipios at this resolution the difference would be serious.

Tip

Rule of thumb: if your polygons are not many times bigger than your cells, use exactextractr.

Two Rasters, One Question

#Step1: Extract elevation the same way
elev <- rast("../data/raster/mexico_elevation.tif")
sta$elev <- terra::extract(elev, vect(sta), fun = mean, na.rm = TRUE)[, 2]

#Step2: How strongly do they move together?
round(cor(sta$tavg, sta$elev), 3)
[1] -0.906
#Step3: How much cooler per 1,000 m?
round(coef(lm(tavg ~ elev, data = st_drop_geometry(sta)))[2] * 1000, 2)
 elev 
-4.14

A correlation of -0.91, and about 4.1 °C cooler per 1,000 metres.

Is That the Right Number?

The physical free-air lapse rate is about 6.5 °C per 1,000 m. We got 4.1. Which is wrong?

Neither. First, check that averaging 32 states did not distort the estimate — refit on all 6,789 raw cells:

xy <- crds(annual, na.rm = FALSE)
cells <- na.omit(data.frame(t = values(annual)[, 1],
                            e = values(elev)[, 1],
                            lat = xy[, 2]))
round(coef(lm(t ~ e, data = cells))[2] * 1000, 2)
    e 
-4.15

-4.15 — essentially identical to the state-level -4.14. Averaging first was harmless here.

Two Different Quantities

The gap is not an error. A surface lapse rate is not a free-air lapse rate.

  • The 6.5 °C/km figure describes air cooling as you rise through the atmosphere at one place
  • Ours describes how surface temperature differs between low and high places — which also depends on how the ground heats the air above it

Surface lapse rates are typically shallower. And it is not latitude doing it:

round(coef(lm(t ~ e + lat, data = cells))["e"] * 1000, 2)
    e 
-4.09

Tip

The habit that matters: when a number disagrees with a textbook, check whether you computed the same quantity before concluding your data is broken. Here, three specifications all give -4.1, which is a strong hint the estimate is fine and the comparison was not.

Part E: From Points to a Raster

Rasterizing Point Counts

Sometimes you want to go the other way: thousands of points, summarised onto a grid.

#Step1: The DENUE health points from Sessions 2 and 4
d   <- read.csv("../data/denue/denue_salud_yucatan.csv")
pts <- vect(st_as_sf(d, coords = c("longitud", "latitud"), crs = 4326))

#Step2: An empty grid to count into -- 0.02 degrees is about 2 km
tmpl <- rast(ext(-90.5, -87.4, 19.9, 21.7), resolution = 0.02, crs = "EPSG:4326")

#Step3: Count the points falling in each cell
counts <- rasterize(pts, tmpl, fun = "count", background = 0)

c(cells = ncell(counts), occupied = sum(values(counts) > 0), max = max(values(counts)))
   cells occupied      max 
   13950      337      564

The Density Grid

ggplot() +
  geom_spatraster(
    data = log1p(counts)) +
  scale_fill_viridis_c(
    option = "magma",
    na.value = NA,
    name = "log(1+n)") +
  theme_void()

Counts run from 0 to 564, so the raw scale shows Mérida and nothing else. log1p — log of one plus the count — compresses that and handles the zeros.

Resolution Changes the Answer

Cell size Grid Occupied cells Busiest cell
~6 km 36 × 62 200 1,478
~2 km 90 × 155 337 564

Same points, same state. The grid you choose decides how concentrated the pattern looks.

Warning

This is the modifiable areal unit problem, and it applies to municipios exactly as much as to grid cells. There is no neutral choice — only a stated one.

When to Rasterize, and When Not To

Do rasterize when you need to combine points with other continuous surfaces, or when your units are arbitrary anyway.

Do not rasterize when a real administrative unit is the thing you care about. A municipio is a budget and a mayor; a grid cell is not.

For the reto: rasterizing is usually a step on the way to an extract, not an output in itself.

Part F: Wrapping Up

What We Have Not Covered

  • Mosaicing — stitching adjacent tiles into one surface
  • Resampling — moving one raster onto another’s grid (near, bilinear, average)
  • Terrain analysis — slope, aspect, hillshade from a DEM
  • Focal operations — moving windows, smoothing, edge detection

All are in terra, and all follow the patterns you have now seen.

Where to Get Rasters

Source What Resolution
WorldClim climate normals 1 km – 20 km
WorldPop population counts 100 m – 1 km
SRTM / Copernicus DEM elevation 30 m
VIIRS night lights economic activity proxy 500 m
INEGI / CONABIO land use, vegetation varies

For the reto, WorldPop is the one to reach for — it gives you the population denominator we kept wishing for in Session 3.

What Have We Learned?

  • A raster is a matrix + extent + CRS; no features, no attribute table, joined only by location
  • Use terra; raster is retired and most online answers are written for it
  • Layers are like columns — twelve months in one file, and mean()/max() work across them
  • crop then mask: rectangle first, shape second
  • extract turns a surface into a column, and then you are back in Session 3
  • Use exactextractr when polygons are not much bigger than cells
  • When a number disagrees with a textbook, check you computed the same quantity — our -4.1 °C/km is a surface lapse rate, not a free-air one
  • Grid resolution changes the pattern: the modifiable areal unit problem

That is the course. You can now read, join, measure, classify, map, and extract — which is most of applied spatial analysis.