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)