ENH: NEXRAD Level 2 nyquist_velocity/unambiguous_range per sweep, mask below-threshold and range-folded codes

docs/history.md

@@ -37,6 +37,7 @@
 * FIX: CfRadial1 export writes string variables (``time_coverage_start``, ``instrument_type``, ``prt_mode``, ...) as char arrays instead of NC_STRING, which Py-ART could not read, and no longer writes ``"None"`` placeholder global attributes or a ``"None: ..."`` history ({issue}`417`) by [@syedhamidali](https://github.com/syedhamidali)
 * FIX: CfRadial1 export writes boolean attributes (e.g. NEXRAD ``mpda_vcp``, ``avset_enabled``) as ``"true"``/``"false"`` instead of failing, as netCDF attributes can't be bool ({issue}`418`) by [@syedhamidali](https://github.com/syedhamidali)
 * FIX: Halo Photonics ``.hpl`` reader: keep the last ray of each sweep (sweeps of one ray were empty and made opening fail, {issue}`303`), number sweeps from 0, take the header length from the ``****`` separator, raise a ``ValueError`` instead of printing on a line-count mismatch, and accept ``pathlib.Path`` and binary file-like objects ({issue}`430`, {pull}`442`) 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.12.0 (2026-04-21)

docs/importers.md

@@ -149,6 +149,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`

tests/io/test_nexrad_level2.py

@@ -758,7 +758,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}"
@@ -835,7 +835,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}"
@@ -867,6 +867,96 @@ 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_open_nexradlevel2_datatree_optional_groups(nexradlevel2_file):
     """Test that optional_groups=True includes metadata subgroups."""
     dtree = open_nexradlevel2_datatree(

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
@@ -69,8 +70,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,
@@ -169,6 +172,12 @@ def _concatenate_chunks(file_list):
     "CFP": "CCORH",
 }
 
+#: special moment codes, ICD 2620002, Table XVII-I (Data Moment
+#: Characteristics and Conversion for Data Names); valid for all moments
+#: and for Message Type 1 data
+BELOW_THRESHOLD = 0
+RANGE_FOLDED = 1
+
 # 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
@@ -1745,6 +1754,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:
@@ -1764,15 +1779,17 @@ 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:
+                data[data == RANGE_FOLDED] = 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(
@@ -1784,29 +1801,37 @@ def __getitem__(self, key):
 
 
 class NexradLevel2Store(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):
@@ -1841,6 +1866,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"
         )
@@ -1905,6 +1932,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),
@@ -1918,6 +1955,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
@@ -2015,6 +2056,7 @@ def open_dataset(
             filename_or_obj,
             group=group,
             lock=lock,
+            mask_and_scale=mask_and_scale,
             loaddata=False,
         )
 
@@ -2103,6 +2145,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.
@@ -2334,6 +2379,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():

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.