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.

Answer 1

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)