make_visual.py

# /// script
# requires-python = ">=3.11"
# dependencies = ["pystac-client", "planetary-computer", "rasterio", "numpy", "scipy", "matplotlib", "requests", "pyproj"]
# ///
"""Trust under uncertainty: MODIS NDVI over the Gimje-Mangyeong rice plain, 2012 vs 2017.

Run: uv run make_visual.py      (keyless; downloads are cached in ./cache)
"""
import datetime as dt
import hashlib
import json
import sys
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import planetary_computer as pc
import pystac_client
import rasterio
import requests
from matplotlib import font_manager
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Rectangle
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.warp import reproject
from rasterio.windows import from_bounds
from scipy.signal import savgol_filter

ROOT = Path(__file__).parent
CACHE, OUT = ROOT / "cache", ROOT / "out"
CACHE.mkdir(exist_ok=True)
OUT.mkdir(exist_ok=True)

INK, ACCENT, PAPER, MUTE = "#22252A", "#C67863", "#F6F4EF", "#8A8D91"
NDVI_CMAP = LinearSegmentedColormap.from_list("calm", ["#E4EBD8", "#B7CFA3", "#7FAE82", "#4C8A66", "#2A5E48", "#17382C"])

BBOX = (126.70, 35.70, 127.10, 35.95)  # lon/lat, Gimje-Mangyeong plain, Jeonbuk
POWER_PT = (126.90, 35.82)
DAMAGE, NORMAL = 2012, 2017
MAP_DOY = 257  # MODIS 16-day composite starting day-of-year used for the maps
UTM = "EPSG:32652"
STAC = "https://planetarycomputer.microsoft.com/api/stac/v1"
TILE = "h28v05"
MIN_VALID = 0.25  # composite is plotted only if this share of cropland pixels is clear
TAG = {DAMAGE: "typhoon year", NORMAL: "reference year"}
EVENTS = {2012: [("Bolaven", dt.date(2012, 8, 28)), ("Sanba", dt.date(2012, 9, 17))]}

MANIFEST = {"generated_utc": None, "aoi_bbox_lonlat": BBOX, "sources": []}


def sha(b: bytes) -> str:
    return hashlib.sha256(b).hexdigest()


def clean(href: str) -> str:
    return href.split("?")[0]


def modis_window(ds):
    t = Transformer.from_crs(4326, ds.crs, always_xy=True)
    pts = [t.transform(lon, lat) for lon in (BBOX[0], BBOX[2]) for lat in (BBOX[1], BBOX[3])]
    xs, ys = [p[0] for p in pts], [p[1] for p in pts]
    return from_bounds(min(xs), min(ys), max(xs), max(ys), ds.transform).round_offsets().round_lengths()


def read_composite(item):
    cached = CACHE / f"{item.id}.npz"
    if not cached.exists():
        with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif"):
            with rasterio.open(pc.sign(item.assets["250m_16_days_NDVI"].href)) as ds:
                win = modis_window(ds)
                ndvi = ds.read(1, window=win)
                tf, crs = ds.window_transform(win), ds.crs.to_wkt()
            with rasterio.open(pc.sign(item.assets["250m_16_days_pixel_reliability"].href)) as ds:
                rel = ds.read(1, window=win)
        np.savez(cached, ndvi=ndvi, rel=rel, tf=np.array(tf)[:6], crs=crs)
    z = np.load(cached)
    return z["ndvi"], z["rel"], rasterio.Affine(*z["tf"]), str(z["crs"])


def utm_grid():
    t = Transformer.from_crs(4326, UTM, always_xy=True)
    pts = [t.transform(lon, lat) for lon in (BBOX[0], BBOX[2]) for lat in (BBOX[1], BBOX[3])]
    x0, x1 = min(p[0] for p in pts) // 250 * 250, max(p[0] for p in pts) // 250 * 250 + 250
    y0, y1 = min(p[1] for p in pts) // 250 * 250, max(p[1] for p in pts) // 250 * 250 + 250
    return rasterio.Affine(250, 0, x0, 0, -250, y1), (int((y1 - y0) / 250), int((x1 - x0) / 250))


def to_utm(arr, tf, crs, fill, dst_tf, shape):
    dst = np.full(shape, fill, arr.dtype)
    reproject(arr, dst, src_transform=tf, src_crs=crs, src_nodata=None, dst_transform=dst_tf, dst_crs=UTM,
              resampling=Resampling.nearest)
    return dst


def fetch_season(client, year):
    items = client.search(collections=["modis-13Q1-061"], bbox=BBOX, datetime=f"{year}-04-01/{year}-11-01").items()
    items = sorted((i for i in items if i.id.startswith("MOD13Q1.") and TILE in i.id), key=lambda i: i.id)
    with ThreadPoolExecutor(8) as ex:
        data = list(ex.map(read_composite, items))
    out = []
    for it, (ndvi, rel, tf, crs) in zip(items, data):
        doy = int(it.id.split(".")[1][5:])
        MANIFEST["sources"].append({
            "kind": "MODIS MOD13Q1 v061 16-day NDVI + pixel_reliability", "year": year, "item_id": it.id,
            "composite_start_doy": doy, "stac_datetime": it.properties.get("start_datetime") or str(it.datetime),
            "ndvi_url": clean(it.assets["250m_16_days_NDVI"].href),
            "reliability_url": clean(it.assets["250m_16_days_pixel_reliability"].href),
            "sha256_ndvi_window_int16": sha(ndvi.tobytes()), "sha256_reliability_window": sha(rel.tobytes()),
            "window_shape": list(ndvi.shape)})
        ndvi, rel = to_utm(ndvi, tf, crs, -3000, *utm_grid()), to_utm(rel, tf, crs, 3, *utm_grid())
        valid = (rel >= 0) & (rel <= 1) & (ndvi > -2000)
        out.append({"doy": doy, "date": dt.date(year, 1, 1) + dt.timedelta(days=doy - 1),
                    "ndvi": np.where(valid, ndvi * 1e-4, np.nan), "rel": rel})
    return out, *utm_grid()


def cropland_fraction(client, shape, tf):
    item = next(client.search(collections=["esa-worldcover"], bbox=BBOX, query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).items())
    href = item.assets["map"].href
    cached = CACHE / "worldcover_crop.npz"
    if not cached.exists():
        with rasterio.open(pc.sign(href)) as ds:
            win = from_bounds(*BBOX, ds.transform).round_offsets().round_lengths()
            lc = ds.read(1, window=win)
            stf = ds.window_transform(win)
        np.savez(cached, lc=lc, tf=np.array(stf)[:6])
    z = np.load(cached)
    lc, stf = z["lc"], rasterio.Affine(*z["tf"])
    frac = np.zeros(shape, np.float32)
    reproject((lc == 40).astype(np.float32), frac, src_transform=stf, src_crs="EPSG:4326",
              dst_transform=tf, dst_crs=UTM, resampling=Resampling.average)
    MANIFEST["sources"].append({"kind": "ESA WorldCover 10 m v200 (2021), cropland class 40 -> share per MODIS pixel",
                                "item_id": item.id, "url": clean(href), "license": "CC BY 4.0",
                                "sha256_window_uint8": sha(lc.tobytes()), "window_shape": list(lc.shape)})
    return frac


def fetch_power(year):
    url = ("https://power.larc.nasa.gov/api/temporal/daily/point?parameters=PRECTOTCORR&community=AG"
           f"&longitude={POWER_PT[0]}&latitude={POWER_PT[1]}&start={year}0401&end={year}1101&format=JSON")
    cached = CACHE / f"power_{year}.json"
    if not cached.exists():
        r = requests.get(url, timeout=60)
        r.raise_for_status()
        cached.write_bytes(r.content)
    raw = cached.read_bytes()
    MANIFEST["sources"].append({"kind": "NASA POWER daily PRECTOTCORR (mm/day), 0.5 deg grid cell", "year": year,
                                "url": url, "sha256_response": sha(raw)})
    series = json.loads(raw)["properties"]["parameter"]["PRECTOTCORR"]
    days = [dt.datetime.strptime(k, "%Y%m%d").date() for k in series]
    return days, np.array([max(v, 0) for v in series.values()])


def smooth(doys, vals):
    doys, vals = np.asarray(doys, float), np.asarray(vals, float)
    ok = ~np.isnan(vals)
    filled = np.interp(doys, doys[ok], vals[ok])
    return savgol_filter(filled, window_length=5, polyorder=2)


def season_stats(comps, crop):
    mask = crop >= 0.5
    rows = []
    for c in comps:
        clear = c["ndvi"][mask]
        share = float(np.mean(~np.isnan(clear)))
        rows.append({"doy": c["doy"], "date": c["date"], "clear": share,
                     "mean": float(np.nanmean(clear)) if share >= MIN_VALID else np.nan})
    sm = smooth([r["doy"] for r in rows], [r["mean"] for r in rows])
    for r, s in zip(rows, sm):
        r["sg"] = float(s)
    return rows


def pick_font():
    have = {f.name for f in font_manager.fontManager.ttflist}
    for name in ("Apple SD Gothic Neo", "AppleGothic", "Noto Sans CJK KR", "NanumGothic", "Malgun Gothic"):
        if name in have:
            return name
    print("warning: no Korean font found; Korean caption will show boxes", file=sys.stderr)
    return "DejaVu Sans"


def style(font):
    plt.rcParams.update({"font.family": [font, "DejaVu Sans"], "text.color": INK, "axes.labelcolor": INK,
                         "xtick.color": INK, "ytick.color": INK, "axes.edgecolor": INK, "figure.facecolor": PAPER,
                         "axes.facecolor": PAPER, "savefig.facecolor": PAPER, "axes.unicode_minus": False})


def draw_map(ax, comp, crop, tf, year, accent):
    h, w = crop.shape
    left, top = tf.c, tf.f
    ext = (left, left + w * tf.a, top + h * tf.e, top)
    rgb = np.ones((h, w, 3)) * np.array(matplotlib.colors.to_rgb("#EAE7E0"))
    ndvi, isc = comp["ndvi"], crop >= 0.5
    ok = isc & ~np.isnan(ndvi)
    rgb[ok] = NDVI_CMAP(np.clip((ndvi[ok] - 0.3) / 0.6, 0, 1))[:, :3]
    ax.imshow(rgb, extent=ext, interpolation="nearest")
    masked = isc & np.isnan(ndvi)
    ov = np.zeros((h, w, 4))
    ov[masked] = (*matplotlib.colors.to_rgb("#B9BBBE"), 1)
    ax.imshow(ov, extent=ext, interpolation="nearest")
    ys, xs = np.nonzero(masked)
    if len(ys):
        ax.contourf(np.linspace(ext[0], ext[1], w), np.linspace(ext[3], ext[2], h), masked.astype(int),
                    levels=[0.5, 1.5], colors="none", hatches=["////"])
    ax.set_xticks([]); ax.set_yticks([])
    for s in ax.spines.values():
        s.set_visible(False)
    km10 = 10_000
    x0, y0 = ext[0] + 0.05 * (ext[1] - ext[0]), ext[2] + 0.07 * (ext[3] - ext[2])
    ax.add_patch(Rectangle((x0 - 600, y0 - 500), km10 + 1200, 3300, color=PAPER, alpha=0.92, zorder=4, linewidth=0))
    ax.add_patch(Rectangle((x0, y0), km10, 450, color=INK, zorder=5))
    ax.text(x0 + km10 / 2, y0 + 1300, "10 km", ha="center", va="bottom", fontsize=9, color=INK, zorder=5)
    clear_share = float(np.mean(~np.isnan(ndvi[isc])))
    ax.set_title(f"{year}", loc="left", fontsize=17, fontweight="bold", color=accent)
    ax.set_title(f"{TAG[year]}  |  composite from {comp['date']:%b} {comp['date'].day}", loc="right", fontsize=8.5, color=MUTE)
    return float(np.nanmean(ndvi[isc])), clear_share


def draw_series(ax, rows, rain, year, accent):
    days, mm = rain
    d0 = dt.date(year, 4, 1)
    weeks = {}
    for d, v in zip(days, mm):
        weeks[(d - d0).days // 7] = weeks.get((d - d0).days // 7, 0) + v
    rx = ax.twinx()
    rx.bar([d0 + dt.timedelta(days=7 * k + 3) for k in weeks], list(weeks.values()), width=6, color=MUTE, alpha=0.35, linewidth=0)
    rx.set_ylim(0, 320); rx.invert_yaxis(); rx.set_yticks([0, 100, 200])
    rx.set_ylabel("rain mm/week", fontsize=8, color=MUTE); rx.tick_params(labelsize=8, colors=MUTE)
    ax.set_zorder(rx.get_zorder() + 1); ax.patch.set_visible(False)
    xs = [r["date"] for r in rows]
    ax.plot(xs, [r["sg"] for r in rows], color=accent, lw=2.6)
    okp = [r for r in rows if not np.isnan(r["mean"])]
    ax.scatter([r["date"] for r in okp], [r["mean"] for r in okp], s=14, color=accent, zorder=4)
    gap = [r for r in rows if np.isnan(r["mean"])]
    ax.scatter([r["date"] for r in gap], [r["sg"] for r in gap], s=22, facecolor=PAPER, edgecolor=accent, zorder=4)
    for name, d in EVENTS.get(year, []):
        ax.axvline(d, color=INK, lw=0.9, ls=":")
        ax.text(d, 0.93, name, rotation=90, va="top", ha="right", fontsize=8, color=INK)
    ax.set_ylim(0.3, 0.9); ax.set_xlim(d0, dt.date(year, 11, 1))
    ax.set_yticks([0.4, 0.6, 0.8]); ax.tick_params(labelsize=8)
    ax.set_xticks([dt.date(year, m, 1) for m in range(4, 12)])
    ax.set_xticklabels(["Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov"])
    ax.set_title(str(year), loc="left", fontsize=13, fontweight="bold", color=accent)
    for s in ("top",):
        ax.spines[s].set_visible(False); rx.spines[s].set_visible(False)
    rx.spines["right"].set_color(MUTE)


CREDIT = ("NDVI: NASA MODIS MOD13Q1 v061 (Terra, 250 m) via Microsoft Planetary Computer. Cropland mask: ESA WorldCover 2021 (CC BY 4.0). "
          "Rain: NASA POWER (0.5° cell). Grey hatch = cloud/low-quality pixels masked; non-cropland pixels pale.")


def add_cbar(fig, rect):
    cax = fig.add_axes(rect) if isinstance(rect, list) else fig.add_subplot(rect)
    cax.imshow(np.linspace(0, 1, 256)[None, :], aspect="auto", cmap=NDVI_CMAP, extent=(0.3, 0.9, 0, 1))
    cax.set_yticks([]); cax.set_xticks([0.3, 0.6, 0.9]); cax.tick_params(labelsize=8)
    cax.set_xlabel("NDVI", fontsize=9)
    for s in cax.spines.values():
        s.set_visible(False)


def main():
    client = pystac_client.Client.open(STAC)
    seasons, tf, shape = {}, None, None
    for y in (DAMAGE, NORMAL):
        seasons[y], tf, shape = fetch_season(client, y)
    crop = cropland_fraction(client, shape, tf)
    rain = {y: fetch_power(y) for y in (DAMAGE, NORMAL)}
    stats = {y: season_stats(seasons[y], crop) for y in seasons}
    font = pick_font()
    style(font)
    col = {DAMAGE: ACCENT, NORMAL: INK}
    comp = {y: next(c for c in seasons[y] if c["doy"] == MAP_DOY) for y in seasons}

    summary = {}
    figm, axm = plt.subplots(1, 2, figsize=(10, 5.4), gridspec_kw={"wspace": 0.04})
    for ax, y in zip(axm, (NORMAL, DAMAGE)):
        mean, clear = draw_map(ax, comp[y], crop, tf, y, col[y])
        summary[y] = {"map_mean_ndvi": mean, "map_clear_share": clear, "map_date": str(comp[y]["date"])}
    figm.suptitle("Gimje-Mangyeong rice plain, MODIS NDVI, mid-September composite", x=0.06, ha="left", fontsize=12, y=0.99)
    add_cbar(figm, [0.4, 0.15, 0.2, 0.03])
    figm.text(0.5, 0.02, CREDIT, ha="center", fontsize=6.3, color=MUTE)
    figm.subplots_adjust(top=0.86, bottom=0.27, left=0.03, right=0.97)
    figm.savefig(OUT / "ndvi_map.png", dpi=200)
    plt.close(figm)

    figs, axs = plt.subplots(1, 2, figsize=(10, 3.4), sharey=True, gridspec_kw={"wspace": 0.38})
    for ax, y in zip(axs, (NORMAL, DAMAGE)):
        draw_series(ax, stats[y], rain[y], y, col[y])
    axs[0].set_ylabel("NDVI (cropland mean, SG-smoothed)", fontsize=8)
    figs.text(0.5, 0.01, "Dots: clear composites. Open circles: composite >75% masked, value interpolated. Savitzky-Golay window 5, order 2. Dotted: typhoon landfall dates (secondary source).",
              ha="center", fontsize=6.3, color=MUTE)
    figs.subplots_adjust(top=0.88, bottom=0.17, left=0.07, right=0.93)
    figs.savefig(OUT / "ndvi_series.png", dpi=200)
    plt.close(figs)

    def val(y, doy):
        return next(r["sg"] for r in stats[y] if r["doy"] == doy)
    d_map = summary[DAMAGE]["map_mean_ndvi"] - summary[NORMAL]["map_mean_ndvi"]
    en = (f"Gimje-Mangyeong rice plain, MODIS NDVI, mid-Sept: {DAMAGE} (typhoons Bolaven, Sanba) "
          f"{summary[DAMAGE]['map_mean_ndvi']:.2f} vs {NORMAL} {summary[NORMAL]['map_mean_ndvi']:.2f} ({d_map:+.2f}); vegetation index only, not an insured loss.")
    ko = (f"김제·만경평야 MODIS NDVI, 9월 중순: {DAMAGE}년(태풍 볼라벤·산바) {summary[DAMAGE]['map_mean_ndvi']:.2f}, "
          f"{NORMAL}년 {summary[NORMAL]['map_mean_ndvi']:.2f} ({d_map:+.2f}). 식생지수 변화일 뿐 보험 손실이 아닙니다.")
    (OUT / "captions.txt").write_text(f"EN: {en}\nKO: {ko}\n", encoding="utf-8")

    fig = plt.figure(figsize=(10, 10.4))
    gs = fig.add_gridspec(4, 1, height_ratios=[5.2, 0.45, 3.0, 1.0], hspace=0.28, left=0.05, right=0.95, top=0.93, bottom=0.03)
    top = gs[0].subgridspec(1, 2, wspace=0.05)
    for j, y in enumerate((NORMAL, DAMAGE)):
        draw_map(fig.add_subplot(top[0, j]), comp[y], crop, tf, y, col[y])
    fig.suptitle("Trust under uncertainty: what a satellite sees in a typhoon year", x=0.05, ha="left", fontsize=14, fontweight="bold", y=0.975)
    cb = gs[1].subgridspec(1, 3, width_ratios=[2, 1, 2])
    add_cbar(fig, cb[0, 1])
    sub = gs[2].subgridspec(1, 2, wspace=0.38)
    sa = [fig.add_subplot(sub[0, j]) for j in range(2)]
    for ax, y in zip(sa, (NORMAL, DAMAGE)):
        draw_series(ax, stats[y], rain[y], y, col[y])
    sa[1].sharey(sa[0]); sa[0].set_ylabel("NDVI (SG-smoothed)", fontsize=8)
    fig.text(0.05, 0.088, en, fontsize=8.6, color=INK)
    fig.text(0.05, 0.063, ko, fontsize=8.6, color=INK)
    fig.text(0.05, 0.012, CREDIT, fontsize=6, color=MUTE)
    fig.savefig(OUT / "trust_sat_hero.png", dpi=200)
    plt.close(fig)

    MANIFEST["generated_utc"] = dt.datetime.now(dt.timezone.utc).isoformat(timespec="seconds")
    MANIFEST["script_sha256"] = sha(Path(__file__).read_bytes())
    MANIFEST["map_composite_doy"] = MAP_DOY
    MANIFEST["summary"] = summary
    MANIFEST["season_stats"] = {str(y): [{**r, "date": str(r["date"])} for r in rows] for y, rows in stats.items()}
    MANIFEST["sources"].sort(key=lambda s: (s["kind"], s.get("year", 0), s.get("item_id", "")))
    (ROOT / "manifest.json").write_text(json.dumps(MANIFEST, indent=1, default=lambda o: None if o != o else float(o)), encoding="utf-8")
    print(en); print(ko)
    for y in stats:
        print(y, " ".join(f"{r['doy']}:{r['mean']:.2f}/{r['clear']:.2f}" for r in stats[y]))


if __name__ == "__main__":
    main()