Skip to content

Instantly share code, notes, and snippets.

@z3tt
Last active November 22, 2022 03:53
Show Gist options
  • Save z3tt/a255ed73ad31601bb418be7f94073596 to your computer and use it in GitHub Desktop.
Save z3tt/a255ed73ad31601bb418be7f94073596 to your computer and use it in GitHub Desktop.
Imperviousness levels in and around Berlin, Germany
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment