How did central Plains infrared cloud patterns change over fifteen minutes?

from datetime import UTC, datetime
from importlib.metadata import version
from pathlib import Path
from tempfile import TemporaryDirectory

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyproj import CRS, Transformer

from usdata import cite_lockfile, pull, verify

# One figure style for every usdata notebook, so previews look alike.
plt.rcParams.update(
    {
        "figure.figsize": (8, 4.5),
        "figure.dpi": 120,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.grid": True,
        "grid.alpha": 0.3,
        "font.size": 10,
    }
)
manifest = Path("dataset.yaml")
print("Executed (UTC):", datetime.now(UTC).isoformat(timespec="seconds"))
print(
    "; ".join(
        f"{name} {version(name)}"
        for name in ("usdata", "xarray", "h5netcdf", "numpy", "pandas", "pyproj")
    )
)

print(manifest.read_text().strip())

result = pull(manifest)
items = sorted(result.fetched, key=lambda item: item.asset.time.start)
assert items and all("CMIPM1-M6C13_G16" in item.asset.id for item in items)
print("\nScenes:", len(items), "· source bytes:", sum(item.provenance.size for item in items))
print("Restored from the committed lockfile:", result.from_lockfile)
print("First source:", items[0].asset.href)
pd.DataFrame(
    {
        "scan_start_utc": item.asset.time.start.strftime("%H:%M:%S.%f")[:-5],
        "bytes": item.provenance.size,
        "retrieved_utc": item.provenance.retrieved_at.isoformat(timespec="seconds"),
        "checksum": item.provenance.checksum[:23] + "...",
    }
    for item in items
).head()

starts = pd.DatetimeIndex([item.asset.time.start for item in items])
ends = pd.DatetimeIndex([item.asset.time.end for item in items])
start = pd.Timestamp("2024-05-06T22:00:00Z")
stop = pd.Timestamp("2024-05-06T22:15:00Z")
assert ((starts >= start) & (starts < stop)).all()
assert not starts.duplicated().any(), "Duplicate scan starts need review."
minute_counts = pd.Series(1, index=starts.floor("min")).groupby(level=0).sum()
expected = pd.date_range(start, stop, freq="min", inclusive="left")
missing = expected.difference(minute_counts.index)
print("Missing one-minute bins:", [time.strftime("%H:%M") for time in missing])
print("Scans per occupied minute:", sorted(minute_counts.unique().tolist()))
cadence = pd.DataFrame(
    {
        "start_utc": starts,
        "end_utc": ends,
        "since_previous_seconds": starts.to_series().diff().dt.total_seconds().to_numpy(),
    }
)
print(
    "Start intervals: "
    f"{cadence.since_previous_seconds.min():.1f} to {cadence.since_previous_seconds.max():.1f} s"
)
cadence

rows, frames, coverage = [], [], []
reference_x = reference_y = reference_projection = reference_extent = None
for item in items:
    with item.open() as scene:
        assert scene.CMI.attrs["units"] == "K" and scene.CMI.dims == ("y", "x")
        assert int(scene.band_id.item()) == 13 and scene.attrs["platform_ID"] == "G16"
        assert pd.Timestamp(scene.attrs["time_coverage_start"]) == pd.Timestamp(
            item.asset.time.start
        )
        assert "good_pixel_qf" in scene.DQF.attrs["flag_meanings"]
        projection = dict(scene.goes_imager_projection.attrs)
        extent = dict(scene.geospatial_lat_lon_extent.attrs)
        if reference_x is None:
            reference_x, reference_y = scene.x.values.copy(), scene.y.values.copy()
            reference_projection, reference_extent = projection, extent
            print("Quality flag values:", scene.DQF.attrs["flag_values"].tolist())
            print("Quality meanings:", scene.DQF.attrs["flag_meanings"])
            transformer = Transformer.from_crs(CRS.from_cf(projection), "EPSG:4326", always_xy=True)
            x, y = np.meshgrid(reference_x, reference_y)
            height = projection["perspective_point_height"]
            lon, lat = transformer.transform(x * height, y * height)
            region = (lon >= -100) & (lon <= -96) & (lat >= 34) & (lat <= 38)
            assert region.any()
            iy, ix = np.where(region)
            ys, xs = slice(iy.min(), iy.max() + 1), slice(ix.min(), ix.max() + 1)
        else:
            assert np.array_equal(scene.x.values, reference_x), "M1 moved: review comparison."
            assert np.array_equal(scene.y.values, reference_y), "M1 moved: review comparison."
            assert projection == reference_projection and extent == reference_extent
        values, flags = scene.CMI.values, scene.DQF.values
        good = region & (flags == 0) & np.isfinite(values)
        accepted = values[good]
        assert accepted.size > 0
        flag_counts = {
            str(int(flag)): int(np.count_nonzero(region & (flags == flag)))
            for flag in np.unique(flags[np.isfinite(flags)])
        }
        rows.append(
            {
                "time": item.asset.time.start,
                "region_pixels": int(region.sum()),
                "accepted_pixels": int(accepted.size),
                "excluded_pixels": int(region.sum() - accepted.size),
                "DQF_counts": flag_counts,
                "missing_DQF": int(np.count_nonzero(region & ~np.isfinite(flags))),
                "p10_K": float(np.quantile(accepted, 0.1)),
                "below_235K_percent": float(100 * np.mean(accepted < 235)),
            }
        )
        coverage.append(
            {
                "time": item.asset.time.start,
                **{
                    key.removeprefix("geospatial_"): float(extent[key])
                    for key in (
                        "geospatial_westbound_longitude",
                        "geospatial_eastbound_longitude",
                        "geospatial_southbound_latitude",
                        "geospatial_northbound_latitude",
                    )
                },
            }
        )
        frames.append(np.where(good, values, np.nan)[ys, xs])
summary = pd.DataFrame(rows).set_index("time")
print("Projection, coordinates, and metadata extent match across all scans.")
print(pd.DataFrame(coverage).head(3).to_string(index=False))
summary[["region_pixels", "accepted_pixels", "excluded_pixels", "p10_K", "below_235K_percent"]]

selected = sorted(set([0, len(items) // 2, len(items) - 1]))
fig, axes = plt.subplots(
    1, len(selected), sharex=True, sharey=True, layout="constrained", squeeze=False
)
for axis, index in zip(axes[0], selected, strict=True):
    mesh = axis.pcolormesh(
        reference_x[xs],
        reference_y[ys],
        frames[index],
        shading="auto",
        cmap="magma",
        vmin=195,
        vmax=295,
    )
    axis.grid(False)
    axis.tick_params(labelsize=7)
    axis.xaxis.set_major_locator(plt.MaxNLocator(3))
    axis.set_title(f"{items[index].asset.time.start:%H:%M:%S} UTC", fontsize=9)
    axis.set_xlabel("East-west scan angle (rad)", fontsize=8)
axes[0, 0].set_ylabel("North-south scan angle (rad)", fontsize=8)
fig.colorbar(
    mesh,
    ax=axes.ravel().tolist(),
    label="Brightness temperature (K)",
    shrink=0.8,
    location="bottom",
)
fig.suptitle("GOES-16 M1 channel 13, 6 May 2024, 34-38°N 100-96°W (DQF = 0)")
plt.show()
plt.close(fig)

fig, axes = plt.subplots(2, 1, sharex=True, layout="constrained")
axes[0].scatter(summary.index, summary.p10_K, color="#23657b", s=24)
axes[0].set_ylabel("10th percentile (K)")
axes[1].scatter(summary.index, summary.below_235K_percent, color="#98442c", s=24)
axes[1].set_ylabel("Pixels <235 K (%)")
axes[1].set_xlabel("Scan-start time (UTC)")
axes[1].xaxis.set_major_formatter(mdates.DateFormatter("%H:%M", tz=UTC))
for axis in axes:
    for time in missing:
        axis.axvspan(time, time + pd.Timedelta(minutes=1), color="gray", alpha=0.2)
fig.suptitle("GOES-16 M1 channel 13 over 34-38°N 100-96°W, 6 May 2024")
plt.show()
plt.close(fig)
print(f"10th percentile: {summary.p10_K.iloc[0]:.2f} → {summary.p10_K.iloc[-1]:.2f} K")
print(
    f"Cold-pixel fraction: {summary.below_235K_percent.iloc[0]:.2f} → "
    f"{summary.below_235K_percent.iloc[-1]:.2f}%"
)

assert verify(manifest) == []
for citation in cite_lockfile(manifest):
    print(citation.as_text())

original = {item.asset.id: item for item in items}
with TemporaryDirectory(prefix="usdata-mesoscale-restore-") as fresh:
    restored = pull(manifest, root=fresh)
    assert restored.from_lockfile and not any(item.from_cache for item in restored.fetched)
    assert restored.lockfile == result.lockfile
    assert {item.asset.id for item in restored.fetched} == set(original)
    assert verify(manifest, root=fresh) == []
    for item in restored.fetched:
        assert item.path.read_bytes() == original[item.asset.id].path.read_bytes()
print(f"All {len(items)} pinned scenes restored byte-for-byte into an empty cache.")