基于Shapely多边形生成无重复包含关系树及嵌套运算问题

我有一组互不重叠的Shapely Polygon，需要根据它们的嵌套层级执行union（并集）或difference（差集）运算，生成代表整个集合的最终MultiPolygon。

举个例子：假设有A-D四个嵌套多边形，需要用A减去B和D，再与C做并集——嵌套层级会交替决定是添加还是移除主多边形的部分。

核心需求是生成一棵无重复的多边形包含关系树（比如A包含B，但不包含C）。我当前的实现用了大量循环：

• 遍历所有多边形，检查当前多边形是否被其他多边形包含，将未被包含的加入列表；
• 这些“未被包含”的多边形是当前层级最外层成员（作为上一层级的子节点），逐个从原列表移除后递归搜索；
• 重复上述步骤直到原列表为空。

当前代码的问题 #

我写的示例代码存在递归问题，输出包含5个多边形而非预期的4个：

from shapely.geometry import Polygon
from shapely import contains
import numpy as np


class Branch:
    def __init__(self, root):
        self.root = root
        self.children = []

    def add_child(self, child):
        self.children.append(child)

    def add_children(self, children):
        for child in children:
            self.add_child(child)

    def __repr__(self) -> str:
        return f"{self.root}:{self.children}"


def get_not_contained(polygons):
    # Find which polygons are not contained within any other polygons
    not_contained = []
    for i in range(len(polygons)):
        con = False
        for j in range(len(polygons)):
            # Don't check against itself
            if polygons[i] == polygons[j]:
                continue

            if contains(polygons[j], polygons[i]):
                con = True
                break

        if not con:
            not_contained.append(polygons[i])

    return not_contained

# a box
coords = np.array([[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]])

# several boxes, some of which are nested (nesting labels match the example image)
polygons = [
    Polygon(coords), # D
    Polygon(coords * 5 - 1), # A
    Polygon(coords + 2), # B
    Polygon(coords * 0.5 + 2.25), # C
]

def get_NC_tree(polygons):
    out = []
    NC = get_not_contained(polygons)
    for poly in NC:
        brnch = Branch(poly)
        out.append(brnch)
        tmp = polygons
        try:
            tmp.remove(poly)
        except ValueError:
            pass
        brnch.add_children(get_NC_tree(tmp))

    return out


print(get_NC_tree(polygons))

更简便的实现方案 #

核心优化思路 #

1. 按面积排序：嵌套的外层多边形面积必然大于内层，先按面积从大到小排序，大幅减少包含检查的次数。
2. 精准查找父节点：每个多边形只需检查比它面积大的多边形是否包含它，找到第一个符合条件的作为父节点，避免重复检查。
3. 避免列表引用错误：原代码中tmp = polygons是引用，修改tmp会影响原列表，新方案用字典映射多边形与节点，彻底解决递归时的列表修改问题。
4. 自动计算最终图形：新增递归函数，根据嵌套层级自动执行union或difference，直接输出最终的MultiPolygon。

修正后的代码 #

from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import unary_union
from shapely import contains
import numpy as np

class Branch:
    def __init__(self, root):
        self.root = root
        self.children = []

    def add_child(self, child):
        self.children.append(child)

    def __repr__(self) -> str:
        # 简化输出，用边界框识别多边形
        return f"Polygon({self.root.bounds}):{self.children}"

def build_containment_tree(polygons):
    # 按面积从大到小排序，外层多边形面积更大
    sorted_polys = sorted(polygons, key=lambda p: -p.area)
    tree = []
    # 用字典记录每个多边形对应的Branch节点
    poly_to_branch = {}

    for poly in sorted_polys:
        branch = Branch(poly)
        poly_to_branch[poly] = branch
        # 寻找当前多边形的父节点（第一个包含它的更大多边形）
        parent_found = False
        for candidate in sorted_polys:
            if candidate == poly:
                continue
            if contains(candidate, poly):
                poly_to_branch[candidate].add_child(branch)
                parent_found = True
                break
        # 无父节点则作为根节点加入树
        if not parent_found:
            tree.append(branch)
    return tree

def compute_final_geometry(tree):
    # 递归计算：根节点层级为0，偶数层执行difference，奇数层执行union
    def recurse(branch, level=0):
        if not branch.children:
            return branch.root
        # 先计算所有子节点的运算结果
        child_results = [recurse(child, level+1) for child in branch.children]
        child_union = unary_union(child_results)
        if level % 2 == 0:
            # 偶数层（根）：当前多边形减去子节点的并集
            return branch.root.difference(child_union)
        else:
            # 奇数层：当前多边形与子节点的并集（可根据需求调整逻辑）
            return branch.root.union(child_union)
    
    # 合并所有根节点的运算结果
    root_results = [recurse(branch) for branch in tree]
    return unary_union(root_results)

# 测试用例
coords = np.array([[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]])
polygons = [
    Polygon(coords),                # D（最小，被A包含）
    Polygon(coords * 5 - 1),        # A（最大，根节点）
    Polygon(coords + 2),            # B（被A包含）
    Polygon(coords * 0.5 + 2.25),   # C（被B包含）
]

# 构建包含关系树
containment_tree = build_containment_tree(polygons)
print("包含关系树：")
print(containment_tree)

# 计算最终MultiPolygon
final_poly = compute_final_geometry(containment_tree)
print("\n最终几何图形类型：", type(final_poly))
print("边界框：", final_poly.bounds)

代码说明 #

• 包含关系树构建：按面积排序后，每个多边形只需要和更大的多边形比较包含关系，找到父节点后加入对应分支，无父节点则作为根节点，确保树结构无重复。
• 最终图形计算：递归处理每个节点，根据层级交替执行差集/并集，最后合并所有根节点的结果，直接得到目标MultiPolygon。

内容的提问来源于stack exchange，提问作者Mandias