Skip to content

Commit 699bca0

Browse files
brendancolbrendancol
authored
Add TPI-based landform classification (#950) (#957)
* Add landforms() for TPI-based terrain classification (#950) Implements the Weiss (2001) 10-class landform classification scheme. Computes TPI at two neighborhood scales using NaN-aware circular convolution, standardizes to z-scores, and classifies each cell based on relative position and slope. Supports all four backends (numpy, cupy, dask+numpy, dask+cupy). * Add landforms to API docs and README feature matrix (#950) * Add landform classification user guide notebook (#950)
1 parent 583883e commit 699bca0

File tree

6 files changed

+602
-1
lines changed

6 files changed

+602
-1
lines changed

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,6 +256,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
256256
| [Terrain Generation](xrspatial/terrain.py) | Generates synthetic terrain from fBm or ridged fractal noise with optional domain warping, Worley blending, and hydraulic erosion | ✅️ | ✅️ | ✅️ | ✅️ |
257257
| [TPI](xrspatial/terrain_metrics.py) | Computes Topographic Position Index (center minus mean of neighbors) | ✅️ | ✅️ | ✅️ | ✅️ |
258258
| [TRI](xrspatial/terrain_metrics.py) | Computes Terrain Ruggedness Index (local elevation variation) | ✅️ | ✅️ | ✅️ | ✅️ |
259+
| [Landforms](xrspatial/terrain_metrics.py) | Classifies terrain into 10 landform types using the Weiss (2001) TPI scheme | ✅️ | ✅️ | ✅️ | ✅️ |
259260
| [Viewshed](xrspatial/viewshed.py) | Determines visible cells from a given observer point on terrain | ✅️ | ✅️ | ✅️ | ✅️ |
260261
| [Min Observable Height](xrspatial/experimental/min_observable_height.py) | Finds the minimum observer height needed to see each cell *(experimental)* | ✅️ | | | |
261262
| [Perlin Noise](xrspatial/perlin.py) | Generates smooth continuous random noise for procedural textures | ✅️ | ✅️ | ✅️ | ✅️ |

‎docs/source/reference/terrain_metrics.rst‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,3 +24,10 @@ Terrain Ruggedness Index (TRI)
2424
:toctree: _autosummary
2525

2626
xrspatial.terrain_metrics.tri
27+
28+
Landform Classification
29+
=======================
30+
.. autosummary::
31+
:toctree: _autosummary
32+
33+
xrspatial.terrain_metrics.landforms

‎examples/user_guide/17_Landform_Classification.ipynb‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Landform Classification\n",
8+
"\n",
9+
"The `landforms()` function classifies each cell in an elevation raster into one of 10 landform categories using the Weiss (2001) TPI-based scheme. It computes Topographic Position Index (TPI) at two neighborhood scales, standardizes both to z-scores, and combines them with slope to assign classes ranging from canyons and valleys to ridges and mountain tops."
10+
]
11+
},
12+
{
13+
"cell_type": "code",
14+
"execution_count": null,
15+
"metadata": {},
16+
"outputs": [],
17+
"source": [
18+
"import numpy as np\n",
19+
"import xarray as xr\n",
20+
"import matplotlib.pyplot as plt\n",
21+
"from matplotlib.colors import ListedColormap, BoundaryNorm\n",
22+
"\n",
23+
"from xrspatial import landforms\n",
24+
"from xrspatial.terrain_metrics import LANDFORM_CLASSES"
25+
]
26+
},
27+
{
28+
"cell_type": "markdown",
29+
"metadata": {},
30+
"source": [
31+
"## 1. Generate synthetic terrain\n",
32+
"\n",
33+
"We build a surface that contains a clear ridge, a valley, and some flat/sloped areas so the classifier has a variety of features to work with."
34+
]
35+
},
36+
{
37+
"cell_type": "code",
38+
"execution_count": null,
39+
"metadata": {},
40+
"outputs": [],
41+
"source": [
42+
"rows, cols = 120, 160\n",
43+
"y = np.linspace(0, 12, rows)\n",
44+
"x = np.linspace(0, 16, cols)\n",
45+
"Y, X = np.meshgrid(y, x, indexing='ij')\n",
46+
"\n",
47+
"# Combine a ridge, a valley, and some rolling terrain\n",
48+
"ridge = 80 * np.exp(-((Y - 3)**2 + (X - 8)**2) / (2 * 2.5**2))\n",
49+
"valley = -40 * np.exp(-((Y - 9)**2 + (X - 5)**2) / (2 * 2**2))\n",
50+
"rolling = 15 * np.sin(Y / 1.5) * np.cos(X / 2)\n",
51+
"base = 200 + 10 * Y # gentle regional slope\n",
52+
"\n",
53+
"elevation = base + ridge + valley + rolling\n",
54+
"\n",
55+
"agg = xr.DataArray(\n",
56+
" elevation,\n",
57+
" dims=['y', 'x'],\n",
58+
" attrs={'res': (0.1, 0.1)},\n",
59+
")\n",
60+
"agg['y'] = np.linspace(y[-1], y[0], rows)\n",
61+
"agg['x'] = x\n",
62+
"\n",
63+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
64+
"im = ax.imshow(elevation, cmap='terrain', aspect='auto')\n",
65+
"ax.set_title('Synthetic elevation')\n",
66+
"ax.set_xlabel('Column')\n",
67+
"ax.set_ylabel('Row')\n",
68+
"fig.colorbar(im, ax=ax, label='Elevation (m)')\n",
69+
"plt.tight_layout()\n",
70+
"plt.show()"
71+
]
72+
},
73+
{
74+
"cell_type": "markdown",
75+
"metadata": {},
76+
"source": [
77+
"## 2. Run landform classification\n",
78+
"\n",
79+
"The default parameters (`inner_radius=3`, `outer_radius=15`) work well for many DEMs. Here we use smaller radii to match the scale of our synthetic terrain."
80+
]
81+
},
82+
{
83+
"cell_type": "code",
84+
"execution_count": null,
85+
"metadata": {},
86+
"outputs": [],
87+
"source": [
88+
"classes = landforms(agg, inner_radius=3, outer_radius=12)\n",
89+
"\n",
90+
"# Build a categorical colormap\n",
91+
"colors = [\n",
92+
" '#1a237e', # 1 Canyon\n",
93+
" '#4a148c', # 2 Midslope drainage\n",
94+
" '#880e4f', # 3 Upland drainage\n",
95+
" '#0d47a1', # 4 U-shaped valley\n",
96+
" '#a5d6a7', # 5 Plain\n",
97+
" '#fff9c4', # 6 Open slope\n",
98+
" '#ffcc80', # 7 Upper slope\n",
99+
" '#e65100', # 8 Local ridge\n",
100+
" '#bf360c', # 9 Midslope ridge\n",
101+
" '#b71c1c', # 10 Mountain top\n",
102+
"]\n",
103+
"cmap = ListedColormap(colors)\n",
104+
"bounds = np.arange(0.5, 11.5, 1)\n",
105+
"norm = BoundaryNorm(bounds, cmap.N)\n",
106+
"\n",
107+
"fig, ax = plt.subplots(figsize=(10, 6))\n",
108+
"im = ax.imshow(classes.values, cmap=cmap, norm=norm, aspect='auto')\n",
109+
"cbar = fig.colorbar(im, ax=ax, ticks=range(1, 11))\n",
110+
"cbar.ax.set_yticklabels(\n",
111+
" [LANDFORM_CLASSES[i] for i in range(1, 11)], fontsize=8\n",
112+
")\n",
113+
"ax.set_title('Weiss (2001) landform classification')\n",
114+
"ax.set_xlabel('Column')\n",
115+
"ax.set_ylabel('Row')\n",
116+
"plt.tight_layout()\n",
117+
"plt.show()"
118+
]
119+
},
120+
{
121+
"cell_type": "markdown",
122+
"metadata": {},
123+
"source": [
124+
"## 3. Effect of neighborhood radii\n",
125+
"\n",
126+
"The inner and outer radii control the scale of features the classifier picks up. Smaller radii detect finer-grained features; larger radii smooth out local variation."
127+
]
128+
},
129+
{
130+
"cell_type": "code",
131+
"execution_count": null,
132+
"metadata": {},
133+
"outputs": [],
134+
"source": [
135+
"configs = [\n",
136+
" (2, 6, 'inner=2, outer=6'),\n",
137+
" (3, 12, 'inner=3, outer=12'),\n",
138+
" (5, 20, 'inner=5, outer=20'),\n",
139+
"]\n",
140+
"\n",
141+
"fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
142+
"for ax, (ir, outr, label) in zip(axes, configs):\n",
143+
" c = landforms(agg, inner_radius=ir, outer_radius=outr)\n",
144+
" ax.imshow(c.values, cmap=cmap, norm=norm, aspect='auto')\n",
145+
" ax.set_title(label)\n",
146+
" ax.axis('off')\n",
147+
"plt.suptitle('Scale sensitivity', y=1.02)\n",
148+
"plt.tight_layout()\n",
149+
"plt.show()"
150+
]
151+
},
152+
{
153+
"cell_type": "markdown",
154+
"metadata": {},
155+
"source": [
156+
"## 4. Slope threshold for plains vs. open slopes\n",
157+
"\n",
158+
"The `slope_threshold` parameter (default 5 degrees) controls whether mid-position cells are labeled as plains or open slopes. A lower threshold classifies more cells as slopes; a higher threshold classifies more as plains."
159+
]
160+
},
161+
{
162+
"cell_type": "code",
163+
"execution_count": null,
164+
"metadata": {},
165+
"outputs": [],
166+
"source": [
167+
"thresholds = [2.0, 5.0, 15.0]\n",
168+
"\n",
169+
"fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
170+
"for ax, thr in zip(axes, thresholds):\n",
171+
" c = landforms(agg, inner_radius=3, outer_radius=12,\n",
172+
" slope_threshold=thr)\n",
173+
" ax.imshow(c.values, cmap=cmap, norm=norm, aspect='auto')\n",
174+
" ax.set_title(f'slope_threshold={thr}')\n",
175+
" ax.axis('off')\n",
176+
"plt.suptitle('Plains vs. open slopes threshold', y=1.02)\n",
177+
"plt.tight_layout()\n",
178+
"plt.show()"
179+
]
180+
},
181+
{
182+
"cell_type": "markdown",
183+
"metadata": {},
184+
"source": [
185+
"## 5. Class distribution\n",
186+
"\n",
187+
"A histogram of class frequencies shows which landform types dominate the study area."
188+
]
189+
},
190+
{
191+
"cell_type": "code",
192+
"execution_count": null,
193+
"metadata": {},
194+
"outputs": [],
195+
"source": [
196+
"classes = landforms(agg, inner_radius=3, outer_radius=12)\n",
197+
"valid = classes.values[~np.isnan(classes.values)].astype(int)\n",
198+
"\n",
199+
"counts = [np.sum(valid == i) for i in range(1, 11)]\n",
200+
"labels = [LANDFORM_CLASSES[i] for i in range(1, 11)]\n",
201+
"\n",
202+
"fig, ax = plt.subplots(figsize=(10, 5))\n",
203+
"ax.barh(range(10), counts, color=colors)\n",
204+
"ax.set_yticks(range(10))\n",
205+
"ax.set_yticklabels(labels, fontsize=9)\n",
206+
"ax.set_xlabel('Cell count')\n",
207+
"ax.set_title('Landform class distribution')\n",
208+
"ax.invert_yaxis()\n",
209+
"plt.tight_layout()\n",
210+
"plt.show()"
211+
]
212+
}
213+
],
214+
"metadata": {
215+
"kernelspec": {
216+
"display_name": "Python 3",
217+
"language": "python",
218+
"name": "python3"
219+
},
220+
"language_info": {
221+
"name": "python",
222+
"version": "3.10.0"
223+
}
224+
},
225+
"nbformat": 4,
226+
"nbformat_minor": 4
227+
}

‎xrspatial/__init__.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,8 @@
7070
from xrspatial.surface_distance import surface_direction # noqa
7171
from xrspatial.surface_distance import surface_distance # noqa
7272
from xrspatial.terrain import generate_terrain # noqa
73+
from xrspatial.terrain_metrics import landforms # noqa
74+
from xrspatial.terrain_metrics import LANDFORM_CLASSES # noqa
7375
from xrspatial.terrain_metrics import roughness # noqa
7476
from xrspatial.terrain_metrics import tpi # noqa
7577
from xrspatial.terrain_metrics import tri # noqa

0 commit comments

Comments
 (0)