#!/usr/bin/env python3
"""Build the bounded browser presentation of the admitted SkyClimate V2.1 atlas.

The authoritative V2.1 atlas stays immutable on R2.  This builder verifies the
admission, manifest and GeoParquet identities before deriving the 45 overall
map layers that are actually present in ``map_summary.parquet``.  It deliberately
does not invent map layers for the 325 conditioned views held by the atlas tile
shards and bounded API.
"""

from __future__ import annotations

import argparse
import gzip
import hashlib
import io
import json
import runpy
import shutil
import subprocess
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path
from typing import Any

import numpy as np
import pyarrow as pa
import pyarrow.parquet as pq

WEB_MERCATOR_LATITUDE_LIMIT = 85.0511287798066
WEB_MERCATOR_MAX_CELL_CENTRE_LAT = 84.75
EXPECTED_PACK_ID = "scga21-4cd927cac60f0a63c7479ad447c69a62bd0f4bae0e8b25df34f6ddee02f57452"
EXPECTED_MAP_SHA256 = "707648c5e62114deb6de6a3041965e95289b82075ab3ac170533b89aeaece29a"
V21_POLICY_ID = "scv21dp1-a7a8aa811cbe404044d6f328a2416febd7b71c5c2cc58a70374367cc091b4b04"
AEROSOL_RELEASE_ROOT = (
    "skyclimate/aerosol-sketch/v2.1/release/"
    "sha256:d90a4531594149e46869cd6d54f7e125a4a28c062c1dc764220b4c15b1252442/"
    "scasp21-9b05685159b13db0136baad07b668800e17568ed445248262216feb5e416c6b1/"
)
AEROSOL_RECEIPT_SHA256 = "1087b5782994a6cca9f571fa1c8be5629fd48cb884a14e1e6273bc9c83965bfc"
AEROSOL_RELEASE_SHA256 = "5ca3ffbed98e2afa2e2c9bf9421f3c914e26f0be51b0c00d17c609fa39edd475"
V4_SUMMARY_ID = "scl137ms1-f8e15c56d9880a9d9ce5b1b391a18f00d87d87cbcb1281ec1bf91874ee78ada9"
V4_REPORT_SHA256 = "ac51ae3da764b6b77098422b0ef81a1c87b8b333693624324f28c0fc4409ec03"
V4_REPORT_KEY = f"skyclimate/turbulence/pressure-level-phase-f/v1/l137-reference-summary-reports/2018-01/{V4_REPORT_SHA256}.json"
V4_FREEZE_SHA256 = "3c341252f9e3f3bf753f1db5c6fca0f235792b20e706e21325df941643f97357"
V4_FREEZE_KEY = (
    "migration-metadata/skyclimate-era5-l137-v4-reference-freeze/"
    "skyclimate-era5-l137-v4-january-reference-2018-01/"
    "c1dd3a6024514a3cfe940aed3ca0e5c5cef720650063dec4df8dff98af9b9198/"
    "era5-l137-v4-reference-freeze-2018-01_2018-03.json"
)
V4_METRICS = {
    "cn2_integral_m_one_third": ("cn2-integral", "Integrated optical turbulence", "m^(1/3)", 1e-12),
    "r0_1550_m": ("r0-1550nm", "Fried parameter · 1550 nm", "m", 1.0),
    "theta0_1550_rad": ("theta0-1550nm", "Isoplanatic angle · 1550 nm", "rad", 1e-4),
    "tau0_1550_s": ("tau0-1550nm", "Coherence time · 1550 nm", "s", 0.05),
    "reference_rytov_variance_1550": ("rytov-1550nm", "Reference Rytov variance · 1550 nm", "1", 1.0),
    "pwv_above_model_surface_kg_m2": ("pwv-model-surface", "Water vapour above model surface", "kg m-2", 60.0),
}
AEROSOL_STATISTICS = {
    "mean": ("mean", "Mean"), "std": ("std", "Population standard deviation"),
    "min": ("min", "Minimum"), "max": ("max", "Maximum"),
    "p05": ("q05", "P05"), "p10": ("q10", "P10"), "p50": ("q50", "Median"),
    "p90": ("q90", "P90"), "p95": ("q95", "P95"),
    "coverage": ("valid_fraction", "Valid three-hourly fraction"),
    "count": ("valid_count", "Valid three-hourly samples"),
}

PALETTES = {
    "cloud": [[10, 23, 38], [21, 64, 78], [31, 120, 126], [127, 193, 176], [239, 244, 222]],
    "water": [[11, 28, 45], [18, 79, 114], [28, 139, 145], [132, 196, 165], [241, 238, 190]],
    "coverage": [[22, 25, 29], [55, 73, 89], [63, 122, 139], [120, 181, 168], [229, 240, 222]],
    "aerosol": [[16, 33, 50], [47, 111, 135], [150, 182, 130], [224, 161, 66], [168, 47, 53]],
}

VARIABLES = {
    "tcc": ("Total cloud cover", "fraction", "total atmospheric column", "cloud"),
    "lcc": ("Low cloud cover", "fraction", "low atmospheric layer", "cloud"),
    "mcc": ("Middle cloud cover", "fraction", "middle atmospheric layer", "cloud"),
    "hcc": ("High cloud cover", "fraction", "high atmospheric layer", "cloud"),
    "tcwv": ("Total column water vapour", "kg m-2", "total atmospheric column", "water"),
}

STATISTICS = {
    "mean": "Mean",
    "std": "Population standard deviation",
    "p05": "P05",
    "p25": "P25",
    "p50": "Median",
    "p75": "P75",
    "p95": "P95",
    "coverage": "Valid hourly fraction",
    "count": "Valid hourly samples",
}


def arguments() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--map-summary", type=Path)
    parser.add_argument("--admission", type=Path)
    parser.add_argument("--manifest", type=Path)
    parser.add_argument("--release-catalogue", type=Path)
    parser.add_argument("--append-aerosol", action="store_true",
                        help="Read existing pinned R2 statistics and append their presentation; no production")
    parser.add_argument("--append-v4-demo", action="store_true",
                        help="Present the retained January V4 summary, without producing or admitting new turbulence")
    parser.add_argument("--package-existing", type=Path,
                        help="Assemble a public release around previously staged native-data indexes; no data production")
    parser.add_argument("--refresh-unpublished-staging", action="store_true",
                        help="Refresh this task's local package metadata before any publication; never use after upload")
    parser.add_argument("--base-catalogue", type=Path, required=True)
    parser.add_argument("--output-root", type=Path, required=True)
    args = parser.parse_args()
    if sum(bool(x) for x in (args.append_aerosol, args.append_v4_demo, args.package_existing)) > 1:
        parser.error("append one existing family at a time")
    if not (args.append_aerosol or args.append_v4_demo or args.package_existing) and any(
        getattr(args, name) is None for name in ("map_summary", "admission", "manifest", "release_catalogue")
    ):
        parser.error("cloud presentation requires --map-summary, --admission, --manifest and --release-catalogue")
    return args


def sha256(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as stream:
        for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def read_json(path: Path) -> dict[str, Any]:
    value = json.loads(path.read_text(encoding="utf-8"))
    if not isinstance(value, dict):
        raise TypeError(f"{path} is not a JSON object")
    return value


def verify_inputs(args: argparse.Namespace) -> tuple[dict[str, Any], dict[str, Any], dict[str, Any]]:
    admission = read_json(args.admission)
    manifest = read_json(args.manifest)
    release_catalogue = read_json(args.release_catalogue)
    if admission.get("status") != "succeeded" or admission.get("pack_id") != EXPECTED_PACK_ID:
        raise ValueError("the atlas admission is not the expected succeeded V2.1 pack")
    if admission.get("provider_readback_verified") is not True:
        raise ValueError("the atlas admission lacks provider read-back verification")
    if manifest.get("pack_id") != EXPECTED_PACK_ID or manifest.get("claim_class") != "source-grid-statistic":
        raise ValueError("the atlas manifest identity or claim class differs")
    if manifest.get("period") != {
        "cadence": "PT1H",
        "end_exclusive": "2026-01-01T00:00:00Z",
        "month_count": 96,
        "start": "2018-01-01T00:00:00Z",
        "valid_time_count": 70128,
    }:
        raise ValueError("the atlas period is not the complete hourly 2018-2025 axis")
    if manifest.get("grid", {}).get("shape") != [721, 1440] or manifest.get("grid", {}).get("source_resolution_degrees") != 0.25:
        raise ValueError("the atlas source grid differs")
    map_record = manifest.get("artifacts", {}).get("map_summary.parquet", {})
    if map_record.get("sha256") != f"sha256:{EXPECTED_MAP_SHA256}" or sha256(args.map_summary) != EXPECTED_MAP_SHA256:
        raise ValueError("map_summary.parquet does not match the admitted atlas")
    catalogue_product = next(
        (item for item in release_catalogue.get("products", []) if item.get("product_id") == "global-atlas-pack"),
        None,
    )
    if not catalogue_product or EXPECTED_PACK_ID not in catalogue_product.get("validation", {}).get("evidence_ids", []):
        raise ValueError("the V2.1 release catalogue does not bind the admitted atlas")
    if catalogue_product.get("rights", {}).get("policy_id") != V21_POLICY_ID:
        raise ValueError("the atlas delivery policy differs")
    return admission, manifest, catalogue_product


def colourise(values: np.ndarray, low: float, high: float, palette: str) -> np.ndarray:
    valid = np.isfinite(values)
    scaled = np.where(valid, np.clip((values - low) / (high - low), 0.0, 1.0), 0.0)
    stops = np.linspace(0.0, 1.0, len(PALETTES[palette]))
    rgba = np.zeros((*values.shape, 4), dtype=np.uint8)
    for channel in range(3):
        rgba[..., channel] = np.interp(
            scaled, stops, [colour[channel] for colour in PALETTES[palette]]
        ).astype(np.uint8)
    rgba[..., 3] = np.where(valid, 224, 0).astype(np.uint8)
    return rgba


def web_mercator_y(latitude_deg: np.ndarray | float) -> np.ndarray:
    latitude = np.clip(np.asarray(latitude_deg, dtype=np.float64), -WEB_MERCATOR_LATITUDE_LIMIT, WEB_MERCATOR_LATITUDE_LIMIT)
    return (1.0 - np.arcsinh(np.tan(np.deg2rad(latitude))) / np.pi) / 2.0


def inverse_web_mercator_y(y: np.ndarray) -> np.ndarray:
    return np.rad2deg(np.arctan(np.sinh(np.pi * (1.0 - 2.0 * y))))


def project(values: np.ndarray, latitudes: np.ndarray, longitudes: np.ndarray) -> tuple[np.ndarray, list[list[float]]]:
    north = WEB_MERCATOR_LATITUDE_LIMIT
    south = -WEB_MERCATOR_LATITUDE_LIMIT
    image_width = len(longitudes) + 1
    x_pixel_span = 1.0 / image_width
    north_y = float(web_mercator_y(north))
    south_y = float(web_mercator_y(south))
    image_height = int(np.ceil((south_y - north_y) / x_pixel_span))
    target_longitudes = -180.0 + (np.arange(image_width, dtype=np.float64) + 0.5) * (360.0 / image_width)
    target_y = north_y + (np.arange(image_height, dtype=np.float64) + 0.5) * ((south_y - north_y) / image_height)
    target_latitudes = inverse_web_mercator_y(target_y)
    longitude_step = float(longitudes[1] - longitudes[0])
    latitude_step = float(latitudes[0] - latitudes[1])
    if longitude_step <= 0 or latitude_step <= 0 or not np.allclose(np.diff(longitudes), longitude_step) or not np.allclose(np.diff(latitudes), -latitude_step):
        raise ValueError("presentation requires a regular, ordered source grid")
    x_indices = np.floor((target_longitudes - longitudes[0]) / longitude_step + 0.5).astype(np.int64) % len(longitudes)
    y_indices = np.clip(np.floor((latitudes[0] - target_latitudes) / latitude_step + 0.5).astype(np.int64), 0, len(latitudes) - 1)
    return values[np.ix_(y_indices, x_indices)], [[-180.0, north], [180.0, north], [180.0, south], [-180.0, south]]


def visual_domain(variable: str, statistic: str) -> tuple[float, float, str]:
    if statistic == "coverage":
        return 0.0, 1.0, "coverage"
    if statistic == "count":
        return 0.0, 70128.0, "coverage"
    if variable != "tcwv":
        return (0.0, 0.5, "cloud") if statistic == "std" else (0.0, 1.0, "cloud")
    return (0.0, 30.0, "water") if statistic == "std" else (0.0, 60.0, "water")


def write_assets(output_root: Path, slug: str, values: np.ndarray, display: np.ndarray, visual_min: float, visual_max: float, palette: str) -> dict[str, Any]:
    from PIL import Image

    finite = values[np.isfinite(values)]
    if finite.size == 0:
        raise ValueError(f"{slug} has no finite values")
    storage_min = float(finite.min())
    storage_max = float(finite.max())
    span = storage_max - storage_min
    encoded = np.full(values.shape, 65535, dtype="<u2")
    valid = np.isfinite(values)
    encoded[valid] = 0 if span == 0.0 else np.rint(np.clip((values[valid] - storage_min) / span, 0.0, 1.0) * 65534.0).astype("<u2")
    layer_path = output_root / "v21" / "layers" / f"{slug}.png"
    value_path = output_root / "v21" / "values" / f"{slug}.u16.bin"
    layer_path.parent.mkdir(parents=True, exist_ok=True)
    value_path.parent.mkdir(parents=True, exist_ok=True)
    Image.fromarray(colourise(display, visual_min, visual_max, palette)).save(layer_path, optimize=True)
    with value_path.open("wb") as raw, gzip.GzipFile(filename="", mode="wb", fileobj=raw, compresslevel=9, mtime=0) as stream:
        stream.write(encoded.tobytes(order="C"))
    max_error = 0.0 if span == 0.0 else span / (2.0 * 65534.0)
    return {
        "imageUrl": f"/data/v21/layers/{layer_path.name}",
        "valueUrl": f"/data/v21/values/{value_path.name}",
        "valueEncoding": "gzip",
        "storageMin": storage_min,
        "storageMax": storage_max,
        "nodataCode": 65535,
        "maximumQuantisationError": max_error,
        "visualMin": visual_min,
        "visualMax": visual_max,
        "palette": PALETTES[palette],
    }


def build_variants(map_path: Path, output_root: Path) -> list[dict[str, Any]]:
    parquet = pq.ParquetFile(map_path)
    required = ["source_latitude_index", "source_longitude_index"] + [
        f"{variable}_{'valid_fraction' if statistic == 'coverage' else 'valid_count' if statistic == 'count' else statistic}"
        for variable in VARIABLES
        for statistic in STATISTICS
    ]
    missing = sorted(set(required) - set(parquet.schema_arrow.names))
    if missing:
        raise ValueError(f"map summary lacks columns: {missing}")
    table = parquet.read(columns=required)
    latitude_index = table["source_latitude_index"].to_numpy()
    longitude_index = table["source_longitude_index"].to_numpy()
    if table.num_rows != 1_038_240 or len(np.unique(latitude_index * 1440 + longitude_index)) != table.num_rows:
        raise ValueError("map summary does not contain one row per global source pixel")
    source_latitudes = np.linspace(90.0, -90.0, 721, dtype=np.float64)
    source_longitudes = np.arange(1440, dtype=np.float64) * 0.25
    latitude_keep = np.flatnonzero(np.abs(source_latitudes) <= WEB_MERCATOR_MAX_CELL_CENTRE_LAT)
    normalised_longitudes = ((source_longitudes + 180.0) % 360.0) - 180.0
    longitude_order = np.argsort(normalised_longitudes)
    latitudes = source_latitudes[latitude_keep]
    longitudes = normalised_longitudes[longitude_order]
    variants: list[dict[str, Any]] = []
    for variable, (label, base_unit, vertical_support, _) in VARIABLES.items():
        for statistic, statistic_label in STATISTICS.items():
            column = f"{variable}_{'valid_fraction' if statistic == 'coverage' else 'valid_count' if statistic == 'count' else statistic}"
            full = np.full((721, 1440), np.nan, dtype=np.float32)
            full[latitude_index, longitude_index] = table[column].to_numpy(zero_copy_only=False).astype(np.float32)
            values = full[np.ix_(latitude_keep, longitude_order)]
            display, coordinates = project(values, latitudes, longitudes)
            visual_min, visual_max, palette = visual_domain(variable, statistic)
            slug = f"scga21-4cd927-{variable}-{statistic}-overall"
            unit = "fraction" if statistic == "coverage" else "hours" if statistic == "count" else base_unit
            variants.append({
                "id": slug,
                "productId": EXPECTED_PACK_ID,
                "variableId": variable,
                "variableLabel": label,
                "statisticId": statistic,
                "statisticLabel": statistic_label,
                "temporalKind": "overall",
                "temporalValue": None,
                "unit": unit,
                "verticalSupport": vertical_support,
                "grid": {
                    "width": 1440,
                    "height": len(latitudes),
                    "north": float(latitudes[0]),
                    "south": float(latitudes[-1]),
                    "west": float(longitudes[0]),
                    "east": float(longitudes[-1]),
                    "latitudeStep": 0.25,
                    "longitudeStep": 0.25,
                    "longitudeCyclic": True,
                    "imageProjection": "web-mercator-nearest",
                    "imageWidth": int(display.shape[1]),
                    "imageHeight": int(display.shape[0]),
                    "imageCoordinates": coordinates,
                },
                **write_assets(output_root, slug, values, display, visual_min, visual_max, palette),
            })
    return variants


def build_product(
    manifest: dict[str, Any],
    catalogue_product: dict[str, Any],
    *,
    manifest_sha256: str,
) -> dict[str, Any]:
    return {
        "productId": EXPECTED_PACK_ID,
        "label": "SkyClimate Global Atlas V2.1",
        "source": "ERA5 hourly cloud layers and total column water vapour",
        "sourceKind": "reanalysis",
        "claimClass": "source-grid-statistic",
        "periodLabel": "2018–2025 complete hourly reference",
        "periodStart": "2018-01-01T00:00:00Z",
        "periodEndExclusive": "2026-01-01T00:00:00Z",
        "cadence": "PT1H · 70,128 hourly instants",
        "resolutionDeg": 0.25,
        # ``rights["grant"]`` records what the upstream Copernicus licence
        # permits Dynamik Orbits to do.  It is not the distribution licence of
        # the combined SkyClimate product.
        "licence": "CC-BY-NC-4.0",
        "deliveryState": "permitted-with-required-notices",
        "sourceUrl": "https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels",
        "datasetDoi": "10.24381/cds.adbb2d47",
        "licenceUrl": "https://creativecommons.org/licenses/by-nc/4.0/",
        "sourceLicenceUrl": "https://creativecommons.org/licenses/by/4.0/",
        "attribution": "Generated using Copernicus Climate Change Service information [2026]",
        "modificationIndication": "Contains modified Copernicus Climate Change Service information [2026]; SkyClimate computes exact source-grid statistics over the complete 2018–2025 hourly axis without spatial interpolation.",
        "disclaimer": "Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data this product contains.",
        "onDemandStudyRequestSha256": None,
        "sourceReleasePackId": None,
        "sourceLineage": None,
        "sourcePackIntegrityVerified": True,
        "validationSourceIds": [],
        "uncertainty": {
            "state": "unavailable",
            "method": None,
            "nominalCoverage": None,
            "reason": "V2.1 engineering admission verifies completeness and source-grid fidelity; it does not provide calibrated local uncertainty.",
        },
        "scientificValidation": None,
        "exactSite": None,
        "siteDossier": None,
        "manifestSha256": f"sha256:{manifest_sha256}",
        "statisticsSha256": manifest["artifacts"]["statistics.nc"]["sha256"],
    }


def fetch_aerosol_summary(output_root: Path) -> tuple[dict[str, Any], Path, str, dict[str, Any]]:
    """Project columns from existing, hash-verified tiles; do not derive AOD."""
    from skyclimate.identity import canonical_json_bytes
    from skyclimate.runner import R2Target

    fetch_runner_publication_metadata = runpy.run_path(
        str(Path(__file__).with_name("fetch-v21-atlas-presentation-source.py"))
    )["fetch_runner_publication_metadata"]

    target = R2Target()
    client = target.client()
    parent = fetch_runner_publication_metadata(
        client, bucket=target.bucket, receipt_key=AEROSOL_RELEASE_ROOT + "receipt.json",
        receipt_sha256=AEROSOL_RECEIPT_SHA256, kernel="aerosol-sketch-release",
        identity_id=AEROSOL_RELEASE_ROOT.rstrip("/").split("/")[-2],
        names=("manifest.json", "release.json"),
    )
    release = parent["documents"]["release.json"]
    if parent["document_references"]["release.json"]["sha256"] != "sha256:" + AEROSOL_RELEASE_SHA256:
        raise ValueError("aerosol publication parent differs")
    release_bytes = canonical_json_bytes(release)
    if hashlib.sha256(release_bytes).hexdigest() != AEROSOL_RELEASE_SHA256:
        raise ValueError("canonical aerosol release bytes differ")
    if (
        release["period"] != {"start": "2018-01-01T00:00:00Z", "end_exclusive": "2026-01-01T00:00:00Z", "cadence": "PT3H"}
        or release["source_grid"]["shape"] != [241, 480]
        or release["wavelengths_nm"] != [500, 550, 1064, 1550, 1640]
        or len(release["tiles"]) != 120
    ):
        raise ValueError("aerosol source coverage differs")
    directory = output_root / "v21" / "sources" / release["release_id"]
    summary = directory / "map_summary.parquet"
    provenance_path = directory / "presentation.json"
    manifest_sha256 = parent["document_references"]["manifest.json"]["sha256"]
    if provenance_path.exists():
        provenance = read_json(provenance_path)
        if (
            provenance.get("parent_release_sha256") != "sha256:" + AEROSOL_RELEASE_SHA256
            or not summary.is_file()
            or sha256(summary) != provenance.get("map_summary_sha256", "").removeprefix("sha256:")
        ):
            raise ValueError("existing aerosol presentation cache differs; not overwriting it")
        # Keep the downloadable parent byte-identical to its pinned R2 hash.
        (directory / "release.json").write_bytes(release_bytes)
        return release, summary, manifest_sha256, provenance
    columns = ["source_pixel_id", "source_grid_row", "source_grid_column", "wavelength_nm"] + [
        f"overall_{field}" for field, _ in AEROSOL_STATISTICS.values()
    ] + ["overall_expected_count"]

    def read_tile(tile: dict[str, Any]) -> pa.Table:
        row = tile["objects"]["statistics.parquet"]
        prefix = tile["manifest_key"].removesuffix("manifest.json")
        if row["key"] != prefix + "statistics.parquet" or not 0 < row["size_bytes"] <= 32 * 1024 * 1024:
            raise ValueError("aerosol statistics reference exceeds the bounded projection")
        body = client.get_object(Bucket=target.bucket, Key=row["key"])["Body"]
        try:
            data = body.read(row["size_bytes"] + 1)
        finally:
            body.close()
        if len(data) != row["size_bytes"] or hashlib.sha256(data).hexdigest() != row["sha256"].removeprefix("sha256:"):
            raise ValueError(f"aerosol tile read-back differs: {tile['shard_id']}")
        return pq.read_table(io.BytesIO(data), columns=columns)

    tables: list[pa.Table] = []
    with ThreadPoolExecutor(max_workers=4) as pool:
        for start in range(0, len(release["tiles"]), 4):
            # Only four compressed objects and their projected columns in flight.
            tables.extend(pool.map(read_tile, release["tiles"][start:start + 4]))
            if (start + 4) % 20 == 0:
                print(json.dumps({"aerosol_statistics_tiles_verified": start + 4, "total": 120}), flush=True)
    table = pa.concat_tables(tables).sort_by([("wavelength_nm", "ascending"), ("source_pixel_id", "ascending")])
    validate_aerosol_summary(table)
    directory.mkdir(parents=True, exist_ok=True)
    pq.write_table(table, summary, compression="zstd")
    (directory / "release.json").write_bytes(release_bytes)
    provenance = {
        "kind": "skyclimate-existing-aerosol-map-projection-v1",
        "parent_release_id": release["release_id"],
        "parent_release_sha256": "sha256:" + AEROSOL_RELEASE_SHA256,
        "parent_receipt": parent["receipt_reference"],
        "parent_documents": parent["document_references"],
        "map_summary_sha256": "sha256:" + sha256(summary),
        "map_summary_size_bytes": summary.stat().st_size,
        "columns": columns,
        "method": "exact-column-selection-no-new-atmospheric-calculation",
        "all_120_parent_statistics_objects_size_and_sha256_verified": True,
    }
    provenance_path.write_text(json.dumps(provenance, indent=2) + "\n", encoding="utf-8")
    return release, summary, manifest_sha256, provenance


def validate_aerosol_summary(table: pa.Table) -> None:
    rows = table["source_grid_row"].to_numpy()
    columns = table["source_grid_column"].to_numpy()
    pixels = table["source_pixel_id"].to_numpy()
    waves = table["wavelength_nm"].to_numpy()
    if (
        table.num_rows != 241 * 480 * 5
        or np.any(rows < 0) or np.any(rows >= 241)
        or np.any(columns < 0) or np.any(columns >= 480)
        or not np.array_equal(pixels, rows.astype(np.int64) * 480 + columns)
        or set(np.unique(waves)) != {500, 550, 1064, 1550, 1640}
    ):
        raise ValueError("aerosol global grid inventory differs")
    for wavelength in (500, 550, 1064, 1550, 1640):
        selected = waves == wavelength
        if not np.array_equal(np.sort(pixels[selected]), np.arange(241 * 480)):
            raise ValueError("aerosol map has duplicate or missing source cells")
    expected = table["overall_expected_count"].to_numpy()
    count = table["overall_valid_count"].to_numpy()
    fraction = table["overall_valid_fraction"].to_numpy()
    if np.any(expected != 23376) or np.any(count < 0) or np.any(count > expected) or not np.allclose(fraction, count / expected, atol=1e-7):
        raise ValueError("aerosol map support differs from the native PT3H axis")


def aerosol_map_coordinates() -> tuple[np.ndarray, np.ndarray]:
    # CAMS temporal admission normalizes longitude to [-180, 180), unlike
    # the cloud presentation's [0, 360) parent. Never reuse its half-world roll.
    return 90 - np.arange(241) * 0.75, -180 + np.arange(480) * 0.75


def append_aerosol(args: argparse.Namespace) -> None:
    base = read_json(args.base_catalogue)
    original_catalogue_sha256 = sha256(args.base_catalogue)
    cloud = next(item for item in base["products"] if item["productId"] == EXPECTED_PACK_ID)
    release, summary, manifest_sha256, provenance = fetch_aerosol_summary(args.output_root)
    table = pq.read_table(summary)
    validate_aerosol_summary(table)
    product_id = release["release_id"]
    source_latitudes, source_longitudes = aerosol_map_coordinates()
    latitude_keep = np.flatnonzero(np.abs(source_latitudes) <= WEB_MERCATOR_MAX_CELL_CENTRE_LAT)
    longitude_order = np.argsort(source_longitudes)
    latitudes, longitudes = source_latitudes[latitude_keep], source_longitudes[longitude_order]
    variants = []
    waves = table["wavelength_nm"].to_numpy()
    rows, columns = table["source_grid_row"].to_numpy(), table["source_grid_column"].to_numpy()
    for wavelength in release["wavelengths_nm"]:
        selected = waves == wavelength
        for statistic, (field, label) in AEROSOL_STATISTICS.items():
            full = np.full((241, 480), np.nan, dtype=np.float32)
            full[rows[selected], columns[selected]] = table[f"overall_{field}"].to_numpy(zero_copy_only=False)[selected]
            values = full[np.ix_(latitude_keep, longitude_order)]
            display, coordinates = project(values, latitudes, longitudes)
            high = 23376.0 if statistic == "count" else 1.0 if statistic == "coverage" else 0.6 if wavelength <= 550 else 0.3
            palette = "coverage" if statistic in {"count", "coverage"} else "aerosol"
            slug = f"scar21-aefc5a-aod-{wavelength}nm-{statistic}-overall"
            variants.append({
                "id": slug, "productId": product_id, "variableId": f"aod-{wavelength}nm",
                "variableLabel": f"Aerosol optical depth · {wavelength} nm",
                "statisticId": statistic, "statisticLabel": label,
                "temporalKind": "overall", "temporalValue": None,
                "unit": "samples" if statistic == "count" else "fraction" if statistic == "coverage" else "1",
                "verticalSupport": "total atmospheric column",
                "grid": {
                    "width": 480, "height": len(latitudes), "north": float(latitudes[0]),
                    "south": float(latitudes[-1]), "west": float(longitudes[0]), "east": float(longitudes[-1]),
                    "latitudeStep": 0.75, "longitudeStep": 0.75, "longitudeCyclic": True,
                    "imageProjection": "web-mercator-nearest", "imageWidth": int(display.shape[1]),
                    "imageHeight": int(display.shape[0]), "imageCoordinates": coordinates,
                },
                **write_assets(args.output_root, slug, values, display, 0.0, high, palette),
            })
    base["products"] = [item for item in base["products"] if item["productId"] != product_id] + [{
        **cloud, "productId": product_id, "label": "SkyClimate five-band aerosol atlas",
        "source": "CAMS EAC4 native and Ångström-derived aerosol optical depth",
        "claimClass": "model-derived-screening", "periodLabel": "2018–2025 · native three-hourly reference",
        "cadence": "PT3H · 23,376 three-hourly instants", "resolutionDeg": 0.75,
        "sourceUrl": "https://ads.atmosphere.copernicus.eu/datasets/cams-global-reanalysis-eac4",
        "datasetDoi": "10.24381/d58bbf47",
        "attribution": "Generated using Copernicus Atmosphere Monitoring Service information [2026]",
        "modificationIndication": "Contains modified Copernicus Atmosphere Monitoring Service information [2026]. AOD at 500 and 1064 nm is spectrally interpolated; 1550 nm is extrapolated from 865/1240 nm; 550 and 1640 nm are native channels. Statistics retain native 0.75° support and PT3H cadence.",
        "manifestSha256": manifest_sha256, "statisticsSha256": provenance["map_summary_sha256"],
        "validationSourceIds": ["aeronet-v3-aod"],
        "scientificValidation": {
            f"aod-{wave}nm": {
                # The existing browser's general 'supported' rule requires
                # BOTH primary and mountain evidence. Do not upgrade it from
                # this family's narrower global result.
                "state": "unsupported",
                "reason": f"Global frozen project gate {item['frozen_project_gate']}; {item['use']}. General/mountain transfer is not validated and no local calibration is claimed. The full family report is linked in provenance.",
            } for wave, item in release["validation"]["bands"].items()
        },
        "uncertainty": {"state": "unavailable", "method": None, "nominalCoverage": None,
                        "reason": "No calibrated local uncertainty; spectral validation is channel-specific and does not establish local station accuracy."},
        "presentationDataUrl": f"/data/v21/sources/{product_id}/map_summary.parquet",
        "presentationManifestUrl": f"/data/v21/sources/{product_id}/presentation.json",
    }]
    base["variants"] = [item for item in base["variants"] if item["productId"] != product_id] + variants
    base["validationSources"] = [item for item in base["validationSources"] if item["serviceId"] != "aeronet-validation"]
    if not any(item["serviceId"] == "aeronet-validation" for item in base["validationSources"]):
        base["validationSources"].append({
            "serviceId": "aeronet-validation", "sourceIds": ["aeronet-v3-aod"],
            "deliveryState": "aggregated-validation-evidence-only", "rawObservationRowsIncluded": False,
            "licenceUrl": "https://aeronet.gsfc.nasa.gov/new_web/data_usage.html",
            "sourceTermsUrl": "https://aeronet.gsfc.nasa.gov/new_web/data_usage.html",
            "notices": {
                "source-citation": "Giles, D. M., et al. (2019), Advancements in the Aerosol Robotic Network (AERONET) Version 3 database – automated near-real-time quality control algorithm with improved cloud screening for Sun photometer aerosol optical depth measurements, Atmospheric Measurement Techniques, 12, 169–209, https://doi.org/10.5194/amt-12-169-2019.",
                "network-acknowledgement": "We thank the AERONET Principal Investigators and their staff for establishing and maintaining the sites used in the aggregated validation evidence. No observation rows are redistributed.",
            },
        })
    audit = base["deliveryPolicy"]["catalogueAudit"]
    audit["product_count"] = len(base["products"])
    audit["dataset_dois"] = sorted({item["datasetDoi"] for item in base["products"]})
    audit["validation_source_count"] = len(base["validationSources"])
    audit["validation_source_ids"] = sorted({source for item in base["validationSources"] for source in item["sourceIds"]})
    limitation = "Aerosol maps expose 55 overall summaries from the complete 2018–2025 native PT3H CAMS axis, not hourly ERA5 data. The five wavelengths have different validation results; none is locally calibrated. Monthly and UTC-hour statistics and temporal sketches remain in the parent family release."
    if limitation not in base["limitations"]:
        base["limitations"].append(limitation)
    if sha256(args.base_catalogue) != original_catalogue_sha256:
        raise ValueError("the base catalogue changed during presentation building; not overwriting it")
    output = args.output_root / "catalog-v21.json"
    output.write_text(json.dumps(base, indent=2) + "\n", encoding="utf-8")
    print(json.dumps({"catalogue": str(output), "aerosol_product_id": product_id, "aerosol_variants": len(variants), "source_projection": str(summary)}))


def validate_v4_demo_lineage(manifest: dict, report: dict, freeze: dict) -> None:
    """A reference demonstration is not an engineering or calibration admission."""
    january = next((item for item in freeze.get("months", []) if item.get("partition_month") == "2018-01"), {})
    if (
        report.get("status") != "succeeded" or report.get("provider_readback_verified") is not True
        or report.get("summary_id") != V4_SUMMARY_ID or manifest.get("summary_id") != V4_SUMMARY_ID
        or manifest.get("kind") != "skyclimate-v21-l137-reference-monthly-summary-v1"
        or manifest.get("claim_class") != "model-derived-screening"
        or manifest.get("release_authority") != "phase-f-evaluation-only"
        or manifest.get("input_pack_id") != january.get("monthly_pack_id")
        or report.get("input_pack_id") != january.get("monthly_pack_id")
        or manifest.get("input_derivation_identity") != january.get("processing_identity")
        or manifest.get("input_derivation_identity", {}).get("turbulence_method_id") != "era5-hmnsp99-monin-obukhov-v4"
        or january.get("durable_publication", {}).get("provider_readback_verified") is not True
        or manifest.get("metric_inventory") != list(V4_METRICS)
        or manifest.get("metric_units") != {name: item[2] for name, item in V4_METRICS.items()}
        or manifest.get("month") != "2018-01" or manifest.get("source_cadence") != "PT3H"
        or manifest.get("valid_time_count") != 248
        or manifest.get("valid_time_start_utc") != "2018-01-01T00:00:00Z"
        or manifest.get("valid_time_stop_exclusive_utc") != "2018-02-01T00:00:00Z"
        or manifest.get("values", {}).get("shape") != [6, 721, 1440]
        or manifest.get("statistic") != "p50-linear-over-independent-finite-values"
        or manifest.get("missingness") != "no-imputation"
        or manifest.get("minimum_finite_count") != 1
        or set(manifest.get("artifacts", {})) != {"p50.npy", "finite_count.npy", "latitude.npy", "longitude.npy"}
    ):
        raise ValueError("V4 demo lineage, method, coverage or support differs from the retained reference")


def fetch_v4_demo_summary(output_root: Path) -> tuple[dict, Path, dict]:
    from skyclimate.runner import R2Target

    client, bucket = R2Target().client(), "optical-ground-siting-derived"
    directory = output_root / "v21" / "sources" / V4_SUMMARY_ID
    directory.mkdir(parents=True, exist_ok=True)

    def verified(key: str, digest: str, size: int | None, name: str) -> bytes:
        path = directory / name
        if path.exists():
            data = path.read_bytes()
        else:
            response = client.get_object(Bucket=bucket, Key=key)
            try:
                if response["ContentLength"] > 50_000_000:
                    raise ValueError("V4 demo read exceeds a precomputed summary object")
                data = response["Body"].read(50_000_001)
            finally:
                response["Body"].close()
        if (size is not None and len(data) != size) or hashlib.sha256(data).hexdigest() != digest.removeprefix("sha256:"):
            raise ValueError(f"V4 demo object size/SHA mismatch: {name}")
        if not path.exists():
            path.write_bytes(data)
        return data

    report = json.loads(verified(V4_REPORT_KEY, V4_REPORT_SHA256, None, "source-report.json"))
    freeze = json.loads(verified(V4_FREEZE_KEY, V4_FREEZE_SHA256, None, "reference-freeze.json"))
    prefix = f"skyclimate/turbulence/pressure-level-phase-f/v1/l137-reference-summaries/{V4_SUMMARY_ID}"
    records = {item["key"]: item for item in report["r2_summary_objects"]}
    record = records[f"{prefix}/manifest.json"]
    manifest = json.loads(verified(record["key"], record["sha256"], record["size_bytes"], "source-manifest.json"))
    validate_v4_demo_lineage(manifest, report, freeze)
    arrays = {}
    for name, artifact in manifest["artifacts"].items():
        record = records[f"{prefix}/{name}"]
        if (record.get("provider_readback_verified") is not True
                or record["size_bytes"] != artifact["size_bytes"]
                or record["sha256"] != artifact["sha256"].removeprefix("sha256:")):
            raise ValueError("V4 summary receipt does not bind every artifact")
        arrays[name] = np.load(io.BytesIO(verified(record["key"], record["sha256"], record["size_bytes"], name)), allow_pickle=False)
    latitude, longitude = arrays["latitude.npy"], arrays["longitude.npy"]
    medians, counts = arrays["p50.npy"], arrays["finite_count.npy"]
    if (not np.array_equal(latitude, 90 - np.arange(721) * 0.25)
            or not np.array_equal(longitude, np.arange(1440) * 0.25)
            or medians.shape != (6, 721, 1440) or counts.shape != medians.shape
            or np.any(counts > 248) or np.any(counts < 0)
            or not np.array_equal(np.isfinite(medians), counts > 0)):
        raise ValueError("V4 summary coordinates or finite support differ")
    summary = directory / "map_summary.parquet"
    columns = {"latitude_deg": np.repeat(latitude, 1440), "longitude_deg": np.tile(longitude, 721)}
    for index, name in enumerate(V4_METRICS):
        columns[name + "_p50"] = medians[index].reshape(-1)
        columns[name + "_finite_count"] = counts[index].reshape(-1)
    pq.write_table(pa.table(columns), summary, compression="zstd")
    provenance = {
        "kind": "skyclimate-existing-v4-reference-presentation-v1",
        "purpose": "demo-presentation-only", "engineering_release_admission": False,
        "claim_class": "model-derived-screening", "local_calibration": False,
        "source_manifest_sha256": "sha256:" + sha256(directory / "source-manifest.json"),
        "source_report_sha256": "sha256:" + V4_REPORT_SHA256,
        "reference_freeze_sha256": "sha256:" + V4_FREEZE_SHA256,
        "summary_id": V4_SUMMARY_ID, "parent_monthly_pack_id": manifest["input_pack_id"],
        "map_period": {"start": manifest["valid_time_start_utc"], "end_exclusive": manifest["valid_time_stop_exclusive_utc"], "cadence": "PT3H", "count": 248},
        "map_summary_sha256": "sha256:" + sha256(summary),
        "map_summary_size_bytes": summary.stat().st_size,
        "method": "exact-existing-monthly-median-and-support-column-selection; no-new-turbulence-calculation",
        "source_release_authority_unchanged": manifest["release_authority"],
        "all_summary_objects_size_and_sha256_verified": True,
        "limitations": ["January 2018 only; not a multi-year or seasonal climatology.",
                        "Median uses each metric's finite samples independently; count is provided and can be as low as one.",
                        "Model-derived V4 screening; surface-layer turbulence can be underestimated and r0 optimistic.",
                        "Reference Rytov variance is not a measured scintillation or a directional link prediction.",
                        "PWV is above the ERA5 model surface, not corrected to a candidate's actual altitude."],
    }
    (directory / "presentation.json").write_text(json.dumps(provenance, indent=2) + "\n", encoding="utf-8")
    return manifest, summary, provenance


def append_v4_demo(args: argparse.Namespace) -> None:
    base, original_sha = read_json(args.base_catalogue), sha256(args.base_catalogue)
    cloud = next(item for item in base["products"] if item["productId"] == EXPECTED_PACK_ID)
    _manifest, summary, provenance = fetch_v4_demo_summary(args.output_root)
    table = pq.read_table(summary)
    latitudes = 90 - np.arange(721) * 0.25
    keep = np.flatnonzero(np.abs(latitudes) <= WEB_MERCATOR_MAX_CELL_CENTRE_LAT)
    latitudes = latitudes[keep]
    source_longitudes = ((np.arange(1440) * 0.25 + 180) % 360) - 180
    order = np.argsort(source_longitudes)
    longitudes = source_longitudes[order]
    variants = []
    for name, (variable, label, unit, upper) in V4_METRICS.items():
        for statistic, suffix, high, palette in [("p50", "p50", upper, "water"), ("count", "finite_count", 248., "coverage")]:
            values = table[name + "_" + suffix].to_numpy().reshape(721, 1440)[np.ix_(keep, order)]
            display, coordinates = project(values, latitudes, longitudes)
            slug = f"scv4demo-2018-01-{variable}-{statistic}"
            variants.append({
                "id": slug, "productId": V4_SUMMARY_ID, "variableId": variable, "variableLabel": label,
                "statisticId": statistic, "statisticLabel": "Median" if statistic == "p50" else "Valid three-hourly samples",
                "temporalKind": "overall", "temporalValue": None, "unit": unit if statistic == "p50" else "samples",
                "verticalSupport": "above ERA5 model surface",
                "grid": {"width": 1440, "height": len(latitudes), "north": float(latitudes[0]), "south": float(latitudes[-1]),
                         "west": -180.0, "east": 179.75, "latitudeStep": 0.25, "longitudeStep": 0.25,
                         "longitudeCyclic": True, "imageProjection": "web-mercator-nearest", "imageWidth": int(display.shape[1]),
                         "imageHeight": int(display.shape[0]), "imageCoordinates": coordinates},
                **write_assets(args.output_root, slug, values, display, 0., high, palette),
            })
    product = {
        **cloud, "productId": V4_SUMMARY_ID, "label": "V4 turbulence demo · January 2018 only",
        "source": "ERA5 L137 · retained V4 model-derived reference summary", "claimClass": "model-derived-screening",
        "periodStart": "2018-01-01T00:00:00Z", "periodEndExclusive": "2018-02-01T00:00:00Z",
        "periodLabel": "January 2018 only · partial V4 screening demo", "cadence": "PT3H · 248 instants; finite support varies by field and cell",
        "sourceUrl": "https://cds.climate.copernicus.eu/datasets/reanalysis-era5-complete", "datasetDoi": "10.24381/cds.143582cf",
        "modificationIndication": "Contains modified Copernicus Climate Change Service information [2026]. Retained V4 monthly medians and per-field finite counts only. Not an engineering admission or a complete climatology.",
        "manifestSha256": provenance["source_manifest_sha256"], "statisticsSha256": provenance["map_summary_sha256"],
        "scientificValidation": {item[0]: {"state": "unsupported", "reason": "V4 screening demo; not locally calibrated. " + " ".join(provenance["limitations"][2:])} for item in V4_METRICS.values()},
        "uncertainty": {"state": "unavailable", "method": None, "nominalCoverage": None, "reason": "No calibrated local uncertainty; missing cells and finite sample counts are preserved."},
        "presentationDataUrl": f"/data/v21/sources/{V4_SUMMARY_ID}/map_summary.parquet",
        "presentationManifestUrl": f"/data/v21/sources/{V4_SUMMARY_ID}/presentation.json",
    }
    base["products"] = [item for item in base["products"] if item["productId"] != V4_SUMMARY_ID] + [product]
    base["variants"] = [item for item in base["variants"] if item["productId"] != V4_SUMMARY_ID] + variants
    audit = base["deliveryPolicy"]["catalogueAudit"]
    audit["product_count"], audit["dataset_dois"] = len(base["products"]), sorted({item["datasetDoi"] for item in base["products"]})
    limitation = "V4 demo maps show January 2018 medians and finite counts only, not full-period turbulence. The retained reference corpus and its original evaluation-only authority are preserved; display does not constitute engineering admission or local calibration."
    if limitation not in base["limitations"]:
        base["limitations"].append(limitation)
    if sha256(args.base_catalogue) != original_sha:
        raise ValueError("the base catalogue changed; not overwriting it")
    (args.output_root / "catalog-v21.json").write_text(json.dumps(base, indent=2) + "\n", encoding="utf-8")
    print(json.dumps({"v4_demo_product": V4_SUMMARY_ID, "variants": len(variants), "new_atmospheric_calculation": False}))


def write_native_cog(path: Path, values: np.ndarray, step: float, tags: dict) -> None:
    """Lossless file-format view: preserve centres, values and missingness."""
    import rasterio
    from rasterio.transform import from_origin

    if values.shape != (round(180 / step) + 1, round(360 / step)):
        raise ValueError("COG must preserve the complete native source grid")
    path.parent.mkdir(parents=True, exist_ok=True)
    profile = {"driver": "COG", "width": values.shape[1], "height": values.shape[0], "count": 1,
               "dtype": values.dtype, "crs": "EPSG:4326", "nodata": np.nan,
               "transform": from_origin(-180 - step / 2, 90 + step / 2, step, step),
               "compress": "DEFLATE", "blocksize": 256, "overview_resampling": "NEAREST"}
    if not path.exists():
        with rasterio.open(path, "w", **profile) as dataset:
            dataset.write(values, 1)
            dataset.update_tags(**tags, AREA_OR_POINT="Point", overview_method="nearest-no-new-information")
    with rasterio.open(path) as dataset:
        if (dataset.tags(ns="IMAGE_STRUCTURE").get("LAYOUT") != "COG"
                or dataset.transform != profile["transform"]
                or not np.array_equal(dataset.read(1), values, equal_nan=True)
                or any(dataset.tags().get(key) != str(value) for key, value in tags.items())):
            raise ValueError("COG read-back differs from its exact source-grid view")


def package_existing(args: argparse.Namespace) -> None:
    """Build distribution files only; source kernels and identities stay untouched."""
    destination = args.package_existing.resolve()
    if not str(destination).startswith("/srv/data/optical-ground-siting/publication-staging/skyclimate/"):
        raise ValueError("publication staging must stay in the bounded data root")
    if (destination / "manifest.json").exists() and not args.refresh_unpublished_staging:
        raise ValueError("package exists; an explicit unpublished-staging refresh or a new release identity is required")
    if (destination.parent / (destination.name + "-publication-receipt.json")).exists():
        raise ValueError("published release cannot be refreshed")
    inventory = read_json(destination / "native-data.json")
    groups = inventory["groups"]
    current_check = destination.parent / (destination.name + "-native-readback-check.json")
    check = read_json(current_check)
    if (check.get("all_sizes_and_stored_sha256_match") is not True
            or check.get("index_inventory_sha256") != sha256(destination / "native-data.json")):
        raise ValueError("current source-object metadata verification is required")
    shutil.copyfile(current_check, destination / "native-source-check.json")
    if (len([g for g in groups if g["family"] == "cloud-tcwv-hourly"]) != 96
            or len([g for g in groups if g["group"].startswith("aerosol-month-")]) != 96):
        raise ValueError("publication requires the existing complete cloud and aerosol temporal axes")
    catalogue = read_json(args.base_catalogue)
    if len(catalogue["products"]) != 3 or len(catalogue["variants"]) != 112:
        raise ValueError("existing-data family presentation differs")
    source_root = args.output_root / "v21/sources"
    cloud_path = Path("/srv/data/optical-ground-siting/skyclimate-v21-presentation-source/map_summary.parquet")
    if sha256(cloud_path) != EXPECTED_MAP_SHA256:
        raise ValueError("cloud presentation source differs")
    aerosol_id = next(p["productId"] for p in catalogue["products"] if p["productId"].startswith("scar21-"))
    paths = {EXPECTED_PACK_ID: cloud_path, aerosol_id: source_root / aerosol_id / "map_summary.parquet",
             V4_SUMMARY_ID: source_root / V4_SUMMARY_ID / "map_summary.parquet"}
    products = {p["productId"]: p for p in catalogue["products"]}
    for product, path in paths.items():
        if product != EXPECTED_PACK_ID:
            proof = read_json(path.parent / "presentation.json")
            if sha256(path) != proof["map_summary_sha256"].removeprefix("sha256:"):
                raise ValueError("presentation projection differs from its verified parent")
            if product == V4_SUMMARY_ID and (proof.get("engineering_release_admission") is not False
                                           or proof.get("source_release_authority_unchanged") != "phase-f-evaluation-only"):
                raise ValueError("V4 demo must not acquire an engineering admission")
    tables = {key: pq.read_table(path) for key, path in paths.items()}
    release_id = destination.name
    base_url = f"https://data.dynamikorbits.com/skyclimate/{release_id}"
    layers, stac_links, cog_checks = [], [], []
    for variant in catalogue["variants"]:
        product = variant["productId"]
        variable, statistic = variant["variableId"], variant["statisticId"]
        table = tables[product]
        step = 0.75 if product == aerosol_id else 0.25
        values = np.full((241, 480) if step == 0.75 else (721, 1440), np.nan, dtype=np.float64)
        if product == EXPECTED_PACK_ID:
            column = variable + "_" + {"coverage": "valid_fraction", "count": "valid_count"}.get(statistic, statistic)
            rows, columns = table["source_latitude_index"].to_numpy(), table["source_longitude_index"].to_numpy()
            values[rows, columns] = table[column].to_numpy(zero_copy_only=False)
            values = np.roll(values, 720, axis=1)
        elif product == aerosol_id:
            wavelength = int(variable.removeprefix("aod-").removesuffix("nm"))
            selected = table["wavelength_nm"].to_numpy() == wavelength
            rows, columns = table["source_grid_row"].to_numpy()[selected], table["source_grid_column"].to_numpy()[selected]
            column = "overall_" + AEROSOL_STATISTICS[statistic][0]
            values[rows, columns] = table[column].to_numpy(zero_copy_only=False)[selected]
        else:
            metric = next(name for name, item in V4_METRICS.items() if item[0] == variable)
            column = metric + ("_finite_count" if statistic == "count" else "_p50")
            values = np.roll(table[column].to_numpy(zero_copy_only=False).reshape(721, 1440), 720, axis=1)
        # Float64 preserves every source value, including the retained V4 medians.
        values = values.astype(np.float64)
        relative = "cog/" + variant["id"] + ".tif"
        source = products[product]
        is_demo = product == V4_SUMMARY_ID
        tags = {"parent_product_id": product, "parent_summary_sha256": sha256(paths[product]),
                "claim_class": source["claimClass"], "units": variant["unit"], "variable": variable,
                "statistic": statistic, "engineering_release_admission": "false" if is_demo else "unchanged-parent-authority"}
        write_native_cog(destination / relative, values, step, tags)
        cog_checks.append({"path": relative, "parent_id": product, "source_column": column,
                           "all_native_values_and_missingness_equal": True, "resolution_degrees": step,
                           "dtype": "float64", "resampling": "none; nearest overviews only"})
        period = {"start": "2018-01-01T00:00:00Z", "end_exclusive": "2018-02-01T00:00:00Z" if is_demo else "2026-01-01T00:00:00Z"}
        name = f"{'DEMO Jan 2018 · ' if is_demo else '2018–2025 · '}{variant['variableLabel']} · {variant['statisticLabel']}"
        item = {"type": "Feature", "stac_version": "1.1.0", "id": variant["id"],
                "geometry": {"type": "Polygon", "coordinates": [[[-180,-90],[180,-90],[180,90],[-180,90],[-180,-90]]]},
                "bbox": [-180,-90,180,90], "properties": {"datetime": None, "start_datetime": period["start"],
                    "end_datetime": "2018-01-31T21:00:00Z" if is_demo else "2025-12-31T21:00:00Z" if product == aerosol_id else "2025-12-31T23:00:00Z",
                    "title": name, "license": "CC-BY-NC-4.0", "skyclimate:parent_id": product,
                    "skyclimate:claim_class": source["claimClass"], "skyclimate:demo_only": is_demo,
                    "skyclimate:unit": variant["unit"], "skyclimate:resolution_degrees": step},
                "links": [{"rel": "root", "href": "../catalog.json", "type": "application/json"}],
                "assets": {"data": {"href": "../" + relative, "type": "image/tiff; application=geotiff; profile=cloud-optimized", "roles": ["data", "visual"]}}}
        item_path = destination / "stac" / (variant["id"] + ".json")
        item_path.parent.mkdir(parents=True, exist_ok=True)
        item_path.write_text(json.dumps(item, indent=2) + "\n")
        stac_links.append({"rel": "item", "href": "./stac/" + item_path.name, "type": "application/json", "title": name})
        if statistic == ("p50" if is_demo else "mean"):
            layers.append({"id": variant["id"], "name": name, "type": "cog", "source": {"type": "raster", "url": base_url + "/" + relative},
                           "sourcePath": base_url + "/" + relative, "visible": variable == "tcc", "opacity": 0.85,
                           "metadata": {"sourceKind": "maplibre-gl-raster", "externalNativeLayer": True,
                               "externalDeckLayer": True, "customLayerType": "raster", "rasterSource": "url",
                               "nativeLayerIds": [variant["id"]], "sourceIds": [], "panelCollapsed": True,
                               "rasterOverlayMode": "interleaved", "bandCount": 1,
                               "rasterState": {"mode": "single", "bands": [1], "nodata": "auto", "colormap": "viridis", "rescale": [[variant["visualMin"], variant["visualMax"]]]},
                               "parent_id": product, "claim_class": source["claimClass"], "demo_only": is_demo}, "style": {}})
    # Source projections and provenance travel with the public package. The
    # original >500 GB corpus remains behind the exact-object alias indexes.
    for product, path in paths.items():
        target = destination / ("reference-demo" if product == V4_SUMMARY_ID else "summaries") / product
        target.mkdir(parents=True, exist_ok=True)
        for source in ([path] if product == EXPECTED_PACK_ID else sorted(path.parent.glob("*.json")) + [path]):
            output = target / source.name
            if not output.exists():
                shutil.copyfile(source, output)
            if sha256(output) != sha256(source):
                raise ValueError("publication copy differs from its verified source")
    def document(name, value):
        (destination / name).write_text(json.dumps(value, indent=2) + "\n", encoding="utf-8")
    reproduction = destination / "reproduction"
    reproduction.mkdir(exist_ok=True)
    for source in (Path(__file__), Path(__file__).with_name("fetch-v21-atlas-presentation-source.py"), Path("uv.lock")):
        shutil.copyfile(source, reproduction / source.name)
    document("catalog.json", {"stac_version": "1.1.0", "type": "Catalog", "id": release_id,
             "description": "Existing SkyClimate data: full-period cloud/TCWV/aerosol data plus a separate January 2018 V4 reference demo.",
             "links": [{"rel": "self", "href": base_url + "/catalog.json"}, *stac_links]})
    document("skyclimate.geolibre.json", {"version": "0.1.0", "name": "SkyClimate · existing data + V4 demo",
             "mapView": {"center": [0,20], "zoom": 1.4, "bearing": 0, "pitch": 0},
             "basemapStyleUrl": "https://tiles.openfreemap.org/styles/positron", "basemapVisible": True,
             "layers": layers, "styles": {}, "metadata": {"release_url": base_url, "licence": "CC-BY-NC-4.0",
                 "all_statistics": base_url + "/catalog.json", "native_temporal_and_sketch_data": base_url + "/native-data.json"}})
    document("provenance.json", {"kind": "skyclimate-existing-data-publication-v1", "release_id": release_id,
             "products": catalogue["products"], "cog_validation": cog_checks,
             "scope": "Existing data only; no new production, no new calibration, no altered source identity or gate",
             "v4": {"public_scope": "January 2018 reference-summary demo", "engineering_release_admission": False,
                    "original_authority": "phase-f-evaluation-only", "other_retained_periods": "Not implied to be part of this public summary demo"},
             "durability": "All aliased exact derived keys are dependencies of this immutable public release and must not be evicted from R2.",
             "git_base_revision": subprocess.check_output(["git", "rev-parse", "HEAD"], text=True).strip(),
             "builder_sources": {str(Path(__file__).relative_to(Path.cwd())): sha256(Path(__file__)),
                "apps/skyclimate-explorer/scripts/fetch-v21-atlas-presentation-source.py": sha256(Path(__file__).with_name("fetch-v21-atlas-presentation-source.py"))}})
    (destination / "README.md").write_text(f"""# SkyClimate: existing-data release and V4 demo

This is the dataset we can publish now, not a claim that the full augmented Atlas is finished.
The parent datasets and their scientific claims are unchanged. No missing month was filled,
no new atmospheric calculation was run, and no calibration was deployed for this publication.

## What is included

Cloud cover (total, low, middle, high) and total column water vapour (TCWV): global ERA5,
2018–2025, hourly, 0.25° source grid. All 96 monthly histories, spatial summaries and the
existing full-catalogue cloud-only and TCWV temporal descriptors are accessible.
TCWV is the integrated atmospheric water-vapour column, not cloud cover and not a
site-altitude-adjusted PWV estimate. Low/middle/high cloud are model layers, not resolved 3D clouds.

Aerosols: global CAMS, 2018–2025, every three hours, 0.75° grid, at 500, 550, 1064,
1550 and 1640 nm. The complete 96 monthly histories, 120 statistical/descriptor tiles
and descriptor normalisation are included. 550 and 1640 nm are native channels;
other wavelengths use the retained spectral-transfer method. AOD is dimensionless
column extinction (scattering plus absorption), not a directional link loss in dB.
The frozen global 1550/1640 tests passed; 500/550/1064 did not, and the mountainous-site
test did not pass at any band. These are not generally validated or locally calibrated products.

V4 turbulence: **January 2018 reference-summary demo only**, 248 three-hourly samples,
global 0.25°. The six medians and their finite sample counts cover integrated Cn²,
r0, theta0, tau0, reference Rytov variance (optical fields at 1550 nm), and PWV above
the model surface. The original authority remains phase-F evaluation-only: this attachment
does not admit a full engineering turbulence release. Surface-layer turbulence is
under-represented; r0 may be optimistic. Reference Rytov is not an arbitrary-terminal
scintillation prediction. This demo is neither eight-year climatology nor a CFLOS product.

## Open and download

Open [the GeoLibre project](./skyclimate.geolibre.json) in GeoSpatial Studio or GeoLibre.
Its 16 initial layers cover all variables; [STAC](./catalog.json) lists all 112
overall statistical views. Each COG preserves the full source-grid values and missingness
with lossless compression; overview pixels are nearest-neighbour display aids.
Cell centres include the poles; the small half-cell edge beyond ±90° is only raster encoding.
No smoothing or resolution improvement is claimed. The Explorer's PNG/u16 map caches are
separate, quantised display files and are not the analytical data.

[native-data.json](./native-data.json) is the index of all hourly and three-hourly packs
and the three separate descriptor families. For each group, read its `index` JSON, then
append the object's relative name to `public_data_prefix` under this release URL:
`{base_url}/`. Each object has its original SHA-256 and byte count. Access is anonymous
HTTPS, GET/HEAD, byte ranges and CORS; there is no bucket-listing or credential endpoint.
The public gateway fails closed if an object's stored checksum metadata or size differs.
The indexes reference existing R2 objects; there is no redundant full-corpus download or copy.
These exact keys are retained public-release dependencies, not disposable derived scratch.
Scripted clients should send an explicit application User-Agent (for example
`SkyClimate-client/1.0`); the existing Cloudflare policy rejects Python urllib's
default bot-like identifier. No access token is required.

The current analytical packs are the existing Parquet/binary temporal layout and NetCDF
statistics, **not newly converted Zarr**. Parquet row-group reads make one spatial tile
available without downloading a whole month. See each parent manifest for dtype, UTC
axis, quantisation and missing-value rules before decoding its binary `payload` column.
For hourly packs use `series/<variable>/tiles-XXXX.parquet` and its explicit tile row;
for aerosols use `series-aod1550.parquet` (or the other band) and the required row group.
Always verify the row's payload SHA-256 before reshaping. The reader implementations are
`r2_hourly_temporal_source.read_r2_hourly_series_tile` and
`aerosol_spectral_dossier._read_aerosol_temporal_tile` in the linked source repository.
Do not treat missing values as zeros or replicate coarse cells into finer observations.
The historic `exact-sites.parquet` pilot is deliberately not exposed; all global grid
cells remain accessible and descriptors use the full retained catalogue, not 186 sites.

The cloud and TCWV descriptors each have 178 nonempty spatial shards, with `candidates`,
`cells`, `descriptors` and `bitsets` files plus a global normalisation. Aerosol descriptors
use 120 native-grid tiles and their own normalisation. Read candidate-to-source-cell
joins and support flags; H3 resolution 5 is a catalogue key, not atmospheric resolution.
These deterministic retrieval descriptors are not learned embeddings, site rankings,
or a universal vector joining unrelated physical families.

## Licence and attribution

The combined SkyClimate distribution is [CC BY-NC 4.0](https://creativecommons.org/licenses/by-nc/4.0/).
Commercial use requires a separate licence from Dynamik Orbits; contact the maintainer
through [the project repository](https://github.com/DynamikOrbits/optical-ground-siting).
Original Copernicus source rights remain distinct and unchanged. Contains modified
Copernicus Climate Change Service information [2026] and Copernicus Atmosphere Monitoring
Service information [2026]. Neither the European Commission nor ECMWF is responsible
for any use that may be made of the Copernicus information or data this product contains.
AERONET is validation context, not republished station observations; acknowledge its
principal investigators and follow its data-use policy for separately obtained station data.

## Deferred, without changing this release

Full 2018–2025 V4 turbulence, RF/rain/wind extensions, further calibration and new paid
compute are deferred. Retained restart plans, budget lessons, provenance and original
method identities remain in the repository. Publishability is not scientific validation.
SkyDome owns directional rays and CFLOS; neither is claimed here.
""", encoding="utf-8")
    artifacts = []
    media = {".json": "application/json", ".md": "text/markdown; charset=utf-8", ".parquet": "application/vnd.apache.parquet", ".tif": "image/tiff"}
    for path in sorted(destination.rglob("*")):
        if not path.is_file() or path == destination / "manifest.json":
            continue
        relative = path.relative_to(destination).as_posix()
        roles = ["visual", "data"] if path.suffix == ".tif" else ["project"] if path.name == "skyclimate.geolibre.json" else ["analysis", "data"] if path.suffix == ".parquet" else ["metadata"]
        artifacts.append({"path": relative, "size_bytes": path.stat().st_size, "sha256": sha256(path), "media_type": media.get(path.suffix, "application/octet-stream"), "roles": roles})
    document("manifest.json", {"$schema": "https://data.dynamikorbits.com/schemas/release-manifest-v1.json", "schema_version": "1.0.0",
             "product_id": "skyclimate", "release_id": release_id, "created_at": "2026-09-06T00:00:00Z", "published_at": None,
             "claim_class": "existing-data-distribution-with-separate-evaluation-demo", "canonical_base_url": base_url,
             "statistics": {"cloud_and_tcwv_months": 96, "aerosol_months": 96, "cog_views": len(cog_checks), "native_objects": sum(g["object_count"] for g in groups),
                            "aliased_existing_bytes": sum(g["existing_source_bytes"] for g in groups), "new_publication_bytes": sum(a["size_bytes"] for a in artifacts)},
             "distributions": {"stac": "catalog.json", "project": "skyclimate.geolibre.json", "native_data": "native-data.json", "readme": "README.md"}, "artifacts": artifacts})
    latest = {"$schema": "https://data.dynamikorbits.com/schemas/latest-v1.json", "schema_version": "1.0.0", "product_id": "skyclimate",
              "release_id": release_id, "manifest": f"skyclimate/{release_id}/manifest.json", "manifest_sha256": sha256(destination / "manifest.json")}
    (destination.parent / "latest.json").write_text(json.dumps(latest, indent=2) + "\n")
    print(json.dumps({"staged_release": str(destination), "artifacts": len(artifacts), "new_bytes": sum(a["size_bytes"] for a in artifacts), "cog_native_value_checks": len(cog_checks)}))


def main() -> None:
    args = arguments()
    if args.package_existing:
        package_existing(args)
        return
    if args.append_aerosol:
        append_aerosol(args)
        return
    if args.append_v4_demo:
        append_v4_demo(args)
        return
    _, manifest, catalogue_product = verify_inputs(args)
    base = read_json(args.base_catalogue)
    variants = build_variants(args.map_summary, args.output_root)
    product = build_product(
        manifest,
        catalogue_product,
        manifest_sha256=sha256(args.manifest),
    )
    # The V2.1 Explorer is a release view, not a museum of superseded pilots.
    # Legacy catalogues remain immutable evidence outside this user-facing view.
    base["products"] = [product]
    base["variants"] = variants
    base["exactAtmosphere"] = None
    base["validationSources"] = []
    audit = base["deliveryPolicy"]["catalogueAudit"]
    dataset_dois = sorted({item["datasetDoi"] for item in base["products"]})
    base["deliveryPolicy"] = {
        "policyId": V21_POLICY_ID,
        "decisionState": "approved-for-engineering-v2.1",
        "decisionDate": "2026-09-02",
        "rawSourceAssetsIncluded": False,
        "aeronetObservationRowsIncluded": False,
        "catalogueAudit": {
            **audit,
            "policy_id": V21_POLICY_ID,
            "decision_state": "approved-for-engineering-v2.1",
            "product_count": len(base["products"]),
            "dataset_dois": dataset_dois,
            "validation_source_count": 0,
            "validation_source_ids": [],
            "accepted": True,
        },
    }
    limitations = [
        "The V2.1 presentation exposes the 45 overall layers physically present in the admitted GeoParquet map summary. The Atlas also contains 325 overall, month, UTC-hour and month-by-UTC-hour views per source cell; those remain available through the bounded API until the shared dynamic map service is admitted.",
        "The V2.1 cloud and TCWV fields are source-grid statistics over ERA5 reanalysis at 0.25 degree. They are not local observations, directional line-of-sight estimates, CFLOS or SkyDome products.",
        "Map images are presentation caches; clicked values come from the companion quantised source-grid arrays.",
        "The Web-Mercator browser map retains native source cells only through ±84.75° latitude; polar source cells remain in the immutable product and require a polar view.",
        "The combined SkyClimate Atlas is licensed for non-commercial use under CC BY-NC 4.0. Commercial use requires a separate licence from Dynamik Orbits.",
    ]
    base["limitations"] = limitations
    output = args.output_root / "catalog-v21.json"
    output.parent.mkdir(parents=True, exist_ok=True)
    output.write_text(json.dumps(base, indent=2, sort_keys=False) + "\n", encoding="utf-8")
    print(json.dumps({
        "catalogue": str(output),
        "product_id": EXPECTED_PACK_ID,
        "variant_count": len(variants),
        "map_summary_sha256": EXPECTED_MAP_SHA256,
    }, indent=2))


if __name__ == "__main__":
    main()