Data Science II

Session 5: Rasters — Continuous Surfaces and Zonal Statistics

Bogdan G. Popescu
bgpopescu@tec.mx

Tec de Monterrey

What You’ll Learn Today

  1. What a raster is — and why it holds the same information as a table, stored differently
  2. Reading, inspecting, and mapping rasters with terra and tidyterra
  3. Raster algebra — arithmetic, masks, and changing the resolution
  4. crop, mask and extract — turning a surface into a column in your table

Today’s data: monthly temperature across Mexico, and the GADM state boundaries.

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)

Two small GeoTIFFs for ../data/raster/, plus the GADM state boundaries you already know from Session 1:

Part A: What a Raster Is

Vector and Raster

  • Spatial data comes as a vector or a raster
  • A river as a vector (left) and as a raster (right)

Sessions 1–4 were all vectors. Today is the other half.

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.

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 3, 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

Already cropped to Mexico. The global download is 37 MB; cropped, the two come to 300 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.

The Same Raster, as a Table

#Step1: any raster can become a data frame -- x, y, one column per layer
tavg_df <- as.data.frame(tavg, xy = TRUE)

dim(tavg_df)
[1] 6789   14
head(tavg_df[, 1:5], 3)
            x     y      Jan      Feb      Mar
211 -115.4167 32.75 12.50225 14.92875 17.32875
212 -115.2500 32.75 12.55275 15.18300 17.78650
213 -115.0833 32.75 13.02050 15.62750 18.26500

One row per cell, one column per layer. The 21,504 cells become 6,789 rows because as.data.frame drops the NA sea cells by default.

#Step2: underneath, one layer is just a matrix -- no x, no y anywhere
as.matrix(tavg[["Jan"]], wide = TRUE)[1:3, 17:21]
        [,1]     [,2]     [,3]     [,4]     [,5]
[1,]     NaN      NaN      NaN      NaN      NaN
[2,]     NaN      NaN 12.50225 12.55275 13.02050
[3,] 13.6095 13.44675 13.37775 12.88275 12.77125

The same three January values, but nothing in there says where they are. Row 2, columns 19–21 is the location — extent and resolution turn an index into a coordinate.

As RStudio Sees Them

Same data, two objects. RStudio offers you tavg_df as a spreadsheet — 6,789 obs. of 14 variables. For tavg it can only say Formal class SpatRaster, because there is no table in there to show you.

Same Picture, Either Way

#Step1: from the raster
ggplot() + geom_spatraster(data = tavg[["Jan"]])

#Step2: from the data frame -- ordinary ggplot2, nothing spatial about it
ggplot(tavg_df, aes(x = x, y = y, fill = Jan)) + geom_raster() + coord_equal()

The same map twice. geom_raster() is plain ggplot2 — x, y, fill. Nothing in the right-hand panel knows it is a map.

So Why Have Rasters at All?

The same 6,789 locations, stored three ways:

what it stores size
sf points a geometry per cell 3,536 KB
data frame x and y columns 771 KB
SpatRaster no coordinates at all extent + resolution

A raster’s cells are on a regular grid, so terra derives any cell’s coordinates from the extent and the resolution. Nothing positional is stored per cell — only the values.

Note

That is the whole trade. Vectors store where each thing is, because things sit anywhere. Rasters store only values, because the grid already says where. Everything else about rasters follows from this.

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.

Most of a Raster Is Nothing

A raster is always a rectangle. 21,504 cells, of which 6,789 carry a value — 68% of this file is NA. Mexico is just the part that happens to have data in it.

mean(values(tavg[["Jan"]]))
[1] NaN
mean(values(tavg[["Jan"]]), na.rm = TRUE)
[1] 15.04444

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.

Your Turn: Map July

We mapped January. Map July instead — then put the two months side by side. Which part of the country changes most?

#Step1: pull one layer out by name -- or by position, tavg[[7]]
jul <- tavg[["Jul"]]

#Step2: map it exactly as before
plot(jul, main = "Mean temperature, July")

#Step3: now both months at once
plot(tavg[[c("Jan", "Jul")]], main = c("January", "July"))

Look hard at the legends on that last plot before you conclude anything.

Mapping: ggplot + tidyterra

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

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

The Same Layer as sf Polygons

#Step1: turn every cell into its own square polygon
jan_sf <- st_as_sf(as.polygons(tavg[["Jan"]], dissolve = FALSE, na.rm = TRUE))

ggplot(jan_sf) + geom_sf(aes(fill = Jan), colour = NA)

The same map — now as 6,789 POLYGON features, each storing its own five corner coordinates. The values take 168 KB; the polygons take 5,150 KB, about 30 times more, to say exactly what extent and resolution already said.

All Twelve Months

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

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_temp <- mean(tavg)
names(annual_temp) <- "tavg"

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

#minmax() reads the range straight off the layer
round(minmax(annual_temp), 1)
    tavg
min  8.0
max 28.2

8.0 °C to 28.2 °C across the country.

Two Derived Surfaces

map_annual_temp <- ggplot() + geom_spatraster(data = annual_temp) +
  scale_fill_viridis_c(option = "inferno", na.value = NA) + ggtitle("Annual mean")
map_temp_range <- ggplot() + geom_spatraster(data = temp_range) +
  scale_fill_viridis_c(option = "mako", na.value = NA) + ggtitle("Seasonal swing")
map_annual_temp + map_temp_range

Two completely different maps from the same twelve layers.

Logical Operations Give You Masks

#Step1: A comparison returns TRUE/FALSE per cell
is_hot <- annual_temp > 25

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.

The Mask, Mapped

#Step1: filter() keeps the matching cells, turns the rest into NA
hot_only <- annual_temp %>%
  filter(tavg > 25)              # = mask(annual_temp, is_hot, ...)

ggplot() + geom_spatraster(data = hot_only) +
  scale_fill_viridis_c(option = "inferno", na.value = NA)

filter() masks by values; masking by a shape still needs mask().

Your Turn: Where Does July Pass 32 °C?

Build the mask, count the cells, then map them. Which part of the country survives?

#Step1: the mask
scorching <- tavg[["Jul"]] > 32

#Step2: count the cells before you draw anything
sum(values(scorching), na.rm = TRUE)

#Step3: keep only those cells, then map as before
hot_july <- tavg[["Jul"]] %>%
  filter(Jul > 32)

Guess the number first. There are 6,789 land cells in total.

Coarsening: aggregate

#Step1: where we start -- 0.17 degrees, 21,504 cells
ncell(annual_temp)
[1] 21504
round(res(annual_temp)[1], 3)
[1] 0.167

Coarsening: aggregate

#Step2: blocks of 3x3 cells, each replaced by their mean
coarse <- aggregate(annual_temp, fact = 3, fun = mean)
ncell(coarse)
[1] 2432
round(res(coarse)[1], 3)
[1] 0.5

Coarsening: aggregate

#Step3: push it to 9x9
coarser <- aggregate(annual_temp, fact = 9, fun = mean)
ncell(coarser)
[1] 286
round(res(coarser)[1], 3)
[1] 1.5

Reprojecting a Raster

projected <- project(annual_temp, "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
states <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_1.shp", quiet = TRUE) %>%
  st_transform(4326) %>%
  rename(state_name = NAME_1)
yuc <- states %>%
  filter(state_name == "Yucatán")

cropped <- crop(annual_temp, 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
tavg_by_state <- terra::extract(annual_temp, vect(states), fun = mean, na.rm = TRUE)
head(tavg_by_state, 3)
  ID     tavg
1  1 16.69675
2  2 18.36143
3  3 21.49858
#Step2: Attach it to the sf object
states$tavg <- tavg_by_state$tavg

What extract Did

Every cell inside a boundary collapses to its mean. On the same colour scale you can see the cost: the surface ran 8.0 to 28.2 °C, the 32 state means only 13.4 to 26.5 °C. Averaging eats both tails.

Hottest States

states %>%
  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

Tabasco at 26.5 °C — and the next three are its neighbours. The whole southeastern lowland, none of it far above sea level.

All Thirty-Two, Ranked

One extract call, and every state has a climate variable. The ordering is altitude, not latitude — the four coolest are all on the central plateau.

A Choropleth of Extracted Values

ggplot(states) +
  geom_sf(aes(fill = tavg),
          colour = "white",
          linewidth = 0.15) +
  scale_fill_viridis_c(option = "inferno",
                       name = "°C")

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:

states$tavg_exact <- exact_extract(annual_temp, states, "mean", progress = FALSE)

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

A maximum difference of 0.165 °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.

The Two Rules, Drawn

Tlaxcala, the smallest state, against 18 km cells. terra::extract keeps the 13 cells whose centre falls inside and ignores the rest. exactextractr uses all 26 cells that touch it, each weighted by its share — and only 3 of those lie entirely within the state.

The Same Polygon, Two Answers

terra::extract(annual_temp, vect(tlax), fun = mean, na.rm = TRUE)$tavg
exact_extract(annual_temp, tlax, "mean")

Same polygon, same raster, two defensible numbers 0.12 °C apart. Nothing is broken — the two functions were asked slightly different questions, and the smaller the polygon gets relative to the cell, the further apart their answers drift.

Part E: Wrapping Up

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. The same information as a table of x, y and values — with the coordinates left out, because the grid already says where
  • 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
  • Most cells are NA, so every summary needs na.rm = TRUE
  • A comparison gives you a mask, and filter() applies it — by value, never by shape
  • Coarsening is cheap and lossy — by fact = 9 the Sierra Madre is gone, but the map still looks like data
  • crop then mask: rectangle first, shape second
  • extract turns a surface into a column and you are back in Session 3 — with exactextractr when the polygons are not much bigger than the cells

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

Teaching Evaluation

forms.gle/3VzPTovn8FxbWSUG9