Roads · compare mapped networks without erasing their differences¶
Context¶
Mapped roads are vector features with provider-specific classes and segment boundaries. Two maps can look similar while answering subtly different questions.
Question¶
What road length is represented inside a small Boulder area, using each source's native meaning?
Pipe¶
pipe(frame) | within_area(area) | length_by_class("EPSG:32613")
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 roads reference. Download this notebook.
Analysis story¶
1. Inspect two native road descriptions¶
Do these sources describe the same analysis area in the same way?
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("roads")
networks = {}
for name in ("overture", "osm"):
body = json.loads((fixture / f"roads_{name}.geojson").read_text())
networks[name] = gpd.GeoDataFrame.from_features(body, crs="EPSG:4326")
networks[name].attrs.update(record["native_metadata"][name])
assert networks[name].source_feature_id.is_unique
# Use one query boundary for both displays; keep full native features.
area = box(*networks["osm"].attrs["requested_bbox"])
def maps(frames, title):
fig, axes = plt.subplots(1, 2, figsize=(10, 4), layout="constrained")
for (name, frame), ax in zip(frames.items(), axes):
frame.plot(ax=ax, color="#236d81", linewidth=1)
west, south, east, north = area.bounds
ax.set(xlim=(west, east), ylim=(south, north), title=name,
xlabel="Longitude (WGS84)", ylabel="Latitude")
ax.ticklabel_format(useOffset=False, style="plain")
ax.xaxis.set_major_locator(MaxNLocator(3))
fig.suptitle(title)
fig.text(.5, -.035, "© OpenStreetMap contributors · Overture Maps Foundation · ODbL", ha="center", fontsize=9)
plt.show()
maps(networks, "Boulder roads · retained native segments")
These are mapped features, not evidence of traffic or road condition. Overture incorporates OSM; apparent agreement is not independent validation.
2. Make the analysis boundary explicit¶
Which parts of the retained features fall inside our area?
def within_area(boundary):
# Project-owned verb: preserve feature attributes and clip geometry explicitly.
def operation(frame):
return frame.clip(boundary).copy()
return operation
clipped = {name: (pipe(frame) | within_area(area)).unwrap()
for name, frame in networks.items()}
assert all(frame.geometry.covered_by(area).all() for frame in clipped.values())
maps(clipped, "Same boundary · explicitly clipped geometries")
Clipping is a deliberate analytical change, not an invisible loader operation. It cannot recover a crossing OSM way omitted by the provider's node-in-bbox query.
3. Measure length without inventing a class crosswalk¶
How is retained road length distributed among each provider's own classes?
def length_by_class(crs):
def operation(frame):
# Geographic degrees are not distances. UTM 13N is local to this Boulder example.
projected = frame.to_crs(crs)
return (projected.assign(length_km=projected.length / 1000)
.groupby("source_classification").length_km.sum().sort_values())
return operation
lengths = {name: (pipe(frame) | within_area(area) | length_by_class("EPSG:32613")).unwrap()
for name, frame in networks.items()}
fig, axes = plt.subplots(1, 2, figsize=(10, 4), layout="constrained")
for (name, result), ax in zip(lengths.items(), axes):
result.plot.barh(ax=ax, color="#236d81")
ax.set(title=f"{name} · native classes", xlabel="Mapped length inside area (km)", ylabel="")
plt.show()
Native classes and segmentation remain different. These lengths describe the retained mapped sample, not completeness, accessibility, routing connectivity, or a source-quality ranking.
Figure¶
Each of the three analytical steps displays its own result.
What the figure tells us¶
Native classes and segmentation remain different. These lengths describe the retained mapped sample, not completeness, accessibility, routing connectivity, or a source-quality ranking.
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 | Overture Maps Foundation and OpenStreetMap contributors |
| Product | Native mapped road segments, Boulder |
| Dates | Overture release 2026-08-19.0; OSM retrieval recorded in native metadata |
| Fixture | tests/fixtures/real_data/source_lessons/roads_overture.geojson |
| Provenance record | tests/fixtures/real_data/source_lessons/roads.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/roads_local_network.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.
See also¶
roads ·