How to divide a complex geometry into n equals subgeometries?

Question

I would like to divide a complex geometry into n subgeometries of the same area.

For instance if I have a rectangle I can do something like this

from shapely.geometry import LineString, MultiPolygon, Polygon
from shapely.ops import split
def splitPolygon(polygon, nx, ny):
    minx, miny, maxx, maxy = polygon.bounds
    dx = (maxx - minx) / nx
    dy = (maxy - miny) / ny

    minx, miny, maxx, maxy = polygon.bounds
    dx = (maxx - minx) / nx  # width of a small part
    dy = (maxy - miny) / ny  # height of a small part
    horizontal_splitters = [LineString([(minx, miny + i*dy), (maxx, miny + i*dy)]) for i in range(ny)]
    vertical_splitters = [LineString([(minx + i*dx, miny), (minx + i*dx, maxy)]) for i in range(nx)]
    splitters = horizontal_splitters + vertical_splitters
    result = polygon
    for splitter in splitters:
        result = MultiPolygon(split(result, splitter))
    return result

myPolygons = splitPolygon(polygon, 5, 5)
import geopandas as gpd
gdfR   = gpd.GeoDataFrame(columns=['geometry'], data=myPolygons.geoms)
f,ax=plt.subplots()
gdfR.boundary.plot(ax=ax, color='red')
polygon.boundary.plot(ax=ax)

I would like to split a complex geometry like following one into n smallest geometries of the same area. Is possible to download the geometry as shapefile here.

Just speaking my mind here: Find the isomorphism that preserves the distance ratios between points from the square to the target geometry? you can frame it as a constrained optimization problem? — Mario Viti, Jul 28 '21 at 08:22
@MarioViti I really do not understand your question. I would like to divide my area in smallest regions, something similar to a Voronoi tessellation — emax, Jul 28 '21 at 08:33
Algorithms like Voronoi Tesselation are an exact solution to an optimization problem. In your case Voronoi (or its dual Delaunay) do not satisfy your constraints of equal areas. So I suggested to frame it as an optimization problem and depending on the shape of the resulting problem using the adequate solver (iterative almost allways work). Check out Voronoi + Lloyd iteartive approach https://www.semanticscholar.org/paper/Lloyd-relaxation-using-analytical-Voronoi-diagram-L-Mouton-B%C3%A9chet/d8c702e0014458dad8ffb04c9c76228a74e90f7a — Mario Viti, Jul 28 '21 at 08:47
There are some github repos e.g. https://github.com/duhaime/lloyd But I haven't tested them :( — Mario Viti, Jul 28 '21 at 09:18

Tim · Answer 1 · 2021-08-14T20:28:35.533

One of the simplest and most efficient solutions is to take the minimum_rotated_rectangle of the geometry.

First, you transform the Polygon or MultiPolygon into a LineString and you get the minimum_rotated_rectangle containing this line.
Then you apply your previous code to this rectangle.
Finally you have to make a second split between the LineString of your polygon and the subgeometries.

A great turnkey solution: https://stackoverflow.com/a/68778560/14864907

How to divide a complex geometry into n equals subgeometries?

1 Answers1