Lab in R

Everything we have mapped so far has been vector data — points, lines and polygons, each with a geometry and a row in a table. Raster data works differently: a regular grid of cells, each holding a single value, with no table of attributes at all.

That difference matters more than it sounds. A raster doesn’t know about “neighbourhoods” or “households” — it just knows that cell (1450, 2013) has the value 847. Much of working with raster data is about bridging that gap: getting from a grid of numbers to something you can join to the social and administrative data you actually care about.

What you’ll be able to do by the end

Load a raster, check its CRS, and reproject it
Crop and mask it down to an area of interest
Style it so it actually communicates something
Extract raster values at point locations — the raster → vector bridge
Compute zonal statistics — one value per administrative area, ready to join and map as a choropleth

Installing Packages

# Provides various utility functions for R programming (we use it to unzip .gz files)
library(R.utils)
# For data manipulation and transformation.
library(dplyr)
# For reshaping data (we use pivot_longer to facet maps by year)
library(tidyr)
# Spatial vector data
library(sf)
# Popular data visualization package in R.  
library(ggplot2)
# Color palettes suitable for data visualization, especially for those with color vision deficiencies.
library(viridis)
# A collection of color palettes for data visualization.
library(RColorBrewer)
# The modern package for working with raster data -- gridded spatial data like satellite imagery or elevation
library(terra)
# Tools for extracting data from raster layers at exact locations, often used in spatial analysis.
library(exactextractr)
# Lets ggplot2 and the tidyverse work directly with terra objects (SpatRaster, SpatVector)
library(tidyterra)

terra, not raster

You may see older tutorials use the raster package. It has been superseded by terra, written by the same author — terra is faster, handles larger files, and is what you should learn now. We use terra throughout this lab. If you hit code online that uses raster, the function names usually map across fairly directly (raster() → rast(), extent() → ext()).

Terrain data

Import raster data

Raster terrain data consists of gridded elevation values that represent the topography of a geographic area. You can download this from the relevant github folder. A good place to download elevation data is Earth Explorer. This video takes you through the download process if you want to try this out yourself.

We first import a raster file for elevation.

elevation <- rast("data/Lebanon/LBN_elevation_w_bathymetry.tif")

elevation

class       : SpatRaster 
size        : 708, 1150, 1  (nrow, ncol, nlyr)
resolution  : 0.0025, 0.0025  (x, y)
extent      : 33.74907, 36.62407, 33.06327, 34.83327  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : LBN_elevation_w_bathymetry.tif 
name        : LBN_elevation_w_bathymetry 
min value   :                      -2082 
max value   :                       3065

Printing the object is worth a moment. Unlike an sf object, there’s no table of rows — instead you get the dimensions (how many cells), the resolution (how much ground each cell covers), the extent (the bounding box), and the CRS. That’s the anatomy of every raster you’ll meet.

Plot it.

plot(elevation)

Note the values below zero: this file includes bathymetry (depth below sea level) as well as elevation, which is why the Mediterranean shows up rather than being blank.

Have a look at the CRS.

crs(elevation)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

Import the Lebanon shapefile

Import the Lebanon shapefile, plot it, and verify its Coordinate Reference System (CRS). Is it the same as the raster’s CRS?

Lebanon_adm1 <- read_sf("data/Lebanon/LBN_adm1.shp")
plot(Lebanon_adm1$geometry)

st_crs(Lebanon_adm1)$input

[1] "Deir ez Zor / Syria Lambert"

Reproject the Raster

We use the terra project() function. We need to define two things:

The object we want to reproject and
The CRS that we want to reproject it to.

elevation <- terra::project(elevation, crs(vect(Lebanon_adm1))) # reproject the elevation data to the crs of the Lebanon shapefile
crs(elevation, describe = TRUE)$name

[1] "Deir ez Zor / Syria Lambert"

Important

Reprojecting a raster is not the same as reprojecting a vector. Vector reprojection just moves coordinates. Raster reprojection has to build a whole new grid and estimate values for the new cells — so it necessarily changes your data slightly. Reproject once, as early as possible, and avoid doing it repeatedly.

Cropping and Masking

Cropping and masking are both spatial operations used to narrow a raster down to an area of interest — but they do different things:

Cropping

Purpose: changes the extent of the raster by cutting it down to a new bounding box. The result is a smaller, rectangular raster.
Typical Use: reducing the size of a raster to focus on a smaller geographic area while retaining all the original values within that area.

Masking

Purpose: sets cells outside a given shape to NA, keeping the extent the same. The result is the same size, but with everything outside your polygon blanked out.
Typical Use: isolating specific areas or features — for example extracting land cover within the boundaries of a protected national park.

Why do both, and in this order?

Cropping is a cheap rectangular operation; masking has to test every cell against a polygon boundary. Crop first, then mask — you shrink the problem before doing the expensive part. On a large raster this can be the difference between seconds and minutes.

elevation_lebanon <- crop(elevation, vect(Lebanon_adm1))
elevation_lebanon_mask <- mask(elevation_lebanon, vect(Lebanon_adm1))

Note vect() — that converts an sf object into terra’s own vector class (SpatVector), which is what terra functions expect.

Note

Compare the two results to see the difference for yourself:

par(mfrow = c(1, 2))
plot(elevation_lebanon, main = "Cropped only")
plot(elevation_lebanon_mask, main = "Cropped + masked")

par(mfrow = c(1, 1))

Plot elevation

plot(elevation_lebanon_mask)
plot(Lebanon_adm1$geometry, col = NA, add = TRUE)

Let’s improve this a bit. Remember that there is a lot we can do with ColorBrewer.

pal <- rev(brewer.pal(6, "Oranges"))
plot(elevation_lebanon_mask, breaks = c(-100, 0, 700, 1200, 1800, 3300), col = pal)
plot(Lebanon_adm1$geometry, col = NA, add = TRUE)

tidyterra gives us geom_spatraster(), which lets a SpatRaster slot straight into the ggplot grammar you already know:

ggplot() +
  geom_spatraster(data = elevation_lebanon_mask) +
  geom_sf(data = Lebanon_adm1, fill = NA, colour = "grey20", linewidth = 0.4) +
  scale_fill_hypso_c(palette = "dem_poster", na.value = NA, name = "Elevation (m)") +
  labs(title = "Elevation in Lebanon") +
  theme_void()

scale_fill_hypso_c() is a tidyterra scale built specifically for elevation — it uses a hypsometric palette, the green-through-brown-to-white convention you see on physical atlases. Much more legible for terrain than a generic sequential ramp.

Questions to ask yourself about how you can improve these maps, going back to geo-visualisation and choropleths:

What are the logical breaks for elevation data? Sea level is a natural one — are the others meaningful, or arbitrary?
Should the colours be changed to standard elevation palettes? (Compare the two tabs above.)
Does the reader know what the units are?

Spatial join with vector data

Here is the bridge between the two data models. You have a raster of elevation, and a set of points — survey households. You want each household to know its own elevation.

# Load some geo-localised survey data 
households <- read_sf("data/Lebanon/random_survey_LBN.shp")

# terra::extract() reads the raster value at each point location.
# 'elevation' is the raster; 'households' the point locations.
household_elevation <- terra::extract(elevation, vect(households))

Warning: [extract] transforming vector data to the CRS of the raster

# Attach elevation at each point to the original households dataframe
households <- cbind(households, household_elevation)

# Check out the data
head(households)

Simple feature collection with 6 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 35.69817 ymin: 33.72841 xmax: 36.17128 ymax: 34.33229
Geodetic CRS:  WGS 84
  id ID LBN_elevation_w_bathymetry                  geometry
1  0  1                   878.1839 POINT (35.90386 33.72841)
2  1  2                  1558.2007 POINT (36.17128 34.20268)
3  2  3                  1096.7252 POINT (35.81425 34.19914)
4  3  4                  1251.2185 POINT (36.03916 34.06442)
5  4  5                   986.2734  POINT (35.69817 34.0405)
6  5  6                  1323.2316 POINT (35.98001 34.33229)

terra::extract() returns a data frame with an ID column (which point) plus one column per raster layer — so after cbind() your points carry their elevation as an ordinary attribute. From here it is just a table: you can model it, summarise it, or map it as a normal point layer.

Important

Make sure all your data is in the same CRS, otherwise the extraction will not work properly.

If you get a column of NAs, this is almost always why. Check with crs(elevation) and st_crs(households), and reproject with terra’s project() before extracting.

Night Lights

This section is a bit more advanced, there are hints along the way to make it simpler.

Nighttime lights are one of the most widely used satellite products in social science — they are a proxy for economic activity that is available consistently, globally, and for decades, in places where reliable statistics simply don’t exist. (Recall Henderson et al. from the lecture.) Here we’ll go from raw imagery to a choropleth of light per administrative region.

Refresher: some useful R functions

This section leans on file manipulation and custom functions. If those are rusty, expand this and work through it first.

# list files 
list.files()

# list files in a specific folder
list.files(file.path("data/Lebanon/Polygons"))

# list files corresponding to a specific pattern ("shp" in the filename)
list.files(file.path("data/Lebanon/Polygons"), pattern = "shp")

# list files corresponding to a specific pattern ("shp" at the end of the filename)
shps <- list.files("data/Lebanon/Polygons", pattern = "*.shp")

# we can also select strings following a pattern inside a list or vector using grepl
shps <- shps[grepl("Lebanon", shps)] 

# let's extract the first element of the list "shps"
file1 <- shps[1]
file1

# how many characters in the filename
nchar(file1)

# let's remove the last 4 characters (the file extension)
file1_short <- substr(file1, 1, nchar(file1) - 4)

# let's add something to the name (concatenate strings) - for example, a new extension ".tif"
paste(file1_short, ".tif", sep = "")

# finally let's create a function MathOperations that first calculates the square and then adds 3
MathOperations <- function(x) {
  sq <- x^2
  z <- sq + 3
  return(z)
}

# try the function on 4, 5, 6
MathOperations(4)
MathOperations(5)
MathOperations(6)

# repeat this operation for the vector 4 to 6 (similar to a loop in STATA)
lapply(4:6, function(x) MathOperations(x))

Download data

Download the Data

We need to download some raster data. NOAA has made nighttime lights data available for 1992 to 2013. It is called the Version 4 DMSP-OLS Nighttime Lights Time Series. The files are cloud-free composites made using all the available archived DMSP-OLS smooth resolution data for calendar years. In cases where two satellites were collecting data, two composites were produced. The products are 30 arc-second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude. We can download the Average, Visible, Stable Lights, & Cloud Free Coverages for 1992 and 2013 and put them in the data/Kenya_Tanzania folder.

Important

If you have trouble downloading from NOAA, a copy of the two years we need is available from here — you need to be logged into your UoL account. Available both as the original tar archives and as ready-to-use TIFs if you’d rather skip the decompression step.

A TAR file is an archive created by tar, a Unix-based utility used to package files together for backup or distribution purposes. It contains multiple files stored in an uncompressed format along with metadata about the archive. TAR archives compressed with GNU Zip compression may become GZ, .TAR.GZ, or .TGZ files. We need to decompress them before using them.

Before you move forward download the data for 1992 and 2013. It is also good practice to create a scratch folder where you do all your unzipping.

In our example, we will only download two years, but generally, you will have to repeat the same cleaning operations many times. Therefore, to speed up the process, we are going to create a new function. The function is going to:

Decompress the files using the untar command,
List the decompressed files using list.files command (notice there are compressed files inside the TAR archive)
Identify the TIF archive files using grepl
Decompress using the gunzip command.

We are then going to run the function on all the TAR files.

Note

You can do these steps manually if you can’t get the below chunk to work.

datafolder <- file.path("./data") # define the location of the data folder

# list downloaded files: they are compressed files using the "tar" format
tars <- list.files(file.path("data/Kenya_Tanzania/scratch"), pattern = "*.tar")

# unzip
UnzipSelect <- function(i) {
  untar(file.path(datafolder,"Kenya_Tanzania/scratch",i), exdir = file.path(datafolder, "Kenya_Tanzania/scratch")) # unzip
  all.files <- list.files(file.path(datafolder,"Kenya_Tanzania/scratch"), pattern = paste0(substr(i, 6, 12), "*")) # list extracted files
  gz <- all.files[grepl("web.stable_lights.avg_vis.tif.gz", all.files)] # select the TIF files

  destfile <- file.path(datafolder, "Kenya_Tanzania", substr(gz, 1, nchar(gz) - 3))

  # If we've already decompressed this file on a previous run, skip it.
  # Without this the chunk only works once, then errors on every re-render.
  if (file.exists(destfile)) return(destfile)

  # gunzip writes to "<destfile>.tmp" first and refuses to start if that
  # already exists -- clear any left behind by an interrupted run.
  tmpfile <- paste0(destfile, ".tmp")
  if (file.exists(tmpfile)) file.remove(tmpfile)

  R.utils::gunzip(filename  = file.path(datafolder,"Kenya_Tanzania/scratch", gz),
                  destname = destfile,
                  overwrite = TRUE) # unzip again
}

# loop over the TAR files
# note that the function returns the last element created - in this example, the TIF files
nl <- lapply(tars, UnzipSelect)

# you can delete the scratch folder with the data we don't need
# unlink(file.path(datafolder,"Kenya_Tanzania/scratch"), recursive = TRUE)

If this chunk errors with “File already exists”

gunzip() decompresses into a temporary .tif.tmp file before renaming it, and it refuses to run if that temp file is already there — which happens if a previous attempt was interrupted part-way. Note that overwrite = TRUE does not cover this: it applies to the final destination, not the temp file.

The guards above handle it automatically. If you’re still stuck, delete the leftovers by hand and re-run:

# remove any stray .tmp files from interrupted runs
file.remove(list.files(file.path(datafolder, "Kenya_Tanzania"),
                       pattern = "\\.tmp$", full.names = TRUE))

We can load and plot the nighttime lights data. When working with many rasters of the same origin, it is faster to hold them together as a single multi-layer object — in terra, rast() does this in one call.

# load the night lights rasters as one multi-layer SpatRaster
# (only possible for rasters of the same extent and resolution)
nl_rasters <- rast(unlist(nl))

# change the layer names
names(nl_rasters) <- c("NL1992", "NL2013")

# plot the result
plot(nl_rasters, axes = FALSE)

Why can’t you see much? Discuss with the person next to you.

Stuck?

These are global rasters, and you’re looking at the whole planet. Nearly all of it is dark, and the few bright pixels are tiny at this scale. Two things are going on: the extent is far larger than our area of interest, and the value distribution is extremely skewed — a handful of very bright cells compress everything else into the bottom of the colour ramp. We’ll fix the first with zonal statistics, and the second with fixed breaks when we map.

Country shapefiles

The second step is to download the shapefiles for Kenya and Tanzania. GADM has made available national and subnational shapefiles for the world. The zips you download, such as gadm36_KEN_shp.zip from GADM should be placed in the Kenya_Tanzania folder. This is the link gadm.

# list country shp that we downloaded from the GADM website
files <- list.files(file.path(datafolder,"Kenya_Tanzania"), pattern = "_shp.zip*", recursive = TRUE, full.names = TRUE)
files

# create a scratch folder
# dir.create(file.path(datafolder,"Kenya_Tanzania/scratch"))

# unzip
lapply(files, function(x) unzip(x, exdir = file.path(datafolder,"Kenya_Tanzania/scratch")))

# GADM has shapefiles for different regional levels (e.g. country, region, district, ward) 
gadm_files <- list.files(file.path(datafolder,"Kenya_Tanzania"), pattern = "gadm*", recursive = TRUE, full.names = TRUE)
gadm_files

# let's select regional level 2
gadm_files_level2 <- gadm_files[grepl("2.shp", gadm_files)]
gadm_files_level2

# load the shapefiles
shps <- lapply(gadm_files_level2, read_sf)
shps

# delete the scratch folder with the data we don't need
# unlink(file.path(datafolder,"Kenya_Tanzania/scratch"), recursive = TRUE)

Zonal statistics

This is the step that turns a raster into something you can treat like any other table. Zonal statistics summarise raster cells within each polygon — giving you one number per administrative region.

We use the package exactextractr, which is fast and, unlike simpler approaches, correctly handles cells that only partly fall inside a polygon by weighting them by the fraction covered.

Again, we use the lapply function to process the two countries successively.

# summarize
ex <- lapply(shps, function(x) exact_extract(nl_rasters, x, c("sum", "mean", "count"), progress = FALSE))
# lapply returns a list of two dataframes, we can use "do.call" to return each element of the list and iterate the function rbind
# the results is a dataframe with the merged rows of the dataframes
ex <- do.call("rbind", ex)

# show first files
head(ex)

  sum.NL1992 sum.NL2013 mean.NL1992 mean.NL2013 count.NL1992 count.NL2013
1          0          0  0.00000000   0.0000000     203.3606     203.3606
2         29        109  0.03748985   0.1409101     773.5427     773.5427
3          0          0  0.00000000   0.0000000    1910.4520    1910.4520
4          0          0  0.00000000   0.0000000    2205.3706    2205.3706
5          0        150  0.00000000   0.1415757    1059.5037    1059.5037
6          0          0  0.00000000   0.0000000    1353.7744    1353.7744

# summary
summary(ex)

   sum.NL1992         sum.NL2013        mean.NL1992         mean.NL2013       
 Min.   :   0.000   Min.   :    0.00   Min.   : 0.000000   Min.   : 0.000000  
 1st Qu.:   0.000   1st Qu.:   34.61   1st Qu.: 0.000000   1st Qu.: 0.008707  
 Median :   1.322   Median :  211.05   Median : 0.000366   Median : 0.105727  
 Mean   : 204.362   Mean   :  629.85   Mean   : 2.180426   Mean   : 3.671596  
 3rd Qu.: 158.505   3rd Qu.:  651.81   3rd Qu.: 0.175598   3rd Qu.: 1.070672  
 Max.   :8505.548   Max.   :17562.79   Max.   :60.147682   Max.   :62.272068  
  count.NL1992        count.NL2013      
 Min.   :    0.002   Min.   :    0.002  
 1st Qu.:  321.951   1st Qu.:  321.951  
 Median :  952.327   Median :  952.327  
 Mean   : 3706.633   Mean   : 3706.633  
 3rd Qu.: 4589.101   3rd Qu.: 4589.101  
 Max.   :47327.062   Max.   :47327.062

Tip

Note we asked for three statistics — sum, mean and count. Which one you want depends on the question. Sum of light is closer to total economic activity; mean controls for the fact that regions differ enormously in size (the same MAUP issue you met with choropleths). We map the mean below — try the sum and see how differently the map reads.

Merge shapefiles

Even though it is not necessary here, we can merge the shapefile to visualize all the regions at once.

Usually, it is easier to process data in small chunks using a function like sapply, lapply, mapply or a loop before merging. For example, when doing zonal statistics, it is faster and easier to process one country at a time and then combine the resulting tables. If you have access to a computer with multiple cores, it is also possible to do “parallel processing” to process each chunk at the same time in parallel.

# merge together
# we select each sf object and merge the rows
# do.call() in R applies a given function to a list as a whole
# rbind() can be used to bind or combine several vectors, matrices, or data frames by rows
tza_ken <- do.call("rbind", shps)

# inspect
str(tza_ken)

sf [484 × 14] (S3: sf/tbl_df/tbl/data.frame)
 $ GID_0    : chr [1:484] "KEN" "KEN" "KEN" "KEN" ...
 $ NAME_0   : chr [1:484] "Kenya" "Kenya" "Kenya" "Kenya" ...
 $ GID_1    : chr [1:484] "KEN.1_1" "KEN.1_1" "KEN.1_1" "KEN.1_1" ...
 $ NAME_1   : chr [1:484] "Baringo" "Baringo" "Baringo" "Baringo" ...
 $ NL_NAME_1: chr [1:484] NA NA NA NA ...
 $ GID_2    : chr [1:484] "KEN.1.1_1" "KEN.1.2_1" "KEN.1.3_1" "KEN.1.4_1" ...
 $ NAME_2   : chr [1:484] "805" "Baringo Central" "Baringo North" "Baringo South" ...
 $ VARNAME_2: chr [1:484] NA NA NA NA ...
 $ NL_NAME_2: chr [1:484] NA NA NA NA ...
 $ TYPE_2   : chr [1:484] "Constituency" "Constituency" "Constituency" "Constituency" ...
 $ ENGTYPE_2: chr [1:484] "Constituency" "Constituency" "Constituency" "Constituency" ...
 $ CC_2     : chr [1:484] "162" "159" "158" "160" ...
 $ HASC_2   : chr [1:484] NA NA NA NA ...
 $ geometry :sfc_MULTIPOLYGON of length 484; first list element: List of 1
  ..$ :List of 1
  .. ..$ : num [1:1017, 1:2] 35.9 35.9 35.9 35.9 35.9 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr [1:13] "GID_0" "NAME_0" "GID_1" "NAME_1" ...

# plot
plot(tza_ken$geometry)

Visualize

Let’s have a first look at our result. Note that from here on, this is just a choropleth — exactly what you built in the vector labs. The raster work is done; the output is an ordinary sf object with numeric columns.

# merge back with shapefile attribute table
# this time instead of merging the rows, we append the columns using cbind
df <- cbind(tza_ken, ex)

ggplot(df) +
  geom_sf(aes(fill = mean.NL1992), colour = "grey60", linewidth = 0.1) +
  scale_fill_viridis_c(name = "Mean 1992", option = "inferno") +
  labs(title = "Mean nighttime lights, 1992",
       subtitle = "Continuous scale — almost everything looks identical") +
  theme_void()

That map is nearly useless, and deliberately so. A handful of very bright regions (Nairobi, Dar es Salaam) stretch the colour scale so far that all the actual variation across the rest of the two countries is squashed into the bottom of the ramp.

The distribution shows why:

ggplot(df, aes(x = mean.NL1992)) +
  geom_histogram(bins = 60, fill = "grey30") +
  labs(x = "Mean nighttime lights, 1992", y = "Number of regions",
       title = "Extremely skewed — most regions are near zero") +
  theme_minimal()

To make the maps tell a story, we need fixed breaks placed where the variation actually is — bunched up near zero, rather than spread evenly from 0 to 63.

# reshape to long format so we can facet the two years side by side
df_long <- df |>
  pivot_longer(
    cols = c(mean.NL1992, mean.NL2013),
    names_to = "year",
    values_to = "mean_nl"
  ) |>
  mutate(
    year = recode(year, "mean.NL1992" = "1992", "mean.NL2013" = "2013"),
    # cut() into explicit classes -- 5 break points give 4 bins,
    # so we need exactly 4 labels and 4 colours below
    nl_class = cut(
      mean_nl,
      breaks = c(0, 0.05, 0.1, 2, 63),
      labels = c("0 – 0.05", "0.05 – 0.1", "0.1 – 2", "2 – 63"),
      include.lowest = TRUE
    )
  )

ggplot(df_long) +
  geom_sf(aes(fill = nl_class), colour = "grey70", linewidth = 0.1) +
  scale_fill_manual(
    name = "Average\nnightlights",
    values = c("#08306B", "#2C7FB8", "#7FCDBB", "#FFFFB2"),
    na.value = "grey90"
  ) +
  facet_wrap(~ year) +
  labs(title = "Nighttime lights, Kenya & Tanzania",
       subtitle = "Fixed breaks — now the variation is visible") +
  theme_void()

Why cut() rather than a continuous scale with breaks

Binning explicitly with cut() and mapping with scale_fill_manual() gives you a discrete variable, so you control exactly how many classes there are and which colour each one gets. Getting this wrong is a classic trap: n break points define n − 1 bins, so a mismatch between the number of breaks and the number of colours silently produces a map where classes blur together — exactly the flat map we were trying to escape.

If you’d rather keep a continuous scale, the equivalent is scale_fill_fermenter() or scale_fill_stepsn() — just count your bins carefully.

Alternative: a log scale

Fixed breaks are one answer to a skewed variable. Another is to transform the scale itself, which avoids having to pick break points by hand:

ggplot(df_long) +
  geom_sf(aes(fill = mean_nl + 0.01), colour = "grey70", linewidth = 0.1) +
  scale_fill_viridis_c(
    name = "Average\nnightlights",
    trans = "log10", option = "inferno",
    labels = scales::label_number(accuracy = 0.01)
  ) +
  facet_wrap(~ year) +
  labs(title = "Nighttime lights on a log scale") +
  theme_void()

Note the + 0.01: log10(0) is undefined, so regions with zero light would drop out entirely. Adding a small constant keeps them in — but be honest about it, since it’s a choice that affects what the reader sees.

Tip

Have a think about what the data is telling you. What’s the story? And — carefully — what isn’t it telling you? Brighter isn’t automatically richer.

Where this fits: rasters in the wild

You have just done, by hand and on a small scale, roughly what a data service does at national scale and industrial volume:

Raw imagery arrives as a noisy grid, partly obscured by cloud
It gets cleaned — many passes combined into a composite
A model turns spectral values into a meaningful indicator — surface temperature, vegetation, cloud probability, light
Pixels are aggregated to administrative units, so the data can be joined to censuses, surveys and policy frameworks

Steps 1–2 we skipped (the NOAA composites arrive pre-cleaned). Step 3 was done for us. Step 4 is exactly what exact_extract() did above.

Want to go further with satellite data?

This is the whole business of Imago, the imagery data service for sustainability, prosperity and wellbeing, part of Smart Data Research UK. They run that pipeline and publish the output as ready-to-use LSOA/MSOA-level statistics — no gigabyte downloads, no remote-sensing algorithms, small-area detail preserved.

From grids to areas — the full version of the four steps above, with the spectral bands / resolution / revisit-cycle background behind them. Start here.
Imago training — free, openly licensed trainings in both R and Python, covering air temperature, precipitation and SPF. Several run directly in your browser with no installation.

Worth knowing: the SPF data you mapped in the choropleths lab was satellite-derived cloud probability that had already been through this whole pipeline. You were doing raster analysis without ever touching a raster — which is rather the point of a service like this.

Resources

Downloading night lights
The package nightlightstats
terra documentation — the reference for everything in this lab
tidyterra — for using rasters inside ggplot
Spatial Data Science with R and terra — free online book