How to smooth adjacent polygons in Python?

Question

I'm looking for a way to smooth polygons such that adjacent/touching polygons remain touching. Individual polygons can be smoothed easily, e.g., with PAEK or Bezier interpolation (https://pro.arcgis.com/en/pro-app/latest/tool-reference/cartography/smooth-polygon.htm), which naturally changes their boundary edge. But how to smooth all polygons such that touching polygons remain that way?

I'm looking for a Python solution ideally, so it can easily be automated. I found an equivalent question for Arcgis (https://gis.stackexchange.com/questions/183718/how-to-smooth-adjacent-polygons), where the top answer outlines a good strategy (converting polygon edges to lines from polygon-junction to junction), smoothing these and then reconstructing the polygons). Perhaps this would the best strategy, but I'm not sure how to convert shared polygon boundaries to individual polylines in Python.

Here is some example code that shows what I'm trying to do for just 2 polygons (but I've created the 'smoothed' polygons by hand):

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');

Update: Using the suggestion from Mountain below and this post, I think the problem could be broken down in the following steps:

• Find boundary edges between each pair of touching polygons (e.g., using this suggestion).
• Transform these into numpy arrays and smooth as per Mountain's bspline suggestion
• Reconstruct polygons using updated/smoothed edges.

Also note that for single (shapely.geometry) polygons, they can be smoothed using: pol.simplify() using Douglas-Peucker algorithm.

Answer 1

After looking through scipy's interpolation functions for a while, I found a 2D interpolation function created by Fnord. You can interpolate ND arrays.

Fnord's Genius: Splines with Python (using control knots and endpoints)

Below is an example of interpolating the points of a circle for a smoother surface. Tried creating a circular array to handle edge effects. The indexing after up-sampling is a bit tricky.

Hopefully this gives a starting point to your quest.

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()

Answer 2

Within the Cartographic toolset in ArcGIS (Pro is my preferred), there are smooth and simplify algorithms. These respect shared feature boundaries. If you are getting unexpected results, I would run Integrate (https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/integrate.htm) on the input dataset before running Smooth/Simplify Polygon. I use these tools in production flood mapping. I keep seeing other resources which assert that these tools do not respect shared boundaries. This is simply untrue.

My in-doc tool runs RepairGeometry->Integrate->Smooth Polygon->Simplify Polygon. For my use-case, it makes the most sense to use PAEK (smooth) and Retain Critical Weighted-Areas (Zhou-Jones). https://pro.arcgis.com/en/pro-app/latest/tool-reference/cartography/smooth-polygon.htm