From b6fbd208726f14833ecfc8f6a1cc2773d6027b5b Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 14 Sep 2026 17:33:35 -0500 Subject: [PATCH] ADD: SVRIMG convective mode benchmark (Cellular vs QLCS) Adds an external, human-labelled benchmark so the labeller can be scored against published ground truth rather than only in-house CSAPR2 labels. SVRIMG (Haberlie, Ashley and Karpinski 2020) provides 1,741 hand-classified GridRad/NEXRAD scenes centred on SPC severe reports. Collapsed to the two convective mode classes, Cellular and QLCS, that is 1,315 scenes across the published split: 997 train, 88 validation, 230 test. lars/preprocessing/svrimg.py mirrors preprocess_radar_data's CSV schema, so the existing inference, Dawid-Skene and kappa tooling runs unchanged; 'unid' and 'split' are added. Two properties of the upstream format drive the implementation: - The distributed PNGs are palette-mode images whose palette is a dummy linear ramp, so they look greyscale in an ordinary viewer. The palette indices are the data. read_svrimg_png reads the indices and refuses any non-palette file, so a silently RGB-converted image fails loudly instead of producing plausible-looking nonsense. - Those indices are reflectivity in dBZ, integer-truncated and floored at zero (upstream GridRad REFC_MAX is itself uint8 with units dBZ). There are no negative values, hence vmin=0 rather than the -20 used for CSAPR2. Thresholds are calibrated on the training split only, not carried over from CODEBOOK_NWSRef.md, whose area fractions were tuned for 1 km CSAPR2 data on a 2048x2048 polar grid and do not transfer to 3.75 km column-max on a 136x136 Cartesian grid. pct_gates_30dBZ is the strongest single discriminator at AUC 0.935, Youden cut 6.48 percent, adopted as 6.5. pct_gates_20/10/40dBZ also separate (AUC 0.929/0.903/0.889) but are omitted to avoid stacking redundant correlated reclassification rules; pct_gates_50dBZ does not discriminate (AUC 0.601). Caveat recorded in the codebook changelog: these coverage fractions largely measure echo extent rather than linearity, so they are expected to misfire on large non-linear convective clusters. The codebook describes classes in reflectivity and topology rather than colour, because COLOR_DBZ_RANGE in lars/nepho/inference.py is calibrated to a ChaseSpectral-like ramp that does not match the NWSRef bands these scenes render with. Note for reviewers: the published split carries a prior shift (train 54% Cellular, test 66% QLCS), so report balanced accuracy and per-class recall rather than raw accuracy. Co-Authored-By: Claude Opus 5 (1M context) --- example_codebooks/CODEBOOK_SVRIMG.md | 146 +++++++ lars/preprocessing/__init__.py | 3 +- lars/preprocessing/labels.py | 3 + lars/preprocessing/svrimg.py | 374 ++++++++++++++++++ .../reports/calibrate_svrimg_thresholds.py | 166 ++++++++ .../reports/svrimg_threshold_calibration.csv | 7 + tests/test_svrimg.py | 278 +++++++++++++ 7 files changed, 976 insertions(+), 1 deletion(-) create mode 100644 example_codebooks/CODEBOOK_SVRIMG.md create mode 100644 lars/preprocessing/svrimg.py create mode 100644 notebooks/reports/calibrate_svrimg_thresholds.py create mode 100644 notebooks/reports/svrimg_threshold_calibration.csv create mode 100644 tests/test_svrimg.py diff --git a/example_codebooks/CODEBOOK_SVRIMG.md b/example_codebooks/CODEBOOK_SVRIMG.md new file mode 100644 index 0000000..3ad50c3 --- /dev/null +++ b/example_codebooks/CODEBOOK_SVRIMG.md @@ -0,0 +1,146 @@ +# Radar Image Labelling Codebook + +A reference guide for annotators labelling radar imagery with LARS. Use this codebook to ensure consistent, reproducible labels across all annotators and sessions. + +--- + +## 1. Overview + +The purpose of this section is to label radar imagery by mode of convection + +- **Radar type:** NEXRAD WSR-88D (GridRad three-dimensional merged mosaic) +- **Data source:** SVRIMG, scenes centred on SPC severe weather reports +- **Geographic scope:** Contiguous United States, 1996-2017 +- **Labelling task:** convective mode classification + +--- + +## 2. Data Description + +### 2.1 Input Fields + +| Field | Units | Description | +|--------------------------|-------|-------------------------------------| +| *reflectivity* | dBZ | Column-maximum intensity of returned radar signal | +| n_gates_30dBZ | percent | The percentage of gates greater than 30 dBZ | + +### 2.2 Image Format + +- **Spatial resolution:** 3.75 km by 3.75 km +- **Temporal resolution:** one scene per severe report, matched to the nearest hourly volume +- **Projection:** Lambert conformal conic, 136x136 grid covering approximately 512 km by 512 km +- **Color scale:** NWSRef colormap with vmin=0 and vmax=80 + +Source values are column-maximum reflectivity stored as integer dBZ and floored +at zero, so no negative reflectivity appears in these scenes. For reference, the +rendered colour bands are: + +| Colour | Approximate reflectivity | +|--------|--------------------------| +| cyan to blue | 5-20 dBZ | +| green | 20-35 dBZ | +| yellow to orange | 35-50 dBZ | +| red | 50-65 dBZ | +| magenta to purple | 65-80 dBZ | + +--- + +## 3. Label Classes + +Each image or region-of-interest must be assigned exactly one primary class. + +### 3.1 Primary Classes + +| Label | Description | +|------------------------|-----------------------------------------------------------------------------| +| Cellular | One or more discrete convective cores, each bounded on every side by weak or absent echo. Cores do not merge into a continuous band, and any echo bridging them is markedly weaker than the cores themselves. Most of the frame is free of echo. The percentage of gates greater than 30 dBZ must not exceed 6.5 percent. If it does exceed 6.5 percent, then classify as QLCS. | +| QLCS | A single continuous, quasi-linear band of strong echo at least roughly 100 km long, whose length is about three times its width or more. The band may be bowed or curved, and may trail a broad region of weaker stratiform echo behind it. Short breaks along the band are acceptable provided the overall structure still reads as one line rather than a field of separate cells. | +| Ambiguous / Uncertain | Cannot be classified with confidence. Use for transitional structures that are neither clearly one line nor clearly separate cells, and for scenes degraded by missing radar coverage. | + +> **Note on enforcement.** Only the quantitative reflectivity thresholds in the +> descriptions above (e.g. "percentage of gates greater than X dBZ must not +> exceed Y percent") are validated automatically, and only when the +> corresponding `pct_gates_*` / `n_gates_*` columns are present in the data. +> Spatial and topological criteria — "bounded on every side by…", "continuous +> band", "length is about three times its width" — are judged by the annotator +> or model and are **not** checked programmatically. + +> **Note on colour criteria.** This codebook deliberately describes classes in +> reflectivity and topology rather than colour. `COLOR_DBZ_RANGE` in +> `lars/nepho/inference.py` is calibrated against a ChaseSpectral-like ramp +> (green 10-30, yellow 30-40, red 40-50, pink 50+), which does not match the +> NWSRef bands tabulated in Section 2.2. Colour phrasings here would therefore +> be validated against the wrong reflectivity bands. + +--- + +## 5. Labelling Procedure + +1. Use :code:`lars.preprocessing.preprocess_svrimg_data` to download scenes and generate images and a .csv file +2. The csv file is populated with the published SVRIMG hand labels in the 'label' column, and a 'split' column giving the published train/validation/test assignment. +3. According to the criteria above, label all images in the 'file_path' column of the .csv file. + +--- + +## 6. Annotator Guidelines + +The bullets in this section are passed verbatim to automated labelling models, +so they must be self-contained for a single image with no external context. + +- Every image is centred on a confirmed severe weather report, so an empty scene is never the right answer; if echo is present but its organisation is unclear, classify as Ambiguous / Uncertain. +- Judge the organisation of the whole image, not only the structures at its centre. +- If a line and separate cells are both present, classify according to whichever structure covers the larger area of the image. +- A broad region of weak echo on its own is not a class here; classify by the arrangement of the strong echo embedded within it. + +### 6.1 Human Annotators Only + +These apply to human annotators and the review process. They are intentionally +kept out of the bullet list above because an automated model labels each image +independently, with no temporal context and no access to the example gallery. + +- Use the provided example gallery (Section 8) to calibrate your judgement. +- Inter-annotator agreement should be checked periodically; raise disagreements with the team lead. + +--- + +## 7. Quality Control + +| Check | Method | +|-------|--------| +| Completeness | All images have a primary label | +| Consistency | Random sample reviewed by second annotator | +| Agreement metric | Cohen's κ computed per annotator pair | +| Class balance | Report balanced accuracy and per-class recall; the published split carries a prior shift (train 54% Cellular, test 66% QLCS) | +| Outlier review | Labels deviating from model predictions flagged for review | + +--- + +## 8. Example Gallery + +Representative training scenes, chosen as the cases nearest each class median +for 30 dBZ coverage. Render them with +:code:`lars.preprocessing.preprocess_svrimg_data`; the source files are +palette-encoded and appear greyscale if opened directly. + +| Class | Representative UNID | pct_gates_30dBZ | Notes | +|--------------------------|-----------------------------------|--------|--------------------------| +| Cellular | 200207250114z000000223 | 3.31 | Class median; discrete cores, mostly echo-free frame | +| Cellular | 200805282245z000000948 | 3.30 | Intense core (71 dBZ) yet still low total coverage | +| QLCS | 200210291000z000000352 | 14.33 | Class median; continuous band | +| QLCS | 201104042340z000292240 | 14.33 | Band with trailing stratiform region | + +--- + +## 9. Changelog + +| Version | Date | Author | Changes | +|---------|------|--------|---------| +| 1.0 | 2026-09-14 | Robert Jackson | Initial release. Adapts CODEBOOK_NWSRef.md to the SVRIMG benchmark: binary Cellular/QLCS collapse (1,315 of 1,741 labelled scenes), vmin 0 rather than -20 because SVRIMG reflectivity is floored at zero, and colour criteria replaced by reflectivity/topology criteria. Thresholds calibrated on the 997-scene training split (see notebooks/reports/calibrate_svrimg_thresholds.py): pct_gates_30dBZ is the strongest single discriminator at AUC 0.935, Youden cut 6.48 percent (sensitivity 0.89, specificity 0.82), adopted as 6.5 percent. Also measured but not adopted, to avoid redundant correlated rules: pct_gates_20dBZ AUC 0.929, pct_gates_10dBZ AUC 0.903, pct_gates_40dBZ AUC 0.889. pct_gates_50dBZ does not discriminate (AUC 0.601) and pct_gates_60dBZ is weak and degenerate (AUC 0.684 at a 0.0 percent cut); neither is used. Caveat: these coverage fractions largely measure echo *extent* rather than linearity, so they are expected to misfire on large non-linear convective clusters. | + +--- + +## 10. References + +- Haberlie, A. M., W. S. Ashley, and M. Karpinski (2020). Mean storms: Composites of radar reflectivity images during two decades of severe thunderstorm events. *International Journal of Climatology*. +- Bowman, K. P., and C. R. Homeyer (2017). GridRad - Three-Dimensional Gridded NEXRAD WSR-88D Radar Data. Research Data Archive at NCAR. +- American Meteorological Society Glossary: https://glossary.ametsoc.org diff --git a/lars/preprocessing/__init__.py b/lars/preprocessing/__init__.py index 3720720..18e28cb 100644 --- a/lars/preprocessing/__init__.py +++ b/lars/preprocessing/__init__.py @@ -1,2 +1,3 @@ from .radar_preprocessing import preprocess_radar_data # noqa: F401 -from .labels import load_labels, save_labels, change_file_path, copy_labels, apply_criteria_to_labels, combine_labels, standardize_labels # noqa: F401 \ No newline at end of file +from .labels import load_labels, save_labels, change_file_path, copy_labels, apply_criteria_to_labels, combine_labels, standardize_labels # noqa: F401 +from .svrimg import (fetch_svrimg_labels, preprocess_svrimg_data, read_svrimg_png, reflectivity_stats, split_for_year) # noqa: F401 diff --git a/lars/preprocessing/labels.py b/lars/preprocessing/labels.py index 1b62f64..e997b0e 100644 --- a/lars/preprocessing/labels.py +++ b/lars/preprocessing/labels.py @@ -8,6 +8,9 @@ "unknown": "Ambiguous / Uncertain", "stratiform": "Stratiform Precipitation", "ambiguous / uncertain": "Ambiguous / Uncertain", + # SVRIMG convective mode classes (see example_codebooks/CODEBOOK_SVRIMG.md). + "cellular": "Cellular", + "qlcs": "QLCS", } diff --git a/lars/preprocessing/svrimg.py b/lars/preprocessing/svrimg.py new file mode 100644 index 0000000..17ed7a1 --- /dev/null +++ b/lars/preprocessing/svrimg.py @@ -0,0 +1,374 @@ +"""Ingest for the SVRIMG severe convective mode benchmark. + +SVRIMG (Haberlie, Ashley and Karpinski 2020) is a collection of hand-labelled +NEXRAD/GridRad column-maximum reflectivity images centred on SPC severe +weather reports. This module downloads those labels and images and emits them +in the same CSV schema as :func:`lars.preprocessing.preprocess_radar_data`, so +the existing inference and inter-annotator agreement tooling runs against +external ground truth without modification. + +Two properties of the upstream format drive the implementation: + +* The distributed PNGs are palette-mode images whose palette is a dummy linear + ramp (index ``i`` maps to RGB ``(3i, 3i+1, 3i+2)``). The *palette indices* + carry the data, which is why the files render as greyscale in an ordinary + image viewer. :func:`read_svrimg_png` reads the indices, not the RGB. +* Those indices are reflectivity in dBZ, integer-truncated and floored at zero. + The upstream GridRad ``REFC_MAX`` field is itself stored as ``uint8`` with + ``units: dBZ``. There are no negative values, so the natural colour scale + bounds are ``vmin=0``, ``vmax=80`` rather than the ``vmin=-20`` used for + CSAPR2 imagery. + +References +---------- +Haberlie, A. M., W. S. Ashley, and M. Karpinski, 2020: Mean storms: Composites +of radar reflectivity images during two decades of severe thunderstorm events. +*International Journal of Climatology*. +""" + +import os + +import numpy as np +import pandas as pd +import requests +from PIL import Image + +SVRIMG_BASE_URL = "https://nimbus.niu.edu/svrimg/data" + +#: Canonical hand-classification table published upstream. +CLASSIFICATION_CSV = "sample_classifications_96-17.csv" + +#: Grid is 136 x 136 at 3.75 km, i.e. roughly a 512 x 512 km box. +GRID_SIZE_PX = 136 + +#: Published train/validation/test split, keyed by inclusive year range. +SPLIT_YEARS = { + "train": (1996, 2011), + "val": (2012, 2013), + "test": (2014, 2017), +} + +#: The two convective mode classes retained by the binary collapse. The +#: upstream taxonomy also contains Tropical, Other, Noise and Missing. +BINARY_CLASSES = ("Cellular", "QLCS") + +_DOWNLOAD_TIMEOUT = 120 + + +def split_for_year(year): + """ + Return the published SVRIMG split a given year belongs to. + + Parameters + ---------- + year (int): Four-digit year. + + Returns + ------- + str or None + 'train', 'val' or 'test', or None if the year falls outside the + 1996-2017 range covered by the dataset. + """ + for split, (first, last) in SPLIT_YEARS.items(): + if first <= year <= last: + return split + return None + + +def _download(url, dest): + """ + Download ``url`` to ``dest`` unless it already exists. + + Parameters + ---------- + url (str): Source URL. + dest (str): Destination path. Parent directories are created as needed. + + Returns + ------- + str + The destination path. + """ + if os.path.exists(dest): + return dest + + os.makedirs(os.path.dirname(dest) or ".", exist_ok=True) + response = requests.get(url, timeout=_DOWNLOAD_TIMEOUT, stream=True) + response.raise_for_status() + + # Write to a temporary name first so an interrupted download does not + # leave a truncated file that later runs would treat as cached. + tmp = dest + ".part" + with open(tmp, "wb") as f: + for chunk in response.iter_content(chunk_size=65536): + f.write(chunk) + os.replace(tmp, dest) + + return dest + + +def read_svrimg_png(path): + """ + Read an SVRIMG PNG and return its reflectivity field in dBZ. + + SVRIMG PNGs are palette-mode images in which the palette is a meaningless + linear ramp and the palette *indices* are the data. Converting such a file + to RGB, as most image readers do by default, destroys the values. This + function reads the indices directly and refuses any file that is not + palette-mode, so a silently converted image fails loudly rather than + producing plausible-looking nonsense. + + Parameters + ---------- + path (str): Path to an SVRIMG .png file. + + Returns + ------- + numpy.ndarray + A (136, 136) uint8 array of reflectivity in dBZ. Zero means no echo. + + Raises + ------ + ValueError + If the image is not palette-mode. + """ + with Image.open(path) as im: + if im.mode != "P": + raise ValueError( + f"Expected a palette-mode ('P') SVRIMG PNG, got mode " + f"'{im.mode}' for {path}. Reading this file as RGB would " + f"discard the reflectivity values." + ) + return np.array(im, dtype=np.uint8) + + +def reflectivity_stats(reflectivity, dbz_thresholds): + """ + Compute reflectivity coverage statistics for a single SVRIMG scene. + + Mirrors the statistics that :func:`preprocess_radar_data` derives from + CSAPR2 sweeps, so the resulting columns are directly comparable. + + Parameters + ---------- + reflectivity (numpy.ndarray): Reflectivity in dBZ, zero meaning no echo. + dbz_thresholds (sequence of float): Thresholds for which gate counts and + coverage percentages are computed. + + Returns + ------- + dict + Mapping with keys 'ref_min', 'ref_max', and one + 'n_gates_dbz'/'pct_gates_dbz' pair per threshold. 'ref_min' and + 'ref_max' are computed over echoing gates only and are NaN for a scene + with no echo at all. + """ + reflectivity = np.asarray(reflectivity) + echo = reflectivity[reflectivity > 0] + + stats = { + "ref_min": float(echo.min()) if echo.size else float("nan"), + "ref_max": float(echo.max()) if echo.size else float("nan"), + } + + total_gates = reflectivity.size + for threshold in dbz_thresholds: + count = int(np.sum(reflectivity > threshold)) + stats[f"n_gates_{threshold}dbz"] = count + stats[f"pct_gates_{threshold}dbz"] = round(count / total_gates * 100, 4) + + return stats + + +def fetch_svrimg_labels(data_dir, haz_type="tor", classes=BINARY_CLASSES, + base_url=SVRIMG_BASE_URL): + """ + Download the SVRIMG hand classifications and join them to scene metadata. + + Parameters + ---------- + data_dir (str): Directory in which downloaded tables are cached. + haz_type (str): Hazard index to join against; 'tor', 'hail' or 'wind'. + Default 'tor'. + classes (sequence of str or None): Class names to retain. Defaults to + ('Cellular', 'QLCS'), the binary convective mode collapse. Pass None + to keep all six upstream classes. + base_url (str): Base URL for the SVRIMG data server. + + Returns + ------- + pd.DataFrame + Indexed by 'unid', with columns 'label', 'radar_time', 'report_time', + 'split', and the upstream scene descriptors 'area', 'convection_area', + 'intense_area', 'max_intensity' and 'mean_intensity'. + """ + class_csv = _download( + f"{base_url}/{CLASSIFICATION_CSV}", + os.path.join(data_dir, CLASSIFICATION_CSV), + ) + index_name = f"96-17_{haz_type}_utc_svrimg_index.csv" + index_csv = _download( + f"{base_url}/{index_name}", os.path.join(data_dir, index_name) + ) + + labels = pd.read_csv(class_csv, index_col="UNID") + labels = labels.rename(columns={"Class Name": "label"}) + labels.index.name = "unid" + + if classes is not None: + labels = labels[labels["label"].isin(classes)] + + index = pd.read_csv(index_csv, index_col="unid") + keep = [ + "radar_time", "report_time", "area", "convection_area", + "intense_area", "max_intensity", "mean_intensity", + ] + joined = labels[["label"]].join(index[keep], how="inner") + + joined["radar_time"] = pd.to_datetime( + joined["radar_time"], format="%m/%d/%Y %H:%M" + ) + joined["report_time"] = pd.to_datetime( + joined["report_time"], format="%m/%d/%Y %H:%M" + ) + joined["split"] = joined["radar_time"].dt.year.map(split_for_year) + + return joined.sort_values("radar_time") + + +def svrimg_image_url(unid, haz_type="tor", base_url=SVRIMG_BASE_URL): + """ + Build the upstream URL for a single SVRIMG scene. + + Parameters + ---------- + unid (str): SVRIMG unique id, e.g. '199601180433z000000018'. The leading + 12 characters are YYYYMMDDHHMM. + haz_type (str): 'tor', 'hail' or 'wind'. Default 'tor'. + base_url (str): Base URL for the SVRIMG data server. + + Returns + ------- + str + Fully qualified URL of the scene's .png file. + """ + year, month = unid[:4], unid[4:6] + return f"{base_url}/{haz_type}/{year}/{month}/{unid}.png" + + +def preprocess_svrimg_data(output_path, data_dir, haz_type="tor", + classes=BINARY_CLASSES, splits=None, + dbz_thresholds=(10, 20, 30, 40, 50), + size_px=256, dpi=150, limit=None, + base_url=SVRIMG_BASE_URL, **kwargs): + """ + Render SVRIMG scenes to labelled .png images for model evaluation. + + Downloads the SVRIMG classifications and the corresponding scenes, decodes + each scene's reflectivity field, and re-renders it with a meteorological + colormap so that a vision model sees imagery consistent with the codebook. + Gates with no echo are masked rather than drawn as 0 dBZ. + The returned DataFrame uses the same schema as + :func:`preprocess_radar_data`, with 'unid' and 'split' added. + + Parameters + ---------- + output_path (str): Directory in which rendered .png images are written. + data_dir (str): Directory used to cache downloaded tables and source PNGs. + haz_type (str): 'tor', 'hail' or 'wind'. Default 'tor'. + classes (sequence of str or None): Class names to retain. Default is the + binary ('Cellular', 'QLCS') collapse. + splits (sequence of str or None): Restrict to these published splits, e.g. + ('test',). Default None keeps all. + dbz_thresholds (sequence of float): Reflectivity thresholds in dBZ for + which gate counts and coverage percentages are computed. + size_px (int): Width and height of the output PNG in pixels. Default 256. + dpi (int): Dots per inch for the saved figure. Default 150. + limit (int or None): Process at most this many scenes. Useful for smoke + tests; default None processes everything. + base_url (str): Base URL for the SVRIMG data server. + + **kwargs: + + Additional keyword arguments are passed to matplotlib's imshow. 'cmap', + 'vmin' and 'vmax' default to 'NWSRef', 0 and 80 respectively, matching + CODEBOOK_SVRIMG.md. + + Returns + ------- + label_df: pd.DataFrame + DataFrame containing labels, paths and times, indexed by 'time'. Note + that several reports can share one hourly radar volume, so the index + is not unique; 'unid' and 'file_path' are the unique keys. + """ + import matplotlib.pyplot as plt + + kwargs.setdefault("vmin", 0) + kwargs.setdefault("vmax", 80) + kwargs.setdefault("cmap", "NWSRef") + if "ax" in kwargs: + raise ValueError("Do not pass in an axis to this function.") + + labels = fetch_svrimg_labels( + data_dir, haz_type=haz_type, classes=classes, base_url=base_url + ) + if splits is not None: + labels = labels[labels["split"].isin(splits)] + if limit is not None: + labels = labels.iloc[:limit] + + os.makedirs(output_path, exist_ok=True) + + dbz_thresholds = list(dbz_thresholds) + count_cols = [f"n_gates_{t}dbz" for t in dbz_thresholds] + pct_cols = [f"pct_gates_{t}dbz" for t in dbz_thresholds] + columns = (["file_path", "time", "unid", "split", "label", + "ref_min", "ref_max"] + count_cols + pct_cols) + rows = [] + + for unid, row in labels.iterrows(): + source = os.path.join( + data_dir, haz_type, unid[:4], unid[4:6], f"{unid}.png" + ) + try: + _download(svrimg_image_url(unid, haz_type, base_url), source) + reflectivity = read_svrimg_png(source) + except (requests.HTTPError, ValueError) as err: + print(f"Skipping {unid}: {err}") + continue + + stats = reflectivity_stats(reflectivity, dbz_thresholds) + + # Zero means "no echo", not "0 dBZ". Mask it so those gates render as + # blank background rather than as the bottom colour of the ramp, which + # would fill the frame with apparent echo. This mirrors the + # `where(field > min_ref)` masking in preprocess_radar_data. + display = np.where(reflectivity > 0, reflectivity, np.nan) + + fig = plt.figure(figsize=(size_px / dpi, size_px / dpi)) + ax = plt.axes() + # Row 0 of the stored array is the southern edge, so flip to draw + # north-up. This matches the upstream svrimg draw_box_plot helper. + ax.imshow(np.flipud(display), **kwargs) + ax.axis("off") + ax.set_title("") + file_name = os.path.join(output_path, f"{unid}.png") + fig.savefig(file_name, dpi=dpi, bbox_inches="tight", pad_inches=0) + plt.close(fig) + + time_str = row["radar_time"].strftime("%Y-%m-%d %H:%M:%S") + rows.append( + { + "file_path": file_name, + "time": time_str, + "unid": unid, + "split": row["split"], + "label": row["label"], + **stats, + } + ) + + out_df = pd.DataFrame(rows, columns=columns) + out_df.set_index("time", inplace=True) + + return out_df.sort_index() diff --git a/notebooks/reports/calibrate_svrimg_thresholds.py b/notebooks/reports/calibrate_svrimg_thresholds.py new file mode 100644 index 0000000..701c96e --- /dev/null +++ b/notebooks/reports/calibrate_svrimg_thresholds.py @@ -0,0 +1,166 @@ +"""Calibrate reflectivity-coverage thresholds for CODEBOOK_SVRIMG.md. + +Cellular vs QLCS is a morphological distinction -- discrete, mutually +separated cores versus a single connected quasi-linear band -- so it is an +open question whether scalar reflectivity coverage fractions discriminate the +two classes at all. This script answers that empirically rather than assuming +either way. + +For each candidate dBZ threshold it computes the coverage fraction over the +SVRIMG *training* split only (1996-2011), then reports ROC AUC and the +Youden-optimal cut point for Cellular vs QLCS. A threshold earns a place in +the codebook only if it clears AUC_GATE; everything else is reported here and +recorded in the codebook changelog as a measured null result. + +Run directly: + + python notebooks/reports/calibrate_svrimg_thresholds.py + +The first run downloads roughly 1,000 scenes (~11 MB) and caches them. +""" +import sys +from unittest.mock import MagicMock +sys.modules.setdefault("pip_system_certs", MagicMock()) +sys.modules.setdefault("pip_system_certs.wrapt_requests", MagicMock()) + +import os + +import numpy as np +import pandas as pd +from sklearn.metrics import roc_auc_score, roc_curve + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..")) + +from lars.preprocessing.svrimg import ( # noqa: E402 + fetch_svrimg_labels, read_svrimg_png, reflectivity_stats, + svrimg_image_url, _download, +) + +SVRIMG_DIR = os.environ.get( + "SVRIMG_DATA_DIR", os.path.expanduser("~/lars_examples/svrimg") +) +OUT_CSV = os.path.join(os.path.dirname(__file__), "svrimg_threshold_calibration.csv") + +DBZ_THRESHOLDS = (10, 20, 30, 40, 50, 60) + +#: Minimum ROC AUC for a threshold to be written into the codebook. +AUC_GATE = 0.65 + +#: QLCS is the positive class throughout. +POSITIVE = "QLCS" + + +def build_feature_table(split="train", haz_type="tor"): + """ + Decode every scene in a split and return its coverage fractions. + + Parameters + ---------- + split (str): Published SVRIMG split to use. Default 'train'. + haz_type (str): Hazard index. Default 'tor'. + + Returns + ------- + pd.DataFrame + One row per scene, with a 'label' column and one + 'pct_gates_dbz' column per threshold in DBZ_THRESHOLDS. + """ + labels = fetch_svrimg_labels(SVRIMG_DIR, haz_type=haz_type) + labels = labels[labels["split"] == split] + print(f"{split}: {len(labels)} scenes " + f"({labels['label'].value_counts().to_dict()})") + + rows = [] + for n, (unid, row) in enumerate(labels.iterrows(), start=1): + dest = os.path.join( + SVRIMG_DIR, haz_type, unid[:4], unid[4:6], f"{unid}.png" + ) + try: + _download(svrimg_image_url(unid, haz_type), dest) + reflectivity = read_svrimg_png(dest) + except Exception as err: # noqa: BLE001 + print(f" skipping {unid}: {err}") + continue + + stats = reflectivity_stats(reflectivity, DBZ_THRESHOLDS) + rows.append({"unid": unid, "label": row["label"], **stats}) + + if n % 200 == 0: + print(f" ...{n}/{len(labels)}") + + return pd.DataFrame(rows).set_index("unid") + + +def evaluate(features): + """ + Score each coverage fraction as a Cellular/QLCS discriminator. + + Parameters + ---------- + features (pd.DataFrame): Output of :func:`build_feature_table`. + + Returns + ------- + pd.DataFrame + One row per threshold with ROC AUC, the orientation of the + relationship, the Youden-optimal cut point and its sensitivity and + specificity, and whether the threshold clears AUC_GATE. + """ + y = (features["label"] == POSITIVE).astype(int).values + results = [] + + for threshold in DBZ_THRESHOLDS: + column = f"pct_gates_{threshold}dbz" + score = features[column].values + + auc = roc_auc_score(y, score) + # An AUC below 0.5 just means coverage runs the other way (higher + # coverage implies Cellular). Orient it so magnitudes compare. + oriented = max(auc, 1 - auc) + direction = "higher=QLCS" if auc >= 0.5 else "higher=Cellular" + + fpr, tpr, cuts = roc_curve(y, score if auc >= 0.5 else -score) + best = int(np.argmax(tpr - fpr)) + cut = cuts[best] if auc >= 0.5 else -cuts[best] + + results.append({ + "threshold_dbz": threshold, + "column": column, + "auc": round(auc, 4), + "oriented_auc": round(oriented, 4), + "direction": direction, + "youden_cut_pct": round(float(cut), 4), + "sensitivity": round(float(tpr[best]), 4), + "specificity": round(float(1 - fpr[best]), 4), + "keep_in_codebook": bool(oriented >= AUC_GATE), + }) + + return pd.DataFrame(results) + + +def main(): + features = build_feature_table() + results = evaluate(features) + + pd.set_option("display.width", 160) + print("\nCellular vs QLCS discrimination on the training split " + f"(n={len(features)}, positive class = {POSITIVE}):\n") + print(results.to_string(index=False)) + + kept = results[results["keep_in_codebook"]] + print(f"\nAUC gate = {AUC_GATE}. " + f"{len(kept)} of {len(results)} thresholds clear it.") + if kept.empty: + print("=> No quantitative coverage criteria belong in the codebook; " + "rely on the topological description.") + else: + for _, r in kept.iterrows(): + print(f"=> {r['column']}: cut at {r['youden_cut_pct']}% " + f"({r['direction']}, AUC {r['oriented_auc']})") + + results.to_csv(OUT_CSV, index=False) + print(f"\nWrote {OUT_CSV}") + + +if __name__ == "__main__": + main() diff --git a/notebooks/reports/svrimg_threshold_calibration.csv b/notebooks/reports/svrimg_threshold_calibration.csv new file mode 100644 index 0000000..cce2de0 --- /dev/null +++ b/notebooks/reports/svrimg_threshold_calibration.csv @@ -0,0 +1,7 @@ +threshold_dbz,column,auc,oriented_auc,direction,youden_cut_pct,sensitivity,specificity,keep_in_codebook +10,pct_gates_10dbz,0.9026,0.9026,higher=QLCS,29.9849,0.9145,0.7468,True +20,pct_gates_20dbz,0.9288,0.9288,higher=QLCS,18.0039,0.8706,0.8318,True +30,pct_gates_30dbz,0.9346,0.9346,higher=QLCS,6.4771,0.8947,0.8244,True +40,pct_gates_40dbz,0.889,0.889,higher=QLCS,2.3519,0.7829,0.8281,True +50,pct_gates_50dbz,0.6006,0.6006,higher=QLCS,0.3136,0.5658,0.6081,False +60,pct_gates_60dbz,0.3161,0.6839,higher=Cellular,0.0,0.739,0.597,True diff --git a/tests/test_svrimg.py b/tests/test_svrimg.py new file mode 100644 index 0000000..52d5201 --- /dev/null +++ b/tests/test_svrimg.py @@ -0,0 +1,278 @@ +import os + +import numpy as np +import pytest +from PIL import Image + +CODEBOOK = os.path.join( + os.path.dirname(__file__), "..", "example_codebooks", "CODEBOOK_SVRIMG.md" +) + +CLASS_CSV = ( + "UNID,Class Code,Class Name\n" + "199604200208z000000206,0,Cellular\n" + "200210291000z000000352,1,QLCS\n" + "201304011200z000000111,0,Cellular\n" + "201504011200z000000222,2,Tropical\n" +) + +INDEX_CSV = ( + "unid,area,convection_area,intense_area,intensity_variance,kind," + "max_intensity,mean_intensity,radar_time,report_time\n" + "199604200208z000000206,1,2,3,4,tor,63,5,4/20/1996 2:00,4/20/1996 2:08\n" + "200210291000z000000352,1,2,3,4,tor,65,5,10/29/2002 10:00,10/29/2002 10:00\n" + "201304011200z000000111,1,2,3,4,tor,55,5,4/1/2013 12:00,4/1/2013 12:00\n" + "201504011200z000000222,1,2,3,4,tor,55,5,4/1/2015 12:00,4/1/2015 12:00\n" +) + + +def _dummy_ramp_palette(): + """The meaningless linear palette SVRIMG PNGs carry.""" + return [v % 256 for i in range(256) for v in (3 * i, 3 * i + 1, 3 * i + 2)] + + +def _write_palette_png(path, array): + """Write ``array`` as a palette-mode PNG, the way SVRIMG distributes them.""" + im = Image.fromarray(np.asarray(array, dtype=np.uint8), mode="P") + im.putpalette(_dummy_ramp_palette()) + os.makedirs(os.path.dirname(path), exist_ok=True) + im.save(path) + return path + + +@pytest.fixture(autouse=True) +def _no_network(monkeypatch): + """Fail loudly if any test tries to fetch from the SVRIMG server.""" + import lars.preprocessing.svrimg as svrimg + + def _offline(url, dest): + if os.path.exists(dest): + return dest + raise AssertionError(f"test attempted a network fetch: {url}") + + monkeypatch.setattr(svrimg, "_download", _offline) + + +def _seed_cache(data_dir): + """Pre-populate the download cache so no test touches the network.""" + os.makedirs(data_dir, exist_ok=True) + with open(os.path.join(data_dir, "sample_classifications_96-17.csv"), "w") as f: + f.write(CLASS_CSV) + with open(os.path.join(data_dir, "96-17_tor_utc_svrimg_index.csv"), "w") as f: + f.write(INDEX_CSV) + return data_dir + + +# -------------------------------------------------------------------------- +# Palette decoding -- the regression guard for the greyscale-looking PNGs +# -------------------------------------------------------------------------- + + +def test_read_svrimg_png_returns_palette_indices_not_rgb(tmp_path): + from lars.preprocessing.svrimg import read_svrimg_png + + expected = np.zeros((136, 136), dtype=np.uint8) + expected[10:20, 10:20] = 55 + expected[30, 30] = 72 + path = _write_palette_png(str(tmp_path / "scene.png"), expected) + + np.testing.assert_array_equal(read_svrimg_png(path), expected) + + +def test_read_svrimg_png_rejects_non_palette_image(tmp_path): + from lars.preprocessing.svrimg import read_svrimg_png + + path = str(tmp_path / "rgb.png") + Image.fromarray(np.zeros((8, 8, 3), dtype=np.uint8), mode="RGB").save(path) + + with pytest.raises(ValueError, match="palette-mode"): + read_svrimg_png(path) + + +# -------------------------------------------------------------------------- +# Coverage statistics +# -------------------------------------------------------------------------- + + +def test_reflectivity_stats_against_hand_computed_array(): + from lars.preprocessing.svrimg import reflectivity_stats + + arr = np.zeros((10, 10), dtype=np.uint8) # 100 gates + arr[0, :10] = 35 # 10 gates > 30 + arr[1, :5] = 55 # 5 gates > 50 + arr[2, :20] = 15 # 10 gates > 10 + + stats = reflectivity_stats(arr, (10, 30, 50)) + + assert stats["n_gates_10dbz"] == 25 # 10 + 5 + 10 + assert stats["n_gates_30dbz"] == 15 # 10 + 5 + assert stats["n_gates_50dbz"] == 5 + assert stats["pct_gates_30dbz"] == 15.0 + assert stats["ref_min"] == 15.0 # echoing gates only, zeros excluded + assert stats["ref_max"] == 55.0 + + +def test_reflectivity_stats_handles_echo_free_scene(): + from lars.preprocessing.svrimg import reflectivity_stats + + stats = reflectivity_stats(np.zeros((5, 5), dtype=np.uint8), (10,)) + + assert np.isnan(stats["ref_min"]) and np.isnan(stats["ref_max"]) + assert stats["n_gates_10dbz"] == 0 + assert stats["pct_gates_10dbz"] == 0.0 + + +def test_split_for_year_matches_published_split(): + from lars.preprocessing.svrimg import split_for_year + + assert split_for_year(1996) == "train" + assert split_for_year(2011) == "train" + assert split_for_year(2012) == "val" + assert split_for_year(2013) == "val" + assert split_for_year(2014) == "test" + assert split_for_year(2017) == "test" + assert split_for_year(2018) is None + + +# -------------------------------------------------------------------------- +# Label fetching +# -------------------------------------------------------------------------- + + +def test_fetch_svrimg_labels_collapses_to_cellular_and_qlcs(tmp_path): + from lars.preprocessing.svrimg import fetch_svrimg_labels + + labels = fetch_svrimg_labels(_seed_cache(str(tmp_path))) + + assert set(labels["label"]) == {"Cellular", "QLCS"} + assert "201504011200z000000222" not in labels.index # Tropical dropped + assert len(labels) == 3 + + +def test_fetch_svrimg_labels_assigns_splits_and_parses_times(tmp_path): + from lars.preprocessing.svrimg import fetch_svrimg_labels + + labels = fetch_svrimg_labels(_seed_cache(str(tmp_path))) + + assert labels.loc["199604200208z000000206", "split"] == "train" + assert labels.loc["201304011200z000000111", "split"] == "val" + assert labels.loc["200210291000z000000352", "radar_time"].hour == 10 + assert labels.loc["199604200208z000000206", "max_intensity"] == 63 + + +def test_fetch_svrimg_labels_can_keep_all_classes(tmp_path): + from lars.preprocessing.svrimg import fetch_svrimg_labels + + labels = fetch_svrimg_labels(_seed_cache(str(tmp_path)), classes=None) + + assert "Tropical" in set(labels["label"]) + assert len(labels) == 4 + + +# -------------------------------------------------------------------------- +# End-to-end preprocessing +# -------------------------------------------------------------------------- + + +def test_preprocess_svrimg_data_matches_radar_schema(tmp_path): + from lars.preprocessing.svrimg import preprocess_svrimg_data + + data_dir = _seed_cache(str(tmp_path / "cache")) + scene = np.zeros((136, 136), dtype=np.uint8) + scene[:40, :40] = 45 # 1600 gates > 30 of 18496 + _write_palette_png( + os.path.join(data_dir, "tor", "1996", "04", + "199604200208z000000206.png"), + scene, + ) + _write_palette_png( + os.path.join(data_dir, "tor", "2002", "10", + "200210291000z000000352.png"), + np.zeros((136, 136), dtype=np.uint8), + ) + + df = preprocess_svrimg_data( + output_path=str(tmp_path / "out"), + data_dir=data_dir, + splits=("train",), + cmap="viridis", # avoid the cmweather dependency + ) + + assert len(df) == 2 # both 1996 and 2002 are train + row = df[df["unid"] == "199604200208z000000206"].iloc[0] + assert row["unid"] == "199604200208z000000206" + assert row["label"] == "Cellular" + assert row["split"] == "train" + assert row["n_gates_30dbz"] == 1600 + assert row["pct_gates_30dbz"] == round(1600 / 18496 * 100, 4) + assert os.path.exists(row["file_path"]) + assert df.index.name == "time" + + # Same schema as preprocess_radar_data, plus unid/split. + base = {"file_path", "label", "ref_min", "ref_max"} + base |= {f"n_gates_{t}dbz" for t in (10, 20, 30, 40, 50)} + base |= {f"pct_gates_{t}dbz" for t in (10, 20, 30, 40, 50)} + assert base <= set(df.columns) + assert {"unid", "split"} <= set(df.columns) + + +def test_preprocess_svrimg_data_rejects_axis_kwarg(tmp_path): + from lars.preprocessing.svrimg import preprocess_svrimg_data + + with pytest.raises(ValueError, match="Do not pass in an axis"): + preprocess_svrimg_data( + output_path=str(tmp_path / "out"), + data_dir=_seed_cache(str(tmp_path / "cache")), + ax="something", + ) + + +# -------------------------------------------------------------------------- +# Codebook wiring +# -------------------------------------------------------------------------- + + +def test_codebook_exposes_binary_convective_mode_classes(): + from lars.nepho.inference import categories_from_codebook + + assert set(categories_from_codebook(CODEBOOK)) == { + "Cellular", "QLCS", "Ambiguous / Uncertain", + } + + +def test_codebook_colormap_uses_zero_floor(): + from lars.nepho.inference import colormap_from_codebook + + # SVRIMG reflectivity is floored at 0, unlike the -20 used for CSAPR2. + assert colormap_from_codebook(CODEBOOK) == { + "colormap": "NWSRef", "vmin": 0, "vmax": 80, + } + + +def test_codebook_criterion_targets_the_emitted_column(): + from lars.nepho.inference import criteria_from_codebook + + criteria = criteria_from_codebook(CODEBOOK) + + assert set(criteria) == {"Cellular"} + criterion = criteria["Cellular"][0] + assert criterion["field"] == "pct_gates_30dbz" + assert criterion["max_value"] == 6.5 + assert criterion["reclassify_as"] == "QLCS" + + +def test_codebook_guidelines_are_non_empty(): + from lars.nepho.inference import guidelines_from_codebook + + assert len(guidelines_from_codebook(CODEBOOK)) >= 3 + + +def test_standardize_labels_canonicalises_svrimg_classes(): + import pandas as pd + from lars.preprocessing.labels import standardize_labels + + df = pd.DataFrame({"label": ["cellular", "QLCS", "qlcs", "Cellular"]}) + + assert list(standardize_labels(df)["label"]) == [ + "Cellular", "QLCS", "QLCS", "Cellular", + ]