From 84871fbe56d499c80d91933443bab9b2f6a7d18f Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Fri, 11 Sep 2026 09:00:46 -0500 Subject: [PATCH] FIX: Stop inflating pluvio gauge variables by 5x in the column match Ground instruments matched with resample="sum" were binned into a hard-coded 5-minute total and then linearly interpolated back onto the column's 1-minute time grid, so every output step received a whole 5-minute bin. Accumulations came out inflated by exactly the ratio of the two steps: a 19.29 mm day was written as 96.45 mm. The sum was also applied indiscriminately to every surviving pluvio variable, including two rain intensities (mm/hr) and the running bucket level, none of which may be summed. sgpradclssC1.c0 was reporting a mean intensity_rt of 139 mm/hr and a bucket_nrt of 3327 mm. - Integrate only true accumulations (accum_nrt, accum_rtnrt). Rates and the bucket level now take the same averaging path as every other instrument. - Regrid accumulations conservatively: integrate to a running total, interpolate that onto the column grid, then difference it back. This preserves the integral on regular and irregular (radar-based) time coordinates alike, and is an exact identity at native resolution. - Keep gaps missing. The old .sum() skipped NaNs with no min_count, so an all-NaN bin reported 0 mm and an outage read as "no rain". - Fix the closed="right" label. Pandas labels bin (t, t+step] as t, a timestamp the bin excludes, shifting data one step early. This also moves met, ldquants, vdisquants, wxt, kazr and sonde timing. - Only close datasets this function opened; _grd_raw was unbound whenever DataSet=True, breaking that match_datasets_act path. Existing c0 files need reprocessing, and the correction is not limited to accum_nrt. Co-Authored-By: Claude Opus 5 (1M context) --- src/radclss/util/column_utils.py | 102 +++++++++++++++++++++++++--- tests/test_column_utils.py | 112 ++++++++++++++++++++++++++++++- 2 files changed, 203 insertions(+), 11 deletions(-) diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index d16c0f9..fd8a39d 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -622,6 +622,74 @@ def subset_points( return ds +#: Variables holding precipitation accumulated over one sampling interval. +#: Only these are re-binned by integrating; every other variable on an +#: accumulating instrument is an instantaneous rate (mm/hr) or a running bucket +#: level, neither of which may be summed. +ACCUMULATION_VARS = frozenset({"accum_nrt", "accum_rtnrt"}) + + +def _column_time_step(column_time, default): + """ + Nominal spacing of a time coordinate, as a ``pd.Timedelta``. + + The median difference is used because radar-based time coordinates are + irregular. Falls back to ``default`` for coordinates too short or too + degenerate to measure. + """ + times = np.asarray(column_time.values).ravel() + if times.size < 2: + return pd.Timedelta(default) + deltas = np.diff(times).astype("timedelta64[ns]").astype(np.int64) + step = pd.Timedelta(int(np.median(deltas)), unit="ns") + if step <= pd.Timedelta(0): + return pd.Timedelta(default) + return step + + +def _accumulate_to_grid(grd_ds, column_time, default_step): + """ + Re-bin per-interval precipitation accumulations onto the column time grid. + + Accumulation variables hold the precipitation that fell during one sampling + interval, stamped at the interval end. Summing them into fixed-width bins + and interpolating the bin totals back onto a finer column grid hands every + output step a whole bin's total, inflating the accumulation by the ratio of + the two steps. Integrating to a running total instead, interpolating that + onto the column grid and differencing it back, yields the accumulation over + each output interval and preserves the integral on regular and irregular + grids alike. + """ + target = np.asarray(column_time.values).ravel() + step = _column_time_step(column_time, default_step) + edges = np.concatenate([[target[0] - step.to_timedelta64()], target]) + + # A zero anchor ahead of the record lets the first output interval + # difference against "nothing accumulated yet" rather than against a NaN. + source_step = _column_time_step(grd_ds["time"], default_step) + anchor = min(grd_ds["time"].values[0] - source_step.to_timedelta64(), edges[0]) + + def _anchored(ds): + head = (ds.isel(time=0) * 0).assign_coords(time=anchor).expand_dims("time") + return xr.concat([head, ds], dim="time").interp(time=edges, method="linear") + + def _running(ds): + # cumsum drops the dimension coordinate, so put it back. + return ds.cumsum("time").assign_coords(time=grd_ds["time"]) + + running = _anchored(_running(grd_ds.fillna(0))) + counted = _anchored(_running(grd_ds.notnull().astype("float64"))) + + # Intervals covered by no valid source sample stay missing instead of + # reporting the zero that a gap in the running total would imply. + matched = running.diff("time").where(counted.diff("time") > 0) + + for name in matched.data_vars: + matched[name].attrs.update(grd_ds[name].attrs) + matched.attrs.update(grd_ds.attrs) + return matched + + def _prepare_match( ground, site, @@ -680,20 +748,32 @@ def _prepare_match( ] grd_ds = grd_ds.drop_vars(non_numeric_vars) - if resample == "mean": - matched = ( - grd_ds.resample(time=resample_time, closed="right") + def _averaged(ds): + # closed="right" bins as (t, t + step]; pandas would otherwise label + # that bin t, a timestamp the bin excludes, shifting data one step early. + return ( + ds.resample(time=resample_time, closed="right", label="right") .mean(keep_attrs=True) .interp(time=column_time, method="linear") ) + + if resample == "mean": + matched = _averaged(grd_ds) elif resample == "skip": matched = grd_ds.interp(time=column_time, method="linear") elif resample == "sum": - matched = ( - grd_ds.resample(time=resample_time, closed="right") - .sum(keep_attrs=True) - .interp(time=column_time, method="linear") - ) + # Only true accumulations may be integrated. The remaining variables on + # a gauge are rain intensities and the running bucket level, which are + # averaged like any other instantaneous measurement. + accum = [v for v in grd_ds.data_vars if v in ACCUMULATION_VARS] + instant = [v for v in grd_ds.data_vars if v not in ACCUMULATION_VARS] + parts = [] + if accum: + parts.append(_accumulate_to_grid(grd_ds[accum], column_time, resample_time)) + if instant: + parts.append(_averaged(grd_ds[instant])) + matched = xr.merge(parts, combine_attrs="override") + matched.attrs.update(grd_ds.attrs) else: raise ValueError( "Invalid resample method. Please choose 'mean', 'sum', or 'skip'." @@ -708,8 +788,10 @@ def _prepare_match( for var in matched.data_vars: matched[var].attrs.update(source=matched.datastream) - grd_ds.close() - _grd_raw.close() + # Only close what this function opened; a caller-supplied dataset is theirs. + if not DataSet: + grd_ds.close() + _grd_raw.close() return site, matched diff --git a/tests/test_column_utils.py b/tests/test_column_utils.py index 44dc56e..0b56d10 100644 --- a/tests/test_column_utils.py +++ b/tests/test_column_utils.py @@ -1,9 +1,15 @@ from unittest.mock import MagicMock, patch import numpy as np +import pandas as pd +import pytest import xarray as xr -from radclss.util.column_utils import get_nexrad_column +from radclss.util.column_utils import ( + _accumulate_to_grid, + _column_time_step, + get_nexrad_column, +) def test_get_nexrad_column(): @@ -157,3 +163,107 @@ def test_get_nexrad_column_integration(): assert "height" in result.dims assert result.dims["station"] == len(input_site_dict) assert "reflectivity" in result.data_vars + + +def _synthetic_gauge(n_minutes=120): + """A 1-minute gauge record with a rain burst, a dry spell and a data gap.""" + time = pd.date_range("2025-06-19T00:01", periods=n_minutes, freq="1min") + accum = np.zeros(n_minutes) + accum[20:60] = np.linspace(0.05, 0.9, 40) # burst + accum[70:80] = 0.2 # second, flatter burst + accum[90:100] = np.nan # instrument gap + return xr.Dataset( + {"accum_nrt": ("time", accum, {"units": "mm", "long_name": "Accum"})}, + coords={"time": time}, + ) + + +@pytest.mark.parametrize("freq", ["1min", "5min", "15min"]) +def test_accumulate_to_grid_conserves_total_on_regular_grids(freq): + """ + Re-binning an accumulation must move rain between time steps, never create + it. Summing into fixed bins and interpolating back onto a finer grid used to + inflate the daily total by the ratio of the two steps. + """ + gauge = _synthetic_gauge() + column_time = xr.DataArray( + pd.date_range(gauge.time.values[0], gauge.time.values[-1], freq=freq), + dims="time", + name="time", + ) + + regridded = _accumulate_to_grid(gauge, column_time, "5Min") + + assert np.isclose( + np.nansum(regridded["accum_nrt"].values), + np.nansum(gauge["accum_nrt"].values), + ) + + +def test_accumulate_to_grid_is_identity_at_native_resolution(): + """On a grid matching the source, the values must come back untouched.""" + gauge = _synthetic_gauge() + + regridded = _accumulate_to_grid(gauge, gauge.time, "5Min") + + np.testing.assert_allclose( + np.nan_to_num(regridded["accum_nrt"].values), + np.nan_to_num(gauge["accum_nrt"].values), + atol=1e-9, + ) + assert regridded["accum_nrt"].attrs["units"] == "mm" + + +def test_accumulate_to_grid_conserves_total_on_irregular_grid(): + """Radar-based time coordinates are irregular; the integral must survive.""" + gauge = _synthetic_gauge() + rng = np.random.default_rng(0) + picks = np.sort(rng.choice(np.arange(1, gauge.sizes["time"]), 40, replace=False)) + column_time = gauge.time.isel(time=picks) + + regridded = _accumulate_to_grid(gauge, column_time, "5Min") + + assert np.isclose( + np.nansum(regridded["accum_nrt"].values), + np.nansum(gauge["accum_nrt"].values), + ) + + +def test_accumulate_to_grid_keeps_gaps_missing(): + """A stretch with no valid samples must stay missing, not report zero rain.""" + gauge = _synthetic_gauge() + column_time = xr.DataArray( + pd.date_range(gauge.time.values[0], gauge.time.values[-1], freq="5min"), + dims="time", + name="time", + ) + + regridded = _accumulate_to_grid(gauge, column_time, "5Min") + + assert np.isnan(regridded["accum_nrt"].values).any() + + +def test_column_time_step_falls_back_when_unmeasurable(): + """Degenerate grids fall back to the supplied default rather than raising.""" + single = xr.DataArray(pd.to_datetime(["2025-06-19T00:00"]), dims="time") + duplicated = xr.DataArray(pd.to_datetime(["2025-06-19T00:00"] * 4), dims="time") + + assert _column_time_step(single, "5Min") == pd.Timedelta("5min") + assert _column_time_step(duplicated, "5Min") == pd.Timedelta("5min") + + +def test_column_time_step_measures_regular_and_irregular_grids(): + regular = xr.DataArray( + pd.date_range("2025-06-19", periods=10, freq="1min"), dims="time" + ) + jittered = xr.DataArray( + pd.to_datetime( + ["2025-06-19T00:00", "2025-06-19T00:05", "2025-06-19T00:11"] + + ["2025-06-19T00:16", "2025-06-19T00:21", "2025-06-19T01:30"] + ), + dims="time", + ) + + assert _column_time_step(regular, "5Min") == pd.Timedelta("1min") + # Median, so the one long outlying gap does not set the step. + assert _column_time_step(jittered, "5Min") == pd.Timedelta("5min")