Dividing Polygon into parts which have equal area using R

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

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)