From a raw radar archive to the flood map on the case page.

"""run_banda_aceh_pipeline.py — stage dispatcher.

One entry point. One env var (ASIA_FLOOD_BASE_DIR). Each stage lives in
its own file so students can run them in isolation and inspect outputs.
"""
import os, sys, subprocess
from pathlib import Path

BASE_DIR = Path(os.environ["ASIA_FLOOD_BASE_DIR"])

STEP_TO_SCRIPT: dict[str, str] = {
    "preprocess":     "preprocess/snap_preprocess_banda_aceh.py",
    "build-dataset":  "dataset/prepare_dataset_from_three_tifs.py",
    "stats":          "infer/calculate_banda_aceh_stats.py",
    "predict":        "infer/predict_banda_aceh_adapted.py",
    "validate":       "validate/validate_predictions.py",
}

# Shorthand macro: build-dataset → predict → validate, without SNAP.
REPRODUCE_NO_SAR: list[str] = ["build-dataset", "predict", "validate"]


def run_step(step: str) -> int:
    """Exec one pipeline stage in a subprocess, propagating env and cwd."""
    ...


def main() -> None:
    """CLI: 'python run_banda_aceh_pipeline.py <step> | reproduce-no-sar'"""
    ...

"""snap_preprocess_banda_aceh.py — raw SAFE → calibrated GeoTIFF.

GPT graph: Apply-Orbit-File → Calibration → Speckle-Filter →
Range-Doppler Terrain-Correction → Subset, with an external SRTM tile
as DEM. A second gdalwarp pass trims black edges.
"""
from pathlib import Path
import subprocess

SNAP_HOME   = Path("/home/yang/snap")
GRAPH_FILE  = "preprocess/grd_preprocessing_external_dem.xml"
DEM_TILE    = "assets/dem/N05E095.tif"
AOI_WKT     = "POLYGON((95.15 5.35, 95.50 5.35, 95.50 5.70, 95.15 5.70, 95.15 5.35))"
FINAL_AOI   = {"lon_min": 95.25, "lon_max": 95.40, "lat_min": 5.45, "lat_max": 5.60}
DATES       = ["20251021", "20251102", "20251126"]


def run_snap(date: str) -> Path:
    """GPT invocation with -PexternalDEMFile pinned for reproducibility."""
    ...


def tighten_with_gdalwarp(src: Path) -> Path:
    """gdalwarp -te <FINAL_AOI> -r bilinear -co COMPRESS=LZW -co TILED=YES."""
    ...


def main() -> None:
    """Run SNAP + gdalwarp for all three acquisition dates."""
    ...

"""prepare_dataset_from_three_tifs.py — 3 scenes → N tiles.

Deterministic non-overlapping 224² grid. Per tile: writes three TIFs
(MS1 = post, SL1 = mid, SL2 = baseline) plus info.json. A sibling script
builds the pickle index that the Dataset class later looks up at load
time.
"""
import math, json, pickle
from pathlib import Path
import rasterio

PATCH_SIZE     = 224
STRIDE         = 224                 # no overlap, deterministic grid
ACT_ID, AOI_ID = 999, 1               # banda_aceh as a custom 'event'


def slice_scene(scene_tif: Path, out_root: Path, role: str) -> list[dict]:
    """Walk the patch grid, write <role>.tif per cell, record info."""
    ...


def build_pickle(records: list[dict], out: Path) -> None:
    """Pickle {(act_id, aoi_id): [{info, path}, ...]} for Dataset lookup."""
    ...


def main() -> None:
    """Run slice_scene on pre-2, pre-1, post. Write pickle index last."""
    ...

"""calculate_banda_aceh_stats.py — per-region, per-clamp normalization.

Input : 3 × processed TIFs (VV, VH bands, > 0 pixels only).
Output: configs/banda_aceh_adapted_configs.json
          { clamp015: { data_mean, data_std, clamp_input }, clamp03: ..., ... }
The four entries (+ 'original' = KuroSiwo defaults) are the four
configs the inference step loops over.
"""
import json
from pathlib import Path
import numpy as np
import rasterio

CLAMPS = [0.15, 0.30, 0.50]


def iter_pixels(tif: Path) -> np.ndarray:
    """Read VV & VH, drop NoData (<= 0), return (N, 2) array."""
    ...


def stats_for_clamp(x: np.ndarray, c: float) -> dict:
    """Per-channel mean/std after clipping, plus truncation ratio."""
    x_c    = np.minimum(x, c)
    trunc  = (x > c).mean(axis=0)
    mean   = x_c.mean(axis=0)
    std    = x_c.std(axis=0)
    return {"clamp_input": c, "data_mean": mean.tolist(),
            "data_std": std.tolist(), "truncated_pct": trunc.tolist()}


def main() -> None:
    pixels  = np.vstack([iter_pixels(p) for p in sorted(TIFS)])
    configs = {f"clamp{int(c * 100):02d}": stats_for_clamp(pixels, c)
               for c in CLAMPS}
    configs["original"] = KUROSIWO_DEFAULTS        # training-time stats
    Path("configs/banda_aceh_adapted_configs.json").write_text(
        json.dumps(configs, indent=2))

"""predict_banda_aceh_adapted.py — 4 configs × one scene.

For each of {original, clamp015, clamp03, clamp05}:
  1. Rebuild KuroSiwoDataset with that config's (clamp, mean, std).
  2. Loop 224² patches, forward, argmax → 3-class mask.
  3. Paste each prediction back at its (row, col) → 2.8-Mpx full map.
  4. Write GeoTIFF + stats.json into prediction_results_adapted_<cfg>/.
"""
import json
from pathlib import Path
import numpy as np
import rasterio
import torch
from torch.utils.data import DataLoader

from models.csmamba           import CSMamba
from dataset.kurosiwo_dataset import KuroSiwoDataset

ADAPTED_CONFIGS = json.loads(
    Path("configs/banda_aceh_adapted_configs.json").read_text()
)


def build_model(device: str) -> CSMamba:
    model = CSMamba(
        img_size=224, in_channels=6, num_classes=3,
        embed_dims=[96, 192, 384, 768], depths=[1, 1, 6, 1], d_state=16,
    ).to(device).eval()
    ckpt = torch.load("assets/checkpoints/CSMamba_FloodFocus_best_model.pt")
    model.load_state_dict(ckpt["model_state_dict"])
    return model


def reassemble(preds, row_idx, col_idx, H_full: int, W_full: int) -> np.ndarray:
    """Paste (B, 224, 224) predictions back into a (H_full, W_full) map."""
    ...


def main() -> None:
    device = "cuda" if torch.cuda.is_available() else "cpu"
    model  = build_model(device)

    for cfg_key in ("original", "clamp015", "clamp03", "clamp05"):
        cfg = ADAPTED_CONFIGS[cfg_key]
        ds  = KuroSiwoDataset(mode="test", configs=cfg, pickle_path=PKL)
        dl  = DataLoader(ds, batch_size=8, num_workers=4, shuffle=False)
        full = np.zeros((H_full, W_full), dtype=np.uint8)

        with torch.no_grad():
            for image, pre1, pre2, _, _, rr, cc in dl:
                x = torch.cat([pre2, pre1, image], dim=1).to(device)      # (B, 6, 224, 224)
                y = model(x).argmax(1).cpu().numpy()                       # (B, 224, 224)
                reassemble_inplace(full, y, rr, cc)

        write_geotiff(full, out_dir / "flood_prediction.tif")
        dump_stats(full, out_dir / "prediction_stats.json")

"""validate_predictions.py — cross-config triangulation.

Inputs : 4 × flood_prediction.tif (one per config).
Outputs: validation_report.json + prediction_comparison.png + 6 diff PNGs.
"""
from itertools import combinations
from pathlib import Path
import json
import numpy as np
import rasterio
from scipy import ndimage


def per_config_diagnostics(pred: np.ndarray) -> dict:
    """Pixel counts per class, flood connected components, boundary ratio."""
    labeled, n = ndimage.label(pred == 2)
    return {"flood_regions": int(n), "flood_pct": float((pred == 2).mean())}


def pairwise_agreement(preds: dict[str, np.ndarray]) -> dict:
    """For every unordered pair, (pred_a == pred_b).mean()."""
    return {f"{a}_vs_{b}": float((preds[a] == preds[b]).mean())
            for a, b in combinations(preds.keys(), 2)}


def render_comparison(preds, vv, vh, out: Path) -> None:
    """5 rows × 3 cols: reference + 4 configs, side-by-side with SAR."""
    ...


def main() -> None:
    preds   = {k: rasterio.open(TIF[k]).read(1) for k in CONFIGS}
    report  = {k: per_config_diagnostics(p) for k, p in preds.items()}
    report |= pairwise_agreement(preds)
    Path("validation_report.json").write_text(json.dumps(report, indent=2))
    render_comparison(preds, vv, vh, Path("prediction_comparison.png"))