Skip to content

Commit d8108a6

Browse files
brendancolbrendancol
committed
Add contour extraction user guide notebook (#964)
1 parent 2bbac27 commit d8108a6

File tree

1 file changed

+271
-0
lines changed

1 file changed

+271
-0
lines changed

‎examples/user_guide/22_Contour_Lines.ipynb‎

Lines changed: 271 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,271 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "title",
6+
"metadata": {},
7+
"source": [
8+
"# Xarray-spatial\n",
9+
"### User Guide: Contour Line Extraction\n",
10+
"-----\n",
11+
"\n",
12+
"Contour lines (isolines) connect points of equal value on a surface. They are one of the most common ways to represent elevation on topographic maps.\n",
13+
"\n",
14+
"The `contours` function uses a marching squares algorithm to trace isolines through a 2D raster at specified elevation values. It supports NumPy, CuPy, Dask, and Dask+CuPy backends.\n",
15+
"\n",
16+
"**Topics covered:**\n",
17+
"- [Basic contour extraction](#Basic-Contour-Extraction)\n",
18+
"- [Automatic level selection](#Automatic-Level-Selection)\n",
19+
"- [GeoDataFrame output](#GeoDataFrame-Output)\n",
20+
"- [Overlaying contours on terrain](#Contours-on-Terrain)"
21+
]
22+
},
23+
{
24+
"cell_type": "code",
25+
"execution_count": null,
26+
"id": "imports",
27+
"metadata": {},
28+
"outputs": [],
29+
"source": [
30+
"import numpy as np\n",
31+
"import xarray as xr\n",
32+
"import matplotlib.pyplot as plt\n",
33+
"\n",
34+
"import xrspatial\n",
35+
"from xrspatial import contours\n",
36+
"from xrspatial.terrain import generate_terrain"
37+
]
38+
},
39+
{
40+
"cell_type": "markdown",
41+
"id": "data-header",
42+
"metadata": {},
43+
"source": [
44+
"## Generate Synthetic Terrain\n",
45+
"\n",
46+
"We start by generating a synthetic elevation raster using `generate_terrain`."
47+
]
48+
},
49+
{
50+
"cell_type": "code",
51+
"execution_count": null,
52+
"id": "gen-terrain",
53+
"metadata": {},
54+
"outputs": [],
55+
"source": [
56+
"W, H = 400, 300\n",
57+
"terrain = xr.DataArray(np.zeros((H, W)))\n",
58+
"terrain = generate_terrain(terrain)\n",
59+
"\n",
60+
"fig, ax = plt.subplots(figsize=(10, 7))\n",
61+
"im = ax.imshow(terrain.values, cmap='terrain', origin='lower')\n",
62+
"plt.colorbar(im, ax=ax, label='Elevation')\n",
63+
"ax.set_title('Synthetic Terrain')\n",
64+
"plt.tight_layout()\n",
65+
"plt.show()"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"id": "basic-header",
71+
"metadata": {},
72+
"source": [
73+
"## Basic Contour Extraction\n",
74+
"\n",
75+
"Extract contour lines at specific elevation levels. The function returns a list of `(level, coordinates)` tuples, where each `coordinates` array has shape `(N, 2)` with `(row, col)` positions."
76+
]
77+
},
78+
{
79+
"cell_type": "code",
80+
"execution_count": null,
81+
"id": "basic-contours",
82+
"metadata": {},
83+
"outputs": [],
84+
"source": [
85+
"# Pick a few specific levels\n",
86+
"vmin, vmax = float(np.nanmin(terrain.values)), float(np.nanmax(terrain.values))\n",
87+
"levels = np.linspace(vmin, vmax, 12)[1:-1] # 10 interior levels\n",
88+
"\n",
89+
"lines = contours(terrain, levels=levels)\n",
90+
"print(f'Extracted {len(lines)} contour polylines across {len(levels)} levels')\n",
91+
"\n",
92+
"# Show the first result\n",
93+
"level, coords = lines[0]\n",
94+
"print(f'First line: level={level:.1f}, {len(coords)} vertices')"
95+
]
96+
},
97+
{
98+
"cell_type": "code",
99+
"execution_count": null,
100+
"id": "plot-basic",
101+
"metadata": {},
102+
"outputs": [],
103+
"source": [
104+
"fig, ax = plt.subplots(figsize=(10, 7))\n",
105+
"im = ax.imshow(terrain.values, cmap='terrain', origin='lower', alpha=0.5)\n",
106+
"\n",
107+
"# Color contours by level\n",
108+
"unique_levels = sorted(set(lvl for lvl, _ in lines))\n",
109+
"cmap = plt.cm.inferno\n",
110+
"norm = plt.Normalize(vmin=min(unique_levels), vmax=max(unique_levels))\n",
111+
"\n",
112+
"for level, coords in lines:\n",
113+
" ax.plot(coords[:, 1], coords[:, 0], color=cmap(norm(level)), linewidth=0.8)\n",
114+
"\n",
115+
"sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
116+
"plt.colorbar(sm, ax=ax, label='Elevation')\n",
117+
"ax.set_title('Contour Lines on Terrain')\n",
118+
"plt.tight_layout()\n",
119+
"plt.show()"
120+
]
121+
},
122+
{
123+
"cell_type": "markdown",
124+
"id": "auto-header",
125+
"metadata": {},
126+
"source": [
127+
"## Automatic Level Selection\n",
128+
"\n",
129+
"When `levels` is not specified, the function automatically picks `n_levels` evenly spaced values between the raster's min and max."
130+
]
131+
},
132+
{
133+
"cell_type": "code",
134+
"execution_count": null,
135+
"id": "auto-levels",
136+
"metadata": {},
137+
"outputs": [],
138+
"source": [
139+
"# Let the function choose 20 levels automatically\n",
140+
"auto_lines = contours(terrain, n_levels=20)\n",
141+
"\n",
142+
"fig, ax = plt.subplots(figsize=(10, 7))\n",
143+
"ax.imshow(terrain.values, cmap='terrain', origin='lower', alpha=0.4)\n",
144+
"\n",
145+
"for level, coords in auto_lines:\n",
146+
" ax.plot(coords[:, 1], coords[:, 0], 'k-', linewidth=0.5)\n",
147+
"\n",
148+
"ax.set_title('Automatic Contour Levels (n_levels=20)')\n",
149+
"plt.tight_layout()\n",
150+
"plt.show()"
151+
]
152+
},
153+
{
154+
"cell_type": "markdown",
155+
"id": "gdf-header",
156+
"metadata": {},
157+
"source": [
158+
"## GeoDataFrame Output\n",
159+
"\n",
160+
"Set `return_type='geopandas'` to get a GeoDataFrame with `level` and `geometry` columns. This is useful for further spatial analysis or export to GIS formats."
161+
]
162+
},
163+
{
164+
"cell_type": "code",
165+
"execution_count": null,
166+
"id": "gdf-output",
167+
"metadata": {},
168+
"outputs": [],
169+
"source": [
170+
"gdf = contours(terrain, levels=levels, return_type='geopandas')\n",
171+
"print(f'GeoDataFrame with {len(gdf)} contour lines')\n",
172+
"gdf.head(10)"
173+
]
174+
},
175+
{
176+
"cell_type": "code",
177+
"execution_count": null,
178+
"id": "gdf-plot",
179+
"metadata": {},
180+
"outputs": [],
181+
"source": [
182+
"fig, ax = plt.subplots(figsize=(10, 7))\n",
183+
"gdf.plot(ax=ax, column='level', cmap='coolwarm', linewidth=0.8, legend=True)\n",
184+
"ax.set_title('Contour GeoDataFrame')\n",
185+
"ax.set_aspect('equal')\n",
186+
"plt.tight_layout()\n",
187+
"plt.show()"
188+
]
189+
},
190+
{
191+
"cell_type": "markdown",
192+
"id": "overlay-header",
193+
"metadata": {},
194+
"source": [
195+
"## Contours on Terrain\n",
196+
"\n",
197+
"Combine hillshade with contour lines for a topographic map look."
198+
]
199+
},
200+
{
201+
"cell_type": "code",
202+
"execution_count": null,
203+
"id": "overlay-plot",
204+
"metadata": {},
205+
"outputs": [],
206+
"source": [
207+
"hillshade = terrain.xrs.hillshade()\n",
208+
"\n",
209+
"fig, ax = plt.subplots(figsize=(10, 7))\n",
210+
"ax.imshow(hillshade.values, cmap='gray', origin='lower')\n",
211+
"ax.imshow(terrain.values, cmap='terrain', origin='lower', alpha=0.35)\n",
212+
"\n",
213+
"# Draw contours with labeled index contours (every 5th level thicker)\n",
214+
"for i, (level, coords) in enumerate(auto_lines):\n",
215+
" lw = 1.2 if i % 5 == 0 else 0.4\n",
216+
" ax.plot(coords[:, 1], coords[:, 0], 'sienna', linewidth=lw, alpha=0.7)\n",
217+
"\n",
218+
"ax.set_title('Topographic Map: Hillshade + Contours')\n",
219+
"ax.set_xlabel('Column')\n",
220+
"ax.set_ylabel('Row')\n",
221+
"plt.tight_layout()\n",
222+
"plt.show()"
223+
]
224+
},
225+
{
226+
"cell_type": "markdown",
227+
"id": "accessor-header",
228+
"metadata": {},
229+
"source": [
230+
"## Using the Accessor\n",
231+
"\n",
232+
"The `contours` function is also available via the `.xrs` accessor on any DataArray."
233+
]
234+
},
235+
{
236+
"cell_type": "code",
237+
"execution_count": null,
238+
"id": "accessor-demo",
239+
"metadata": {},
240+
"outputs": [],
241+
"source": [
242+
"# Equivalent to contours(terrain, levels=[500])\n",
243+
"accessor_lines = terrain.xrs.contours(levels=levels[:3])\n",
244+
"print(f'{len(accessor_lines)} contour lines via .xrs accessor')"
245+
]
246+
},
247+
{
248+
"cell_type": "markdown",
249+
"id": "references",
250+
"metadata": {},
251+
"source": [
252+
"### References\n",
253+
"- Marching squares algorithm: https://en.wikipedia.org/wiki/Marching_squares\n",
254+
"- Contour lines in cartography: https://en.wikipedia.org/wiki/Contour_line"
255+
]
256+
}
257+
],
258+
"metadata": {
259+
"kernelspec": {
260+
"display_name": "Python 3",
261+
"language": "python",
262+
"name": "python3"
263+
},
264+
"language_info": {
265+
"name": "python",
266+
"version": "3.10.0"
267+
}
268+
},
269+
"nbformat": 4,
270+
"nbformat_minor": 5
271+
}

0 commit comments

Comments
 (0)