FIX: Stop inflating pluvio gauge variables by 5x in the column match

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
 
 

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")