How to detect inner polygons from a multipolygon shapely object?

Question

I would like to detect inner polygons from a multipolygon shapely object. Great lakes, Black Sea and Caspian sea should be inner polygons and not be filled.

How to do this properly with shapefile ?

Please find the script bellow for investigating.

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

PBrockmann · Answer 1 · 2023-04-21T12:49:43.333

So the strategy is to separate the different self intersected polygons from their duplicate points into multiple polygons. Then build a path with eventually holes, with option closed=True. Then a patch object. Then no overlay, the filling of each polygons is correct whatever the projection.

Note that for the path object, exterior ring will be oriented counter-clockwise and the interior rings (holes) will be oriented clockwise.

The first idea was to orient correctly initial self intersected polygons but shapely.geometry.polygon.orient does not give a correct orientation with such objects. You have first to separate them into individual polygons.

See: See: Separate vtk self intersected polydata into multiple polygons from duplicate points? and https://discourse.vtk.org/t/create-separated-regions-from-polydata

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

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

with open('./polys_03.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)
    
#-----------------------------------------
def polygon2path(poly):
    path = Path.make_compound_path(
               Path(poly.exterior.coords, closed=True),
               *[Path(ring.coords) for ring in poly.interiors])
    return path
#-----------------------------------------
fig = plt.figure(figsize=(10,5))

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

transform = ccrs.Geodetic()

paths = []
holesNumber = []
for n,polygonA in enumerate(polygons):
    holes = []
    for m,polygonB in enumerate(polygons):
        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 = shapely.geometry.Polygon(polygonA.exterior.coords, holes=holes)
    polygonNew = shapely.geometry.polygon.orient(polygonNew)   # Orient to be oriented counter-clockwise
    
    path = polygon2path(polygonNew)
    paths.append(path)

    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = PathPatch(path, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)

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

How to detect inner polygons from a multipolygon shapely object?

1 Answers1