Durchschnittliches Alter und Einkommen auf Gemeindeebene (https://twitter.com/matschnetzer/status/1606233928745988096)

ageinc.R

library(tidyverse)
library(sf)
library(msthemes)
library(raster)
library(rmapshaper)
library(smoothr)
library(patchwork)
library(showtext)
select <- dplyr::select

# Schriftarten hinzufügen
font_add_google(name = "Playfair Display", family = "Playfair Display")
font_add_google(name = "Roboto Condensed", family = "Roboto")
showtext_auto()
showtext_opts(dpi = 600)

# Daten-Download von Statistik Austria und data.gv.at
# Gemeindegrenzen: https://data.statistik.gv.at/web/meta.jsp?dataset=OGDEXT_GEM_1
# Dauersiedlungsraum: https://data.statistik.gv.at/web/meta.jsp?dataset=OGDEXT_DSR_1
# Gewässer: https://www.data.gv.at/katalog/dataset/ce50ffa6-5032-4771-90a2-1c48d6a0ac85 
# Alter und Einkommen: https://www.statistik.at/atlas/  


# Karte mit Siedlungsraum erstellen
gem <- read_sf("STATISTIK_AUSTRIA_GEM_20220101.shp") |>
  ms_simplify(keep = 0.1)

sied <- read_sf("STATISTIK_AUSTRIA_DSR_20111031.shp") |>
  filter(ID %in% 2:3) |>
  ms_simplify(keep = 0.01) |>
  st_buffer(dist = 1000, endCapStyle = "SQUARE",
            joinStyle = "MITRE", nQuadSegs = 2) |>
  summarise(geometry = st_union(geometry)) |>
  fill_holes(threshold = units::set_units(200, km^2))

wat <- read_sf("stehendeGewaesser.shp") |>
  st_simplify(dTolerance = 100)

siedmap <- st_intersection(gem, sied)

# Bundesgrenze erstellen
autborder <- gem |> st_union() |> st_as_sf()

# Originaldatei (zusätzlich extern bearbeitet): https://www.isticktoit.net/?p=483 
relief <- raster("relief.tif") |>
  raster::mask(autborder) |> 
  as("SpatialPixelsDataFrame") |>
  as.data.frame()

# Einkommens- und Alteresdaten laden
inc <- read.csv("lohnsteuerstatistik-jahresbruttobezug.csv", skip = 7, sep = ";") |> 
  select(id = ID, name = Name, inc = "X.") |> 
  mutate(id = ifelse(name == "Matrei am Brenner", 70370, id))
age <- read.csv("bevölkerung_nach_alter.csv", skip = 7, sep = ";", dec = ",") |> 
  select(id = ID, name = Name, age = Wert)
df <- left_join(age, inc)

# Bundesländergrenzen
bldborder <- gem |> 
  mutate(bl = as.factor(substr(id, 1, 1))) |> 
  group_by(bl) |> summarise()

# Farbpalette erstellen
colorscale <- tribble(
  ~group, ~fill,
  "3 - 3", "#413079",
  "2 - 3", "#49708d",
  "1 - 3", "#48a065", 
  "3 - 2", "#685891",
  "2 - 2", "#6f8ba0", 
  "1 - 2", "#85b798",
  "3 - 1", "#8e82ab", 
  "2 - 1", "#93a6b4",
  "1 - 1", "#a3c5af"
)

# Variablen erstellen und Farbpalette hinzufügen
plotmap <- siedmap |> 
  left_join(df |> mutate(id = as.character(id)), by = "id") |> 
  drop_na() |> 
  mutate(incg = base::cut(inc, breaks = c(0,40000,50000,100000)),
         ageg = base::cut(age, breaks = c(30,40,45,60)),
         group = paste(as.numeric(incg), "-", as.numeric(ageg))) |> 
  left_join(colorscale)

# Karte erstellen
map <- ggplot() +
  geom_raster(data=relief, interpolate=T, aes(x = x, y = y, fill = "white", alpha = layer)) + 
  scale_alpha_continuous(name = "", range = c(0.7,0.9), guide = F) +
  geom_sf(data = plotmap, aes(fill = fill), color = NA) +
  geom_sf(data = bldborder, aes(fill = NA), color = "black", size = 0.3) +
  geom_sf(data = wat, fill = "#D6F1FF", color = "transparent") +
  scale_fill_identity() +
  coord_sf(datum = NA) +
  theme_ms(alttf = T, dark = T) +
  labs(title="Einkommen und Alter",
       subtitle="Durchschnittswerte auf Gemeindeebene in Österreich, 2021",
       caption="Daten: Statistik Austria, Lohnsteuerstatistik. Grafik: @matschnetzer") +
  theme(axis.title = element_blank(),
        legend.position = "none",
        plot.background = element_rect(fill="gray10", color=NA),
        panel.background = element_rect(fill="gray10", color=NA),
        panel.border = element_blank(),
        plot.title = element_text(hjust=0.5),
        plot.subtitle = element_text(hjust=0.5, family = "Playfair Display"),
        plot.caption = element_text(size = 6, family = "Roboto")) +
  expand_limits(x = 800000, y = 260000) +
  
  # Text und Pfeile 
  annotate("text",label="Dunkelgrün bedeutet niedrige Einkommen\nund hohes Durchschnittsalter, wie in der\nSüdost-Steiermark oder im Waldviertel", size=2, family="Roboto", hjust=0.5, x=720000, y=320000, color = "white") +
  annotate("text",label="Mittlere Einkommen im Inntal\n und niedrige Einkommen in den\n Seitentälern Tirols.", size=2, family="Roboto", hjust=0.5, x=200000, y=280000, color = "white") +
  annotate("text",label="Dunkles Violett bedeutet hohe\n Einkommen und hohes Alter,\n z.B. im Wiener Speckgürtel.", size=2, family="Roboto", hjust=0.5, x=750000, y=510000, color = "white") +
  geom_curve(aes(x=240000,xend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Innsbruck"]))[1], y=300000,yend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Innsbruck"]))[2]), curvature = 0.1, ncp=8, linewidth=0.1, color = "white",
             arrow=arrow(length=unit(0.01, "npc"), type="closed")) +
  geom_curve(aes(x=640000,xend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Straden"]))[1], y=320000,yend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Straden"]))[2]), curvature = -0.1, ncp=8, linewidth=0.1, color = "white",
             arrow=arrow(length=unit(0.01, "npc"), type="closed")) +
  geom_curve(aes(x=690000,xend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Klosterneuburg"]))[1], y=510000,yend=st_coordinates(st_centroid(siedmap$geometry[siedmap$name == "Klosterneuburg"]))[2]), curvature = 0.1, ncp=8, linewidth=0.1, color = "white",
             arrow=arrow(length=unit(0.01, "npc"), type="closed")) 

# Legende erstellen
collegend <- colorscale |> 
  separate(group, into = c("inc", "age"), sep = " - ") |> 
  mutate(across(c(inc,age), as.numeric))
legend <- ggplot() +
  geom_tile(data = collegend,
            mapping = aes(x = inc, y = age, fill = fill)) +
  scale_fill_identity() +
  scale_x_continuous(breaks= 1:3,labels=c("<40k €","",">50k €")) +
  scale_y_continuous(breaks= 1:3,labels=c("<40","",">45")) +
  labs(x = "Einkommen",
       y = "Alter") +
  theme_ms(dark = T, grid = F) +
  theme(
    plot.background = element_rect(fill = "gray10", color = NA),
    panel.background = element_rect(fill = "gray10", color = NA),
    axis.title = element_text(size = 6.5, hjust = 0.5, family = "Roboto"),
    axis.text.x = element_text(size=5.5, hjust=0.5, margin = margin(r = 0), 
                               family = "Roboto"),
    axis.text.y = element_text(size=5.5, hjust=0.5, margin = margin(r = 0), 
                               family = "Roboto", angle=90)) +
  coord_fixed()

# Karte und Legende zusammenfügen
map + inset_element(legend, left = 0, bottom = 0.6, right = 0.35, top = 0.95) &
  plot_annotation(theme = theme(plot.background = element_rect(fill ="gray10", color = NA)))
  
ggsave("incagemap.png", width = 8, height = 4.21, dpi=600)