Add polygonize simplification user guide notebook (#1151)

Author: brendancol

examples/user_guide/50_Polygonize_Simplification.ipynb

@@ -0,0 +1,213 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Polygonize with Geometry Simplification\n",
+    "\n",
+    "`polygonize()` converts raster regions into vector polygons. On high-resolution rasters the result can have thousands of vertices per polygon, which slows rendering and inflates file size.\n",
+    "\n",
+    "The `simplify_tolerance` parameter simplifies polygon boundaries during polygonization. Two algorithms are available:\n",
+    "- `simplify_method=\"douglas-peucker\"` (default): distance-based, removes vertices that deviate less than the tolerance from the simplified line\n",
+    "- `simplify_method=\"visvalingam-whyatt\"`: area-based, removes vertices that form triangles smaller than the tolerance threshold\n",
+    "\n",
+    "Adjacent polygons share identical simplified boundaries, so no gaps or overlaps are introduced."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib.patches import Polygon as MplPolygon\n",
+    "from matplotlib.collections import PatchCollection\n",
+    "\n",
+    "from xrspatial import polygonize"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate a classified raster\n",
+    "\n",
+    "A synthetic land-cover raster with irregular region boundaries."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rng = np.random.default_rng(42)\n",
+    "shape = (80, 120)\n",
+    "\n",
+    "from scipy.ndimage import gaussian_filter\n",
+    "noise = rng.standard_normal(shape)\n",
+    "smooth = gaussian_filter(noise, sigma=8)\n",
+    "classified = np.digitize(smooth, bins=[-0.5, 0.0, 0.5]) + 1\n",
+    "\n",
+    "raster = xr.DataArray(classified.astype(np.int32))\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "raster.plot(ax=ax, cmap=\"Set2\", add_colorbar=True)\n",
+    "ax.set_title(\"Classified raster (4 land-cover types)\")\n",
+    "ax.set_aspect(\"equal\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Polygonize without simplification"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "col_raw, pp_raw = polygonize(raster)\n",
+    "total_verts_raw = sum(len(r) for rings in pp_raw for r in rings)\n",
+    "print(f\"Polygons: {len(pp_raw)}, Total vertices: {total_verts_raw}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_polygons(polygon_points, column, title, ax):\n",
+    "    cmap = plt.cm.Set2\n",
+    "    vals = sorted(set(column))\n",
+    "    val_to_color = {v: cmap(i / max(len(vals) - 1, 1)) for i, v in enumerate(vals)}\n",
+    "\n",
+    "    patches = []\n",
+    "    colors = []\n",
+    "    for val, rings in zip(column, polygon_points):\n",
+    "        ext = rings[0]\n",
+    "        patches.append(MplPolygon(ext[:, :2], closed=True))\n",
+    "        colors.append(val_to_color[val])\n",
+    "\n",
+    "    pc = PatchCollection(patches, facecolors=colors, edgecolors=\"black\",\n",
+    "                         linewidths=0.3)\n",
+    "    ax.add_collection(pc)\n",
+    "    ax.set_xlim(0, 120)\n",
+    "    ax.set_ylim(0, 80)\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_title(title)\n",
+    "    ax.invert_yaxis()\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "plot_polygons(pp_raw, col_raw, f\"Raw polygons ({total_verts_raw} vertices)\", ax)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Douglas-Peucker simplification\n",
+    "\n",
+    "Increasing tolerance values show the trade-off between fidelity and vertex count. The tolerance is in coordinate units (pixels here)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(18, 6))\n",
+    "\n",
+    "for ax, tol in zip(axes, [0.5, 1.5, 3.0]):\n",
+    "    col, pp = polygonize(raster, simplify_tolerance=tol)\n",
+    "    n_verts = sum(len(r) for rings in pp for r in rings)\n",
+    "    reduction = 100 * (1 - n_verts / total_verts_raw)\n",
+    "    plot_polygons(pp, col,\n",
+    "                  f\"tolerance={tol}  ({n_verts} verts, {reduction:.0f}% reduction)\",\n",
+    "                  ax)\n",
+    "\n",
+    "plt.suptitle(\"Douglas-Peucker simplification\", fontsize=14, y=1.02)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visvalingam-Whyatt simplification\n",
+    "\n",
+    "Area-based simplification removes vertices that contribute the least area change. The tolerance is a minimum triangle area threshold."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(18, 6))\n",
+    "\n",
+    "for ax, tol in zip(axes, [0.25, 1.0, 3.0]):\n",
+    "    col, pp = polygonize(raster, simplify_tolerance=tol,\n",
+    "                         simplify_method=\"visvalingam-whyatt\")\n",
+    "    n_verts = sum(len(r) for rings in pp for r in rings)\n",
+    "    reduction = 100 * (1 - n_verts / total_verts_raw)\n",
+    "    plot_polygons(pp, col,\n",
+    "                  f\"tolerance={tol}  ({n_verts} verts, {reduction:.0f}% reduction)\",\n",
+    "                  ax)\n",
+    "\n",
+    "plt.suptitle(\"Visvalingam-Whyatt simplification\", fontsize=14, y=1.02)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## GeoDataFrame output\n",
+    "\n",
+    "`simplify_tolerance` works with all return types including GeoDataFrame."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf = polygonize(raster, simplify_tolerance=1.5, return_type=\"geopandas\",\n",
+    "                 column_name=\"landcover\")\n",
+    "print(gdf.head(10))\n",
+    "print(f\"\\nAll geometries valid: {gdf.geometry.is_valid.all()}\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}