Elevation · read a landscape at its native scale¶
Context¶
A terrain surface is useful before it becomes a time series. Start with a small real 3DEP window near Boulder and ask what can be learned without changing its grid.
Question¶
How does elevation vary within this hillside, and what is lost when we summarize it?
Pipe¶
pipe(terrain) | v.mean(dim="y", keep_dim=False)
Acquisition and input checks come before the analytical sentence. The cells below use frozen real loader outputs, so running the lesson does not depend on provider availability. For live acquisition, see the elevation reference. Download this notebook.
Analysis story¶
1. Read the landscape in its native cells¶
What is high and low within this small Boulder hillside?
from pathlib import Path
import hashlib
import json
import numpy as np
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from shapely.geometry import box
from cubedynamics import pipe, verbs as v
# Frozen outputs of the real bounded loaders, not generated measurements.
root = next(p for p in (Path.cwd(), *Path.cwd().parents)
if (p / "tests/fixtures/real_data/source_lessons").is_dir())
fixture = root / "tests/fixtures/real_data/source_lessons"
def verify_input(name):
record = json.loads((fixture / f"{name}.provenance.json").read_text())
assert record["is_synthetic"] is False
for relative, expected in record["files"].items():
assert hashlib.sha256((fixture / relative).read_bytes()).hexdigest() == expected
return record
record = verify_input("elevation")
with xr.open_dataarray(fixture / "elevation.nc", engine="scipy") as source:
terrain = source.load() # Only the retained 99 by 99 window.
assert terrain.dims == ("y", "x") and terrain.attrs["units"] == "m"
fig, ax = plt.subplots(figsize=(7, 4), layout="constrained")
terrain.plot(ax=ax, cmap="terrain", cbar_kwargs={"label": "Elevation (m; native vertical datum)"})
ax.set(title="3DEP · Boulder hillside · native cells", xlabel="Longitude (EPSG:4269)", ylabel="Latitude")
ax.ticklabel_format(useOffset=False, style="plain")
ax.xaxis.set_major_locator(MaxNLocator(4))
plt.show()
A real 3DEP window, north up. This is static terrain: no invented time axis, resampling, or vertical-datum conversion.
2. Describe relief relative to this window¶
Where is terrain above or below the local window mean?
def center_on_window_mean(surface):
# A spatial baseline, not the time-based anomaly verb or a sea-level reference.
return surface - surface.mean(("y", "x"))
relief = (pipe(terrain) | center_on_window_mean).unwrap()
np.testing.assert_allclose(relief, terrain - terrain.mean(), rtol=0, atol=1e-9)
fig, ax = plt.subplots(figsize=(7, 4), layout="constrained")
relief.plot(ax=ax, cmap="RdBu_r", center=0, cbar_kwargs={"label": "Departure from window mean (m)"})
ax.set(title="Same cells · a local relief baseline", xlabel="Longitude (EPSG:4269)", ylabel="Latitude")
ax.ticklabel_format(useOffset=False, style="plain")
ax.xaxis.set_major_locator(MaxNLocator(4))
plt.show()
The pipe changes the reference level, not the terrain. Negative departures are below the selected window mean; changing the window changes that baseline.
3. Reduce the map to a west–east profile¶
How does average elevation vary across the window?
# Collapse only y. This is a cell mean, not an area-weighted regional statistic.
profile = (pipe(terrain) | v.mean(dim="y", keep_dim=False)).unwrap()
np.testing.assert_allclose(profile, terrain.mean("y"))
fig, ax = plt.subplots(figsize=(7, 3.5), layout="constrained")
profile.plot(ax=ax, color="#236d81")
ax.set(title="Mean elevation by longitude", xlabel="Longitude (EPSG:4269)", ylabel="Elevation (m)")
ax.ticklabel_format(useOffset=False, style="plain")
ax.xaxis.set_major_locator(MaxNLocator(4))
plt.show()
A simpler profile makes the broad gradient legible but removes north–south structure. It is not a road grade, slope map, or watershed delineation.
Figure¶
Each of the three analytical steps displays its own result.
What the figure tells us¶
A simpler profile makes the broad gradient legible but removes north–south structure. It is not a road grade, slope map, or watershed delineation.
Input integrity and numerical checks are reproducible; these teaching results do not establish broader scientific suitability. Preserve the scope and source metadata when reusing them.
Data used¶
| Field | Frozen analysis input |
|---|---|
| Provider | USGS 3DEP |
| Product | 1/3 arc-second native elevation window |
| Dates | Static terrain; tile version USGS_13_n40w106_20260630 |
| Fixture | tests/fixtures/real_data/source_lessons/elevation.nc |
| Provenance record | tests/fixtures/real_data/source_lessons/elevation.provenance.json |
The source reference describes current supported scope, quality checks, and limits; the fixture record above identifies the observations used here. Data validation documents checksums and acceptance checks. The analytical baseline and thresholds belong to this story, not the provider.
Reproduce¶
Clone the repository, then run these commands from its root:
python -m pip install -e ".[vignettes]"
python scripts/run_vignettes.py docs/vignettes/elevation_landscape.ipynb
No network is needed after installation. Open the downloaded notebook in Jupyter and run all cells to see the same figures. The website executes these cells during its strict build. Environment setup and the vignette contract explain the workflow.