如何用Python双层循环填充多边形邻接关系二维数组

我是Python新手，目前遇到技术问题：我有一个名为triangles的列表，包含24个Shapely多边形，部分多边形彼此相邻。我需要构建一个二维数组（矩阵），存储每个多边形对应的邻接多边形索引（而非仅布尔值），例如多边形0邻接1、2、10，矩阵第一行对应位置需填入这些索引。我已有双层循环的大致框架，但需要具体实现方法。生成这些三角形的代码如下：

import matplotlib.pyplot as plt
from shapely.geometry import Point, LineString
from shapely.ops import unary_union, polygonize
from matplotlib.pyplot import cm
import numpy as np

def plot_coords(coords, color):
    pts = list(coords)
    x, y = zip(*pts)
    # print(color)
    plt.plot(x,y, color='k', linewidth=1)
    plt.fill_between(x, y, facecolor=color)


def plot_polys(polys, color):
    for poly, color in zip(polys, color):
        plot_coords(poly.exterior.coords, color)

x = 0
y = 0
h = 1.73205080757

points = [# center
          Point(x, y),
          #  first ring
          Point((x + 2), y),
          Point((x - 2), y),
          Point((x + 1), (y + h)),
          Point((x - 1), (y + h)),
          Point((x + 1), (y - h)),
          Point((x - 1), (y - h)),
          # second ring
          Point((x + 3), h),
          Point((x - 3), h),
          Point((x + 3), -h),
          Point((x - 3), -h),
          Point((x + 2), (h + h)),
          Point((x - 2), (h + h)),
          Point((x + 2), (-h + -h)),
          Point((x - 2), (-h + -h)),
          Point((x + 4), y),
          Point((x - 4), y),
          Point(x, (h + h)),
          Point(x, (-h + -h)),
          #third ring
          Point((x + 4), (h + h)),
          Point((x - 4), (h + h)),
          Point((x + 4), (-h + -h)),
          Point((x - 4), (-h + -h)),
          Point((x + 1), (h + h + h)),
          Point((x - 1), (h + h + h)),
          Point((x + 1), (-h + -h + -h)),
          Point((x - 1), (-h + -h + -h)),
          Point((x + 5), h),
          Point((x - 5), h),
          Point((x + 5), -h),
          Point((x - 5), -h)]

# buffer points to create circle polygons

circles = []
for point in points:
    circles.append(point.buffer(2))


# unary_union and polygonize to find overlaps

rings = [LineString(list(pol.exterior.coords)) for pol in circles]
union = unary_union(rings)
result_polys = [geom for geom in polygonize(union)]

# remove tiny sliver polygons
threshold = 0.01
filtered_polys = [polygon for polygon in result_polys if polygon.area > threshold]

# remove outer circle fragments
complete_polys = [polygon for polygon in filtered_polys if (polygon.centroid.x**2 + polygon.centroid.y**2 < 4**2)]

print("total polygons = " + str(len(result_polys)))
print("filtered polygons = " + str(len(filtered_polys)))
print("complete polygons = " + str(len(complete_polys)))

centroidCoords = [
    Point(0.9999999999999997, -0.5773635180985456),
    Point(1.5000000000000004, -0.866025403785),
    Point(1.9999999999999987, -1.154687289471454),
    Point(0, -1.1546872894714544),
    Point(0, -1.7320508075699992),
    Point(-0.0000000000000002, -2.3094143256685467),
    Point(-1.0000000000000004, -0.577363518098546),
    Point(-1.5, -0.8660254037850001),
    Point(-1.9999999999999996, -1.154687289471454),
    Point(-1, 0.5773635180985459),
    Point(-1.5000000000000004, 0.8660254037850001),
    Point(-1.9999999999999993, 1.1546872894714537),
    Point(-0.0000000000000001, 1.1546872894714542),
    Point(0, 1.73205080757),
    Point(0.0000000000000001, 2.309414325668546),
    Point(1, 0.5773635180985456),
    Point(1.5000000000000002, 0.8660254037849999),
    Point(1.9999999999999993, 1.1546872894714542),
    Point(2.9999999999999987, -0.5773635180985457),
    Point(3.5, -0.8660254037849997),
    Point(2, -1.73205080757),
    Point(1.9999999999999996, -2.309414325668546),
    Point(0.5, -0.8660254037850001),
    Point(0.5, 0.866025403785),
    Point(2, 1.7320508075700005),
    Point(2.0000000000000004, 2.3094143256685467),
    Point(3.0000000000000004, 0.5773635180985458),
    Point(3.4999999999999996, 0.8660254037849999),
    Point(-0.5000000000000001, -0.866025403785),
    Point(-2.0000000000000004, -1.7320508075699999),
    Point(-2.000000000000001, -2.3094143256685475),
    Point(-3.000000000000001, -0.5773635180985459),
    Point(-3.500000000000001, -0.866025403785),
    Point(-3.0000000000000018, 0.5773635180985461),
    Point(-3.4999999999999987, 0.8660254037849999),
    Point(-2.0000000000000004, 1.73205080757),
    Point(-1.9999999999999996, 2.3094143256685467),
    Point(-0.5, 0.8660254037849999),
    Point(2.5, 0.866025403785),
    Point(1, 0),
    Point(-0.5000000000000001, 2.5980762113550004),
    Point(-0.9999999999999994, 2.886738097041452),
    Point(1.0000000000000009, 2.8867380970414565),
    Point(0.9999999999999997, 3.4641016151399993),
    Point(2.5000000000000004, 2.598076211355),
    Point(-0.9999999999999997, 0),
    Point(-2.5, 0.8660254037850001),
    Point(-2.500000000000001, 2.598076211355001),
    Point(-0.9999999999999999, 3.4641016151399993),
    Point(0.4999999999999998, 2.5980762113549987),
    Point(2.5000000000000004, -2.598076211355),
    Point(1.0000000000000002, -2.8867380970414547),
    Point(1.0000000000000002, -3.4641016151400006),
    Point(-0.5, -2.598076211355001),
    Point(-0.9999999999999997, -2.8867380970414533),
    Point(2.5000000000000004, -0.8660254037850001),
    Point(0.5, -2.598076211355),
    Point(-0.9999999999999992, -3.464101615139999),
    Point(-2.5, -2.598076211355001),
    Point(-2.5000000000000004, -0.8660254037850001),
    Point(3, 0),
    Point(1.5000000000000002, 2.5980762113550004),
    Point(-3.0000000000000004, 0),
    Point(-1.4999999999999998, 2.5980762113549987),
    Point(1.4999999999999998, -2.598076211355001),
    Point(-1.5, -2.598076211355)]


# separating into petals and triangles

limit = 0.66

triangles = [polygon for polygon in complete_polys if polygon.area < limit]
petals = [polygon for polygon in complete_polys if polygon.area > limit]


for i in range(len(triangles)):
    plot_coords(triangles[i].exterior.coords, triangles_colors[i])
    centroidCoords = triangles[i].centroid
    plt.scatter(centroidCoords.x,centroidCoords.y, marker="$"+str(i)+"$")
    
print(len(triangles))
ax.set_aspect('equal')
plt.show()

邻接矩阵实现方案 #

直接基于你已有的双层循环框架，结合Shapely的几何判断方法即可完成：

核心代码

# 初始化邻接矩阵，每个元素是对应多边形的邻接索引列表
adjacency_matrix = []

n = len(triangles)
for i in range(n):
    neighbors = []
    for j in range(n):
        # 跳过自身，判断两个多边形是否邻接（共享边界且不重叠）
        if i != j and triangles[i].touches(triangles[j]):
            neighbors.append(j)
    adjacency_matrix.append(neighbors)

# 输出示例：查看每个多边形的邻接索引
for idx, neighbor_list in enumerate(adjacency_matrix):
    print(f"多边形{idx}的邻接索引：{neighbor_list}")

精度优化说明

如果touches()方法因为浮点数精度问题出现误判（比如微小间隙或重叠导致邻接判断错误），可以改用交集长度判断：

# 替代touches()的判断逻辑，通过交集线段长度确认真正邻接
if i != j:
    intersection = triangles[i].intersection(triangles[j])
    # 若交集是线段且长度大于极小值，视为邻接
    if intersection.geom_type == "LineString" and intersection.length > 1e-6:
        neighbors.append(j)

复杂度说明

对于24个多边形，O(n²)的双层循环完全足够，无需额外优化。如果后续多边形数量大幅增加，可以考虑提前提取每个多边形的边界线段，减少重复计算。

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