用匀称的方式填充 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

我找到了这篇文章：

Convert Multipolygon to Polygon in 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))