在具有多边形形状的MultiPolygon中填充孔-荷兰2位邮政编码[重复]

Question

此问题已经在这里有了答案：

Convert Multipolygon to Polygon in Python [closed]（2个答案）

关闭2分钟前。

我在荷兰找到了一个用于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))