diff --git a/pyproject.toml b/pyproject.toml index 2b6bbb6..e0a8be4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "radclss" -version = "2026.9.11" +version = "2026.9.12" description = "Extracted Radar Columns and In Situ Sensors" readme = "README.md" requires-python = ">=3.10" diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index fd8a39d..f6e5f27 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -574,21 +574,28 @@ def subset_points( valid = np.isfinite(da["height"]) n_valid = int(valid.sum()) interpolated = False - dvars = da.data_vars + # time_offset is bookkeeping, not a measurement: gate_time + # below is computed from it. A column outside the radar's + # coverage has no gate times at all, so purging it along with + # the all-NaN fields would leave nothing for gate_time to read, + # and letting it drive dropna would discard every height. + dvars = [v for v in da.data_vars if v != "time_offset"] for v in dvars: if np.all(np.isnan(da[v].values)): - da = da.drop(v) + da = da.drop_vars(v) + measured = [v for v in da.data_vars if v != "time_offset"] if n_valid > 0: - da_clean = da.dropna("height").sortby("height") + da_clean = da.dropna("height", subset=measured).sortby("height") if da_clean.sizes.get("height", 0) > 0: try: da = da_clean.interp(height=height_bins) except pd.errors.InvalidIndexError: da_clean = da_clean.drop_duplicates("height", keep="first") da = da_clean.interp(height=height_bins) - time_offset = time_offset.drop_duplicates( - "height", keep="first" - ) + if "rhi" in radar.scan_type: + time_offset = time_offset.drop_duplicates( + "height", keep="first" + ) interpolated = True if not interpolated: diff --git a/tests/test_column_utils.py b/tests/test_column_utils.py index 0b56d10..114bd0e 100644 --- a/tests/test_column_utils.py +++ b/tests/test_column_utils.py @@ -9,6 +9,7 @@ _accumulate_to_grid, _column_time_step, get_nexrad_column, + subset_points, ) @@ -267,3 +268,107 @@ def test_column_time_step_measures_regular_and_irregular_grids(): 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") + + +def _out_of_coverage_column(heights): + """ + What Py-ART hands back for a site the radar never sampled. + + ``column_vertical_profile`` averages an empty set of rays into NaN for + every field, and ``time_offset`` -- the mean gate time -- comes back NaN + right along with them. + """ + n = len(heights) + return xr.Dataset( + { + "reflectivity": (["height"], np.full(n, np.nan)), + "time_offset": (["height"], np.full(n, np.nan)), + "base_time": np.datetime64("2025-06-19T00:00:00"), + }, + coords={"height": (["height"], heights)}, + ).set_coords(["base_time"]) + + +def test_subset_points_survives_a_column_outside_radar_coverage(): + """ + A site the radar never sampled yields an all-NaN column, ``time_offset`` + included. Purging every all-NaN variable used to take ``time_offset`` with + it and leave the gate_time calculation reading a variable that no longer + existed, raising ``KeyError: 'time_offset'`` and losing the whole scan -- + every co-located site with it, not just the one out of coverage. + """ + input_site_dict = {"M1": (34.34525, -87.33842, 293)} + height_bins = np.arange(500, 8500, 250) + + mock_radar = MagicMock() + mock_radar.scan_type = "ppi" + mock_radar.metadata = {"scan_mode": "ppi", "facility_id": "M1"} + mock_radar.time = {"data": np.arange(10.0)} + mock_radar.sweep_start_ray_index = {"data": np.ma.array([0])} + mock_radar.sweep_end_ray_index = {"data": np.ma.array([9])} + + column = _out_of_coverage_column(np.arange(500.0, 8500.0, 100.0)) + + with ( + patch("radclss.util.column_utils.pyart.io.read", return_value=mock_radar), + patch( + "radclss.util.column_utils.pyart.util.columnsect.column_vertical_profile", + return_value=column, + ), + ): + result = subset_points( + "bnfcsapr2cfrS3.a1.20250619.000000.nc", + input_site_dict, + height_bins=height_bins, + ) + + assert result is not None + assert "time_offset" in result + assert result.sizes["station"] == 1 + np.testing.assert_array_equal(result["height"].values, height_bins) + # No gate times to report, but the column still comes back and says so. + assert np.all(np.isnat(result["time_offset"].values)) + assert np.all(np.isnat(result["gate_time"].values)) + + +def test_subset_points_keeps_heights_a_partly_missing_field_would_drop(): + """ + ``time_offset`` must not decide which heights survive ``dropna``. Gates the + radar did sample stay in the column even where the gate time is missing. + """ + input_site_dict = {"M1": (34.34525, -87.33842, 293)} + height_bins = np.arange(500, 2500, 250) + heights = np.arange(500.0, 2500.0, 100.0) + + offsets = np.zeros(len(heights)) + offsets[::2] = np.nan # gate times missing on half the gates + column = xr.Dataset( + { + "reflectivity": (["height"], np.linspace(10.0, 40.0, len(heights))), + "time_offset": (["height"], offsets), + "base_time": np.datetime64("2025-06-19T00:00:00"), + }, + coords={"height": (["height"], heights)}, + ).set_coords(["base_time"]) + + mock_radar = MagicMock() + mock_radar.scan_type = "ppi" + mock_radar.metadata = {"scan_mode": "ppi", "facility_id": "M1"} + mock_radar.time = {"data": np.arange(10.0)} + mock_radar.sweep_start_ray_index = {"data": np.ma.array([0])} + mock_radar.sweep_end_ray_index = {"data": np.ma.array([9])} + + with ( + patch("radclss.util.column_utils.pyart.io.read", return_value=mock_radar), + patch( + "radclss.util.column_utils.pyart.util.columnsect.column_vertical_profile", + return_value=column, + ), + ): + result = subset_points( + "bnfcsapr2cfrS3.a1.20250619.000000.nc", + input_site_dict, + height_bins=height_bins, + ) + + assert np.isfinite(result["reflectivity"].values).all()