Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
82 changes: 51 additions & 31 deletions src/radclss/vis/quicklooks.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,4 @@
import datetime
import sys
from datetime import timedelta

import act
import matplotlib.pyplot as plt
Expand All @@ -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",
Expand Down Expand Up @@ -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}")

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
135 changes: 135 additions & 0 deletions tests/test_vis.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
Loading