Add polygon clipping function (#1144)

README.md

@@ -460,6 +460,7 @@ Same-CRS tiles skip reprojection entirely and are placed by direct coordinate al
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
 | [Apply](xrspatial/zonal.py) | Applies a custom function to each zone in a classified raster | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Clip Polygon](xrspatial/polygon_clip.py) | Clips a raster to an arbitrary polygon with masking | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Crop](xrspatial/zonal.py) | Extracts the bounding rectangle of a specific zone | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Regions](xrspatial/zonal.py) | Identifies connected regions of non-zero cells | Standard (CCL) | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Trim](xrspatial/zonal.py) | Removes nodata border rows and columns from a raster | Standard | ✅️ | ✅️ | ✅️ | ✅️ |

docs/source/reference/zonal.rst

@@ -16,6 +16,13 @@ Apply
 
     xrspatial.zonal.apply
 
+Clip Polygon
+============
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.polygon_clip.clip_polygon
+
 Crop
 ====
 .. autosummary::

examples/user_guide/50_Polygon_Clipping.ipynb

@@ -0,0 +1,311 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a1",
+   "metadata": {},
+   "source": [
+    "# Polygon Clipping\n",
+    "\n",
+    "Extracting raster data for an irregular boundary is one of the most common GIS operations. `clip_polygon()` masks pixels that fall outside an arbitrary polygon to a nodata value (NaN by default) and optionally trims the result to the polygon's bounding box.\n",
+    "\n",
+    "This is different from `crop()`, which only extracts a rectangular bounding box. When your study area is a watershed, administrative boundary, or any non-rectangular shape, `clip_polygon()` is what you want."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a2",
+   "metadata": {},
+   "source": [
+    "### What you'll see\n",
+    "\n",
+    "1. Generate a synthetic terrain raster\n",
+    "2. Clip to a single polygon\n",
+    "3. Compare crop=True vs crop=False\n",
+    "4. Clip with a custom nodata value\n",
+    "5. Use all_touched to include boundary pixels\n",
+    "6. Clip with a list of coordinate pairs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3",
+   "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 shapely.geometry import Polygon\n",
+    "\n",
+    "from xrspatial import generate_terrain\n",
+    "from xrspatial.polygon_clip import clip_polygon"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a4",
+   "metadata": {},
+   "source": [
+    "## Generate terrain\n",
+    "\n",
+    "A 200x300 synthetic elevation raster gives us something realistic to clip."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "template = xr.DataArray(np.zeros((200, 300)))\n",
+    "terrain = generate_terrain(template, x_range=(0, 30), y_range=(0, 20))\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "terrain.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
+    "ax.set_title('Synthetic terrain')\n",
+    "ax.set_aspect('equal')\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "print(f'Shape: {terrain.shape}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a6",
+   "metadata": {},
+   "source": [
+    "## Define a clipping polygon\n",
+    "\n",
+    "An irregular polygon representing a study area boundary. We'll draw it on top of the terrain to see what we're about to extract."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "study_area = Polygon([\n",
+    "    (5, 4), (8, 2), (14, 3), (20, 5),\n",
+    "    (22, 10), (18, 16), (12, 17),\n",
+    "    (6, 14), (3, 9),\n",
+    "])\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "terrain.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
+    "patch = MplPolygon(\n",
+    "    list(study_area.exterior.coords),\n",
+    "    closed=True, fill=False, edgecolor='red', linewidth=2,\n",
+    ")\n",
+    "ax.add_patch(patch)\n",
+    "ax.set_title('Terrain with study area boundary')\n",
+    "ax.set_aspect('equal')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a8",
+   "metadata": {},
+   "source": [
+    "## Basic clip\n",
+    "\n",
+    "With `crop=True` (the default), the result is trimmed to the polygon's bounding box. Pixels outside the polygon are set to NaN."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "clipped = clip_polygon(terrain, study_area)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
+    "\n",
+    "terrain.plot(ax=axes[0], cmap='terrain', add_colorbar=False)\n",
+    "axes[0].set_title('Original')\n",
+    "axes[0].set_aspect('equal')\n",
+    "\n",
+    "clipped.plot(ax=axes[1], cmap='terrain', add_colorbar=False)\n",
+    "axes[1].set_title('Clipped (crop=True)')\n",
+    "axes[1].set_aspect('equal')\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "print(f'Original shape: {terrain.shape}')\n",
+    "print(f'Clipped shape:  {clipped.shape}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a10",
+   "metadata": {},
+   "source": [
+    "## crop=False: keep the original extent\n",
+    "\n",
+    "When `crop=False`, the output has the same shape as the input. Pixels outside the polygon are still masked to NaN, but the spatial extent is unchanged. This is useful when you need the result to align pixel-for-pixel with other rasters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "clipped_full = clip_polygon(terrain, study_area, crop=False)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "clipped_full.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
+    "ax.set_title('Clipped (crop=False) — same extent, masked outside polygon')\n",
+    "ax.set_aspect('equal')\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "print(f'Shape matches original: {clipped_full.shape == terrain.shape}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a12",
+   "metadata": {},
+   "source": [
+    "## Custom nodata value\n",
+    "\n",
+    "By default, masked pixels get NaN. You can set a different value with the `nodata` parameter. This is helpful for integer rasters or when downstream tools expect a specific sentinel value."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "clipped_nd = clip_polygon(terrain, study_area, nodata=-9999, crop=False)\n",
+    "\n",
+    "# Count cells with each status\n",
+    "vals = clipped_nd.values\n",
+    "n_nodata = np.sum(vals == -9999)\n",
+    "n_valid = np.sum(vals != -9999)\n",
+    "print(f'Valid pixels:  {n_valid}')\n",
+    "print(f'Nodata pixels: {n_nodata}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a14",
+   "metadata": {},
+   "source": [
+    "## all_touched: boundary pixel inclusion\n",
+    "\n",
+    "By default, only pixels whose centre falls inside the polygon are kept. With `all_touched=True`, any pixel that the polygon boundary touches is also included. This gives a slightly larger footprint."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "clipped_default = clip_polygon(terrain, study_area, crop=False)\n",
+    "clipped_touched = clip_polygon(terrain, study_area, crop=False, all_touched=True)\n",
+    "\n",
+    "n_default = np.count_nonzero(np.isfinite(clipped_default.values))\n",
+    "n_touched = np.count_nonzero(np.isfinite(clipped_touched.values))\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
+    "\n",
+    "clipped_default.plot(ax=axes[0], cmap='terrain', add_colorbar=False)\n",
+    "axes[0].set_title(f'all_touched=False ({n_default} pixels)')\n",
+    "axes[0].set_aspect('equal')\n",
+    "\n",
+    "clipped_touched.plot(ax=axes[1], cmap='terrain', add_colorbar=False)\n",
+    "axes[1].set_title(f'all_touched=True ({n_touched} pixels)')\n",
+    "axes[1].set_aspect('equal')\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "print(f'all_touched adds {n_touched - n_default} boundary pixels')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a16",
+   "metadata": {},
+   "source": [
+    "## Coordinate array input\n",
+    "\n",
+    "If you don't have shapely installed in your workflow, you can pass the polygon as a list of (x, y) coordinate pairs. `clip_polygon()` builds the polygon internally."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coords = [(5, 4), (8, 2), (14, 3), (20, 5),\n",
+    "          (22, 10), (18, 16), (12, 17), (6, 14), (3, 9)]\n",
+    "\n",
+    "clipped_coords = clip_polygon(terrain, coords)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "clipped_coords.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
+    "ax.set_title('Clipped from coordinate list')\n",
+    "ax.set_aspect('equal')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a18",
+   "metadata": {},
+   "source": [
+    "## Accessor syntax\n",
+    "\n",
+    "All xrspatial functions are also available through the `.xrs` accessor on DataArrays."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import xrspatial  # registers the .xrs accessor\n",
+    "\n",
+    "clipped_acc = terrain.xrs.clip_polygon(study_area)\n",
+    "print(f'Accessor result shape: {clipped_acc.shape}')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.14.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

xrspatial/__init__.py

@@ -120,6 +120,7 @@
 from xrspatial.terrain_metrics import tri  # noqa
 from xrspatial.hydro import twi  # noqa: unified wrapper
 from xrspatial.hydro import twi_d8  # noqa
+from xrspatial.polygon_clip import clip_polygon  # noqa
 from xrspatial.polygonize import polygonize  # noqa
 from xrspatial.viewshed import viewshed  # noqa
 from xrspatial.visibility import cumulative_viewshed  # noqa

xrspatial/accessor.py

@@ -362,6 +362,10 @@ def crop(self, zones, zones_ids, **kwargs):
         from .zonal import crop
         return crop(zones, self._obj, zones_ids, **kwargs)
 
+    def clip_polygon(self, geometry, **kwargs):
+        from .polygon_clip import clip_polygon
+        return clip_polygon(self._obj, geometry, **kwargs)
+
     def trim(self, **kwargs):
         from .zonal import trim
         return trim(self._obj, **kwargs)