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

You can do this with the topojson library.

Sample script:

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

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:
topo = topojson.Topology(df_pols_coarse)
topo_smooth = topo.toposimplify(1)
df_pols_smooth = topo_smooth.to_gdf()

## 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")
plt.show()

Answer 2

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