Skip to content

Commit ce4f0d8

Browse files
brendancolbrendancol
authored
Add polygon clipping function (#1144) (#1173)
`clip_polygon()` masks pixels outside a polygon to nodata, with optional bbox crop. Supports all four backends. Accepts shapely geometries, GeoDataFrames, or coordinate lists. Adds 20 tests, API docs, user guide notebook, and README entry.
1 parent 871a51c commit ce4f0d8

File tree

7 files changed

+845
-0
lines changed

7 files changed

+845
-0
lines changed

‎README.md‎

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

‎docs/source/reference/zonal.rst‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,13 @@ Apply
1616

1717
xrspatial.zonal.apply
1818

19+
Clip Polygon
20+
============
21+
.. autosummary::
22+
:toctree: _autosummary
23+
24+
xrspatial.polygon_clip.clip_polygon
25+
1926
Crop
2027
====
2128
.. autosummary::

‎examples/user_guide/50_Polygon_Clipping.ipynb‎

Lines changed: 311 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,311 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "a1",
6+
"metadata": {},
7+
"source": [
8+
"# Polygon Clipping\n",
9+
"\n",
10+
"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",
11+
"\n",
12+
"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."
13+
]
14+
},
15+
{
16+
"cell_type": "markdown",
17+
"id": "a2",
18+
"metadata": {},
19+
"source": [
20+
"### What you'll see\n",
21+
"\n",
22+
"1. Generate a synthetic terrain raster\n",
23+
"2. Clip to a single polygon\n",
24+
"3. Compare crop=True vs crop=False\n",
25+
"4. Clip with a custom nodata value\n",
26+
"5. Use all_touched to include boundary pixels\n",
27+
"6. Clip with a list of coordinate pairs"
28+
]
29+
},
30+
{
31+
"cell_type": "code",
32+
"execution_count": null,
33+
"id": "a3",
34+
"metadata": {},
35+
"outputs": [],
36+
"source": [
37+
"import numpy as np\n",
38+
"import xarray as xr\n",
39+
"import matplotlib.pyplot as plt\n",
40+
"from matplotlib.patches import Polygon as MplPolygon\n",
41+
"from shapely.geometry import Polygon\n",
42+
"\n",
43+
"from xrspatial import generate_terrain\n",
44+
"from xrspatial.polygon_clip import clip_polygon"
45+
]
46+
},
47+
{
48+
"cell_type": "markdown",
49+
"id": "a4",
50+
"metadata": {},
51+
"source": [
52+
"## Generate terrain\n",
53+
"\n",
54+
"A 200x300 synthetic elevation raster gives us something realistic to clip."
55+
]
56+
},
57+
{
58+
"cell_type": "code",
59+
"execution_count": null,
60+
"id": "a5",
61+
"metadata": {},
62+
"outputs": [],
63+
"source": [
64+
"template = xr.DataArray(np.zeros((200, 300)))\n",
65+
"terrain = generate_terrain(template, x_range=(0, 30), y_range=(0, 20))\n",
66+
"\n",
67+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
68+
"terrain.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
69+
"ax.set_title('Synthetic terrain')\n",
70+
"ax.set_aspect('equal')\n",
71+
"plt.tight_layout()\n",
72+
"plt.show()\n",
73+
"\n",
74+
"print(f'Shape: {terrain.shape}')"
75+
]
76+
},
77+
{
78+
"cell_type": "markdown",
79+
"id": "a6",
80+
"metadata": {},
81+
"source": [
82+
"## Define a clipping polygon\n",
83+
"\n",
84+
"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."
85+
]
86+
},
87+
{
88+
"cell_type": "code",
89+
"execution_count": null,
90+
"id": "a7",
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"study_area = Polygon([\n",
95+
" (5, 4), (8, 2), (14, 3), (20, 5),\n",
96+
" (22, 10), (18, 16), (12, 17),\n",
97+
" (6, 14), (3, 9),\n",
98+
"])\n",
99+
"\n",
100+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
101+
"terrain.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
102+
"patch = MplPolygon(\n",
103+
" list(study_area.exterior.coords),\n",
104+
" closed=True, fill=False, edgecolor='red', linewidth=2,\n",
105+
")\n",
106+
"ax.add_patch(patch)\n",
107+
"ax.set_title('Terrain with study area boundary')\n",
108+
"ax.set_aspect('equal')\n",
109+
"plt.tight_layout()\n",
110+
"plt.show()"
111+
]
112+
},
113+
{
114+
"cell_type": "markdown",
115+
"id": "a8",
116+
"metadata": {},
117+
"source": [
118+
"## Basic clip\n",
119+
"\n",
120+
"With `crop=True` (the default), the result is trimmed to the polygon's bounding box. Pixels outside the polygon are set to NaN."
121+
]
122+
},
123+
{
124+
"cell_type": "code",
125+
"execution_count": null,
126+
"id": "a9",
127+
"metadata": {},
128+
"outputs": [],
129+
"source": [
130+
"clipped = clip_polygon(terrain, study_area)\n",
131+
"\n",
132+
"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
133+
"\n",
134+
"terrain.plot(ax=axes[0], cmap='terrain', add_colorbar=False)\n",
135+
"axes[0].set_title('Original')\n",
136+
"axes[0].set_aspect('equal')\n",
137+
"\n",
138+
"clipped.plot(ax=axes[1], cmap='terrain', add_colorbar=False)\n",
139+
"axes[1].set_title('Clipped (crop=True)')\n",
140+
"axes[1].set_aspect('equal')\n",
141+
"\n",
142+
"plt.tight_layout()\n",
143+
"plt.show()\n",
144+
"\n",
145+
"print(f'Original shape: {terrain.shape}')\n",
146+
"print(f'Clipped shape: {clipped.shape}')"
147+
]
148+
},
149+
{
150+
"cell_type": "markdown",
151+
"id": "a10",
152+
"metadata": {},
153+
"source": [
154+
"## crop=False: keep the original extent\n",
155+
"\n",
156+
"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."
157+
]
158+
},
159+
{
160+
"cell_type": "code",
161+
"execution_count": null,
162+
"id": "a11",
163+
"metadata": {},
164+
"outputs": [],
165+
"source": [
166+
"clipped_full = clip_polygon(terrain, study_area, crop=False)\n",
167+
"\n",
168+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
169+
"clipped_full.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
170+
"ax.set_title('Clipped (crop=False) — same extent, masked outside polygon')\n",
171+
"ax.set_aspect('equal')\n",
172+
"plt.tight_layout()\n",
173+
"plt.show()\n",
174+
"\n",
175+
"print(f'Shape matches original: {clipped_full.shape == terrain.shape}')"
176+
]
177+
},
178+
{
179+
"cell_type": "markdown",
180+
"id": "a12",
181+
"metadata": {},
182+
"source": [
183+
"## Custom nodata value\n",
184+
"\n",
185+
"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."
186+
]
187+
},
188+
{
189+
"cell_type": "code",
190+
"execution_count": null,
191+
"id": "a13",
192+
"metadata": {},
193+
"outputs": [],
194+
"source": [
195+
"clipped_nd = clip_polygon(terrain, study_area, nodata=-9999, crop=False)\n",
196+
"\n",
197+
"# Count cells with each status\n",
198+
"vals = clipped_nd.values\n",
199+
"n_nodata = np.sum(vals == -9999)\n",
200+
"n_valid = np.sum(vals != -9999)\n",
201+
"print(f'Valid pixels: {n_valid}')\n",
202+
"print(f'Nodata pixels: {n_nodata}')"
203+
]
204+
},
205+
{
206+
"cell_type": "markdown",
207+
"id": "a14",
208+
"metadata": {},
209+
"source": [
210+
"## all_touched: boundary pixel inclusion\n",
211+
"\n",
212+
"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."
213+
]
214+
},
215+
{
216+
"cell_type": "code",
217+
"execution_count": null,
218+
"id": "a15",
219+
"metadata": {},
220+
"outputs": [],
221+
"source": [
222+
"clipped_default = clip_polygon(terrain, study_area, crop=False)\n",
223+
"clipped_touched = clip_polygon(terrain, study_area, crop=False, all_touched=True)\n",
224+
"\n",
225+
"n_default = np.count_nonzero(np.isfinite(clipped_default.values))\n",
226+
"n_touched = np.count_nonzero(np.isfinite(clipped_touched.values))\n",
227+
"\n",
228+
"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
229+
"\n",
230+
"clipped_default.plot(ax=axes[0], cmap='terrain', add_colorbar=False)\n",
231+
"axes[0].set_title(f'all_touched=False ({n_default} pixels)')\n",
232+
"axes[0].set_aspect('equal')\n",
233+
"\n",
234+
"clipped_touched.plot(ax=axes[1], cmap='terrain', add_colorbar=False)\n",
235+
"axes[1].set_title(f'all_touched=True ({n_touched} pixels)')\n",
236+
"axes[1].set_aspect('equal')\n",
237+
"\n",
238+
"plt.tight_layout()\n",
239+
"plt.show()\n",
240+
"\n",
241+
"print(f'all_touched adds {n_touched - n_default} boundary pixels')"
242+
]
243+
},
244+
{
245+
"cell_type": "markdown",
246+
"id": "a16",
247+
"metadata": {},
248+
"source": [
249+
"## Coordinate array input\n",
250+
"\n",
251+
"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."
252+
]
253+
},
254+
{
255+
"cell_type": "code",
256+
"execution_count": null,
257+
"id": "a17",
258+
"metadata": {},
259+
"outputs": [],
260+
"source": [
261+
"coords = [(5, 4), (8, 2), (14, 3), (20, 5),\n",
262+
" (22, 10), (18, 16), (12, 17), (6, 14), (3, 9)]\n",
263+
"\n",
264+
"clipped_coords = clip_polygon(terrain, coords)\n",
265+
"\n",
266+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
267+
"clipped_coords.plot(ax=ax, cmap='terrain', add_colorbar=True)\n",
268+
"ax.set_title('Clipped from coordinate list')\n",
269+
"ax.set_aspect('equal')\n",
270+
"plt.tight_layout()\n",
271+
"plt.show()"
272+
]
273+
},
274+
{
275+
"cell_type": "markdown",
276+
"id": "a18",
277+
"metadata": {},
278+
"source": [
279+
"## Accessor syntax\n",
280+
"\n",
281+
"All xrspatial functions are also available through the `.xrs` accessor on DataArrays."
282+
]
283+
},
284+
{
285+
"cell_type": "code",
286+
"execution_count": null,
287+
"id": "a19",
288+
"metadata": {},
289+
"outputs": [],
290+
"source": [
291+
"import xrspatial # registers the .xrs accessor\n",
292+
"\n",
293+
"clipped_acc = terrain.xrs.clip_polygon(study_area)\n",
294+
"print(f'Accessor result shape: {clipped_acc.shape}')"
295+
]
296+
}
297+
],
298+
"metadata": {
299+
"kernelspec": {
300+
"display_name": "Python 3",
301+
"language": "python",
302+
"name": "python3"
303+
},
304+
"language_info": {
305+
"name": "python",
306+
"version": "3.14.2"
307+
}
308+
},
309+
"nbformat": 4,
310+
"nbformat_minor": 5
311+
}

‎xrspatial/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,6 +121,7 @@
121121
from xrspatial.terrain_metrics import tri # noqa
122122
from xrspatial.hydro import twi # noqa: unified wrapper
123123
from xrspatial.hydro import twi_d8 # noqa
124+
from xrspatial.polygon_clip import clip_polygon # noqa
124125
from xrspatial.polygonize import polygonize # noqa
125126
from xrspatial.viewshed import viewshed # noqa
126127
from xrspatial.visibility import cumulative_viewshed # noqa

‎xrspatial/accessor.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -362,6 +362,10 @@ def crop(self, zones, zones_ids, **kwargs):
362362
from .zonal import crop
363363
return crop(zones, self._obj, zones_ids, **kwargs)
364364

365+
def clip_polygon(self, geometry, **kwargs):
366+
from .polygon_clip import clip_polygon
367+
return clip_polygon(self._obj, geometry, **kwargs)
368+
365369
def trim(self, **kwargs):
366370
from .zonal import trim
367371
return trim(self._obj, **kwargs)

0 commit comments

Comments
 (0)