Algorithm to subdivide a polygon in smaller polygons

Question

I have a polygon made of successive edges on a plane, and would like to subdivide it in sub-polygons being triangles or rectangles. Where can I find an algorithm to do this ? Thanks !

Answer 1

shapely.delaunay_triangles

import shapely
import matplotlib.pyplot as plt
from shapely.plotting import plot_polygon

#Create a polygon geometry from a WKT string
wkt_polygon = 'Polygon ((515931 6429487, 515612 6430222, 516062 6430294, 516140 6430782, 516491 6430997, 517682 6430470, 517077 6429356, 516869 6428706, 516634 6428478, 516237 6428855, 516413 6429695, 515931 6429487))'
poly = shapely.wkt.loads(wkt_polygon)

#Triangulate it
triangulated = shapely.delaunay_triangles(geometry=poly)
triangulated = list(triangulated.geoms) #From a geometry collection to a list of polygons

#Plot
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10,5))
plot_polygon(polygon=poly, ax=ax[0])
for triangle in triangulated:
    plot_polygon(polygon=triangle, ax=ax[1])

Answer 2

Stumbled across this after many searches.

Thanks @Aditya Chhabra for your submission, it works great but get_squares_from_rect is very slow for small side lengths due to iterative clips.

We can do this instantaneously if we combine all LineStrings into a single collection, then clip and polygonize in one step, which I found in in this question.

Previously side lengths of 0.0001 (EPSG:4326) took > 1 minute, now it takes no time.

from shapely.ops import unary_union, polygonize, linemerge
from shapely.geometry import LineString
import numpy as np


def get_squares_from_rect_faster(RectangularPolygon, side_length=0.0025):
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)

    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]

    lines = horizontal_splitters + vertical_splitters
    lines.append(RectangularPolygon.boundary)
    lines = unary_union(lines)
    lines = linemerge(lines)

    return list(polygonize(lines))

Answer 3

I was looking for an answer for this myself but couldn't find one. Tried to stitch together several pieces and here's the result.

This is not necessarily the most optimal routine but it did the job for me. If you want to increase performance, try experimenting with the code.

A brief description of the algorithm:

1. Using the boundaries of the original geometry itself, and the boundaries of its convex hull, and its minimum rotated rectangle, derive all possible rectangles.

2. Divide all rectangles into smaller squares of specified side length.

3. Drop duplicates using a rounded off centroid. (r: round off param)

4. Retain either those squares 'within' the geometry, or those that 'intersect' the geometry, depending on whichever is closer to the total number of required squares.

EDITED

New Solution

#### Python script for dividing any shapely polygon into smaller equal sized polygons

import numpy as np
from shapely.ops import split
import geopandas
from shapely.geometry import MultiPolygon, Polygon

def rhombus(square):
    """
    Naively transform the square into a Rhombus at a 45 degree angle
    """
    coords = square.boundary.coords.xy
    xx = list(coords[0])
    yy = list(coords[1])
    radians = 1
    points = list(zip(xx, yy))
    Rhombus = Polygon(
        [
            points[0],
            points[1],
            points[3],
            ((2 * points[3][0]) - points[2][0], (2 * points[3][1]) - points[2][1]),
            points[4],
        ]
    )
    return Rhombus


def get_squares_from_rect(RectangularPolygon, side_length=0.0025):
    """
    Divide a Rectangle (Shapely Polygon) into squares of equal area.

    `side_length` : required side of square

    """
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = RectangularPolygon
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)

    return square_polygons


def split_polygon(G, side_length=0.025, shape="square", thresh=0.9):
    """
    Using a rectangular envelope around `G`, creates a mesh of squares of required length.
    
    Removes non-intersecting polygons. 
            

    Args:
    
    - `thresh` : Range - [0,1]

        This controls - the number of smaller polygons at the boundaries.
        
        A thresh == 1 will only create (or retain) smaller polygons that are 
        completely enclosed (area of intersection=area of smaller polygon) 
        by the original Geometry - `G`.
        
        A thresh == 0 will create (or retain) smaller polygons that 
        have a non-zero intersection (area of intersection>0) with the
        original geometry - `G` 

    - `side_length` : Range - (0,infinity)
        side_length must be such that the resultant geometries are smaller 
        than the original geometry - `G`, for a useful result.

        side_length should be >0 (non-zero positive)

    - `shape` : {square/rhombus}
        Desired shape of subset geometries. 


    """
    assert side_length>0, "side_length must be a float>0"
    Rectangle    = G.envelope
    squares      = get_squares_from_rect(Rectangle, side_length=side_length)
    SquareGeoDF  = geopandas.GeoDataFrame(squares).rename(columns={0: "geometry"})
    Geoms        = SquareGeoDF[SquareGeoDF.intersects(G)].geometry.values
    if shape == "rhombus":
        Geoms = [rhombus(g) for g in Geoms]
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    elif shape == "square":
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    return geoms

# Reading geometric data

geo_filepath = "/data/geojson/pc_14.geojson"
GeoDF        = geopandas.read_file(geo_filepath)

# Selecting random shapely-geometry

G = np.random.choice(GeoDF.geometry.values)


squares   = split_polygon(G,shape='square',thresh=0.5,side_length=0.025)
rhombuses = split_polygon(G,shape='rhombus',thresh=0.5,side_length=0.025)

Previous Solution:

import numpy as np
import geopandas
from shapely.ops import split
from shapely.geometry import MultiPolygon, Polygon, Point, MultiPoint

def get_rect_from_geom(G, r=2):
    """
    Get rectangles from a geometry.
    r = rounding factor.
    small r ==> more rounding off ==> more rectangles
    """
    coordinate_arrays = G.exterior.coords.xy
    coordinates = list(
        zip(
            [np.round(c, r) for c in coordinate_arrays[0]],
            [np.round(c, r) for c in coordinate_arrays[1]],
        )
    )
    Rectangles = []
    for c1 in coordinates:
        Coords1 = [a for a in coordinates if a != c1]
        for c2 in Coords1:
            Coords2 = [b for b in Coords1 if b != c2]
            x1, y1 = c1[0], c1[1]
            x2, y2 = c2[0], c2[1]
            K1 = [k for k in Coords2 if k == (x1, y2)]
            K2 = [k for k in Coords2 if k == (x2, y1)]
            if (len(K1) > 0) & (len(K2) > 0):
                rect = [list(c1), list(K1[0]), list(c2), list(K2[0])]
                Rectangles.append(rect)
    return Rectangles


def get_squares_from_rect(rect, side_length=0.0025):
    """
    Divide a rectangle into equal area squares
    side_length = required side of square
    """
    y_list = [r[1] for r in rect]
    x_list = [r[0] for r in rect]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells, ycells = int(np.round(width / side_length)), int(
        np.round(height / side_length)
    )

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = Polygon(rect)
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)
    return [np.stack(SQPOLY.exterior.coords.xy, axis=1) for SQPOLY in square_polygons]


def round_centroid(g, r=10):
    """
    Get Centroids.
    Round off centroid coordinates to `r` decimal points.
    """
    C = g.centroid.coords.xy
    return (np.round(C[0][0], r), np.round(C[1][0], r))


def subdivide_polygon(g, side_length=0.0025, r=10):
    """
    1. Create all possible rectangles coordinates from the geometry, its minimum rotated rectangle, and its convex hull.
    2. Divide all rectangles into smaller squares.

    small r ==> more rounding off ==> fewer overlapping squares. (these are dropped as duplicates)
    large r ==> may lead to a few overlapping squares.
    """
    #   Number of squares required.
    num_squares_reqd = g.area // (side_length ** 2)

    # Some of these combinations can be dropped to improve performance.
    Rectangles = []
    Rectangles.extend(get_rect_from_geom(g))
    Rectangles.extend(get_rect_from_geom(g.minimum_rotated_rectangle))
    Rectangles.extend(get_rect_from_geom(g.convex_hull))

    Squares = []
    for r in range(len(Rectangles)):
        rect = Rectangles[r]
        Squares.extend(get_squares_from_rect(rect, side_length=side_length))
    SquarePolygons = [Polygon(square) for square in Squares]
    GDF = geopandas.GeoDataFrame(SquarePolygons).rename(columns={0: "geometry"})
    GDF.loc[:, "centroid"] = GDF.geometry.apply(round_centroid, r=r)
    GDF = GDF.drop_duplicates(subset=["centroid"])

    wgeoms = GDF[GDF.within(g)].geometry.values
    igeoms = GDF[GDF.intersects(g)].geometry.values

    w = abs(num_squares_reqd - len(wgeoms))
    i = abs(num_squares_reqd - len(igeoms))

    print(w, i)
    if w <= i:
        return wgeoms
    else:
        return igeoms



geoms = subdivide(g)

Answer 4

In computational geometry, the problem you want to solve is called triangulation.

There are algorithms to solve this problem, giving triangulations with different properties. You will need to decide which one is the best fit.