Perimeter of Voronoi Cells

I am working with the R programming language. I have the following dataset:

# simulate data
set.seed(123)

df1 <- data.frame(longitude = rnorm(100,10,1),
                 latitude = rnorm(100,10,1), color = "red")

df2 <- data.frame(longitude = rnorm(100,8,1),
                 latitude = rnorm(100,8,1), color = "blue")

df = rbind(df1, df2)

ggplot(df, aes(x=longitude, y=latitude, color=color)) +
  geom_point()

Using this data, I made Voronoi Cells for this data (based on the colors):

library(ggvoronoi) 

ggplot(df, aes(x = longitude, y = latitude, fill = color)) +   geom_voronoi() +   geom_point() +   scale_fill_manual(values = c("red", "blue"))

My Question: It it somehow possible to find out the perimeter of the red and blue shapes - and then store these perimeters in a shapefile?

Based on some references I read, I tried to adapt them for this problem - but I am not sure if this is correct?

library(ggvoronoi)
library(sf)

vor_spdf <- voronoi_polygon(data = df, x = "longitude", y = "latitude")
vor_df <- fortify_voronoi(vor_spdf)

vor_sf <- st_as_sf(vor_spdf)

perimeter <- st_length(vor_sf)

vor_sf$perimeter <- perimeter

st_write(vor_sf, "voronoi.shp")

Can someone please show me how to do this correctly?

Thanks!

sf provides st_voronoi() function, an interface for GEOS implementation for creating Voronoi tesselations so that you can work with sf / sfc / sfg objects right from the start:

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)

# simulate data
set.seed(123)
df1 <- data.frame(longitude = rnorm(100,10,1),
                  latitude = rnorm(100,10,1), color = "red")
df2 <- data.frame(longitude = rnorm(100,8,1),
                  latitude = rnorm(100,8,1), color = "blue")
df = rbind(df1, df2)

# create sf from data.frame
points_sf <- st_as_sf(df, coords = c("longitude", "latitude"))

# create Voronoi tesselation from union of points, extract polygons, 
# turn geometry set to sf
voronoi_sf <-
  points_sf |> 
  st_union() |> 
  st_voronoi() |> 
  st_collection_extract("POLYGON") |>
  st_sf(geometry = _)

# chck bbox
st_bbox(voronoi_sf)
#>      xmin      ymin      xmax      ymax 
#> -1.463466 -2.172836 19.894271 20.947978

# create a fake boundary polygon around points for more cmopact output
fake_boundary_sf <- 
  points_sf |>
  st_union() |>
  st_convex_hull() |>
  st_buffer(1, joinStyle = "MITRE") 

# to use more compact shape instead of default envelope of st_voronoi,
# cut resulting sf object into a shape of fake_boundary_sf
voronoi_fbound_sf <- st_intersection(voronoi_sf, fake_boundary_sf)

# spatial join to add color attribute of points_sf  
voronoi_fbound_sf <- st_join(voronoi_fbound_sf, points_sf)

ggplot() +
  geom_sf(data = voronoi_sf, fill = NA, color = "grey") +
  geom_sf(data = voronoi_fbound_sf, aes(fill = color), alpha = .4) +
  geom_sf(data = points_sf, aes(color = color), size = .5 ) +
  scale_fill_identity() +
  scale_color_identity() +
  coord_sf(xlim = c(5,15), ylim = c(5,15)) +
  theme_minimal()

Resulting sf object:

voronoi_fbound_sf
#> Simple feature collection with 200 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5.252011 ymin: 4.533196 xmax: 13.14487 ymax: 14.26178
#> CRS:           NA
#> First 10 features:
#>    color                       geometry
#> 1   blue POLYGON ((6.143967 7.55766,...
#> 2   blue POLYGON ((7.407921 12.20419...
#> 3   blue POLYGON ((6.664117 6.580956...
#> 4   blue POLYGON ((6.664117 6.580956...
#> 5   blue POLYGON ((6.583638 8.390378...
#> 6   blue POLYGON ((6.607826 7.792981...
#> 7   blue POLYGON ((6.908212 9.287342...
#> 8   blue POLYGON ((7.158668 10.04188...
#> 9   blue POLYGON ((6.583638 8.390378...
#> 10  blue POLYGON ((6.304645 9.325164...

Write and read:

st_write(voronoi_fbound_sf, "voronoi_fbound.shp")
#> writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRS
#> Writing layer `voronoi_fbound' to data source 
#>   `voronoi_fbound.shp' using driver `ESRI Shapefile'
#> Writing 200 features with 1 fields and geometry type Polygon.
st_layers("voronoi_fbound.shp")
#> Driver: ESRI Shapefile 
#> Available layers:
#>       layer_name geometry_type features fields
#> 1 voronoi_fbound       Polygon      200      1
#>                                    crs_name
#> 1 Undefined Cartesian SRS with unknown unit
read_sf("voronoi_fbound.shp") |> plot(pal =c("blue", "red"))

Edit: dissolve individual Vorogonoi polygons by color attribute into (2) multipolygons; optionally transform to (2) polygons (keep largest sub-poly + remove holes)

# union polygons by "color", 200 polygons are turned into 2 MULTIPOLYGONs (polygon collections)
voronoi_multipolygon <- 
  voronoi_fbound_sf |>
  group_by(color) |>
  summarise() 

voronoi_multipolygon
#> Simple feature collection with 2 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5.252011 ymin: 4.533196 xmax: 13.14487 ymax: 14.26178
#> Projected CRS: Undefined Cartesian SRS with unknown unit
#> # A tibble: 2 × 2
#>   color                                                                 geometry
#>   <chr>                                                           <MULTIPOLYGON>
#> 1 blue  (((9.402573 4.904994, 9.394385 4.901992, 8.721431 4.7367, 7.892898 4.53…
#> 2 red   (((13.06049 8.0261, 11.99858 6.95691, 11.61793 6.573656, 10.37808 7.740…
plot(voronoi_multipolygon, pal =c("blue", "red"))

# split MULTIPOLYGONs to POLYGONs ("explode", 2 rows to 10),
# group by color, from each group select polygon with greatest area 
# (resulting polygons include holes),
# remove holes
voronoi_main_no_holes <-
  voronoi_multipolygon |>
  st_cast("POLYGON") |>
  group_by(color) |>
  slice_max(st_area(geometry)) |>
  sfheaders::sf_remove_holes()

voronoi_main_no_holes
#> Simple feature collection with 2 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5.252011 ymin: 4.533196 xmax: 13.14487 ymax: 14.26178
#> Projected CRS: Undefined Cartesian SRS with unknown unit
#> # A tibble: 2 × 2
#> # Groups:   color [2]
#>   color                                                                 geometry
#> * <chr>                                                                <POLYGON>
#> 1 blue  ((9.402573 4.904994, 9.394385 4.901992, 8.721431 4.7367, 7.892898 4.533…
#> 2 red   ((13.06049 8.0261, 11.99858 6.95691, 11.61793 6.573656, 10.37808 7.7400…
plot(voronoi_main_no_holes, pal =c("blue", "red"))

^{Created on 2023-08-08 with reprex v2.0.2}