Creating Heatmaps in R

Heatmaps are a popular way to represent how data is distributed across space. You’ll see them used to visualize population density, disease burden, service coverage, crime patterns, and many other spatial patterns, because they make “where things are concentrated” easy to grasp at a glance.

In this tutorial, we’ll build a simple population density heatmap in R using Uganda’s most recent population census outputs from the Uganda Bureau of Statistics. We’ll pull district boundaries via the rgeoboundaries package, then combine those boundaries with UBOS district-level totals to create a grid-based heatmap for five districts in Eastern Uganda: Mbale, Sironko, Bududa, Manafwa, and Tororo.

Note

This work is inspired by the broader R community, especially R-Ladies chapters that regularly publish approachable mapping and spatial-data learning materials. If you’re looking for more spatial learning pathways, the R-Ladies Ottawa community is a great place to start, and the community-maintained R Spatial Ladies resources list is packed with references and tools.

As we began working with the data, we ran into a limitation, which we of course knew before hand, our dataset contains only one population value per district, yet we wanted a heatmap that feels more continuous and visually informative. With district totals alone, we cannot observe true within-district variation, we can only represent an overall total (or an average once we convert totals into density).

For this tutorial, we use a practical workaround: we compute population density per district and then assign that same density value to every grid cell inside the district boundary. This produces a heatmap-like surface while staying consistent with what the data can support.

In cases where more detailed data are available such as sub-county or parish level counts, we can move beyond uniform assignment and generate a much more realistic surface, where density varies within districts in a way that better reflects how population is actually distributed.

Overview

We shall follow the steps below:

1. Create a regular grid of cells over the mapped area.

2. Clip the grid to district boundaries.

3. Calculate population density for each district.

4. Apply that density uniformly to all grid cells within each district

Let’s break down each step.

Step 1: Load Libraries

library(tidyverse)      # for data manipulation
library(sf)             # for spatial operations
library(scales)         # for formatting numbers
library(rgeoboundaries) # for getting administrative boundaries

Step 2: Get the Data

First, we need district boundaries and population data for the selected districts:

# population data
ugpop <- readxl::read_xlsx('ug_population_2024.xlsx')

# get district-level boundaries for uganda
uganda_districts <- geoboundaries(country = "Uganda", adm_lvl = 3)

# filter for our 5 districts of interest
selected_districts <- uganda_districts |>
  filter(shapeName %in% c("Mbale", "Sironko", "Bududa", "Manafwa", "Tororo"))

# join with population data
selected_districts_pop <- selected_districts |>
  left_join(ugpop, by = c("shapeName" = "district"))

Step 3: Create a the Grid

The st_make_grid() function creates a regular grid of square cells. By default, it creates a 10×10 grid.

# create a regular grid (10×10 cells)
grid <- st_make_grid(selected_districts)

ggplot() +
  geom_sf(data = selected_districts) +
  geom_sf(
    data = grid,
    alpha = 0.2) +
  theme_void()

However, to make our heatmap, we’d like to have more cells. So let’s adjust our code to make a 100 by 100 grid:

# grid with 100x100 cells
grid <- st_make_grid(selected_districts, n = c(100, 100))

ggplot() +
  geom_sf(data = selected_districts) +
  geom_sf(data = grid, alpha = 0.2) +
  theme_void()

This creates 10,000 small square cells covering the entire bounding box of our five districts.

Step 4: Clip the Grid to District Boundaries

When we create a regular grid over an area, many grid cells will inevitably fall outside the districts we’re interested in (especially along edges). If we mapped those cells as-is, we’d end up visualizing values in places that are not part of our map. The solution is to clip (trim) the grid so we keep only the portions that fall inside our selected district boundaries.

st_intersection() overlays the district polygons and the grid, returning only the grid pieces that intersect the districts (i.e., the grid clipped to the district shapes), st_make_valid() fixes any invalid geometries that can cause spatial operations to fail. This is especially useful in Uganda’s administrative boundaries where complex shapes may include geometry issues that trigger errors during intersection. mutate(grid_id = row_number()) creates a unique ID for each resulting grid cell/piece. We’ll use this grid_id later when joining attributes (like density values) back onto the grid.

# clip grid to district boundaries
grid_clip <- selected_districts |>
  st_intersection(grid) |>
  st_make_valid() |>
  mutate(grid_id = row_number()) |>
  dplyr::select(grid_id)


ggplot() +
  geom_sf(data = selected_districts) +
  geom_sf(
    data = grid_clip,
    alpha = 0.2) +
  theme_void()

After trimming the grid to the district boundaries, we can plot the result to confirm everything worked as expected. The clipped grid now contains only the cells that fall inside the study area, a total of 5,097 grid cells within the five selected districts give us a clean base for assigning and visualizing population density.

Step 5: Calculate Population Density

Now we compute population density (people per km²) for each district and attach that value to the clipped grid. This step matters because sf::st_area() returns area in square meters, but population density is usually reported in people per square kilometer so we must convert units.

districts_grid_pop <- grid_clip |>
  st_join(selected_districts_pop) |>
  group_by(shapeName) |>
  mutate(
    # calculate total district area in km² (not m²)
    district_area_km2 = sum(as.numeric(st_area(geometry))) / 1e6,
  
    # calculate district density (people per km²)
    pop_density = total[1] / district_area_km2
  ) |>
  ungroup() |>
  group_by(grid_id) |>
  summarize(
    pop_density = mean(pop_density, na.rm = TRUE),
    .groups = "drop"
  )

Step 6: Create the Visualization

ggplot() +
  # district boundaries
  geom_sf(
    data = selected_districts,
    fill = NA,
    color = "gray20",
    linewidth = 0.8) +
  
  # heatmap
  geom_sf(
    data = districts_grid_pop,
    aes(fill = pop_density),
    color = NA) +
  
  # color gradient
  scale_fill_gradient(
    low = "#FFF7BC",
    high = "#D7301F",
    na.value = "gray95",
    labels = comma,
    name = "Population Density (people per km²)",
    guide = guide_colorbar(
      barwidth = 15,
      barheight = 0.8,
      title.position = "top",
      title.hjust = 0.5)) +
  
  # district labels
  geom_sf_text(
    data = selected_districts,
    aes(label = shapeName),
    size = 3.5,
    fontface = "bold",
    color = "gray20") +
  
  # styling
  labs(
    title = "Population Density Heatmap",
    subtitle = "Eastern Uganda Districts") +
  theme_void() +
  theme(
    plot.title = element_text(
      size = 16,
      face = "bold",
      hjust = 0.5,
      margin = margin(t = 10, b = 5)),
    plot.subtitle = element_text(
      size = 10,
      hjust = 0.5,
      color = "gray30",
      margin = margin(b = 15)),
    plot.caption = element_text(
      size = 8,
      color = "gray50",
      hjust = 1,
      margin = margin(t = 10)),
    legend.position = "bottom",
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9),
    plot.background = element_rect(fill = "white", color = NA),
    plot.margin = margin(15, 15, 15, 15))

The final map produces a clean, heatmap-style view of population density across the five districts. Each grid cell is filled using the district’s computed density value (people per km²).

While we have chosen a few districts, you can use this same process to make heatmaps for other data and for different locations, depending on the map data that you have. With more data availability, we can further improve the analysis and produce more realistic density patterns.

Good luck making some interesting heatmaps for your work, and feel free to build on this workflow as shared.