如何将质心（标记）的数据指定给它所属的 Voronoi/Thiessen 多边形？（R）

Question

在他的文章中，Kyle Walker展示了一种在Leaflet中制作Voronoi多边形的方法。他通过以下代码绘制了Fort Worth每家星巴克咖啡馆周围的Voronoi多边形：

library(leaflet);    library(rgeos)
library(rgdal);    library(spatstat)
library(maptools)

starbucks <- read.csv('starbucks.csv')
fw <- subset(starbucks, City == 'Fort Worth')
coords <- cbind(fw$Longitude, fw$Latitude)
## Spatial points w/the WGS84 datum
sp_fw <- SpatialPointsDataFrame(coords = coords, data = fw, 
                proj4string = CRS("+proj=longlat +datum=WGS84"))
sp_fw_proj <- spTransform(sp_fw, CRS("+init=epsg:26914"))
fw_coords <- sp_fw_proj@coords
## Create the window for the polygons
window <- owin(range(fw_coords[,1]), range(fw_coords[,2]))
## Create the polygons
d <- dirichlet(as.ppp(fw_coords, window))
## Convert to a SpatialPolygonsDataFrame and calculate an "area" field.  
dsp <- as(d, "SpatialPolygons")
dsp_df <- SpatialPolygonsDataFrame(dsp,
                                   data = data.frame(id =  1:length(dsp@polygons)))
proj4string(dsp_df) <- CRS("+init=epsg:26914")
dsp_df$area <- round((gArea(dsp_df, byid = TRUE) / 1000000), 1)
dsp_xy <- spTransform(dsp_df, CRS("+proj=longlat +datum=WGS84"))
## Map it!
leaflet() %>%
  addMarkers(data = fw, 
         lat = ~ Latitude, 
         lng = ~ Longitude, 
         popup = fw$Name) %>%
  addPolygons(data = dsp_xy, 
          color = "green", 
          fill = "green", 
          popup = paste0("Area: ", 
                         as.character(dsp_xy$area), 
                         " square km")) %>%
    addTiles()

我想给他的地图添加一个额外的功能：我想为一个多边形分配一个特定的颜色。该颜色取决于最近标记（质心）的特征。

例如，用星巴克的质心“绿色”着色每个多边形，并用唐恩·肯德基的质心“紫色”着色。（假设starbucks.csv也包括唐恩·肯德基的坐标）

换句话说，我想将质心（“fw”）的数据与其所属的多边形（“dsp_xy”）的数据合并。

有人能帮我解决这个问题吗？

- Zigaar

---

我们不是来帮你完成作业的。你目前尝试了什么？你卡在哪里了？具体来说，你遇到了什么问题？ - Adam

可能在这里 sp::over 会有用。 - user3603486

Answer 1

“voronoi”函数来自“dismo”包，适用于编程相关内容。我还将使用此帖子演示R的新“sf”包。

让我们生成一个可重复的虚假数据集，包含星巴克和唐恩都乐的位置：

library(leaflet)
library(sf)
library(dismo)
library(sp)

set.seed(1983)

# Get some sample data

long <- sample(seq(-118.4, -118.2, 0.001), 50, replace = TRUE)

lat <- sample(seq(33.9, 34.1, 0.001), 50, replace = TRUE)

type <- sample(c("Starbucks", "Dunkin"), 50, replace = TRUE)

接下来，让我们从数据中创建一个名为sf的数据框，并查看它：

points <- data.frame(long = long, lat = lat, type = type) %>%
  st_as_sf(crs = 4326, coords = c("long", "lat"))

plot(points)

接下来，我们使用dismo包中的voronoi函数创建Voronoi多边形，非常简单，然后将其与我们的点具有相同的坐标系。在实际工作流程中，您应该使用投影坐标系，但我只是为了说明而使用WGS84（这些操作将假定为平面）。同时请注意，我在sf和sp类之间来回切换；R世界将随着时间的推移完全支持sf，但在过渡期间强制转换很简单。

polys <- points %>%
  as("Spatial") %>%
  voronoi() %>%
  st_as_sf() %>%
  st_set_crs(., 4326)

plot(polys)

现在，使用Leaflet将其可视化，并使用您想要的颜色：

pal <- colorFactor(c("purple", "green"), polys$type)

polys %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(fillColor = ~pal(type), weight = 0.5, color = "grey") %>%
  addCircleMarkers(data = points, label = ~type, color = ~pal(type))

我们这里不需要它，但是你也需要了解 sf 中的一个函数 st_join，它可以无缝处理空间连接，并且适用于您最初提出的覆盖类型。

- kwalkertcu