Splitting polygons using triangulation

triangulation.py

from collections import deque
from shapely.geometry import Polygon, Point, LinearRing
from shapely import force_2d, has_z, has_m, constrained_delaunay_triangles
from shapely.validation import make_valid
from typing import List, Generator

def triangulate_polygon(polygon: Polygon) -> List[Polygon]:
    """Triangulate a given polygon."""
    return [tri for tri in constrained_delaunay_triangles(polygon).geoms if not tri.is_empty and tri.is_valid]

def split_triangle(triangle: Polygon, area_threshold: float) -> Generator[Polygon, None, None]:
    """Recursively split a triangle into two smaller triangles if its area exceeds the threshold."""
    queue = deque([triangle])

    while queue:
        current_triangle = queue.popleft()  # deque provides faster popping from the left
        if current_triangle.area > area_threshold:
            # Retrieve the coordinates of the triangle's vertices
            coords = list(current_triangle.exterior.coords)[:3]
            
            # Determine the longest edge using squared length for efficiency
            longest_edge = max(
                ((coords[i], coords[(i + 1) % 3]) for i in range(3)),
                key=lambda pair: (pair[1][0] - pair[0][0]) ** 2 + (pair[1][1] - pair[0][1]) ** 2
            )

            # Midpoint of the longest edge
            midpoint = (
                (longest_edge[0][0] + longest_edge[1][0]) / 2,
                (longest_edge[0][1] + longest_edge[1][1]) / 2
            )

            # Find the remaining vertex not on the longest edge
            remaining_vertex = next(pt for pt in coords if pt not in longest_edge)

            # Create two new triangles using the remaining vertex and the midpoint
            new_triangles = [
                Polygon([longest_edge[0], midpoint, remaining_vertex]),
                Polygon([longest_edge[1], midpoint, remaining_vertex]),
            ]

            for tri in new_triangles:
                if tri.area > area_threshold:
                    queue.append(tri)  # Add to queue for further splitting
                else:
                   yield tri
        else:
            yield current_triangle

def subdivide_polygon(geometry: Polygon, threshold: float, check_2D: bool = True) -> Generator[Polygon, None, None]:
    """
    Subdivide the input polygon into smaller triangles based on an area threshold.
    Invalid input `geometry` will silently be made valid (if possible).
    Any LinearRings will be converted to Polygons.
    """
    if geometry is None:
        return []
    if isinstance(geometry, LinearRing):
        geometry = Polygon(geometry)
    if check_2D and (has_z(geometry) or has_m(geometry)):
        geometry = force_2d(geometry)
    if not geometry.is_valid:
        geometry = make_valid(geometry)
        if geometry.geom_type == "GeometryCollection":
            geometry.normalize()
        geometry = geometry.buffer(0)
    
    if geometry.area < threshold:
        yield geometry
        return
    
    initial_triangles = triangulate_polygon(geometry)

    for triangle in initial_triangles:
        if not triangle.is_valid or triangle.is_empty:
            continue
        yield from split_triangle(triangle, threshold)

if __name__ == '__main__':
    from shapely import wkt
    from shapely.geometry import GeometryCollection

    polygon_a = wkt.loads(
        'POLYGON ((-62.490234 37.09024, -87.363281 30.524413, -86.484375 21.616579, '
        '-81.5625 15.029686, -64.248047 16.972741, -64.599609 25.562265, '
        '-70.576172 25.641526, -59.941406 31.052934, -52.470703 35.317366, '
        '-62.490234 37.09024))'
    )
    polygon_b = wkt.loads('POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))')
    polygon_c = wkt.loads('MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),(30 20, 20 15, 20 25, 30 20)))')
    polygon_d = wkt.loads('POLYGON ((-42.978516 31.952162, -65.478516 29.382175, -47.460937 26.509905, -72.509766 19.311143, -48.779297 20.632784, -47.109375 11.436955, -44.033203 21.207459, -34.980469 16.467695, -40.869141 23.483401, -27.333984 26.037042, -35.068359 28.304381, -43.242188 27.605671, -42.978516 31.952162))')

    area_threshold = 0.5

    def timing(f):
        @wraps(f)
        def wrap(*args, **kw):
            ts = time()
            result = f(*args, **kw)
            te = time()
            print('func:%r args:[%r, %r] took: %2.4f sec' % \
            (f.__name__, args, kw, te-ts))
            return result
        return wrap

    @timing
    def timed_test_execution(*args):
        for geom in [polygon_a, polygon_b, polygon_c, polygon_d]:
            collection = list(subdivide_polygon(geom, *args))
            # print(GeometryCollection(collection), '\n')

    timed_test_execution(area_threshold)