Working with Geodata in R

E1493 - Data Journalism

Author

Simon Munzert

Published

March 23, 2026

NoteBefore you start

Install the packages used in this workshop:

install.packages(c(
  "tidyverse", "remotes", "sf", "terra", "stars",
  "rnaturalearth", "rnaturalearthdata",
  "nominatimlite", "tidygeocoder", "osmdata",
  "patchwork", "ggrepel", "ggspatial", "prettymapr", "leaflet",
  "ggimage", "elevatr", "tidyterra", "hexbin",
  "mapview", "mapedit", "mapdeck"
))
remotes::install_github("ropensci/rnaturalearthhires")

This script accompanies the Session 6 slides on Mapping and Geodata. It’s designed as a hands-on complement — the slides cover concepts like spatial data types, projections, and distance; this script lets you work with real spatial data in R. We use the German breweries dataset throughout as a running example.

The script builds on concepts from Pebesma & Bivand (2023), Spatial Data Science: With Applications in R, freely available online. If you want to go deeper, that book is the definitive reference for sf and stars.

The running example: Franconian brewery density

For those of you who took Introduction to Data Science, this dataset will feel familiar. We are going back to the instructor’s Franconian roots — a region in northern Bavaria famous for its extraordinarily dense network of small breweries. Back then, we had scraped brewery locations from biermap24.de.

The central question we explore throughout this script is: does the data actually confirm that Franconia is the region with the highest brewery density in Germany? It is one thing to claim this as local pride; it is another to test it with spatial data.

We will start by loading administrative boundaries, geocode brewery locations, perform spatial joins to assign each brewery to a federal state, compute density metrics, and finally visualise the results in a range of map styles. By the end, we will have a quantitative answer.

Part 1: Importing Spatial Data

The two main types of spatial data — vector (points, lines, polygons) and raster (gridded pixels) — each have their own packages and workflows in R. The sf package handles vector data; terra (and stars) handle raster data. Both integrate with the tidyverse.

Vector data formats

Before diving into code, it helps to know the landscape of common spatial file formats you will encounter:

Common spatial data formats
Format Extension Pros Cons
Shapefile .shp + 3–5 sidecar files Universal legacy standard; accepted everywhere Multiple files must stay together; 10-character column name limit; 2 GB file size cap
GeoJSON .geojson or .json Human-readable text; web-native; single file Large file sizes for dense geometry; not ideal for big datasets
GeoPackage .gpkg Modern single-file format; no column name limits; large files OK; open standard Less familiar to non-GIS users
KML/KMZ .kml / .kmz Google Earth native; easy to share Limited attribute support; verbose XML
CSV with coords .csv Dead simple; works with any tool No CRS metadata; no native geometry type support
TipShapefile vs. GeoPackage

Shapefiles are the legacy standard — they consist of 4–6 files (.shp, .dbf, .shx, .prj, …) that must travel together. GeoPackage (.gpkg) is the modern replacement: a single file, no column name truncation, no file size limits. Use GeoPackage when you can.

Vector data with sf

The workhorse function is st_read(), which can read shapefiles, GeoJSON, GeoPackage, KML, and dozens of other formats via GDAL. Let’s start by loading a polygon map of Germany from the rnaturalearth package:

library(sf)
library(rnaturalearth)

# Country boundary as an sf object
germany <- ne_countries(
  scale = "medium",
  country = "Germany",
  returnclass = "sf"
)

# What does an sf object look like?
class(germany)
[1] "sf"         "data.frame"

An sf object is just a data frame with a special geometry column. This means all your familiar dplyr verbs work on it — filter(), mutate(), select(), left_join(), you name it. The geometry just comes along for the ride.

# The geometry column
st_geometry(germany)
Geometry set for 1 feature 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 5.85752 ymin: 47.27881 xmax: 15.0166 ymax: 55.05874
Geodetic CRS:  WGS 84
# CRS (coordinate reference system)
st_crs(germany)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

The CRS tells you how coordinates map onto the Earth’s surface. Here we have WGS 84 (EPSG:4326) — the standard longitude/latitude system used by GPS.

Reading files from disk

In real projects, you’ll often read shapefiles or GeoJSON files. The syntax is the same regardless of format — st_read() detects the driver automatically:

# Reading a shapefile (hypothetical example)
# districts <- st_read("data/german_districts.shp")

# Reading GeoJSON
# districts <- st_read("data/german_districts.geojson")

# Reading GeoPackage (modern, recommended format)
# districts <- st_read("data/german_districts.gpkg")

Getting German state boundaries

For our brewery analysis, we’ll need state-level polygons. Let’s get them:

# German federal states (Bundesländer)
germany_states <- ne_states(country = "Germany", returnclass = "sf")

ggplot(germany_states) +
  geom_sf(fill = "gray95", color = "gray40") +
  theme_minimal() +
  labs(title = "German federal states (Bundesländer)")

Modifying and enriching sf objects

Because sf objects are data frames, all standard dplyr verbs work on them. The geometry column is automatically kept or dropped as appropriate.

# (a) mutate() works just like on a plain data frame
# Add a column with the name in upper case
germany_states <- germany_states |>
  mutate(name_upper = toupper(name))

# (b) Compute area of each state with st_area()
germany_states <- germany_states |>
  mutate(area_km2 = as.numeric(st_area(geometry) / 1e6))

germany_states |>
  st_drop_geometry() |>
  select(name, area_km2) |>
  arrange(desc(area_km2))
# (c) Join with an external non-spatial data frame
# Approximate 2023 population figures for German federal states
state_population <- tibble(
  name = c(
    "Bayern", "Niedersachsen", "Baden-Württemberg", "Nordrhein-Westfalen",
    "Brandenburg", "Mecklenburg-Vorpommern", "Sachsen-Anhalt",
    "Thüringen", "Sachsen", "Hessen", "Rheinland-Pfalz",
    "Schleswig-Holstein", "Saarland", "Hamburg", "Berlin",
    "Bremen"
  ),
  population = c(
    13.4, 8.1, 11.2, 18.1,
    2.5, 1.6, 2.2,
    2.1, 4.1, 6.4, 4.2,
    3.0, 1.0, 1.9, 3.7,
    0.7
  )
)

# Left-join: geometry stays in germany_states, population comes in
germany_states <- germany_states |>
  left_join(state_population, by = "name")

# (d) Note: the geometry column comes along automatically — no special handling needed
germany_states |>
  st_drop_geometry() |>
  select(name, area_km2, population) |>
  head(5)
NoteWhere to find spatial data

Good sources for administrative boundaries and other geodata include:

Part 2: Projections

A coordinate reference system (CRS) defines how coordinates map onto the Earth. Getting this right matters — wrong projections lead to distorted maps and incorrect distance/area calculations.

The key distinction: geographic vs. projected

Geographic CRS (e.g., WGS 84 / EPSG:4326): Coordinates are in degrees of longitude and latitude on an ellipsoidal model of the Earth. This is what GPS uses.

Projected CRS (e.g., UTM Zone 32N / EPSG:32632): Coordinates are in meters (or feet) on a flat 2D plane. Every projection introduces some distortion — the question is which property you want to preserve: area, shape, distance, or direction.

When does the choice matter more?

  • Large areas and high latitudes: distortion is small for a map of Berlin but substantial for all of Europe. At 60° N latitude, one degree of longitude is barely half as wide as at the equator.
  • Area comparisons: the Mercator projection (still the default for many web maps) dramatically inflates the apparent size of Greenland, Russia, and Canada relative to countries near the equator. For data journalism stories involving area or density, this is a real problem.
  • Distance and buffer operations: any distance measured in degrees is geometrically meaningless. You need a metric CRS before calling st_buffer() or st_distance().

The following figure illustrates how dramatically projection choice affects the perception of area and shape at the global scale:

world <- ne_countries(scale = "medium", returnclass = "sf")

p_mercator <- ggplot(world) +
  geom_sf(fill = "gray80", color = "white", linewidth = 0.1) +
  coord_sf(crs = "+proj=merc", ylim = c(-15000000, 18000000)) +
  theme_void(base_size = 11) +
  labs(title = "Mercator (EPSG:3857)",
       subtitle = "Distorts area")

p_eqearth <- ggplot(world) +
  geom_sf(fill = "gray80", color = "white", linewidth = 0.1) +
  coord_sf(crs = "+proj=eqearth") +
  theme_void(base_size = 11) +
  labs(title = "Equal Earth",
       subtitle = "Preserves relative areas faithfully")

wrap_elements(p_mercator) + wrap_elements(p_eqearth)

NoteWhere do CRS numbers come from? The EPSG Registry

The EPSG codes (e.g., EPSG:4326, EPSG:25832) come from the EPSG Geodetic Parameter Dataset, originally maintained by the European Petroleum Survey Group (hence EPSG) and now managed by the International Association of Oil & Gas Producers (IOGP). The registry provides a standardised, uniquely numbered catalogue of coordinate reference systems used worldwide.

Key EPSG codes
EPSG Code Name Use case
4326 WGS 84 GPS coordinates; universal standard for data exchange
3857 Web Mercator (Pseudo-Mercator) Google Maps, OpenStreetMap tiles, web mapping
25832 ETRS89 / UTM Zone 32N Official German reference system; metric distances
3035 LAEA Europe Equal-area maps of Europe; statistics, Eurostat

The website epsg.io is the easiest way to look up any code, understand its properties, and find the right projection for your region.

Seeing the world as a globe

An orthographic projection centres the map on a point on the Earth’s surface, giving the visual impression of viewing a globe. It is a striking presentation choice for magazine and news graphics:

# orthographic projection centered on central Europa
ortho_crs <- "+proj=ortho +lat_0=48 +lon_0=10"

# globe boundary
globe_disc <- st_sfc(st_point(c(0, 0)), crs = ortho_crs) |>
  st_buffer(6370000)  # mean Earth radius in metres

ggplot(world) +
    geom_sf(data = globe_disc, fill = "#c8dff0", color = "grey55",  # ocean + rim
          linewidth = 0.6) +
  geom_sf(fill = "lightgreen", color = "black", linewidth = 0.15) +
  coord_sf(crs = ortho_crs) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "navy", color = NA)
    )

Checking and setting CRS

# Check the CRS of our Germany polygon
st_crs(germany)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
# The EPSG code is the shorthand
st_crs(germany)$epsg
[1] 4326

Reprojecting data

Use st_transform() to reproject vector data. For Germany, a good projected CRS is ETRS89 / UTM Zone 32N (EPSG:25832) — this is the official German reference system:

# Reproject to ETRS89 / UTM zone 32N
germany_proj <- st_transform(germany, crs = 25832)

# Now coordinates are in meters, not degrees
st_bbox(germany)       # in degrees
    xmin     ymin     xmax     ymax 
 5.85752 47.27881 15.01660 55.05874 
st_bbox(germany_proj)  # in meters
     xmin      ymin      xmax      ymax 
 279659.3 5236826.5  919717.1 6101489.7 
# Visual comparison
p1 <- ggplot(germany) +
  geom_sf(fill = "gray90") +
  labs(title = "WGS 84 (geographic)") +
  theme_minimal()

p2 <- ggplot(germany_proj) +
  geom_sf(fill = "gray90") +
  labs(title = "ETRS89 / UTM 32N (projected)") +
  theme_minimal()

p1 + p2

The visual difference for Germany is subtle, but the implications for calculations are real: area and distance computations on projected data use planar geometry (fast, exact), while geographic coordinates require spherical/ellipsoidal calculations.

When does the projection matter?

For display (making maps), the projection choice is mostly aesthetic. For analysis (calculating distances, areas, buffers), it’s critical. A common workflow:

  1. Store and exchange data in WGS 84 (EPSG:4326) — the universal standard.
  2. Reproject to an appropriate local CRS before doing any spatial calculations.
  3. Reproject back to WGS 84 for web mapping or data exchange if needed.
# Area of Germany in different CRS
st_area(germany) |> units::set_units("km^2")        # spherical calculation
356166.9 [km^2]
st_area(germany_proj) |> units::set_units("km^2")    # planar calculation
357442.1 [km^2]
ImportantA common trap

If you compute distances or areas on unprojected (lon/lat) data, sf will use spherical geometry by default (via the s2 library), which is correct but slower. If you accidentally treat lon/lat as planar coordinates (e.g., by disabling s2), you’ll get wrong results — especially at higher latitudes where one degree of longitude is much less than 111 km.

Part 3: Transformations and Spatial Overlay

Spatial transformations let you combine datasets to answer real journalistic questions. Here are a few concrete use cases that arise regularly:

  • “Which constituency does this polling station fall in?” — assign each point to its containing polygon.
  • “How many hospitals are within 30 minutes of each neighbourhood?” — buffer + intersect.
  • “Which reporters cover which regions?” — spatial join to assign coverage zones.
  • “Do renewable energy plants cluster near or far from population centres?” — distance calculations and density analysis.

This is where spatial data becomes genuinely powerful for journalism beyond mapping.

Loading our brewery data

Let’s load the geocoded brewery coordinates from the scraping exercise:

# Load the pre-geocoded brewery data
load("data/coords_breweries.RData")
head(coords_df)
# Clean: filter to plausible German coordinates
breweries <- coords_df |>
  filter(
    !is.na(lon), !is.na(lat),
    lon >= 5.5, lon <= 15.5,
    lat >= 47, lat <= 55.5
  )

nrow(breweries)
[1] 656

Converting a data frame to an sf object

The brewery data is just a regular data frame with lon/lat columns. To make it spatial, we use st_as_sf():

breweries_sf <- st_as_sf(
  breweries,
  coords = c("lon", "lat"),
  crs = 4326  # WGS 84
)

# It's now a spatial object
class(breweries_sf)
[1] "sf"         "data.frame"
head(breweries_sf)
# Quick check: plot breweries on top of Germany using projected CRS
ggplot() +
  geom_sf(data = st_transform(germany, 25832), fill = "gray95") +
  geom_sf(data = st_transform(breweries_sf, 25832),
          color = "red", size = 0.4, alpha = 0.5) +
  coord_sf(crs = 25832) +
  theme_minimal(base_size = 14) +
  labs(title = "Breweries in Germany")

Spatial join: which state is each brewery in?

The classic overlay operation: assign each point to the polygon it falls within. This is st_join() — the spatial equivalent of left_join():

# Join breweries to states: each brewery gets the attributes of the state it's in
breweries_with_state <- st_join(breweries_sf, germany_states["name"])

head(breweries_with_state)
# Count breweries per state
breweries_per_state <- breweries_with_state |>
  st_drop_geometry() |>
  count(name, sort = TRUE) |>
  rename(state = name, n_breweries = n)

breweries_per_state

Bavaria comes out on top by a large margin — but how much of this reflects the fact that Bavaria is also the largest German state by area? The introductory question was specifically about density, not raw counts. Franconia (the northern part of Bavaria, covering districts like Oberfranken, Mittelfranken, and Unterfranken) is where the real concentration lies. We will examine this more carefully when we compute brewery density per km² in the exercises.

# Merge counts back to the state polygons
states_with_breweries <- germany_states |>
  left_join(breweries_per_state, by = c("name" = "state")) |>
  mutate(n_breweries = replace_na(n_breweries, 0))

# Use discretised bins for cleaner color mapping
states_with_breweries <- states_with_breweries |>
  mutate(brewery_bin = cut(
    n_breweries,
    breaks = c(-Inf, 10, 50, 100, 200, Inf),
    labels = c("0–10", "11–50", "51–100", "101–200", "201+"),
    right = TRUE
  ))

ggplot(states_with_breweries) +
  geom_sf(aes(fill = brewery_bin), color = "white") +
  scale_fill_brewer(
    palette = "YlOrRd",
    name = "Breweries",
    na.value = "gray90"
  ) +
  coord_sf(crs = 25832) +
  theme_minimal(base_size = 14) +
  labs(
    title = "Number of breweries by German state",
    caption = "Source: biermap24.de"
  )

Buffering, intersection, and other geometric operations

TipWhy you must reproject before buffering

st_buffer() takes a dist argument — but in what unit? If your data is in geographic coordinates (EPSG:4326, degrees of lon/lat), then dist = 50000 means 50,000 degrees, which is meaningless. Even dist = 0.5 (half a degree) produces a buffer that is not a true circle: a half-degree in longitude is much narrower than a half-degree in latitude, especially at higher latitudes like Germany’s (~50° N).

The solution is to reproject to a metric CRS (e.g., EPSG:25832 for Germany) before buffering. Then dist = 50000 means exactly 50,000 metres = 50 km, and the buffer is a geometrically correct circle.

# Reproject for metric operations
breweries_proj <- st_transform(breweries_sf, 25832)
germany_states_proj <- st_transform(germany_states, 25832)

# Buffer: 50 km circle around Berlin
berlin_point <- st_sfc(st_point(c(13.405, 52.52)), crs = 4326) |>
  st_transform(25832)
berlin_buffer <- st_buffer(berlin_point, dist = 50000)  # 50 km in meters

# Which breweries fall within 50 km of Berlin?
breweries_near_berlin <- st_intersection(breweries_proj, berlin_buffer)
cat("Breweries within 50 km of Berlin:", nrow(breweries_near_berlin), "\n")
Breweries within 50 km of Berlin: 6 
# Visualize the buffer
ggplot() +
  geom_sf(data = germany_states_proj, fill = "gray95", color = "gray60") +
  geom_sf(data = berlin_buffer, fill = "steelblue", alpha = 0.2, color = "steelblue") +
  geom_sf(data = breweries_proj, color = "red", size = 0.3, alpha = 0.3) +
  geom_sf(data = breweries_near_berlin, color = "darkred", size = 1.5) +
  coord_sf(
    crs = 25832,
    xlim = st_bbox(berlin_buffer)[c(1, 3)] + c(-100000, 100000),
    ylim = st_bbox(berlin_buffer)[c(2, 4)] + c(-100000, 100000)
  ) +
  theme_minimal(base_size = 14) +
  labs(
    title = "Breweries within 50 km of Berlin",
    subtitle = "Blue circle = 50 km buffer around city center"
  )

Key transformation functions at a glance

Here are the sf functions you’ll use most often for spatial analysis:

# Predicates (return TRUE/FALSE or index lists)
st_intersects(x, y)     # do geometries intersect?
st_within(x, y)         # is x fully within y?
st_contains(x, y)       # does x fully contain y?
st_is_within_distance(x, y, dist)  # within distance threshold?

# Operations (return new geometries)
st_intersection(x, y)   # clip x to y (overlay)
st_difference(x, y)     # x minus y
st_union(x)             # merge all geometries into one
st_buffer(x, dist)      # expand/contract by distance
st_centroid(x)          # center point of polygons

# Measurements
st_distance(x, y)       # pairwise distances
st_area(x)              # polygon areas
st_length(x)            # line lengths

# Joins
st_join(x, y)           # spatial left join (like dplyr::left_join, but by geometry)
st_nearest_feature(x, y) # index of nearest feature in y for each x
NoteFor data journalists

The spatial join (st_join) is your single most useful function. It lets you answer questions like: “Which district is this event in?”, “How many incidents occurred in each neighborhood?”, “Which schools are in the flood zone?” — all without manual matching or lookup tables.

Part 4: Geocoding

Geocoding turns addresses or place names into coordinates. It’s essential when your source data contains location names but no geometry — think lists of schools, hospitals, crime incidents, or… breweries.

Option 2: nominatimlite (lightweight, sf-native)

If you want geocoded results directly as sf objects (instead of lat/lon columns), nominatimlite is a nice lightweight option:

library(nominatimlite)

# Geocode and get sf geometry directly
munich_sf <- geo_lite_sf("München, Germany", limit = 1)
munich_sf
# Geocode multiple places
brewery_cities <- c("Bamberg, Germany", "Kulmbach, Germany", "Bayreuth, Germany")
brewery_cities_sf <- geo_lite_sf(brewery_cities, limit = 1)

  |                                                        
  |                                                  |   0%
  |                                                        
  |=================                                 |  33%
  |                                                        
  |=================================                 |  67%
  |                                                        
  |==================================================| 100%
brewery_cities_sf

Geocoding the brewery dataset

In the scraping-breweries.R script, we geocoded brewery cities using the Nominatim API. Here’s what that workflow looks like with tidygeocoder — a cleaner, more modern approach:

# If you needed to re-geocode the brewery cities from scratch:
# (Don't run this — it takes a while due to rate limits)

# unique_cities <- unique(cities_from_scrape)
#
# geocoded <- tibble(city = unique_cities) |>
#   geocode(
#     address = city,
#     method = "osm",
#     custom_query = list(countrycodes = "de")  # restrict to Germany
#   )
ImportantGeocoding gotchas
  1. Ambiguity. “Frankfurt” could be Frankfurt am Main or Frankfurt an der Oder. Always add country context (e.g., countrycodes = "de").
  2. Rate limits. Nominatim allows 1 request/second. For 500+ cities, that’s 8+ minutes. tidygeocoder handles rate limiting automatically.
  3. Quality varies. Geocoding works well for cities and street addresses, but poorly for vague locations (“near the old bridge”). Always sanity-check results on a map.
  4. Cache your results. Geocoding is slow and API-dependent. Save results to a file (save() or write_csv()) so you don’t re-geocode every time.

Reverse geocoding

Sometimes you have coordinates and need addresses — that’s reverse geocoding:

# What's at these coordinates?
coords <- tibble(
  lat = c(52.5128, 49.4263, 40.6892),
  lon = c(13.3866, 11.1208, -74.0471)
)

coords |>
  reverse_geocode(
    lat = lat, long = lon,
    method = "osm",
    address = "address_found"
  )

Reverse geocoding is useful when you have GPS traces, satellite data, or geotagged photos and need to identify the location in human terms.

Part 5: Plotting Maps

Map design matters in data journalism — a well-crafted map can reveal patterns that tables and bar charts can’t, but a bad map can mislead. This section covers practical map-making in R with an eye toward journalistic standards.

The basic recipe: ggplot2 + geom_sf()

Since sf objects are data frames, they plug directly into ggplot2 via geom_sf(). This is the single most important thing to know about mapping in R:

ggplot() +
  geom_sf(data = st_transform(germany_states, 25832),
          fill = "gray95", color = "gray40") +
  geom_sf(data = st_transform(breweries_sf, 25832),
          color = "#C2001A", size = 0.3, alpha = 0.4) +
  coord_sf(crs = 25832) +
  theme_void(base_size = 14) +
  labs(
    title = "Breweries in Germany",
    caption = "Source: biermap24.de"
  )

Already not bad. But let’s make it better.

Choropleth maps: encoding data in color

The most common journalistic map type. The key decisions are: what variable to show, which color scale, and how to bin the data. Using discrete bins (rather than a continuous gradient) makes the map easier to read and interpret:

states_map <- states_with_breweries |>
  mutate(
    n_breweries = replace_na(n_breweries, 0),
    brewery_bin = cut(
      n_breweries,
      breaks = c(-Inf, 10, 50, 100, 200, Inf),
      labels = c("0–10", "11–50", "51–100", "101–200", "201+"),
      right = TRUE
    )
  )

ggplot(states_map) +
  geom_sf(aes(fill = brewery_bin), color = "white", linewidth = 0.3) +
  scale_fill_brewer(
    palette = "YlOrRd",
    name = "Number of\nbreweries",
    na.value = "gray90"
  ) +
  guides(fill = guide_legend(nrow = 1, title.position = "top")) +
  coord_sf(crs = 25832) +
  theme_void(base_size = 14) +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(color = "gray40")
  ) +
  labs(
    title = "Germany's brewery landscape",
    subtitle = "Bavaria dominates with more than half of all breweries",
    caption = "Source: biermap24.de"
  )

WarningThe choropleth trap

Choropleth maps encode data by area, which means large but sparsely populated regions dominate visually. A map of “total breweries per state” will be dominated by Bavaria partly because Bavaria is big. Always consider normalizing: breweries per capita, per km², or using a different map type entirely.

Kernel density heatmap: understanding spatial concentration

When you have thousands of points, a kernel density map reveals spatial concentration more clearly than individual dots. Instead of every point being drawn independently, the surface shows where points cluster, smoothed into a continuous field. This is particularly useful for understanding the Franconian cluster without getting lost in overlapping dots:

ggplot() +
  geom_sf(data = germany_states, fill = "#f5f5f0", color = "gray70",
          linewidth = 0.2) +
  stat_density_2d(
    data = breweries,
    aes(x = lon, y = lat, fill = after_stat(level)),
    geom = "polygon",
    alpha = 0.7,
    bins = 10,
    h = c(1.2, 1.2)
  ) +
  scale_fill_viridis_c(
    option = "magma",
    direction = -1,
    name = "Density"
  ) +
  coord_sf(crs = 4326, xlim = c(5.5, 15.5), ylim = c(47, 55.5)) +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.margin = margin(10, 10, 10, 10)
  ) +
  labs(
    title = "Brewery density across Germany",
    subtitle = "Kernel density estimate — darker = higher concentration of breweries",
    caption = "Source: biermap24.de"
  )

The Franconian cluster in northern Bavaria shows up strikingly here. Unlike a choropleth aggregated at the state level, the kernel density surface does not hide the sub-state geography: you can see that the high density is concentrated in the northern part of Bavaria, not spread evenly across the whole state.

Combining layers: points + polygons + labels

Journalistic maps often combine multiple layers. Let’s highlight Bavaria and add context:

# Isolate Bavaria
bavaria <- germany_states |> filter(name == "Bayern")

# Major cities for context
major_cities <- tibble(
  name = c("München", "Berlin", "Hamburg", "Bamberg", "Nürnberg"),
  lon = c(11.58, 13.41, 9.99, 10.89, 11.08),
  lat = c(48.14, 52.52, 53.55, 49.89, 49.45)
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# ggrepel needs plain x/y coordinates matching the active CRS (25832).
# Extract projected coords as a plain data frame — do NOT pass sf geometry
# to geom_text_repel, and use meter-scale nudge_y (not degrees).
major_cities_coords <- major_cities |>
  st_transform(25832) |>
  mutate(
    x = st_coordinates(geometry)[, 1],
    y = st_coordinates(geometry)[, 2]
  ) |>
  st_drop_geometry()

ggplot() +
  # Base: all states in light gray
  geom_sf(data = germany_states, fill = "#f0f0f0", color = "gray70",
          linewidth = 0.2) +
  # Highlight Bavaria
  geom_sf(data = bavaria, fill = "#FFF3CD", color = "#B8860B",
          linewidth = 0.4) +
  # All breweries
  geom_sf(data = breweries_sf, color = "#C2001A", size = 0.15,
          alpha = 0.35) +
  # Major cities as labeled points (geom_sf handles reprojection automatically)
  geom_sf(data = major_cities, color = "black", size = 2, shape = 21,
          fill = "white") +
  # Labels use the projected x/y data frame (nudge in metres, not degrees)
  ggrepel::geom_text_repel(
    data = major_cities_coords,
    aes(x = x, y = y, label = name),
    size = 3.2, fontface = "bold",
    nudge_y = 50000, segment.size = 0.3
  ) +
  coord_sf(crs = 25832) +
  theme_void(base_size = 14) +
  theme(plot.title = element_text(face = "bold", size = 16)) +
  labs(
    title = "Germany's breweries cluster in Bavaria and Franconia",
    caption = "Source: biermap24.de"
  )

Hex bin maps: handling overplotting

When you have thousands of points, individual dots overlap. Hex bins aggregate into a grid while preserving the spatial pattern.

Note: geom_hex() and coord_sf() do not mix well — the hex grid is computed in the coordinate space of the data before coord_sf() applies its transformation, leading to misaligned or empty maps. The fix is to project first, extract plain x/y coordinates, and then use coord_sf(crs = 25832) for the axis labels only:

# geom_hex() requires the hexbin package
library(hexbin)

# Project breweries to UTM 25832 for metric hex bins
# (geom_hex + coord_sf conflict: compute hex grid in projected coords,
#  then let coord_sf handle axis labels only)
breweries_proj <- st_transform(breweries_sf, 25832)

# Extract plain numeric coordinates as a data frame
brew_xy <- as.data.frame(st_coordinates(breweries_proj))

# Project state polygons too
germany_states_proj <- st_transform(germany_states, 25832)
germany_proj        <- st_transform(germany, 25832)

ggplot() +
  geom_sf(data = germany_proj, fill = NA, color = "gray40") +
  geom_hex(
    data = brew_xy,
    aes(x = X, y = Y),
    bins = 30, alpha = 0.85
  ) +
  scale_fill_viridis_c(
    option = "inferno",
    direction = -1,
    name = "Breweries\nper hex"
  ) +
  coord_sf(crs = 25832) +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = c(0.15, 0.4)
  ) +
  labs(
    title = "Brewery density across Germany",
    subtitle = "Hexagonal bins aggregate point locations into a density surface",
    caption = "Source: biermap24.de"
  )

Topographic backdrop: elevation data

For stories where terrain matters (e.g., showing how breweries follow river valleys, or mapping Alpine vs. lowland regions), adding an elevation backdrop adds a compelling layer. The elevatr package downloads elevation rasters on demand, and tidyterra integrates terra rasters with ggplot2.

# Requires internet access and additional packages:
# install.packages(c("elevatr", "tidyterra"))
# Note: downloading elevation data may take 30–60 seconds.

library(elevatr)
library(tidyterra)
library(terra)

# Get Bavaria boundary in WGS84 for elevation download
bavaria_wgs84 <- germany_states |>
  filter(name == "Bayern") |>
  st_transform(4326)

# Download elevation raster for Bavaria (z = 7 is ~1 km resolution)
bavaria_elev <- get_elev_raster(bavaria_wgs84, z = 7, clip = "bbox")

# Convert to terra SpatRaster and mask to Bavaria boundary
bavaria_elev_terra <- rast(bavaria_elev)
bavaria_elev_masked <- mask(bavaria_elev_terra, vect(bavaria_wgs84))

# Breweries clipped to Bavaria
breweries_bavaria <- st_intersection(
  st_transform(breweries_sf, 4326),
  bavaria_wgs84
)

ggplot() +
  tidyterra::geom_spatraster(data = bavaria_elev_masked) +
  scale_fill_hypso_c(
    palette = "dem_screen",
    na.value = NA,
    name = "Elevation (m)"
  ) +
  geom_sf(data = bavaria_wgs84, fill = NA, color = "gray30", linewidth = 0.4) +
  geom_sf(data = breweries_bavaria, color = "#C2001A", size = 0.6, alpha = 0.6) +
  theme_void(base_size = 14) +
  theme(plot.title = element_text(face = "bold", size = 16)) +
  labs(
    title = "Breweries in Bavaria against an elevation backdrop",
    subtitle = "Most Franconian breweries sit in the relatively flat northern part of Bavaria",
    caption = "Elevation: elevatr / Amazon Web Services Terrain Tiles | Breweries: biermap24.de"
  )

Hybrid maps with OpenStreetMap tiles

For maps that need geographic context (street names, rivers, settlements), raster tile basemaps from OpenStreetMap provide a familiar reference. The ggspatial package adds tiles and map accessories to ggplot2:

# Requires internet access and:
# install.packages(c("ggspatial", "prettymapr"))
# prettymapr is a backend dependency of annotation_map_tile()

library(prettymapr)
library(ggspatial)

# Clip breweries to Franconia / northern Bavaria bounding box
breweries_franconia <- breweries |>
  filter(lon >= 9.5, lon <= 12.5, lat >= 49.0, lat <= 51.0)

ggplot() +
  annotation_map_tile(type = "osm", zoom = 8) +
  geom_point(
    data = breweries_franconia,
    aes(x = lon, y = lat),
    color = "#C2001A", size = 1.2, alpha = 0.7
  ) +
  annotation_scale(location = "bl", width_hint = 0.2) +
  annotation_north_arrow(
    location = "tr",
    style = north_arrow_fancy_orienteering()
  ) +
  coord_sf(crs = 4326, xlim = c(9.5, 12.5), ylim = c(49.0, 51.0)) +
  theme_void(base_size = 13) +
  labs(
    title = "Breweries in Franconia",
    subtitle = "OpenStreetMap tile basemap with scale bar and north arrow",
    caption = "Tiles: © OpenStreetMap contributors | Breweries: biermap24.de"
  )

Replacing dots with symbols: beer icons

For a more editorial look, you can replace point markers with emoji icons. ggimage::geom_image() renders an image file (local path or URL) at each data coordinate. Using a CDN-hosted Twemoji PNG is the most reliable cross-platform approach — no emoji font dependency, consistent rendering everywhere. Use sparingly: best for small samples or magazine-style illustrations.

# Requires internet at render time to fetch the emoji PNG from the CDN.
# install.packages("ggimage")

library(ggimage)

# Use a random sample to minimize overplotting.
# Add 'icon' column pointing to locally stored beer mug png
set.seed(42)
breweries_sample <- breweries |>
  slice_sample(n = 150) |>
  mutate(icon = "data/beermug.png")

ggplot() +
  geom_sf(data = germany_states, fill = "#f5f5f0", color = "gray70",
          linewidth = 0.2) +
  geom_image(
    data = breweries_sample,
    aes(x = lon, y = lat, image = icon),
    size = 0.025
  ) +
  coord_sf(crs = 4326, xlim = c(5.5, 15.5), ylim = c(47, 55.5)) +
  theme_void(base_size = 14) +
  labs(
    title = "150 randomly sampled German breweries",
    caption = "Source: biermap24.de"
  )

Interactive maps with leaflet

For web publication, interactive maps let readers explore the data themselves. The leaflet package creates HTML widgets that work in RMarkdown, Quarto, and Shiny. This fully-featured example includes tile switching, marker clustering, popups, a minimap, and a scale bar:

# Produces an HTML widget — works in Quarto HTML output with embed-resources: true

library(leaflet)

leaflet(breweries) |>
  setView(lng = 10.5, lat = 51.2, zoom = 6) |>
  # Tile layer providers
  addProviderTiles("CartoDB.Positron",  group = "CartoDB Positron") |>
  addProviderTiles("Esri.WorldImagery", group = "Satellite") |>
  addProviderTiles("OpenStreetMap",     group = "OpenStreetMap") |>
  # Brewery markers with clustering
  addCircleMarkers(
    lng = ~lon, lat = ~lat,
    radius = 4,
    color = "#C2001A",
    fillOpacity = 0.7,
    stroke = FALSE,
    popup = ~city,
    group = "Breweries",
    clusterOptions = markerClusterOptions()
  ) |>
  # Layer control
  addLayersControl(
    baseGroups    = c("CartoDB Positron", "Satellite", "OpenStreetMap"),
    overlayGroups = "Breweries",
    options       = layersControlOptions(collapsed = FALSE)
  ) |>
  addMiniMap(toggleDisplay = TRUE) |>
  addScaleBar(position = "bottomright")

Best practices for journalistic maps

Here are concrete guidelines drawn from data journalism practice:

1. Choose the right map type for the story. Point maps for discrete events (crimes, businesses, incidents). Choropleth maps for rates or aggregated data (unemployment rate, vote share). Avoid choropleths for raw counts — normalize first.

2. Use sequential color scales for quantitative data. Single-hue scales (e.g., light to dark blue) work for “more vs. less.” Diverging scales (e.g., blue–white–red) work when there’s a meaningful midpoint (e.g., above/below average). Never use rainbow scales — they have no perceptual ordering.

3. Keep it minimal. Strip away grid lines, axis labels, and unnecessary borders. Use theme_void() as your starting point. Add labels only where they serve the reader.

4. Mind the visual area bias. Large regions attract attention regardless of their data values. Consider alternatives like cartograms, hex maps, or dot density maps when geographic area could mislead.

5. Show your sources. Always include a data source caption. If the boundaries come from a specific year or administrative revision, say so.

6. Consider interactivity. For web publication, leaflet or mapview let readers zoom and click. But a well-designed static map is often better for print and social sharing. See the leaflet section above.

Bonus: Modern Geospatial Packages Worth Knowing

The sf + ggplot2 workflow covers most data journalism needs, but a growing ecosystem of specialised packages extends what is possible. Here is a curated overview of well-maintained packages worth knowing about.

mapview — one-line interactive map exploration

mapview wraps Leaflet into a single function call. Pass any sf object and get a fully interactive map instantly — great for exploratory data analysis when you just want to see the data.

library(mapview)

# One line to get an interactive map with popups, zoom, and tile switcher
mapview(breweries_sf, zcol = NULL, col.regions = "#C2001A", cex = 2)

Good for: interactive output, quick EDA. Not suited for publication-quality design.

mapedit — draw and edit spatial features interactively

mapedit opens a Shiny-based editor in the browser where you can draw points, lines, and polygons by hand, or edit existing features. The result is returned as an sf object directly in your R session.

library(mapedit)
library(mapview)

# Launch interactive editor; draw features; result comes back as sf
my_features <- editMap(mapview(germany_states))
# my_features$drawn is an sf object of what you drew

Good for: interactive input/editing. Useful for defining regions of interest, manually cleaning spatial data, or teaching.

mapdeck — WebGL-powered big-data maps

mapdeck uses the deck.gl JavaScript library (from Uber) via a Mapbox token to render millions of data points in the browser using WebGL. It supports H3 hexagon layers, arc layers, animated trips, and more.

library(mapdeck)

# Requires a free Mapbox access token: https://account.mapbox.com/
set_token("YOUR_MAPBOX_TOKEN")

mapdeck(style = mapdeck_style("dark")) |>
  add_scatterplot(
    data    = breweries,
    lon     = "lon",
    lat     = "lat",
    radius  = 3000,
    fill_colour = "#C2001A",
    tooltip = "city"
  )

Good for: interactive output. Ideal when you have hundreds of thousands of points and need GPU-accelerated rendering.

ggspatial — scale bars, north arrows, and tile basemaps

ggspatial extends ggplot2 with three map accessories that are otherwise painful to add manually: scale bars, north arrows, and raster tile basemaps (from OSM, Stamen, ESRI, and others).

library(ggspatial)

ggplot() +
  annotation_map_tile(type = "cartolight", zoom = 7) +
  geom_sf(data = st_transform(breweries_sf, 4326),
          color = "#C2001A", size = 0.8, alpha = 0.5) +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tr",
                          style = north_arrow_minimal()) +
  coord_sf(crs = 4326, xlim = c(5.5, 15.5), ylim = c(47, 55.5)) +
  theme_void()

Good for: static output. The tile basemap requires internet access at render time.

osmdata — query OpenStreetMap’s Overpass API

osmdata provides a direct interface to OSM’s Overpass API. You can query for any OSM feature — streets, buildings, parks, cafes, hospitals, breweries — anywhere in the world, and get the result as an sf object.

library(osmdata)

# Query OSM for breweries in Bavaria
bavaria_breweries_osm <- opq("Bavaria, Germany") |>
  add_osm_feature(key = "craft", value = "brewery") |>
  osmdata_sf()

# Points are in $osm_points, polygons (buildings) in $osm_polygons
head(bavaria_breweries_osm$osm_points)

ggplot() +
  geom_sf(data = st_transform(bavaria, 25832), fill = "gray95") +
  geom_sf(data = bavaria_breweries_osm$osm_points,
          color = "#C2001A", size = 1) +
  coord_sf(crs = 25832) +
  theme_void() +
  labs(title = "Bavarian breweries from OpenStreetMap")

Good for: static or interactive output. The definitive source for granular, locally-contributed geographic features.

Beyond R: other tools worth knowing

R is not the only tool in the data journalist’s spatial toolkit. Knowing when to reach for something else saves time.

QGIS is the leading open-source desktop GIS application. It provides a graphical interface for loading, reprojecting, editing, and visualising geodata — all without writing code. QGIS is indispensable for exploring unfamiliar datasets, manually digitising features, and producing print-quality map layouts. It reads every format sf supports (shapefiles, GeoPackage, PostGIS, WFS, etc.) and has a plugin ecosystem covering everything from field data collection to network routing. Many newsrooms use QGIS alongside R: R for analysis and automation, QGIS for visual QA and final layout polish.

ArcGIS / ArcGIS Online (Esri) is the commercial industry standard in government and large organisations. ArcGIS Online lets you publish and share interactive web maps without any programming. For data journalists it is most relevant as a source of authoritative boundary data and a target format — many public agencies publish data only as Esri File Geodatabases or via ArcGIS REST services, both of which sf can read.

Google Earth Engine is a cloud platform for planetary-scale analysis of satellite imagery and raster time series. If your story involves global deforestation, flood mapping, or any analysis requiring terabytes of satellite data, Earth Engine is the right tool. It cannot be replicated on a laptop — the compute runs in Google’s data centres, and you query it via a JavaScript or Python API.

Mapshaper is a fast, browser-based tool for simplifying, inspecting, and converting vector data. It is the quickest way to reduce a bloated shapefile to a web-friendly size or convert between formats without writing any code. Indispensable for web mapping preparation.

Exercises

NoteExercise 1: European capitals — who’s closest to whom?

Use tidygeocoder::geocode() to geocode 10 European capital cities of your choice. Convert the result to an sf object with st_as_sf().

(a) Plot the capitals as labeled points on a base map of Europe (ne_countries(continent = "Europe", returnclass = "sf")). Use ggrepel for the labels.

(b) Compute a full pairwise distance matrix with st_distance(). Which two capitals in your list are furthest apart? Which are nearest?

(c) Highlight the furthest pair by drawing a connecting line on the map. Hint: st_union() the two points, then st_cast("LINESTRING") to get a line geometry.

NoteExercise 2: Coffee culture — an OSM snapshot

Pick any major European city you like. Use osmdata to download all locations tagged amenity = "cafe".

(a) Create a kernel density heatmap (stat_density_2d) of café locations on a plain city base map.

(b) Geocode the city’s main square or central train station. Buffer it by 1 km (remember to reproject to a metric CRS first). How many cafés fall within walking distance?

(c) Plot the result on an OpenStreetMap tile basemap using ggspatial::annotation_map_tile(). Add a scale bar.

NoteExercise 3: Projections and global inequality

Install gapminder (install.packages("gapminder")). Join the gapminder dataset filtered to the year 2007 to ne_countries(scale = "medium", returnclass = "sf") world polygons.

(a) Map lifeExp as a choropleth using an equal-area projection (try crs = "+proj=eqearth").

(b) Remake the same map using the Mercator projection (crs = "+proj=merc").

(c) Write 2–3 sentences: which regions draw the eye differently under each projection? Does the choice of projection change the story the map appears to tell?

NoteExercise 4: Real-time geodata — earthquakes

The USGS publishes live earthquake data as a GeoJSON feed you can read directly into R:

earthquakes <- st_read(
  "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_month.geojson"
)

(a) Plot all earthquakes from the past month on a world map, sizing each point by mag (magnitude). Use an equal-area or orthographic projection.

(b) Filter to M ≥ 6 events and add a 500 km buffer around each (don’t forget to reproject!). Which world regions show the highest overlap of buffers?

(c) Perform a spatial join with ne_countries() to count how many M ≥ 5 earthquakes occurred in each country. Which country tops the list?

Further reading

sessionInfo()
R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.3.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] hexbin_1.28.5           prettymapr_0.2.5        tidyterra_1.1.0        
 [4] elevatr_0.99.1          ggimage_0.3.5           leaflet_2.2.3          
 [7] ggspatial_1.1.10        ggrepel_0.9.7           patchwork_1.3.2        
[10] tidygeocoder_1.0.6      nominatimlite_0.4.3     rnaturalearthdata_1.0.0
[13] rnaturalearth_1.2.0     terra_1.9-1             sf_1.1-0               
[16] lubridate_1.9.5         forcats_1.0.1           stringr_1.6.0          
[19] dplyr_1.2.0             purrr_1.2.1             readr_2.2.0            
[22] tidyr_1.3.2             tibble_3.3.1            ggplot2_4.0.2          
[25] tidyverse_2.0.0        

loaded via a namespace (and not attached):
 [1] DBI_1.3.0                     s2_1.1.9                     
 [3] rlang_1.1.7                   magrittr_2.0.4               
 [5] otel_0.2.0                    e1071_1.7-17                 
 [7] compiler_4.5.3                png_0.1-8                    
 [9] systemfonts_1.3.2             vctrs_0.7.1                  
[11] pkgconfig_2.0.3               wk_0.9.5                     
[13] crayon_1.5.3                  fastmap_1.2.0                
[15] magick_2.9.1                  labeling_0.4.3               
[17] rmarkdown_2.30                tzdb_0.5.0                   
[19] xfun_0.56                     rnaturalearthhires_1.0.0.9000
[21] jsonlite_2.0.0                progress_1.2.3               
[23] rosm_0.3.1                    prettyunits_1.2.0            
[25] R6_2.6.1                      stringi_1.8.7                
[27] RColorBrewer_1.1-3            jquerylib_0.1.4              
[29] Rcpp_1.1.1                    knitr_1.51                   
[31] pacman_0.5.1                  leaflet.providers_2.0.0      
[33] timechange_0.4.0              tidyselect_1.2.1             
[35] abind_1.4-8                   rstudioapi_0.18.0            
[37] yaml_2.3.12                   codetools_0.2-20             
[39] curl_7.0.0                    lattice_0.22-9               
[41] plyr_1.8.9                    withr_3.0.2                  
[43] S7_0.2.1                      evaluate_1.0.5               
[45] gridGraphics_0.5-1            units_1.0-1                  
[47] proxy_0.4-29                  isoband_0.3.0                
[49] pillar_1.11.1                 KernSmooth_2.23-26           
[51] ggfun_0.2.0                   generics_0.1.4               
[53] sp_2.2-1                      hms_1.1.4                    
[55] scales_1.4.0                  class_7.3-23                 
[57] glue_1.8.0                    gdtools_0.5.0                
[59] tools_4.5.3                   data.table_1.18.2.1          
[61] ggiraph_0.9.6                 fs_1.6.7                     
[63] grid_4.5.3                    slippymath_0.3.1             
[65] crosstalk_1.2.2               raster_3.6-32                
[67] cli_3.6.5                     rappdirs_0.3.4               
[69] fontBitstreamVera_0.1.1       viridisLite_0.4.3            
[71] gtable_0.3.6                  yulab.utils_0.2.4            
[73] digest_0.6.39                 fontquiver_0.2.1             
[75] progressr_0.18.0              classInt_0.4-11              
[77] ggplotify_0.1.3               htmlwidgets_1.6.4            
[79] farver_2.1.2                  htmltools_0.5.9              
[81] lifecycle_1.0.5               httr_1.4.8                   
[83] fontLiberation_0.1.0          MASS_7.3-65