python 如何从多多边形形状物体中检测出内部多边形？

ckocjqey 于 2022-12-10 发布在 Python

关注(0)|答案(1)|浏览(252)

我想从一个多多边形形状的物体中检测内部多边形。五大湖，黑海和里海应该是内部多边形，并且不被填充。
如何使用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

来源：https://stackoverflow.com/questions/74710358/how-to-detect-inner-polygons-from-a-multipolygon-shapely-object

1条答案

按热度按时间

l7wslrjt1#

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

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()

五大湖、黑海和里海现在都没有填满。
但注意南极画得不对。不明白为什么？