如何从多边形形状对象中检测内部多边形？答案

【问题标题】：How to detect inner polygons from a multipolygon shapely object?如何从多边形形状对象中检测内部多边形？
【发布时间】：2022-12-09 18:25:05
【问题描述】：

我想从多多边形形状对象中检测内部多边形。大湖、黑海和里海应该是内部多边形而不是被填充。

如何使用 shapefile 正确执行此操作？

请找到下面的脚本进行调查。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely import geometry
import random
import pickle

! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle

with open('./polys.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson(10))

transform = ccrs.Geodetic()

for polygon in polygons.geoms:
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    x = polygon.exterior.coords.xy[0]
    y = polygon.exterior.coords.xy[1]
    ax.fill(x, y, transform=transform, color=random_color, lw=0.5, edgecolor="black")

ax.set_global()
ax.gridlines()
plt.show()

【问题讨论】：

标签： python shapely

【解决方案1】：

这是一个使用补丁的解决方案，因为您可以处理带孔的多边形。请注意，外环将逆时针定向，内环（孔）将顺时针定向。

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle

with open('./polys.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    path = Path.make_compound_path(
           Path(poly.exterior.coords),
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(-10)
#map_proj = ccrs.Orthographic(-10, -60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
    holes = []
    for m,polygonB in enumerate(polygons.geoms):
        if (m == n): continue     
        if polygonA.contains(polygonB):
            holes.append(polygonB.exterior.coords)
            holesNumber.append(m)
    if n in holesNumber: continue  # n is a hole
    polygonNew = geometry.Polygon(polygonA.exterior.coords, holes=holes)
    polygonNew = shapely.geometry.polygon.orient(polygonNew)   # Orient to be oriented counter-clockwise
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

大湖、黑海和里海现在还没有被填满。

但请注意，南极洲的绘制不正确。不明白为什么？

更明显的是正射投影填充不正确的问题。看到紫罗兰直填满非洲。在这一点上也欢迎任何帮助吗？

已经发布了一个关于这个的问题：https://github.com/SciTools/cartopy/issues/2111

【讨论】：