From d9c1aab9a6faf2d5fa32d25ce4f8a230dd154678 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 4 May 2026 15:27:29 -0500 Subject: [PATCH 1/8] ENH: VPT support in subset_points - return ray time series for co-located site, NaNs otherwise Co-Authored-By: Claude Sonnet 4.6 --- src/radclss/util/column_utils.py | 141 +++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index b8488e6..57bdcdb 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -126,6 +126,137 @@ def _grab_90_degree_rays(radar): return column +def _vpt_to_column_timeseries(radar, height_bins): + """Convert all rays of a VPT scan to a time series of columns. + + For VPT scans elevation ≈ 90°, so gate range ≈ height above radar. + Each ray becomes one time step; all rays are stacked into a (time, height) dataset. + """ + zgate = radar.range["data"] + radar.altitude["data"][0] + + base_vol_time = np.datetime64( + pyart.util.datetime_from_radar(radar).isoformat(), "ns" + ) + ray_time_offsets = pd.to_timedelta(radar.time["data"], unit="s") + abs_times = (base_vol_time + ray_time_offsets).astype("datetime64[ns]") + + da_meta = [ + "units", + "standard_name", + "long_name", + "valid_max", + "valid_min", + "coordinates", + ] + data_vars = {} + for key in radar.fields: + arr = np.ma.filled(radar.fields[key]["data"], np.nan).astype(float) + attrs = { + tag: radar.fields[key][tag] for tag in da_meta if tag in radar.fields[key] + } + data_vars[key] = xr.DataArray(arr, dims=["time", "height"], attrs=attrs) + + ds = xr.Dataset(data_vars, coords={"height": zgate, "time": abs_times}) + + valid_h = np.isfinite(ds["height"]) + if int(valid_h.sum()) > 0: + try: + ds = ds.dropna("height").sortby("height").interp(height=height_bins) + except pd.errors.InvalidIndexError: + ds = ( + ds.drop_duplicates("height", keep="first") + .dropna("height") + .sortby("height") + .interp(height=height_bins) + ) + else: + ds = ds.reindex(height=xr.DataArray(height_bins, dims="height", name="height")) + + abs_times_s = abs_times.astype("datetime64[s]") + ds["base_time"] = xr.DataArray(abs_times_s, dims="time") + ds["time_offset"] = xr.DataArray( + ray_time_offsets.values.astype("timedelta64[s]"), dims="time" + ) + ds["gate_time"] = xr.DataArray(abs_times_s, dims="time") + + height_des = ( + "Height Above Sea Level [in meters] for the Center of Each" + + " Radar Gate Above the Target Location" + ) + ds.height.attrs.update( + long_name="Height of Radar Beam", + units="m", + standard_name="height", + description=height_des, + ) + ds["base_time"].attrs.update(long_name="UTC Reference Time", units="seconds") + ds["time_offset"].attrs.update( + long_name="Time in Seconds Since Volume Start", units="seconds" + ) + ds.attrs["distance_from_radar"] = "0 km" + ds.attrs["latitude_of_location"] = str(radar.latitude["data"][0]) + " degrees" + ds.attrs["longitude_of_location"] = str(radar.longitude["data"][0]) + " degrees" + return ds + + +def _vpt_nan_fill(radar, height_bins): + """Return a NaN-filled dataset matching the VPT column time-series shape. + + Used for stations that are not co-located with the VPT radar site so that + the time axis still aligns with the radar-site columns. + """ + base_vol_time = np.datetime64( + pyart.util.datetime_from_radar(radar).isoformat(), "ns" + ) + ray_time_offsets = pd.to_timedelta(radar.time["data"], unit="s") + abs_times = (base_vol_time + ray_time_offsets).astype("datetime64[ns]") + abs_times_s = abs_times.astype("datetime64[s]") + + nrays = radar.nrays + n_heights = len(height_bins) + da_meta = [ + "units", + "standard_name", + "long_name", + "valid_max", + "valid_min", + "coordinates", + ] + data_vars = {} + for key in radar.fields: + attrs = { + tag: radar.fields[key][tag] for tag in da_meta if tag in radar.fields[key] + } + data_vars[key] = xr.DataArray( + np.full((nrays, n_heights), np.nan), + dims=["time", "height"], + attrs=attrs, + ) + + ds = xr.Dataset(data_vars, coords={"height": height_bins, "time": abs_times}) + ds["base_time"] = xr.DataArray(abs_times_s, dims="time") + ds["time_offset"] = xr.DataArray( + ray_time_offsets.values.astype("timedelta64[s]"), dims="time" + ) + ds["gate_time"] = xr.DataArray(abs_times_s, dims="time") + + height_des = ( + "Height Above Sea Level [in meters] for the Center of Each" + + " Radar Gate Above the Target Location" + ) + ds.height.attrs.update( + long_name="Height of Radar Beam", + units="m", + standard_name="height", + description=height_des, + ) + ds["base_time"].attrs.update(long_name="UTC Reference Time", units="seconds") + ds["time_offset"].attrs.update( + long_name="Time in Seconds Since Volume Start", units="seconds" + ) + return ds + + def get_nexrad_column( rad_time, site, @@ -373,6 +504,16 @@ def subset_points( # Make sure we are interpolating from the radar's location above sea level # NOTE: interpolating throughout Troposphere to match sonde to in the future + if "vpt" in radar.scan_type: + if radar.metadata.get("facility_id", "") == site: + da = _vpt_to_column_timeseries(radar, height_bins) + else: + da = _vpt_nan_fill(radar, height_bins) + da["lat"] = lat + da["lon"] = lon + column_list.append(da) + continue + if "rhi" not in radar.scan_type: da = pyart.util.columnsect.column_vertical_profile(radar, lat, lon) else: From 44469355d1d0f76219ce8bec52ee42ffc759b627 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 4 May 2026 15:56:37 -0500 Subject: [PATCH 2/8] FIX: Drop duplicate time indices in VPT column timeseries and NaN fill Co-Authored-By: Claude Sonnet 4.6 --- src/radclss/util/column_utils.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index 57bdcdb..10d794a 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -157,6 +157,7 @@ def _vpt_to_column_timeseries(radar, height_bins): data_vars[key] = xr.DataArray(arr, dims=["time", "height"], attrs=attrs) ds = xr.Dataset(data_vars, coords={"height": zgate, "time": abs_times}) + ds = ds.drop_duplicates("time", keep="first") valid_h = np.isfinite(ds["height"]) if int(valid_h.sum()) > 0: @@ -172,12 +173,11 @@ def _vpt_to_column_timeseries(radar, height_bins): else: ds = ds.reindex(height=xr.DataArray(height_bins, dims="height", name="height")) - abs_times_s = abs_times.astype("datetime64[s]") - ds["base_time"] = xr.DataArray(abs_times_s, dims="time") - ds["time_offset"] = xr.DataArray( - ray_time_offsets.values.astype("timedelta64[s]"), dims="time" - ) - ds["gate_time"] = xr.DataArray(abs_times_s, dims="time") + dedup_times_s = ds["time"].values.astype("datetime64[s]") + dedup_offsets = (ds["time"].values - base_vol_time).astype("timedelta64[s]") + ds["base_time"] = xr.DataArray(dedup_times_s, dims="time") + ds["time_offset"] = xr.DataArray(dedup_offsets, dims="time") + ds["gate_time"] = xr.DataArray(dedup_times_s, dims="time") height_des = ( "Height Above Sea Level [in meters] for the Center of Each" @@ -210,7 +210,6 @@ def _vpt_nan_fill(radar, height_bins): ) ray_time_offsets = pd.to_timedelta(radar.time["data"], unit="s") abs_times = (base_vol_time + ray_time_offsets).astype("datetime64[ns]") - abs_times_s = abs_times.astype("datetime64[s]") nrays = radar.nrays n_heights = len(height_bins) @@ -234,11 +233,13 @@ def _vpt_nan_fill(radar, height_bins): ) ds = xr.Dataset(data_vars, coords={"height": height_bins, "time": abs_times}) - ds["base_time"] = xr.DataArray(abs_times_s, dims="time") - ds["time_offset"] = xr.DataArray( - ray_time_offsets.values.astype("timedelta64[s]"), dims="time" - ) - ds["gate_time"] = xr.DataArray(abs_times_s, dims="time") + ds = ds.drop_duplicates("time", keep="first") + + dedup_times_s = ds["time"].values.astype("datetime64[s]") + dedup_offsets = (ds["time"].values - base_vol_time).astype("timedelta64[s]") + ds["base_time"] = xr.DataArray(dedup_times_s, dims="time") + ds["time_offset"] = xr.DataArray(dedup_offsets, dims="time") + ds["gate_time"] = xr.DataArray(dedup_times_s, dims="time") height_des = ( "Height Above Sea Level [in meters] for the Center of Each" From 2f60141f2c399d28969696166de5a9df3db2b955 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 4 May 2026 20:58:32 +0000 Subject: [PATCH 3/8] FIX: Correct attribute --- src/radclss/util/column_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index 10d794a..82d611b 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -505,7 +505,7 @@ def subset_points( # Make sure we are interpolating from the radar's location above sea level # NOTE: interpolating throughout Troposphere to match sonde to in the future - if "vpt" in radar.scan_type: + if "vpt" in radar.metadata["scan_mode"]: if radar.metadata.get("facility_id", "") == site: da = _vpt_to_column_timeseries(radar, height_bins) else: From 4135f29d5cfcc4dad3096842a67eec74b31b3bf9 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 4 May 2026 16:02:47 -0500 Subject: [PATCH 4/8] FIX: Use .flat[0] for base_time to handle variable-length VPT time arrays Co-Authored-By: Claude Sonnet 4.6 --- src/radclss/core/radclss_core.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/radclss/core/radclss_core.py b/src/radclss/core/radclss_core.py index 7c70b46..3f8b057 100644 --- a/src/radclss/core/radclss_core.py +++ b/src/radclss/core/radclss_core.py @@ -232,7 +232,9 @@ def radclss( max_times = {} for k in columns.keys(): if "radar" in k and len(columns[k]) > 0: - times = np.array([x["base_time"].values[0] for x in columns[k]]) + times = np.array( + [x["base_time"].values.flat[0] for x in columns[k] if x is not None] + ) min_times[k] = np.min(times) max_times[k] = np.max(times) if verbose: From bd0896df45f2d7f0160f6b93f5ec3fe5caa515e4 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 4 May 2026 16:57:51 -0500 Subject: [PATCH 5/8] FIX: Normalize time dimension before concat to handle mixed VPT/non-VPT files Co-Authored-By: Claude Sonnet 4.6 --- src/radclss/core/radclss_core.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/radclss/core/radclss_core.py b/src/radclss/core/radclss_core.py index 3f8b057..4a5d17e 100644 --- a/src/radclss/core/radclss_core.py +++ b/src/radclss/core/radclss_core.py @@ -386,9 +386,18 @@ def _get_nexrad_wrapper(time_str): if verbose: print(f" Processing {k}...") valid = [data for data in columns[k] if data is not None] + # Non-VPT datasets have no 'time' dimension; VPT datasets do. + # Expand non-VPT datasets so all entries have a consistent 'time' + # coordinate before concatenation. + expanded = [] + for ds in valid: + if "time" not in ds.dims: + bt = np.atleast_1d(ds["base_time"].values.flat[0]) + ds = ds.expand_dims("time").assign_coords(time=bt) + expanded.append(ds) chunks = [ - xr.concat(valid[i : i + _CONCAT_CHUNK], dim="time") - for i in range(0, len(valid), _CONCAT_CHUNK) + xr.concat(expanded[i : i + _CONCAT_CHUNK], dim="time") + for i in range(0, len(expanded), _CONCAT_CHUNK) ] ds_concat[k] = xr.concat(chunks, dim="time") if len(chunks) > 1 else chunks[0] if verbose: From a33ed9b4f671046edd9806fd10dd33ebb9575158 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Tue, 5 May 2026 17:19:36 +0000 Subject: [PATCH 6/8] FIX: Height bins for nexrad --- src/radclss/core/radclss_core.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/radclss/core/radclss_core.py b/src/radclss/core/radclss_core.py index 4a5d17e..808eaf0 100644 --- a/src/radclss/core/radclss_core.py +++ b/src/radclss/core/radclss_core.py @@ -303,6 +303,7 @@ def _get_nexrad_wrapper(time_str): output_config["site"], input_site_dict, nexrad_radar=nexrad_site, + height_bins=height_bins ) results = current_client.map(_get_nexrad_wrapper, time_list) @@ -339,6 +340,7 @@ def _get_nexrad_wrapper(time_str): output_config["site"], input_site_dict, nexrad_radar=nexrad_site, + height_bins=height_bins ) successful_count = 0 @@ -434,6 +436,7 @@ def _get_nexrad_wrapper(time_str): print("=" * 80) print(f" Time coordinate method: {time_coords}") + ds_concat[k] = ds_concat[k].drop_duplicates('time') if "radar" in time_coords: if verbose: print(f" Reindexing all datasets to {time_coords} time coordinates") From 2e854c6009fb867da61867d4d2cfa6ee47e16319 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Wed, 6 May 2026 22:13:45 +0000 Subject: [PATCH 7/8] FIX: VPT scans are now processed correctly. --- src/radclss/util/column_utils.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/radclss/util/column_utils.py b/src/radclss/util/column_utils.py index 82d611b..5eebd05 100644 --- a/src/radclss/util/column_utils.py +++ b/src/radclss/util/column_utils.py @@ -41,7 +41,7 @@ def _read_nexrad_cache(path): def _grab_90_degree_rays(radar): """Special case for column right over the radar in an RHI""" # Get the rays within 0.5 degrees of 90 degrees - ray = np.argmin(radar.elevation["data"] - 90.0) + ray = np.argmin(np.abs(radar.elevation["data"] - 90.0)) moment = {key: [] for key in radar.fields.keys()} # Determine the center of each gate for the subsetted rays. rhi_z = radar.range["data"] @@ -132,7 +132,7 @@ def _vpt_to_column_timeseries(radar, height_bins): For VPT scans elevation ≈ 90°, so gate range ≈ height above radar. Each ray becomes one time step; all rays are stacked into a (time, height) dataset. """ - zgate = radar.range["data"] + radar.altitude["data"][0] + zgate = radar.range["data"] base_vol_time = np.datetime64( pyart.util.datetime_from_radar(radar).isoformat(), "ns" @@ -157,22 +157,21 @@ def _vpt_to_column_timeseries(radar, height_bins): data_vars[key] = xr.DataArray(arr, dims=["time", "height"], attrs=attrs) ds = xr.Dataset(data_vars, coords={"height": zgate, "time": abs_times}) - ds = ds.drop_duplicates("time", keep="first") - + #ds = ds.dropna("height") valid_h = np.isfinite(ds["height"]) + print(ds["reflectivity"], np.sum(np.isnan(ds["reflectivity"].values))) if int(valid_h.sum()) > 0: try: - ds = ds.dropna("height").sortby("height").interp(height=height_bins) + ds = ds.sortby("height").interp(height=height_bins) except pd.errors.InvalidIndexError: ds = ( ds.drop_duplicates("height", keep="first") - .dropna("height") .sortby("height") .interp(height=height_bins) ) else: ds = ds.reindex(height=xr.DataArray(height_bins, dims="height", name="height")) - + print(ds["reflectivity"][:, :150]) dedup_times_s = ds["time"].values.astype("datetime64[s]") dedup_offsets = (ds["time"].values - base_vol_time).astype("timedelta64[s]") ds["base_time"] = xr.DataArray(dedup_times_s, dims="time") @@ -232,7 +231,8 @@ def _vpt_nan_fill(radar, height_bins): attrs=attrs, ) - ds = xr.Dataset(data_vars, coords={"height": height_bins, "time": abs_times}) + ds = xr.Dataset(data_vars, + coords={"height": height_bins, "time": abs_times}) ds = ds.drop_duplicates("time", keep="first") dedup_times_s = ds["time"].values.astype("datetime64[s]") From 48a17444a5672ad605486b7477c135508da0a237 Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Wed, 6 May 2026 22:14:13 +0000 Subject: [PATCH 8/8] VER: 2025.5.6 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index bbc16c1..fb8409d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "radclss" -version = "2026.4.20" +version = "2026.5.6" description = "Extracted Radar Columns and In Situ Sensors" readme = "README.md" requires-python = ">=3.10"