Required packages

library(tmap)
library(tmap.mapgl)
library(sf)
library(terra)
library(dplyr)

Exercise 1 — World data: view mode vs MapLibre

1.1 — Proportional symbol map

The World dataset contains the variable HPI (Happy Planet Index), a measure of sustainable wellbeing combining life expectancy, wellbeing, inequality, and ecological footprint.

  1. Set tmap to "view" mode. Make a proportional symbol map of World where:
    • fill colour encodes HPI
    • size encodes pop_est (population estimate)
tmap_mode("view")
data("World")

tm_shape(World) +
  tm_symbols(fill = "HPI", size = "pop_est")
  1. Can you improve the map? Try enlarging the bubbles and applying a diverging colour palette for HPI.
tm_shape(World) +
  tm_symbols(fill = "HPI", size = "pop_est",
             fill.scale = tm_scale_intervals(values = "pu_gn"),
             size.scale = tm_scale_continuous(values.scale = 3))

values.scale = 3 multiplies the default symbol sizes by 3. The "pu_gn" palette is a purple–green diverging scheme — low HPI maps to purple, high HPI to green. Alternatives can be explored via cols4all::c4a_gui() under the diverging tab.

  1. Add a mouse hover label showing the country name.
tm_shape(World) +
  tm_symbols(fill = "HPI", size = "pop_est",
             fill.scale = tm_scale_intervals(values = "pu_gn"),
             size.scale = tm_scale_continuous(values.scale = 3),
             hover = "name")
  1. Choose informative popup variables and pass them to popup.vars. Which did you choose and why?
tm_shape(World) +
  tm_symbols(
    fill       = "HPI",
    size       = "pop_est",
    fill.scale = tm_scale_intervals(values = "pu_gn"),
    size.scale = tm_scale_continuous(values.scale = 3),
    hover      = "name",
    popup.vars = c("name", "HPI", "life_exp", "well_being",
                   "footprint", "gdp_cap_est", "pop_est")
  )

life_exp and well_being are HPI components; footprint is the ecological cost; gdp_cap_est provides economic context. Including name keeps the popup self-contained.

1.2 — Switch to MapLibre

  1. Switch to "maplibre" mode and reproduce the map from 1.1d.

    You may need to increase values.scale — MapLibre renders symbols at a different scale than view mode. Also note: in the globe projection, bubbles are flat circles projected onto a sphere, which can make them appear to “float” when rotating. Mapping to sphere volume rather than circle area gives better results on a globe: replace values.scale = ... with values = tm_seq(0, 1, power = 1/3).

tmap_mode("maplibre")

tm_shape(World) +
  tm_symbols(
    fill       = "HPI",
    size       = "pop_est",
    fill.scale = tm_scale_intervals(values = "pu_gn"),
    size.scale = tm_scale_continuous(values = tm_seq(0, 1, power = 1/3),
                                     values.scale = 6),
    hover      = "name",
    popup.vars = c("name", "HPI", "life_exp", "well_being",
                   "footprint", "gdp_cap_est", "pop_est")
  )

values.scale = 6 is needed in MapLibre because its default symbol scale is smaller than in Leaflet view mode. tm_seq(0, 1, power = 1/3) applies a cube-root transformation to the size mapping, which corresponds to volume scaling for spheres (radius ∝ value^(1/3)), reducing the visual dominance of the very largest bubbles and making the globe feel more natural when rotating.

  1. Use tmap_providers() to inspect the available basemaps. For diverging palettes, high contrast backgrounds work best — pure white or black perform better than any shade of grey. "ofm.bright" (OpenFreeMap) is a clean, freely available bright basemap.
tmap_providers()

tm_basemap("ofm.bright") +
tm_shape(World) +
  tm_symbols(
    fill       = "HPI",
    size       = "pop_est",
    fill.scale = tm_scale_intervals(values = "pu_gn"),
    size.scale = tm_scale_continuous(values = tm_seq(0, 1, power = 1/3),
                                     values.scale = 6),
    hover      = "name",
    popup.vars = c("name", "HPI", "life_exp", "well_being",
                   "footprint", "gdp_cap_est", "pop_est")
  )

A bright white basemap keeps the full range of the purple–green palette visible. Grey basemaps reduce apparent saturation of mid-range values, making it harder to distinguish countries near the HPI median. Dark basemaps have the same problem at the other end of the scale.

Exercise 2 — 3D population density map with WorldPop data

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.

2.1 — Download the global population raster

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_global

2.2 — Crop to a country or continent

Note: very small countries may not crop well at 1km resolution — larger countries work best.

A note on border cells: The World boundaries 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 to NA. 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, and means population totals near borders are not precisely attributed to a single country.

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

2.3 — Aggregate to a workable resolution

aggregate_pop <- function(r, fact = 8) {
  r_agg        <- terra::aggregate(r, fact = fact, fun = sum, na.rm = TRUE)
  names(r_agg) <- "pop"
  r_agg
}

# Germany: fact = 8  → ~8km cells
# Africa:  fact = 64 → ~64km cells
pop_agg <- aggregate_pop(pop_cropped, fact = 8)
pop_agg

2.4 — Visualise as a 3D extruded map

Note 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.

  1. Flat 2D check map. tm_maplibre(pitch = 45) already tilts the view to give a sense of perspective before adding extrusion.
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) +
  1. 3D extruded map.
tm_shape(pop_agg) +
  tm_polygons_3d(
    fill         = "pop",
    fill.scale   = tm_scale_continuous(values = "brewer.yl_or_rd"),
    fill.legend  = tm_legend("Population (2025)"),
    height       = "pop",
    height.scale = tm_scale_continuous(values.scale = 3),
    options      = opt_tm_polygons_3d(height.max = "10%", height.min = "0.1%")
  ) +
  tm_maplibre(pitch = 45) +

height.max = "10%" caps the tallest bar at 10% of the viewport, preventing dense coastal cities from dwarfing everything else. tm_maplibre(pitch = 45) sets the initial tilt so the 3D effect is immediately visible.

  1. Compare Germany and Africa.
# 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 = 64)

tm_shape(pop_agg2) +
  tm_polygons_3d(
    fill         = "pop",
    fill.scale   = tm_scale_continuous(values = "brewer.yl_or_rd"),
    fill.legend  = tm_legend("Population (2025)"),
    height       = "pop",
    height.scale = tm_scale_continuous(values.scale = 3),
    options      = opt_tm_polygons_3d(height.max = "10%", height.min = "0.1%")
  ) +
  tm_maplibre(pitch = 45) +

Germany shows a polycentric pattern — recognisable peaks at Berlin, Hamburg, Munich, Frankfurt, Cologne, and the Rhine-Ruhr conurbation. Africa at the continental scale reveals very different urban structures: dense peaks around Lagos and the Niger delta, the Nile delta, Kinshasa, and the East African highlands, with large sparsely populated stretches of the Sahara and Kalahari in between.