如何在 Python 中平滑相邻的多边形？

Question

我正在寻找一种方法来平滑多边形，使相邻/接触的多边形保持接触。可以轻松平滑单个多边形，例如使用 PAEK 或贝塞尔插值 (https://pro.arcgis.com/en/pro-app/latest/tool-reference/cartography/smooth-polygon.htm)，这自然改变他们的边界边缘。但是如何平滑所有多边形以使接触的多边形保持原样？

我正在理想地寻找 Python 解决方案，因此它可以轻松实现自动化。我发现了一个 Arcgis 的等效问题 (https://gis.stackexchange.com/questions/183718/how-to-smooth-adjacent-polygons)，其中最佳答案概述了一个很好的策略（将多边形边转换为线多边形连接到连接），平滑这些然后重建多边形）。也许这是最好的策略，但我不确定如何在 Python 中将共享的多边形边界转换为单独的多段线。

这里是一些示例代码，显示了我试图对 2 个多边形执行的操作（但我已经手动创建了“平滑”多边形）：

import matplotlib.pyplot as plt
import geopandas as gpd
from shapely import geometry

x_min, x_max, y_min, y_max = 0, 20, 0, 20

## Create original (coarse) polygons:
staircase_points = [[(ii, ii), (ii, ii + 1)] for ii in range(x_max)]
staircase_points_flat = [coord for double_coord in staircase_points for coord in double_coord] + [(x_max, y_max)]

list_points = {1: staircase_points_flat + [(x_max, y_min)],
               2: staircase_points_flat + [(x_min, y_max)]}
pols_coarse = {}
for ind_pol in [1, 2]:
    list_points[ind_pol] = [geometry.Point(x) for x in list_points[ind_pol]]
    pols_coarse[ind_pol] = geometry.Polygon(list_points[ind_pol])

df_pols_coarse = gpd.GeoDataFrame({'geometry': pols_coarse.values(), 'id': pols_coarse.keys()})

## Create smooth polygons (manually):
pols_smooth = {1: geometry.Polygon([geometry.Point(x) for x in [(x_min, y_min), (x_max, y_min), (x_max, y_max)]]),
               2: geometry.Polygon([geometry.Point(x) for x in [(x_min, y_min), (x_min, y_max), (x_max, y_max)]])}
df_pols_smooth = gpd.GeoDataFrame({'geometry': pols_smooth.values(), 'id': pols_smooth.keys()})

## Plot
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
df_pols_coarse.plot(column='id', ax=ax[0])
df_pols_smooth.plot(column='id', ax=ax[1])
ax[0].set_title('Original polygons')
ax[1].set_title('Smoothed polygons');

Answer 1

在 ArcGIS 的制图工具集中（Pro 是我的首选），有平滑和简化的算法。这些尊重共享特征边界。如果您得到意想不到的结果，我会在运行之前对输入数据集运行 Integrate (https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/integrate.htm)平滑/简化多边形。 我在生产洪水映射中使用这些工具。 我不断看到其他资源断言这些工具不尊重共享边界。这完全是不真实的。

我的文档内工具运行 RepairGeometry->Integrate->Smooth Polygon->Simplify Polygon。对于我的用例，使用 PAEK（平滑）和保留关键加权区域 (Zhou-Jones) 最有意义。 https://pro.arcgis.com/en/pro-app/latest/tool-reference/cartography/smooth-polygon.htm

Answer 2

在浏览了一段时间的 scipy 插值函数后，我发现了 Fnord 创建的 2D 插值函数。您可以插入 ND 数组。

Fnord 的天才：Python 样条曲线（使用控制节点和端点）

下面是一个插值圆的点以获得更平滑表面的示例。尝试创建一个圆形阵列来处理边缘效应。上采样后的索引有点棘手。

希望这能为您的探索提供一个起点。

import numpy as np
import geopandas as gpd
from shapely import Polygon
from matplotlib import pyplot as plt
import scipy.interpolate as si

def bspline(cv, n=100, degree=3):
    # this guy is a boss
    # https://stackoverflow.com/questions/28279060/splines-with-python-using-control-knots-and-endpoints
    """ Calculate n samples on a bspline

        cv :      Array ov control vertices
        n  :      Number of samples to return
        degree:   Curve degree
    """
    cv = np.asarray(cv)
    count = cv.shape[0]

    # Prevent degree from exceeding count-1, otherwise splev will crash
    degree = np.clip(degree,1,count-1)

    # Calculate knot vector
    kv = np.array([0]*degree + list(range(count-degree+1)) + [count-degree]*degree,dtype='int')

    # Calculate query range
    u = np.linspace(0,(count-degree),n)

    # Calculate result
    return np.array(si.splev(u, (kv,cv.T,degree))).T

angle = np.linspace(0, 2*np.pi, 10, endpoint=False)
points = np.vstack((np.cos(angle), np.sin(angle))).T

# upsample
overlap = round(len(points)/2)
upsample = 5
extended_points = np.vstack((points[-overlap:], points, points[:overlap]))
new_points = bspline(extended_points, n=len(extended_points)*upsample)[(overlap-2)*upsample:(-overlap+1)*upsample]

p1 = gpd.GeoSeries(Polygon(points))
p1.plot()

p2 = gpd.GeoSeries(Polygon(new_points))
p2.plot()

plt.show()