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)
[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: enter image description here

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


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.


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)
Units: [m^2]
[1] 2735339585 2880268727

polygon cutted by half


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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.