Python : shapely, 一个多边形内的级联交叉点

Question

我想将一个多边形拆分为一个多边形列表，该列表对应于与其他多边形的所有交集(以及它们之间的交集)。

from shapely.geometry import Point

circleA = Point((0, 0)).buffer(1)
circleB = Point((1, 0)).buffer(1)
circleC = Point((1, 1)).buffer(1)

def cascaded_intersections(poly1, lst_poly):
    # ???
    return result

result = cascaded_intersections(circleA, (circleB, circleC))

结果应该是 4 个多边形的列表，对应于 A 的 4 个互补部分(上图:[AC!B、ABC、AB!C、A 的其余部分])。

这个问题与将多边形从覆盖线串列表中分割成最小的部分相同。

如何编写cascaded_intersections？

Answer 1

我的一位同事 Pascal L. 找到了一个解决方案:

#!/usr/bin/env python
# -*- coding:utf-8 -*-

from shapely.geometry import MultiPolygon, Polygon, Point, GeometryCollection
from shapely.ops import cascaded_union

EMPTY = GeometryCollection()

def partition(poly_a, poly_b):
    """
    Splits polygons A and B into their differences and intersection.
    """
    if not poly_a.intersects(poly_b):
        return poly_a, poly_b, EMPTY
    only_a = poly_a.difference(poly_b)
    only_b = poly_b.difference(poly_a)
    inter  = poly_a.intersection(poly_b)
    return only_a, only_b, inter


def eliminate_small_areas(poly, small_area):
    """
    Eliminates tiny parts of a MultiPolygon (or Polygon)
    """
    if poly.area < small_area:
        return EMPTY
    if isinstance(poly, Polygon):
        return poly
    assert isinstance(poly, MultiPolygon)
    l = [p for p in poly if p.area > small_area]
    if len(l) == 0:
        return EMPTY
    if len(l) == 1:
        return l[0]
    return MultiPolygon(l)


def cascaded_intersections(poly1, lst_poly):
    """
    Splits Polygon poly1 into intersections of/with list of other polygons.
    """

    result = [(lst_poly[0], (0,))]

    for i, poly in enumerate(lst_poly[1:], start=1):

        current = []

        while result:
            result_geometry, result_indexes = result.pop(0)
            only_result, only_poly, inter = partition(result_geometry, poly)
            for geometry, indexes in ((only_result, result_indexes), (inter, result_indexes + (i,))):
                if not geometry.is_empty:
                    current.append((geometry, indexes))
        current_union = cascaded_union([elt[0] for elt in current])
        only_poly = poly.difference(current_union)
        if not only_poly.is_empty:
            current.append((only_poly, (i,)))
        result = current

    for r in range(len(result)-1, -1, -1):
        geometry, indexes = result[r]
        if poly1.intersects(geometry):
            inter = poly1.intersection(geometry)
            result[r] = (inter, indexes)
        else:
            del result[r]

    only_poly1 = poly1.difference(cascaded_union([elt[0] for elt in result]))
    only_poly1 = eliminate_small_areas(only_poly1, 1e-16*poly1.area)
    if not only_poly1.is_empty:
        result.append((only_poly1, None))

    return [r[0] for r in result]

a=Point(0,0).buffer(1)
b1=Point(0,1).buffer(1)
b2=Point(1,0).buffer(1)
b3=Point(1,1).buffer(1)

result = cascaded_intersections(a, (b1,b2,b3))