San Juan Riparian Watch: how the pipeline actually works

01Satellite ETL: making data selection queryable

The premise of the whole project is a phenological one. In a semi-arid basin, riparian phreatophytes — cottonwood, willow, and the invasive tamarisk and Russian olive — stay green into the dry season because their roots reach groundwater, while the surrounding upland browns out. That contrast is the signal. Capturing it requires a time series, not a scene, which is why the entry point is a STAC query rather than a file download.

The searcher opens the Planetary Computer catalog with planetary_computer.sign_inplace as a pystac-client modifier, so every asset href returned by the catalog is signed lazily at access time. The consequence is that the data-selection step — area of interest, time window, cloud threshold — is a declarative record in code rather than an untracked manual download.

class PlanetaryComputerSearcher:
    def __init__(self, stac_url: str = STAC_API_URL) -> None:
        self._catalog = pystac_client.Client.open(
            stac_url, modifier=planetary_computer.sign_inplace,
        )

# Sentinel-2: cloud cover is a *server-side* filter, not a post-hoc drop
items = searcher.search(
    SENTINEL_2_COLLECTION, request.bbox, request.date_range,
    query={"eo:cloud_cover": {"lt": request.max_cloud_cover}},
)
cube = stac_load(
    items, bands=bands, bbox=request.bbox, crs=STORAGE_CRS,
    resolution=_deg_resolution(request.resolution_m, request.bbox),
    chunks={},              # dask-backed, lazy until .compute()
    groupby="solar_day",   # mosaic same-day overlapping orbits
)
if "SCL" in cube:
    cube = _apply_scl_mask(cube)

Two masking layers matter. The eo:cloud_cover < 20 query drops whole scenes server-side; _apply_scl_mask then drops individual bad pixels using Sentinel-2's scene classification layer, treating no-data, saturated, cloud-shadow, medium- and high-probability cloud, cirrus and snow as invalid:

SCL_INVALID = frozenset({0, 1, 3, 8, 9, 10, 11})

def _apply_scl_mask(reflectance, scl):
    valid = ~scl.isin(list(SCL_INVALID))
    return reflectance.where(valid)   # invalid → NaN, propagates to stats

Sentinel-1 RTC is loaded through the same path with query=None — radar is cloud-independent, so applying an optical cloud filter to it would be a category error. The land-cover grids (ESA WorldCover, io-lulc-9-class) are loaded with like=cube, which snaps them onto the Sentinel-2 geobox so every product shares one pixel grid. No resampling bugs downstream because there is no downstream resampling.

02Feature engineering: 22 numbers per pixel

A single NDVI value cannot separate riparian woodland from an irrigated alfalfa field — both are green in July. What separates them is behaviour over the season, canopy water content, and spatial texture. The feature stack encodes all three.

Five spectral indices are computed per timestep, each chosen for a physical reason rather than because it was available: NDVI for greenness and density, NDMI for canopy water using the SWIR band, NDRE for chlorophyll via the red edge, EVI for a soil- and aerosol-resistant greenness, and kNDVI as a kernel transform that resists saturation in dense canopy.

def ndvi(nir, red):        return _safe_ratio(nir - red, nir + red)
def ndmi(nir, swir1):      return _safe_ratio(nir - swir1, nir + swir1)
def ndre(nir, rededge):    return _safe_ratio(nir - rededge, nir + rededge)

def evi(nir, red, blue):
    denom = nir + 6.0 * red - 7.5 * blue + 1.0
    return 2.5 * _safe_ratio(nir - red, denom)

def kndvi(nir, red):       return np.tanh(ndvi(nir, red) ** 2)

Each index is then collapsed over the time axis into four statistics. The amplitude (p90 − p10) is the phenology term — it is what tells a groundwater-fed cottonwood stand apart from a field that greens up and is then cut. The p10/p90 percentiles are used instead of min/max because a single unmasked cloud edge would dominate an extremum.

for name, da in index_cube.data_vars.items():
    median = da.median(dim="time", skipna=True)
    p10 = da.quantile(0.10, dim="time", skipna=True).drop_vars("quantile")
    p90 = da.quantile(0.90, dim="time", skipna=True).drop_vars("quantile")
    out[f"{name}_median"]    = median
    out[f"{name}_p10"]       = p10
    out[f"{name}_p90"]       = p90
    out[f"{name}_amplitude"] = p90 - p10   # seasonal swing

That yields 5 indices × 4 statistics = 20 features, plus two texture bands computed on the median NDVI with a 3×3 window (local standard deviation and local range, a cheap scipy.ndimage proxy for GLCM texture). Riparian gallery forest is structurally rough; an alfalfa pivot is smooth. Total: 22 features, assembled into a FeatureStack where a pixel is valid only if every one of its features is finite.

03Labels without a labeller, and a terrain prior

There is no hand-digitised training set for this basin. The pipeline manufactures one from the intersection of independent global products, then constrains the model's output with a physical terrain envelope so that the machine learning cannot invent riparian vegetation on a hillside.

Weak labels: woody ∧ near-water

The operating definition — and it is worth being explicit, because it drives every number that follows — is that riparian means woody vegetation growing near water, which is not the same thing as wetland. Positives are woody cover within 100 m of water in either land cover product, or a wetland class. Negatives must be upland and beyond 200 m from water. The 100–200 m band is deliberately left unlabelled: it is exactly where the two classes blur, and forcing a label there teaches the model noise.

def near_water_mask(water, resolution_m, dist_m):
    dist_px = ndimage.distance_transform_edt(~water)
    return (dist_px * resolution_m) <= dist_m

water = _isin(wc, WC_WATER, shape) | _isin(io, IO_WATER, shape)
near  = near_water_mask(water, resolution_m, NEAR_WATER_M)   # 100 m
far   = ~near_water_mask(water, resolution_m, FAR_WATER_M)   # 200 m

wc_hit  = (_isin(wc, WC_WOODY, shape) & near) | _isin(wc, WC_WETLAND, shape)
io_hit  = (_isin(io, IO_WOODY, shape) & near) | _isin(io, IO_WETLAND, shape)
nwi_hit = nwi_mask if nwi_mask is not None else np.zeros(shape, dtype=bool)

# agreement (0–3) doubles as a confidence score
agreement = wc_hit.astype(np.int8) + io_hit.astype(np.int8) + nwi_hit.astype(np.int8)
positive  = agreement >= 1
upland    = _isin(wc, WC_UPLAND, shape) | _isin(io, IO_UPLAND, shape)
negative  = (~positive) & upland & far

label = np.full(shape, LABEL_EXCLUDE, dtype=np.int8)   # -1 = don't train here
label[negative] = LABEL_NEGATIVE
label[positive] = LABEL_POSITIVE

The HAND envelope

Height Above Nearest Drainage is a terrain-derived measure: for each cell, the vertical distance to the nearest cell on the drainage network, following flow direction. Low HAND means hydrologically connected valley bottom. It is computed from the 3DEP/Copernicus DEM through a standard pysheds conditioning chain (fill pits → fill depressions → resolve flats → flow direction → accumulation), loaded with like=cube so it lands on the same grid.

This is the structural replacement for the fixed hydrology buffer the project started with. Rather than asserting "riparian is 30 m from the centreline," it asks the terrain where water can actually reach, and uses that as a container. Cells with HAND ≤ 8 m are candidate corridor; the model's output is intersected with it:

riparian = (prob_grid >= threshold) & valid_mask
if envelope is not None:
    riparian &= envelope    # upland pixels the model scored high are dropped

04The model, and the validation that makes it honest

The classifier is deliberately boring — a 300-tree random forest with balanced class weights, because riparian is the minority class by a wide margin. The interesting engineering is not the model. It is the cross-validation.

clf = RandomForestClassifier(
    n_estimators=n_estimators, max_depth=max_depth,
    random_state=random_state,
    class_weight="balanced",   # riparian is the minority class
    n_jobs=-1,
)
clf.fit(features, labels.astype(int))

Why a random split would lie

Satellite pixels are spatially autocorrelated: a pixel and its neighbour are nearly the same observation. Under a random train/test split, almost every test pixel has a near-duplicate in the training set, and the model scores brilliantly by memorising. The fix is to hold out whole spatial blocks. Pixels are binned into ~2 km tiles by integer-flooring their coordinates, and those tile IDs become groups for GroupKFold:

def assign_spatial_folds(lats, lons, block_deg=0.02):   # ~0.02 deg ≈ 2 km
    row = np.floor(lats / block_deg).astype(np.int64)
    col = np.floor(lons / block_deg).astype(np.int64)
    return row * 100_000 + col      # block id → GroupKFold group

Every reported metric — precision, recall, F1, PR-AUC, ROC-AUC — comes from folds where an entire geographic block was unseen during training. _safe_auc returns NaN rather than a misleading number when a fold happens to contain a single class.

05Independent truth: NMRipMap

Weak labels validate against weak labels, which proves nothing. The project pulls NMRipMap v2.0 Plus — a professionally mapped riparian inventory from the New Mexico Natural Heritage program — live from an ArcGIS MapServer, rasterises it onto the same grid, and compares pixel-wise.

inter = int((reference & prediction).sum())
union = int((reference | prediction).sum())
tp = inter
fp = (prediction & ~reference).sum()
fn = (reference & ~prediction).sum()

iou       = inter / union
precision = tp / (tp + fp)
recall    = tp / (tp + fn)
f1        = 2 * precision * recall / (precision + recall)

Running that comparison exposed the weak-label problem. On the Malpais tile the weak-label model reached F1 0.71; on the Animas tile — an agricultural valley where irrigated fields sit within 100 m of water and satisfy "woody-near-water" — it reached approximately zero. The labels were teaching the model to find farms.

The response was to make the label source a parameter and train directly on the reference polygons where they exist, writing a distinct model_version so the two models are never confused in the database:

if label_source == "nmripmap":
    grid = _nmripmap_label_grid(...)
    model_version = "rf-nmripmap-v1"

def _nmripmap_label_grid(...):
    ref_mask = rasterize_mask(fetch_nmripmap(bbox), ...)
    # raises if the tile has no NMRipMap coverage — CO tiles fall back to weak labels
    return replace(grid, label=ref_mask.astype(np.int8))

06The foundation model — a negative result, then the fair test

OlmoEarth-v1-Nano encodes the Sentinel-2 datacube into learned embeddings, which then feed the same random-forest head through the same spatial-CV harness. Identical labels, identical folds — the only variable is the representation.

Getting the encoder to run required resolving a token-mask convention: the model expects a per-band-set mask over a 5-D (batch, H, W, T, band_sets) tensor, where Sentinel-2 L2A is tokenised into three band sets by native resolution. Valid pixels are marked for the online encoder, invalid ones as missing. Patch tokens are then averaged over time and band sets to give one embedding vector per patch.

valid = np.isfinite(s2).all(axis=-1)                        # (1, H, W, T)
mask_pix = np.where(valid, MaskValue.ONLINE_ENCODER.value,
                           MaskValue.MISSING.value)
mask = np.repeat(mask_pix[..., np.newaxis], S2_NUM_BAND_SETS, axis=-1)

output_dict = model.encoder(sample, patch_size=patch_size)
latent, _pooled, _kwargs = unpack_encoder_output(output_dict)

tokens = latent.sentinel2_l2a         # (B, P_H, P_W, T, Band_Sets, D)
emb = tokens.mean(dim=(3, 4)).squeeze(0)   # (P_H, P_W, D)

The AOI is cropped to a square whose side is a multiple of 64, because the 2-D sincos positional encoding requires it. Patch-level labels are assigned by majority vote at ≥ 0.5, the block size for spatial CV shrinks to 0.004° to keep enough folds at patch resolution, and np.kron upsamples patch probabilities back to pixels before reusing the same vectorise-and-write path as the baseline.

Three defects, all found after publishing:

1. The ground truth was ~45% wrong. fetch_nmripmap rasterised every NMRipMap polygon as riparian. On this AOI 341 polygons are mapped but only 189 are woody riparian — the rest are developed, agriculture, upland and water, all of which were being taught to both models as riparian.
2. The foundation model's time axis was averaged away before the classifier ever saw it (the mean-pooling described below).
3. The labels and the imagery were four years apart. NMRipMap v2.0 Plus was photo-interpreted from NAIP 2020; the run used Sentinel-2 2024. Corridors move in four years — defoliation, floods, channel migration — so every model was fed label noise we introduced ourselves.

The re-run below fixes the labels and the pooling, holds everything else constant (same AOI, same patch grid, same spatial folds, same RandomForest head), and varies only the representation:

The baseline won — but this is not a fair test, and comparing against Ai2's own published recipe shows why. Ai2's olmoearth_projects ships a mangrove project that is a near-exact analog of this task — segment a woody vegetation class near water from a Sentinel-2 time series, validated against an authoritative reference map (Global Mangrove Watch there, NMRipMap here). It reports 97.6% accuracy. Its config differs from the run above in four load-bearing ways:

The obvious suspect was the last row. Mean-pooling across the time axis discards the phenology signal — and phenology is the entire thesis of this project (phreatophytes stay green into the dry season because their roots reach groundwater, while the uplands brown out). The model was denied the one thing it exists to exploit, and then judged for failing. That defect is real, and it is now pinned by a unit test that builds two classes with identical seasonal means and different trajectories: under mean-pooling they collapse to the same vector and no classifier can separate them; with the trajectory preserved they separate perfectly.

So the honest statement is narrower than either the original claim or its correction. Every arm above is a frozen encoder feeding a scikit-learn RandomForest — a configuration Ai2 endorses nowhere. What remained genuinely untested was OlmoEarth as Ai2 actually uses it: a fine-tuned BASE with a per-pixel segmentation head, scaffolded at experiments/riparian_extent/. That run has since been done — on a GPU, against the RF on an honest leave-one-reach-out test — and it changes the verdict.

The verdict is GO — the foundation model clears the pre-registered +0.04 macro-mean bar (by +0.074). But the per-fold column tells the real story: on the three river reaches the RF already handled (0.845–0.905) the two tie — the ±0.01–0.03 differences sit inside the metric-definition mismatch above, so "ties on the rivers" is the defensible read, not "RF slightly wins." Every bit of the macro-mean advantage comes from one fold: rescuing the desert arroyo, 0.557 → 0.889 (+0.332) — a single gap larger than the entire net macro gain. OlmoEarth is not a uniformly better model; it is a specialist for the hard, under-represented morphology the context-free per-pixel RF is blind to. So the RF ships for plain extent (no GPU), and the foundation model earns its keep exactly where the training distribution runs out. Full write-up: the LORO result.

07From raster to polygon

Thresholded probability becomes a mask, connected components become polygons, and two cartographic filters run before anything is written: a 500 m² minimum mapping unit that discards single-pixel speckle, and a topology-preserving simplify at three-quarters of a pixel.

if int(region_mask.sum()) * resolution_m ** 2 < MIN_MAPPING_UNIT_M2:
    continue                                    # drop speckle (< 500 m²)

mean_prob = round(float(prob_grid[region_mask].mean()), 4)
poly = shapely_shape(geom).simplify(
    resolution_m * _DEG_PER_M * 0.75, preserve_topology=True)

The write is idempotent by construction — a delete scoped to (method, model_version, huc12) followed by a batch insert — so re-running a tile replaces exactly that tile's output for exactly that model, and never touches another model's rows.

INSERT INTO silver.riparian_extent
    (method, model_version, is_riparian, riparian_probability,
     cell_size_m, huc12, geom)
VALUES (:method, :model_version, :is_riparian, :riparian_probability,
        :cell_size_m, :huc12, ST_SetSRID(ST_GeomFromText(:wkt), 4269))

08The ingest and scoring pipeline

Running alongside the delineation track is the original buffer-centric ETL, which still supplies the vector context layers, the NDVI time series, and the condition score.

Paginated ArcGIS ingest

Every ArcGIS REST layer enforces maxRecordCount = 2000. The client paginates on resultOffset, retries 5xx with exponential backoff, and requests outSR=4269 so nothing is reprojected client-side. The loop terminates on a short page, not only an empty one — a detail that saves one wasted round trip per layer and, more importantly, avoids an infinite loop when a server returns a short final page followed by an empty one.

offset = 0
while True:
    batch = self._client.query(
        url=url, where=where, geometry_filter=envelope,
        result_offset=offset, result_record_count=ARCGIS_BATCH_SIZE,
    )
    if batch.empty:
        break
    batches.append(batch)
    offset += ARCGIS_BATCH_SIZE
    if len(batch) < ARCGIS_BATCH_SIZE:
        break        # short page = last page

The pipeline order is bronze first (watershed boundary loads first because it is the spatial envelope every other query filters against, then NHDPlus, parcels and NWI concurrently on a three-worker pool), then silver spatial processing, then optional enrichment, then gold aggregation. Buffers are generated in a single SQL statement that casts to geography — never buffering in degrees:

INSERT INTO silver.riparian_buffers (stream_id, buffer_distance_m, area_sq_m, geom)
SELECT s.id, :buffer_distance,
    ST_Area(ST_Buffer(s.geom::geography, :buffer_distance)),
    ST_SetSRID(ST_Buffer(s.geom::geography, :buffer_distance)::geometry, 4269)
FROM bronze.streams s

NDVI: scene-first zonal statistics

The naive implementation opens the raster once per buffer. This one inverts the loop — read each scene once, then extract every intersecting buffer's statistics from the in-memory array. Writes are idempotent on (buffer_id, acquisition_date, satellite).

mask = geometry_mask([mapping(geom)], out_shape=ndvi_array.shape,
                     transform=transform, invert=True)
buffer_pixels = ndvi_array[mask]
mean_val, min_val, max_val = compute_ndvi_stats(buffer_pixels)   # valid ∈ [-1, 1]

-- writer
INSERT INTO silver.vegetation_health (buffer_id, acquisition_date, mean_ndvi, ...)
VALUES (...) ON CONFLICT (buffer_id, acquisition_date, satellite) DO NOTHING

Condition scoring

The SMP composite is an 80/10/10 weighting of vegetation structure, connectivity, and contributing area, decomposed into ten sub-scores on a 0–10 scale (NDVI, vertical complexity, species composition, shrub layer, patchiness, native regeneration, native cover) that roll up into an A–F grade. It is computed in Python and materialised into gold.buffer_health_score — the SQL migrations are pure DDL, with no spatial computation hidden in them.

09The data model

Bronze holds raw ingest, silver holds spatial processing, gold holds aggregates. Every geometry column is EPSG:4269, every table carries an audit timestamp, and every spatial table carries a GiST index. Data flows one direction only.

CREATE TABLE IF NOT EXISTS silver.riparian_extent (
    id                   BIGSERIAL PRIMARY KEY,
    method               TEXT          NOT NULL,   -- 'rf' | 'olmoearth'
    model_version        TEXT          NOT NULL,
    is_riparian          BOOLEAN       NOT NULL,
    riparian_probability NUMERIC(5, 4) NOT NULL,
    cell_size_m          NUMERIC(6, 2) NOT NULL,
    geom                 geometry(Polygon, 4269) NOT NULL,
    processed_at         TIMESTAMPTZ   NOT NULL DEFAULT now(),
    CONSTRAINT ck_riparian_extent_method CHECK (method IN ('rf', 'olmoearth')),
    CONSTRAINT ck_riparian_extent_probability
        CHECK (riparian_probability >= 0 AND riparian_probability <= 1)
);
CREATE INDEX IF NOT EXISTS idx_riparian_extent_geom
    ON silver.riparian_extent USING GIST (geom);

Storing method and model_version as columns rather than as separate tables is what makes the RF-versus-foundation-model comparison a query rather than a migration. The same map endpoint serves either by filtering on one column, and a disagreement map between the two is a single self-join away.

10Serving: Aspire, Dapper, and vector tiles

The C# side is a thin, strictly layered read API over PostGIS. Endpoints hold no SQL, services hold no data access, and the repository holds no business logic.

Orchestration

.NET Aspire wires the four services together and injects connection strings by service discovery. Two details are load-bearing for a project living on an external drive: the PostGIS container is bind-mounted rather than using a Docker volume (volumes are lost when Docker restarts on removable media), and the ETL defaults to incremental so an orchestrator restart cannot silently wipe the NDVI table.

var ripariandb = builder.AddPostgres("postgres")
    .WithImage("postgis/postgis").WithImageTag("16-3.4")
    .WithLifetime(ContainerLifetime.Persistent)
    .WithDataBindMount("../pgdata")      // survives Docker restarts
    .AddDatabase("ripariandb");

var api = builder.AddProject<Projects.RiparianPoc_Api>("api")
    .WithReference(ripariandb).WaitFor(ripariandb)
    .WithExternalHttpEndpoints();

var etl = builder.AddDockerfile("etl", "../python-etl")
    .WithReference(ripariandb).WaitFor(ripariandb)
    // incremental by default — a restart must not truncate NDVI
    .WithEnvironment("ETL_MODE", builder.Configuration["ETL_MODE"] ?? "incremental");

The three layers

Interfaces are segregated by return shape, not by table: ISpatialQueryService returns GeoJSON feature collections and MVT byte arrays; IComplianceDataService returns typed records. A handler does one thing.

// endpoint: inject, call, return. No SQL here, ever.
private static async Task<IResult> GetBuffersWithHealth(
    ISpatialQueryService spatialService, CancellationToken ct)
{
    var fc = await spatialService.GetBuffersWithHealthAsync(ct);
    return TypedResults.Ok(fc);
}

// service: owns the SQL and the span
using var activity = Source.StartActivity("SpatialQuery.GetBuffersWithHealth");
const string sql = """
    SELECT rb.id, rb.stream_id, rb.buffer_distance_m, rb.area_sq_m,
           s.gnis_name AS stream_name,
           vh.mean_ndvi, vh.health_category, vh.acquisition_date,
           ST_AsGeoJSON(rb.geom) AS geojson
    FROM silver.riparian_buffers rb
    JOIN bronze.streams s ON s.id = rb.stream_id
    LEFT JOIN LATERAL (
        SELECT mean_ndvi, health_category, acquisition_date
        FROM silver.vegetation_health
        WHERE buffer_id = rb.id AND season_context = 'peak_growing'
        ORDER BY acquisition_date DESC
        LIMIT 1
    ) vh ON TRUE
    """;
var fc = await _repository.QueryGeoJsonAsync(sql, null, ct);
activity?.SetTag(FeatureCountTag, fc.Count);

The LEFT JOIN LATERAL … LIMIT 1 is the right tool here: it fetches the most recent peak-growing reading per buffer in one pass, keeps buffers with no reading at all (which render in the default colour rather than vanishing), and avoids both a correlated subquery per row and a window function over the entire history table.

Vector tiles with an index-backed pre-filter

Every tile query follows one pattern. ST_TileEnvelope builds the tile bounds in Web Mercator; the join condition transforms that envelope back to 4269 and uses the && bounding-box operator, which is what lets PostgreSQL hit the GiST index before doing any expensive geometry work. Reprojection happens only on the handful of rows that survive.

WITH tile AS (
    SELECT ST_TileEnvelope(@z, @x, @y) AS envelope
),
mvt_geom AS (
    SELECT b.id,
           ST_AsMVTGeom(ST_Transform(b.geom, 3857), tile.envelope) AS geom,
           b.buffer_distance_m,
           COALESCE(h.score_grade, 'Unknown') as grade,
           h.composite_score
    FROM silver.riparian_buffers b
    JOIN tile ON b.geom && ST_Transform(tile.envelope, 4269)  -- GiST pre-filter
    LEFT JOIN LATERAL (
        SELECT score_grade, composite_score
        FROM gold.buffer_health_score s
        WHERE s.buffer_id = b.id ORDER BY s.id DESC LIMIT 1
    ) h ON true
)
SELECT ST_AsMVT(mvt_geom.*, 'buffers', 4096, 'geom') FROM mvt_geom;

That change alone made tile queries 10–40× faster. The repository returns the protobuf blob straight from ExecuteScalarAsync<byte[]> — the tile is encoded in the database and never materialised as .NET geometry objects.

11The map

React 18 with react-map-gl/maplibre over MapLibre GL 3.6.2, Vite 6, Tailwind 3. Layers are vector tiles wherever volume justifies it and GeoJSON where it does not. Styling is pushed into MapLibre paint expressions so the GPU does the work, not React.

Categorical data uses match; continuous data uses interpolate. Model probability is continuous, so riparian extent ramps rather than steps:

<Source id="riparian-extent-source" type="geojson" data={riparianExtent ?? EMPTY_FC}>
  <Layer id="riparian-extent-fill" type="fill"
    layout={{ visibility: showRiparianExtent ? 'visible' : 'none' }}
    paint={{
      'fill-color': [
        'interpolate', ['linear'], ['get', 'riparian_probability'],
        0.5,  '#a7f3d0',
        0.75, '#34d399',
        1.0,  '#059669',
      ],
      'fill-opacity': 0.6,
      'fill-outline-color': '#047857',
    }} />
</Source>

The timelapse works by rebuilding the tile URL template when the selected date changes, and keying the <Source> on that date so MapLibre remounts and refetches. Layer visibility is toggled through the layout.visibility property rather than by unmounting components — the tiles stay in MapLibre's cache, so toggling a layer back on is instant.

const ndviTileUrl = useMemo(() => {
  if (selectedDate)
    return `${API_URL}/api/tiles/buffers-ndvi/${selectedDate}/{z}/{x}/{y}.pbf`;
  return `${API_URL}/api/tiles/buffers-ndvi/{z}/{x}/{y}.pbf`;
}, [selectedDate]);

// key forces a source remount when the date changes
<Source key={`ndvi-tiles-${selectedDate ?? 'latest'}`}
        id="buffer-ndvi-source" type="vector" tiles={[ndviTileUrl]}>

Session correlation is end-to-end. A UUID is minted per page load and attached to every request — including MapLibre's own tile fetches, via transformRequest, which is easy to forget and would otherwise leave the tile traffic untraceable.

const SESSION_ID = crypto.randomUUID();

async function fetchJson<T>(url: string): Promise<T> {
  const res = await fetch(url, { headers: { 'X-Session-Id': SESSION_ID } });
  const correlationId = res.headers.get('X-Correlation-Id');
  if (correlationId) console.debug(`[API] ${url} correlation=${correlationId}`);
  if (!res.ok) { /* parse ApiErrorResponse, surface correlationId */ }
  return res.json() as Promise<T>;
}

// MapLibre tile requests carry the session too
const transformRequest = useCallback((url: string) =>
  url.includes('/api/') ? { url, headers: { 'X-Session-Id': SESSION_ID } } : { url }, []);

12Observability and failure behaviour

Three custom ActivitySources are registered additively on top of Aspire's defaults, giving a trace hierarchy of HTTP request → service span → repository span → Npgsql span. Every repository call tags its duration and row count; a failed query marks the span red before rethrowing.

builder.Services.AddOpenTelemetry().WithTracing(t => t
    .AddSource("RiparianPoc.Api.Repository")
    .AddSource("RiparianPoc.Api.SpatialQuery")
    .AddSource("RiparianPoc.Api.ComplianceData"));

// one place decides how a failure becomes a status code
var (statusCode, message) = exception switch
{
    NpgsqlException            => (503, "Database temporarily unavailable"),
    _ when exception.InnerException is NpgsqlException
                               => (503, "Database temporarily unavailable"),
    ArgumentException ex       => (400, ex.Message),      // 4xx: safe to echo
    KeyNotFoundException ex    => (404, ex.Message),
    OperationCanceledException => (504, "Request timed out"),
    _                          => (500, "An unexpected error occurred"),  // 5xx: never echo
};

Client messages are asymmetric on purpose: 4xx responses echo the exception message because the caller caused it and can fix it; 5xx responses return a generic string, with the exception detail exposed only in Development. Every response carries an X-Correlation-Id, which CORS explicitly exposes so the browser can read it — the reason a user-reported error can be traced to one span in the Aspire dashboard.

Where a reviewer should attack this

Collected in one place, ordered by how much they should change your confidence in the result. The first three are science; the rest are engineering.

Worth saying plainly: the reason this list can be written at all is that the pipeline was built to produce falsifiable numbers. Spatial cross-validation against an independent reference is what surfaced the weak-label failure, and it is what made the foundation model's loss visible instead of flattering.

What building this actually required

Mapped against the four qualification bullets this document was written to test. Three hold up; one does not, and two large areas are missing from the list entirely.

Read as written, the original four bullets describe someone who could build the map and wire up the ETL. They would not describe the person who wrote validate.py — and that file is the reason any of the numbers above can be believed.