From 6116f00f8fb241cbd1b4e46cd3d97be8769b55db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9rgio=20Costa?= <noreply@example.com>
Date: Fri, 2 Oct 2026 01:43:13 +0000
Subject: [PATCH 1/2] Add the `osm` source: OpenStreetMap ways through Overpass

type = "osm" reads ways (roads, rivers, ...) by Overpass statements, clipped to
the grid plus a margin, and registers them as a vector source with checksum and
provenance, ready for the `distance` and `count` operators.

- The answer is cached by request (statements, bbox, date): a pipeline reads the
  same bytes, hence the same spec_hash, on every run and workspace.
- `date` asks for an older snapshot ([date:"..."]); with `years`,
  date = "{year}-07-01" gives one per year.
- Several public Overpass servers are tried in turn with retries; HTTP 406 stops
  at once and points to OSM_CONTACT; an empty answer is an error and not cached.
- Provenance: request, server, time retrieved, timestamp_osm_base, features,
  checksum, ODbL attribution.
- Tests run against a local fake Overpass server (no network); guide
  docs/guides/osm.md.

Ports 02b_fetch_osm.py of lab15-reconstruction. Not exercised against the real
Overpass API (unreachable from the environment this was written in).

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                  |  10 ++
 README.md                     |   1 +
 disscube/pipeline/runner.py   |  16 +-
 disscube/pipeline/schema.py   |  27 +++-
 disscube/sources/osm.py       | 247 ++++++++++++++++++++++++++++++
 docs/guides/osm.md            |  79 ++++++++++
 docs/guides/pipeline_files.md |   3 +-
 mkdocs.yml                    |   1 +
 tests/test_osm.py             | 276 ++++++++++++++++++++++++++++++++++
 9 files changed, 653 insertions(+), 7 deletions(-)
 create mode 100644 disscube/sources/osm.py
 create mode 100644 docs/guides/osm.md
 create mode 100644 tests/test_osm.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58a60cb..dc9ce15 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+- **`type = "osm"` source**: ways from OpenStreetMap through the Overpass API, clipped to the grid
+  plus a `margin`, registered as a vector source with checksum and provenance (request, server,
+  time retrieved, `timestamp_osm_base`, ODbL). The answer is cached by request, so a pipeline reads
+  the same data (and gets the same `spec_hash`) on every run; `date` asks for an older snapshot
+  (`date = "{year}-07-01"` with `years`); several public servers are tried with retries; HTTP 406
+  points to `OSM_CONTACT`. Replaces the OpenStreetMap download script of the Lab15 reconstruction.
+
+## [0.4.0] - 2026-09-30
+
 ### Changed
 - **DisSModel is now optional.** DisSCube no longer depends on it: install
   `disscube[dissmodel]` only to hand a cube to a model. `CubeClient.to_lucc_data()`
diff --git a/README.md b/README.md
index a0255a5..8643c84 100644
--- a/README.md
+++ b/README.md
@@ -246,6 +246,7 @@ disscube/
 │   ├── _categorical.py  legends (.qml/.json/.csv) and reclassification
 │   ├── mapbiomas.py  MapBiomas annual land-cover maps (Collection 11, 10 m series)
 │   ├── prodes.py     PRODES deforestation (download + cache, legend from the .qml)
+│   ├── osm.py        OpenStreetMap ways via Overpass (cached, with date snapshots)
 │   └── classified.py any classified map with its legend, e.g. from SITS
 └── utils.py          Checksums (sha256_file) and BDC tile importer (import_bdc_grids)
 ```
diff --git a/disscube/pipeline/runner.py b/disscube/pipeline/runner.py
index 339173f..c576101 100644
--- a/disscube/pipeline/runner.py
+++ b/disscube/pipeline/runner.py
@@ -23,6 +23,7 @@
     FileSource,
     GridConfig,
     MapbiomasSource,
+    OsmSource,
     PipelineConfig,
     ProdesSource,
     UnionSource,
@@ -56,7 +57,7 @@ def base_dir(self) -> Path:
 @dataclass
 class PlannedSource:
     id: str
-    config: FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | UnionSource
+    config: FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | UnionSource
     year: int | None = None
 
 
@@ -412,7 +413,7 @@ def _register_source(
 ):
     c = s.config
     # ── Validação defensiva: apenas fontes dinâmicas em janela na nuvem precisam de bbox_geo ──
-    if isinstance(c, (BdcSource, MapbiomasSource, ProdesSource, ClassifiedSource)):
+    if isinstance(c, (BdcSource, MapbiomasSource, ProdesSource, ClassifiedSource, OsmSource)):
         if bbox_geo is None:
             raise PipelineError(
                 f"source {c.id!r} ({c.type}): windowed cloud sources require an"
@@ -448,6 +449,17 @@ def _register_source(
                 files=files,
                 name=c.name,
             )
+        if isinstance(c, OsmSource):
+            from disscube.sources import osm
+
+            try:
+                return osm.register_osm_source(
+                    cube, c.id, c.query, bbox_geo, raw, margin=c.margin, date=c.date,
+                    endpoints=c.endpoints, cache=(base / c.cache) if c.cache else None,
+                    time=c.time, name=c.name,
+                )
+            except (osm.OsmError, ValueError) as exc:
+                raise PipelineError(f"source {c.id!r}: {exc}") from None
         if isinstance(c, ClassifiedSource):
             from disscube.sources.classified import register_classified_map
 
diff --git a/disscube/pipeline/schema.py b/disscube/pipeline/schema.py
index 2ce5252..bd1c577 100644
--- a/disscube/pipeline/schema.py
+++ b/disscube/pipeline/schema.py
@@ -131,6 +131,25 @@ class ClassifiedSource(_SourceBase):
     producer: str | None = None
 
 
+class OsmSource(_SourceBase):
+    """Ways from OpenStreetMap (Overpass), clipped to the grid plus ``margin``.
+
+    ``query`` holds Overpass statements without the bounding box (``way["highway"]["ref"~"BR-163"]``),
+    separated by ``;`` or given as a list. ``margin`` is in degrees: a ``distance`` needs
+    features beyond the grid, a ``count`` does not. ``date`` asks for the map as it was on
+    that day; with ``years``, ``date = "{year}-07-01"`` gives one snapshot per year. The answer
+    is cached, so a pipeline reads the same data until the cache is cleared. Set ``OSM_CONTACT``.
+    """
+
+    type: Literal["osm"]
+    query: str | list[str]
+    margin: float = Field(default=0.0, ge=0)
+    date: str | None = None
+    endpoints: list[str] | None = None
+    cache: str | None = None
+    time: int | None = None
+
+
 class UnionSource(_SourceBase):
     """The features of several vector sources of this file as one source."""
 
@@ -139,7 +158,7 @@ class UnionSource(_SourceBase):
 
 
 Source = Annotated[
-    FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | UnionSource,
+    FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | UnionSource,
     Field(discriminator="type"),
 ]
 
@@ -174,7 +193,7 @@ class PipelineConfig(_Strict):
     grid: GridConfig | None = None
     extent: list[float] | None = Field(default=None, min_length=4, max_length=4)
     """``[min_lon, min_lat, max_lon, max_lat]`` (WGS84) a sources-only file reads
-    windowed sources over (classified maps, BDC, MapBiomas, PRODES)."""
+    windowed sources over (classified maps, BDC, MapBiomas, PRODES, OpenStreetMap)."""
     sources_from_catalog: bool = False
     """Let [[derive]] blocks use sources this file does not declare, registered
     in the workspace's catalog by another pipeline file. Off by default, so a
@@ -190,14 +209,14 @@ def _grid_or_extent(self):
 
         # Apenas fontes dinâmicas em janela na nuvem precisam de extent.
         # Fontes 'file' e 'union' já têm extensão definida por seus arquivos locais.
-        windowed_types = {"bdc", "mapbiomas", "prodes", "classified"}
+        windowed_types = {"bdc", "mapbiomas", "prodes", "classified", "osm"}
         needs_extent = any(
             getattr(s, "type", None) in windowed_types for s in self.source
         )
 
         if self.grid is None and self.extent is None and needs_extent:
             raise ValueError(
-                "a sources-only file with windowed sources (BDC, MapBiomas, PRODES)"
+                "a sources-only file with windowed sources (BDC, MapBiomas, PRODES, OpenStreetMap)"
                 " needs `extent` (or a [grid])"
             )
         return self
diff --git a/disscube/sources/osm.py b/disscube/sources/osm.py
new file mode 100644
index 0000000..6dd6503
--- /dev/null
+++ b/disscube/sources/osm.py
@@ -0,0 +1,247 @@
+"""Vector features from OpenStreetMap, through the Overpass API.
+
+A pipeline asks for ways by Overpass statements and gets them clipped to the grid
+(plus a margin) as a vector source with a checksum and provenance::
+
+    [[source]]
+    id     = "roads"
+    type   = "osm"
+    query  = 'way["highway"]["ref"~"BR-163"]'
+    margin = 0.7                       # degrees beyond the grid: distances need far features
+
+The statements are Overpass QL without the bounding box, which is added to each one
+(``way[...](south,west,north,east)``), so ``query`` can hold several, separated by ``;``
+or given as a list. Only ways are read; they come back as ``LineString`` in EPSG:4326
+(a closed way, such as a reservoir, stays as its outline).
+
+**Reproducibility.** OpenStreetMap is a moving target. The answer is saved once in a
+cache (``$DISSCUBE_CACHE/osm``, or ``~/.cache/disscube/osm``) keyed by what was asked
+(statements, bounding box, date), and every later run reuses it: the same pipeline gives
+the same bytes, hence the same ``spec_hash``, until the cache is deleted. Its checksum and
+the time it was retrieved go to the provenance. ``date`` asks Overpass for the map as it
+was on that day (``[date:"..."]``), which keeps later roads out of an earlier period.
+
+**Etiquette.** Overpass answers HTTP 406 to clients that do not identify themselves, so set
+``OSM_CONTACT`` (an e-mail or URL) in the environment; it goes in the ``User-Agent``. Several
+public servers are tried in turn, with retries. Data © OpenStreetMap contributors (ODbL).
+"""
+
+from __future__ import annotations
+
+import hashlib
+import json
+import logging
+import os
+import time
+import urllib.error
+import urllib.parse
+import urllib.request
+from collections.abc import Callable, Sequence
+from datetime import UTC, datetime
+from pathlib import Path
+from typing import Any
+
+from disscube.utils import sha256_file
+
+log = logging.getLogger("disscube.sources.osm")
+
+#: Public Overpass servers, tried in turn (the main one often answers 504 when busy).
+DEFAULT_ENDPOINTS: tuple[str, ...] = (
+    "https://overpass-api.de/api/interpreter",
+    "https://overpass.kumi.systems/api/interpreter",
+    "https://overpass.private.coffee/api/interpreter",
+    "https://maps.mail.ru/osm/tools/overpass/api/interpreter",
+)
+#: Passes over the endpoints before giving up.
+ROUNDS = 3
+#: Seconds the Overpass query may run server-side.
+QUERY_TIMEOUT = 180
+LICENSE = "ODbL"
+ATTRIBUTION = "© OpenStreetMap contributors"
+#: Tags kept as feature properties.
+KEPT_TAGS = ("name", "ref", "highway", "waterway", "natural", "railway", "boundary")
+
+# Seconds to wait between attempts; replaced in tests.
+_sleep: Callable[[float], None] = time.sleep
+
+
+class OsmError(RuntimeError):
+    """The OpenStreetMap data could not be obtained."""
+
+
+def default_cache_dir() -> Path:
+    """``$DISSCUBE_CACHE/osm``, or ``~/.cache/disscube/osm``."""
+    root = os.environ.get("DISSCUBE_CACHE") or Path.home() / ".cache" / "disscube"
+    return Path(root) / "osm"
+
+
+def user_agent() -> str:
+    contact = os.environ.get("OSM_CONTACT", "").strip()
+    if not contact:
+        log.warning("OSM_CONTACT is not set: public Overpass servers reject clients without a contact (HTTP 406)")
+    return f"disscube ({contact or 'no contact given'})"
+
+
+def statements_of(query: str | Sequence[str]) -> list[str]:
+    """The Overpass statements of ``query`` (a string split at ``;``, or a list), stripped."""
+    parts = query.split(";") if isinstance(query, str) else list(query)
+    statements = [p.strip().rstrip(";").strip() for p in parts if p and p.strip()]
+    if not statements:
+        raise ValueError("an osm source needs a non-empty 'query'")
+    for st in statements:
+        if st.split("[", 1)[0].strip() != "way":
+            raise ValueError(f"an osm statement must filter ways, e.g. way[\"highway\"]: {st!r}")
+    return statements
+
+
+def padded(bbox_geo: Sequence[float], margin: float) -> list[float]:
+    """``[west, south, east, north]`` widened by ``margin`` degrees, clipped to the globe."""
+    w, s, e, n = (float(v) for v in bbox_geo)
+    return [max(w - margin, -180.0), max(s - margin, -90.0), min(e + margin, 180.0), min(n + margin, 90.0)]
+
+
+def build_query(statements: Sequence[str], bbox: Sequence[float], date: str | None = None,
+                timeout: int = QUERY_TIMEOUT) -> str:
+    """The Overpass QL text asking for ``statements`` inside ``bbox`` (``[w, s, e, n]``)."""
+    w, s, e, n = bbox
+    header = f"[out:json][timeout:{timeout}]"
+    if date:
+        header += f'[date:"{_overpass_date(date)}"]'
+    body = "".join(f"{st}({s},{w},{n},{e});" for st in statements)
+    return f"{header};({body});out geom;"
+
+
+def _overpass_date(date: str) -> str:
+    """``2019-01-01`` or a full ISO timestamp, as Overpass wants it (``YYYY-MM-DDThh:mm:ssZ``)."""
+    try:
+        parsed = datetime.fromisoformat(date)
+    except ValueError:
+        raise ValueError(f"osm 'date' must be ISO 8601 (e.g. 2019-01-01): {date!r}") from None
+    if parsed.tzinfo is not None:
+        parsed = parsed.astimezone(UTC)
+    return parsed.strftime("%Y-%m-%dT%H:%M:%SZ")
+
+
+def _require_http(url: str) -> str:
+    if urllib.parse.urlparse(url).scheme not in ("http", "https"):
+        raise ValueError(f"Overpass endpoint must be http(s): {url!r}")
+    return url
+
+
+def overpass(query: str, endpoints: Sequence[str] = DEFAULT_ENDPOINTS, timeout: int = QUERY_TIMEOUT,
+             rounds: int = ROUNDS) -> tuple[dict[str, Any], str]:
+    """Run ``query`` on the first endpoint that answers; returns ``(result, endpoint)``.
+
+    Failures (HTTP 5xx, timeouts, a server-side ``remark`` with no elements) move on to the
+    next endpoint; after a full pass it waits and tries again, ``rounds`` times. HTTP 406
+    stops at once: it means the client is not identified (set ``OSM_CONTACT``).
+    """
+    data = urllib.parse.urlencode({"data": query}).encode()
+    agent = user_agent()
+    errors: list[str] = []
+    for rnd in range(rounds):
+        for url in endpoints:
+            req = urllib.request.Request(_require_http(url), data=data, headers={"User-Agent": agent})
+            host = urllib.parse.urlparse(url).netloc
+            try:
+                with urllib.request.urlopen(req, timeout=timeout + 30) as resp:  # nosec B310
+                    result = json.load(resp)
+                if "remark" in result and not result.get("elements"):
+                    raise OsmError(str(result["remark"]))  # server-side timeout / out of memory
+                return result, url
+            except urllib.error.HTTPError as exc:
+                if exc.code == 406:
+                    raise OsmError("HTTP 406 from Overpass: set OSM_CONTACT (an e-mail or URL) so it can identify you") from exc
+                errors.append(f"{host}: HTTP {exc.code}")
+            except Exception as exc:  # noqa: BLE001 — any failure of one server: try the next
+                errors.append(f"{host}: {type(exc).__name__}: {exc}")
+            log.warning("%s", errors[-1])
+            _sleep(3)
+        if rnd < rounds - 1:
+            _sleep(10 * (rnd + 1))
+    raise OsmError("Overpass did not answer: " + "; ".join(errors[-len(endpoints):]))
+
+
+def to_geojson(result: dict[str, Any]) -> dict[str, Any]:
+    """Overpass ``out geom`` ways as a GeoJSON FeatureCollection of ``LineString`` (EPSG:4326)."""
+    features = []
+    for el in result.get("elements", []):
+        geom = el.get("geometry")
+        if el.get("type") != "way" or not geom or len(geom) < 2:
+            continue
+        tags = el.get("tags", {})
+        features.append({
+            "type": "Feature",
+            "properties": {"osm_id": el["id"], **{k: tags[k] for k in KEPT_TAGS if k in tags}},
+            "geometry": {"type": "LineString", "coordinates": [[p["lon"], p["lat"]] for p in geom]},
+        })
+    return {"type": "FeatureCollection", "crs": {"type": "name", "properties": {"name": "EPSG:4326"}},
+            "features": features}
+
+
+def cache_key(statements: Sequence[str], bbox: Sequence[float], date: str | None) -> str:
+    """What was asked, as a short hash: the same request finds the same cached answer."""
+    request = {"statements": list(statements), "bbox": [round(float(v), 6) for v in bbox],
+               "date": _overpass_date(date) if date else None}
+    return hashlib.sha256(json.dumps(request, sort_keys=True).encode()).hexdigest()[:16]
+
+
+def fetch(query: str | Sequence[str], bbox_geo: Sequence[float], margin: float = 0.0, date: str | None = None,
+          endpoints: Sequence[str] | None = None, cache: str | Path | None = None) -> tuple[Path, dict[str, Any]]:
+    """The GeoJSON file for ``query`` over ``bbox_geo`` (+ ``margin``) and what is known about it.
+
+    Returns ``(path, info)``. A request already in the cache is not repeated. An empty answer
+    raises :class:`OsmError` and is not cached: it almost always means a wrong query or area.
+    """
+    statements = statements_of(query)
+    bbox = padded(bbox_geo, margin)
+    key = cache_key(statements, bbox, date)
+    folder = Path(cache) if cache else default_cache_dir()
+    path, meta_path = folder / f"osm-{key}.geojson", folder / f"osm-{key}.json"
+    if path.exists() and meta_path.exists():
+        log.info("OpenStreetMap: reusing %s", path)
+        return path, json.loads(meta_path.read_text(encoding="utf-8"))
+
+    result, endpoint = overpass(build_query(statements, bbox, date), tuple(endpoints or DEFAULT_ENDPOINTS))
+    collection = to_geojson(result)
+    if not collection["features"]:
+        raise OsmError(f"OpenStreetMap returned no ways for {statements} in {bbox}"
+                       + (f" as of {date}" if date else "") + ": check the query and the area")
+
+    folder.mkdir(parents=True, exist_ok=True)
+    tmp = path.with_suffix(".tmp")
+    tmp.write_text(json.dumps(collection), encoding="utf-8")
+    tmp.replace(path)
+    info = {"statements": statements, "bbox": bbox, "margin": margin, "date": date, "endpoint": endpoint,
+            "retrieved": datetime.now(UTC).isoformat(timespec="seconds"), "features": len(collection["features"]),
+            "checksum": sha256_file(path), "license": LICENSE, "attribution": ATTRIBUTION,
+            "osm_base": result.get("osm3s", {}).get("timestamp_osm_base")}
+    meta_path.write_text(json.dumps(info, indent=2, ensure_ascii=False), encoding="utf-8")
+    return path, info
+
+
+def register_osm_source(cube, source_id: str, query: str | Sequence[str], bbox_geo: Sequence[float],
+                        out_dir: str | Path, margin: float = 0.0, date: str | None = None,
+                        endpoints: Sequence[str] | None = None, cache: str | Path | None = None,
+                        time: int | None = None, name: str | None = None):
+    """Fetch (or reuse) the OpenStreetMap ways and register them as a vector source.
+
+    The GeoJSON is copied to ``out_dir/<source_id>.geojson`` next to a
+    ``<source_id>.provenance.json`` with the request, the server, the time retrieved, the
+    ``timestamp_osm_base`` of the data and the checksum.
+    """
+    from disscube.models import SpatialSource
+
+    cached, info = fetch(query, bbox_geo, margin=margin, date=date, endpoints=endpoints, cache=cache)
+    out = Path(out_dir)
+    out.mkdir(parents=True, exist_ok=True)
+    path = out / f"{source_id}.geojson"
+    path.write_bytes(cached.read_bytes())
+    prov_path = out / f"{source_id}.provenance.json"
+    provenance = {"type": "osm", "file": str(path), "crs": "EPSG:4326", "format": "vector", **info}
+    prov_path.write_text(json.dumps(provenance, indent=2, ensure_ascii=False), encoding="utf-8")
+    src = SpatialSource(id=source_id, name=name or source_id, format="vector", asset_url=str(path),
+                        crs="EPSG:4326", time=time, checksum=info["checksum"],
+                        tags=["osm", f"provenance:{prov_path}"])
+    cube.register_spatial_source(src)
+    return src
diff --git a/docs/guides/osm.md b/docs/guides/osm.md
new file mode 100644
index 0000000..1fbe6dd
--- /dev/null
+++ b/docs/guides/osm.md
@@ -0,0 +1,79 @@
+# OpenStreetMap
+
+`type = "osm"` takes roads, rivers and other **ways** from OpenStreetMap (Overpass API) and
+registers them as a vector source, clipped to the grid plus a margin, with a checksum and
+provenance. It is the usual input of the `distance` and `count` operators.
+
+```toml
+[grid]
+name = "area"
+bbox = [-54.842, -3.587, -54.459, -3.168]
+resolution = 500
+
+[[source]]
+id     = "br163"
+type   = "osm"
+query  = 'way["highway"]["ref"~"BR-163"]'
+margin = 0.7                              # degrees; distances need features beyond the grid
+
+[[derive]]
+target   = "dist_br"
+source   = "br163"
+operator = "distance"
+```
+
+```bash
+export OSM_CONTACT="you@example.org"      # Overpass answers HTTP 406 to anonymous clients
+disscube run area.toml --workspace outputs/area
+```
+
+## Fields
+
+| Field | Meaning |
+|---|---|
+| `query` | Overpass statements **without** the bounding box, as a string separated by `;` or a list. Only `way[...]` filters are read. |
+| `margin` | Degrees added around the grid (default `0`). A `distance` to a road outside the grid needs it; a `count` of roads inside does not. |
+| `date` | The map as it was on that day (`2019-01-01` or a full ISO timestamp), through Overpass's `[date:"…"]`. |
+| `endpoints` | Overpass servers to try, in order (default: four public ones). Use it for your own instance. |
+| `cache` | Cache folder, relative to the pipeline file (default `$DISSCUBE_CACHE/osm` or `~/.cache/disscube/osm`). |
+| `time` | Time of the source; set from `years` when you use `{year}`. |
+
+Features come back as `LineString` in EPSG:4326 with a few tags (`name`, `ref`, `highway`,
+`waterway`, `natural`, `railway`, `boundary`). A closed way, such as a reservoir, stays as its outline.
+
+## Reproducibility
+
+OpenStreetMap changes every minute, so DisSCube fetches **once** and keeps the answer in the
+cache, keyed by what was asked (statements, bounding box, date). Every later run, in any
+workspace, reuses it: the same bytes, the same checksum, the same `spec_hash`. To get fresh
+data, delete the cache file. `raw/<id>.provenance.json` records the request, the server, the
+time it was retrieved, the `timestamp_osm_base` of the data, the number of features, the
+checksum and the licence (ODbL, © OpenStreetMap contributors).
+
+An empty answer is an error and is not cached: it almost always means a wrong query or area.
+
+## Maps of the past
+
+OpenStreetMap is today's map by default, so roads opened after the period you study leak into
+the drivers. With `date` the query asks for an older state; with `years` you get one snapshot
+per year:
+
+```toml
+[[source]]
+id     = "roads_{year}"
+type   = "osm"
+query  = 'way["highway"]'
+date   = "{year}-07-01"
+years  = [2010, 2015, 2020]
+margin = 0.1
+```
+
+OpenStreetMap is sparse in many regions before the mid-2000s, so an old snapshot is a
+lower bound, not the road network of the time.
+
+## Etiquette
+
+The public Overpass servers are shared. Keep queries small, set `OSM_CONTACT`, and rely on the
+cache instead of running a pipeline repeatedly with a cleared one. Failures (HTTP 5xx,
+timeouts) move on to the next server and retry; HTTP 406 stops at once with a message about
+`OSM_CONTACT`.
diff --git a/docs/guides/pipeline_files.md b/docs/guides/pipeline_files.md
index c475eea..5333701 100644
--- a/docs/guides/pipeline_files.md
+++ b/docs/guides/pipeline_files.md
@@ -71,7 +71,7 @@ resolution = 300
 
 [[source]]                    # one block per source; see the types below
 id = "…"
-type = "file" | "bdc" | "mapbiomas" | "prodes" | "classified" | "union"
+type = "file" | "bdc" | "mapbiomas" | "prodes" | "classified" | "osm" | "union"
 
 [[derive]]                    # one block per derived variable
 target = "urban_pct"
@@ -97,6 +97,7 @@ relative to the pipeline file, so a folder with its TOML and data is portable.
 | `mapbiomas` | `year`, `collection` (11), `resolution` (30), `url` | `disscube.sources.mapbiomas` |
 | `prodes` | `year`, `url`, `cache` | `disscube.sources.prodes` |
 | `classified` | `path`, `legend` (table or `.qml`/`.json`/`.csv` file), `time`, `nodata`, `producer` | `disscube.sources.classified` |
+| `osm` | `query`, `margin`, `date`, `endpoints`, `cache`, `time` (see [OpenStreetMap](osm.md)) | `disscube.sources.osm` |
 | `union` | `of` (ids of vector sources declared before it) | their features in one GeoPackage under `raw/` |
 
 Every source block also accepts `name` and `years`.
diff --git a/mkdocs.yml b/mkdocs.yml
index 612d65a..662be32 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -43,4 +43,5 @@ nav:
       - BDC Integration: guides/bdc.md
       - MapBiomas: guides/mapbiomas.md
       - "Land-cover maps: PRODES and SITS": guides/landcover.md
+      - "OpenStreetMap": guides/osm.md
       - Custom Grids: guides/grids.md
diff --git a/tests/test_osm.py b/tests/test_osm.py
new file mode 100644
index 0000000..9b27957
--- /dev/null
+++ b/tests/test_osm.py
@@ -0,0 +1,276 @@
+"""OpenStreetMap source: query building, caching, failover and the pipeline type.
+
+A local HTTP server plays Overpass, so nothing here touches the network.
+"""
+
+from __future__ import annotations
+
+import json
+import threading
+import urllib.parse
+from http.server import BaseHTTPRequestHandler, HTTPServer
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from disscube.pipeline import PipelineError, load, plan, run
+from disscube.sources import osm
+
+# A straight road across the middle of the 3 km grid below (EPSG:31983 ≈ lon -46.3, lat -23.6 at the origin).
+BBOX = [-46.30, -23.60, -46.27, -23.57]
+ROAD = {"type": "way", "id": 11, "tags": {"highway": "trunk", "ref": "BR-1", "surface": "asphalt"},
+        "geometry": [{"lat": -23.585, "lon": -46.31}, {"lat": -23.585, "lon": -46.26}]}
+
+
+class _Overpass:
+    """A fake Overpass server. ``script`` holds one status per request; the last one repeats."""
+
+    def __init__(self, elements=None, script=(200,), remark=None):
+        self.requests: list[dict] = []
+        self.script = list(script)
+        outer = self
+
+        class Handler(BaseHTTPRequestHandler):
+            def do_POST(self):
+                body = self.rfile.read(int(self.headers["Content-Length"]))
+                outer.requests.append({"query": urllib.parse.parse_qs(body.decode())["data"][0],
+                                       "agent": self.headers.get("User-Agent")})
+                status = outer.script[min(len(outer.requests) - 1, len(outer.script) - 1)]
+                payload = {"osm3s": {"timestamp_osm_base": "2026-09-30T00:00:00Z"},
+                           "elements": [ROAD] if elements is None else elements}
+                if remark:
+                    payload = {"remark": remark, "elements": []}
+                self.send_response(status)
+                self.send_header("Content-Type", "application/json")
+                self.end_headers()
+                self.wfile.write(json.dumps(payload).encode())
+
+            def log_message(self, *args):
+                pass
+
+        self.server = HTTPServer(("127.0.0.1", 0), Handler)
+        self.url = f"http://127.0.0.1:{self.server.server_port}/api/interpreter"
+        threading.Thread(target=lambda: self.server.serve_forever(poll_interval=0.01), daemon=True).start()
+
+    def close(self):
+        self.server.shutdown()
+        self.server.server_close()
+
+
+@pytest.fixture
+def overpass():
+    servers = []
+
+    def make(**kw):
+        servers.append(_Overpass(**kw))
+        return servers[-1]
+
+    yield make
+    for s in servers:
+        s.close()
+
+
+@pytest.fixture(autouse=True)
+def _env(tmp_path, monkeypatch):
+    monkeypatch.setenv("DISSCUBE_CACHE", str(tmp_path / "cache"))
+    monkeypatch.setenv("OSM_CONTACT", "tests@example.org")
+    monkeypatch.setattr(osm, "_sleep", lambda s: None)
+
+
+# --- query building -------------------------------------------------------------------------
+
+def test_statements_are_split_cleaned_and_checked():
+    assert osm.statements_of('way["highway"]; way["waterway"] ;') == ['way["highway"]', 'way["waterway"]']
+    assert osm.statements_of(['way["a"]', 'way["b"]']) == ['way["a"]', 'way["b"]']
+    for bad in ("", "  ;  ", 'node["place"]', 'rel["route"]'):
+        with pytest.raises(ValueError):
+            osm.statements_of(bad)
+
+
+def test_bbox_is_added_to_each_statement_and_padded():
+    assert osm.padded([-46.3, -23.6, -46.27, -23.57], 0.1) == pytest.approx([-46.4, -23.7, -46.17, -23.47])
+    assert osm.padded([-179.9, -89.95, 179.9, 89.95], 1) == [-180.0, -90.0, 180.0, 90.0]
+    q = osm.build_query(['way["a"]', 'way["b"]'], [1, 2, 3, 4])
+    assert 'way["a"](2,1,4,3);way["b"](2,1,4,3);' in q and q.endswith("out geom;")
+    assert "[date:" not in q
+
+
+def test_date_gives_a_snapshot_query_and_must_be_iso():
+    assert '[date:"2019-01-01T00:00:00Z"]' in osm.build_query(['way["a"]'], [1, 2, 3, 4], "2019-01-01")
+    assert '[date:"2019-07-01T12:30:00Z"]' in osm.build_query(['way["a"]'], [1, 2, 3, 4], "2019-07-01T12:30:00Z")
+    with pytest.raises(ValueError, match="ISO 8601"):
+        osm.build_query(['way["a"]'], [1, 2, 3, 4], "last year")
+
+
+def test_geojson_keeps_lines_and_a_few_tags():
+    fc = osm.to_geojson({"elements": [ROAD, {"type": "node", "id": 1, "lat": 0, "lon": 0},
+                                      {"type": "way", "id": 2, "geometry": [{"lat": 0, "lon": 0}]}]})
+    assert [f["properties"]["osm_id"] for f in fc["features"]] == [11]
+    props = fc["features"][0]["properties"]
+    assert props == {"osm_id": 11, "highway": "trunk", "ref": "BR-1"}  # `surface` is not kept
+    assert fc["features"][0]["geometry"]["coordinates"][0] == [-46.31, -23.585]
+
+
+# --- fetching --------------------------------------------------------------------------------
+
+def test_fetch_caches_the_answer_and_identifies_itself(overpass):
+    srv = overpass()
+    path, info = osm.fetch('way["highway"]', BBOX, margin=0.05, endpoints=[srv.url])
+    assert len(srv.requests) == 1 and srv.requests[0]["agent"] == "disscube (tests@example.org)"
+    assert info["features"] == 1 and info["license"] == "ODbL" and info["osm_base"] == "2026-09-30T00:00:00Z"
+    assert info["bbox"] == pytest.approx([-46.35, -23.65, -46.22, -23.52])
+
+    again, info2 = osm.fetch('way["highway"]', BBOX, margin=0.05, endpoints=[srv.url])
+    assert again == path and info2 == info and len(srv.requests) == 1  # no second request
+
+    osm.fetch('way["highway"]', BBOX, margin=0.06, endpoints=[srv.url])  # another margin is another request
+    osm.fetch('way["highway"]', BBOX, margin=0.05, date="2015-01-01", endpoints=[srv.url])
+    assert len(srv.requests) == 3
+
+
+def test_failover_to_the_next_server_and_retry_rounds(overpass):
+    broken, ok = overpass(script=(504,)), overpass()
+    osm.fetch('way["highway"]', BBOX, endpoints=[broken.url, ok.url])
+    assert len(broken.requests) == 1 and len(ok.requests) == 1
+
+    flaky = overpass(script=(503, 503, 200))
+    osm.fetch('way["waterway"]', BBOX, endpoints=[flaky.url])  # fails twice, then answers
+    assert len(flaky.requests) == 3
+
+
+def test_gives_up_after_all_rounds_and_says_why(overpass):
+    down = overpass(script=(504,))
+    with pytest.raises(osm.OsmError, match="did not answer.*HTTP 504"):
+        osm.fetch('way["highway"]', BBOX, endpoints=[down.url])
+    assert len(down.requests) == osm.ROUNDS
+
+
+def test_406_stops_at_once_and_points_to_osm_contact(overpass):
+    srv = overpass(script=(406,))
+    with pytest.raises(osm.OsmError, match="OSM_CONTACT"):
+        osm.fetch('way["highway"]', BBOX, endpoints=[srv.url, srv.url])
+    assert len(srv.requests) == 1
+
+
+def test_server_side_remark_counts_as_a_failure(overpass):
+    srv = overpass(remark="runtime error: Query timed out", script=(200,))
+    with pytest.raises(osm.OsmError, match="Query timed out"):
+        osm.fetch('way["highway"]', BBOX, endpoints=[srv.url])
+
+
+def test_empty_answer_is_an_error_and_is_not_cached(overpass):
+    srv = overpass(elements=[])
+    for _ in range(2):
+        with pytest.raises(osm.OsmError, match="no ways"):
+            osm.fetch('way["highway"]', BBOX, endpoints=[srv.url])
+    assert len(srv.requests) == 2 and not list(osm.default_cache_dir().glob("*.geojson"))
+
+
+def test_endpoint_must_be_http():
+    with pytest.raises(ValueError, match="http"):
+        osm.fetch('way["highway"]', BBOX, endpoints=["file:///etc/passwd"])
+
+
+# --- in a pipeline ---------------------------------------------------------------------------
+
+def _pipeline(tmp_path, srv, extra="", source_extra=""):
+    path = tmp_path / "osm.toml"
+    path.write_text(f"""
+schema = 1
+[grid]
+name = "g"
+crs = "EPSG:31983"
+bbox = [333000, 7388000, 336000, 7391000]
+resolution = 300
+{extra}
+[[source]]
+id = "road"
+type = "osm"
+query = 'way["highway"]["ref"~"BR-1"]'
+endpoints = ["{srv.url}"]
+{source_extra}
+[[derive]]
+target = "dist_road"
+source = "road"
+operator = "distance"
+""", encoding="utf-8")
+    return path
+
+
+def test_pipeline_fetches_registers_and_derives_a_distance(tmp_path, overpass):
+    srv = overpass()
+    report = run(_pipeline(tmp_path, srv), workspace=tmp_path / "ws")
+    assert report.derived and len(srv.requests) == 1
+    prov = json.loads((tmp_path / "ws" / "raw" / "road.provenance.json").read_text())
+    assert prov["type"] == "osm" and prov["format"] == "vector" and prov["features"] == 1
+    assert prov["checksum"].startswith("sha256:") and prov["statements"] == ['way["highway"]["ref"~"BR-1"]']
+    assert prov["attribution"] == "© OpenStreetMap contributors" and Path(prov["pipeline"]["file"]).name == "osm.toml"
+    assert (tmp_path / "ws" / "raw" / "road.geojson").exists()
+
+
+def test_same_pipeline_in_another_workspace_reuses_the_cache_and_the_spec_hash(tmp_path, overpass):
+    from disscube import CubeClient
+
+    srv = overpass()
+    path = _pipeline(tmp_path, srv)
+    run(path, workspace=tmp_path / "a")
+    run(path, workspace=tmp_path / "b")
+    assert len(srv.requests) == 1
+
+    def spec_hash(ws):
+        cube = CubeClient(str(ws / "catalog.db"), str(ws / "store"))
+        return [d.spec_hash for d in cube.catalog.search_derived_variables() if d.name == "dist_road"]
+
+    assert spec_hash(tmp_path / "a") == spec_hash(tmp_path / "b") != []
+
+
+def test_years_give_one_dated_snapshot_each(tmp_path, overpass):
+    srv = overpass()
+    path = tmp_path / "years.toml"
+    path.write_text(f"""
+schema = 1
+extent = [-46.30, -23.60, -46.27, -23.57]
+[[source]]
+id = "road_{{year}}"
+type = "osm"
+query = 'way["highway"]'
+date = "{{year}}-07-01"
+endpoints = ["{srv.url}"]
+years = [2010, 2020]
+""", encoding="utf-8")
+    p = plan(path)
+    assert [(s.id, s.config.time, s.config.date) for s in p.sources] == [
+        ("road_2010", 2010, "2010-07-01"), ("road_2020", 2020, "2020-07-01")]
+    run(path, workspace=tmp_path / "ws")
+    assert ['[date:"2010-07-01T00:00:00Z"]' in srv.requests[0]["query"],
+            '[date:"2020-07-01T00:00:00Z"]' in srv.requests[1]["query"]] == [True, True]
+
+
+def test_a_failed_download_is_a_pipeline_error_naming_the_source(tmp_path, overpass):
+    srv = overpass(script=(406,))
+    with pytest.raises(PipelineError, match=r"source 'road'.*OSM_CONTACT"):
+        run(_pipeline(tmp_path, srv), workspace=tmp_path / "ws")
+
+
+def test_validation_without_network(tmp_path, overpass):
+    srv = overpass()
+    load(_pipeline(tmp_path, srv, source_extra="margin = 0.5\ndate = \"2015-01-01\""))
+    assert not srv.requests
+    bad = tmp_path / "bad.toml"
+    bad.write_text('schema = 1\n[[source]]\nid = "r"\ntype = "osm"\nquery = "way[x]"\n', encoding="utf-8")
+    with pytest.raises(PipelineError, match="extent"):  # a sources-only osm file needs an area
+        load(bad)
+    neg = _pipeline(tmp_path, srv, source_extra="margin = -1")
+    with pytest.raises(PipelineError, match="margin"):
+        load(neg)
+
+
+def test_distance_to_the_fetched_road_is_right(tmp_path, overpass):
+    from disscube import CubeClient
+
+    srv = overpass()
+    run(_pipeline(tmp_path, srv, source_extra="margin = 0.2"), workspace=tmp_path / "ws")
+    cube = CubeClient(str(tmp_path / "ws" / "catalog.db"), str(tmp_path / "ws" / "store"))
+    d = np.asarray(cube.load("dist_road").values)
+    assert np.isfinite(d).all() and d.min() >= 0 and d.max() > d.min()  # a gradient away from the road

From 0c5a5130889f77207bde9d764237660158c8e94e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9rgio=20Costa?= <noreply@example.com>
Date: Fri, 2 Oct 2026 06:54:47 +0000
Subject: [PATCH 2/2] Add the `dem` source: elevation and slope from SRTM,
 Copernicus or TOPODATA

type = "dem" reads a digital elevation model over the grid (plus a margin) and
registers elevation, or slope in degrees or percent, as a raster source with
checksum and provenance, ready for the zonal operators.

- dem = srtm | copernicus | topodata, or `tiles` with your own GeoTIFFs. Only the
  window over the grid is read (HTTP range requests for SRTM and Copernicus);
  TOPODATA sheets are downloaded once into $DISSCUBE_CACHE/dem.
- Slope is computed in the source, because an operator does not see the
  neighbours of a pixel: bilinear to the UTM zone of the grid's centre at
  `resolution` metres, then central differences (numpy.gradient).
- Voids and values outside -1000..9000 m become NaN. A tile that cannot be
  opened is skipped with a warning; none readable is an error.
- Asking for a Copernicus slope logs a warning: GLO-30 includes the forest
  canopy, so slopes at forest edges are overestimated.
- Tests use synthetic tiles (a plane of known slope, a gzip-ed HGT, a ZIP served
  by a local HTTP server); guide docs/guides/dem.md.

Ports 02c_fetch_dem.py of lab15-reconstruction. Builds on the `osm` commit
(same schema, runner and docs insertion points). Not exercised against the real
tile servers (unreachable from the environment this was written in).

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                  |   6 +
 README.md                     |   1 +
 disscube/pipeline/runner.py   |  17 +-
 disscube/pipeline/schema.py   |  32 +++-
 disscube/sources/dem.py       | 252 ++++++++++++++++++++++++++++++
 docs/guides/dem.md            |  66 ++++++++
 docs/guides/pipeline_files.md |   3 +-
 mkdocs.yml                    |   1 +
 tests/test_dem.py             | 285 ++++++++++++++++++++++++++++++++++
 9 files changed, 656 insertions(+), 7 deletions(-)
 create mode 100644 disscube/sources/dem.py
 create mode 100644 docs/guides/dem.md
 create mode 100644 tests/test_dem.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index dc9ce15..222120e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **`type = "dem"` source**: elevation, or slope in degrees or percent, from SRTM, Copernicus GLO-30
+  or TOPODATA (or your own `tiles`), read over the grid plus a `margin` and registered as a raster
+  source with checksum and provenance. The slope is computed on a metric (UTM) grid by central
+  differences; asking for a Copernicus slope logs a warning, since that model includes the forest
+  canopy. TOPODATA sheets are downloaded once into a cache. Replaces the DEM script of the Lab15
+  reconstruction.
 - **`type = "osm"` source**: ways from OpenStreetMap through the Overpass API, clipped to the grid
   plus a `margin`, registered as a vector source with checksum and provenance (request, server,
   time retrieved, `timestamp_osm_base`, ODbL). The answer is cached by request, so a pipeline reads
diff --git a/README.md b/README.md
index 8643c84..fc92ac6 100644
--- a/README.md
+++ b/README.md
@@ -247,6 +247,7 @@ disscube/
 │   ├── mapbiomas.py  MapBiomas annual land-cover maps (Collection 11, 10 m series)
 │   ├── prodes.py     PRODES deforestation (download + cache, legend from the .qml)
 │   ├── osm.py        OpenStreetMap ways via Overpass (cached, with date snapshots)
+│   ├── dem.py        Elevation and slope from SRTM, Copernicus GLO-30 or TOPODATA
 │   └── classified.py any classified map with its legend, e.g. from SITS
 └── utils.py          Checksums (sha256_file) and BDC tile importer (import_bdc_grids)
 ```
diff --git a/disscube/pipeline/runner.py b/disscube/pipeline/runner.py
index c576101..e1b978d 100644
--- a/disscube/pipeline/runner.py
+++ b/disscube/pipeline/runner.py
@@ -18,6 +18,7 @@
 from disscube.pipeline.schema import (
     BdcSource,
     ClassifiedSource,
+    DemSource,
     DeriveConfig,
     ExportConfig,
     FileSource,
@@ -57,7 +58,7 @@ def base_dir(self) -> Path:
 @dataclass
 class PlannedSource:
     id: str
-    config: FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | UnionSource
+    config: FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | DemSource | UnionSource
     year: int | None = None
 
 
@@ -413,7 +414,7 @@ def _register_source(
 ):
     c = s.config
     # ── Validação defensiva: apenas fontes dinâmicas em janela na nuvem precisam de bbox_geo ──
-    if isinstance(c, (BdcSource, MapbiomasSource, ProdesSource, ClassifiedSource, OsmSource)):
+    if isinstance(c, (BdcSource, MapbiomasSource, ProdesSource, ClassifiedSource, OsmSource, DemSource)):
         if bbox_geo is None:
             raise PipelineError(
                 f"source {c.id!r} ({c.type}): windowed cloud sources require an"
@@ -460,6 +461,18 @@ def _register_source(
                 )
             except (osm.OsmError, ValueError) as exc:
                 raise PipelineError(f"source {c.id!r}: {exc}") from None
+        if isinstance(c, DemSource):
+            from disscube.sources import dem
+
+            try:
+                return dem.register_dem_source(
+                    cube, c.id, bbox_geo, raw, dem=c.dem,
+                    tiles=[_resolve(base, t) for t in c.tiles] if c.tiles else None,
+                    product=c.product, margin=c.margin, resolution=c.resolution,
+                    cache=(base / c.cache) if c.cache else None, name=c.name,
+                )
+            except (dem.DemError, ValueError) as exc:
+                raise PipelineError(f"source {c.id!r}: {exc}") from None
         if isinstance(c, ClassifiedSource):
             from disscube.sources.classified import register_classified_map
 
diff --git a/disscube/pipeline/schema.py b/disscube/pipeline/schema.py
index bd1c577..3b61829 100644
--- a/disscube/pipeline/schema.py
+++ b/disscube/pipeline/schema.py
@@ -150,6 +150,30 @@ class OsmSource(_SourceBase):
     time: int | None = None
 
 
+class DemSource(_SourceBase):
+    """Elevation or slope from a DEM (SRTM, Copernicus GLO-30, TOPODATA), over the grid plus ``margin``.
+
+    Give ``dem`` to download one, or ``tiles`` (GeoTIFFs, paths relative to the pipeline file or URLs
+    GDAL can read) to use your own. ``product`` is ``elevation`` (m), ``slope_deg`` or ``slope_pct``;
+    a slope is computed on a metric (UTM) grid of ``resolution`` metres. Prefer SRTM or TOPODATA to
+    Copernicus for slope: Copernicus includes the forest canopy.
+    """
+
+    type: Literal["dem"]
+    dem: Literal["srtm", "copernicus", "topodata"] | None = None
+    tiles: list[str] | None = None
+    product: Literal["elevation", "slope_deg", "slope_pct"] = "elevation"
+    margin: float = Field(default=0.02, ge=0)
+    resolution: float = Field(default=30.0, gt=0)
+    cache: str | None = None
+
+    @model_validator(mode="after")
+    def _one_input(self):
+        if (self.dem is None) == (self.tiles is None):
+            raise ValueError(f"source {self.id!r}: give exactly one of 'dem' or 'tiles'")
+        return self
+
+
 class UnionSource(_SourceBase):
     """The features of several vector sources of this file as one source."""
 
@@ -158,7 +182,7 @@ class UnionSource(_SourceBase):
 
 
 Source = Annotated[
-    FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | UnionSource,
+    FileSource | BdcSource | MapbiomasSource | ProdesSource | ClassifiedSource | OsmSource | DemSource | UnionSource,
     Field(discriminator="type"),
 ]
 
@@ -193,7 +217,7 @@ class PipelineConfig(_Strict):
     grid: GridConfig | None = None
     extent: list[float] | None = Field(default=None, min_length=4, max_length=4)
     """``[min_lon, min_lat, max_lon, max_lat]`` (WGS84) a sources-only file reads
-    windowed sources over (classified maps, BDC, MapBiomas, PRODES, OpenStreetMap)."""
+    windowed sources over (classified maps, BDC, MapBiomas, PRODES, OpenStreetMap, DEM)."""
     sources_from_catalog: bool = False
     """Let [[derive]] blocks use sources this file does not declare, registered
     in the workspace's catalog by another pipeline file. Off by default, so a
@@ -209,14 +233,14 @@ def _grid_or_extent(self):
 
         # Apenas fontes dinâmicas em janela na nuvem precisam de extent.
         # Fontes 'file' e 'union' já têm extensão definida por seus arquivos locais.
-        windowed_types = {"bdc", "mapbiomas", "prodes", "classified", "osm"}
+        windowed_types = {"bdc", "mapbiomas", "prodes", "classified", "osm", "dem"}
         needs_extent = any(
             getattr(s, "type", None) in windowed_types for s in self.source
         )
 
         if self.grid is None and self.extent is None and needs_extent:
             raise ValueError(
-                "a sources-only file with windowed sources (BDC, MapBiomas, PRODES, OpenStreetMap)"
+                "a sources-only file with windowed sources (BDC, MapBiomas, PRODES, OpenStreetMap, DEM)"
                 " needs `extent` (or a [grid])"
             )
         return self
diff --git a/disscube/sources/dem.py b/disscube/sources/dem.py
new file mode 100644
index 0000000..64e609b
--- /dev/null
+++ b/disscube/sources/dem.py
@@ -0,0 +1,252 @@
+"""Elevation and slope from a digital elevation model (SRTM, Copernicus GLO-30, TOPODATA).
+
+A pipeline names the DEM and the product and gets a raster source with a checksum and
+provenance, ready for the zonal operators (``mean``, ``min``...)::
+
+    [[source]]
+    id      = "slope"
+    type    = "dem"
+    dem     = "srtm"            # srtm | copernicus | topodata
+    product = "slope_deg"       # elevation | slope_deg | slope_pct
+
+    [[derive]]
+    target   = "media_decl"
+    source   = "slope"
+    operator = "mean"
+
+Only the window over the grid (plus ``margin``) is read: SRTM and Copernicus tiles are read
+straight from public cloud storage, TOPODATA sheets are downloaded once into a cache
+(``$DISSCUBE_CACHE/dem``, or ``~/.cache/disscube/dem``). ``tiles`` takes your own GeoTIFFs
+instead.
+
+**Slope** needs the neighbours of every pixel, which an operator that aggregates a cell does
+not have, so it is computed here: the window is resampled (bilinear) to the UTM zone of the
+grid at ``resolution`` metres, so the slope is metric, and taken by central differences
+(``numpy.gradient``), in degrees or percent. That is not Horn's 3 × 3 formula, so values differ
+slightly from ``gdaldem slope``.
+
+**Which DEM.** SRTM (C-band radar, February 2000) and TOPODATA (INPE's refinement of it)
+measure close to the ground; Copernicus GLO-30 (X-band, 2011–2015) is a *surface* model that
+includes the forest canopy, so every forest/pasture edge becomes a ~30 m "cliff" and the slope
+is overestimated exactly where deforestation happens. For slope prefer SRTM or TOPODATA.
+"""
+
+from __future__ import annotations
+
+import logging
+import math
+import os
+import shutil
+import urllib.parse
+import urllib.request
+import zipfile
+from collections.abc import Sequence
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+import rasterio
+from rasterio.crs import CRS
+from rasterio.errors import RasterioIOError
+from rasterio.transform import array_bounds
+from rasterio.warp import Resampling, calculate_default_transform, reproject
+
+from disscube.sources._raster import Window2D, mosaic, read_window, register_raster
+
+log = logging.getLogger("disscube.sources.dem")
+
+#: Copernicus DEM GLO-30, public Cloud-Optimized GeoTIFFs on AWS (no account), named by the SW corner.
+COPERNICUS_URL = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{ns}{lat:02d}_00_{ew}{lon:03d}_00_DEM/"
+                  "Copernicus_DSM_COG_10_{ns}{lat:02d}_00_{ew}{lon:03d}_00_DEM.tif")
+#: SRTM 1 arc-second as republished in the AWS Terrain Tiles ("skadi", public, no key), gzip-ed HGT.
+SRTM_URL = "https://s3.amazonaws.com/elevation-tiles-prod/skadi/{ns}{lat:02d}/{ns}{lat:02d}{ew}{lon:03d}.hgt.gz"
+#: TOPODATA (INPE): 1° × 1.5° sheets named by their north-west corner, e.g. ``03S555`` = 3–4°S, 55.5–54°W. ZN = altitude.
+TOPODATA_URL = "http://www.dsr.inpe.br/topodata/data/geotiff/{sheet}ZN.zip"
+
+DATASETS: dict[str, dict[str, str]] = {
+    "srtm": {"description": "SRTM 1 arc-second (Feb 2000, C-band radar), AWS Terrain Tiles (skadi)",
+             "license": "public domain (NASA/USGS)", "surface": "near ground (canopy partly penetrated)"},
+    "copernicus": {"description": "Copernicus DEM GLO-30 (2011-2015, X-band): a surface model",
+                   "license": "Copernicus DEM, © DLR e.V. and Airbus Defence and Space GmbH; free licence",
+                   "surface": "surface: includes the forest canopy"},
+    "topodata": {"description": "TOPODATA (INPE): SRTM refined to 1 arc-second by kriging",
+                 "license": "INPE, free access", "surface": "near ground (refined SRTM)"},
+}
+PRODUCTS: dict[str, str] = {"elevation": "m", "slope_deg": "degrees", "slope_pct": "percent"}
+#: Plausible elevations in metres; anything else is a void (SRTM marks them -32768).
+VALID_RANGE = (-1000.0, 9000.0)
+
+
+class DemError(RuntimeError):
+    """The elevation data could not be obtained or does not cover the area."""
+
+
+def default_cache_dir() -> Path:
+    """``$DISSCUBE_CACHE/dem``, or ``~/.cache/disscube/dem``."""
+    root = os.environ.get("DISSCUBE_CACHE") or Path.home() / ".cache" / "disscube"
+    return Path(root) / "dem"
+
+
+def utm_epsg(bbox_geo: Sequence[float]) -> int:
+    """EPSG code of the WGS84 / UTM zone of the centre of ``bbox_geo`` (``[w, s, e, n]``)."""
+    lon, lat = (bbox_geo[0] + bbox_geo[2]) / 2, (bbox_geo[1] + bbox_geo[3]) / 2
+    zone = min(max(int((lon + 180) // 6) + 1, 1), 60)
+    return (32600 if lat >= 0 else 32700) + zone
+
+
+def padded(bbox_geo: Sequence[float], margin: float) -> list[float]:
+    w, s, e, n = (float(v) for v in bbox_geo)
+    return [max(w - margin, -180.0), max(s - margin, -90.0), min(e + margin, 180.0), min(n + margin, 90.0)]
+
+
+def _degree_range(lo: float, hi: float) -> range:
+    """Integer degrees of the 1° tiles (named by their lower edge) that ``[lo, hi]`` touches."""
+    first = math.floor(lo)
+    return range(first, max(math.ceil(hi) - 1, first) + 1)
+
+
+def tile_urls(bbox_geo: Sequence[float], dem: str) -> list[str]:
+    """Where the tiles of ``dem`` covering ``bbox_geo`` (``[w, s, e, n]``, WGS84) live."""
+    if dem not in DATASETS:
+        raise ValueError(f"unknown DEM {dem!r}; choose from {sorted(DATASETS)}")
+    w, s, e, n = bbox_geo
+    if dem == "topodata":
+        sheets: set[str] = set()
+        for lat in range(math.ceil(s), math.ceil(n) + 1):  # north edge of the sheet
+            for k in range(math.floor(w / 1.5), math.floor(e / 1.5) + 1):
+                west = -k * 1.5  # west edge, degrees W
+                sheets.add(f"{abs(lat):02d}{'S' if lat <= 0 else 'N'}{round(west * 10):03d}")
+        return [TOPODATA_URL.format(sheet=sh) for sh in sorted(sheets)]
+    template = COPERNICUS_URL if dem == "copernicus" else SRTM_URL
+    return [template.format(ns="N" if lat >= 0 else "S", lat=abs(lat), ew="E" if lon >= 0 else "W", lon=abs(lon))
+            for lat in _degree_range(s, n) for lon in _degree_range(w, e)]
+
+
+def _require_http(url: str) -> str:
+    if urllib.parse.urlparse(url).scheme not in ("http", "https"):
+        raise ValueError(f"DEM URL must be http(s): {url!r}")
+    return url
+
+
+def _download(url: str, folder: Path, timeout: float = 300) -> Path:
+    """Fetch ``url`` once into ``folder`` (atomically); an existing file is reused."""
+    folder.mkdir(parents=True, exist_ok=True)
+    target = folder / url.rsplit("/", 1)[-1]
+    if target.exists():
+        return target
+    log.info("downloading %s", url)
+    tmp = target.with_suffix(target.suffix + ".part")
+    try:
+        with urllib.request.urlopen(_require_http(url), timeout=timeout) as resp, open(tmp, "wb") as out:  # nosec B310
+            shutil.copyfileobj(resp, out, length=1 << 20)
+    except OSError as exc:  # URLError, HTTPError and timeouts are all OSError
+        tmp.unlink(missing_ok=True)
+        raise DemError(f"could not download {url}: {exc}") from exc
+    tmp.replace(target)
+    return target
+
+
+def openable(url: str, dem: str, cache: str | Path | None = None) -> str:
+    """A path GDAL can open for a tile: the COG itself, gzip over HTTP, or a downloaded ZIP."""
+    if dem == "copernicus":
+        return url
+    if dem == "srtm":
+        return f"/vsigzip//vsicurl/{url}" if "://" in url else f"/vsigzip/{url}"
+    archive = _download(url, Path(cache) if cache else default_cache_dir())
+    with zipfile.ZipFile(archive) as zf:
+        tif = next((nm for nm in zf.namelist() if nm.lower().endswith(".tif")), None)
+    if tif is None:
+        raise DemError(f"{archive.name} holds no GeoTIFF")
+    return f"/vsizip/{archive}/{tif}"
+
+
+def read_dem(bbox_geo: Sequence[float], dem: str | None = None, tiles: Sequence[str] | None = None,
+             cache: str | Path | None = None) -> tuple[Window2D, list[str]]:
+    """Elevation (metres, NaN at voids) over ``bbox_geo``, and the tiles it came from.
+
+    ``tiles`` (GeoTIFFs readable by GDAL) replace the download of ``dem``. A tile that cannot be
+    opened is skipped with a warning (the sea has none), but at least one must be read.
+    """
+    if (dem is None) == (tiles is None):
+        raise ValueError("give exactly one of 'dem' or 'tiles'")
+    sources = list(tiles) if tiles is not None else tile_urls(bbox_geo, dem)  # type: ignore[arg-type]
+    windows: list[Window2D] = []
+    used: list[str] = []
+    failures: list[str] = []
+    for src in sources:
+        try:
+            href = src if tiles is not None else openable(src, dem, cache)  # type: ignore[arg-type]
+            win = read_window(href, bbox_geo)
+        except (RasterioIOError, rasterio.errors.WindowError) as exc:
+            failures.append(f"{src}: {exc}")
+            log.warning("tile skipped: %s", failures[-1])
+            continue
+        if win.data.size == 0:
+            continue
+        data = win.data
+        data[(data < VALID_RANGE[0]) | (data > VALID_RANGE[1])] = np.nan
+        windows.append(win)
+        used.append(src)
+    if not windows:
+        raise DemError("no elevation data over the area: " + ("; ".join(failures) or "no tile intersects it"))
+    return mosaic(windows), used
+
+
+def to_utm(window: Window2D, epsg: int, resolution: float) -> Window2D:
+    """Resample ``window`` (bilinear) to the UTM zone ``epsg`` with square ``resolution``-metre pixels."""
+    rows, cols = window.data.shape
+    dst_crs = CRS.from_epsg(epsg)
+    transform, width, height = calculate_default_transform(
+        window.crs, dst_crs, cols, rows, *array_bounds(rows, cols, window.transform), resolution=resolution)
+    out = np.full((height, width), np.nan, dtype="float32")
+    reproject(window.data, out, src_transform=window.transform, src_crs=window.crs, src_nodata=np.nan,
+              dst_transform=transform, dst_crs=dst_crs, dst_nodata=np.nan, resampling=Resampling.bilinear)
+    return Window2D(data=out, transform=transform, crs=dst_crs)
+
+
+def slope(window: Window2D) -> tuple[np.ndarray, np.ndarray]:
+    """Slope in degrees and in percent of a DEM in metric coordinates (central differences)."""
+    dzdy, dzdx = np.gradient(window.data, abs(window.transform.e), abs(window.transform.a))
+    rise = np.hypot(dzdx, dzdy)
+    return np.degrees(np.arctan(rise)).astype("float32"), (100 * rise).astype("float32")
+
+
+def make_product(elevation: Window2D, product: str, bbox_geo: Sequence[float], resolution: float) -> Window2D:
+    """The raster of ``product`` from an elevation window (see :data:`PRODUCTS`)."""
+    if product not in PRODUCTS:
+        raise ValueError(f"unknown DEM product {product!r}; choose from {sorted(PRODUCTS)}")
+    if product == "elevation":
+        return elevation
+    metric = to_utm(elevation, utm_epsg(bbox_geo), resolution)
+    degrees, percent = slope(metric)
+    return Window2D(data=degrees if product == "slope_deg" else percent, transform=metric.transform, crs=metric.crs)
+
+
+def register_dem_source(cube, source_id: str, bbox_geo: Sequence[float], out_dir: str | Path, *,
+                        dem: str | None = None, tiles: Sequence[str] | None = None, product: str = "elevation",
+                        margin: float = 0.02, resolution: float = 30.0, cache: str | Path | None = None,
+                        name: str | None = None):
+    """Read the DEM over the grid, make ``product`` and register it as a raster source.
+
+    The GeoTIFF is ``<out_dir>/<source_id>.tif`` (float32, NaN nodata; elevation in EPSG:4326,
+    slopes in the UTM zone of the grid) with a ``<source_id>.provenance.json`` holding the DEM,
+    its licence, the tiles, the window, the processing and the checksum.
+    """
+    if product not in PRODUCTS:
+        raise ValueError(f"unknown DEM product {product!r}; choose from {sorted(PRODUCTS)}")
+    if dem == "copernicus" and product != "elevation":
+        log.warning("%s: Copernicus GLO-30 includes the forest canopy, so slopes at forest edges are overestimated; "
+                    "SRTM or TOPODATA are closer to the ground", source_id)
+    bbox = padded(bbox_geo, margin)
+    elevation, used = read_dem(bbox, dem, tiles, cache)
+    layer = make_product(elevation, product, bbox, resolution)
+
+    info: dict[str, Any] = DATASETS.get(dem or "", {"description": "user-supplied tiles", "license": "see the data's provider"})
+    processing = ("elevation as read, EPSG:4326" if product == "elevation" else
+                  f"bilinear to UTM {layer.crs.to_string()} at {resolution:g} m, central-difference slope (numpy.gradient)")
+    provenance = {"type": "dem", "dem": dem or "tiles", "description": info["description"], "license": info["license"],
+                  "product": product, "unit": PRODUCTS[product], "tiles": used, "bbox_geo": bbox, "margin": margin,
+                  "processing": processing, "crs": layer.crs.to_string()}
+    if "surface" in info:
+        provenance["surface"] = info["surface"]
+    return register_raster(cube, source_id, layer, out_dir, provenance, name=name, tags=["dem", dem or "tiles", product])
diff --git a/docs/guides/dem.md b/docs/guides/dem.md
new file mode 100644
index 0000000..dbb86ce
--- /dev/null
+++ b/docs/guides/dem.md
@@ -0,0 +1,66 @@
+# Elevation and slope (DEM)
+
+`type = "dem"` reads a digital elevation model over the grid and gives **elevation** or **slope**
+as a raster source, with a checksum and provenance. The usual next step is a zonal operator such
+as `mean`.
+
+```toml
+[grid]
+name = "area"
+bbox = [-54.842, -3.587, -54.459, -3.168]
+resolution = 500
+
+[[source]]
+id      = "slope"
+type    = "dem"
+dem     = "srtm"            # srtm | copernicus | topodata
+product = "slope_deg"       # elevation | slope_deg | slope_pct
+
+[[derive]]
+target   = "mean_slope"
+source   = "slope"
+operator = "mean"
+```
+
+## Fields
+
+| Field | Meaning |
+|---|---|
+| `dem` | `srtm`, `copernicus` or `topodata`. Give this **or** `tiles`. |
+| `tiles` | Your own GeoTIFFs (paths relative to the pipeline file, or URLs GDAL reads) instead of a download. |
+| `product` | `elevation` (m, default), `slope_deg` or `slope_pct`. |
+| `margin` | Degrees added around the grid (default `0.02`): a slope needs the neighbours of the edge pixels. |
+| `resolution` | Metres of the grid on which the slope is computed (default `30`). |
+| `cache` | Folder for downloaded TOPODATA sheets (default `$DISSCUBE_CACHE/dem` or `~/.cache/disscube/dem`). |
+
+## The three DEMs
+
+| `dem` | What | How it is read |
+|---|---|---|
+| `srtm` | SRTM 1 arc-second, February 2000 (C-band radar), as republished in the AWS Terrain Tiles. No account. | Only the window over the grid, through HTTP range requests. |
+| `copernicus` | Copernicus GLO-30, 2011–2015 (X-band). Free licence, © DLR/Airbus. | Only the window, from the public Cloud-Optimized GeoTIFFs on AWS. |
+| `topodata` | INPE's TOPODATA: SRTM refined to 1 arc-second by kriging; Brazil only. | The 1° × 1.5° ZIP sheet is downloaded once into the cache. |
+
+**For slope, prefer SRTM or TOPODATA.** Copernicus GLO-30 is a *surface* model: it includes the
+forest canopy, so every forest/pasture edge looks like a ~30 m cliff and the slope is overestimated
+exactly where deforestation happens. DisSCube logs a warning when you ask for a Copernicus slope.
+In the Lab15 reconstruction, moving from Copernicus to SRTM raised the correlation of the mean slope
+with the original from 0.82 to 0.90.
+
+## How the slope is made
+
+An operator aggregates the pixels of a cell and does not see their neighbours, so the slope is
+computed in the source: the window is resampled (bilinear) to the UTM zone of the grid's centre at
+`resolution` metres, so the slope is metric, and taken by central differences (`numpy.gradient`).
+This is not Horn's 3 × 3 formula, so values differ slightly from `gdaldem slope`. Voids (SRTM marks
+them −32768) and anything outside −1 000 to 9 000 m become NaN, and the slope is NaN around them.
+
+The GeoTIFF is float32 with NaN as nodata: elevation in EPSG:4326, slopes in UTM.
+`raw/<id>.provenance.json` records the DEM, its licence, the tiles, the window, the processing and
+the checksum.
+
+## Tiles are not always there
+
+SRTM and Copernicus have no tile over the open sea. A tile that cannot be opened is skipped with a
+warning; if none can be read, the run stops with the reasons. Check the warning if the grid reaches
+the coast.
diff --git a/docs/guides/pipeline_files.md b/docs/guides/pipeline_files.md
index 5333701..890522f 100644
--- a/docs/guides/pipeline_files.md
+++ b/docs/guides/pipeline_files.md
@@ -71,7 +71,7 @@ resolution = 300
 
 [[source]]                    # one block per source; see the types below
 id = "…"
-type = "file" | "bdc" | "mapbiomas" | "prodes" | "classified" | "osm" | "union"
+type = "file" | "bdc" | "mapbiomas" | "prodes" | "classified" | "osm" | "dem" | "union"
 
 [[derive]]                    # one block per derived variable
 target = "urban_pct"
@@ -98,6 +98,7 @@ relative to the pipeline file, so a folder with its TOML and data is portable.
 | `prodes` | `year`, `url`, `cache` | `disscube.sources.prodes` |
 | `classified` | `path`, `legend` (table or `.qml`/`.json`/`.csv` file), `time`, `nodata`, `producer` | `disscube.sources.classified` |
 | `osm` | `query`, `margin`, `date`, `endpoints`, `cache`, `time` (see [OpenStreetMap](osm.md)) | `disscube.sources.osm` |
+| `dem` | `dem` or `tiles`, `product`, `margin`, `resolution`, `cache` (see [Elevation and slope](dem.md)) | `disscube.sources.dem` |
 | `union` | `of` (ids of vector sources declared before it) | their features in one GeoPackage under `raw/` |
 
 Every source block also accepts `name` and `years`.
diff --git a/mkdocs.yml b/mkdocs.yml
index 662be32..eaee26c 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -44,4 +44,5 @@ nav:
       - MapBiomas: guides/mapbiomas.md
       - "Land-cover maps: PRODES and SITS": guides/landcover.md
       - "OpenStreetMap": guides/osm.md
+      - "Elevation and slope (DEM)": guides/dem.md
       - Custom Grids: guides/grids.md
diff --git a/tests/test_dem.py b/tests/test_dem.py
new file mode 100644
index 0000000..05cdc58
--- /dev/null
+++ b/tests/test_dem.py
@@ -0,0 +1,285 @@
+"""DEM source: tiles, voids, mosaics, slope on a known plane, TOPODATA download, and the pipeline type.
+
+Everything uses synthetic tiles on disk and a local HTTP server; nothing touches the network.
+"""
+
+from __future__ import annotations
+
+import gzip
+import json
+import threading
+import zipfile
+from http.server import BaseHTTPRequestHandler, HTTPServer
+
+import numpy as np
+import pytest
+import rasterio
+from rasterio.transform import from_origin
+
+from disscube import CubeClient
+from disscube.pipeline import PipelineError, load, run
+from disscube.sources import dem
+
+VALID_MAX = dem.VALID_RANGE[1]
+LAB15 = [-54.842, -3.587, -54.459, -3.168]
+M_PER_DEG_LAT = 110_580.0  # one degree of latitude near 3.4°S
+RISE = 0.05                # the test plane climbs 5 cm per metre northward (stays under 9 000 m over a degree)
+EXPECTED_DEG = float(np.degrees(np.arctan(RISE)))  # 2.862°
+
+
+def _write_tif(path, west, north, size_deg, pixel, z):
+    """A 1-band float32 GeoTIFF in EPSG:4326 whose upper-left corner is (west, north)."""
+    rows, cols = z.shape
+    with rasterio.open(path, "w", driver="GTiff", height=rows, width=cols, count=1, dtype="float32",
+                       crs="EPSG:4326", transform=from_origin(west, north, pixel, pixel), nodata=np.nan) as d:
+        d.write(z.astype("float32"), 1)
+    return str(path)
+
+
+def _plane(west, south, east, north, pixel):
+    """Elevation that climbs northward by RISE m per metre, in a tile covering the given box."""
+    lat = north - (np.arange(round((north - south) / pixel)) + 0.5) * pixel
+    cols = round((east - west) / pixel)
+    return np.tile(((lat - south) * M_PER_DEG_LAT * RISE)[:, None], (1, cols))
+
+
+@pytest.fixture
+def plane_tile(tmp_path):
+    box = (-54.95, -3.70, -54.35, -3.05)
+    return _write_tif(tmp_path / "plane.tif", box[0], box[3], 0.0, 0.0005, _plane(*box, 0.0005))
+
+
+@pytest.fixture
+def cube(tmp_path):
+    return CubeClient(str(tmp_path / "cat.db"), str(tmp_path / "store"))
+
+
+# --- naming and geometry ---------------------------------------------------------------------
+
+def test_tile_names_for_the_lab15_area():
+    box = dem.padded(LAB15, 0.02)
+    assert dem.tile_urls(box, "srtm") == ["https://s3.amazonaws.com/elevation-tiles-prod/skadi/S04/S04W055.hgt.gz"]
+    assert dem.tile_urls(box, "copernicus")[0].endswith("Copernicus_DSM_COG_10_S04_00_W055_00_DEM.tif")
+    assert dem.tile_urls(box, "topodata") == ["http://www.dsr.inpe.br/topodata/data/geotiff/03S555ZN.zip"]
+
+
+def test_tiles_across_boundaries_and_exact_edges():
+    names = [u.rsplit("/", 1)[1] for u in dem.tile_urls([-55.2, -4.5, -53.8, -2.5], "srtm")]
+    assert sorted(names) == sorted(f"S{la:02d}W{lo:03d}.hgt.gz" for la in (5, 4, 3) for lo in (56, 55, 54))
+    assert len(dem.tile_urls([-55.0, -4.0, -54.0, -3.0], "srtm")) == 1  # edges on whole degrees: one tile
+    assert [u.rsplit("/", 1)[1] for u in dem.tile_urls([10.2, 50.2, 10.8, 50.8], "srtm")] == ["N50E010.hgt.gz"]
+    with pytest.raises(ValueError, match="unknown DEM"):
+        dem.tile_urls(LAB15, "aster")
+
+
+def test_utm_zone_of_the_grid():
+    assert dem.utm_epsg(LAB15) == 32721
+    assert dem.utm_epsg([10, 50, 11, 51]) == 32632
+    assert dem.utm_epsg([-70, 40, -69, 41]) == 32619
+
+
+# --- reading ---------------------------------------------------------------------------------
+
+def test_slope_of_a_known_plane(plane_tile):
+    elevation, used = dem.read_dem(dem.padded(LAB15, 0.02), tiles=[plane_tile])
+    assert used == [plane_tile]
+    deg = dem.make_product(elevation, "slope_deg", LAB15, 30)
+    pct = dem.make_product(elevation, "slope_pct", LAB15, 30)
+    assert deg.crs.to_epsg() == 32721 and abs(deg.transform.a) == 30
+    inner = (slice(40, -40), slice(40, -40))  # away from the edges of the window
+    assert np.isfinite(elevation.data).all() and np.nanmax(elevation.data) < VALID_MAX  # the plane is not read as voids
+    assert np.nanmean(deg.data[inner]) == pytest.approx(EXPECTED_DEG, abs=0.05)
+    assert np.nanmean(pct.data[inner]) == pytest.approx(100 * RISE, abs=0.15)
+    valid = deg.data[inner][np.isfinite(deg.data[inner])]
+    assert np.mean(np.abs(valid - EXPECTED_DEG) > 0.1) < 0.01  # a plane has one slope (edge pixels aside)
+
+
+def test_elevation_stays_in_wgs84_and_voids_become_nan(tmp_path):
+    z = np.full((400, 400), 120.0)
+    z[100:110, 100:110] = -32768.0
+    z[200, 200] = 12000.0
+    path = _write_tif(tmp_path / "v.tif", -54.9, -3.1, 0.0, 0.0005, z)
+    win, _ = dem.read_dem([-54.88, -3.28, -54.72, -3.12], tiles=[path])
+    assert win.crs.to_epsg() == 4326
+    assert np.isnan(win.data).sum() > 0 and np.nanmax(win.data) == 120.0
+
+
+def test_two_tiles_make_one_continuous_window(tmp_path):
+    left = _write_tif(tmp_path / "l.tif", -55.0, -3.0, 0.0, 0.001, np.full((1000, 1000), 10.0))
+    right = _write_tif(tmp_path / "r.tif", -54.0, -3.0, 0.0, 0.001, np.full((1000, 1000), 20.0))
+    win, used = dem.read_dem([-54.1, -3.6, -53.9, -3.4], tiles=[left, right])  # straddles -54°
+    assert used == [left, right] and set(np.unique(win.data[~np.isnan(win.data)])) == {10.0, 20.0}
+    assert not np.isnan(win.data).any()
+
+
+def test_a_missing_tile_is_skipped_with_a_warning_but_all_missing_is_an_error(tmp_path, plane_tile, caplog):
+    with caplog.at_level("WARNING", logger="disscube.sources.dem"):
+        _, used = dem.read_dem(dem.padded(LAB15, 0.02), tiles=[str(tmp_path / "nope.tif"), plane_tile])
+    assert used == [plane_tile] and "tile skipped" in caplog.text
+    with pytest.raises(dem.DemError, match="no elevation data"):
+        dem.read_dem(LAB15, tiles=[str(tmp_path / "nope.tif")])
+
+
+def test_exactly_one_of_dem_or_tiles():
+    with pytest.raises(ValueError, match="exactly one"):
+        dem.read_dem(LAB15)
+    with pytest.raises(ValueError, match="exactly one"):
+        dem.read_dem(LAB15, dem="srtm", tiles=["a.tif"])
+
+
+# --- SRTM (gzip-ed HGT) and Copernicus names, through the real download paths -----------------
+
+def test_srtm_hgt_gz_is_read_through_vsigzip(tmp_path, monkeypatch):
+    n = 1201
+    z = np.full((n, n), 300, dtype="int16")
+    z[500:520, 500:520] = -32768  # a void
+    hgt = tmp_path / "S04W055.hgt"
+    with rasterio.open(hgt, "w", driver="SRTMHGT", height=n, width=n, count=1, dtype="int16", crs="EPSG:4326",
+                       transform=from_origin(-55, -3, 1 / 1200, 1 / 1200)) as d:
+        d.write(z, 1)
+    (tmp_path / "S04W055.hgt.gz").write_bytes(gzip.compress(hgt.read_bytes()))
+    monkeypatch.setattr(dem, "SRTM_URL", str(tmp_path) + "/{ns}{lat:02d}{ew}{lon:03d}.hgt.gz")
+    win, used = dem.read_dem(dem.padded(LAB15, 0.02), dem="srtm")
+    assert used[0].endswith("S04W055.hgt.gz")
+    assert np.nanmax(win.data) == 300 and np.isnan(win.data).sum() > 0 and win.data.shape[0] > 100
+
+
+def test_copernicus_warns_about_the_canopy_when_asked_for_slope(tmp_path, plane_tile, cube, monkeypatch, caplog):
+    monkeypatch.setattr(dem, "COPERNICUS_URL", plane_tile.replace("plane.tif", "{ns}{lat:02d}{ew}{lon:03d}.tif"))
+    (tmp_path / "S04W055.tif").write_bytes((tmp_path / "plane.tif").read_bytes())
+    with caplog.at_level("WARNING", logger="disscube.sources.dem"):
+        dem.register_dem_source(cube, "slope", LAB15, tmp_path / "raw", dem="copernicus", product="slope_deg")
+    assert "canopy" in caplog.text
+    prov = json.loads((tmp_path / "raw" / "slope.provenance.json").read_text())
+    assert prov["dem"] == "copernicus" and "canopy" in prov["surface"]
+    caplog.clear()
+    with caplog.at_level("WARNING", logger="disscube.sources.dem"):
+        dem.register_dem_source(cube, "elev", LAB15, tmp_path / "raw", dem="copernicus", product="elevation")
+    assert "canopy" not in caplog.text  # the warning is about slope only
+
+
+# --- TOPODATA: a ZIP downloaded once ---------------------------------------------------------
+
+@pytest.fixture
+def topodata_server(tmp_path, monkeypatch):
+    sheet = tmp_path / "03S555ZN.tif"
+    _write_tif(sheet, -55.5, -3.0, 0.0, 0.0005, _plane(-55.5, -4.0, -54.0, -3.0, 0.0005))
+    archive = tmp_path / "03S555ZN.zip"
+    with zipfile.ZipFile(archive, "w") as zf:
+        zf.write(sheet, "03S555ZN.tif")
+    served: list[str] = []
+
+    class Handler(BaseHTTPRequestHandler):
+        def do_GET(self):
+            served.append(self.path)
+            if self.path.endswith("03S555ZN.zip"):
+                body = archive.read_bytes()
+                self.send_response(200)
+                self.send_header("Content-Length", str(len(body)))
+                self.end_headers()
+                self.wfile.write(body)
+            else:
+                self.send_error(404)
+
+        def log_message(self, *a):
+            pass
+
+    server = HTTPServer(("127.0.0.1", 0), Handler)
+    threading.Thread(target=lambda: server.serve_forever(poll_interval=0.01), daemon=True).start()
+    monkeypatch.setattr(dem, "TOPODATA_URL", f"http://127.0.0.1:{server.server_port}/{{sheet}}ZN.zip")
+    yield served
+    server.shutdown()
+    server.server_close()
+
+
+def test_topodata_zip_is_downloaded_once_and_cached(tmp_path, topodata_server):
+    cache = tmp_path / "cache"
+    win, _ = dem.read_dem(dem.padded(LAB15, 0.02), dem="topodata", cache=cache)
+    assert np.nanmax(win.data) > 0 and (cache / "03S555ZN.zip").exists()
+    dem.read_dem(dem.padded(LAB15, 0.02), dem="topodata", cache=cache)
+    assert len(topodata_server) == 1
+
+
+def test_topodata_download_failure_is_a_clear_error(tmp_path, topodata_server, monkeypatch):
+    monkeypatch.setattr(dem, "TOPODATA_URL", dem.TOPODATA_URL.replace("{sheet}ZN", "{sheet}XX"))
+    with pytest.raises(dem.DemError, match="could not download"):
+        dem.read_dem(dem.padded(LAB15, 0.02), dem="topodata", cache=tmp_path / "c2")
+    assert not list((tmp_path / "c2").glob("*.part"))
+
+
+def test_download_only_accepts_http():
+    with pytest.raises(ValueError, match="http"):
+        dem._download("file:///etc/passwd", dem.default_cache_dir())
+
+
+# --- in a pipeline ---------------------------------------------------------------------------
+
+def _pipeline(tmp_path, source_extra=""):
+    path = tmp_path / "dem.toml"
+    path.write_text(f"""
+schema = 1
+[grid]
+name = "g"
+crs = "EPSG:4326"
+bbox = [-54.8415, -3.5865, -54.459, -3.168]
+resolution = 0.0045
+[[source]]
+id = "slope"
+type = "dem"
+tiles = ["plane.tif"]
+product = "slope_deg"
+{source_extra}
+[[derive]]
+target = "media_decl"
+source = "slope"
+operator = "mean"
+""", encoding="utf-8")
+    return path
+
+
+def test_pipeline_derives_the_mean_slope_per_cell(tmp_path, plane_tile):
+    report = run(_pipeline(tmp_path), workspace=tmp_path / "ws")
+    assert report.derived
+    cube = CubeClient(str(tmp_path / "ws" / "catalog.db"), str(tmp_path / "ws" / "store"))
+    values = np.asarray(cube.load("media_decl").values, dtype=float)
+    assert np.isfinite(values).mean() > 0.99
+    assert np.nanmean(values) == pytest.approx(EXPECTED_DEG, abs=0.1)
+    prov = json.loads((tmp_path / "ws" / "raw" / "slope.provenance.json").read_text())
+    assert prov["type"] == "dem" and prov["product"] == "slope_deg" and prov["unit"] == "degrees"
+    assert prov["crs"] == "EPSG:32721" and prov["checksum"].startswith("sha256:")
+    assert "central-difference" in prov["processing"] and prov["pipeline"]["file"].endswith("dem.toml")
+
+
+def test_same_tiles_same_spec_hash_and_a_new_resolution_changes_it(tmp_path, plane_tile):
+    def spec_hash(ws):
+        cube = CubeClient(str(ws / "catalog.db"), str(ws / "store"))
+        return next(d.spec_hash for d in cube.catalog.search_derived_variables() if d.name == "media_decl")
+
+    path = _pipeline(tmp_path)
+    run(path, workspace=tmp_path / "a")
+    run(path, workspace=tmp_path / "b")
+    run(_pipeline(tmp_path, source_extra="resolution = 60"), workspace=tmp_path / "c")
+    assert spec_hash(tmp_path / "a") == spec_hash(tmp_path / "b") != spec_hash(tmp_path / "c")
+
+
+def test_a_missing_tile_file_is_a_pipeline_error_naming_the_source(tmp_path):
+    with pytest.raises(PipelineError, match=r"source 'slope'.*no elevation data"):
+        run(_pipeline(tmp_path), workspace=tmp_path / "ws")
+
+
+def test_validation(tmp_path, plane_tile):
+    load(_pipeline(tmp_path))  # no network, no tile read
+    bad = [
+        ('type = "dem"\ntiles = ["plane.tif"]', 'type = "dem"'),                    # neither dem nor tiles
+        ('tiles = ["plane.tif"]', 'tiles = ["plane.tif"]\ndem = "srtm"'),           # both
+        ('product = "slope_deg"', 'product = "aspect"'),                            # unknown product
+        ('tiles = ["plane.tif"]', 'dem = "aster"'),                                 # unknown dem
+        ('tiles = ["plane.tif"]', 'tiles = ["plane.tif"]\nmargin = -1'),           # negative margin
+        ('tiles = ["plane.tif"]', 'tiles = ["plane.tif"]\nresolution = 0'),        # non-positive resolution
+    ]
+    text = _pipeline(tmp_path).read_text()
+    for old, new in bad:
+        broken = tmp_path / "broken.toml"
+        broken.write_text(text.replace(old, new), encoding="utf-8")
+        with pytest.raises(PipelineError):
+            load(broken)