Imperviousness levels in and around Berlin, Germany

imp_berlin.R

library(tidyverse)
library(sf)
library(terra)
library(stars)
library(ggspatial)
library(systemfonts)
library(patchwork)

register_variant(
  name = "Input Mono Light",
  family = "Input Mono", 
  weight = "light"
)

register_variant(
  name = "Cabinet Grotesk ExtraBold",
  family = "Cabinet Grotesk", 
  weight = "ultrabold"
)

theme_set(theme_void(base_family = "Cabinet Grotesk", base_size = 15))
theme_update(
  legend.position = "top",  
  plot.title = element_text(
    size = rel(1.67), family = "Cabinet Grotesk ExtraBold", 
    face = "bold", margin = margin(5, 0, 10, 0)
  ),
  plot.title.position = "plot",
  plot.caption = element_text(
    size = rel(.8), color = "grey40", margin = margin(10, 0, -3, 0)
  ),
  plot.tag = element_text(
    face = "bold", margin = margin(0, 0, -10, 0), size = rel(1.25)
  ),
  axis.text = element_text(
    size = rel(.8), family = "Input Mono Light", 
    color = "grey40", margin = margin(rep(5, 4))
  ),
  axis.ticks.length = unit(.4, "lines"), axis.ticks = element_line(color = "grey75"),
  plot.margin = margin(rep(10, 4))
)

imperv <- rast(here::here("data", "imp_bb_mv_b_20m_3035.tif"))
b <- as(extent(4520000, 4585000, 3246000, 3301000), 'SpatialPolygons')
crs(b) <- crs(imperv)
imperv_crop <- st_as_stars(crop(imperv, b))
imperv_crop_200m <- st_as_stars(aggregate(crop(imperv, b), fact = 10, fun = "mean"))
imperv_crop_500m <- st_as_stars(aggregate(crop(imperv, b), fact = 25, fun = "mean"))
imperv_crop_1000m <- st_as_stars(aggregate(crop(imperv, b), fact = 50, fun = "mean"))

## imperviousness map
map_imp <- 
  ggplot() +
  geom_stars(data = imperv_crop) +
  geom_sf(data = d6berlin::sf_berlin, fill = "transparent", color = "grey80") +
  scale_fill_scico(
    palette = "batlow", begin = .1, limits = c(0, 100),
    name = "Imperviousness", breaks = 0:10*10, labels = function(x) paste0(x, "%"),
    guide = guide_colorsteps(barwidth = unit(28, "lines"), barheight = unit(.6, "lines"),
                             title.position = "top", title.hjust = .5, show.limits = TRUE)
  ) +
  annotation_scale(
    location = 'tr', text_family = "Input Mono", text_cex = 1, text_col = "white"
  ) +
  coord_sf(
    expand = FALSE, crs = crs(imperv)
  ) +
  labs(x = NULL, y = NULL)


## overview map
sf_world <- 
  st_as_sf(rworldmap::getMap(resolution = "low")) %>% 
  st_transform(crs = st_crs(imperv)) %>% 
  st_buffer(dist = 0) %>% 
  dplyr::select(ISO_A2, SOVEREIGNT, LON, continent) %>% 
  mutate(area = st_area(.))

map_europe <- 
  ggplot(sf_world) +
  geom_sf(fill = "grey80", color = "grey95", lwd = .1) +
  geom_rect(
    xmin = 4500000, xmax = 4600000, ymin = 3230000, ymax = 3310000,
    color = "#212121", fill = "#b28c32", size = .7
  ) +
  geom_sf_text(
    data = filter(sf_world, ISO_A2 == "DE")),
    aes(label = ISO_A2),
    family = "Open Sans", color = "grey40", fontface = "bold", size = 4.5,
    nudge_x = 20000, nudge_y = -10000
  ) +
  annotation_scale(
    location = 'tr', text_family = "Input Mono", text_cex = 1
  ) +
  coord_sf(xlim = c(2650000, 5150000), ylim = c(1650000, 5100000)) +
  scale_x_continuous(expand = c(0, 0), breaks = seq(-10, 30, by = 10)) +
  theme(panel.ontop = FALSE,
        panel.grid.major = element_line(color = "grey75", linetype = "15", size = .3)) +
  labs(x = NULL, y = NULL)


## globe
map_globe <- d6berlin::globe(col_earth = "grey80", col_water = "grey95", bg = TRUE)


### combined map
map_overview <- map_europe + labs(tag = "A.") + inset_element(map_globe, .14, .75, .56, 1.2, align_to = "plot")

m <- map_overview + (map_imp + labs(tags = "B.")) +
  plot_annotation(
    title = "Percentages of sealed soil in and around Berlin, Germany, on a 20 m scale",
    caption = "Source: Copernicus Land Monitoring Service"
  )

ggsave("berlin_imp.pdf", width = 12, height = 7.7, bg = "white", device = cairo_pdf)