Generate a basemap for Germany and combine it with travel distance data

distanceplot.R

# Generate a basemap for Germany and combine it with travel distance data

# Init --------------------------------------------------------------------

library(tidyverse)
library(sf)
library(eurostat)

# Constants ---------------------------------------------------------------

# input and output paths
paths <- list()
paths$input <- list(
  distancedata.rds = "distancedata.rds"
)
paths$output <- list(
  distanceplot.rds = "distanceplot.rds",
  distanceplot.png = "distanceplot.png"
)

# list containers for analysis artifacts
distanceplot <- list()

# bounding box which crops to Germany and surrounding countries
# coordinates of bounding box expressed in Pseudo-Mercator coordinates
# https://epsg.io/3857
distanceplot$bounding_box <- c(
  xmin =  633500, ymin = 5944900,
  xmax = 1700000, ymax = 7399100
)

# the crs we transform all spatial data to
# project to Pseudo Mercator as used by OSM or Google Maps
# https://epsg.io/3857
distanceplot$target_crs <- 3857

# Create region layer ---------------------------------------------

# district outlines in SF format
distanceplot$regions <-
  # just to be explicit we state that the coordinates are expressed
  # in latitude-longitude format, i.e. https://epsg.io/4326
  get_eurostat_geospatial(
    nuts_level = "1", # German Bundesländer
    crs = 4326,
    resolution = 03
  ) |>
  # select only German regions
  filter(
    CNTR_CODE %in% "DE"
  ) |>
  # project to target CRS
  st_transform(crs = st_crs(distanceplot$target_crs)) |>
  # crop to Germany and surrounding countries
  # coordinates of bounding box expressed in Mercator coordinates
  st_crop(y = distanceplot$bounding_box)

ggplot(distanceplot$regions) +
  geom_sf()

# extract internal borders
distanceplot$internal_borders <-
  distanceplot$regions |>
  st_boundary() |>
  st_union() |>
  st_difference(
    distanceplot$regions |>
      st_union() |>
      st_boundary()
  ) |>
  st_collection_extract("LINESTRING")

ggplot(distanceplot$internal_borders) +
  geom_sf()

# Create background layer for map of Germany ------------------------------

distanceplot$background <-
  get_eurostat_geospatial(
    nuts_level = "0", # nation states
    crs = 4326,
    resolution = 03
  ) |>
  # project to target CRS
  st_transform(crs = st_crs(distanceplot$target_crs)) |>
  # crop to Germany and surrounding countries
  # coordinates of bounding box expressed in Mercator coordinates
  st_crop(y = distanceplot$bounding_box)

ggplot(distanceplot$background) +
  geom_sf()

# Create outline of Germany -----------------------------------------------

# make a union of all regions belonging to France and return
# the outline of the resulting polygon
distanceplot$outline <-
  distanceplot$regions |>
  summarise(id = "de") |>
  st_union()

ggplot(distanceplot$outline) +
  geom_sf(fill = NA, lwd = 1, color = "black") +
  theme_void()

# Merge basemap with actual data ------------------------------------------

# cut travel distances into intervals
distanceplot$colorscale <- tribble(
  ~cut, ~color, ~label,
  0, "#b6f9d9", "under 10",
  10, "#a6e17c", "10 to <20",
  20, "#cda149", "20 to <30",
  30, "#da5b68", "30 to <40",
  40, "#993a8f", "over 40"
)

# travel time data
distanceplot$distances <-
  readRDS(paths$input$distancedata.rds) |>
  # project to pseudo Mercator
  st_transform(crs = st_crs(distanceplot$target_crs)) |>
  select(distance = dist_2021) |>
  mutate(
    distance = cut(distance, c(distanceplot$colorscale$cut, Inf), right = FALSE)
  )


# plot travel distance
distanceplot$fig <-
  ggplot() +
  # background
  geom_sf(
    data = distanceplot$background,
    color = "grey50"
  ) +
  # distance raster
  geom_sf(
    aes(fill = distance, color = distance),
    data = distanceplot$distances
  ) +
  # internal state borders only
  geom_sf(
    color = "black", lwd = 0.2,
    data = distanceplot$internal_borders
  ) +
  # germany outline
  geom_sf(
    fill = NA, lwd = 0.2, color = "black",
    data = distanceplot$outline
  ) +
  scale_color_manual(
    values = distanceplot$colorscale$color,
    labels = distanceplot$colorscale$label,
    aesthetics = c("colour", "fill")
  ) +
  theme_void()
distanceplot$fig

# Export ----------------------------------------------------------

saveRDS(distanceplot, paths$output$distanceplot.rds, compress = "xz")
ggsave(
  filename = paths$output$distanceplot.png,
  plot = distanceplot$fig,
  device = "png",
  width = 175,
  dpi = 300,
  units = "mm",
  bg = "white"
)