Skip to main content

Module 6 lab: putting data on a map

Work through this at your own pace. Each step ends with a check – something you should be able to see on your screen before moving on. If a check fails, that is the point to put your hand up rather than pushing on.

This lab is all in R. There is no Excel half, because Excel does not read a shapefile and will not reproject anything.

The shape of the day: read boundaries, draw them, colour them by a number you joined on, then ask the map questions it can answer that a table cannot – how big, how far, what is next to what, and what is inside what.

Plan on about an hour and a quarter for the ten steps, and half an hour more for the Practice boxes. It covers 19  Spatial Data, 20  Reading and Mapping Spatial Data and 21  Spatial Operations and Maps. If you run out of time, Step 8 onwards is the part to finish at home.

Step 1: Install the spatial packages and get the data

Open Positron and open your AREC_261 folder, the one from the Module 2 lab. Every read_csv("data/...") and st_read("data/...") below is relative to the folder Positron has open.

Spatial work needs packages you have not installed yet. In the console – not the script, since you only do this once:

install.packages(c("sf", "terra", "tidyterra", "exactextractr"))

This takes a few minutes and prints a great deal. On Mac it may compile from source and show warnings in red that are not errors.

Check: library(sf) runs and prints something like Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1. Those three are the C libraries that do the actual geometry – sf is a wrapper around them. If this line errors, stop here; nothing below will work.

Download these into data/:

Make a new script in code/ called module6_lab.R:

# ---
# Title: Module 6 lab
# Author: Your Name
# Date: 2026-10-21
# Description:
#   Reading, mapping and querying spatial data.
# ---

# Load packages
library(tidyverse)
library(sf)
library(terra)
library(tidyterra)
library(janitor)

Step 2: Read the RM boundaries and look at what you got

# 296 rural municipalities. quiet = TRUE stops st_read printing the file path.
rms <- st_read("data/sask_rms.geojson", quiet = TRUE) |>
  clean_names()

rms

Check: 296 rows and 3 columns – rmno, rmnm and geometry. The print-out above the table says Geometry type: MULTIPOLYGON and Geodetic CRS: WGS 84.

Printing an sf object gives you a header the ordinary data frame print does not: geometry type, bounding box, and CRS. Read that header every time you open a spatial file. It is where the problems announce themselves.

Now look at the geometry column itself:

rms$geometry[[1]]

Check: a long list of coordinate pairs around -106 and 49 – degrees of longitude and latitude. That is one RM’s outline, stored as the corners of its polygon.

An sf object is an ordinary data frame with one extra column holding the shape. Everything you know from Module 2 still works on it:

# the dplyr verbs work unchanged
rms |>
  filter(rmnm == "ABERDEEN")

Check: one row, RM number 373, and the geometry column comes along with it.

Step 3: Draw it

ggplot(rms) +
  geom_sf()

Check: the outline of Saskatchewan, filled with 296 polygons. It takes a second or two to draw.

That is the whole trick: geom_sf() knows to look for the geometry column, so you do not name an x or a y. Tidy it up:

ggplot(rms) +
  geom_sf(fill = "grey95", colour = "grey70", linewidth = 0.2) +
  theme_void()

theme_void() removes the axes. Latitude and longitude gridlines around a provincial map are clutter – nobody reads an RM map off its coordinates.

Check: a clean provincial map, pale fill, thin grey borders, no axes or grey background.

Highlight one RM by drawing it as a second layer on top:

aberdeen <- filter(rms, rmnm == "ABERDEEN")

ggplot() +
  geom_sf(data = rms, fill = "grey95", colour = "grey70", linewidth = 0.2) +
  geom_sf(data = aberdeen, fill = "darkgreen") +
  theme_void()

Check: one dark green RM, north-east of Saskatoon.

Note the change: ggplot() is now empty and each geom_sf() names its own data. With two different datasets on one map there is no single dataset to put in ggplot().

Step 4: Join a number onto the shapes and colour by it

A map of outlines says where things are. A map becomes useful when the colour carries a number.

yields <- read_csv("data/rm_yields_1990_2025.csv", show_col_types = FALSE) |>
  clean_names()

glimpse(yields)

Check: 71,104 rows and 5 columns – year, rm, crop, yield, unit. This is long data: one row per RM per crop per year.

Cut it to one crop in one year, then join:

canola_2024 <- yields |>
  filter(year == 2024, crop == "Canola")

nrow(canola_2024)

Check: 293 rows. There are 296 RMs, so three are already missing.

# rmno is text in the geojson and rm is a number in the CSV
rms <- mutate(rms, rmno = as.numeric(rmno))

rm_canola <- left_join(rms, canola_2024, by = c("rmno" = "rm"))

sum(is.na(rm_canola$yield))

Check: 296 rows out, 3 with a missing yield.

Two things worth noticing. The join key needed a type conversion – rmno arrived as text and rm is numeric, and without the as.numeric() the join fails with a type error. And left_join with the sf object on the left keeps the geometry; put the CSV first and you get a plain data frame back with the shapes dropped.

ggplot(rm_canola) +
  geom_sf(aes(fill = yield), colour = NA) +
  scale_fill_viridis_c(option = "magma", na.value = "grey90") +
  labs(title = "Canola yield by RM, 2024", fill = "bu/ac") +
  theme_void()

Check: a filled map, pale in the south-west and bright through the centre and east. Three grey RMs where the yield was missing.

na.value = "grey90" is the line that matters. Without it the missing RMs are drawn in the same background grey as the panel and look like holes in the province rather than like absent data.

Practice

Change colour = NA to colour = "white" and redraw. Then try scale_fill_viridis_c() without the option argument. Which version reads better from across the room?

White borders separate the RMs and make individual ones findable, at the cost of some ink; with 296 small polygons they start to dominate. colour = NA suits a map read as a surface, white borders suit one where a reader looks up their own RM.

The default viridis runs dark blue to yellow, magma runs black through red to pale. Both are colourblind-safe and both are perceptually even, which the old rainbow palettes were not – a rainbow makes some equal steps look larger than others, which is why it has fallen out of use for exactly this kind of map.

Step 5: Measure the shapes

The geometry column is not just for drawing. It can be measured.

st_area(aberdeen)

Check: about 714,217,471 square metres, printed with [m^2] after it.

Square metres is a useless unit for an RM, and the number carries its unit around with it, which gets in the way of arithmetic. Convert:

as.numeric(st_area(aberdeen)) / 1e6

Check: 714.2175 – square kilometres.

as.numeric() strips the unit. Do it to every RM at once:

rms <- mutate(rms, km2 = as.numeric(st_area(geometry)) / 1e6)

median(rms$km2)
sum(rms$km2 > 714.2175)

Check: the median RM is 853.03 km², and 268 of the 296 are bigger than Aberdeen.

So Aberdeen is a small RM – bottom tenth by area. You could not have read that off the map; the eye is poor at comparing areas, especially irregular ones.

NoteWhy the area is right even though the CRS is in degrees

The file is in WGS 84, whose coordinates are degrees, and a degree is not a fixed distance – it narrows toward the poles. Measuring area in degrees would be meaningless.

st_area() knows this. For an unprojected CRS it computes the area on the ellipsoid and returns metres anyway. Not every function is so careful, which is why 21  Spatial Operations and Maps reprojects before doing geometry.

Step 6: Ask what is next to what

st_touches(aberdeen, rms)

Check: a list with one element, containing a single number.

One neighbour. Look at the map from Step 3 and that is plainly wrong – Aberdeen has several RMs around it. Try the other predicate:

neighbours <- st_intersects(aberdeen, rms)[[1]]
rms$rmnm[neighbours]

Check: five names – GRANT, ROSTHERN, BLUCHER, ABERDEEN, CORMAN PARK.

Five, including Aberdeen itself, so four genuine neighbours against the one that st_touches() found.

The difference is the definition. st_touches() is true only when two polygons share boundary points exactly and their interiors do not overlap. st_intersects() is true when they share any point at all. These boundaries were digitised from paper maps, so adjacent RMs are a metre or two apart in places, or overlap slightly. Exact touching is a strict test and real data rarely passes it.

The predicate you want is almost always st_intersects(), and the self-match is the thing to remember to drop:

neighbours <- setdiff(neighbours, which(rms$rmnm == "ABERDEEN"))

ggplot() +
  geom_sf(data = rms, fill = "grey95", colour = "grey70", linewidth = 0.2) +
  geom_sf(data = rms[neighbours, ], fill = "lightblue") +
  geom_sf(data = aberdeen, fill = "darkgreen") +
  theme_void()

Check: Aberdeen in green with four light blue RMs around it.

Practice

Find the neighbours of CORMAN PARK – the RM that surrounds Saskatoon – and count them. Then compare st_touches() and st_intersects() for it.

Corman Park is one of the largest RMs in the province and borders a great many others, so the count is well into double figures with st_intersects().

The gap between the two predicates is not a fixed number – it depends on how each pair of boundaries happened to be digitised. That is the real lesson: a result that depends on exact coordinate agreement is a result that depends on who drew the file.

Step 7: Reproject, and see why it matters

Every layer on one map must be in the same CRS. Read the risk zones:

zones <- st_read("data/sk_risk_zones.geojson", quiet = TRUE)

st_crs(rms)$Name
st_crs(zones)$Name

Check: both say WGS 84.

These two agree, so they can be drawn together:

ggplot() +
  geom_sf(data = rms, fill = NA, colour = "grey85", linewidth = 0.15) +
  geom_sf(data = zones, fill = NA, colour = "firebrick", linewidth = 0.5) +
  theme_void()

Check: the RM mesh in pale grey with 23 heavier red zone outlines over it.

Now reproject one of them and look at the shape:

rms_utm <- st_transform(rms, 32613)      # UTM zone 13N, metres

ggplot(rms_utm) +
  geom_sf(fill = "grey95", colour = "grey70", linewidth = 0.2) +
  theme_void()

Check: Saskatchewan again, but visibly narrower at the top and leaning – the province is no longer a neat trapezoid.

Both pictures are correct. WGS 84 plotted directly stretches the north; UTM 13N is accurate near its central meridian and distorts as you move east and west of it. There is no undistorted flat map of a curved surface, and the question is only which distortion you can live with.

Check: st_crs(rms_utm)$units returns "m". Coordinates are now metres, so distances and areas come out in sensible units without the ellipsoid calculation.

Step 8: Read a raster

Everything so far has been vector – points, lines, polygons. A raster is a grid of cells.

aci <- rast("data/aci_slice_2023.tif")

aci

Check: the print-out reports size: 2997, 2745, 1, resolution: 30.02163, 30.02163, and coord. ref. : WGS 84 / UTM zone 13N.

Read those three lines. The grid is 2,997 rows by 2,745 columns – 8.2 million cells. Each cell is about 30 m by 30 m on the ground, which is the nominal resolution of the satellite this was built from. And the CRS is not the one the RMs are in.

ggplot() +
  geom_spatraster(data = aci) +
  theme_void()

Check: a square of blocky colour. It takes several seconds, and the legend is a continuous colour bar running to about 200.

The legend is wrong, and it is wrong in an instructive way. The cell values are codes – 146 means spring wheat, 153 means canola – and a continuous colour bar implies 153 is more than 146 in some meaningful sense. It is not. The file:

classes <- read_csv("data/aci_crop_classes.csv", show_col_types = FALSE)

head(classes)

Check: 75 rows with value, crop and colour columns.

Step 9: Cut the raster to one RM

The raster is a square slab; Aberdeen is an irregular polygon. Two operations, in order.

The CRS has to match first:

aberdeen_utm <- st_transform(aberdeen, crs(aci))

ab_crop <- crop(aci, vect(aberdeen_utm))
ab_mask <- mask(ab_crop, vect(aberdeen_utm))

ab_mask

Check: dimensions now 651, 506, 1 – down from 8.2 million cells to about 330,000.

Three things happened. st_transform() put the polygon into the raster’s CRS, because crop() will not reconcile them for you and either errors or returns an empty raster. vect() converted the sf polygon into terra’s own format – sf and terra are separate packages with separate object types. And then crop() cut the bounding box while mask() set everything outside the polygon to missing.

crop() alone leaves a rectangle. mask() alone keeps all 8.2 million cells and blanks most of them. Both, in that order, is the pattern.

ggplot() +
  geom_spatraster(data = ab_mask) +
  geom_sf(data = aberdeen_utm, fill = NA, colour = "white", linewidth = 0.4) +
  theme_void()

Check: the crop pattern inside Aberdeen’s outline, with blank space in the corners of the panel.

Now count what is in there:

counts <- freq(ab_mask) |>
  as_tibble() |>
  left_join(classes, by = "value") |>
  arrange(desc(count))

head(counts, 5)

Check: spring wheat first with 26,447 cells, then canola 14,767, barley 11,745, wetland 5,639, water 4,364.

Cells are about 30 m square, so roughly 900 m² each. Turn the counts into area:

counts |>
  mutate(km2 = round(count * 900 / 1e6, 1)) |>
  select(crop, km2) |>
  head(5)

Check: spring wheat about 23.8 km², canola about 13.3.

Compare that with Step 5: Aberdeen is 714 km², and the classified cells here total about 63 km². The slice in this file covers only a corner of the RM. A number computed from a clipped raster is a number about the part of the polygon the raster covered, and nothing warns you – the code runs and returns a confident answer either way.

Practice

Run the same crop-mask-count sequence for BLUCHER, Aberdeen’s neighbour to the south-east. How much of it does the raster cover?

Compare the total classified area against st_area() for the same RM. Whichever RM you pick, do that comparison before quoting any figure from a raster clip – it is the check that catches the error, and it is two lines.

Step 10: Make one map worth handing in

Put the pieces together. One crop, one year, properly labelled:

rm_canola |>
  ggplot() +
  geom_sf(aes(fill = yield), colour = "white", linewidth = 0.1) +
  geom_sf(data = zones, fill = NA, colour = "grey20", linewidth = 0.4) +
  scale_fill_viridis_c(option = "magma", na.value = "grey90") +
  labs(
    title = "Canola yield by rural municipality, 2024",
    subtitle = "Crop insurance risk zones outlined in dark grey",
    fill = "bu/ac",
    caption = "Source: Saskatchewan Crop Insurance"
  ) +
  theme_void()

ggsave("output/canola_2024_map.png", width = 8, height = 6, dpi = 300)

Check: output/canola_2024_map.png exists and opens. If you get Error: Cannot find directory output, make it with dir.create("output").

Then read your own map the way Step 9 of the Module 4 lab had you read a chart:

  1. Does the title say what is mapped, where and when?
  2. Is the legend labelled with a unit?
  3. Can a reader tell grey-for-missing from the low end of the colour scale?
  4. Is the extra layer earning its ink, or is it just lines?
  5. Would a table of the top ten RMs answer the question better?

Check: you found at least one thing to change. Change it, re-run, re-save.

The last question matters more for maps than for charts. A map is the right answer when the spatial pattern is the finding – when where something happens is what you are trying to show. If the reader’s real question is “which RM had the best canola”, that is a sorted table, and a map makes them hunt for it.

Practice on your own

The test bank for this module is here, and it is the best use of the remaining time.

Two other things if you want a break from the bank.

Map a different crop and year from rm_yields_1990_2025.csv and put the two side by side. Durum against canola is a good pair – they are grown in different parts of the province, and the two maps make that obvious in a way the table does not.

Then take the Step 10 map and find out how much of it survives being printed in black and white. Add scale_fill_viridis_c(option = "magma") and convert the PNG to greyscale, or just print it. Viridis palettes are designed to stay monotonic in greyscale, which most palettes do not – it is the reason to prefer them beyond the colourblind-safety.

Before you leave

Save your script. You should have one that reads a vector file and a raster, joins a table onto the shapes, measures areas, finds neighbours, reprojects, clips a raster to a polygon, and saves a finished map to output/.

The habit to keep is reading the header. Every spatial file announces its geometry type, its extent and its CRS when you print it, and nearly every spatial bug is one of those three not being what you assumed.