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

Question

我正在寻找一种平滑多边形的方法，使相邻/接触的多边形保持接触。单个多边形可以很容易地进行平滑处理，例如使用PAEK或Bezier插值（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[1:-1] + [(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');

更新： 使用Mountain的建议和这篇文章，我认为问题可以分解成以下步骤：

• 找到每对相邻多边形之间的边界边缘（例如，使用这个建议）。
• 将它们转换为numpy数组并根据Mountain的bspline建议进行平滑处理。
• 使用更新/平滑的边缘重构多边形。

还要注意，对于单个（shapely.geometry）多边形，可以使用pol.simplify()使用Douglas-Peucker算法进行平滑处理。

- Thijs

---

我不熟悉这些库，但如果是2D几何图形，你不能使用numpy来生成几何图形和平滑吗？ - DGKang

我所知道的唯一能够实现这一点的库是JavaScript的topojson。如果你在Python中找到了什么东西，我很乐意听听！话虽如此，不幸的是，在这里要求库推荐并不是主题。 - Michael Delgado

是的，完全正确。很遗憾，我不知道这些库中有任何开箱即用的方法来完成这个任务。 - Michael Delgado

对于问题中的多边形情况，也许可以取多边形中相互接触的中间点，然后进行线性回归以获得二次方程来绘制直线？ - DGKang

是的，没错 - 但我正在寻找适用于N个不仅仅是直线的多边形的通用解决方案。只是用这两个多边形作为简单的例子，抱歉没有表述清楚。 - Thijs

Answer 1

在查看了一段时间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://dev59.com/kV4c5IYBdhLWcg3wJnQ-
    """ 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()

- Mountain

---

嗨，这是一个很好的建议，非常感谢！我正在使用原始帖子中概述的策略来实现它，等我有更新时会在这里发布。 - Thijs

恐怕您可能会很难让相邻的多边形与此匹配得很好，但期待看到结果如何！ - AKX

Answer 2

在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

- Kevin