Bivariate map for vaccinations and covid cases in Austria (https://twitter.com/matschnetzer/status/1478750626339303424)

covidmap_bivariate.R

############# 
#### ACKNOWLEDGMENTS
#############

# Thanks for the original idea and the helpful code to @grssnbchr and @benjazehr: https://timogrossenbacher.ch/2019/04/bivariate-maps-with-ggplot2-and-sf/

############# 
#### LOAD PACKAGES
#############

library(tidyverse)
library(sf)
library(rmapshaper)
library(smoothr)
library(msthemes)
library(patchwork)

############# 
#### GET DATA
#############

## District map (https://data.statistik.gv.at/web/meta.jsp?dataset=OGDEXT_POLBEZ_1)
bez <- read_sf("STATISTIK_AUSTRIA_POLBEZ_20210101.shp") %>%
        ms_simplify(keep = 0.1)

## Water areas (https://www.data.gv.at/katalog/dataset/2217e33c-6bd3-4fbf-8344-51dfb76f5ff6)
wat <- read_sf("stehendeGewaesser.shp") %>%
  st_simplify(dTolerance = 100)

## Populated areas (https://data.statistik.gv.at/web/meta.jsp?dataset=OGDEXT_DSR_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)) 

## Get Covid data
cov <- read_csv2("https://covid19-dashboard.ages.at/data/CovidFaelle_GKZ.csv") %>%
  mutate(inzidenz = AnzahlFaelle7Tage/AnzEinwohner*100000)

## Get vaccination data
impf <- read_csv2("https://info.gesundheitsministerium.gv.at/data/COVID19_vaccination_municipalities.csv") %>%
  mutate(id = substr(municipality_id, 1, 3),
         id = ifelse(str_starts(id,"9"), 900, id)) %>%
  group_by(id) %>% 
  summarise(totpop = sum(municipality_population), 
            totimpf = sum(valid_certificates)) %>%
  mutate(impfquote = totimpf/totpop*100,
         id = as.numeric(id))


############# 
#### RELIEF (THAT'S THE HARDEST PART!!)
#############

# Found a Relief Map for Austria here: https://www.isticktoit.net/?p=483
# Convert the jpeg into GeoTiff with gdal in your Terminal:
### gdalwarp -s_srs EPSG:3035 -t_srs EPSG:31287 -r bilinear -srcnodata 0,0,0 -dstalpha shade25AT1.jpg relief_t.tif
### gdal_translate -of GTiff -co tiled=yes -outsize 3000 0 -scale 150 200 0 255 relief_t.tif relief.tif
# I'm not sure if this is the best way, but it worked
# Then go back to R, read the tif, intersect with Austrian border, and save as data frame

autborder <- bez %>% st_union() %>% as_Spatial()
relief <- raster("relief.tif") %>%
  raster::mask(autborder) %>%
  as("SpatialPixelsDataFrame") %>% 
  as.data.frame()


############# 
#### PREPARATIONS
#############

# Create states borders
bldborder <- bez %>% mutate(bl = as.factor(substr(id,1,1))) %>% group_by(bl) %>% split(.$bl) %>% 
  lapply(st_union) %>% do.call(c,.) %>% st_cast()

# Replace districts of Vienna with city borders
vienna <- bez %>% filter(str_starts(id, "9")) %>% summarise(geometry = st_union(geometry)) %>% 
  mutate(id = "900", name = "Wien")
bez <- bind_rows(bez %>% filter(!str_starts(id, "9")), vienna)

# Intersect districts map and populated areas
siedbez <- st_intersection(bez, sied)

# Combine map with data
plotmap <- siedbez %>% mutate(id = as.numeric(id)) %>%
  left_join(cov %>% dplyr::select(GKZ,inzidenz), by=c("id"="GKZ")) %>%
  left_join(impf) %>%
  mutate(impfq = base::cut(impfquote, breaks = c(0,65,70,100), labels = 1:3),
         covq = base::cut(inzidenz, breaks = c(0,200,300,Inf), labels = 3:1))

# Create colorscale
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")

# Assign colors to data
plotmap %<>% mutate(group = paste(impfq, "-", covq)) %>%
  left_join(colorscale, by="group")


############# 
#### CREATE MAP
#############

# Map
map <- ggplot() +
  geom_raster(data=relief, interpolate=T, aes(x=x, y=y, alpha=relief)) + 
  scale_alpha_continuous(name="", range = c(0, 0.3), guide = F) +
  geom_sf(data=plotmap, aes(fill=fill), color="grey98", size=0.07) +
  geom_sf(data=bldborder, aes(fill=NA), color="black", size=0.3) +
  geom_sf(data=wat, fill="#D6F1FF", color="transparent") +
  annotate("text",label="Anzahl der\n Bezirke", size=1.8, family="IBMPlexSans", 
           hjust=0.5, vjust=1, x=320000, y=530000) +
  geom_curve(aes(x=315000,xend=280000, y=535000,yend=535000), curvature = 0.2, 
             ncp=8, size=0.1, arrow=arrow(length=unit(0.01, "npc"), type="closed")) +
  scale_fill_identity() +
  coord_sf(datum = NA) +
  theme_ms(alttf = T) +
  labs(title="Covid-19 in Österreich",
       subtitle="Impfquote und 7-Tages-Inzidenz nach Bezirken, 2022",
       caption="Daten: AGES (3.1.2022). Grafik: @matschnetzer") +
  theme(axis.title = element_blank(),
        legend.position = "none",
        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 = "IBMPlexSans"))


# Number of districts for legend
anzahl <- as.data.frame(plotmap) %>% group_by(group) %>% count()
colorscale <- colorscale %>% left_join(anzahl) %>%
  separate(group, into = c("impf", "cov"), sep = " - ") %>%
  mutate(across(c(impf, cov), as.numeric))
  
# Legend
legend <- ggplot(data = colorscale, aes(x = impf, y = cov)) +
  geom_tile(aes(fill = fill)) +
  geom_text(aes(label=n), color="white", size = 2.5, family = "IBMPlexSans")+
  scale_fill_identity() +
  scale_y_continuous(breaks= 1:3,labels=c(">300","","<200")) +
  scale_x_continuous(breaks= 1:3,labels=c("<65%","",">70%")) +
  labs(x = "Höhere Impfquote →",
       y = "Niedrigere Inzidenz →") +
  theme_ms(grid=F) +
  theme(
    plot.background = element_rect(fill="transparent"),
    axis.title = element_text(size = 5, hjust = 0, family = "IBMPlexSans"),
    axis.text.y = element_text(angle=90, size=4, hjust=0.5, margin = margin(r = 0)),
    axis.text.x = element_text(size=4, hjust=0.5, margin = margin(r = 0)),
    axis.line = element_blank(),
    axis.ticks = element_blank()
  ) +
  coord_fixed()

# Assemble plot with patchwork
map + inset_element(legend, left = 0.15, right = 0.45, bottom = 0.6, top = 0.90)

# Save plot
ggsave("covidmap_bivariate.png", width = 8.29, height = 4.21, dpi = 600)