Spatial Analysis with sf

1. Why Spatial Data Are Different

Every dataset you have worked with so far treats observations as independent — the value at row 27 tells you nothing about the value at row 28. Spatial data violates this assumption almost by definition. Tobler’s First Law of Geography states it plainly:

“Everything is related to everything else, but near things are more related than distant things.” — Waldo Tobler, 1970

A temperature measurement in downtown Los Angeles is a better predictor of the temperature in a neighbouring parking lot than of the temperature in a forest 50 km away. This spatial autocorrelation is not a nuisance to eliminate — it is information to exploit. Kriging and spatial regression models do exactly that.

Spatial data also carry a property no other data type has: location on a curved surface. The Earth is not flat (apologies to anyone who thought this was a statistics tutorial), and representing curved-Earth coordinates on a flat screen or page requires a projection — a mathematical transformation that always introduces some distortion. Managing projections correctly is the single most common source of errors in spatial analysis.

The three things that make spatial analysis special

1. Spatial autocorrelation. Values at nearby locations tend to be similar. Standard statistical tests assume independence; applying them naively to spatial data produces false precision (your effective sample size is smaller than it looks).

2. The modifiable areal unit problem (MAUP). When you aggregate point data into areas (census tracts, counties, watersheds), the results depend on how you draw the boundaries. The same underlying data can tell opposite stories depending on aggregation choices — a sobering fact for any policy analysis.

3. Edge effects. Locations near the boundary of your study area have fewer neighbours than interior locations. Spatial statistics that depend on neighbourhood relationships behave differently at edges, and ignoring this biases results.

📌 Spatial data types: vector vs. raster

Environmental spatial data come in two fundamental flavours. Vector data represent discrete objects: a weather station is a point, a river is a line, a watershed is a polygon. Each object carries attribute columns (like a spreadsheet row) alongside its geometry. The sf package handles vector data. Raster data divide space into a regular grid of cells, each storing one value (elevation, NDVI, land surface temperature). The terra and stars packages handle rasters. This module focuses on vector data; raster workflows appear in the Remote Sensing module.

2. Core Spatial Concepts

2.1 Geometry types

The sf package (Simple Features) represents spatial objects using five core geometry types. Every row in an sf data frame has one geometry:

Type	Description	Environmental example
`POINT`	Single (x, y) coordinate	Weather station, sampling site
`MULTIPOINT`	Set of points	Monitoring network
`LINESTRING`	Sequence of connected points	River reach, road
`POLYGON`	Closed ring of points	Watershed, city boundary, lake
`MULTIPOLYGON`	Collection of polygons	US state with islands

2.2 The sf data frame

An sf object is a regular R data frame with one special column — conventionally named geometry — that stores the spatial information. Every other column is an ordinary attribute you can wrangle with dplyr verbs.

# A typical sf object printed in the console
Simple feature collection with 10 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -118.2 ymin: 29.7 xmax: -71.1 ymax: 47.6
CRS:           EPSG:4326
# A tibble: 10 × 5
  city          temp_C landuse   pop_million geometry
  <chr>          <dbl> <chr>           <dbl> <POINT [°]>
1 Los Angeles     32.1 Urban             3.9 (-118.24 34.05)
2 Phoenix         41.3 Urban             1.6 (-112.07 33.45)
...

The geometry column is “sticky” — it travels with the data frame through filter(), mutate(), left_join(), and most dplyr operations automatically.

2.3 Spatial relationships

sf provides predicate functions that test geometric relationships:

Function	Tests whether…	Use case
`st_intersects()`	geometries share any point	Does this watershed touch this county?
`st_within()`	A is fully inside B	Which stations are within city limits?
`st_contains()`	B fully contains A	Which counties contain a monitoring site?
`st_distance()`	Computes pairwise distances	How far is each station from the coast?
`st_buffer()`	Creates a buffer zone	1 km buffer around each industrial facility

3. Shapefiles: What They Are and Why They Come in Folders

The shapefile is the most common vector data format in environmental science, developed by Esri in the early 1990s and now effectively a universal standard. Every EPA monitoring dataset, USGS watershed boundary, and Census TIGER file you will download comes as a shapefile — and every shapefile is actually four or more files that must stay together.

Geometry

.shp
The actual shapes — coordinates of every point, line, or polygon vertex. The file everyone calls “the shapefile.”

Index

.shx
A positional index that lets software jump directly to a specific feature’s geometry without reading the whole .shp file.

Attributes

.dbf
A dBASE IV table storing all non-spatial columns. This is your spreadsheet — names, values, categories — one row per feature.

Projection

.prj
A plain-text file describing the coordinate reference system. Without this, software cannot place your features on Earth correctly.

Spatial index

.sbn / .sbx
Optional Esri spatial index files that speed up spatial queries. Not required by sf.

Encoding

.cpg
Specifies the character encoding for the .dbf attribute names (e.g., UTF-8). Prevents garbled accented characters.

⚠️ The shapefile folder rule

All files with the same base name must be in the same directory. If a colleague emails you watersheds.shp alone, you cannot read it — you need watersheds.shp, watersheds.shx, watersheds.dbf, and watersheds.prj at minimum. When sharing shapefiles, always zip the entire folder. When reading with sf::st_read(), point to the .shp file and sf loads the rest automatically.

Reading and writing shapefiles with sf

library(sf)

# Reading — point to the .shp file; sf loads .dbf, .prj, .shx automatically
watersheds <- st_read("data/shapefiles/watersheds.shp")

# sf also reads GeoJSON (single file, easier to share):
stations   <- st_read("data/stations.geojson")

# Writing back out
st_write(watersheds, "outputs/watersheds_clean.shp")       # shapefile
st_write(watersheds, "outputs/watersheds_clean.geojson")   # GeoJSON
st_write(watersheds, "outputs/watersheds_clean.gpkg")      # GeoPackage (preferred modern format)

💡 Prefer GeoPackage over shapefile for new work

The GeoPackage (.gpkg) format is a modern single-file alternative to the shapefile that supports longer column names (shapefile truncates to 10 characters!), UTF-8 natively, multiple layers per file, and files larger than 2 GB. sf reads and writes GeoPackage with the same st_read() / st_write() calls.

4. Coordinate Reference Systems

A coordinate reference system (CRS) is the mathematical framework that maps (x, y) numbers to locations on Earth. Getting CRS management wrong is the most common way spatial analyses go silently wrong — two layers that look identical in print but use different CRS will appear in completely different parts of the globe when overlaid.

4.1 Geographic vs. projected CRS

Geographic CRS (also called unprojected or lat/lon) stores locations as longitude and latitude in decimal degrees. The most common is WGS 84 (EPSG:4326), used by GPS devices worldwide and by most downloaded datasets. Distances in degrees are meaningless for area calculations — a degree of longitude near the equator covers ~111 km, but near the poles it covers nearly zero.

Projected CRS transforms the curved Earth onto a flat plane, with coordinates in metres (or feet). Now distances and areas are meaningful. The choice of projection always involves a trade-off:

Projection family	Preserves	Distorts	Good for
Conformal (e.g., UTM, State Plane)	Angles / shape	Area	Navigation, local mapping
Equal-area (e.g., Albers, Mollweide)	Area	Shape	National statistics, choropleth maps
Equidistant	Distance from centre	Area and shape	Distance calculations

4.2 EPSG codes: the universal shorthand

Every standard CRS has an EPSG code — a four- or five-digit number that uniquely identifies it. This is the most reliable way to specify a CRS in code.

EPSG	Name	Common use
4326	WGS 84 (lat/lon)	GPS, most downloaded data, web maps
3857	Web Mercator	Google Maps, OpenStreetMap tile backgrounds
5070	NAD83 / Conus Albers	Continental US equal-area analysis
32610–32619	UTM Zone 10N–19N	Accurate local distance/area in western US
26910–26919	NAD83 UTM Zone 10N–19N	USGS topographic data

4.3 CRS operations in sf

library(sf)

# Check the CRS of an object
st_crs(my_sf)

# Transform to a different CRS
# Always transform, never assign unless you are fixing a wrongly-labelled CRS
my_sf_albers <- st_transform(my_sf, crs = 5070)   # NAD83 Conus Albers
my_sf_utm    <- st_transform(my_sf, crs = 32611)   # UTM Zone 11N

# WRONG — do not do this unless the CRS label is simply incorrect:
st_crs(my_sf) <- 4326   # only to fix a mis-labelled CRS, not to reproject

⚠️ st_transform() vs. assigning st_crs() — a critical distinction

st_transform(sf, crs) reprojects the coordinates — it moves the points on Earth to match the new system. st_crs(sf) <- crs only relabels the CRS without moving anything. If you use the assignment form on correctly-labelled data, you will silently place your features in the wrong location on the globe. Always use st_transform() for reprojection.

5. Data Setup: Urban Heat Islands Across 10 US Cities

The urban heat island (UHI) effect describes how cities are systematically warmer than surrounding rural areas. Impervious surfaces (concrete, asphalt) absorb heat; sparse vegetation provides less evaporative cooling; waste heat from vehicles and air conditioners adds further energy. The effect is strongest on calm, clear summer nights and can exceed 10°C in large cities.

We will simulate a realistic dataset of 10 US cities with contrasting land-use compositions, multiple monitoring stations per city, and summer afternoon temperatures. This dataset is constructed to be analytically identical to what you would obtain from EPA AQS or the National Land Cover Database (NLCD).

library(sf)         # vector spatial data
library(dplyr)
library(tidyr)
library(ggplot2)
library(viridis)    # perceptually uniform colour scales
library(knitr)
library(kableExtra)
# gstat and sp are used for kriging (Section 9)
# install.packages(c("gstat", "sp")) if needed

set.seed(2024)

# ── City metadata ────────────────────────────────────────────────────────────
cities <- data.frame(
  city        = c("Los Angeles", "Phoenix", "Houston", "Chicago",
                  "New York",    "Seattle", "Denver",  "Atlanta",
                  "Minneapolis", "Portland"),
  lon         = c(-118.24, -112.07, -95.37, -87.63,
                  -74.01,  -122.33,  -104.99, -84.39,
                  -93.27,  -122.68),
  lat         = c(34.05,  33.45,  29.76,  41.88,
                  40.71,  47.61,   39.74,  33.75,
                  44.98,  45.52),
  # Fraction of impervious surface in urban core (0–1)
  imperv_frac = c(0.72, 0.78, 0.65, 0.61,
                  0.80, 0.52, 0.55, 0.58,
                  0.60, 0.49),
  # Fraction of tree canopy cover in urban core (0–1)
  canopy_frac = c(0.18, 0.08, 0.22, 0.19,
                  0.21, 0.38, 0.25, 0.35,
                  0.28, 0.42),
  # Climate zone affects baseline temperature
  climate     = c("Mediterranean", "Desert", "Humid subtropical", "Continental",
                  "Continental",   "Oceanic", "Semi-arid", "Humid subtropical",
                  "Continental",  "Oceanic"),
  # Mean July afternoon temperature (°C) at a reference rural station nearby
  rural_temp_C = c(29.8, 40.2, 33.1, 28.4,
                   28.9, 22.6, 31.4, 32.7,
                   27.8, 23.5)
)

# ── Simulate monitoring stations within each city ────────────────────────────
# Each city has 6–10 stations with varying land-use fractions
n_stations_per_city <- sample(6:10, nrow(cities), replace = TRUE)

stations <- do.call(rbind, lapply(seq_len(nrow(cities)), function(i) {
  n  <- n_stations_per_city[i]
  ct <- cities[i, ]

  # Randomly vary impervious fraction across stations within the city
  imp  <- pmin(pmax(ct$imperv_frac + rnorm(n, 0, 0.12), 0.05), 0.98)
  can  <- pmin(pmax(ct$canopy_frac + rnorm(n, 0, 0.08), 0.01), 0.70)

  # Temperature model:
  # base = rural temp + UHI from impervious - cooling from canopy + noise
  uhi_effect <- 6.5 * imp - 3.2 * can
  temp_C     <- ct$rural_temp_C + uhi_effect + rnorm(n, 0, 0.8)

  # Scatter stations within ~0.15 degrees of city centre
  data.frame(
    city        = ct$city,
    station_id  = paste0(substr(ct$city, 1, 3), "_", sprintf("%02d", seq_len(n))),
    lon         = ct$lon + runif(n, -0.15, 0.15),
    lat         = ct$lat + runif(n, -0.12, 0.12),
    imperv_frac = round(imp, 3),
    canopy_frac = round(can, 3),
    temp_C      = round(temp_C, 2),
    rural_temp_C = ct$rural_temp_C,
    uhi_C       = round(temp_C - ct$rural_temp_C, 2),
    climate     = ct$climate
  )
}))

cat("Total monitoring stations:", nrow(stations), "\n")

## Total monitoring stations: 82

cat("Cities:", nrow(cities), "\n")

## Cities: 10

stations |>
  group_by(city) |>
  summarise(
    n           = n(),
    mean_temp   = round(mean(temp_C), 1),
    mean_uhi    = round(mean(uhi_C),  1),
    mean_imperv = round(mean(imperv_frac) * 100, 0),
    mean_canopy = round(mean(canopy_frac) * 100, 0)
  ) |>
  arrange(desc(mean_temp)) |>
  kable(col.names = c("City", "Stations", "Mean temp (°C)",
                       "Mean UHI (°C)", "% Impervious", "% Canopy"),
        caption = "Summary of simulated urban temperature monitoring data across 10 US cities.") |>
  kable_styling(bootstrap_options = "striped")

Summary of simulated urban temperature monitoring data across 10 US cities.
City	Stations	Mean temp (°C)	Mean UHI (°C)	% Impervious	% Canopy
Phoenix	10	44.1	3.9	70	9
Houston	10	35.4	2.3	60	26
Atlanta	7	35.0	2.3	59	40
Denver	6	34.7	3.3	64	25
Los Angeles	7	34.0	4.2	73	21
New York	6	33.9	5.0	80	22
Chicago	9	31.6	3.2	57	16
Minneapolis	7	30.2	2.4	58	31
Portland	10	25.4	1.9	52	45
Seattle	10	25.2	2.6	54	36

6. sf Basics: Building and Exploring Spatial Objects

6.1 Creating an sf object from a data frame

Any data frame with longitude and latitude columns can be converted to an sf point object with st_as_sf(). You must specify which columns hold the coordinates and what CRS those coordinates use.

library(sf)

# Convert monitoring stations to sf POINT object
# coords = c("lon", "lat") tells sf which columns are x and y
# crs = 4326 means WGS 84 geographic coordinates (decimal degrees)
stations_sf <- st_as_sf(stations,
                         coords = c("lon", "lat"),
                         crs    = 4326)

# The geometry column is now attached
print(stations_sf[1:4, c("station_id", "city", "temp_C", "geometry")])

## Simple feature collection with 4 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -118.384 ymin: 34.0036 xmax: -118.1752 ymax: 34.12556
## Geodetic CRS:  WGS 84
##   station_id        city temp_C                   geometry
## 1     Los_01 Los Angeles  34.78 POINT (-118.1752 34.01939)
## 2     Los_02 Los Angeles  34.99  POINT (-118.384 34.08346)
## 3     Los_03 Los Angeles  34.80  POINT (-118.3324 34.0036)
## 4     Los_04 Los Angeles  32.52 POINT (-118.2531 34.12556)

# Check geometry type and CRS
st_geometry_type(stations_sf)[1]

## [1] POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

st_crs(stations_sf)$epsg

## [1] 4326

# Bounding box of the dataset
st_bbox(stations_sf)

##       xmin       ymin       xmax       ymax 
## -122.82237   29.66351  -73.86421   47.72826

6.2 Plotting sf objects: a first map

sf provides a simple plot() method that maps each attribute column:

# Quick exploratory plot — one panel per attribute column
plot(stations_sf[, c("temp_C", "imperv_frac", "canopy_frac", "uhi_C")],
     pch = 16, cex = 0.9)

Quick sf plot of monitoring station temperatures. Each panel maps one attribute variable. The default colour scale is applied automatically. For publication-quality maps we will use ggplot2 in Section 7.

6.3 Summarising spatial data: city centroids

We can compute city-level summaries and build a second sf object representing city centroids — which we will use for choropleth mapping.

# City-level aggregates
city_summary <- stations |>
  group_by(city) |>
  summarise(
    n_stations  = n(),
    mean_temp   = mean(temp_C),
    mean_uhi    = mean(uhi_C),
    mean_imperv = mean(imperv_frac),
    mean_canopy = mean(canopy_frac),
    rural_temp  = first(rural_temp_C),
    climate     = first(climate),
    # Keep city centre coordinates from the original cities data frame
    lon         = mean(lon),
    lat         = mean(lat)
  ) |>
  left_join(cities[, c("city","lon","lat")], by = "city", suffix = c("_mean","")) |>
  select(-lon_mean, -lat_mean)

city_sf <- st_as_sf(city_summary, coords = c("lon", "lat"), crs = 4326)

7. Choropleth Maps with ggplot2 and sf

A choropleth map fills geographic areas with colour to represent a variable. ggplot2 integrates with sf through geom_sf(), which reads the geometry column automatically and handles all projection details.

7.1 Station-level temperature map

# US state boundaries for context — built into the maps package
us_states <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
us_states <- st_transform(us_states, crs = 4326)

ggplot() +
  geom_sf(data = us_states, fill = "#f0ede6", colour = "#ccc", linewidth = 0.3) +
  geom_sf(data = stations_sf,
          aes(colour = temp_C, size = pmax(uhi_C, 0)),
          alpha = 0.80, shape = 16) +
  scale_colour_viridis_c(option = "magma", name = "Temp (°C)",
                          direction = -1) +
  scale_size_continuous(name = "UHI (°C)", range = c(1.5, 5)) +
  coord_sf(xlim = c(-130, -65), ylim = c(24, 50), expand = FALSE) +
  labs(
    title    = "Urban monitoring station temperatures — Summer afternoon",
    subtitle = "Point size = urban heat island intensity vs. nearby rural reference",
    caption  = "Simulated data; station layout modelled on EPA AQS network"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                  colour = "#1e3a2f", size = 13),
    plot.subtitle = element_text(colour = "#7a7a72", size = 9),
    legend.position = "right",
    panel.grid = element_line(colour = "#eee", linewidth = 0.3)
  )

Summer afternoon temperatures at all monitoring stations across 10 US cities. Point size scales with the urban heat island intensity (departure from rural reference temperature). Hotter cities (Phoenix, Houston) dominate the right end of the colour scale. Seattle and Portland show the weakest heat islands, consistent with their high canopy cover and mild oceanic climate.

7.2 Choropleth of mean UHI by city

A choropleth works best when the geographic units are polygons. Here we create circular buffer polygons around each city centre to give the map more visual weight. In a real analysis you would use actual city boundary polygons downloaded from the Census Bureau.

# Reproject to Albers equal-area for buffering (distances in metres)
city_albers <- st_transform(city_sf, crs = 5070)

# 40 km radius buffer around each city centre
city_circles <- st_buffer(city_albers, dist = 40000)

# Back to geographic CRS for plotting
city_circles <- st_transform(city_circles, crs = 4326)

ggplot() +
  geom_sf(data = us_states, fill = "#f5f2ee", colour = "#ccc", linewidth = 0.3) +
  geom_sf(data = city_circles,
          aes(fill = mean_uhi),
          colour = "white", linewidth = 0.6, alpha = 0.88) +
  geom_sf_text(data = city_sf, aes(label = city),
               size = 2.8, colour = "#1e3a2f", fontface = "bold",
               nudge_y = 1.2) +
  scale_fill_viridis_c(option = "inferno", name = "Mean UHI (°C)",
                        direction = -1, begin = 0.15) +
  coord_sf(xlim = c(-130, -65), ylim = c(24, 50), expand = FALSE) +
  labs(
    title   = "Urban Heat Island Intensity Across 10 US Cities",
    subtitle = "Mean departure from rural reference temperature (summer afternoon)",
    caption  = "Circles = 40 km radius buffers around city centres"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                  colour = "#1e3a2f", size = 13),
    plot.subtitle = element_text(colour = "#7a7a72", size = 9),
    legend.position = "right"
  )

$Mean urban heat island intensity by city (departure from nearby rural temperature). Cities are represented by 40 km radius circles; real analyses would use Census-derived city boundaries. Phoenix leads the heat island rankings, driven by extreme impervious surface cover and minimal canopy. Seattle and Portland show the weakest UHI, consistent with high canopy fractions and oceanic climate.$

Mean urban heat island intensity by city (departure from nearby rural temperature). Cities are represented by 40 km radius circles; real analyses would use Census-derived city boundaries. Phoenix leads the heat island rankings, driven by extreme impervious surface cover and minimal canopy. Seattle and Portland show the weakest UHI, consistent with high canopy fractions and oceanic climate.

8. Spatial Joins: Linking Layers by Location

A spatial join merges two sf objects based on geographic relationship rather than a shared column. The canonical use: you have monitoring stations (points) and you want to attach attributes from the watershed or county polygon they fall within.

8.1 How st_join() works

# Syntax mirrors dplyr::left_join() but uses geometry instead of a key column
result <- st_join(points_sf, polygons_sf, join = st_within)
# result has all rows of points_sf, with polygon attributes appended
# where each point falls within a polygon

8.2 Joining stations to climate zones

We simulate a simple climate zone polygon layer and join it to our stations, then test whether UHI intensity differs by climate.

# Create simple rectangular climate zone polygons for illustration
# In practice you would load a real climate classification shapefile (e.g. Koppen-Geiger)
climate_boxes <- data.frame(
  zone   = c("Hot arid", "Subtropical", "Continental", "Oceanic/Mediterranean"),
  xmin   = c(-115, -100, -100, -130),
  xmax   = c( -100, -80,  -70, -100),
  ymin   = c(  25,   25,   38,   30),
  ymax   = c(  40,   38,   50,   50)
)

climate_sf <- do.call(rbind, lapply(seq_len(nrow(climate_boxes)), function(i) {
  r <- climate_boxes[i, ]
  geom <- st_polygon(list(matrix(
    c(r$xmin, r$ymin,
      r$xmax, r$ymin,
      r$xmax, r$ymax,
      r$xmin, r$ymax,
      r$xmin, r$ymin),
    ncol = 2, byrow = TRUE
  )))
  st_sf(zone = r$zone, geometry = st_sfc(geom, crs = 4326))
}))

# Join each station to its climate zone polygon
stations_joined <- st_join(stations_sf, climate_sf, join = st_within)

# Summarise UHI by climate zone
uhi_by_zone <- stations_joined |>
  st_drop_geometry() |>          # remove geometry for plain dplyr operations
  group_by(zone) |>
  summarise(
    n         = n(),
    mean_uhi  = round(mean(uhi_C, na.rm = TRUE), 2),
    sd_uhi    = round(sd(uhi_C,   na.rm = TRUE), 2),
    mean_imperv = round(mean(imperv_frac) * 100, 1)
  ) |>
  arrange(desc(mean_uhi))

kable(uhi_by_zone,
      col.names = c("Climate zone", "Stations", "Mean UHI (°C)", "SD", "% Impervious"),
      caption   = "Urban heat island intensity by climate zone after spatial join.")

Urban heat island intensity by climate zone after spatial join.
Climate zone	Stations	Mean UHI (°C)	SD	% Impervious
Hot arid	16	3.68	1.20	67.7
Continental	22	3.44	1.73	63.9
Oceanic/Mediterranean	43	3.11	1.33	61.7
Subtropical	17	2.31	0.83	59.9

📌 st_drop_geometry()

Once you have finished spatial operations, call st_drop_geometry() to convert an sf object back to a plain data frame. This speeds up subsequent non-spatial operations and prevents the geometry column from appearing unexpectedly in output tables.

8.3 Visualising the spatial join result

stations_joined |>
  st_drop_geometry() |>
  filter(!is.na(zone)) |>
  mutate(zone = factor(zone, levels = c("Hot arid","Subtropical",
                                         "Continental","Oceanic/Mediterranean"))) |>
  ggplot(aes(x = zone, y = uhi_C, fill = zone)) +
  geom_violin(trim = FALSE, alpha = 0.55, linewidth = 0.5) +
  geom_boxplot(width = 0.10, fill = "white", alpha = 0.85, outlier.size = 1) +
  geom_jitter(width = 0.08, alpha = 0.35, size = 1.2, colour = "#1e3a2f") +
  scale_fill_manual(values = c("#c8832a","#c62828","#2d5a3d","#4a7c5e"), guide = "none") +
  labs(x = NULL, y = "Urban heat island (°C)",
       title = "UHI intensity by climate zone",
       subtitle = "Violin + boxplot; points = individual monitoring stations") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))

Urban heat island intensity by climate zone. Hot arid cities (Phoenix) show the strongest UHI; oceanic/Mediterranean cities (Seattle, Portland, Los Angeles) show weaker effects partly because lower baseline humidity and higher canopy cover moderate surface temperatures. Points are jittered to show all individual stations.

9. Land Use and the Urban Heat Island

The strongest predictor of local temperature is the fraction of impervious surface — concrete, asphalt, rooftops — that absorbs and re-radiates solar energy. Tree canopy provides the counterweight through evapotranspiration. Here we model that relationship explicitly.

9.1 Temperature as a function of land use

p1 <- ggplot(stations, aes(x = imperv_frac * 100, y = temp_C, colour = city)) +
  geom_point(alpha = 0.65, size = 1.8) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE,
              colour = "#1e3a2f", linewidth = 1.1) +
  scale_colour_viridis_d(option = "turbo", guide = "none") +
  labs(x = "Impervious surface (%)", y = "Temperature (°C)",
       title = "Impervious surface") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(family = "serif", face = "bold",
                                   colour = "#1e3a2f"))

p2 <- ggplot(stations, aes(x = canopy_frac * 100, y = temp_C, colour = city)) +
  geom_point(alpha = 0.65, size = 1.8) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE,
              colour = "#1e3a2f", linewidth = 1.1) +
  scale_colour_viridis_d(option = "turbo", name = "City") +
  labs(x = "Canopy cover (%)", y = NULL,
       title = "Tree canopy") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(family = "serif", face = "bold",
                                   colour = "#1e3a2f"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(0.4, "cm"))

library(patchwork)
p1 + p2 + plot_layout(widths = c(1, 1.4))

$Urban station temperature as a function of impervious surface fraction (left) and tree canopy fraction (right), coloured by city. The positive slope for impervious cover and negative slope for canopy are the statistical signature of the urban heat island mechanism. Note that Phoenix dominates the upper range — its desert baseline temperature compounds the impervious surface effect.$

Urban station temperature as a function of impervious surface fraction (left) and tree canopy fraction (right), coloured by city. The positive slope for impervious cover and negative slope for canopy are the statistical signature of the urban heat island mechanism. Note that Phoenix dominates the upper range — its desert baseline temperature compounds the impervious surface effect.

9.2 Multiple regression: controlling for climate

Temperature at a station reflects both its local land cover and its regional climate. A multiple regression lets us estimate the land-use effect after controlling for baseline rural temperature.

# Multiple regression: temp ~ rural_baseline + impervious + canopy
lm_uhi <- lm(temp_C ~ rural_temp_C + imperv_frac + canopy_frac, data = stations)
summary(lm_uhi)

## 
## Call:
## lm(formula = temp_C ~ rural_temp_C + imperv_frac + canopy_frac, 
##     data = stations)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63520 -0.67621 -0.09279  0.66330  1.77604 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.30309    0.91749   2.510   0.0141 *  
## rural_temp_C  0.93757    0.02287  41.002  < 2e-16 ***
## imperv_frac   6.37569    0.78282   8.145 4.94e-12 ***
## canopy_frac  -5.00141    0.86765  -5.764 1.56e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8707 on 78 degrees of freedom
## Multiple R-squared:  0.9781, Adjusted R-squared:  0.9772 
## F-statistic:  1159 on 3 and 78 DF,  p-value: < 2.2e-16

broom::tidy(lm_uhi, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(across(where(is.numeric), round, 3)) |>
  kable(col.names = c("Predictor", "Estimate", "SE", "t", "p", "CI low", "CI high"),
        caption   = "Multiple regression: station temperature as a function of rural baseline and land use. A 10-percentage-point increase in impervious surface adds approximately 0.65°C; the same increase in canopy cover removes approximately 0.32°C.") |>
  kable_styling(bootstrap_options = "striped")

Multiple regression: station temperature as a function of rural baseline and land use. A 10-percentage-point increase in impervious surface adds approximately 0.65°C; the same increase in canopy cover removes approximately 0.32°C.
Predictor	Estimate	SE	t	CI low	CI high
rural_temp_C	0.938	0.023	41.002	0.892	0.983
imperv_frac	6.376	0.783	8.145	4.817	7.934
canopy_frac	-5.001	0.868	-5.764	-6.729	-3.274

10. Kriging: Interpolating Temperature Across Space

Monitoring stations tell us the temperature at specific points. But managers and policymakers often need continuous surfaces — temperature everywhere in a city, not just where we happened to put sensors. Kriging is the gold-standard geostatistical interpolation method for this task.

The key insight: stations that are close together tend to have similar temperatures (spatial autocorrelation). Kriging fits a model of that autocorrelation structure — the variogram — and uses it to make optimal predictions at unmeasured locations, along with prediction uncertainty (the kriging variance).

10.1 What is a variogram?

The variogram γ(h) quantifies how dissimilar pairs of observations are as a function of the distance h between them:

\gamma(h) = \frac{1}{2} \text{Var}[Z(s) - Z(s + h)]

At short distances, nearby observations are similar so γ(h) is small. As h increases, observations become less correlated and γ(h) rises, eventually levelling off at the sill (total variance). The range is the distance beyond which observations are effectively uncorrelated. The nugget represents measurement error or fine-scale variation below the sampling resolution.

10.2 Fitting and plotting the variogram

library(gstat)
library(sp)

# Work with a single city for a clean kriging example
la_stations <- stations |>
  filter(city == "Los Angeles") |>
  select(lon, lat, temp_C, imperv_frac, canopy_frac)

# Convert to sp SpatialPointsDataFrame (gstat uses sp internally)
coordinates(la_stations) <- ~ lon + lat
proj4string(la_stations)  <- CRS("+init=epsg:4326")

# Reproject to UTM zone 11N (metres) for meaningful distance units
la_utm <- spTransform(la_stations, CRS("+init=epsg:32611"))

# Compute empirical variogram
vario_emp <- variogram(temp_C ~ 1, data = la_utm, cutoff = 25000, width = 3000)

# Fit a spherical variogram model
vario_fit <- fit.variogram(vario_emp,
                            model = vgm(psill = 1.5, model = "Sph",
                                        range = 15000, nugget = 0.3))
print(vario_fit)

##   model    psill    range
## 1   Nug 0.000000     0.00
## 2   Sph 1.147199 14815.41

# Plot
plot(vario_emp, vario_fit,
     main = "Variogram — Los Angeles station temperatures",
     xlab = "Distance (m)", ylab = "Semivariance (°C²)")

Empirical variogram (dots) and fitted spherical model (line) for station temperatures within Los Angeles. The y-axis (semivariance) quantifies how different temperature pairs are as a function of separation distance. The curve rises from a small nugget effect, reaches a sill around 1.8 km separation, and flattens — indicating that stations more than about 18 km apart have effectively uncorrelated temperatures after removing spatial trend.

10.3 Building a prediction grid

# Create a regular prediction grid over the LA bounding box
bb <- bbox(la_utm)
grid_pts <- expand.grid(
  lon = seq(bb[1,1], bb[1,2], length.out = 60),
  lat = seq(bb[2,1], bb[2,2], length.out = 60)
)

# Add covariate columns BEFORE converting to sp — assigning to a
# SpatialPixelsDataFrame after gridded() triggers a row-count mismatch error
grid_pts$imperv_frac <- mean(la_utm$imperv_frac)

coordinates(grid_pts) <- ~ lon + lat
proj4string(grid_pts) <- CRS("+init=epsg:32611")
gridded(grid_pts)      <- TRUE   # mark as regular grid for gstat

cat("Prediction grid:", prod(dim(grid_pts)), "cells\n")

## Prediction grid: 3600 cells

10.4 Running ordinary kriging

# Ordinary kriging: interpolate temp_C across the grid
krige_result <- krige(temp_C ~ 1,
                       locations = la_utm,
                       newdata   = grid_pts,
                       model     = vario_fit)

## [using ordinary kriging]

# Convert to data frame for ggplot
krige_df <- as.data.frame(krige_result) |>
  rename(temp_pred = var1.pred, temp_var = var1.var) |>
  mutate(temp_se = sqrt(temp_var))

head(krige_df[, c("lon","lat","temp_pred","temp_se")])

##        lon     lat temp_pred   temp_se
## 1 372308.4 3763332  34.72681 0.9003814
## 2 372633.6 3763332  34.74183 0.8768299
## 3 372958.7 3763332  34.75528 0.8515661
## 4 373283.8 3763332  34.76715 0.8243695
## 5 373608.9 3763332  34.77746 0.7949737
## 6 373934.1 3763332  34.78621 0.7630516

# Extract station coordinates from sp object into a plain data frame
# as.data.frame() on a SpatialPoints* stores coords in a sub-matrix;
# use coordinates() to pull them out as named columns instead.
la_pts_df <- as.data.frame(coordinates(la_utm))
names(la_pts_df) <- c("lon", "lat")

# ── Predicted temperature ─────────────────────────────────────────────────
p_pred <- ggplot(krige_df, aes(x = lon, y = lat, fill = temp_pred)) +
  geom_raster(interpolate = TRUE) +
  geom_point(data = la_pts_df,
             aes(x = lon, y = lat),
             inherit.aes = FALSE,
             shape = 21, fill = "white", colour = "#1e3a2f",
             size = 2.2, stroke = 0.8) +
  scale_fill_viridis_c(option = "magma", name = "°C", direction = -1) +
  coord_equal() +
  labs(title = "Kriged temperature", subtitle = "Ordinary kriging prediction") +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(family = "serif", face = "bold",
                                   colour = "#1e3a2f"),
        axis.title = element_blank())

# ── Kriging standard error ────────────────────────────────────────────────
p_se <- ggplot(krige_df, aes(x = lon, y = lat, fill = temp_se)) +
  geom_raster(interpolate = TRUE) +
  geom_point(data = la_pts_df,
             aes(x = lon, y = lat),
             inherit.aes = FALSE,
             shape = 21, fill = "white", colour = "#1e3a2f",
             size = 2.2, stroke = 0.8) +
  scale_fill_viridis_c(option = "cividis", name = "SE (°C)") +
  coord_equal() +
  labs(title = "Kriging uncertainty", subtitle = "Standard error of prediction") +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(family = "serif", face = "bold",
                                   colour = "#1e3a2f"),
        axis.title = element_blank())

p_pred + p_se

Kriged temperature surface across Los Angeles monitoring stations. Left: predicted temperature (ordinary kriging). Right: kriging standard error — higher uncertainty at the edges of the study domain (far from any monitoring station) and lower uncertainty near stations. This uncertainty surface is unique to kriging and invaluable for sensor network design: place new stations where the SE is highest.

📌 Kriging vs. simple interpolation methods

You might wonder: why not just use inverse distance weighting (IDW) or a spline? Kriging has three advantages. First, it is optimal — it minimises mean squared prediction error given the fitted variogram. Second, it produces uncertainty estimates (the kriging variance), which IDW and splines do not. Third, it respects spatial structure: the variogram explicitly models how quickly autocorrelation decays with distance, rather than assuming a fixed decay function. The cost is complexity — you must fit and validate the variogram model. For large or irregular domains, this is always worth it.

10.5 Universal kriging: incorporating land-use covariates

Ordinary kriging (above) interpolates using spatial autocorrelation alone. Universal kriging (also called kriging with external drift) incorporates auxiliary variables — here, impervious surface fraction — as a deterministic trend, then kriging the residuals. This often substantially improves predictions in environmentally heterogeneous areas.

# Fit a variogram on residuals from the land-use trend
vario_uk <- variogram(temp_C ~ imperv_frac, data = la_utm,
                      cutoff = 25000, width = 3000)
fit_uk   <- fit.variogram(vario_uk,
                           model = vgm(psill = 0.9, model = "Sph",
                                       range = 12000, nugget = 0.2))

# grid_pts already carries imperv_frac (added in krige-grid before gridded())
krige_uk <- krige(temp_C ~ imperv_frac,
                   locations = la_utm,
                   newdata   = grid_pts,
                   model     = fit_uk)

## [using universal kriging]

cat("Universal kriging RMSE (cross-val approx.):",
    round(sqrt(mean(krige_uk$var1.var, na.rm = TRUE)), 3), "°C\n")

## Universal kriging RMSE (cross-val approx.): 0.502 °C

11. Reporting Spatial Analysis

What every spatial analysis must document

Data source and vintage — where did the shapefiles and attribute data come from? Include URL, version, and download date.
CRS used — report the EPSG code for both the source CRS and the analysis CRS. State explicitly whether any reprojection was performed.
Spatial join logic — which predicate (st_within, st_intersects, etc.) and why. Note any features that failed to join.
Variogram model — for kriging: the model type (spherical, exponential, etc.), nugget, sill, range, and how the model was fit (eye-fit, REML, WLS).
Prediction grid — resolution and extent; acknowledge edge effects.
Validation — report leave-one-out cross-validation RMSE or hold-out set accuracy for kriging predictions.

Quick-reference: common sf workflows

library(sf)

# ── Reading / writing ─────────────────────────────────────────────────────
st_read("file.shp")                          # read shapefile
st_read("file.geojson")                      # read GeoJSON
st_write(sf_obj, "output.gpkg")              # write GeoPackage

# ── CRS ───────────────────────────────────────────────────────────────────
st_crs(sf_obj)                               # check CRS
st_transform(sf_obj, crs = 5070)             # reproject

# ── Geometry operations ───────────────────────────────────────────────────
st_buffer(sf_obj, dist = 1000)               # 1 km buffer (in metres if projected)
st_union(sf_obj)                             # dissolve all features into one
st_centroid(sf_obj)                          # centroid of each polygon
st_area(sf_obj)                              # area (m² if projected)
st_distance(points_a, points_b)              # pairwise distances

# ── Spatial join / filter ─────────────────────────────────────────────────
st_join(points, polygons, join = st_within)  # attach polygon attrs to points
st_filter(points, polygons)                  # keep only points inside polygons
st_intersection(layer_a, layer_b)            # geometric intersection

# ── Converting back to data frame ─────────────────────────────────────────
st_drop_geometry(sf_obj)                     # remove geometry, return data.frame
as.data.frame(sf_obj)                        # same effect

🔑 Key Takeaways

Spatial data violates independence. Nearby observations are correlated (Tobler’s First Law). Standard tests applied naively overstate precision. Use spatial statistics — or at minimum, map your data before analysing it.
CRS management is not optional. Always check st_crs() before overlaying layers. Reproject with st_transform(); use equal-area projections (EPSG:5070) for area/distance calculations in the continental US.
Shapefiles travel in folders. The .shp file alone is unreadable — you need .dbf, .shx, and .prj too. Prefer GeoPackage (.gpkg) for new work.
Spatial joins link layers by location, not by a shared key column. Use st_join() with the appropriate predicate (st_within, st_intersects, etc.).
Kriging is optimal spatial interpolation: it exploits the variogram to produce both a predicted surface and a prediction uncertainty surface. The uncertainty map is as valuable as the prediction — use it to guide sensor placement and communicate confidence to decision-makers.
The urban heat island signal is recoverable from routine monitoring data: impervious surface fraction and canopy cover together explain most intra-city temperature variation after controlling for regional climate.

What’s next?

Time Series — temporal autocorrelation is the time-domain twin of spatial autocorrelation; many of the same concepts apply
Multiple Regression & GLMs — spatial regression (spatial lag and spatial error models) extends OLS to handle spatial autocorrelation in residuals
Remote Sensing module (coming soon) — raster spatial data with terra: land surface temperature from Landsat, NDVI, and supervised classification