Interpolating and plotting the votes for the BoerBurgerBeweging in Utrecht

bbb_votes.R

# Plot & interpolate BBB vote percentage in Utrecht
# using open data from data.utrecht.nl
# last edited 20230320 by 
library(tidyverse)
library(sf)
library(curl)
library(glue)
library(jsonlite)
library(pbapply)
library(ggspatial)
library(gstat)
library(firatheme)

# Download the raw data
data_url <- "https://data.utrecht.nl/sites/default/files/open-data/osv4-3_telling_ps2023_utrecht_gemeente_utrecht.csv"
df <- read_delim(data_url, delim = ";", skip = 5)

# Geocode locations using openstreetmaps ----
# Get stemlokaal locations from raw data
locations <- 
  df[2,] |> 
  select(-Lijstnummer:-Totaal) |> 
  pivot_longer(
    cols = everything(), 
    names_to = "stemlokaal", 
    values_to = "postcode"
  ) |>
  mutate(
    adres = str_extract(stemlokaal, "(^.*? [0-9]+)"), 
    .keep = "unused"
  ) |> 
  distinct()

# Function to geocode
geo_code <- function(adres, postcode) {
  adres_esc <- curl_escape(adres)
  postc_esc <- curl_escape(postcode)
  result <- 
    glue("https://nominatim.openstreetmap.org/search?street={adres_esc}&postalcode={postc_esc}&limit=1&format=json") |> 
    read_lines(progress = FALSE) |> 
    fromJSON()
  
  if (length(result) == 0L) return(list("lat" = NA, "lon" = NA))
  
  result |> select("lat", "lon")
}

# Perform geocoding (with progress bar)
latlon <- pbsapply(
  X = 1:nrow(locations), 
  FUN = function(i) geo_code(locations[i, "adres"], locations[i, "postcode"])
)

# Add geocode to locations
locations$lat <- as.numeric(latlon[1,])
locations$lon <- as.numeric(latlon[2,])

# Get BBB votes per location ----
# get votes per party per location
party_votes <-
  df[19:nrow(df),] |> 
  filter(!is.na(Aanduiding)) |> 
  select(-Lijstnummer, -Volgnummer, -`Naam kandidaat`, -Totaal) |> 
  pivot_longer(-Aanduiding, names_to = "stemlokaal") |>
  mutate(adres = str_extract(stemlokaal, "(^.*? [0-9]+)")) |> 
  mutate(value = parse_integer(value)) |> 
  group_by(adres, Aanduiding) |> 
  summarize(votes = sum(value))

# compute percentage BBB per location
bbb_votes <- 
  party_votes |> 
  group_by(adres) |> 
  mutate(perc = votes / sum(votes) * 100) |> 
  filter(Aanduiding == "BBB") |> 
  select(adres, perc)


# Kriging and mapping ----
# create spatial dataset with the right CRS
geo_votes <- 
  left_join(bbb_votes, locations) |>
  ungroup() |> 
  st_as_sf(coords = c("lon","lat"), crs = 4326) |> 
  st_transform(28992)

# create prediction grid for kriging
grid_sfc <- st_make_grid(st_buffer(geo_votes, 500), cellsize = 30)

# Create an empirical variogram 
vdist <- variogram(perc ~ 1, geo_votes, cutoff = 2000, width = 30)

# Estimate a spherical model for the variogram
mdist <- fit.variogram(vdist, vgm(6, "Gau", 2000, 1), fit.method = 6)


# kriging: generate predictions for grid
grid_pred <- krige(
  formula   = perc ~ 1,
  locations = geo_votes, 
  newdata   = st_centroid(grid_sfc), 
  model     = mdist
)

pred_plot <- 
  grid_pred |> 
  mutate(geometry = grid_sfc) |> 
  ggplot() +
  annotation_map_tile(zoom = 14) +
  geom_sf(aes(fill = var1.pred), colour = NA) +
  geom_sf(data = geo_votes, aes(size = perc), colour = "#343434") +
  scale_size_continuous(
    range = c(.8, 4), 
    breaks = c(0, 5, 10, 15, 20, 25, 30),
    labels = c(0, 5, 10, 15, 20, 25, 30)
  ) +
  scale_fill_gradient(
    low = "#93c01f00", 
    high = "#93c01fff", 
    guide = "none", 
  ) +
  labs(
    title = "BoerBurgerBeweging in Utrecht",
    subtitle = "Verschil tussen binnenstad en buitengebied",
    size = "% BBB"
  ) + 
  theme_fira() +
  theme(
    axis.text = element_blank(), 
    axis.line = element_blank(), 
    panel.grid.major = element_blank()
  )

ggsave("bbb_plot.png", plot = pred_plot, width = 12, height = 8, bg = "white")