Use shapely.Polygon.simplify() to reduce ice wedge polygon vertices while maintaining geometry fidelity

simplify-polygons.py

import fiona
from pprint import pprint
from shapely.geometry import Point, Polygon
from descartes import PolygonPatch
import matplotlib.pyplot as plt
from pathlib import Path

data_folder = os.path.join(
    Path.home(),
    "Documents",
    "Development",
    "csjx",
    "permafrost-ice-wedge-layer",
    "resources", 
    "WV02_20120911234552_103001001B2C9200_12SEP11234552-M1BS-052903555090_01_P003_u16rf3413_pansh")

path = os.path.join(data_folder, "WV02_20120911234552_103001001B2C9200_12SEP11234552-M1BS-052903555090_01_P003_u16rf3413_pansh.shp")
with fiona.open(path, "r") as source:
    # The lon/lat bounds stated in the original tif file

    count = 0
    tolerance = 0.1 # The simplification tolerance level
    for feature in source:

        # Create a figure
        fig = plt.figure(1, figsize=[20, 10], dpi=72)
        fig.set_frameon(True)
        # Add the left side subplot
        ax = fig.add_subplot(121)
        # Plot the current ice wedge polygon
        polygon = Polygon(feature['geometry']['coordinates'][0])
        xmin, ymin, xmax, ymax = polygon.bounds
        ax.set_xlim(xmin - 1, xmax+ 1)
        ax.set_ylim(ymin - 1, ymax + 1)
        patch = PolygonPatch(polygon, facecolor="skyblue", edgecolor="blue", alpha=0.5, zorder=2)
        ax.add_patch(patch)

        # Plot the exterior ring vertices
        for p in polygon.exterior.coords:
            point = Point(p)
            x, y = point.xy
            ax.plot(x, y, markersize="2", marker="o", color="red")

        # Simplify the polygon to reduce vertices based on the tolerance level
        polygon2 = polygon.simplify(tolerance=tolerance, preserve_topology=True)
        # And add the right-side subplot to the figure
        ax = fig.add_subplot(122)
        ax.set_xlim(xmin - 1, xmax+ 1)
        ax.set_ylim(ymin - 1, ymax + 1)
        patch = PolygonPatch(polygon2, facecolor="skyblue", edgecolor="blue", alpha=0.5, zorder=2)
        ax.add_patch(patch)

        # Plot the simplified exterior ring vertices
        for p in polygon2.exterior.coords:
            point = Point(p)
            x, y = point.xy
            ax.plot(x, y, markersize="2", marker="o", color="red")

        # Display the results
        plt.show()
        poly_length = len(polygon.exterior.coords)
        poly2_length = len(polygon2.exterior.coords)
        percent = 100 - ((poly2_length/poly_length)*100)
        print(f"polygon2 = polygon.simplify(tolerance={tolerance:0.1f}, preserve_topology=True)")
        print(f"Poly 1 vertices: {str(poly_length)}")
        print(f"Poly 2 vertices: {str(len(polygon2.exterior.coords))}")
        print(f"{percent:.2f}% decrease in size\n")
        
        count += 1

        # Show deformed polygons at higher tolerance levels
        if count == 2:
            tolerance = 0.3
            continue

        if count == 3:
            tolerance = 0.5
            continue

        if count == 4:
            tolerance = 0.7
            continue

        # Limit the number of figures to 20 polygons
        if count > 20:
            break

        # After 4 plots, plot the rest at the minimum tolerance again
        if count > 4:
            tolerance = 0.1
            continue