Dividing Polygon into parts which have equal area using R

Question

I have a Multipolygon and want to split it into parts with an equal area.

How can I do it using R?

my data:

pol <-readOGR("/Users/Desktop/test.shp")
pol_wkt <- wicket::sp_convert(pol)
print(pol_wkt)
[1] "MULTIPOLYGON(((8.23806 48.6899,8.25024 48.6919,8.25993 48.6867,8.25993 48.6867,8.27037 48.6693,8.25819 48.6638,8.28901 48.6504,8.34146 48.6486,8.33624 48.6185,8.26142 48.6203,8.23259 48.6551,8.23806 48.6899)),
((8.20773 48.6031,8.25247 48.6041,8.25247 48.6041,8.31959 48.5932,8.33922 48.5688,8.24427 48.5731,8.23134 48.5542,8.28578 48.5539,8.31089 48.5186,8.23681 48.5221,8.20773 48.6031)))"

how the polygon looks like:

I want divide them into n parts, so that the parts are near to an equal area

I answered this for sp objects here using rgeos and a root-finding algorithm to find the level to split a polygon into equal (or variable) area slices: https://groups.google.com/g/r-help-archive/c/V87xEntSpQc?pli=1 - if anyone cares to convert that to sf.... — Spacedman, Sep 29 '20 at 18:51

score 5 · Accepted Answer · edited Jan 28 '24 at 17:00

As I needed something similar, I wrote a small function for R following the steps of this question on "Splitting polygon into equal area polygons using QGIS" I hope this helps somebody, even though this is a late answer. This function takes an sf polygon:

library(sf)
library(mapview)
library(tidyverse) 
library(dismo)
library(osmdata)  
library(mapview)
split_poly <- function(sf_poly, n_areas) {
    # Create random points
    points_rnd <- st_sample(sf_poly, size = 10000)
    # k-means clustering
    points <- do.call(rbind, st_geometry(points_rnd)) %>%
        as_tibble() %>% setNames(c("lon","lat"))
    k_means <- kmeans(points, centers = n_areas)
    # Create voronoi polygons
    voronoi_polys <- dismo::voronoi(k_means$centers, ext = sf_poly)
    # Clip to sf_poly
    crs(voronoi_polys) <- crs(sf_poly)
    voronoi_sf <- st_as_sf(voronoi_polys)
    equal_areas <- st_intersection(voronoi_sf, sf_poly)
    equal_areas$area <- st_area(equal_areas)
    return(equal_areas)
}

Using the same data as @Elio Diaz earlier in this thread:

pol <- osmdata::getbb("aguascalientes", format_out = "sf_polygon") 
pol_areas <- split_poly(pol, 20)
mapview(pol_areas)

The areas are not completely equal, but it is close enough for my case. See pol_areas$area for m² per polygon.

score 2 · Answer 2 · answered Sep 29 '20 at 17:38

This is a script that approximates the fractions, it has a great field for optimization. It only does horizontal cutting, not in an oriented bounding box. In the porcientos argument you may put as many values as you like, it is not only for halves (c(.5,.5), this means c(.4, .3, .2, .1) would be a valid vector as well.

library(units)
library(sf)
library(dplyr)
library(osmdata)
pol <- osmdata::getbb("aguascalientes", format_out = "sf_polygon") 
porcientos <- c(.5,.5) # the half argument
polycent <- function(poly, porcientos) {
  df   <- st_sf(id = 1:length(porcientos), crs = 4326, # empty sf for populating
                geometry = st_sfc(lapply(1:length(porcientos), function(x) st_multipolygon())))
  area1 <- st_area(poly) %>% sum() # st_area takes multipolygons as one; # area1 is constant
  poly2 <- poly # duplicating for the final cut
  for(j in  seq_along(porcientos[-length(porcientos)])) { 
    bb = st_bbox(poly2)
    top <- bb['ymax']
    bot <- bb['ymin']
    steps <- seq(bot, top, by = (top - bot) / 80)
    for(i in steps[2:length(steps)]) {  # 2:n because 1:n renders a line "ymax" = "ymin"
      bf <- bb
      bf['ymax'] = i
      temp <- st_intersection(poly, st_as_sfc(bf, 4326))
      area2 <- st_area(temp) %>% sum()           # con get(.., i) coz st_area prints rounded
      if(drop_units(area2)/drop_units(area1) >= porcientos[j]) break
      df$geometry[j] <- st_geometry(temp)
    }
    poly2 <- st_difference(poly2, st_union(df))
  }
  df$geometry[length(porcientos)] <- st_geometry(st_difference(poly, st_union(df)))
  poly <- df
}
ea = polycent(pol, porcientos)
plot(rbind(ea[1,], ea[2,]), graticule = T, axes = T)
st_area(ea)
Units: [m^2]
[1] 2735339585 2880268727

robin loche · Answer 3 · 2020-09-29T18:07:33.323

From my understanding, this is a very difficult question.

You can divide your geometries into a number of vertices max easily, which most of the time can look like that, with for exemple this function, but I'm guessing this is not really for an objective of performance (this function is mainly used to simplify computation).

If it's really the area that you care about, I don't think there is an exact algorithm to do that, but you can get approximation. What you can do is to rasterize your geometries (see for exemple here), then count the number of pixel that you have, and divide it into the number that you want, and transform it back to polygon. This is not perfect but you can control the precision by controlling the size of the pixels while rasterizing, and you can control the form of your sub-part with the way you select the pixels in each sub-part.

The result will probably not be very "organic" (it will look squared) depending on how you select the pixels that goes into each sub-part, but doing differently is probably a little more difficult.

Maybe, if you want so, you can try to generate a lot of random points uniformly inside your polygons and try to apply a clustering algorithm like k-means on it? (if you have multiple parts separated, like in your exemple, and you want a small number of sub-parts, it will maybe not be what you're looking for, because the clustering will probably always separate the polygons first)

Dividing Polygon into parts which have equal area using R

3 Answers3

Linked