library(tmap)
library(tmap.mapgl)
library(sf)
library(terra)
library(dplyr)The World dataset contains the variable HPI
(Happy Planet Index), a measure of sustainable wellbeing combining life
expectancy, wellbeing, inequality, and ecological footprint.
"view" mode. Make a proportional symbol map
of World where:
fill colour encodes HPIsize encodes pop_est (population
estimate)tmap_mode("view")
data("World")
tm_shape(World) +
tm_symbols(fill = "HPI", size = "pop_est")size.scale = tm_scale_continuous(values.scale = ...)fill.scale = tm_scale_intervals(values = "...")cols4all::c4a_gui() — look at
the diverging palettes. HPI ranges from low (bad) to high (good), so a
green–red or purple–green diverging palette works well.# your code here — build on (a)hover argument.# your code here — build on (b)?World to see all available variables. Choose a
selection that you think is informative for popups and pass them to
popup.vars inside tm_symbols(). Which
variables did you choose and why?# your code here — build on (c)
# hint: popup.vars = c("name", "HPI", "life_exp", ...)Switch to "maplibre" mode using
tmap_mode("maplibre") and reproduce the map from 1.1d.
You may need to increase values.scale compared to view
mode — MapLibre renders symbols at a different scale. Adjust until the
bubble sizes feel right.
Also note: in MapLibre’s globe projection, bubbles are flat circles
projected onto a sphere, which can make them appear to “float” slightly
above the surface when rotating the globe. The standard
size.scale maps data values to circle area, which
is perceptually correct for flat maps. On a globe, mapping to sphere
volume instead gives better results — replace
values.scale = ... with
values = tm_seq(0, 1, power = 1/3) inside
tm_scale_continuous().
tmap_mode("maplibre")
# your code here — try adjusting values.scale, then try values = tm_seq(0, 1, power = 1/3)tmap_providers() to inspect the available basemaps.
For a diverging colour palette, high contrast backgrounds work best —
pure white or black, rather than any grey. Try
"ofm.bright". Which basemap do you prefer and why?tmap_providers()
tm_basemap("ofm.bright") +
# your map hereNote: Some basemap providers require an API key. The
tmap_providers()output indicates which ones are freely available.
In this exercise you will download the 2025 global population raster from WorldPop, crop it to a country or continent, aggregate it to a workable resolution, and visualise it as a 3D extruded map in MapLibre mode.
About the raster cell size
The WorldPop data are stored in a geographic coordinate system (WGS 84 longitude/latitude). This means all cells have the same angular size — approximately 0.0083° × 0.0083° at 1 km resolution — but their physical size varies with latitude. Cell height (north–south) is constant at roughly 59 km after aggregation by a factor of 64. Cell width (east–west) shrinks as you move away from the equator, because lines of longitude converge toward the poles: at the equator a cell is ~59 km wide, while near the poles it can be as narrow as ~6 km. This is an inherent property of geographic rasters, not a data artefact.
In practice, when you aggregate by a fixed factor the resulting cells represent a larger physical area near the equator than near the poles. Keep this in mind when comparing population counts across different latitudes.
The WorldPop 2025 global 1km constrained population dataset is a single global file (~280 MB). The code below downloads it once and caches it locally.
worldpop_url <- paste0(
"https://data.worldpop.org/GIS/Population/",
"Global_2015_2030/R2025A/2025/0_Mosaicked/v1/1km_ua/constrained/",
"global_pop_2025_CN_1km_R2025A_UA_v1.tif"
)
local_file <- "data/global_pop_2025_1km.tif"
dir.create("data", showWarnings = FALSE)
if (!file.exists(local_file)) {
message("Downloading WorldPop 2025 (~280 MB) — this will take a while...")
options(timeout = 2000)
download.file(worldpop_url, destfile = local_file, mode = "wb")
}
pop_global <- rast(local_file)
pop_globalUse the World sf dataset to select a boundary and crop
the global raster to it. Because the raster is at 1km resolution, very
small countries may not crop well — larger countries work best.
A note on border cells: The
Worldboundaries are simplified polygons. Where a grid cell straddles a country border,terra::mask()decides by cell centre: if the centre falls inside the polygon the whole cell is kept; if outside, the whole cell is set toNA. This means border cells may contain population from the neighbouring country, or conversely some border cells from within the country may be dropped. The effect is small for large areas but can produce a slightly ragged edge at the boundary.
data("World")
# --- Option A: Germany ---
area <- World |> filter(name == "Germany")
# --- Option B: Africa ---
# area <- World |> filter(continent == "Africa") |> st_union()
area_reproj <- st_transform(area, crs(pop_global))
pop_cropped <- crop(pop_global, area_reproj) |> mask(area_reproj)
names(pop_cropped) <- "pop"
pop_cropped#' Aggregate a population raster by a given factor
#'
#' @param r SpatRaster with population counts
#' @param fact Aggregation factor (e.g. 8 = 8×8 cells → ~8km resolution)
#' @return Aggregated SpatRaster with layer named "pop"
aggregate_pop <- function(r, fact = 8) {
r_agg <- terra::aggregate(r, fact = fact, fun = sum, na.rm = TRUE)
names(r_agg) <- "pop"
r_agg
}
# Suggested factors:
# Germany: fact = 8 → ~8km cells
# Africa: fact = 64 → ~32km cells
pop_agg <- aggregate_pop(pop_cropped, fact = 8)
pop_aggNote on CRS: For large areas like Germany or Africa the default globe CRS works well and has no artefacts.
tm_crs(3857)(Web Mercator) is only needed for small regions like the Netherlands where 3D polygon artefacts can appear.
tm_maplibre(pitch = 45) to tilt the view — this already
gives a sense of the 3D effect we will add next.tmap_mode("maplibre")
tm_shape(pop_agg) +
tm_raster(col = "pop",
col.scale = tm_scale_continuous(values = "brewer.yl_or_rd"),
col.legend = tm_legend("Population (2025)")) +
tm_maplibre(pitch = 45) +tm_polygons_3d().
Set fill = "pop" for colour and height = "pop"
for extrusion height. Experiment with values.scale in
height.scale and with
height.max/height.min in
opt_tm_polygons_3d().tm_shape(pop_agg) +
tm_polygons_3d(
fill = "pop",
fill.scale = tm_scale_continuous(values = "brewer.yl_or_rd"),
height = "pop",
height.scale = tm_scale_continuous(values.scale = ___),
options = opt_tm_polygons_3d(height.max = "10%", height.min = "0.1%")
) +
tm_maplibre(pitch = 45) +# Option A: Germany
area2 <- World |> filter(name == "Germany")
# Option B: Africa
# area2 <- World |> filter(continent == "Africa") |> st_union()
area2_reproj <- st_transform(area2, crs(pop_global))
pop_cropped2 <- crop(pop_global, area2_reproj) |> mask(area2_reproj)
names(pop_cropped2) <- "pop"
pop_agg2 <- aggregate_pop(pop_cropped2, fact = ___) # 8 for Germany, 64 for Africa
tm_shape(pop_agg2) +
tm_polygons_3d(
fill = "pop",
fill.scale = tm_scale_continuous(values = "brewer.yl_or_rd"),
height = "pop",
height.scale = tm_scale_continuous(values.scale = ___),
options = opt_tm_polygons_3d(height.max = "10%", height.min = "0.1%")
) +
tm_maplibre(pitch = 45) +