r - 使用来自 OpenStreetMap 的数据在 R 中构建 Voronoi 图

Question

在此示例中，我尝试根据与村中两座教堂之一的接近程度为英国韦特旺村的 map 着色。这是没有 Voronoi 图的基本图的代码，它将村庄和两座教堂的多边形显示为点。

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()

这是输出:

我正在尝试使用 st_voronoi 来构建这些 Voronoi 图。但是，它似乎不起作用:

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

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

这会失败，因为

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

如果不传递信封参数，我会得到一张 map ，其中图表并未延伸到 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()

那么我做错了什么？

Answer 1

一种可能的解决方案:

# 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))