pythongispolygongeopandasshapely

How to smooth adjacent polygons in Python?

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

Expected outcome smoothed touching polygons

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

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

Solution

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