diff --git a/docs/history.md b/docs/history.md index d17ee0ec..4d61c6ce 100644 --- a/docs/history.md +++ b/docs/history.md @@ -24,6 +24,7 @@ * FIX: the calibration values of Furuno SCNX files are scaled (the header gives antenna gain, transmit power, radar constant and noise power in tenths) and have ``units``; the transmit power is ``xmit_power_h/v`` in dBm and is written to CfRadial1 files ({issue}`538`, {pull}`543`) by [@syedhamidali](https://github.com/syedhamidali) * FIX: the CfRadial1 reader gave each sweep the rays of another sweep when the sweeps of the file were not in time order (e.g. written by ``to_cfradial1`` from ``sweep=[5, 0]``); the sweeps are sliced in file order and sorted by time within each sweep ({issue}`550`, {pull}`551`) by [@syedhamidali](https://github.com/syedhamidali) * ENH: ``.xradar.reindex_angle`` and ``.xradar.reindex_range`` for Dataset (a sweep) and DataTree (every sweep), the same arguments as ``xradar.util.reindex_angle`` and ``reindex_range`` ({issue}`539`, {pull}`544`) by [@syedhamidali](https://github.com/syedhamidali) +* ENH: NEXRAD Level 2 reader exposes ``nyquist_velocity`` and ``unambiguous_range`` per sweep (MSG_31 radial data block / MSG_1 header) and masks the special moment codes below threshold (``0``) and range folded (``1``) to ``NaN`` via ``_FillValue`` when ``mask_and_scale=True``; raw codes stay available with ``mask_and_scale=False`` ({pull}`457`) by [@syedhamidali](https://github.com/syedhamidali) ## 0.13.0 (2026-10-21) diff --git a/docs/importers.md b/docs/importers.md index eadf14c2..8b4cac08 100644 --- a/docs/importers.md +++ b/docs/importers.md @@ -190,6 +190,20 @@ metadata. Finally, the xarray machinery returns a {py:class}`xarray:xarray.Datas with wanted group (eg. ``0``). Depending on the used backend kwargs several more functions are applied on that {py:class}`xarray:xarray.Dataset`. +Each sweep carries the ``nyquist_velocity`` (m s-1) and ``unambiguous_range`` +(meters) of the elevation cut, read from the radial data block (MSG_31) or the +digital radar data header (MSG_1). + +The special moment codes below threshold (``0``) and range folded (``1``) are +decoded to ``NaN`` with ``mask_and_scale=True`` (default). To tell them apart, +open the data with ``mask_and_scale=False`` (globally or per variable, e.g. +``mask_and_scale={"VRADH": False}``) and check the raw codes: + +```python +raw = xr.open_dataset(filename, engine="nexradlevel2", group="sweep_1", mask_and_scale=False) +range_folded = raw.VRADH == 1 +``` + ### open_nexradlevel2_datatree With {func}`~xradar.io.backends.nexrad_level2.open_nexradlevel2_datatree` diff --git a/tests/io/test_nexrad_level2.py b/tests/io/test_nexrad_level2.py index 3d80f66f..61b8c834 100644 --- a/tests/io/test_nexrad_level2.py +++ b/tests/io/test_nexrad_level2.py @@ -759,7 +759,7 @@ def test_open_nexradlevel2_datatree(nexradlevel2_file): # Verify a sample variable in one of the sweep groups (adjust as needed based on expected variables) sample_sweep = sweep_groups[0] - assert len(dtree[sample_sweep].data_vars) == 9 + assert len(dtree[sample_sweep].data_vars) == 11 assert ( "DBZH" in dtree[sample_sweep].data_vars ), f"DBZH should be a data variable in {sample_sweep}" @@ -836,7 +836,7 @@ def test_open_nexradlevel2_msg1_datatree(nexradlevel2_msg1_file): # Verify a sample variable in one of the sweep groups (adjust as needed based on expected variables) sample_sweep = sweep_groups[0] - assert len(dtree[sample_sweep].data_vars) == 6 + assert len(dtree[sample_sweep].data_vars) == 8 assert ( "DBZH" in dtree[sample_sweep].data_vars ), f"DBZH should be a data variable in {sample_sweep}" @@ -868,6 +868,110 @@ def test_open_nexradlevel2_msg1_datatree(nexradlevel2_msg1_file): assert dtree.attrs["scan_name"] == "VCP-0" +@pytest.mark.parametrize( + "sweep, nyquist, unambiguous_range", + [(0, 8.81, 466000.0), (1, 35.09, 117000.0), (9, 29.97, 137000.0)], +) +def test_nexradlevel2_nyquist_unambiguous_range( + nexradlevel2_file, sweep, nyquist, unambiguous_range +): + """Nyquist velocity and unambiguous range from the MSG_31 RAD block.""" + ds = open_dataset(nexradlevel2_file, engine="nexradlevel2", group=f"sweep_{sweep}") + assert ds.nyquist_velocity.dims == () + assert ds.nyquist_velocity.attrs == { + "standard_name": "nyquist_velocity", + "units": "m s-1", + } + np.testing.assert_allclose(ds.nyquist_velocity, nyquist) + assert ds.unambiguous_range.dims == () + assert ds.unambiguous_range.attrs["units"] == "meters" + np.testing.assert_allclose(ds.unambiguous_range, unambiguous_range) + with NEXRADLevel2File(nexradlevel2_file, loaddata=False) as nex: + rad = nex.msg_31_data_header[sweep]["msg_31_data_header"]["RAD"] + np.testing.assert_allclose(ds.nyquist_velocity, rad["nyquist_vel"] / 100.0) + + dtree = open_nexradlevel2_datatree(nexradlevel2_file, sweep=[sweep]) + np.testing.assert_allclose(dtree[f"sweep_{sweep}"].nyquist_velocity, nyquist) + + +def test_nexradlevel2_msg1_nyquist_unambiguous_range(nexradlevel2_msg1_file): + """Nyquist velocity and unambiguous range from the MSG_1 header.""" + dtree = open_nexradlevel2_datatree(nexradlevel2_msg1_file, sweep=[0, 1, 9]) + # surveillance-only cut has no Nyquist velocity + assert np.isnan(dtree["sweep_0"].nyquist_velocity) + np.testing.assert_allclose(dtree["sweep_0"].unambiguous_range, 466000.0) + np.testing.assert_allclose(dtree["sweep_1"].nyquist_velocity, 25.37) + np.testing.assert_allclose(dtree["sweep_1"].unambiguous_range, 148000.0) + np.testing.assert_allclose(dtree["sweep_9"].nyquist_velocity, 27.41) + np.testing.assert_allclose(dtree["sweep_9"].unambiguous_range, 137000.0) + + +@pytest.mark.parametrize("moment", ["DBZH", "VRADH", "WRADH"]) +def test_nexradlevel2_mask_special_codes(nexradlevel2_file, moment): + """Below threshold (0) and range folded (1) are masked by default.""" + raw = open_dataset( + nexradlevel2_file, engine="nexradlevel2", group="sweep_1", mask_and_scale=False + )[moment] + ds = open_dataset(nexradlevel2_file, engine="nexradlevel2", group="sweep_1")[moment] + # raw codes are kept with mask_and_scale=False + assert raw.dtype == np.uint8 + assert raw.attrs["_FillValue"] == 0 + below = (raw == 0).values + folded = (raw == 1).values + assert below.any() and folded.any() + + # both codes decode to NaN, all other values are valid + np.testing.assert_array_equal(ds.isnull().values, below | folded) + assert ds.encoding["_FillValue"] == 0 + scale, offset = raw.attrs["scale_factor"], raw.attrs["add_offset"] + valid = ~(below | folded) + np.testing.assert_allclose( + ds.values[valid], raw.values[valid] * scale + offset, rtol=1e-12 + ) + # decoded special code values do not show up as data + assert float(ds.min()) > offset + scale + + +def test_nexradlevel2_mask_special_codes_dual_pol(nexradlevel2_file): + """Special codes are masked for 16 bit moments (PHIDP) and RHOHV, ZDR.""" + raw = open_dataset( + nexradlevel2_file, engine="nexradlevel2", group="sweep_0", mask_and_scale=False + ) + ds = open_dataset(nexradlevel2_file, engine="nexradlevel2", group="sweep_0") + for moment in ["ZDR", "PHIDP", "RHOHV"]: + codes = raw[moment].isin([0, 1]).values + assert codes.any() + np.testing.assert_array_equal(ds[moment].isnull().values, codes) + assert raw["PHIDP"].dtype == np.uint16 + + +def test_nexradlevel2_mask_and_scale_mapping(nexradlevel2_file): + """Per-variable mask_and_scale keeps raw codes only for that variable.""" + ds = open_dataset( + nexradlevel2_file, + engine="nexradlevel2", + group="sweep_1", + mask_and_scale={"VRADH": False}, + ) + assert ds.VRADH.dtype == np.uint8 + assert (ds.VRADH == 1).any() + assert ds.DBZH.isnull().any() + + +def test_nexradlevel2_mask_cfp_flags(closing): + """CFP codes 0-7 are flags, not clutter power (ICD Table XVII-I, note 30).""" + fname = DATASETS.fetch("KILX20230629_154426_V06") + raw = closing( + open_dataset( + fname, engine="nexradlevel2", group="sweep_0", mask_and_scale=False + ) + ) + ds = closing(open_dataset(fname, engine="nexradlevel2", group="sweep_0")) + assert (raw.CCORH == 2).any() + np.testing.assert_array_equal(ds.CCORH.isnull().values, (raw.CCORH < 8).values) + assert float(ds.CCORH.min()) >= 0.0 + + def test_open_nexradlevel2_datatree_optional_groups(nexradlevel2_file): """Test that optional_groups=True includes metadata subgroups.""" dtree = open_nexradlevel2_datatree( diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 323977f2..5d83b53b 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -41,6 +41,7 @@ import struct import warnings from collections import OrderedDict, defaultdict +from collections.abc import Mapping import numpy as np import xarray as xr @@ -86,8 +87,10 @@ get_elevation_attrs, get_latitude_attrs, get_longitude_attrs, + get_nyquist_velocity_attrs, get_range_attrs, get_time_attrs, + get_unambiguous_range_attrs, moment_attrs, radar_calibration_subgroup, radar_parameters_subgroup, @@ -186,6 +189,15 @@ def _concatenate_chunks(file_list): "CFP": "CCORH", } +#: special moment codes, ICD 2620002, Table XVII-I (Data Moment +#: Characteristics and Conversion for Data Names), note 21; also for +#: Message Type 1 data +BELOW_THRESHOLD = 0 +RANGE_FOLDED = 1 +#: CFP codes 0-7 are flags (filter not applied, point filter, dual pol +#: only, reserved), its data start at 8 (Table XVII-I, note 30) +CFP_FIRST_DATA = 8 + # NEXRAD Level II file structures and sizes # The deails on these structures are documented in: # "Interface Control Document for the Achive II/User" RPG Build 12.0 @@ -1765,6 +1777,12 @@ def __init__(self, datastore, name, var): width = {8: 1, 16: 2}[word_size] self.dtype = np.dtype(f">u{width}") self.shape = (nrays, nbins) + mask_and_scale = datastore.mask_and_scale + if isinstance(mask_and_scale, Mapping): + mask_and_scale = mask_and_scale.get(nexrad_mapping.get(name, name), True) + # range folded is merged into below threshold (``_FillValue``) when + # masking, so both special codes decode to NaN + self.mask_range_folded = mask_and_scale is None or bool(mask_and_scale) def _getitem(self, key): with self.datastore.lock: @@ -1784,15 +1802,18 @@ def _getitem(self, key): x = np.uint16(0x7FF) else: x = np.uint8(0xFF) - if len(data[0]) < self.shape[1]: - return np.pad( - np.vstack(data) & x, - ((0, 0), (0, self.shape[1] - len(data[0]))), + data = np.vstack(data) & x + if self.mask_range_folded: + first_data = CFP_FIRST_DATA if self.name == "CFP" else RANGE_FOLDED + 1 + data[data < first_data] = BELOW_THRESHOLD + if data.shape[1] < self.shape[1]: + data = np.pad( + data, + ((0, 0), (0, self.shape[1] - data.shape[1])), mode="constant", - constant_values=0, - )[key] - else: - return (np.vstack(data) & x)[key] + constant_values=BELOW_THRESHOLD, + ) + return data[key] def __getitem__(self, key): return indexing.explicit_indexing_adapter( @@ -1804,29 +1825,37 @@ def __getitem__(self, key): class NexradLevel2Store(_ManagedStoreMixin, AbstractDataStore): - def __init__(self, manager, group=None, lock=NEXRADL2_LOCK): + def __init__(self, manager, group=None, lock=NEXRADL2_LOCK, mask_and_scale=True): self._manager = manager self._group = int(group[6:]) self._filename = self.filename self.lock = ensure_lock(lock) + self.mask_and_scale = mask_and_scale @classmethod - def open(cls, filename, mode="r", group=None, lock=None, **kwargs): + def open( + cls, filename, mode="r", group=None, lock=None, mask_and_scale=True, **kwargs + ): if lock is None: lock = NEXRADL2_LOCK manager = CachingFileManager( NEXRADLevel2File, filename, mode=mode, kwargs=kwargs ) - return cls(manager, group=group, lock=lock) + return cls(manager, group=group, lock=lock, mask_and_scale=mask_and_scale) @classmethod - def open_groups(cls, filename, groups, mode="r", lock=None, **kwargs): + def open_groups( + cls, filename, groups, mode="r", lock=None, mask_and_scale=True, **kwargs + ): if lock is None: lock = NEXRADL2_LOCK manager = CachingFileManager( NEXRADLevel2File, filename, mode=mode, kwargs=kwargs ) - return {group: cls(manager, group=group, lock=lock) for group in groups} + return { + group: cls(manager, group=group, lock=lock, mask_and_scale=mask_and_scale) + for group in groups + } @property def filename(self): @@ -1861,6 +1890,8 @@ def open_store_variable(self, name, var): attrs = {key: mapping[key] for key in moment_attrs if key in mapping} attrs["scale_factor"] = 1.0 / var["scale"] attrs["add_offset"] = -var["offset"] / var["scale"] + # below threshold (and range folded, see NexradLevel2ArrayWrapper) + attrs["_FillValue"] = BELOW_THRESHOLD attrs["coordinates"] = ( "elevation azimuth range latitude longitude altitude time" ) @@ -1925,6 +1956,16 @@ def open_store_coordinates(self): fixed_angle = self.root.msg_31_header[header_idx][0]["elevation_angle"] fixed_angle *= angle_scale + # nyquist velocity and unambiguous range, constant over the sweep + # MSG_31 RAD block (Table XVII-H) or MSG_1 header (Table III) + # scaled by 0.01 m s-1 and 0.1 km + rad = self.ds.get("sweep_constant_data", {}).get("RAD", msg_31_header[0]) + nyquist = rad.get("nyquist_vel", 0) / 100.0 + unambiguous_range = rad.get("unambig_range", 0) * 100.0 + # surveillance-only MSG_1 cuts carry zero values + nyquist = nyquist if nyquist > 0 else np.nan + unambiguous_range = unambiguous_range if unambiguous_range > 0 else np.nan + coords = { "azimuth": Variable((dim,), azimuth, get_azimuth_attrs(), encoding), "elevation": Variable((dim,), elevation, get_elevation_attrs(), encoding), @@ -1938,6 +1979,10 @@ def open_store_coordinates(self): "prt_mode": Variable((), prt_mode), "follow_mode": Variable((), follow_mode), "sweep_fixed_angle": Variable((), fixed_angle), + "nyquist_velocity": Variable((), nyquist, get_nyquist_velocity_attrs()), + "unambiguous_range": Variable( + (), unambiguous_range, get_unambiguous_range_attrs() + ), } return coords @@ -2080,6 +2125,7 @@ def open_dataset( filename_or_obj, group=group, lock=lock, + mask_and_scale=mask_and_scale, loaddata=False, ) @@ -2347,6 +2393,9 @@ def open_nexradlevel2_datatree( mask_and_scale : bool, optional If True, replaces values in the dataset that match `_FillValue` with NaN and applies scale and offset adjustments. Default is True. + The special moment codes below threshold (0) and range folded (1) are + both decoded to NaN. Use ``mask_and_scale=False`` to get the raw codes, + e.g. to keep range folded gates distinguishable. decode_times : bool, optional If True, decodes time variables according to CF conventions. Default is True. @@ -2475,6 +2524,7 @@ def open_sweeps_as_dict( filename=filename_or_obj, lock=lock, groups=sweeps, + mask_and_scale=mask_and_scale, ) groups_dict = {} for path_group, store in stores.items(): diff --git a/xradar/model.py b/xradar/model.py index 1616ab0e..7b37a37d 100644 --- a/xradar/model.py +++ b/xradar/model.py @@ -37,6 +37,7 @@ "get_nyquist_velocity_attrs", "get_range_attrs", "get_time_attrs", + "get_unambiguous_range_attrs", "moment_attrs", "sweep_vars_mapping", "get_azimuth_dataarray", @@ -754,6 +755,14 @@ def get_nyquist_velocity_attrs(units="m s-1"): return nyquist_attrs +def get_unambiguous_range_attrs(units="meters"): + unambiguous_range_attrs = { + "long_name": "unambiguous_range", + "units": units, + } + return unambiguous_range_attrs + + def get_moment_attrs(moment): """Get Radar Moment CF attributes.