From 101246c0f9a07f63e4c8753950080e24e008f08f Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Thu, 10 Sep 2026 16:31:28 -0500 Subject: [PATCH] FIX: Handle midnight wrap-around in the column plotting suite A RadCLss file holds the radar volumes for a single processing day, but the first volume of that day routinely starts a few minutes before midnight, so time[0] falls on the previous calendar day. Both plotting functions anchored their 24 hour window on time[0] floored to midnight, which selected the wrong day entirely. In create_radclss_columns the window is a .sel(time=slice(...)), so the data was discarded rather than merely clipped: nsaradclssC1.c0.20260713 .235803.nc plotted 2 of its 1440 times. 5 of the 12 RadCLss files on hand were affected, at both NSA and BNF. It went unnoticed because the test fixture starts exactly at 00:00:00. Replace the anchoring with _daily_time_window(), which takes the day from the timestamp the bulk of the samples fall on and then widens the window to keep the pre-midnight leader and any post-midnight tail, so no data is dropped. The derived day agrees with the authoritative -b processing day recorded in each file's command_line attribute for all 12 files tested. Also give the reflectivity panel of the rainfall timeseries the same xlim as the rain rate and accumulation panels, which previously had none and so never aligned, and report the product day in its suptitle. Adds six regression tests built on a synthetic wrap-around dataset that mirrors the NSA file; the two end-to-end tests were confirmed to fail against the pre-fix code. Co-Authored-By: Claude Opus 5 (1M context) --- pyproject.toml | 2 +- src/radclss/vis/quicklooks.py | 82 +++++++++++++-------- tests/test_vis.py | 135 ++++++++++++++++++++++++++++++++++ 3 files changed, 187 insertions(+), 32 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index c099636..57c0dc0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "radclss" -version = "2026.7.22" +version = "2026.9.10" description = "Extracted Radar Columns and In Situ Sensors" readme = "README.md" requires-python = ">=3.10" diff --git a/src/radclss/vis/quicklooks.py b/src/radclss/vis/quicklooks.py index 7472b09..13b29eb 100644 --- a/src/radclss/vis/quicklooks.py +++ b/src/radclss/vis/quicklooks.py @@ -1,6 +1,4 @@ -import datetime import sys -from datetime import timedelta import act import matplotlib.pyplot as plt @@ -9,6 +7,46 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable +def _daily_time_window(times): + """ + Determine the 24 hour window a RadCLss file represents. + + A RadCLss file holds the radar volumes belonging to a single processing + day, but the first volume of that day routinely *starts* a few minutes + before midnight (the NSA XSAPR volume beginning 23:58:03 UTC belongs to + the following day's file). Anchoring on ``time[0]`` therefore selects + the previous calendar day and discards nearly every sample, so the day is + taken from the timestamp the bulk of the samples fall on and then widened + to keep any pre-midnight leader and post-midnight tail. + + Input + ----- + times : array-like of numpy.datetime64 + Time coordinate of the RadCLss dataset. + + Output + ------ + start : numpy.datetime64 + Start of the plotting window. + end : numpy.datetime64 + End of the plotting window. + ref_day : numpy.datetime64 + Midnight of the day the file represents, as a ``datetime64[D]``. + + """ + times = np.asarray(times, dtype="datetime64[s]") + if times.size == 0: + raise ValueError( + "\nERROR - (_daily_time_window):" + + " \n\tRadCLss dataset has an empty time coordinate.\n" + ) + days, counts = np.unique(times.astype("datetime64[D]"), return_counts=True) + ref_day = days[np.argmax(counts)] + start = min(ref_day.astype("datetime64[s]"), times.min()) + end = max((ref_day + 1).astype("datetime64[s]"), times.max()) + return start, end, ref_day + + def create_radclss_columns( radclss, field="corrected_reflectivity", @@ -85,22 +123,16 @@ def create_radclss_columns( fig, axarr = plt.subplots(nrows, ncols, figsize=(width, height)) plt.subplots_adjust(hspace=0.8) - # Define the time of the radar file we are plotting against - radar_time = datetime.datetime.strptime( - np.datetime_as_string(ds["time"].data[0], unit="s"), "%Y-%m-%dT%H:%M:%S" - ).replace(tzinfo=datetime.timezone.utc) - final_time = radar_time + timedelta(days=1) + # Define the window of the radar file we are plotting against + start_time, end_time, _ = _daily_time_window(ds["time"].data) for i, station in enumerate(stations): row = i // 2 col = i % 2 if len(axarr.shape) == 1: axarr = np.expand_dims(axarr, axis=0) - ds[field].sel(station=station).sel( - time=slice( - radar_time.strftime("%Y-%m-%dT00:00:00"), - final_time.strftime("%Y-%m-%dT00:00:00"), - ) - ).plot(y="height", ax=axarr[row, col], vmin=vmin, vmax=vmax, **kwargs) + ds[field].sel(station=station).sel(time=slice(start_time, end_time)).plot( + y="height", ax=axarr[row, col], vmin=vmin, vmax=vmax, **kwargs + ) long_name = ds[field].attrs.get("long_name", field) axarr[row, col].set_title(f"{station} {long_name}") @@ -193,11 +225,8 @@ def create_radclss_rainfall_timeseries( print("\n") return - # Define the time of the radar file we are plotting against - radar_time = datetime.datetime.strptime( - np.datetime_as_string(ds["time"].data[0], unit="s"), "%Y-%m-%dT%H:%M:%S" - ).replace(tzinfo=datetime.timezone.utc) - final_time = radar_time + timedelta(days=1) + # Define the window of the radar file we are plotting against + start_time, end_time, ref_day = _daily_time_window(ds["time"].data) # ----------------------------------------------- # Side Plot A - Display the RadClss Radar Field @@ -214,6 +243,7 @@ def create_radclss_rainfall_timeseries( ) ax2.set_ylabel("Height [m]") ax2.set_xlabel("Time [UTC]") + ax2.set_xlim([start_time, end_time]) # -------------------------------------- # Side Plot B - Display the Rain Rates @@ -237,12 +267,7 @@ def create_radclss_rainfall_timeseries( ax3.set_title(" ") ax3.set_ylabel("Precipitation Rate \n[mm/hr]") ax3.set_xlabel("Time [UTC]") - ax3.set_xlim( - [ - radar_time.strftime("%Y-%m-%dT00:00:00"), - final_time.strftime("%Y-%m-%dT00:00:00"), - ] - ) + ax3.set_xlim([start_time, end_time]) ax3.legend(loc="upper right") ax3.grid(True) ax3.set_ylim(rr_min, rr_max) @@ -282,12 +307,7 @@ def create_radclss_rainfall_timeseries( ax4.set_xlabel("Time [UTC]") ax4.legend(loc="upper left") ax4.grid(True) - ax4.set_xlim( - [ - radar_time.strftime("%Y-%m-%dT00:00:00"), - final_time.strftime("%Y-%m-%dT00:00:00"), - ] - ) + ax4.set_xlim([start_time, end_time]) ax4.set_ylim(cum_min, cum_max) # Add a blank space next to the subplot to shape it as the above plot divider = make_axes_locatable(ax4) @@ -302,7 +322,7 @@ def create_radclss_rainfall_timeseries( if title_flag is True: plt.suptitle( "BNF Extracted Radar Columns and In-Situ Sensors (RadCLss) \n" - + radar_time.strftime("%Y-%m-%d") + + str(ref_day) ) # Clean up this function diff --git a/tests/test_vis.py b/tests/test_vis.py index b05768b..0ad595b 100644 --- a/tests/test_vis.py +++ b/tests/test_vis.py @@ -1,8 +1,14 @@ +import datetime + import arm_test_data +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import numpy as np import pytest import xarray as xr import radclss +from radclss.vis import quicklooks @pytest.mark.mpl_image_compare @@ -156,3 +162,132 @@ def test_create_radclss_timeseries_no_title(): assert fig is not None assert axarr is not None return fig + + +def _wraparound_dataset(start="2026-07-13T23:58:00", periods=1440): + """ + Build a RadCLss-style dataset whose first volume starts before midnight. + + Mirrors ``nsaradclssC1.c0.20260713.235803.nc``: the file is the 2026-07-14 + product, but because the first radar volume of the day starts at 23:58 the + first timestamp falls on 2026-07-13. + """ + time = np.datetime64(start) + np.arange(periods) * np.timedelta64(1, "m") + height = np.arange(0.0, 3000.0, 500.0) + station = ["M1", "S30"] + + ramp = np.linspace(0.0, 1.0, time.size) + field = ( + ramp[:, None, None] + * np.ones((1, height.size, 1)) + * np.ones((1, 1, len(station))) + ) * 40.0 - 10.0 + rain = np.tile(ramp[:, None], (1, len(station))) + + return xr.Dataset( + { + "corrected_reflectivity": ( + ("time", "height", "station"), + field, + {"long_name": "Equivalent reflectivity factor", "units": "dBZ"}, + ), + "rain_rate_A": ( + ("time", "height", "station"), + field * 0.0 + 1.0, + {"long_name": "Rain rate", "units": "mm/hr"}, + ), + "intensity_rtnrt": ( + ("time", "station"), + rain, + {"long_name": "Pluvio2 rain rate", "units": "mm/hr"}, + ), + "ldquants_rain_rate": ( + ("time", "station"), + rain, + {"long_name": "LDQUANTS rain rate", "units": "mm/hr"}, + ), + }, + coords={"time": time, "height": height, "station": station}, + ) + + +def test_daily_time_window_wraparound(): + """The day is taken from the bulk of the samples, not from time[0].""" + ds = _wraparound_dataset() + start, end, ref_day = quicklooks._daily_time_window(ds["time"].data) + + # time[0] is on 07-13, but the file is the 07-14 product + assert ref_day == np.datetime64("2026-07-14", "D") + # the pre-midnight leader and the last sample are both kept + assert start == np.datetime64("2026-07-13T23:58:00", "s") + assert end == np.datetime64("2026-07-15T00:00:00", "s") + + +def test_daily_time_window_no_wraparound(): + """A file starting exactly at midnight is anchored on its own day.""" + ds = _wraparound_dataset(start="2026-07-14T00:00:00") + start, end, ref_day = quicklooks._daily_time_window(ds["time"].data) + + assert ref_day == np.datetime64("2026-07-14", "D") + assert start == np.datetime64("2026-07-14T00:00:00", "s") + assert end == np.datetime64("2026-07-15T00:00:00", "s") + + +def test_daily_time_window_partial_day(): + """A short file that wraps is still anchored on the majority day.""" + ds = _wraparound_dataset(start="2026-07-17T23:54:00", periods=184) + _start, _end, ref_day = quicklooks._daily_time_window(ds["time"].data) + + assert ref_day == np.datetime64("2026-07-18", "D") + + +def test_daily_time_window_empty(): + with pytest.raises(ValueError, match="empty time coordinate"): + quicklooks._daily_time_window(np.array([], dtype="datetime64[s]")) + + +def test_create_radclss_columns_wraparound_keeps_all_times(): + """Columns spanning midnight must not be discarded by the day slice.""" + ds = _wraparound_dataset() + fig, axarr = radclss.vis.create_radclss_columns( + ds, field="corrected_reflectivity", vmin=-10, vmax=30 + ) + + # every sample in the file is drawn, not just the two before midnight + mesh = axarr[0, 0].collections[0] + n_plotted = mesh.get_coordinates().shape[1] - 1 + assert n_plotted == ds.sizes["time"] + + # and the axis covers the 07-14 product day rather than 07-13 + x_min, x_max = axarr[0, 0].get_xlim() + span = mdates.num2date(x_max) - mdates.num2date(x_min) + assert span > datetime.timedelta(hours=23) + assert ( + mdates.num2date(x_min) + < datetime.datetime(2026, 7, 14, 12, tzinfo=datetime.timezone.utc) + < mdates.num2date(x_max) + ) + plt.close(fig) + + +def test_create_radclss_timeseries_wraparound_panels_aligned(): + """All three timeseries panels share the corrected day window.""" + ds = _wraparound_dataset() + fig, axarr = radclss.vis.create_radclss_rainfall_timeseries( + ds, field="corrected_reflectivity", rheight=1000 + ) + + limits = [ax.get_xlim() for ax in axarr] + assert limits[0] == limits[1] == limits[2] + + x_min, x_max = limits[0] + assert mdates.num2date(x_min) == datetime.datetime( + 2026, 7, 13, 23, 58, tzinfo=datetime.timezone.utc + ) + assert mdates.num2date(x_max) == datetime.datetime( + 2026, 7, 15, tzinfo=datetime.timezone.utc + ) + + # the suptitle reports the product day, not the day of time[0] + assert fig._suptitle.get_text().endswith("2026-07-14") + plt.close(fig)