R – How to Divide a Polygon into Equal Area Parts

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

Answer 1

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