Repository navigation
Expand file tree
/
Copy pathmapvectorR.qmd
More file actions
420 lines (308 loc) · 20.8 KB
/
Copy pathmapvectorR.qmd
File metadata and controls
420 lines (308 loc) · 20.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
# Lab in R {#sec-vector-data-R .unnumbered}
## Choropleths
In this session, we will build on all we have learnt so far about loading and manipulating (spatial) data and apply it to one of the most commonly used forms of spatial analysis: choropleths. Remember these are maps that display the spatial distribution of a variable encoded in a color scheme, also called palette. Although there are many ways in which you can convert the values of a variable into a specific color, we will focus in this context only on a handful of them, in particular:
- Unique values
- Equal interval
- Quantiles
- Fisher-Jenks
## Installing Packages
Before all this mapping fun, let us get the importing of libraries and data loading out of the way:
```{r message=FALSE}
# Load the 'sf' library, which stands for Simple Features, used for working with spatial data.
library(sf)
# Load the 'tidyverse' library, a collection of packages for data manipulation and visualization.
library(tidyverse)
# The 'readr' library provides a fast and user-friendly way to read data from common formats like CSV.
library(readr)
# RColorBrewer library for creating visually appealing color schemes for plots and data visualizations
library(RColorBrewer)
# Working with class intervals and classification methods, esp in the context of spatial data analysis.
library(classInt)
```
## Data
We will be using data from the [**Imago Data Service**](https://data.imago.ac.uk) for this section — specifically the **Sun Probability Framework (SPF)**, a UK-wide dataset of annual cloud probability estimates for every small area in the country. Higher values mean cloudier conditions (a lower probability of direct sunlight); lower values mean clearer conditions. We'll be using the 2025 release.
::: callout-note
## Download the data
Download the SPF 2025 GeoPackage from the Imago Data Service:
<https://data.imago.ac.uk/datasets/cloud-probability-statistics-per-small-area-in-2025-version-2-0>
Place it inside a `data/` directory. This file already contains both the small-area geometries *and* the SPF value in one GeoPackage — no separate join needed to get started, and `sf` reads GeoPackages natively with `read_sf()`.
:::
```{r}
# A GeoPackage can bundle multiple layers/tables -- check what's actually
# in the file before reading, since read_sf() silently picks the first
# layer if you don't specify one, which may not be the spatial layer
st_layers("data/Imago_spf/Cloud probability statistics per small area in 2025 (GeoPackage).gpkg")
```
```{r}
spf_25 <- read_sf("data/Imago_spf/Cloud probability statistics per small area in 2025 (GeoPackage).gpkg",
layer = "SPF_LSOA_level") # replace with the actual layer name from st_layers() above
```
Don't forget that before you go further, you want to check the `CRS` of the `sf` object as well as the `dataframe`.
::: {.callout-tip collapse="true"}
## Answer
```{r crs}
st_crs(spf_25)
head(spf_25)
```
:::
The principal variables we'll use throughout are:
- `data_zone_code`: a unique small-area identifier, harmonised across England, Wales, Scotland and Northern Ireland — the first letter tells you which nation an area belongs to (`E` = England, `W` = Wales, `S` = Scotland, `N` = Northern Ireland)
- `cloud_probability`: the annual average SPF value
- `geometry`: the small-area boundary
Since our "Unique values" example needs a genuinely categorical variable, and SPF itself is continuous, we'll derive one from the `data_zone_code` prefix — which UK nation each small area belongs to:
```{r}
spf_25 <- spf_25 %>%
mutate(nation = case_when(
str_starts(data_zone_code, "E") ~ "England",
str_starts(data_zone_code, "W") ~ "Wales",
str_starts(data_zone_code, "S") ~ "Scotland",
str_starts(data_zone_code, "N") ~ "Northern Ireland",
TRUE ~ "Unknown"
))
```
Now we are fully ready to map!
We will be using `ggplot2` throughout the module — see the documentation [here](https://ggplot2.tidyverse.org). There are other good options for thematic mapping in R, such as `tmap` (see [here](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html)) and `mapsf` (see [here](https://riatelab.github.io/mapsf/)), but we'll stick to `ggplot2` here for consistency with the rest of the course.
## Unique values
A choropleth for categorical variables simply assigns a different color to every potential value in the series. Variables could be both nominal or ordinal.
**Nominal**: Nominal variables represent categories or labels *without* any inherent order or ranking. The categories are distinct and do not have a natural progression or hierarchy, such as "apple," "banana," and "orange" for fruit types.
**Ordinal** : Ordinal variables represent categories or labels with a meaningful *order* or *ranking*. The relative order or hierarchy among the categories is significant, indicating a clear progression from lower to higher values, such as "low," "medium," and "high" for satisfaction levels.
In `R`, creating categorical choropleths is possible with one line of code. To demonstrate this, we can plot which UK nation each small area belongs to (`nation`, which we derived above). `nation` is nominal — there's no inherent ranking between England, Scotland, Wales and Northern Ireland. `geom_sf` is calling the geometric object and `fill` is defining what values we want to fill the polygons with.
```{r}
ggplot(data = spf_25) +
geom_sf(aes(fill = nation), color = NA) +
theme_void()
```
## Equal Interval
If, instead of categorical variables, we want to display the geographical distribution of a continuous phenomenon, we need to select a way to encode each value into a color. One potential solution is applying what is usually called "equal intervals". The intuition of this method is to split the range of the distribution, the difference between the minimum and maximum value, into equally large segments and to assign a different color to each of them according to a palette that reflects the fact that values are ordered.
Creating the choropleth is relatively straightforward in `R`. For example, to create an equal interval map of `cloud_probability`.
First we need to prepare the data, going back to our [data wrangling](openscienceR.html).
```{r}
#| message: false
#| warning: false
spf_filtered <- spf_25 %>%
# Step 1: Filter out rows where 'cloud_probability' is missing (i.e., remove NA values)
filter(!is.na(cloud_probability)) %>%
# Step 2: Round to a whole number for a cleaner legend
mutate(cloud_probability = round(cloud_probability))
```
Now let's map using equal intervals.
Mapping in `ggplot` can be a bit tricky at the beginning. You will want to take a look at the package `classInt` [here](https://www.rdocumentation.org/packages/classInt/versions/0.4-10/topics/classIntervals).
**Step 1**: Calculate equal interval breaks with function `classIntervals` and store them.
```{r}
#| warning: false
e_breaks <- classIntervals(spf_filtered$cloud_probability, n = 7, style = "equal")
# Assign the class breaks to the data
spf_filtered$e_breaks <- cut(spf_filtered$cloud_probability, e_breaks$brks)
```
- `e_breaks`: This is a variable name that you are assigning to store the result of the class intervals
7438
calculation.
- `classIntervals()`: This is a function that calculates class intervals for a given numeric vector. It is provided by the `classInt` package.
- `spf_filtered$cloud_probability`: This is the data vector that you want to create class intervals for. `n = 7`: This argument specifies that you want to divide the data into 7 classes. `style = "equal"`: This argument specifies that you want to use equal interval classification, which means that the data range will be divided into equal-sized intervals. These class intervals can be useful for creating data visualizations like **choropleth maps** or **histograms**.
**Step 2**: Finally we can map! as you see a bit of extra work was needed, but you ultimately have more control.
```{r}
num_bins <- 7
cmap <- brewer.pal(num_bins, "Blues")
ggplot() +
geom_sf(data = spf_filtered, aes(fill = e_breaks), color = NA) +
scale_fill_manual(
values = cmap,
name = "SPF (cloud probability)", # Improved legend title
labels = gsub("[,]", "-", paste0(gsub("[\\[\\]()]", " ", levels(spf_filtered$e_breaks), perl = TRUE))) # Replace comma with hyphen and remove brackets/parentheses from labels
) +
labs(
title = "SPF 2025: Equal Interval Classification",
fill = NULL # Remove the fill label
) +
theme_void() +
theme(legend.position = "bottom")
```
It is important to understand that **equal intervals** can first and foremost be visualised on the data distribution. We have already created these intervals with the function `classIntervals` in Step 1 of the `ggplot` code above. Here we need a couple of extra steps to collect the break values and plot them in histogram form.
::: callout-tip
## Class Interval styles
The function `classIntervals` has the following styles: "fixed", "sd", "equal", "pretty", "quantile", "kmeans", "hclust", "bclust", "fisher", "jenks", "dpih", "headtails", "maximum", or "box".
:::
```{r}
#| warning: false
# Same step as above
e_breaks <- classIntervals(spf_filtered$cloud_probability, n = 7, style = "equal")
spf_filtered$e_breaks <- cut(spf_filtered$cloud_probability, e_breaks$brks)
# Collect the values of the breaks
e_break_values <- e_breaks$brks
# Place the values in a dataframe
e_break_values_df <- data.frame(BreakValues = e_break_values)
# Create a ggplot2 visualization with 'spf_filtered' dataset as the data source
# and 'cloud_probability' as the variable for the x-axis.
ggplot(spf_filtered, aes(x = cloud_probability)) +
# Add a density plot to the visualization with fill color set to dark blue
# and transparency (alpha) set to 0.4.
geom_density(fill = "darkblue", alpha = 0.4) +
# Add a rug plot (small tick marks) along the x-axis with transparency (alpha) set to 0.5.
geom_rug(alpha = 0.5) +
# Add vertical lines to the plot based on the 'e_break_values_df' dataset
# with x-intercepts specified by the 'BreakValues' variable.
# The color of these lines is set to orange.
geom_vline(data = e_break_values_df, aes(xintercept = BreakValues), color = "darkorange") +
# Apply the 'theme_minimal()' theme to the plot for a minimalistic appearance.
theme_minimal() +
# Modify the x-axis label.
labs(x = "SPF (cloud probability)")
```
Technically speaking, the figure is created by overlaying a KDE plot with vertical bars for each of the break points. This makes much more explicit the issue highlighted by which the middle bins contain a large amount of observations while the ones at either extreme only encompass a handful of them.
## Quantiles
One solution to obtain a more balanced classification scheme is using quantiles. This, by definition, assigns the same amount of values to each bin: the entire series is laid out in order and break points are assigned in a way that leaves exactly the same amount of observations between each of them. This "observation-based" approach contrasts with the "value-based" method of equal intervals and, although it can obscure the magnitude of extreme values, it can be more informative in cases with skewed distributions.
The code required to create the choropleth mirrors that needed above for equal intervals. As before, first we create the intervals, in this case quantiles.
```{r}
#| warning: false
# Find quantile breaks for data segmentation into four groups.
qt_breaks <- classIntervals(spf_filtered$cloud_probability, n = 4, style = "quantile")
# Assign the class breaks to the data
spf_filtered$qt_breaks <- cut(spf_filtered$cloud_probability, qt_breaks$brks)
```
Then we map the data:
```{r}
num_bins <- 4
# Define a color palette for visualizing data.
cmap <- brewer.pal(num_bins, "Blues")
# plot
ggplot() +
geom_sf(data = spf_filtered, aes(fill = qt_breaks), color = NA) +
theme_void() + # remove x and y axis
scale_fill_manual(
values = cmap,
name = "SPF (cloud probability)", # Improved legend title
labels = gsub("[,]", "-", paste0(gsub("[\\[\\]()]", " ", levels(spf_filtered$qt_breaks), perl = TRUE)))) + # Replace comma with hyphen and remove brackets/parentheses from labels
labs(
title = "SPF 2025: Quantile Classification",
fill = NULL # Remove the fill label
) +
theme_void() +
theme(legend.position = "bottom")
```
As we are dealing with small areas across the whole UK, it is easier to see how the data is being divided differently in the histogram.
```{r}
#| warning: false
qt_breaks <- classIntervals(spf_filtered$cloud_probability, n = 4, style = "quantile")
spf_filtered$qt_breaks <- cut(spf_filtered$cloud_probability, qt_breaks$brks)
# Collect the values of the breaks
qt_break_values <- qt_breaks$brks
# Place the values in a dataframe
qt_break_values_df <- data.frame(BreakValues = qt_break_values)
# Create a ggplot2 visualization
ggplot(spf_filtered, aes(x = cloud_probability)) +
# Density plot
geom_density(fill = "darkblue", alpha = 0.4) +
# Add a rug plot (small tick marks)
geom_rug(alpha = 0.5) +
# Add vertical lines at 'BreakValues'
geom_vline(data = qt_break_values_df, aes(xintercept = BreakValues), color = "darkorange") +
theme_minimal() +
labs(x = "SPF (cloud probability)")
```
## Fisher-Jenks
Equal interval and quantiles are only two examples of very many classification schemes to encode values into colors. As an example of a more sophisticated one, let us create a Fisher-Jenks choropleth. As before, first we create the intervals, in this case fisher jenks
```{r}
#| warning: false
# Find fisher breaks for data segmentation into 7 groups.
fish_breaks <- classIntervals(spf_filtered$cloud_probability, n = 7, style = "fisher")
# Assign the class breaks to the data
spf_filtered$fish_breaks <- cut(spf_filtered$cloud_probability, fish_breaks$brks)
```
Then we map the data:
```{r}
num_bins <- 7
# Define a color palette for visualizing data.
cmap <- brewer.pal(num_bins, "Blues")
# plot
ggplot() +
geom_sf(data = spf_filtered, aes(fill = fish_breaks), color = NA) +
theme_void() + # remove x and y axis
scale_fill_manual(
values = cmap,
name = "SPF (cloud probability)", # Improved legend title
labels = gsub("[,]", "-", paste0(gsub("[\\[\\]()]", " ", levels(spf_filtered$fish_breaks), perl = TRUE)))) + # Replace comma with hyphen and remove brackets/parentheses from labels
labs(
title = "SPF 2025: Fisher-Jenks Classification",
fill = NULL # Remove the fill label
) +
theme_void() +
theme(legend.position = "bottom")
```
Now let's look at the density plot
```{r}
#| warning: false
fish_breaks <- classIntervals(spf_filtered$cloud_probability, n = 7, style = "fisher")
spf_filtered$fish_breaks <- cut(spf_filtered$cloud_probability, fish_breaks$brks)
# Collect the values of the breaks
fish_break_values <- fish_breaks$brks
# Place the values in a dataframe
fish_break_values_df <- data.frame(BreakValues = fish_break_values)
# Create a ggplot2 visualization
ggplot(spf_filtered, aes(x = cloud_probability)) +
# Density plot
geom_density(fill = "darkblue", alpha = 0.4) +
# Add a rug plot (small tick marks)
geom_rug(alpha = 0.5) +
# Add vertical lines at 'BreakValues'
geom_vline(data = fish_break_values_df, aes(xintercept = BreakValues), color = "darkorange") +
theme_minimal() +
labs(x = "SPF (cloud probability)")
```
For example, the bins at the extremes of the distribution cover a much wider span than those in the middle, because there are fewer small areas in those value ranges.
You will notice a lot cooler difference once you play around with a larger dataset.
## Zooming into the map
A general map of an entire region, or urban area, can sometimes obscure local patterns because they happen at a much smaller scale that cannot be perceived in the global view. One way to solve this is by providing a focus of a smaller part of the map in a separate figure. Although there are many ways to do this in `R`, the most straightforward one is to define the bounding box.
As an example, let us zoom into Wales — which stands out on the SPF map above with generally higher cloud probability than the rest of the UK. Rather than hardcoding coordinates (which depend on knowing the exact CRS of the data), we compute the bounding box directly from a `nation` subset of the data:
### Zoom into full map
```{r}
# Get the bounding box of Wales directly from the data --
# this works regardless of which CRS the file happens to be in
wales_bbox <- spf_25 %>%
filter(nation == "Wales") %>%
st_bbox()
```
We use the function `coord_sf` to zoom at the desired level.
```{r}
ggplot(data = spf_25) +
geom_sf(aes(fill = cloud_probability), color = NA) +
scale_fill_distiller(palette = "Blues", direction = 1) +
theme_void() +
coord_sf(xlim = c(wales_bbox["xmin"], wales_bbox["xmax"]),
ylim = c(wales_bbox["ymin"], wales_bbox["ymax"])) +
labs(fill = "SPF", title = "SPF 2025: Zoomed into Wales")
```
## Putting it all together: a publication-ready map
The four maps above were deliberately rough — the point was to show the classification mechanics, not to produce a finished product. A map intended for a report or publication needs a few more cartographic elements: a clear title, a colour scheme chosen on purpose, a north arrow, a scale bar, and a source credit.
For the colour scheme, we'll switch to `viridis` — a perceptually uniform, colourblind-safe palette, and a deliberate change from the `Blues` used throughout this lab, to show that the choice of palette is independent of the choice of classification. We'll also use `scale_fill_viridis_c()` for a continuous gradient fill rather than manually cutting the data into bins, which is a quicker route to a clean map when you don't need the fine control of a manual classification.
For the north arrow and scale bar, we'll use the [`ggspatial`](https://paleolimbot.github.io/ggspatial/) package's `annotation_north_arrow()` and `annotation_scale()` — both are CRS-aware, meaning they work out the correct orientation and distance regardless of what projection `spf_filtered` happens to be in, the same "don't hardcode it, derive it from the data" principle we used for the zoom above.
```{r}
#| message: false
library(ggspatial)
ggplot(data = spf_filtered) +
geom_sf(aes(fill = cloud_probability), color = NA) + # continuous fill -- no manual breaks here, a full gradient rather than discrete bins
scale_fill_viridis_c(name = "SPF (cloud probability)", direction = -1) +
annotation_north_arrow(location = "tr", which_north = "true", style = north_arrow_fancy_orienteering) +
annotation_scale(location = "br", width_hint = 0.3) +
labs(
title = "Sun Probability Framework, UK, 2025",
caption = "Source: Imago Data Service, Sun Probability Framework 2025"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(face = "bold", size = 16),
plot.caption = element_text(size = 8, color = "grey40")
)
```
::: callout-tip
`ggspatial` isn't installed by default — `install.packages("ggspatial")` if you don't have it. `which_north = "true"` points the arrow at true north from wherever it sits on the map, rather than just "up" (`"grid"`), which matters more the further you get from a map's central meridian.
:::
## Additional resources
- On Drawing beautiful maps with `sf` and `ggplot` see [here](https://r-spatial.org/r/2018/10/25/ggplot2-sf.html)
- If you want to have a look at **Choropleths in Python** have a look at the chapter on choropleth mapping by [Rey, Arribas-Bel and Wolf](https://geographicdata.science/book/notebooks/05_choropleth.html)
::: callout-tip
## Want to go further with Imago data?
The SPF dataset we used in this lab is one of several openly available products from [**Imago**](https://imago.ac.uk), the imagery data service for sustainability, prosperity and wellbeing, part of [Smart Data Research UK](https://www.sdruk.ukri.org/). They publish their own training materials, free and openly licensed:
- [**Imago training**](https://imago-sdruk.github.io/Imago_training/) — the full set, covering air temperature, precipitation and SPF, in both `R` and `Python`. Several come with notebooks you can run straight in your browser, no installation needed.
- [**From grids to areas**](https://imago-sdruk.github.io/Imago_training/satellites_to_areas.html) — start here if you're curious where a dataset like SPF actually *comes from*. It walks through how raw satellite imagery becomes neighbourhood-level statistics: spectral bands, resolution and revisit cycles, then the four-stage pipeline from noisy cloud-covered grid → cleaned composite → modelled indicator → values aggregated to LSOAs.
- [**SPF: Sun Probability Framework**](https://imago-sdruk.github.io/Imago_training/spf.html) and [**Comparing 2024 and 2025**](https://imago-sdruk.github.io/Imago_training/spf_1.html) — more on the specific dataset you've just been mapping, including a worked example across multiple years.
Worth a look if you're thinking about using Earth observation data in your own assignment or dissertation.
:::