Create voronoi cells within each polygon seperately

Question

Data: In the data below I have clusters, which are 2 large groupings of the data. Within each cluster are 5 districts. I use the points within each cluster to create a polygon for the cluster.

Problem: I'm attempting to calculate voronoi for each district within each cluster. So each of the 2 cluster polygons should have 5 voronoi cells within it. How can I create 5 voronoi cells bounded by each cluster polygon?

library(dbplyr)
library(sf)
library(purrr)
library(concaveman)

# Create raw data
df <- data.frame(
  "cluster" = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2),
  "district" = c('1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_1','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_2','1_3','1_3','1_3','1_3','1_3','1_3','1_3','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_4','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','1_5','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_1','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_2','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_3','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_4','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5','2_5'),
  "mx" = c(0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,0.973009090909091,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,1.10794444444444,0.983014285714286,0.983014285714286,0.983014285714286,0.983014285714286,0.983014285714286,0.983014285714286,0.983014285714286,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.20296666666667,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,1.31765526315789,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.560226315789474,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.377593333333333,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.7201725,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.215625,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571,0.369878571428571),
  "my" = c(-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.27477272727273,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.49177777777778,-1.07938571428571,-1.07938571428571,-1.07938571428571,-1.07938571428571,-1.07938571428571,-1.07938571428571,-1.07938571428571,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.0937,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.15119473684211,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.02930526315789,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-1.10028666666667,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.9094475,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.642783333333333,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286,-0.795914285714286),
  "ID" = 1:200,
  "X" = c(1.0083,0.9068,1.0232,1.0005,0.8388,0.8655,1.0133,1.0106,1.0139,1.0537,0.8759,1.0063,1.0187,1.0241,1.0004,0.8886,1.0803,0.8518,0.9998,1.0154,0.8851,1.0252,1.0926,1.064,1.1015,1.0354,1.1511,1.1074,1.202,1.1063,1.0173,1.137,1.1156,1.0776,1.1315,1.2281,0.9974,1.1487,1.0098,1.2197,0.9598,0.9695,0.9268,1.0008,1.0827,0.9331,1.0084,1.1311,1.1856,1.1932,1.2464,1.2331,1.1944,1.1846,1.2203,1.1753,1.2245,1.2396,1.2682,1.2359,1.1918,1.1205,1.3166,1.3072,1.2984,1.2842,1.3451,1.3197,1.2996,1.3482,1.28,1.3051,1.4471,1.315,1.3177,1.3387,1.32,1.3508,1.2747,1.3681,1.2735,1.325,1.3093,1.3244,1.2626,1.3123,1.2819,1.2619,1.3639,1.3099,1.2783,1.3313,1.3895,1.3559,1.3521,1.2589,1.4204,1.2303,1.2808,1.3125,0.55,0.5483,0.6329,0.5006,0.5233,0.566,0.5673,0.4864,0.5658,0.5636,0.542,0.6146,0.5648,0.6092,0.5329,0.5726,0.574,0.62,0.51,0.3787,0.3769,0.3579,0.3881,0.3806,0.403,0.3962,0.409,0.3422,0.4361,0.3853,0.3568,0.3105,0.3674,0.3752,0.6694,0.6492,0.6568,0.6773,0.6237,0.7265,0.7201,0.7596,0.7049,0.7699,0.6555,0.7105,0.731,0.6376,0.8865,0.754,0.7983,0.699,0.7223,0.7214,0.6496,0.7907,0.7418,0.7825,0.7417,0.8978,0.7875,0.6874,0.7761,0.6189,0.706,0.7037,0.7149,0.7059,0.687,0.7888,0.6514,0.7271,0.6679,0.7067,0.2631,0.2701,0.0822,0.069,0.2196,0.2848,0.2661,0.2343,0.2905,0.0684,0.2874,0.252,0.3301,0.5509,0.3343,0.343,0.2951,0.2524,0.5442,0.3187,0.3143,0.3731,0.3352,0.2711,0.455,0.4609),
  "Y" = c(-1.2547,-1.2297,-1.3071,-1.237,-1.362,-1.3776,-1.2552,-1.2354,-1.2396,-1.3493,-1.3019,-1.2484,-1.2435,-1.233,-1.217,-1.2715,-1.3396,-1.34,-1.233,-1.2511,-1.3333,-1.1851,-1.5461,-1.5043,-1.5452,-1.5412,-1.4937,-1.5425,-1.4155,-1.523,-1.549,-1.5077,-1.5458,-1.369,-1.5033,-1.3805,-1.5183,-1.4288,-1.5429,-1.3952,-1.0349,-1.0838,-1.0615,-1.0702,-1.0446,-1.1367,-1.124,-1.0626,-1.0958,-1.0808,-1.0775,-1.0499,-1.0963,-1.0341,-1.0348,-1.0838,-1.162,-1.0487,-1.0924,-1.1537,-1.2107,-1.1224,-1.1499,-1.1803,-1.2877,-1.1151,-1.1339,-1.1431,-1.1521,-1.1675,-1.1407,-1.1916,-1.1229,-1.1308,-1.1154,-1.183,-1.1214,-1.0793,-1.1857,-1.0679,-1.1633,-1.075,-1.1354,-1.1494,-1.162,-1.1582,-1.18,-1.1234,-1.1077,-1.1144,-1.1305,-1.1482,-1.2058,-1.1685,-1.2152,-1.1439,-1.1252,-1.2113,-1.1632,-1.1965,-0.9917,-0.9861,-1.064,-0.9898,-0.9799,-1.1048,-1.0085,-0.9395,-1.0425,-1.0806,-1.0132,-1.0785,-1.1109,-1.0632,-0.945,-1.009,-1.1005,-1.0759,-0.9732,-1.1279,-1.1746,-1.1957,-1.0738,-0.9896,-1.0601,-0.9735,-1.0953,-1.0849,-1.0406,-1.1202,-1.0781,-1.2002,-1.157,-1.1328,-0.7416,-0.9323,-0.9372,-0.7013,-0.9191,-0.9356,-0.9838,-0.9302,-0.9407,-1.044,-0.6478,-0.9147,-0.9688,-0.9272,-0.9494,-1.0817,-0.9915,-0.9329,-0.8514,-0.9665,-0.783,-0.9601,-0.8996,-0.7717,-1.0067,-0.9839,-1.0594,-0.9705,-1.01,-0.9163,-1.0049,-0.6829,-0.7918,-0.7353,-0.9295,-0.9944,-0.9524,-0.9257,-0.936,-0.7661,-0.5825,-0.5989,-0.7375,-0.7262,-0.594,-0.6145,-0.571,-0.7069,-0.6377,-0.7865,-0.5962,-0.5615,-0.7729,-0.7873,-0.7713,-0.7774,-0.816,-0.8865,-0.7689,-0.8591,-0.7913,-0.7231,-0.7859,-0.8542,-0.7447,-0.8042)
)

#Create cluster polygons from data
sf <- st_as_sf(df, coords = c("X","Y"), crs = 4326)
shapes <- map(unique(sf$cluster),
              ~ concaveman(sf[sf$cluster %in% .,])
) %>% 
  map2(unique(sf$cluster), ~ mutate(.x, cluster = .y)) %>% 
  reduce(rbind)
  
###################################################
# OK, here is where I start running into problems #
###################################################

# Attempt to calculate voronoi cells within each cluster polygon
bbox_polygon <- function(x) {
  bb <- sf::st_bbox(x)
  
  p <- matrix(
    c(bb["xmin"], bb["ymin"], 
      bb["xmin"], bb["ymax"],
      bb["xmax"], bb["ymax"], 
      bb["xmax"], bb["ymin"], 
      bb["xmin"], bb["ymin"]),
    ncol = 2, byrow = T
  )
  
  sf::st_polygon(list(p))
}

sf_district <- df %>%
  select(district, mx, my) %>%
  st_as_sf(coords = c("mx","my"), crs = 4326)
sfbox <- st_sfc(bbox_polygon(sf_district))

v <- st_voronoi(st_union(sf_district), sfbox)

# Plot to see what it looks like
plot(st_intersection(st_cast(v), st_union(shapes)), col = 0)

If you look at the picture, you can see that I'm doing something wrong because the northwest polygon has 6 voronoi cells when there are only 5 districts in the data. I think the problem is that my script is just creating one big voronoi tesselation and then laying the polygon shapes on top of it, rather than calculating the voronoi for each cluster separately and then bounding them within each cluster polygon.

See the minimal part of the minimal reproducible example guidance. If all you need help with is getting the voronoi calculation by group, just give us a representative sample of the data at that point—no need for us to have the packages necessary to plot or compute hulls if that's not specifically part of the problem. Since sf works with dplyr verbs, can you just use a group_by at some point? — camille
@camille, thank you for the direction. I've simplified my example as much as I could. The first half is just setting up the data. I then put a comment in where I run into trouble. I want to loop through each cluster and calculate the voronoi for that cluster, bounding the voronoi to the cluster polygon. So I need to calculate the voronoi for each cluster AND calculate a bounding box for each cluster using the "shapes" geometry. — Obed

Michael Dorman Michael Dorman · Accepted Answer · 2020-11-08T20:28:43

To get separate voronoi polygons for each "cluster", you can run a for loop. Instead of the expression v <- st_voronoi(st_union(sf_district), sfbox), as follows:

sf_district = st_join(sf_district, shapes)
result = list()
for(i in unique(sf_district$cluster)) {
    u = sf_district[sf_district$cluster == i, ]
    v = st_voronoi(st_union(u), sfbox)
    v = st_intersection(st_cast(v), shapes[shapes$cluster == i, ])
    result[[i]] = v
}
result = do.call(c, result)

This ensures that the voronoi polygons in each cluster are unaffected by other points, and that each cluster has as many polygons as there are points.

Here is a plot of the result:

plot(result, col = hcl.colors(length(result), "Set 2"))
plot(st_geometry(sf_district), add = TRUE)

Create voronoi cells within each polygon seperately

1 Answers