# Petrel Flow Accumulation Layer v0.1 — Methodology ## Data pipeline ``` Copernicus GLO-30 DEM (© DLR/Airbus, provided under COPERNICUS by the EU + ESA) │ ▼ read each 10° land tile with a 1° (~110 km) halo of neighbour DEM │ ▼ pysheds: fill_pits → fill_depressions → resolve_flats → flowdir │ ▼ D8 flow accumulation (upstream cell count) ── accumulation (intermediate) │ ┌───────────────┴───────────────┐ ▼ ▼ channels = accumulation > 1000 cells TWI = ln( (area+1)·cell_area / tan slope ) │ │ ▼ ▼ Euclidean distance transform → metres ── twi │ ▼ ── distance_to_drainage │ ▼ crop halo → tile core; write per-tile sub-layers, mosaic to global COG │ ▼ coast_clip via petrelml.utils.coast_clip (gshhg_f, all_touched) │ ▼ compact uint16 reencode with PETREL_SCALE tag (cap distance 50 km, TWI 35) │ ▼ GeoTIFF COG w/ overviews ``` **Single-input, Petrel-owned.** The entire derivation starts and ends with the Copernicus GLO-30 DEM. No HydroSHEDS, HydroRIVERS, or other share-alike / EULA-flow-down drainage source is touched at any stage. This is deliberate — see [`data_sources.md`](data_sources.md) and [`model_card.md`](model_card.md) for why single-input provenance is what makes the Layer both giftable and commercially clean. ## Flow derivation — pysheds D8 Per DEM window, the drainage grid is computed with the standard pysheds depression-handling chain: 1. **`fill_pits`** — remove single-cell sinks. 2. **`fill_depressions`** — flood multi-cell depressions to their spill level so flow can leave every closed basin. 3. **`resolve_flats`** — impose a gradient across flat areas so flow directions are defined. 4. **`flowdir`** — D8 flow direction (each cell drains to its steepest of 8 neighbours). 5. **`accumulation`** — upstream contributing cell count per cell. Accumulation is expressed in **cells** (each cell's own contribution counts as 1). At 9 arcsec the physical cell area is latitude-scaled from the grid resolution; TWI converts the cell count to a contributing **area** before taking the log. ## Halo tiling and the cross-tile caveat Flow accumulation is inherently global: water crosses tile boundaries, and a tile computed in isolation under-counts flow entering from outside. v0.1 uses **halo-tiled** derivation — each 10° tile is read together with a **1° (~110 km) halo** of surrounding DEM, the flow grid is derived on the halo'd array, and the result is cropped back to the tile core. Consequences, stated honestly: - **Correct** for any catchment whose headwaters lie within the halo, and for the **drainage network structure** everywhere — channel positions, distance-to-drainage, and TWI depend on the local topology the halo captures, not on the absolute upstream count. - **Under-counted** absolute accumulation for major **trans-halo** rivers at tile interiors far from their source (the source lies beyond the 110 km halo). This is exactly why raw accumulation is held back as an intermediate in v0.1 and only the topology-driven derivatives ship. - A **Barnes-style global two-pass accumulation** (2016) — which reconciles per-tile partial sums across the whole globe in a second pass — is the flagged v1 upgrade that removes this limitation. Lahar routing (the other consumer of this core) is unaffected, since lahar reach (<80 km) is far shorter than the halo. ## Distance to drainage The channel network is the boolean mask `accumulation > 1000 cells`. Distance to drainage is the **Euclidean distance transform** of the complement of that mask (`scipy.ndimage.distance_transform_edt`), converted from pixels to **metres** using a latitude-scaled metres-per-pixel factor (`res° × 111 320 × √cos(lat)`). Cells on a channel read 0; the value grows with straight-line distance to the nearest channel. The far tail (deserts, endorheic interiors, any tile with no channel) is capped at **50 km** (`cap_physical = 50000 m`) for encoding. ## Topographic wetness index (TWI) ``` TWI = ln( (accumulation + 1) · cell_area / tan(slope) ) ``` - `(accumulation + 1) · cell_area` is the upslope contributing area in m². - `slope` is the local gradient magnitude of the DEM (`np.gradient` / metres-per-pixel), floored at 1e-4 rad to avoid a divide-by-zero on flats. - TWI is unitless and typically spans ~0–30; it is capped at **35** (`cap_physical = 35`) for encoding. High TWI marks flat, water-gathering, saturation-prone terrain; low TWI marks steep, well-drained slopes. ## Coastline masking The canonical `coast_clip` step (GSHHG-f full-resolution shoreline, `all_touched=True`) removes ocean cells after the global raster is assembled. Any cell touching real land keeps its value; open ocean and — by Petrel standard — Antarctica (GSHHG L1 only) return NoData. Because the product ships at its native derivation grid (no coarse-to-fine downscale), there is no Carolina-Beach-class coastal-hole regression here; the audit still verifies coastal land cells resolve to data and ocean / Antarctica resolve to NoData. ## Sanity checks the COG must pass (promote audit) 1. **PETREL_COAST_CLIP + PETREL_SCALE tags present** — enforces the canonical mask + encoding steps ran. 2. **Structure sanity** — distance-to-drainage reads ~0 along major river lines and high on ridgelines; TWI reads high on valley floors and low on steep slopes. 3. **Mask correctness** — Antarctica and open ocean are NoData; coastal land cells resolve to data. 4. **Range** — distance within [0, 50 000] m; TWI within [~0, 35]. ## Versioning - `v0.1` — pysheds D8 on the Copernicus GLO-30 DEM, halo-tiled at 9 arcsec, shipping distance-to-drainage + TWI. First shipped. - `v1` — Barnes-style global two-pass accumulation (removes the trans-halo under-count); adds a log-encoded raw-accumulation sub-layer. - Later — 3 arcsec (90 m) regional derivations where the DEM and downstream hazards justify the resolution.