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