用 shapely 填充 MultiPolygon 中的洞 - 荷兰 2 位数邮政编码

Question

我在这里找到了荷兰 4 位邮政编码的 shapefile：

https://public.opendatasoft.com/explore/dataset/openpostcodevlakkenpc4/export/?location=9,52.16172,5.56595&basemap=jawg.streets

我想做的是将共享前两位数字的邮政编码组合起来，然后将它们绘制在地图上。

不幸的是，数据似乎有问题——荷兰的一些地区似乎没有被覆盖。结果，当我组合多边形以获取 2 位数多边形（从 4 位数多边形）时，生成的多边形中存在孔。

我需要填补这些漏洞。我看过类似问题的帖子，但似乎没有什么能完全满足我的需要。特别是，我看到一篇关于使用凹壳的帖子，但这在这里似乎有点矫枉过正。

我已经设法修复了几个工件（例如“条子”），但漏洞仍然存在。

这是我到目前为止所拥有的：

from shapely.ops import cascaded_union
from shapely.geometry import JOIN_STYLE, Polygon, MultiPolygon

import os
import folium
import pandas as pd
import geopandas as gpd

GEOM_LOC = r"PATH_TO_FILE_ABOVE\openpostcodevlakkenpc4.shx"

# Get the data
geom = gpd.read_file(GEOM_LOC)

# Remove empty or nan records
is_empty = geom.geometry.is_empty
is_nan = geom.geometry.isna()
geom = geom[~(is_empty | is_nan)]

# Add field for 2 digit postcode
geom["digit"] = geom.pc4.apply(lambda x: x[0:2])
geom = geom[["digit", "geometry"]]

# Make dataframe for 2-digit zones
tags = list(set(geom["digit"]))
df = pd.DataFrame(tags, columns=["tag"])

# Function returning union of geometries
def combine_borders(*geoms):
    return cascaded_union([
        geom if geom.is_valid else geom.buffer(0) for geom in geoms
    ])

# Wrapping function above for application to dataframe
def combine(tag):
    sub_geom = geom[geom["digit"] == tag]
    bounds = list(sub_geom.geometry)
    return(combine_borders(*bounds))

# Apply the function
df["boundary"] = df.tag.apply(lambda x: combine(x))

# The polygons generated above do not fit perfectly together,
# resulting in artifacts. We fix that now
eps = 0.00001

def my_buffer(area):
    return(
        area.buffer(eps, 1, join_style=JOIN_STYLE.mitre).buffer(-eps, 1, join_style=JOIN_STYLE.mitre)
    )

df["boundary"] = df.boundary.apply(lambda x: my_buffer(x))

# Have a look at one of the holes
test = geom[geom.digit == "53"]

test = df.loc[df.tag == "53", "boundary"].item()
m = folium.Map(location=[51.8, 5.4], zoom_start=10, tiles="CartoDB positron")
m.choropleth(test)
m

如您所见，我已经在两位数区域“53”上进行了测试。有一个坑需要填补：

我意识到底层几何结构相当复杂，但是有没有一种简单的方法来填充这个“洞”？

非常感谢您的帮助。

编辑 - 2020-04-25-16.14

作为参考，这里是来自维基百科的 2 位邮政编码区域“53”：

所以这并不是说邮政编码内有邮政编码——邮政编码飞地。

编辑 - 2020-04-27

我找到了这个帖子：

在 Python 中将多边形转换为多边形

应用代码

no_holes = MultiPolygon(Polygon(p.exterior) for p in m)

我的多边形封闭了这个洞：

Answer 1

这封闭了 GeoDataFrame 的多边形几何形状中的简单漏洞。

def close_holes(poly: Polygon) -> Polygon:
        """
        Close polygon holes by limitation to the exterior ring.
        Args:
            poly: Input shapely Polygon
        Example:
            df.geometry.apply(lambda p: close_holes(p))
        """
        if poly.interiors:
            return Polygon(list(poly.exterior.coords))
        else:
            return poly

df = df.geometry.apply(lambda p: close_holes(p))