feat(cloud_io): Dask / chunked-tile backend for continental DEMs (P30)

[MAfarrag]

Context

Continental-scale DEMs are 10⁹–10¹⁰ cells; even at int32 that's 4–40 GB, exceeding typical workstation memory.
The Phase 1–3 algorithms assume the full DEM fits in RAM. Phase 4 unlocks tiled operation by introducing
window-based chunking primitives and (eventually) a dask.array-backed DEM that processes chunks lazily.

Problem / Current Behaviour

• No chunking primitive existed pre-Phase-4 — every operation loaded the full raster.
• Phase 1–3 inner kernels are in-memory only; large rasters trip workstation memory limits.

Affected locations

File	Symbol	Notes
src/digitalrivers/cloud_io.py	tile_windows(...)	shipped — window iterator.
src/digitalrivers/cloud_io.py	dask_backend(...)	pending — currently a stub.

Motivation Example

from digitalrivers import DEM
from digitalrivers.cloud_io import tile_windows

dem = DEM.read_file("continental_30m.tif")

# Shipped: iterate over tile windows for window-based processing.
for window in tile_windows(dem, tile_size_cells=2048):
    tile = dem.read_array(window=window)
    ...

# Pending: full Dask-backed pipeline.
# da_dem = DEM.from_dask(dask_array, transform=dem.transform, crs=dem.crs, nodata=-9999)
# filled = da_dem.fill_depressions()   # routes through tile-aware kernels

Proposed Solution

from digitalrivers.cloud_io import tile_windows, dask_backend

# Shipped — window iterator over a Dataset.
tile_windows(dataset, tile_size_cells: int = 2048) -> Iterator[Window]

# Pending — Dask backend.
class DEM(Dataset):
    @classmethod
    def from_dask(cls, da: "dask.array.Array", transform, crs, nodata) -> "DEM": ...
    def to_dask(self, chunks: tuple[int, int] | int = 2048) -> "dask.array.Array": ...

Implementation strategy:

• Embarrassingly-parallel kernels (slope, hillshade, fine-resolution flow direction, building burn,
  hydroflatten of a single polygon) dispatch on input type via dask.array.overlap.map_overlap.
• Topological-order kernels (flow accumulation, watershed, basin labelling, HAND, Strahler) require
  tile ghosting + cross-tile boundary-flow exchange (Wallis 2009 reduction). Implement via map_blocks
  per tile + iterative outer convergence loop on the boundaries.
• Cap each tile at tile_size_cells = 2048**2 (16 MB float32).
• Optional digitalrivers.config.set_chunk_size(int) / set_n_workers(int) runtime config.

References

• Dask documentation: https://docs.dask.org/en/stable/
• xarray-spatial for reference patterns (map_overlap stencil kernels, _boundary_store.py for cross-tile
  exchange).
• rioxarray open_rasterio for lazy chunked GeoTIFF loading.
• Wallis et al. 2009 — iterative method for adaptive parallelisation of flow accumulation.
• Barnes 2014 — parallel Priority-Flood for trillion-cell DEMs.

Out of Scope

• Cloud-streamed I/O (S3/GCS reads) — separate issue P31 (feat(cloud_io): cloud-optimised raster I/O (COG / Zarr / S3 / GCS) (P31) #45); builds on this Dask backend.
• GPU kernels (CuPy / Dask-CuPy) — Phase 5+.

Effort Estimate

Size: XL (8–12 days)
Rationale: Topological-order tile stitching is the hard part. Flow accumulation alone is 3–5 days.

Definition of Done

• cloud_io.tile_windows(...) ships and is exercised by tests.
• Pending: DEM.from_dask(...) / to_dask(...) constructors.
• Pending: 10000×10000 float32 DEM end-to-end (fill → routing → accumulation → streams) in < 5 minutes
  on an 8-core workstation, peak RSS <= 4 GB.
• Pending: 100000×100000 DEM processable on a 16-core, 64 GB workstation in < 1 hour via Dask local
  scheduler.
• Pending: Cross-tile consistency — Dask-backed output bitwise-equal to in-memory output.
• Pending: dask.distributed cluster support; CI smoke test.
• Pending: Coello (4000×4000) regression — Dask backend no more than 30% slower than in-memory.

Status: Partial on branch feat/phase-4. Shipped: tile_windows window iterator. Pending: full Dask
graph construction for parallel kernels. See tests/test_cloud_io.py.