Constructing a Voronoi diagram in R with data from OpenStreetMap

October 25, 2022

For this example, I am trying to color the map of the UK village of Wetwang based on proximity to one of the two churches in the village. Here is the code for the basic plot without Voronoi diagrams, which shows the polygon of the village and the two churches as points.

library(osmdata)
library(sf)
library(tidyverse)

bb <- getbb("Wetwang", featuretype = "settlement", format_out = "polygon")

wetwang <- getbb("Wetwang", featuretype = "settlement") %>% 
  opq() %>% 
  add_osm_feature(key = "boundary", value = "administrative") %>% 
  osmdata_sf() %>% 
  trim_osmdata(bb)

churches <- getbb("Wetwang", featuretype = "settlement") %>% 
  opq() %>% 
  add_osm_feature(key = "building", value = "church") %>% 
  osmdata_sf() %>% 
  trim_osmdata(bb)

ggplot() +
  geom_sf(data = wetwang$osm_multipolygons) +
  geom_sf(data = churches$osm_polygons %>% 
            st_centroid()) +
  theme_void()

This is the output:

I am trying to use st_voronoi to construct these Voronoi diagrams. However, it doesn’t seem to be working:

envelope <- st_cast(wetwang$osm_multipolygons, "POLYGON")

church_areas <- churches$osm_polygons %>% 
  st_centroid() %>% 
  st_union() %>% 
  st_voronoi(envelope = envelope)

This fails because

Error in vapply(lst, class, rep(NA_character_, 3)) : 
  values must be length 3,
 but FUN(X[[1]]) result is length 2

Without passing the envelope argument, I get a map where the diagram doesn’t extend over the whole area of Wetwang:

church_areas <- churches$osm_polygons %>% 
  st_centroid() %>% 
  st_union() %>% 
  st_voronoi() %>% 
  st_collection_extract(., "POLYGON") %>% 
  st_as_sf()

ggplot() +
  geom_sf(data = wetwang$osm_multipolygons) +
  geom_sf(data = churches$osm_points) +
  geom_sf(data = church_areas) +
            theme_void()

So what am I doing wrong?

>Solution :

One possible solution:

# create sf objects to use. churches$osm_points has many points for each church
#  finding the centroid of churches$osm_polygons solves this
churches_single_point <- churches$osm_polygons %>% st_centroid()
churches_polygon <- churches$osm_polygons
wetwang_sf <- wetwang$osm_multipolygons

church_voronoi <- churches_single_point %>%
  st_geometry() %>%
  st_union() %>%
  st_voronoi(envelope = st_geometry(wetwang_sf)) %>%
  st_cast() %>%
  st_as_sf()

# use st_intersection to "crop" the voronoi & join it with churches_single_point
#  plot & fill by name.y
church_voronoi %>% 
  st_intersection(wetwang_sf) %>%
  st_join(churches_single_point) %>%
  ggplot() + 
  geom_sf(aes(fill = name.y))