diff --git a/era5_analysis/README.md b/era5_analysis/README.md new file mode 100644 index 0000000..80e244b --- /dev/null +++ b/era5_analysis/README.md @@ -0,0 +1,89 @@ +# ERA5 buoyancy-topology analysis + +Population study of lifted-parcel buoyancy topologies and LNB-convention +sensitivity over 569,478 random ERA5 ocean columns (1980–2024), for the three +parcels used by the tcpyPI potential-intensity calculation. Companion to +`discontinuity_analysis/` (the issue #77 root-cause study): this directory +measures, at climatological scale, how often the multi-crossing buoyancy +structures behind that failure occur and how much the LNB definition matters. + +## Data + +Input sample (public ARCO-ERA5, no credentials needed to rebuild): + + s3://gridded-data-dev/bmares/2026-07-18-era5-tcpypi-sample/ + tcpypi_inputs_1000h_seed20260718.parquet (~40 MB, ~142k columns) + tcpypi_inputs_4000h_seed20260719.parquet (~154 MB, 569,478 columns) + +Rebuild from scratch with `build_tcpypi_era5_sample.py` (see its docstring for +the exact commands and the conversion recipe; schema metadata inside the +parquet carries units, levels, sampling parameters, and caveats). + +## Pipeline + +All scripts read/write a work directory: `export ERA5_SCRATCH=...` (default: +`./work`). They import tcpyPI from this repository's `src/` (numba required; +scientific results in the writeups used the max-work cape() of the +`lnb-max-work` branch — the scan scripts implement both conventions +internally, so the checked-out kernel only matters for `validate_scan.py` +and `pi_conv.py` cross-checks). + +1. `convert.py` — parquet -> `profiles_converted.npz` (units, bottom-up + level order, r = q/(1-q)); expects the parquet at + `$ERA5_SCRATCH/profiles.parquet`. +2. `buoyancy_scan.py` — environmental-parcel topology/crossings/partial-sums + scan under four LNB conventions -> `buoyancy_scan_results.npz`. +3. `validate_scan.py` — cross-check E_top/E_max against real legacy/max-work + `cape()` on random profiles. +4. `export_stats.py` — per-profile results parquet + population statistics. +5. `pi_conv.py` — full pi() under four LNB conventions (2.28M solves) + -> `pi_conv_results.npz` (+ flat parquet); includes the wild-population + legacy non-convergence census. +6. `scan_final.py` — the definitive three-parcel scan (A environmental, + B eyewall, C saturated core at each column's converged max-work P_M), + perturbation-convention topology labels, candidate-based conventions + -> `parcel_topology_final.npz`. +7. `scan_pi_parcels.py` — earlier two-parcel scan + hybrid-vs-ln(p) + quadrature comparison (kept for the quadrature numbers). +8. `gen_curves.py` + `plot_curves.py` — spaghetti figures per topology class + (`figures/curves_parcel_{A,B,C}.png`). +9. `prep_artifact_data.py` + `sounding_scope_template.html` — the interactive + "Sounding Scope" flipbook (inject `artifact_data.json` into the template's + `__DATA__` placeholder to obtain a single self-contained HTML file). + +## Headline results (TC-relevant subset: SST >= 26 C, sp >= 1000 hPa) + +- Topology frequencies (perturbation convention: leading sign taken + infinitesimally above launch): parcel A is a rich mixture (`+-+-` 32%, + `-+-` 17%, `+-+-+-` 16%, `+-` 14%, ...); parcel B is 84% `+-` with a 16% + multi-hump tail; parcel C is 98.4% `+-` plus 1.24% `+` (buoyant at the + 70 hPa retained-column top). +- LNB-convention sensitivity is inversely proportional to parcel energy: + E_top clamps to 0 while E_max > 0 on 15.4% (A), 1.75% (B), 0.00% (C) of + columns; |E_top - E_max| > 1 J/kg on 13.1% / 5.7% / 0.0%. +- Full-pi() convention comparison over 476k columns (SST > 5 C): + legacy topmost-level LNB fails to converge on 8 columns (~1 in 60,000; + median max-work VMAX 64 m/s — real storm environments returned missing); + max-work fails on 0; first-crossing and lifted-ballistic conventions fail + ~20x more often than legacy. Where PI matters (SST >= 26), top vs max + agree to p99 = 0.53 m/s. +- Column ceiling: the topology chain (`scan_final.py` onward) retains + levels down to 30 hPa (ptop=20) -- the minimal ceiling observing every + buoyancy crossing in the sample. Highest crossing: 47.6 hPa (saturated + core parcel, interpolated between the 50 and 30 hPa nodes; parcel B max + 65.0, parcel A max 87.5). At this ceiling zero profiles clip AND zero + entropy solves fail (retaining the 30-20 hPa layer caused 21 failures). + tcpyPI's own ptop=50 convention (used by the pi()/CAPE scans above) + truncates at 70 hPa and clips 1.24% of TC-relevant core parcels, whose + true crossings sit at median 68.6 hPa. The solver's `ENEW > P-1` guard + prevents extension beyond ~10-20 hPa without modification. Buoyancy-plot + x-ranges are set by the deepest interior dip (the most negative value a + curve reaches before its last positive level) plus 10% — the terminal + stratospheric plunge never sets the range. +- Strict from-rest ballistic CAPE is degenerate for surface-launched + environmental parcels (b(launch) = 0 by construction and surface CIN or + neutral layering is universal); the lifted-to-LFC variant is the + meaningful ballistic convention. + +Outputs referenced above (results parquets, npz files) are regenerated by +the pipeline; the two published to date live next to the input sample on S3. diff --git a/era5_analysis/build_tcpypi_era5_sample.py b/era5_analysis/build_tcpypi_era5_sample.py new file mode 100644 index 0000000..f13ce0e --- /dev/null +++ b/era5_analysis/build_tcpypi_era5_sample.py @@ -0,0 +1,225 @@ +#!/usr/bin/env python3 +"""Assemble generic tcpyPI input samples from ARCO-ERA5. + +What it does +------------ +For N random hours in [1980, 2025), it reads four ERA5 fields from the public +Analysis-Ready Cloud-Optimized ERA5 store on GCS (no credentials needed): + + sea_surface_temperature, surface_pressure (single level) + temperature, specific_humidity (37 pressure levels) + +For each hour it draws M cosine-weighted (equal-area) random lat/lon points, +snaps them to the 0.25 deg grid, discards land points (NaN SST), and emits one +record per surviving ocean point containing latitude, longitude, time, the two +surface scalars, and the two 37-level vertical profiles. Output is a single +zstd-compressed parquet file whose schema metadata carries the pressure levels, +units, the exact tcpyPI conversion recipe, the sampling parameters, and known +data caveats. + +Provenance / vendoring +---------------------- +The ARCO-ERA5 URL and the variable set are taken from Climate Central's +attribution pipeline (attribution/hurricanes/aggregate.py, which uses the same +store for hurricane potential-intensity climatology) but are inlined here so +this script stands entirely on its own. + +Requirements (all pip-installable, no attribution package): + pip install numpy pandas pyarrow xarray zarr gcsfs + +Reading the ARCO store is anonymous/public; writing a local parquet needs no +credentials. (Uploading the result to S3 is out of scope for this script.) + +Reproduce the two published samples +----------------------------------- + python build_tcpypi_era5_sample.py --n-hours 1000 --seed 20260718 \ + --threads 8 --out tcpypi_inputs_1000h_seed20260718.parquet + python build_tcpypi_era5_sample.py --n-hours 4000 --seed 20260719 \ + --threads 8 --out tcpypi_inputs_4000h_seed20260719.parquet + +Convert a stored record to tcpyPI inputs +---------------------------------------- + SST_C = sst_K - 273.15 + MSL_hPa = sp_Pa / 100.0 # over ocean surface_pressure ~= MSL + P_hPa = pressure_level_hPa # from file metadata (ascending index) + T_C = temperature_K - 273.15 + R_gkg = 1000 * q / (1 - q) # q = specific_humidity_kgkg +""" + +import argparse +import json +import time as _time +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +import pandas as pd +import pyarrow as pa +import pyarrow.parquet as pq +import xarray as xr + +# --- vendored constants (see attribution/hurricanes/aggregate.py) ------------- +ARCO_URL = "gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3" +SURFACE_VARS = ("sea_surface_temperature", "surface_pressure") +PROFILE_VARS = ("temperature", "specific_humidity") + +# ERA5 0.25 deg grid geometry: latitude[i] = 90 - 0.25*i (i in 0..720), +# longitude[j] = 0.25*j (j in 0..1439). +LAT0, DLAT, NLAT = 90.0, 0.25, 721 +DLON, NLON = 0.25, 1440 +TIME_START, TIME_STOP = "1980-01-01", "2025-01-01" # [start, stop) + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + p.add_argument("--n-hours", type=int, default=1000, + help="number of distinct random hours to sample") + p.add_argument("--n-points", type=int, default=200, + help="random lat/lon points drawn per hour (before land filter)") + p.add_argument("--seed", type=int, default=20260718) + p.add_argument("--threads", type=int, default=8, + help="concurrent hours (I/O bound; each holds ~0.3 GB)") + p.add_argument("--out", type=str, required=True, help="output .parquet path") + return p.parse_args() + + +def open_arco(): + """Open the public ARCO-ERA5 store lazily (anonymous access).""" + return xr.open_zarr(ARCO_URL, chunks=None, storage_options=dict(token="anon")) + + +def sample_grid_indices(rng, n): + """Draw n cosine-weighted (equal-area) points; return nearest grid indices. + + Latitude density proportional to cos(lat) is achieved by lat = arcsin(U), + U ~ Uniform(-1, 1); longitude is uniform on [0, 360). + """ + lat = np.degrees(np.arcsin(rng.uniform(-1.0, 1.0, size=n))) + lon = rng.uniform(0.0, 360.0, size=n) + lat_idx = np.rint((LAT0 - lat) / DLAT).astype(int).clip(0, NLAT - 1) + lon_idx = np.rint(lon / DLON).astype(int) % NLON + return lat_idx, lon_idx + + +def process_hour(ds, it, timestamp, rng, n_points): + """Read one hour's global fields, sample points, keep ocean points.""" + t_arr = ds.temperature.isel(time=it).values # (37, 721, 1440) K + q_arr = ds.specific_humidity.isel(time=it).values # (37, 721, 1440) kg/kg + sp_arr = ds.surface_pressure.isel(time=it).values # (721, 1440) Pa + sst_arr = ds.sea_surface_temperature.isel(time=it).values # (721, 1440) K + + lat_idx, lon_idx = sample_grid_indices(rng, n_points) + sst_pt = sst_arr[lat_idx, lon_idx] + keep = np.isfinite(sst_pt) + lat_idx, lon_idx = lat_idx[keep], lon_idx[keep] + if lat_idx.size == 0: + return None + + return dict( + latitude=(LAT0 - DLAT * lat_idx).astype("float32"), + longitude=(((DLON * lon_idx + 180.0) % 360.0) - 180.0).astype("float32"), + time=np.repeat(np.datetime64(timestamp, "ns"), lat_idx.size), + sst_K=sst_pt[keep].astype("float32"), + sp_Pa=sp_arr[lat_idx, lon_idx].astype("float32"), + temperature_K=t_arr[:, lat_idx, lon_idx].T.astype("float32"), # (n, 37) + specific_humidity_kgkg=q_arr[:, lat_idx, lon_idx].T.astype("float32"), + ) + + +def build_metadata(levels, args): + return { + b"source": ARCO_URL.encode(), + b"pressure_level_hPa": json.dumps(levels.tolist()).encode(), + b"units": json.dumps({ + "sst_K": "K", "sp_Pa": "Pa", "temperature_K": "K", + "specific_humidity_kgkg": "kg/kg", "pressure_level": "hPa", + "latitude": "degrees_north", "longitude": "degrees_east (-180..180)", + }).encode(), + b"profile_order": ( + b"ascending pressure_level index " + b"(level[0]=1 hPa ... level[-1]=1000 hPa)" + ), + b"tcpyPI_conversion": json.dumps({ + "SST_C": "sst_K - 273.15", + "MSL_hPa": "sp_Pa / 100.0 (over ocean surface_pressure ~= MSL)", + "P_hPa": "pressure_level_hPa", + "T_C": "temperature_K - 273.15", + "R_gkg": "1000 * q/(1-q) with q = specific_humidity_kgkg", + }).encode(), + b"sampling": json.dumps({ + "n_hours": args.n_hours, "n_points_per_hour": args.n_points, + "seed": args.seed, "time_window": f"[{TIME_START}, {TIME_STOP})", + "lat_weighting": "cosine (equal-area via arcsin)", + "land_filter": "dropped points with NaN sea_surface_temperature", + }).encode(), + b"known_caveats": json.dumps({ + "inland_lakes": ( + "ERA5 sea_surface_temperature is defined over large/high-elevation " + "inland lakes (Victoria ~885 hPa, Titicaca ~633 hPa, Urmia, etc.), so " + "~0.1% of records are lakes, not open ocean. Filter on sp_Pa " + "(e.g. >= 95000) if only open ocean is desired." + ), + "negative_specific_humidity": ( + "A tiny fraction of stratospheric q values are slightly negative " + "(~-4e-6 kg/kg), an ERA5 interpolation artifact; negligible for tcpyPI." + ), + }).encode(), + } + + +def main(): + args = parse_args() + rng = np.random.default_rng(args.seed) + + ds = open_arco() + levels = ds.level.values.astype("int64") # hPa + + times = pd.DatetimeIndex(ds.time.values) + valid = np.where((times >= TIME_START) & (times < TIME_STOP))[0] + chosen = np.sort(rng.choice(valid, size=args.n_hours, replace=False)) + print(f"sampling {args.n_hours} hours x {args.n_points} pts; " + f"valid hour pool={valid.size}; threads={args.threads}", flush=True) + + # Independent per-hour RNGs => result is identical regardless of thread order. + hour_rngs = [np.random.default_rng([args.seed, int(it)]) for it in chosen] + + def _run(k): + it = int(chosen[k]) + t0 = _time.time() + rec = process_hour(ds, it, times[it], hour_rngs[k], args.n_points) + kept = 0 if rec is None else rec["latitude"].size + print(f" hour {k + 1}/{args.n_hours} {times[it]} " + f"kept={kept}/{args.n_points} ({_time.time() - t0:.1f}s)", flush=True) + return rec + + t_start = _time.time() + if args.threads > 1: + with ThreadPoolExecutor(max_workers=args.threads) as ex: + recs = list(ex.map(_run, range(args.n_hours))) + else: + recs = [_run(k) for k in range(args.n_hours)] + recs = [r for r in recs if r is not None] + + cat = {key: np.concatenate([r[key] for r in recs]) for key in recs[0]} + elapsed = _time.time() - t_start + print(f"total ocean records: {cat['latitude'].size} ({elapsed:.1f}s, " + f"{elapsed / args.n_hours:.2f}s/hour)", flush=True) + + table = pa.table({ + "latitude": pa.array(cat["latitude"]), + "longitude": pa.array(cat["longitude"]), + "time": pa.array(cat["time"]), + "sst_K": pa.array(cat["sst_K"]), + "sp_Pa": pa.array(cat["sp_Pa"]), + "temperature_K": pa.array( + list(cat["temperature_K"]), type=pa.list_(pa.float32()) + ), + "specific_humidity_kgkg": pa.array( + list(cat["specific_humidity_kgkg"]), type=pa.list_(pa.float32()) + ), + }).replace_schema_metadata(build_metadata(levels, args)) + pq.write_table(table, args.out, compression="zstd") + print(f"wrote {args.out} rows={table.num_rows}", flush=True) + + +if __name__ == "__main__": + main() diff --git a/era5_analysis/buoyancy_scan.py b/era5_analysis/buoyancy_scan.py new file mode 100644 index 0000000..a36b466 --- /dev/null +++ b/era5_analysis/buoyancy_scan.py @@ -0,0 +1,271 @@ +"""Buoyancy-topology scan over ERA5 profiles (tcpyPI thermodynamic conventions). + +For each profile, lift the lowest-level parcel (cape()'s environmental-parcel +semantics, ascent_flag=0, ptop=50 -> 28 retained levels 1000..70 hPa), compute +the buoyancy (density-temperature difference) at each level, and extract: +sign topology, interpolated crossings (log-pressure), per-region unsigned CAPE +contributions, running signed partial sums, argmax, and CAPE under the four +LNB conventions (legacy topmost / max-work / ballistic-reachable / first-LNB). + +Outputs padded arrays in buoyancy_scan_results.npz (MAXR regions). +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import os +import sys +import time + +sys.path.insert(0, _SRC) + +import numpy as np +from numba import njit, prange + +from tcpyPI import constants, utilities +from tcpyPI.pi import solve_temperature_from_entropy + + +MAXR = 28 # = retained level count; regions can never exceed this +PTOP = 50.0 + +RD = constants.RD + + +@njit(cache=True) +def scan_one(TC, Rgkg, P, nlvl, srow, + region_abs, region_peak, psums, cross_lnp): + """Scan one profile, filling srow = [iflag, n, first_pos, argmax, clipped, + clamped, E_top, E_max, E_reach, E_first, LNB_top, LNB_max, LNB_reach, + LNB_first, lcl, b_low, b_top] and the padded per-region rows.""" + T = TC + 273.15 + R = Rgkg * 0.001 + + TP = T[0] + RP = R[0] + PP = P[0] + if (RP < 1e-6) or (TP < 200.0): + srow[0] = 0.0 + srow[3] = -1.0 + return + + TPC = utilities.T_ktoC(TP) + ESP = utilities.es_cc(TPC) + EVP = utilities.ev(RP, PP) + RH = min(EVP / ESP, 1.0) + S = utilities.entropy_S(TP, RP, PP) + PLCL = utilities.e_pLCL(TP, RH, PP) + + # buoyancy profile (identical formulas to cape()) + b = np.empty(nlvl) + for j in range(nlvl): + if P[j] >= PLCL: + TG = TP * (P[j] / PP) ** (constants.RD / constants.CPD) + RG = RP + TLVR = utilities.Trho(TG, RG, RG) + TVENV = utilities.Trho(T[j], R[j], R[j]) + b[j] = TLVR - TVENV + else: + TG, RG, IFLAG = solve_temperature_from_entropy(S, P[j], RP, T[j]) + if IFLAG == 2: + srow[0] = 2.0 + srow[3] = -1.0 + srow[14] = PLCL + return + TLVR = utilities.Trho(TG, RP, RG) # reversible ascent (ascent_flag=0) + TENV = utilities.Trho(T[j], R[j], R[j]) + b[j] = TLVR - TENV + + # --- segmentation into alternating sign regions (b > 0 is '+') --- + first_pos = 1 if b[0] > 0.0 else 0 + k = 0 # current region index + W = 0.0 # running signed work + region_abs[0] = 0.0 + region_peak[0] = abs(b[0]) + cross_p = np.empty(MAXR) # crossing pressures (region ends), local + for j in range(1, nlvl): + b0 = b[j - 1] + b1 = b[j] + pos0 = b0 > 0.0 + pos1 = b1 > 0.0 + if pos0 != pos1: + # region boundary: interpolated crossing (node itself if b hits 0) + if b0 * b1 < 0.0: + PC = (P[j] * b0 - P[j - 1] * b1) / (b0 - b1) + else: + PC = P[j] if b1 == 0.0 else P[j - 1] + piece0 = RD * b0 * (P[j - 1] - PC) / (P[j - 1] + PC) + piece1 = RD * b1 * (PC - P[j]) / (PC + P[j]) + region_abs[k] += abs(piece0) + cross_p[k] = PC + cross_lnp[k] = np.log(PC) + W += piece0 + psums[k] = W + k += 1 + region_abs[k] = abs(piece1) + region_peak[k] = abs(b1) + W += piece1 + else: + piece = RD * (b0 + b1) * (P[j - 1] - P[j]) / (P[j] + P[j - 1]) + region_abs[k] += abs(piece) + if abs(b1) > region_peak[k]: + region_peak[k] = abs(b1) + W += piece + n = k + 1 + psums[k] = W + cross_p[k] = P[nlvl - 1] # last region "ends" at the profile top + nstore = n + + # --- diagnostics from the partial sums --- + # argmax of partial sums (-1 if all <= 0) + best = -np.inf + aidx = -1 + for kk in range(nstore): + if psums[kk] > best: + best = psums[kk] + aidx = kk + if best <= 0.0: + aidx = -1 + + # sign of region kk: first_pos alternating + top_sign_pos = ((n - 1) % 2 == 0) == (first_pos == 1) + clipped = 1 if top_sign_pos else 0 + + # E_max (max-work) and its LNB + E_max = best if best > 0.0 else 0.0 + LNB_max = cross_p[aidx] if aidx >= 0 else 0.0 + + # E_top (legacy pcmin: signed integral to the topmost positive level, + # clamped at zero) + clamped = 0 + if n == 1 and first_pos == 0: + E_top = 0.0 + LNB_top = 0.0 + else: + if top_sign_pos: + raw = psums[nstore - 1] + LNB_top = cross_p[nstore - 1] + else: + raw = psums[nstore - 2] + LNB_top = cross_p[nstore - 2] + if raw < 0.0: + clamped = 1 + E_top = 0.0 + else: + E_top = raw + + # E_reach_strict (ballistic from rest at launch: max W over the first + # connected component of {W>0}; W monotone within regions, so region-end + # checks suffice). NOTE: b[0]=0 by construction and surface CIN is nearly + # universal, so this is ~always 0 for launch-level parcels. + E_reach = 0.0 + LNB_reach = 0.0 + for kk in range(nstore): + if psums[kk] > E_reach: + E_reach = psums[kk] + LNB_reach = cross_p[kk] + if psums[kk] <= 0.0: + break + + # index of the first positive region (the parcel's LFC-to-LNB layer) + r1 = 0 if first_pos == 1 else 1 + + # E_first (first-LNB convention: signed integral through any leading CIN + # up to the FIRST down-crossing), clamped at zero + if r1 < nstore: + E_first = psums[r1] if psums[r1] > 0.0 else 0.0 + LNB_first = cross_p[r1] + else: + E_first = 0.0 + LNB_first = 0.0 + + # E_reach_lfc (lifted-parcel ballistic: externally lifted to the LFC with + # zero kinetic energy, then coasting; kinetic energy at p is W(p) - floor, + # floor = W at LFC entry; stalls where W returns to the floor) + E_reach_lfc = 0.0 + LNB_reach_lfc = 0.0 + if r1 < nstore: + floor = psums[r1 - 1] if r1 >= 1 else 0.0 + best = -np.inf + for kk in range(r1, nstore): + if psums[kk] > best: + best = psums[kk] + LNB_reach_lfc = cross_p[kk] + if psums[kk] <= floor and kk > r1: + break + if best > 0.0: + E_reach_lfc = best + else: + LNB_reach_lfc = 0.0 + + srow[0] = 1.0 + srow[1] = n + srow[2] = first_pos + srow[3] = aidx + srow[4] = clipped + srow[5] = clamped + srow[6] = E_top + srow[7] = E_max + srow[8] = E_reach + srow[9] = E_first + srow[10] = LNB_top + srow[11] = LNB_max + srow[12] = LNB_reach + srow[13] = LNB_first + srow[14] = PLCL + srow[15] = b[0] + srow[16] = b[nlvl - 1] + srow[17] = E_reach_lfc + srow[18] = LNB_reach_lfc + return + + +@njit(parallel=True, cache=True) +def scan_all(TCs, Rs, P, nlvl, out_scalar, region_abs, region_peak, psums, cross_lnp): + nprof = TCs.shape[0] + for i in prange(nprof): + scan_one(TCs[i], Rs[i], P, nlvl, out_scalar[i], + region_abs[i], region_peak[i], psums[i], cross_lnp[i]) + + +def main(): + d = np.load(f"{SCRATCH}/profiles_converted.npz") + P_full = d["P"] + nlvl = int((P_full > PTOP).sum()) + P = P_full[:nlvl].copy() + TCs = d["TC"][:, :nlvl].copy() + Rs = d["R"][:, :nlvl].copy() + nprof = TCs.shape[0] + print(f"{nprof} profiles, {nlvl} retained levels ({P[0]:.0f}..{P[-1]:.0f} hPa)") + + out_scalar = np.full((nprof, 19), np.nan) + region_abs = np.full((nprof, MAXR), np.nan) + region_peak = np.full((nprof, MAXR), np.nan) + psums = np.full((nprof, MAXR), np.nan) + cross_lnp = np.full((nprof, MAXR - 1), np.nan) + + t0 = time.time() + scan_all(TCs[:100], Rs[:100], P, nlvl, out_scalar[:100], + region_abs[:100], region_peak[:100], psums[:100], cross_lnp[:100]) + print(f"compile+warmup: {time.time()-t0:.1f}s") + t0 = time.time() + scan_all(TCs, Rs, P, nlvl, out_scalar, region_abs, region_peak, psums, cross_lnp) + print(f"full scan: {time.time()-t0:.1f}s") + + np.savez_compressed( + f"{SCRATCH}/buoyancy_scan_results.npz", + scalar=out_scalar, region_abs=region_abs, region_peak=region_peak, + psums=psums, cross_lnp=cross_lnp, + scalar_names=np.array( + ["iflag", "n", "first_pos", "argmax", "clipped", "clamped", + "E_top", "E_max", "E_reach", "E_first", + "LNB_top", "LNB_max", "LNB_reach", "LNB_first", + "lcl_hPa", "b_low", "b_top", "E_reach_lfc", "LNB_reach_lfc"]), + ) + print("saved buoyancy_scan_results.npz") + + +if __name__ == "__main__": + main() diff --git a/era5_analysis/convert.py b/era5_analysis/convert.py new file mode 100644 index 0000000..f0c2915 --- /dev/null +++ b/era5_analysis/convert.py @@ -0,0 +1,52 @@ +"""parquet -> npz: units and ordering for the tcpyPI-convention scanner. + +Input arrays are top-down on the 37 standard ERA5 pressure levels. +Output: bottom-up (decreasing pressure) float64 arrays, T in Celsius, +mixing ratio r = q/(1-q) in g/kg, MSL-ish surface pressure in hPa. +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import numpy as np +import pandas as pd + + + +# standard ERA5 37 pressure levels, top-down (hPa) +PLEV_TD = np.array( + [1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200, 225, 250, + 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 775, 800, 825, 850, + 875, 900, 925, 950, 975, 1000], dtype=np.float64) + +df = pd.read_parquet(f"{SCRATCH}/profiles.parquet") +n = len(df) +print("profiles:", n) + +t = np.stack(df["temperature_K"].to_numpy()).astype(np.float64) # (n, 37) top-down +q = np.stack(df["specific_humidity_kgkg"].to_numpy()).astype(np.float64) + +# sanity: top-down ordering means q tiny at index 0 +assert np.nanmedian(q[:, 0]) < 1e-4 < np.nanmedian(q[:, -1]), "ordering assumption violated" + +TC = t[:, ::-1] - 273.15 # bottom-up, Celsius +r = q[:, ::-1] +r = r / (1.0 - r) * 1000.0 # specific humidity -> mixing ratio, g/kg +P = PLEV_TD[::-1].copy() # bottom-up: 1000 ... 1 + +np.savez( + f"{SCRATCH}/profiles_converted.npz", + P=P, + TC=np.ascontiguousarray(TC), + R=np.ascontiguousarray(r), + sst_C=df["sst_K"].to_numpy().astype(np.float64) - 273.15, + sp_hPa=df["sp_Pa"].to_numpy().astype(np.float64) / 100.0, + lat=df["latitude"].to_numpy().astype(np.float64), + lon=df["longitude"].to_numpy().astype(np.float64), + time=df["time"].to_numpy().astype("datetime64[s]").astype(np.int64), +) +print("saved profiles_converted.npz") +print("NaN rows (any T NaN):", int(np.isnan(TC).any(axis=1).sum())) +print("sp_hPa range:", df["sp_Pa"].min() / 100, df["sp_Pa"].max() / 100) diff --git a/era5_analysis/export_stats.py b/era5_analysis/export_stats.py new file mode 100644 index 0000000..4302cd2 --- /dev/null +++ b/era5_analysis/export_stats.py @@ -0,0 +1,124 @@ +"""Build the per-profile results parquet and print population statistics.""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import numpy as np +import pandas as pd + + + +d = np.load(f"{SCRATCH}/profiles_converted.npz") +res = np.load(f"{SCRATCH}/buoyancy_scan_results.npz") +sc = res["scalar"] +names = [str(x) for x in res["scalar_names"]] +ix = {n: i for i, n in enumerate(names)} +N = sc.shape[0] + +iflag = sc[:, ix["iflag"]].astype(np.int64) +n = np.where(np.isnan(sc[:, ix["n"]]), 0, sc[:, ix["n"]]).astype(np.int64) +first_pos = sc[:, ix["first_pos"]].astype(np.int64) + +# topology string: strictly alternating -> determined by (first sign, n) +topo = np.empty(N, dtype=object) +for i in range(N): + if iflag[i] != 1: + topo[i] = "" + continue + c0 = "+" if first_pos[i] == 1 else "-" + s = ("+-" if c0 == "+" else "-+") * ((n[i] + 1) // 2) + topo[i] = s[: n[i]] + +# ragged list columns, truncated to each row's n +cross = res["cross_lnp"] +rabs = res["region_abs"] +ps = res["psums"] +crossings, regions, partials = [], [], [] +for i in range(N): + k = n[i] + if iflag[i] != 1 or k == 0: + crossings.append([]); regions.append([]); partials.append([]) + else: + crossings.append(cross[i, : k - 1].tolist()) + regions.append(rabs[i, : k].tolist()) + partials.append(ps[i, : k].tolist()) + +df = pd.DataFrame({ + "latitude": d["lat"], "longitude": d["lon"], + "time": d["time"].astype("datetime64[s]"), + "sst_C": d["sst_C"], "sp_hPa": d["sp_hPa"], + "iflag": iflag, + "topology": topo, + "n_regions": n, + "crossings_lnp": crossings, + "region_cape_abs": regions, + "partial_sums": partials, + "argmax_partial": sc[:, ix["argmax"]].astype(np.int64), + "E_top": sc[:, ix["E_top"]], "E_max": sc[:, ix["E_max"]], + "E_reach_strict": sc[:, ix["E_reach"]], "E_first": sc[:, ix["E_first"]], + "E_reach_lfc": sc[:, ix["E_reach_lfc"]], + "LNB_top_hPa": sc[:, ix["LNB_top"]], "LNB_max_hPa": sc[:, ix["LNB_max"]], + "LNB_reach_strict_hPa": sc[:, ix["LNB_reach"]], + "LNB_reach_lfc_hPa": sc[:, ix["LNB_reach_lfc"]], + "LNB_first_hPa": sc[:, ix["LNB_first"]], + "clipped_top": sc[:, ix["clipped"]].astype(np.float64), + "clamped_E_top": sc[:, ix["clamped"]].astype(np.float64), + "lcl_hPa": sc[:, ix["lcl_hPa"]], + "b_lowest_K": sc[:, ix["b_low"]], "b_top_K": sc[:, ix["b_top"]], + "subsurface_launch": d["sp_hPa"] < 1000.0, +}) +df.to_parquet(f"{SCRATCH}/buoyancy_topology_4000h.parquet", index=False) +print(f"wrote buoyancy_topology_4000h.parquet ({len(df)} rows)\n") + +# ---------------- statistics ---------------- +ok = iflag == 1 +print(f"iflag: ok={ok.sum()} improper-parcel={(iflag==0).sum()} " + f"entropy-nonconv={(iflag==2).sum()}") +print(f"subsurface launch (sp<1000): {(d['sp_hPa']<1000).sum()} " + f"({100*(d['sp_hPa']<1000).mean():.1f}%)\n") + +v = df[ok] +print("=== topology distribution (n = number of sign regions) ===") +tc = v.groupby("n_regions").size() +for k, c in tc.items(): + print(f" n={k:2d}: {c:7d} ({100*c/len(v):6.2f}%)") +print(f" multi-crossing (n>=3): {(v.n_regions>=3).sum()} " + f"({100*(v.n_regions>=3).mean():.2f}%)") +print(" most common topologies:") +print(v.topology.value_counts().head(8).to_string()) + +print("\n=== convention disagreements (J/kg) ===") +for a, b in [("E_top", "E_max"), ("E_reach_lfc", "E_max"), ("E_first", "E_max"), + ("E_top", "E_reach_lfc"), ("E_reach_strict", "E_max")]: + dd = (v[a] - v[b]).abs() + print(f" |{a}-{b}|: >0.01: {(dd>0.01).sum():6d} ({100*(dd>0.01).mean():5.2f}%) " + f">10: {(dd>10).sum():6d} ({100*(dd>10).mean():5.2f}%) " + f"max={dd.max():9.2f} mean={dd.mean():7.3f}") + +pos = v[v.E_max > 0] +print(f"\nprofiles with positive max-work CAPE: {len(pos)} ({100*len(pos)/len(v):.1f}%)") +dd = (pos.E_top - pos.E_max).abs() +rel = dd / pos.E_max +print(f" among those, E_top disagrees >1%: {(rel>0.01).sum()} ({100*(rel>0.01).mean():.2f}%), " + f">50%: {(rel>0.5).sum()} ({100*(rel>0.5).mean():.2f}%)") +print(f" E_top clamped to 0 while E_max>0: {((pos.E_top==0)).sum()} " + f"({100*(pos.E_top==0).mean():.2f}%) (the issue-77 configuration)") +lnbd = (pos.LNB_top_hPa - pos.LNB_max_hPa).abs() +print(f" LNB_top vs LNB_max differ >50 hPa: {(lnbd>50).sum()} ({100*(lnbd>50).mean():.2f}%)") + +print(f"\nclamped E_top (raw signed integral < 0): {(v.clamped_E_top==1).sum()} " + f"({100*(v.clamped_E_top==1).mean():.2f}%)") +print(f"clipped at profile top (b>0 at {70:.0f} hPa): {(v.clipped_top==1).sum()} " + f"({100*(v.clipped_top==1).mean():.2f}%)") +print(f"argmax=-1 (all partial sums <=0): {(v.argmax_partial==-1).sum()} " + f"({100*(v.argmax_partial==-1).mean():.2f}%)") + +oc = v[~v.subsurface_launch] +print(f"\n=== same stats, ocean-valid launches only (sp>=1000; {len(oc)} rows) ===") +print(f" multi-crossing (n>=3): {100*(oc.n_regions>=3).mean():.2f}%") +p2 = oc[oc.E_max > 0] +print(f" E_max>0: {100*len(p2)/len(oc):.1f}%; E_top==0 among them: " + f"{100*(p2.E_top==0).mean():.2f}%; |E_top-E_max|>10 J/kg: " + f"{100*((p2.E_top-p2.E_max).abs()>10).mean():.2f}%") diff --git a/era5_analysis/figures/curves_parcel_A.png b/era5_analysis/figures/curves_parcel_A.png new file mode 100644 index 0000000..018c4d7 Binary files /dev/null and b/era5_analysis/figures/curves_parcel_A.png differ diff --git a/era5_analysis/figures/curves_parcel_B.png b/era5_analysis/figures/curves_parcel_B.png new file mode 100644 index 0000000..c316ccb Binary files /dev/null and b/era5_analysis/figures/curves_parcel_B.png differ diff --git a/era5_analysis/figures/curves_parcel_C.png b/era5_analysis/figures/curves_parcel_C.png new file mode 100644 index 0000000..97ba5dc Binary files /dev/null and b/era5_analysis/figures/curves_parcel_C.png differ diff --git a/era5_analysis/gen_curves.py b/era5_analysis/gen_curves.py new file mode 100644 index 0000000..5695c55 --- /dev/null +++ b/era5_analysis/gen_curves.py @@ -0,0 +1,77 @@ +"""Sample buoyancy curves per (parcel, topology class) for spaghetti plots.""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import sys + +import numpy as np + +sys.path.insert(0, SCRATCH) +sys.path.insert(0, _SRC) + +from scan_final import buoyancy, EPS, PTOP # njit'd +from tcpyPI import utilities + +NSAMP = 150 + +d = np.load(f"{SCRATCH}/profiles_converted.npz") +z = np.load(f"{SCRATCH}/parcel_topology_final.npz") +P_full = d["P"] +nlvl = int((P_full > PTOP).sum()) +P = P_full[:nlvl].copy() +tc = (d["sst_C"] >= 26) & (d["sp_hPa"] >= 1000) & np.isfinite(z["PM"]) + +def topo_str(fs, n): + c0 = "+" if fs > 0 else "-" + return (("+-" if c0 == "+" else "-+") * ((int(n) + 1) // 2))[: int(n)] + +rng = np.random.default_rng(3) +out = {} +CLASSES = {"A": 6, "B": 4, "C": 3} +for pname, nclass in CLASSES.items(): + o = z[pname] + ok = tc & np.isfinite(o[:, 0]) & (o[:, 0] != 0) + keys = [(fs, n) for fs, n in {(o[i, 0], o[i, 1]) for i in np.where(ok)[0]}] + counts = {k: 0 for k in keys} + lab = np.full(len(o), -1) + keymap = {} + for i in np.where(ok)[0]: + k = (o[i, 0], o[i, 1]) + if k not in keymap: + keymap[k] = len(keymap) + lab[i] = keymap[k] + counts[k] += 1 + top = sorted(counts, key=lambda k: -counts[k])[:nclass] + Ntot = ok.sum() + for k in top: + rows = np.where(ok & (lab == keymap[k]))[0] + sel = rng.choice(rows, size=min(NSAMP, len(rows)), replace=False) + curves = np.empty((len(sel), nlvl)) + b = np.empty(nlvl) + for m, i in enumerate(sel): + T = d["TC"][i, :nlvl] + 273.15 + R = d["R"][i, :nlvl] * 0.001 + if pname == "A": + f = buoyancy(T[0], R[0], P[0], T, R, P, nlvl, b) + else: + PP = min(z["PM"][i], 1000.0) + if pname == "B": + MSL = d["sp_hPa"][i] + RP = EPS * R[0] * MSL / (PP * (EPS + R[0]) - R[0] * MSL) + f = buoyancy(T[0], RP, PP, T, R, P, nlvl, b) + else: + SSTK = d["sst_C"][i] + 273.15 + ES0 = utilities.es_cc(d["sst_C"][i]) + f = buoyancy(SSTK, utilities.rv(ES0, PP), PP, T, R, P, nlvl, b) + curves[m] = b if f == 1 else np.nan + t = topo_str(*k) + share = 100.0 * counts[k] / Ntot + out[f"{pname}|{t}|{share:.2f}"] = curves + print(f"parcel {pname} {t:14s} {share:6.2f}% sampled {len(sel)}") + +np.savez_compressed(f"{SCRATCH}/topology_curves.npz", P=P, + **{k: v for k, v in out.items()}) +print("saved topology_curves.npz") diff --git a/era5_analysis/pi_conv.py b/era5_analysis/pi_conv.py new file mode 100644 index 0000000..af86647 --- /dev/null +++ b/era5_analysis/pi_conv.py @@ -0,0 +1,257 @@ +"""Potential intensity under four LNB conventions, over all ERA5 profiles. + +Conventions (terminal-point rule for every CAPE evaluation inside pi()): + 0 = top legacy pcmin: signed W at the crossing above the topmost + positive grid level (or profile top), clamped at 0 + 1 = max max-W: maximum of W over all region-end candidates (>=0) + 2 = first signed W at the first down-crossing, clamped at 0 + 3 = reach_lfc lifted-ballistic: best W in the component reachable from the + first positive region with floor = W at its entry + +Quadrature: running signed work with trapezoids split at every crossing +(<=1% from cape()'s unsplit-interval arithmetic; validated separately). +Outputs pi_conv_results.npz: (nprof, 4) arrays VMAX, PMIN, IFL, TO, OTL, NIT. +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import sys +import time + +sys.path.insert(0, _SRC) + +import numpy as np +from numba import njit, prange + +from tcpyPI import constants, utilities +from tcpyPI.pi import solve_temperature_from_entropy + + +PTOP = 50.0 +RD = constants.RD +EPS = constants.EPS +CKCD = 0.9 +V_REDUC = 0.8 + + +@njit(cache=True) +def cape_conv(TP, RP, PP, T, R, P, nlvl, conv): + """CAPE/TOB/LNB for one parcel under one convention. + Returns (E, TOB, LNB, flag): flag 1 ok, 0 improper parcel, 2 no-converge.""" + if (RP < 1e-6) or (TP < 200.0): + return 0.0, T[0], 0.0, 0 + + TPC = utilities.T_ktoC(TP) + ESP = utilities.es_cc(TPC) + EVP = utilities.ev(RP, PP) + RH = min(EVP / ESP, 1.0) + S = utilities.entropy_S(TP, RP, PP) + PLCL = utilities.e_pLCL(TP, RH, PP) + + b = np.empty(nlvl) + for j in range(nlvl): + if P[j] >= PLCL: + TG = TP * (P[j] / PP) ** (constants.RD / constants.CPD) + TLVR = utilities.Trho(TG, RP, RP) + b[j] = TLVR - utilities.Trho(T[j], R[j], R[j]) + else: + TG, RG, IFLAG = solve_temperature_from_entropy(S, P[j], RP, T[j]) + if IFLAG == 2: + return 0.0, T[0], P[0], 2 + TLVR = utilities.Trho(TG, RP, RG) + b[j] = TLVR - utilities.Trho(T[j], R[j], R[j]) + + # running signed work with per-crossing splits; region-end candidates + W = RD * (PP - P[0]) / (PP + P[0]) * b[0] # surface partial term + maxreg = nlvl + 1 + ps = np.empty(maxreg) # W at region ends + cp = np.empty(maxreg) # region-end pressures (crossing or top) + ct = np.empty(maxreg) # temperature at region ends (interp / top) + k = 0 + for j in range(1, nlvl): + b0 = b[j - 1] + b1 = b[j] + if (b0 > 0.0) != (b1 > 0.0): + if b0 * b1 < 0.0: + PC = (P[j] * b0 - P[j - 1] * b1) / (b0 - b1) + else: + PC = P[j] if b1 == 0.0 else P[j - 1] + W += RD * b0 * (P[j - 1] - PC) / (P[j - 1] + PC) + ps[k] = W + cp[k] = PC + ct[k] = (T[j - 1] * (PC - P[j]) + T[j] * (P[j - 1] - PC)) / (P[j - 1] - P[j]) + k += 1 + W += RD * b1 * (PC - P[j]) / (PC + P[j]) + else: + W += RD * (b0 + b1) * (P[j - 1] - P[j]) / (P[j] + P[j - 1]) + ps[k] = W + cp[k] = P[nlvl - 1] + ct[k] = T[nlvl - 1] + n = k + 1 + first_pos = b[0] > 0.0 + top_pos = b[nlvl - 1] > 0.0 + + E = 0.0 + TOB = T[0] + LNB = 0.0 + if conv == 0: + # legacy: topmost positive grid level among j>=1 (pcmin's INB loop + # excludes the launch node) + has_pos = False + for j in range(1, nlvl): + if b[j] > 0.0: + has_pos = True + break + if has_pos: + kk = n - 1 if top_pos else n - 2 + if kk >= 0 and ps[kk] > 0.0: + E = ps[kk] + TOB = ct[kk] + LNB = cp[kk] + elif conv == 1: + best = 0.0 + for kk in range(n): + if ps[kk] > best: + best = ps[kk] + TOB = ct[kk] + LNB = cp[kk] + E = best + if E == 0.0: + TOB = T[0] + LNB = 0.0 + elif conv == 2: + r1 = 0 if first_pos else 1 + if r1 < n and ps[r1] > 0.0: + E = ps[r1] + TOB = ct[r1] + LNB = cp[r1] + else: + r1 = 0 if first_pos else 1 + if r1 < n: + floor = ps[r1 - 1] if r1 >= 1 else 0.0 + best = -np.inf + bi = -1 + for kk in range(r1, n): + if ps[kk] > best: + best = ps[kk] + bi = kk + if ps[kk] <= floor and kk > r1: + break + if best > 0.0: + E = best + TOB = ct[bi] + LNB = cp[bi] + return E, TOB, LNB, 1 + + +@njit(cache=True) +def pi_conv(SSTC, MSL, T, R, P, nlvl, conv): + """pi() minimum-pressure loop under one CAPE convention. + Returns (VMAX, PMIN, IFL, TO, OTL, NP).""" + SSTK = utilities.T_Ctok(SSTC) + if SSTC <= 5.0 or SSTC > 100.0: + return np.nan, np.nan, 0, np.nan, np.nan, 0 + if np.min(T) <= 100.0 or np.max(T) - 273.15 > 100.0: + return np.nan, np.nan, 0, np.nan, np.nan, 0 + ES0 = utilities.es_cc(SSTC) + + IFL = 1 + CAPEA, _, _, fA = cape_conv(T[0], R[0], P[0], T, R, P, nlvl, conv) + if fA != 1: + IFL = fA + + NP = 0 + PM = 970.0 + PMOLD = PM + PNEW = 0.0 + TO = np.nan + OTL = np.nan + CAPEM = 0.0 + CAPEMS = 0.0 + RAT = 1.0 + TVAV = 300.0 + while np.abs(PNEW - PMOLD) > 0.5: + PP = min(PM, 1000.0) + RP = EPS * R[0] * MSL / (PP * (EPS + R[0]) - R[0] * MSL) + CAPEM, _, _, fM = cape_conv(T[0], RP, PP, T, R, P, nlvl, conv) + if fM != 1: + IFL = fM + RPS = utilities.rv(ES0, PP) + CAPEMS, TOMS, LNBS, fS = cape_conv(SSTK, RPS, PP, T, R, P, nlvl, conv) + if fS != 1: + IFL = fS + TO = TOMS + OTL = LNBS + RAT = SSTK / TO + TV0 = utilities.Trho(T[0], R[0], R[0]) + TVSST = utilities.Trho(SSTK, RPS, RPS) + TVAV = 0.5 * (TV0 + TVSST) + CAT = max((CAPEM - CAPEA) + 0.5 * CKCD * RAT * (CAPEMS - CAPEM), 0.0) + PNEW = MSL * np.exp(-CAT / (RD * TVAV)) + PMOLD = PM + PM = PNEW + NP += 1 + if (NP > 200) or (PM < 400.0): + return np.nan, np.nan, 2, np.nan, np.nan, NP + + CATFAC = 0.5 * (1.0 + 1.0 / constants.b) + CAT = max((CAPEM - CAPEA) + CKCD * RAT * CATFAC * (CAPEMS - CAPEM), 0.0) + PMIN = MSL * np.exp(-CAT / (RD * TVAV)) + FAC = max(0.0, CAPEMS - CAPEM) + VMAX = V_REDUC * np.sqrt(CKCD * RAT * FAC) + return VMAX, PMIN, IFL, TO, OTL, NP + + +@njit(parallel=True, cache=True) +def run_all(SST, MSL, TCs, Rs, P, nlvl, VMAX, PMIN, IFL, TO, OTL, NIT): + nprof = TCs.shape[0] + for i in prange(nprof): + T = TCs[i] + 273.15 + R = Rs[i] * 0.001 + for conv in range(4): + v, p, f, t, o, np_ = pi_conv(SST[i], MSL[i], T, R, P, nlvl, conv) + VMAX[i, conv] = v + PMIN[i, conv] = p + IFL[i, conv] = f + TO[i, conv] = t + OTL[i, conv] = o + NIT[i, conv] = np_ + + +def main(): + d = np.load(f"{SCRATCH}/profiles_converted.npz") + P_full = d["P"] + nlvl = int((P_full > PTOP).sum()) + P = P_full[:nlvl].copy() + TCs = d["TC"][:, :nlvl].copy() + Rs = d["R"][:, :nlvl].copy() + SST = d["sst_C"] + MSL = d["sp_hPa"] + nprof = TCs.shape[0] + + VMAX = np.full((nprof, 4), np.nan) + PMIN = np.full((nprof, 4), np.nan) + IFL = np.zeros((nprof, 4), dtype=np.int64) + TO = np.full((nprof, 4), np.nan) + OTL = np.full((nprof, 4), np.nan) + NIT = np.zeros((nprof, 4), dtype=np.int64) + + t0 = time.time() + run_all(SST[:50], MSL[:50], TCs[:50], Rs[:50], P, nlvl, + VMAX[:50], PMIN[:50], IFL[:50], TO[:50], OTL[:50], NIT[:50]) + print(f"warmup: {time.time()-t0:.1f}s") + t0 = time.time() + run_all(SST, MSL, TCs, Rs, P, nlvl, VMAX, PMIN, IFL, TO, OTL, NIT) + print(f"full run ({nprof} x 4 conventions): {time.time()-t0:.1f}s") + + np.savez_compressed(f"{SCRATCH}/pi_conv_results.npz", + VMAX=VMAX, PMIN=PMIN, IFL=IFL, TO=TO, OTL=OTL, NIT=NIT, + conv_names=np.array(["top", "max", "first", "reach_lfc"])) + print("saved pi_conv_results.npz") + + +if __name__ == "__main__": + main() diff --git a/era5_analysis/plot_curves.py b/era5_analysis/plot_curves.py new file mode 100644 index 0000000..4e2e274 --- /dev/null +++ b/era5_analysis/plot_curves.py @@ -0,0 +1,95 @@ +"""Spaghetti plots of buoyancy curves conditioned on topology.""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + + +BLUE = "#2a78d6" +DBLUE = "#104281" +INK = "#0b0b0b" +INK2 = "#52514e" +GRID = "#e5e4e0" + +plt.rcParams.update({ + "font.size": 9, "font.family": "serif", + "axes.edgecolor": INK2, "axes.labelcolor": INK, + "axes.titlesize": 10, "axes.titleweight": "bold", + "axes.spines.top": False, "axes.spines.right": False, + "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, + "xtick.color": INK2, "ytick.color": INK2, + "figure.dpi": 150, "savefig.bbox": "tight", +}) + +z = np.load(f"{SCRATCH}/topology_curves.npz") +P = z["P"] +YT = [1000, 850, 700, 500, 400, 300, 200, 150, 100, 70, 50, 30] + +groups = {} +for k in z.files: + if k == "P": + continue + pname, topo, share = k.split("|") + groups.setdefault(pname, []).append((float(share), topo, z[k])) + +TITLES = {"A": "Parcel A — environmental (ambient launch)", + "B": "Parcel B — eyewall (moisture-enriched, at converged $P_M$)", + "C": "Parcel C — saturated core (at converged $P_M$)"} + +for pname, items in groups.items(): + items.sort(key=lambda x: -x[0]) + ncol = min(3, len(items)) + nrow = int(np.ceil(len(items) / ncol)) + fig, axs = plt.subplots(nrow, ncol, figsize=(3.5 * ncol, 4.1 * nrow), + sharey=True, squeeze=False) + # x lower limit: most negative interior dip (before the curve's last + # positive level); the terminal stratospheric plunge never returns to + # positive and should not set the range. 10% margin. + lo = 0.0 + hi = 0.0 + for _, _, curves_ in items: + for c in curves_: + fin = np.isfinite(c) + pos = np.where(fin & (c > 0))[0] + if len(pos): + lo = min(lo, np.nanmin(c[: pos[-1] + 1])) + hi = max(hi, np.nanmax(c[fin]) if fin.any() else 0.0) + hi = min(hi, np.nanpercentile( + np.concatenate([c.ravel() for _, _, c in items]), 99.95)) + lo *= 1.1 + hi *= 1.05 + pad = 0.0 + for ax, (share, topo, curves) in zip(axs.flat, items): + for c in curves: + ax.plot(c, P, color=BLUE, alpha=0.10, lw=0.7, zorder=3) + ax.plot(np.nanmedian(curves, axis=0), P, color=DBLUE, lw=1.8, zorder=5, + label="pointwise median") + ax.axvline(0, color=INK, lw=0.9, zorder=4) + ax.set_yscale("log") + ax.set_ylim(1010, 28) + ax.set_yticks(YT) + ax.set_yticklabels([str(t) for t in YT]) + ax.minorticks_off() + ax.set_xlim(lo - pad, hi + pad) + ax.set_title(f"{topo} ({share:.2f}%)") + ax.set_xlabel("buoyancy [K]") + for ax in axs[:, 0]: + ax.set_ylabel("pressure [hPa], log scale (up = higher)") + for ax in axs.flat[len(items):]: + ax.set_visible(False) + axs.flat[0].legend(loc="upper left", frameon=False, fontsize=8) + fig.suptitle(f"{TITLES[pname]} — {sum(len(c) for _, _, c in items)} sampled buoyancy profiles by topology" + " (x-axis clipped to the interior-dip range)", + fontsize=11, y=1.005) + fig.tight_layout() + fig.savefig(f"{SCRATCH}/curves_parcel_{pname}.png", dpi=140) + fig.savefig(f"{SCRATCH}/curves_parcel_{pname}.pdf") + plt.close(fig) + print(f"saved curves_parcel_{pname}.png/.pdf") diff --git a/era5_analysis/prep_artifact_data.py b/era5_analysis/prep_artifact_data.py new file mode 100644 index 0000000..f0da802 --- /dev/null +++ b/era5_analysis/prep_artifact_data.py @@ -0,0 +1,266 @@ +"""Assemble the JSON payload for the sounding-scope artifact. + +The per-profile statistics table covers the FULL population (every column, +binary-packed little-endian, gzipped, base64) so that counts, class shares, +densities and the dVMAX scatters are exact under any marginalization. The +flip-book player carries curves for a ~20k stratified sample (curves are the +expensive part; the sample only has to be dense enough to animate). +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import base64 +import gzip +import json +import sys + +import numpy as np + +sys.path.insert(0, SCRATCH) +sys.path.insert(0, _SRC) + +from scan_final import buoyancy, EPS, PTOP +from tcpyPI import utilities + +# materialize everything: NpzFile decompresses the whole array per access +_d = np.load(f"{SCRATCH}/profiles_converted.npz") +d = {k: _d[k] for k in _d.files} +_z = np.load(f"{SCRATCH}/parcel_topology_final.npz") +z = {k: _z[k] for k in ("A", "B", "C", "PM")} +_pc = np.load(f"{SCRATCH}/pi_conv_results.npz") +PI = {k: _pc[k] for k in ("VMAX", "PMIN", "IFL")} # (n, 4): top,max,first,reach +P_full = d["P"] +nlvl = int((P_full > PTOP).sum()) +P = P_full[:nlvl].copy() +# per-parcel populations: A is defined everywhere; B/C wherever the +# max-work pressure iteration converged (SST > 5 C) +mA = np.isfinite(z["A"][:, 0]) +mBC = np.isfinite(z["PM"]) & np.isfinite(z["B"][:, 0]) & np.isfinite(z["C"][:, 0]) +popmask = {"A": mA, "B": mBC, "C": mBC} +pop = {p: np.where(popmask[p])[0] for p in "ABC"} +print("populations:", {p: len(pop[p]) for p in "ABC"}) + +t64 = d["time"].astype("datetime64[s]") +years = t64.astype("datetime64[Y]").astype(int) + 1970 +doy = (t64.astype("datetime64[D]") + - t64.astype("datetime64[Y]").astype("datetime64[D]")).astype(int) + 1 + + +def topo_str(fs, n): + if fs == 0: + return "0" + c0 = "+" if fs > 0 else "-" + return (("+-" if c0 == "+" else "-+") * ((int(n) + 1) // 2))[: int(n)] + + +# per-parcel class tables (top 7 by count within tc, then "other") +tables = {} +labels = {} +allstrs = {} +for pname in "ABC": + o = z[pname] + rows = pop[pname] + strs = np.array([topo_str(o[i, 0], o[i, 1]) for i in rows]) + allstrs[pname] = strs + uniq, cnt = np.unique(strs, return_counts=True) + order = np.argsort(-cnt) + top = [str(uniq[k]) for k in order[:7]] + tables[pname] = {t: i for i, t in enumerate(top)} + labels[pname] = [[str(uniq[k]), int(cnt[k]), + float(round(100.0 * cnt[k] / len(rows), 2))] + for k in order[:7]] + other = int(len(rows) - cnt[order[:7]].sum()) + labels[pname].append(["other", other, float(round(100.0 * other / len(rows), 2))]) + + +def cls(pname, i): + if not popmask[pname][i]: + return -1 + return tables[pname].get(topo_str(z[pname][i, 0], z[pname][i, 1]), 7) + + +rng = np.random.default_rng(42) + +# ---- player pool: ensure every top class of every parcel is represented ---- +pool = set(rng.choice(pop["A"], 14000, replace=False).tolist()) +# force-include the wild legacy non-convergence columns (IFL_top=2, IFL_max=1) +legfail_rows = np.where((PI["IFL"][:, 0] == 2) & (PI["IFL"][:, 1] == 1) & mA)[0] +pool.update(legfail_rows.tolist()) +print("legacy-failure columns force-included:", len(legfail_rows)) +# force flip-book coverage of the extreme convention gaps (capped) +with np.errstate(invalid="ignore"): + _dv = PI["VMAX"][:, 1] - PI["VMAX"][:, 0] +for _m in (np.isfinite(_dv) & (_dv > 8.0) & mA, + np.isfinite(_dv) & (_dv < -2.0) & mA): + _rows = np.where(_m)[0] + if len(_rows) > 400: + _rows = rng.choice(_rows, 400, replace=False) + pool.update(_rows.tolist()) +for pname in "ABC": + for t in tables[pname]: + rows = pop[pname][allstrs[pname] == t] + take = min(600, len(rows)) + pool.update(rng.choice(rows, take, replace=False).tolist()) +pool = np.array(sorted(pool)) +print("player pool:", len(pool)) + +# curve rectangles: (3, npool, nlvl) int16 in centi-J/kg (0.01 precision, +# same as the old JSON rounding); zeros where B/C undefined -- gzip erases +# the constant runs, and the JS side rebuilds [] from the class sentinel +npool = len(pool) +curvebuf = np.zeros((3, npool, nlvl), np.float64) +b = np.empty(nlvl) +for k, i in enumerate(pool): + T = d["TC"][i, :nlvl] + 273.15 + R = d["R"][i, :nlvl] * 0.001 + buoyancy(T[0], R[0], P[0], T, R, P, nlvl, b) + curvebuf[0, k] = b + if mBC[i]: + PP = min(z["PM"][i], 1000.0) + MSL = d["sp_hPa"][i] + RP = EPS * R[0] * MSL / (PP * (EPS + R[0]) - R[0] * MSL) + buoyancy(T[0], RP, PP, T, R, P, nlvl, b) + curvebuf[1, k] = b + ES0 = utilities.es_cc(d["sst_C"][i]) + buoyancy(d["sst_C"][i] + 273.15, utilities.rv(ES0, PP), PP, T, R, P, nlvl, b) + curvebuf[2, k] = b +curves = {pn: curvebuf[j] for j, pn in enumerate("ABC")} +curve_i2 = np.round(curvebuf * 100).astype(" 180, d["lon"][pool] - 360, d["lon"][pool]) +player = { + "lat": [round(float(x), 2) for x in d["lat"][pool]], + "lon": [round(float(x), 2) for x in lonw], + "year": [int(x) for x in years[pool]], + "doy": [int(x) for x in doy[pool]], + "hour": [int(x) for x in ((d["time"][pool] // 3600) % 24)], + "sst": [round(float(x), 1) for x in d["sst_C"][pool]], +} +for pname in "ABC": + m = popmask[pname] + o = z[pname] + player[f"c{pname}"] = [int(cls(pname, i)) for i in pool] + for tag, col in (("Etop", 2), ("Emax", 3), ("Ereach", 5)): + player[f"{tag}{pname}"] = [int(round(o[i, col])) if m[i] else 0 for i in pool] + for tag, col in (("Ltop", 8), ("Lmax", 9), ("Lreach", 11)): + player[f"{tag}{pname}"] = [round(float(o[i, col]), 1) if m[i] else 0 for i in pool] + player[f"topo{pname}"] = [topo_str(o[i, 0], o[i, 1]) if m[i] else "" + for i in pool] +# PI per convention for the column (ptop=50 tcpyPI convention; -1 = missing) +for k, tag in zip((0, 1, 3), ("top", "max", "reach")): + player[f"V{tag}"] = [round(float(PI["VMAX"][i, k]), 1) + if np.isfinite(PI["VMAX"][i, k]) else -1 for i in pool] +player["Ptop"] = [round(float(PI["PMIN"][i, 0]), 1) + if np.isfinite(PI["PMIN"][i, 0]) else -1 for i in pool] +player["Pmax"] = [round(float(PI["PMIN"][i, 1]), 1) + if np.isfinite(PI["PMIN"][i, 1]) else -1 for i in pool] +player["legfail"] = [1 if (PI["IFL"][i, 0] == 2 and PI["IFL"][i, 1] == 1) else 0 + for i in pool] + +# ---- statistics table: FULL population, binary-packed ---- +N = len(mA) +lon_all = np.where(d["lon"] > 180, d["lon"] - 360, d["lon"]) +hours_all = (d["time"] // 3600) % 24 + +clsv, gv = {}, {} +for pname in "ABC": + o = z[pname] + ci = np.full(N, -1, np.int8) + lut = tables[pname] + ci[pop[pname]] = np.array([lut.get(s, 7) for s in allstrs[pname]], np.int8) + clsv[pname] = ci + # disagreement bitmask: 1 = clamp-flip (E_top=0 while E_max>0); + # 2 = |E_top-E_max| > 1 J/kg; 4 = |E_reach-E_max| > 1 J/kg + with np.errstate(invalid="ignore"): + g = ((o[:, 2] == 0.0) & (o[:, 3] > 0.0)).astype(np.uint8) + g |= (np.abs(o[:, 2] - o[:, 3]) > 1.0).astype(np.uint8) << 1 + g |= (np.abs(o[:, 5] - o[:, 3]) > 1.0).astype(np.uint8) << 2 + g[~popmask[pname]] = 0 + gv[pname] = g + +# VMAX per convention, deci-m/s; -10 encodes missing (decodes to -1.0) +vq = {} +for tag, k in (("top", 0), ("max", 1), ("reach", 3)): + V = PI["VMAX"][:, k] + vq[tag] = np.where(np.isfinite(V), np.round(V * 10), -10).astype(" 8.0) + dv_lo = np.isfinite(dv_all) & (dv_all < -2.0) +gpi_all = np.where(legf, 1, np.where(dv_hi, 2, + np.where(dv_lo, 3, 0))).astype(np.uint8) +print("dv>8:", int(dv_hi.sum()), " dv<-2:", int(dv_lo.sum())) + +# int16 columns first so every offset stays 2-byte aligned +bincols = [ + ("lat", "i2", 100, np.round(d["lat"] * 100).astype(" " + f"{len(mapbin['b64'])/1e6:.1f} MB b64gz; " + f"curves: {curve_i2.nbytes/1e6:.1f} MB raw -> " + f"{len(curvebin['b64'])/1e6:.1f} MB b64gz") + +# per-parcel x-limits for the scope (tropospheric percentiles) +# x-limits: lower bound = the most negative any curve goes BEFORE its last +# positive level (interior dips only; the terminal stratospheric plunge never +# returns to positive and should not set the range), with a 10% margin. +xlims = {} +poolBC = mBC[pool] +for pname in "ABC": + arr = curves[pname] if pname == "A" else curves[pname][poolBC] + lo = 0.0 + for c in arr: + pos = np.where(c > 0)[0] + if len(pos) == 0: + continue + m = c[: pos[-1] + 1].min() + if m < lo: + lo = m + hi = float(np.percentile(arr, 99.9)) + xlims[pname] = [round(float(lo * 1.1), 1), round(hi * 1.05 + 0.5, 1)] + +payload = { + "P": [float(x) for x in P], + "player": player, + "mapbin": mapbin, + "curvebin": curvebin, + "classes": labels, + "xlims": xlims, + "years": [int(years[pop["A"]].min()), int(years[pop["A"]].max())], + "npop": {p: int(len(pop[p])) for p in "ABC"}, + "ntc": int(len(pop["A"])), +} +js = json.dumps(payload, separators=(",", ":")) +open(f"{SCRATCH}/artifact_data.json", "w").write(js) +print(f"payload: {len(js)/1e6:.2f} MB; years {payload['years']}") +for pname in "ABC": + print(pname, labels[pname][:4]) diff --git a/era5_analysis/scan_final.py b/era5_analysis/scan_final.py new file mode 100644 index 0000000..e659535 --- /dev/null +++ b/era5_analysis/scan_final.py @@ -0,0 +1,321 @@ +"""Final topology/convention scan for all three PI parcels. + +Topology convention: leading sign = sign of buoyancy infinitesimally above +launch (= sign of first nonzero b for the PWL interpolant); no zero-width +leading regions. Conventions are computed candidate-based (down-crossings of +the running work W), so they are unaffected by the labeling change: + E_top : W at the last down-crossing / top (legacy; requires a positive + grid level above the launch node), clamped at 0 + E_max : max of W over all candidates (>= 0) + E_first : W at the first down-crossing, clamped at 0 + E_reach : lifted-ballistic (floor = W entering the first positive region) + +Parcel A at ambient; parcels B, C at the converged max-work P_M. +Output: parcel_topology_final.npz +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import sys +import time + +sys.path.insert(0, _SRC) + +import numpy as np +from numba import njit, prange + +from tcpyPI import constants, utilities +from tcpyPI.pi import solve_temperature_from_entropy + +PTOP = 20.0 # retain levels with P > 20 hPa -> top 30 hPa: the minimal +# ceiling observing every buoyancy crossing in the 569k-column sample (highest: +# 47.6 hPa, saturated core parcel, interpolated between the 50 and 30 hPa +# nodes). At this ceiling zero profiles clip and zero entropy solves fail (the +# 30-20 hPa layer caused 21 failures when retained). +RD = constants.RD +EPS = constants.EPS +CKCD = 0.9 + + +@njit(cache=True) +def buoyancy(TP, RP, PP, T, R, P, nlvl, b): + if (RP < 1e-6) or (TP < 200.0): + return 0 + TPC = utilities.T_ktoC(TP) + RH = min(utilities.ev(RP, PP) / utilities.es_cc(TPC), 1.0) + S = utilities.entropy_S(TP, RP, PP) + PLCL = utilities.e_pLCL(TP, RH, PP) + for j in range(nlvl): + if P[j] >= PLCL: + TG = TP * (P[j] / PP) ** (constants.RD / constants.CPD) + b[j] = utilities.Trho(TG, RP, RP) - utilities.Trho(T[j], R[j], R[j]) + else: + TG, RG, IFLAG = solve_temperature_from_entropy(S, P[j], RP, T[j]) + if IFLAG == 2: + return 2 + b[j] = utilities.Trho(TG, RP, RG) - utilities.Trho(T[j], R[j], R[j]) + return 1 + + +@njit(cache=True) +def summarize(b, P, PP, nlvl, out): + """out = [first_sign(+1/-1/0), n, E_top, E_max, E_first, E_reach, + clamped, clipped, LNB_top, LNB_max, LNB_first, LNB_reach] + (LNB pressures in hPa; 0 where the convention returns zero CAPE)""" + # ---- topology under the perturbation convention ---- + # effective sign sequence: skip zero values for the leading sign; a zero + # value elsewhere adopts the '-' side as before (b>0 is '+') + j0 = 0 + while j0 < nlvl and b[j0] == 0.0: + j0 += 1 + if j0 == nlvl: + out[0] = 0.0 + out[1] = 1.0 + return + first_pos = b[j0] > 0.0 + n = 1 + prev = first_pos + for j in range(j0 + 1, nlvl): + cur = b[j] > 0.0 + if cur != prev: + n += 1 + prev = cur + out[0] = 1.0 if first_pos else -1.0 + out[1] = n + + # ---- running work; candidates at down-crossings (+ -> -) and top ---- + W = RD * (PP - P[0]) / (PP + P[0]) * b[0] + ncand = 0 + cand_W = np.empty(nlvl + 1) + cand_P = np.empty(nlvl + 1) + up_W = 0.0 # W entering the first positive region + seen_pos = b[0] > 0.0 + if seen_pos: + up_W = 0.0 # launch is inside a positive region (parcels B, C) + have_up = seen_pos + for j in range(1, nlvl): + b0 = b[j - 1] + b1 = b[j] + if (b0 > 0.0) != (b1 > 0.0): + if b0 * b1 < 0.0: + PC = (P[j] * b0 - P[j - 1] * b1) / (b0 - b1) + else: + PC = P[j] if b1 == 0.0 else P[j - 1] + W += RD * b0 * (P[j - 1] - PC) / (P[j - 1] + PC) + if b0 > 0.0: + cand_W[ncand] = W # down-crossing candidate + cand_P[ncand] = PC + ncand += 1 + else: + if not have_up: + up_W = W # first entry into a + region + have_up = True + W += RD * b1 * (PC - P[j]) / (PC + P[j]) + else: + W += RD * (b0 + b1) * (P[j - 1] - P[j]) / (P[j] + P[j - 1]) + top_pos = b[nlvl - 1] > 0.0 + if top_pos: + cand_W[ncand] = W # profile-top candidate (clipped) + cand_P[ncand] = P[nlvl - 1] + ncand += 1 + out[7] = 1.0 if top_pos else 0.0 + + if ncand == 0: + return + + # E_max over all candidates (launch, W=0, implicit) + best = 0.0 + for kk in range(ncand): + if cand_W[kk] > best: + best = cand_W[kk] + out[9] = cand_P[kk] # LNB_max + out[3] = best + + # E_top: legacy requires a positive grid level above the launch node + has_pos_above = False + for j in range(1, nlvl): + if b[j] > 0.0: + has_pos_above = True + break + if has_pos_above: + raw = cand_W[ncand - 1] + if raw < 0.0: + out[6] = 1.0 + else: + out[2] = raw + out[8] = cand_P[ncand - 1] # LNB_top + + # E_first + if cand_W[0] > 0.0: + out[4] = cand_W[0] + out[10] = cand_P[0] # LNB_first + + # E_reach (lifted-ballistic): walk candidates; between candidate k and + # k+1 the parcel traverses a negative region whose minimum W is the entry + # of the next + region; we approximate the stall check with the next + # up-crossing W (monotone within regions makes this exact). + # Reconstruct up-crossing Ws in a second pass for correctness. + floor = up_W + bestr = 0.0 + W2 = RD * (PP - P[0]) / (PP + P[0]) * b[0] + entered = b[0] > 0.0 + stalled = False + for j in range(1, nlvl): + b0 = b[j - 1] + b1 = b[j] + if (b0 > 0.0) != (b1 > 0.0): + if b0 * b1 < 0.0: + PC = (P[j] * b0 - P[j - 1] * b1) / (b0 - b1) + else: + PC = P[j] if b1 == 0.0 else P[j - 1] + W2 += RD * b0 * (P[j - 1] - PC) / (P[j - 1] + PC) + if b0 > 0.0: + if entered and not stalled and W2 > bestr: + bestr = W2 + out[11] = PC # LNB_reach + else: + if not entered: + entered = True + elif entered and W2 <= floor: + stalled = True + W2 += RD * b1 * (PC - P[j]) / (PC + P[j]) + else: + W2 += RD * (b0 + b1) * (P[j - 1] - P[j]) / (P[j] + P[j - 1]) + if top_pos and entered and not stalled and W2 > bestr: + bestr = W2 + out[11] = P[nlvl - 1] + if bestr > 0.0: + out[5] = bestr + else: + out[5] = 0.0 + out[11] = 0.0 + + +@njit(cache=True) +def pm_converged(SSTC, MSL, T, R, P, nlvl): + """Max-work pi loop; returns converged PM (RMW-level pressure) or NaN.""" + if SSTC <= 5.0 or SSTC > 100.0: + return np.nan + SSTK = utilities.T_Ctok(SSTC) + ES0 = utilities.es_cc(SSTC) + b = np.empty(nlvl) + tmp = np.empty(12) + + # E_max-only cape via summarize (out[3]) + def _capemax(TP, RP, PP): + for q in range(12): + tmp[q] = 0.0 + f = buoyancy(TP, RP, PP, T, R, P, nlvl, b) + if f != 1: + return 0.0, T[0], -1.0 + summarize(b, P, PP, nlvl, tmp) + return tmp[3], 0.0, 1.0 + + CAPEA, _, _ = _capemax(T[0], R[0], P[0]) + NP = 0 + PM = 970.0 + PMOLD = PM + PNEW = 0.0 + while np.abs(PNEW - PMOLD) > 0.5: + PP = min(PM, 1000.0) + RP = EPS * R[0] * MSL / (PP * (EPS + R[0]) - R[0] * MSL) + CAPEM, _, _ = _capemax(T[0], RP, PP) + RPS = utilities.rv(ES0, PP) + # saturated call: need TOB -> compute buoyancy then find argmax temp: + f = buoyancy(SSTK, RPS, PP, T, R, P, nlvl, b) + if f != 1: + return np.nan + for q in range(12): + tmp[q] = 0.0 + summarize(b, P, PP, nlvl, tmp) + CAPEMS = tmp[3] + # outflow temperature: T at coldest buoyancy-weighted... use simple + # proxy: T at the level nearest the LNB is not tracked here; RAT via + # tropopause-min temperature of the environment as in BE02 spirit. + # For PM convergence only the CAT magnitude matters; use env min T. + TOmin = T[0] + for j in range(nlvl): + if T[j] < TOmin: + TOmin = T[j] + RAT = SSTK / TOmin + TV0 = utilities.Trho(T[0], R[0], R[0]) + TVSST = utilities.Trho(SSTK, RPS, RPS) + TVAV = 0.5 * (TV0 + TVSST) + CAT = max((CAPEM - CAPEA) + 0.5 * CKCD * RAT * (CAPEMS - CAPEM), 0.0) + PNEW = MSL * np.exp(-CAT / (RD * TVAV)) + PMOLD = PM + PM = PNEW + NP += 1 + if (NP > 200) or (PM < 400.0): + return np.nan + return PM + + +@njit(parallel=True, cache=True) +def run(SST, MSL, TCs, Rs, P, nlvl, PMs, outA, outB, outC): + nprof = TCs.shape[0] + for i in prange(nprof): + T = TCs[i] + 273.15 + R = Rs[i] * 0.001 + b = np.empty(nlvl) + # parcel A at ambient + f = buoyancy(T[0], R[0], P[0], T, R, P, nlvl, b) + if f == 1: + summarize(b, P, P[0], nlvl, outA[i]) + else: + outA[i, 0] = np.nan + # converged PM (max-work) + PM = pm_converged(SST[i], MSL[i], T, R, P, nlvl) + PMs[i] = PM + if not np.isfinite(PM): + outB[i, 0] = np.nan + outC[i, 0] = np.nan + continue + PP = min(PM, 1000.0) + RP = EPS * R[0] * MSL[i] / (PP * (EPS + R[0]) - R[0] * MSL[i]) + f = buoyancy(T[0], RP, PP, T, R, P, nlvl, b) + if f == 1: + summarize(b, P, PP, nlvl, outB[i]) + else: + outB[i, 0] = np.nan + SSTK = SST[i] + 273.15 + ES0 = utilities.es_cc(SST[i]) + RPS = utilities.rv(ES0, PP) + f = buoyancy(SSTK, RPS, PP, T, R, P, nlvl, b) + if f == 1: + summarize(b, P, PP, nlvl, outC[i]) + else: + outC[i, 0] = np.nan + + +def main(): + d = np.load(f"{SCRATCH}/profiles_converted.npz") + P_full = d["P"] + nlvl = int((P_full > PTOP).sum()) + P = P_full[:nlvl].copy() + TCs = d["TC"][:, :nlvl].copy() + Rs = d["R"][:, :nlvl].copy() + nprof = TCs.shape[0] + PMs = np.full(nprof, np.nan) + outA = np.zeros((nprof, 12)) + outB = np.zeros((nprof, 12)) + outC = np.zeros((nprof, 12)) + + t0 = time.time() + run(d["sst_C"][:50], d["sp_hPa"][:50], TCs[:50], Rs[:50], P, nlvl, + PMs[:50], outA[:50], outB[:50], outC[:50]) + run(d["sst_C"], d["sp_hPa"], TCs, Rs, P, nlvl, PMs, outA, outB, outC) + print(f"scan: {time.time()-t0:.1f}s") + np.savez_compressed(f"{SCRATCH}/parcel_topology_final.npz", + A=outA, B=outB, C=outC, PM=PMs, + cols=np.array(["first_sign", "n", "E_top", "E_max", + "E_first", "E_reach", "clamped", "clipped", + "LNB_top", "LNB_max", "LNB_first", "LNB_reach"])) + print("saved parcel_topology_final.npz") + + +if __name__ == "__main__": + main() diff --git a/era5_analysis/scan_pi_parcels.py b/era5_analysis/scan_pi_parcels.py new file mode 100644 index 0000000..e04c049 --- /dev/null +++ b/era5_analysis/scan_pi_parcels.py @@ -0,0 +1,236 @@ +"""Topology scan for the two PI parcels (eyewall, saturated core) at each +profile's converged max-work P_M, plus a quadrature-scheme comparison +(current hybrid vs self-consistent linear-in-ln(p)) for the env parcel. +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import sys +import time + +sys.path.insert(0, _SRC) + +import numpy as np +from numba import njit, prange + +from tcpyPI import constants, utilities +from tcpyPI.pi import solve_temperature_from_entropy + + +PTOP = 50.0 +RD = constants.RD +EPS = constants.EPS +CKCD = 0.9 + + +@njit(cache=True) +def buoyancy(TP, RP, PP, T, R, P, nlvl, b): + """Fill b; return 1 ok / 0 improper / 2 entropy no-converge.""" + if (RP < 1e-6) or (TP < 200.0): + return 0 + TPC = utilities.T_ktoC(TP) + RH = min(utilities.ev(RP, PP) / utilities.es_cc(TPC), 1.0) + S = utilities.entropy_S(TP, RP, PP) + PLCL = utilities.e_pLCL(TP, RH, PP) + for j in range(nlvl): + if P[j] >= PLCL: + TG = TP * (P[j] / PP) ** (constants.RD / constants.CPD) + b[j] = utilities.Trho(TG, RP, RP) - utilities.Trho(T[j], R[j], R[j]) + else: + TG, RG, IFLAG = solve_temperature_from_entropy(S, P[j], RP, T[j]) + if IFLAG == 2: + return 2 + b[j] = utilities.Trho(TG, RP, RG) - utilities.Trho(T[j], R[j], R[j]) + return 1 + + +@njit(cache=True) +def summarize(b, P, PP, nlvl, lnp_mode): + """Topology summary for one buoyancy profile. + lnp_mode: 0 = current hybrid quadrature, 1 = linear-in-ln(p) consistent. + Returns (n, n_plus, E_top, E_max, clamped, clipped, first_pos).""" + if lnp_mode == 0: + W = RD * (PP - P[0]) / (PP + P[0]) * b[0] + else: + W = RD * np.log(PP / P[0]) * 0.5 * (b[0] + b[0]) + ps = np.empty(nlvl + 1) + k = 0 + for j in range(1, nlvl): + b0 = b[j - 1] + b1 = b[j] + if (b0 > 0.0) != (b1 > 0.0): + if lnp_mode == 0: + if b0 * b1 < 0.0: + PC = (P[j] * b0 - P[j - 1] * b1) / (b0 - b1) + else: + PC = P[j] if b1 == 0.0 else P[j - 1] + W += RD * b0 * (P[j - 1] - PC) / (P[j - 1] + PC) + ps[k] = W + W += RD * b1 * (PC - P[j]) / (PC + P[j]) + else: + l0 = np.log(P[j - 1]) + l1 = np.log(P[j]) + if b0 * b1 < 0.0: + lc = (l1 * b0 - l0 * b1) / (b0 - b1) + else: + lc = l1 if b1 == 0.0 else l0 + W += RD * 0.5 * b0 * (l0 - lc) + ps[k] = W + W += RD * 0.5 * b1 * (lc - l1) + k += 1 + else: + if lnp_mode == 0: + W += RD * (b0 + b1) * (P[j - 1] - P[j]) / (P[j] + P[j - 1]) + else: + W += RD * 0.5 * (b0 + b1) * np.log(P[j - 1] / P[j]) + ps[k] = W + n = k + 1 + first_pos = 1 if b[0] > 0.0 else 0 + top_pos = b[nlvl - 1] > 0.0 + n_plus = (n + (1 if first_pos == 1 else -1)) // 2 + (1 if (first_pos == 0 and top_pos) else 0) + # count '+' regions robustly: alternating from first_pos over n regions + n_plus = 0 + for kk in range(n): + pos = (kk % 2 == 0) == (first_pos == 1) + if pos: + n_plus += 1 + best = 0.0 + for kk in range(n): + if ps[kk] > best: + best = ps[kk] + E_max = best + clamped = 0 + if top_pos: + raw = ps[n - 1] + elif n >= 2: + raw = ps[n - 2] + else: + raw = 0.0 + has_pos_above = False + for j in range(1, nlvl): + if b[j] > 0.0: + has_pos_above = True + break + if not has_pos_above: + E_top = 0.0 + elif raw < 0.0: + E_top = 0.0 + clamped = 1 + else: + E_top = raw + clipped = 1 if top_pos else 0 + return n, n_plus, E_top, E_max, clamped, clipped, first_pos + + +@njit(parallel=True, cache=True) +def run(SST, MSL, PMs, TCs, Rs, P, nlvl, out_eye, out_core, out_env_cmp): + nprof = TCs.shape[0] + for i in prange(nprof): + T = TCs[i] + 273.15 + R = Rs[i] * 0.001 + b = np.empty(nlvl) + SSTK = SST[i] + 273.15 + # env-parcel quadrature comparison (hybrid vs lnp-consistent) + f = buoyancy(T[0], R[0], P[0], T, R, P, nlvl, b) + if f == 1: + n0, np0, Et0, Em0, cl0, _, _ = summarize(b, P, P[0], nlvl, 0) + n1, np1, Et1, Em1, cl1, _, _ = summarize(b, P, P[0], nlvl, 1) + out_env_cmp[i, 0] = Em0 + out_env_cmp[i, 1] = Em1 + out_env_cmp[i, 2] = Et0 + out_env_cmp[i, 3] = Et1 + if not np.isfinite(PMs[i]): + out_eye[i, 0] = -1.0 + out_core[i, 0] = -1.0 + continue + PP = min(PMs[i], 1000.0) + # eyewall parcel: enriched boundary-layer air at PP + RP = EPS * R[0] * MSL[i] / (PP * (EPS + R[0]) - R[0] * MSL[i]) + f = buoyancy(T[0], RP, PP, T, R, P, nlvl, b) + if f == 1: + n, npl, Et, Em, cl, cp, fp = summarize(b, P, PP, nlvl, 0) + out_eye[i, 0] = n + out_eye[i, 1] = npl + out_eye[i, 2] = Et + out_eye[i, 3] = Em + out_eye[i, 4] = cl + out_eye[i, 5] = cp + out_eye[i, 6] = fp + else: + out_eye[i, 0] = -float(f + 1) + # saturated core parcel + ES0 = utilities.es_cc(SST[i]) + RPS = utilities.rv(ES0, PP) + f = buoyancy(SSTK, RPS, PP, T, R, P, nlvl, b) + if f == 1: + n, npl, Et, Em, cl, cp, fp = summarize(b, P, PP, nlvl, 0) + out_core[i, 0] = n + out_core[i, 1] = npl + out_core[i, 2] = Et + out_core[i, 3] = Em + out_core[i, 4] = cl + out_core[i, 5] = cp + out_core[i, 6] = fp + else: + out_core[i, 0] = -float(f + 1) + + +def main(): + d = np.load(f"{SCRATCH}/profiles_converted.npz") + r = np.load(f"{SCRATCH}/pi_conv_results.npz") + P_full = d["P"] + nlvl = int((P_full > PTOP).sum()) + P = P_full[:nlvl].copy() + TCs = d["TC"][:, :nlvl].copy() + Rs = d["R"][:, :nlvl].copy() + # converged eye pressure under max-work as the P_M proxy; PMIN is the eye + # pressure -- close enough to RMW pressure for a population scan, and + # available per profile. Only rows converged under max (IFL==1). + PMs = np.where(r["IFL"][:, 1] == 1, r["PMIN"][:, 1], np.nan) + + nprof = TCs.shape[0] + out_eye = np.full((nprof, 7), np.nan) + out_core = np.full((nprof, 7), np.nan) + out_env_cmp = np.full((nprof, 4), np.nan) + + t0 = time.time() + run(d["sst_C"][:50], d["sp_hPa"][:50], PMs[:50], TCs[:50], Rs[:50], P, nlvl, + out_eye[:50], out_core[:50], out_env_cmp[:50]) + run(d["sst_C"], d["sp_hPa"], PMs, TCs, Rs, P, nlvl, out_eye, out_core, out_env_cmp) + print(f"scan: {time.time()-t0:.1f}s") + + np.savez_compressed(f"{SCRATCH}/pi_parcel_topology.npz", + eye=out_eye, core=out_core, env_cmp=out_env_cmp) + + # ---- quadrature comparison ---- + ok = np.isfinite(out_env_cmp[:, 0]) + dEm = np.abs(out_env_cmp[ok, 0] - out_env_cmp[ok, 1]) + base = np.maximum(out_env_cmp[ok, 0], 50.0) + print("\n=== hybrid vs linear-in-ln(p) quadrature (env parcel, E_max) ===") + print(f" |dE_max|: median {np.median(dEm):.3f} p99 {np.percentile(dEm,99):.2f} " + f"max {dEm.max():.2f} J/kg; rel max {np.max(dEm/base)*100:.2f}%") + + # ---- PI-parcel topology stats (TC-relevant subset) ---- + tc = (d["sst_C"] >= 26) & (d["sp_hPa"] >= 1000) & (r["IFL"][:, 1] == 1) + for name, out in [("eyewall", out_eye), ("core(sat)", out_core)]: + v = out[tc] + okv = v[:, 0] > 0 + n = v[okv, 0].astype(np.int64) + npl = v[okv, 1].astype(np.int64) + print(f"\n=== {name} parcel at converged P_M (TC-relevant, n={okv.sum()}) ===") + print(f" launches buoyant (b_surface>0): {100*v[okv,6].mean():.1f}%") + print(f" candidate-LNB count (+ regions): 1: {100*(npl==1).mean():.1f}% " + f"2: {100*(npl==2).mean():.1f}% >=3: {100*(npl>=3).mean():.1f}% 0: {100*(npl==0).mean():.2f}%") + print(f" E_top clamped to 0 while E_max>0: " + f"{100*((v[okv,4]==1)&(v[okv,3]>0)).mean():.2f}%") + dd = np.abs(v[okv, 2] - v[okv, 3]) + print(f" |E_top-E_max|: >1 J/kg {100*(dd>1).mean():.2f}% >10 {100*(dd>10).mean():.2f}% " + f"max {dd.max():.1f} J/kg") + print(f" clipped at profile top (b>0 at 70 hPa): {100*v[okv,5].mean():.2f}%") + + +if __name__ == "__main__": + main() diff --git a/era5_analysis/sounding_scope_template.html b/era5_analysis/sounding_scope_template.html new file mode 100644 index 0000000..7ae8193 --- /dev/null +++ b/era5_analysis/sounding_scope_template.html @@ -0,0 +1,1086 @@ +Sounding Scope — buoyancy topologies + + +
+

Sounding Scope

+ buoyancy-profile flipbook over global ERA5 ocean columns · parcels B/C on the with a converged storm pressure (SST > 5 °C) · conditioned on sign topology +
+ +
+ + + + + + + + +
+
+
+
+ +
+
+

buoyancy vs log-pressure

+ +
+
+
+

position · density of filtered population

+ +
+
+
+

day of year

+ +
+
+

hour of day (UTC)

+ +
+
+
+

year

+ +
+
+
+
+

ΔVMAX: max − top

+ +
+
+

ΔVMAX: reach − top

+ +
+
+
+
+ +
+ Each flash is one real reanalysis column: its lifted-parcel buoyancy curve (left), + its location, and its date. The shaded layers show where and when the currently + filtered topology lives: a 2°-binned density on the map, a 10-day rose on the + ring, and yearly counts on the track. Curve and markers take the color of the + profile's topology class. LNB conventions, named consistently throughout as + top (the legacy pcmin rule: integrate the signed work to the topmost + positive level, clamp at zero), max (maximum buoyancy work), and + reach (ballistic coasting from the LFC): the readout turns red when top + returns zero for a convectively viable profile. Parcels: A = ambient air lifted + as-is; B = boundary-layer air moisture-enriched by the storm's converged central + pressure; C = air saturated at the sea-surface temperature. VMAX² ∝ + (SST/TO)(EC − EB), so with B or C + selected the other of the pair is drawn thin and dashed: the PI winds live in + the gap between those two curves. Densities, class percentages, matching counts, + and the ΔVMAX scatters are computed over the full population + (every column); only the animated flip-book draws from a ~20k stratified + sample. Space bar toggles play. Click any marginal to filter the samples: + the map selects a position (±1,500 km), the rings a day-of-year window + (±15 days) or an hour-of-day window (±3 h UTC), the track a single year; each marginal keeps showing the + density under the other filters, with its own selection dashed. + Click the × chips or press Escape to clear; on the ΔVMAX panels, + clicking inside the current selection rectangle also clears it. The condition chips restrict + everything to profiles where the LNB conventions disagree: The OTL arrow on the + scope's pressure axis marks the storm outflow level — parcel C's LNB under + max, the one that sets TO in the PI formulas — for every profile, + whichever parcel is selected. T/M/R lanes + mark the selected parcel's LNB under top, max, and reach (the first-crossing + convention is omitted — it is materially wrong too often to be + interesting); CAPE and potential intensity per + convention appear in the readout (PI computed at tcpyPI's own ptop = 50). + The eight real columns where the top convention fails to converge are + included — select "top non-convergent" to see them. +
+ + + diff --git a/era5_analysis/validate_scan.py b/era5_analysis/validate_scan.py new file mode 100644 index 0000000..2e3c6bc --- /dev/null +++ b/era5_analysis/validate_scan.py @@ -0,0 +1,76 @@ +"""Cross-validate the scanner against the real cape() implementations: +E_top vs legacy (modernize_eg) cape(), E_max vs max-work (lnb-max-work) cape(). +""" + +import os as _os +from pathlib import Path as _Path +SCRATCH = _os.environ.get("ERA5_SCRATCH", str(_Path(__file__).resolve().parent / "work")) +_SRC = str(_Path(__file__).resolve().parents[1] / "src") + +import importlib.util +import sys + +import numpy as np + + +OLD_SCRATCH = _os.environ.get("LEGACY_TREE", SCRATCH) + + +def load_pkg(name, src_root): + spec = importlib.util.spec_from_file_location( + name, f"{src_root}/tcpyPI/__init__.py", + submodule_search_locations=[f"{src_root}/tcpyPI"]) + mod = importlib.util.module_from_spec(spec) + sys.modules[name] = mod + # make relative imports resolve under the alias + sys.path.insert(0, src_root) + spec.loader.exec_module(mod) + sys.path.pop(0) + return mod + + +d = np.load(f"{SCRATCH}/profiles_converted.npz") +res = np.load(f"{SCRATCH}/buoyancy_scan_results.npz") +sc = res["scalar"] +names = list(res["scalar_names"]) +ix = {n: i for i, n in enumerate(names)} + +P = d["P"] +rng = np.random.default_rng(7) +ok_rows = np.where(sc[:, ix["iflag"]] == 1)[0] +sel = rng.choice(ok_rows, size=400, replace=False) +# ensure multi-region rows are represented +multi = ok_rows[sc[ok_rows, ix["n"]] >= 4] +if len(multi): + sel = np.concatenate([sel, rng.choice(multi, size=min(200, len(multi)), replace=False)]) + +# import the NEW (max-work) cape from the lnb-max-work checkout +sys.path.insert(0, _SRC) +from tcpyPI.pi import cape as cape_new # noqa: E402 + +# legacy cape via the archived modernize_eg tree (pure-python instrumented copy) +sys.path.insert(0, OLD_SCRATCH) +from gmap import cape_instr as cape_legacy # noqa: E402 (pcmin convention) + +bad_top = bad_max = 0 +worst_top = worst_max = 0.0 +for i in sel: + T = d["TC"][i] + 273.15 + R = d["R"][i] * 0.001 + out_new = cape_new(T[0], R[0], P[0], T, R, P, 0, 50, 1) + o_leg = cape_legacy(T[0], R[0], P[0], T, R, P, 0, 50) + e_max_scan = sc[i, ix["E_max"]] + e_top_scan = sc[i, ix["E_top"]] + dmax = abs(out_new[0] - e_max_scan) + dtop = abs(o_leg["CAPED"] - e_top_scan) + worst_max = max(worst_max, dmax) + worst_top = max(worst_top, dtop) + if dmax > 1e-8 * max(1.0, e_max_scan): + bad_max += 1 + if dtop > 1e-8 * max(1.0, e_top_scan): + bad_top += 1 + +print(f"checked {len(sel)} profiles " + f"({int((sc[sel, ix['n']] >= 4).sum())} with n>=4 regions)") +print(f"E_max vs new cape(): mismatches {bad_max}, worst |d| = {worst_max:.3e}") +print(f"E_top vs legacy cape(): mismatches {bad_top}, worst |d| = {worst_top:.3e}")