ENH: Add a reader for nexrad level2 files

.github/workflows/ci.yml

@@ -59,6 +59,7 @@ jobs:
     - name: Install xradar
       run: |
         python -m pip install . --no-deps
+        python setup.py build_ext --inplace
     - name: Version Info
       run: |
         python -c "import xradar; print(xradar.version.version)"
@@ -104,6 +105,7 @@ jobs:
     - name: Install xradar
       run: |
         python -m pip install . --no-deps
+        python setup.py build_ext --inplace
     - name: Version Info
       run: |
         python -c "import xradar; print(xradar.version.version)"

.github/workflows/upload_pypi.yml

@@ -30,4 +30,5 @@ jobs:
           TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
         run: |
           python -m build
+          python setup.py build_ext --inplace
           twine upload dist/*

MANIFEST.in

@@ -8,4 +8,5 @@ recursive-include tests *
 recursive-exclude * __pycache__
 recursive-exclude * *.py[co]
 
+global-include *.pyx *pxd
 recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif

ci/notebooktests.yml

@@ -7,6 +7,7 @@ dependencies:
   - cmweather
   - codecov
   - coverage
+  - cython
   - dask
   - h5netcdf
   - h5py

ci/unittests.yml

@@ -6,6 +6,7 @@ dependencies:
   - codecov
   - cmweather
   - coverage
+  - cython
   - dask
   - h5netcdf
   - h5py

pyproject.toml

@@ -39,12 +39,15 @@ gamic = "xradar.io.backends:GamicBackendEntrypoint"
 iris = "xradar.io.backends:IrisBackendEntrypoint"
 odim = "xradar.io.backends:OdimBackendEntrypoint"
 rainbow = "xradar.io.backends:RainbowBackendEntrypoint"
+nexradlevel2 = "xradar.io.backends:NexradLevel2BackendEntrypoint"
 
 [build-system]
 requires = [
     "setuptools>=45",
     "wheel",
     "setuptools_scm[toml]>=7.0",
+    "cython",
+    "numpy"
 ]
 build-backend = "setuptools.build_meta"
 

setup.py

@@ -0,0 +1,24 @@
+from Cython.Build import cythonize
+from Cython.Compiler import Options
+from setuptools import Extension, setup
+
+# These are optional
+Options.docstrings = True
+Options.annotate = False
+
+# Modules to be compiled and include_dirs when necessary
+extensions = [
+    Extension(
+        "xradar.interpolate._nexrad_interpolate",
+        sources=["xradar/interpolate/_nexrad_interpolate.pyx"],
+    ),
+]
+
+
+# This is the function that is executed
+setup(
+    # external to be compiled
+    ext_modules=cythonize(
+        extensions, compiler_directives={"language_level": "3", "cpow": True}
+    ),
+)

tests/conftest.py

@@ -58,3 +58,8 @@ def iris0_file():
 @pytest.fixture(scope="session")
 def iris1_file():
     return DATASETS.fetch("SUR210819000227.RAWKPJV")
+
+
+@pytest.fixture(scope="session")
+def nexradlevel2_file():
+    return DATASETS.fetch("KATX20130717_195021_V06")

tests/io/backends/test_common.py

@@ -0,0 +1,9 @@
+from xradar.io.backends import common
+
+
+def test_lazy_dict():
+    d = common.LazyLoadDict({"key1": "value1", "key2": "value2"})
+    assert d["key1"] == "value1"
+    lazy_func = lambda: 999
+    d.set_lazy("lazykey1", lazy_func)
+    assert d["lazykey1"] == 999

tests/io/test_nexrad_level2.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+# Copyright (c) 2024, openradar developers.
+# Distributed under the MIT License. See LICENSE for more info.
+
+"""Tests for `xradar.io.nexrad_archive` module."""
+
+import xarray as xr
+
+from xradar.io.backends import open_nexradlevel2_datatree
+
+
+def test_open_nexradlevel2_datatree(nexradlevel2_file):
+    dtree = open_nexradlevel2_datatree(nexradlevel2_file)
+    ds = dtree["sweep_0"]
+    assert ds.attrs["instrument_name"] == "KATX"
+    assert ds.attrs["nsweeps"] == 16
+    assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters"
+    assert ds["DBZH"].shape == (719, 1832)
+    assert ds["DBZH"].dims == ("azimuth", "range")
+    assert int(ds.sweep_number.values) == 0
+
+
+def test_open_nexrad_level2_backend(nexradlevel2_file):
+    ds = xr.open_dataset(nexradlevel2_file, engine="nexradlevel2")
+    assert ds.attrs["instrument_name"] == "KATX"
+    assert ds.attrs["nsweeps"] == 16
+    assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters"
+    assert ds["DBZH"].shape == (719, 1832)
+    assert ds["DBZH"].dims == ("azimuth", "range")
+    assert int(ds.sweep_number.values) == 0

xradar/__init__.py

@@ -20,6 +20,7 @@
 # import subpackages
 from . import accessors  # noqa
 from . import georeference  # noqa
+from . import interpolate  # noqa
 from . import io  # noqa
 from . import model  # noqa
 from . import util  # noqa

xradar/interpolate/__init__.pxd

@@ -0,0 +1,11 @@
+#!/usr/bin/env python
+# Copyright (c) 2024, openradar developers.
+# Distributed under the MIT License. See LICENSE for more info.
+
+"""
+XRadar Interpolation
+====================
+
+"""
+
+from xradar.interpolate._nexrad_interpolate cimport * # noqa

xradar/interpolate/__init__.py

@@ -0,0 +1,8 @@
+#!/usr/bin/env python
+# Copyright (c) 2024, openradar developers.
+# Distributed under the MIT License. See LICENSE for more info.
+
+"""
+XRadar Interpolation
+====================
+"""

xradar/interpolate/_nexrad_interpolate.pyx

@@ -0,0 +1,110 @@
+"""
+Interpolation of NEXRAD moments from 1000 meter to 250 meter gate spacing.
+
+"""
+
+def _fast_interpolate_scan_4(
+        float[:, :] data, float[:] scratch_ray, float fill_value,
+        int start, int end, int moment_ngates, int linear_interp):
+    """ Interpolate a single NEXRAD moment scan from 1000 m to 250 m. """
+    # This interpolation scheme is only valid for NEXRAD data where a 4:1
+    # (1000 m : 250 m) interpolation is needed.
+    #
+    # The scheme here performs a linear interpolation between pairs of gates
+    # in a ray when the both of the gates are not masked (below threshold).
+    # When one of the gates is masked the interpolation changes to a nearest
+    # neighbor interpolation. Nearest neighbor is also performed at the end
+    # points until the new range bin would be centered beyond half of the range
+    # spacing of the original range.
+    #
+    # Nearest neighbor interpolation is performed when linear_interp is False,
+    # this is equivalent to repeating each gate four times in each ray.
+    #
+    # No transformation of the raw data is performed prior to interpolation, so
+    # reflectivity will be interpolated in dB units, velocity in m/s, etc,
+    # this may not be the best method for interpolation.
+    #
+    # This method was adapted from Radx and Py-ART
+    cdef int ray_num, i, interp_ngates
+    cdef float gate_val, next_val, delta
+
+    interp_ngates = 4 * moment_ngates  # number of gates interpolated
+
+    for ray_num in range(start, end+1):
+
+        # repeat each gate value 4 times
+        for i in range(moment_ngates):
+            gate_val = data[ray_num, i]
+            scratch_ray[i*4 + 0] = gate_val
+            scratch_ray[i*4 + 1] = gate_val
+            scratch_ray[i*4 + 2] = gate_val
+            scratch_ray[i*4 + 3] = gate_val
+
+        if linear_interp:
+            # linear interpolate
+            for i in range(2, interp_ngates - 4, 4):
+                gate_val = scratch_ray[i]
+                next_val = scratch_ray[i+4]
+                if gate_val == fill_value or next_val == fill_value:
+                    continue
+                delta = (next_val - gate_val) / 4.
+                scratch_ray[i+0] = gate_val + delta * 0.5
+                scratch_ray[i+1] = gate_val + delta * 1.5
+                scratch_ray[i+2] = gate_val + delta * 2.5
+                scratch_ray[i+3] = gate_val + delta * 3.5
+
+        for i in range(interp_ngates):
+            data[ray_num, i] = scratch_ray[i]
+
+
+def _fast_interpolate_scan_2(
+        float[:, :] data, float[:] scratch_ray, float fill_value,
+        int start, int end, int moment_ngates, int linear_interp):
+    """ Interpolate a single NEXRAD moment scan from 300 m to 150 m. """
+    # This interpolation scheme is only valid for NEXRAD TWDR data where a 2:1
+    # (300 m : 150 m) interpolation is needed.
+    #
+    # The scheme here performs a linear interpolation between pairs of gates
+    # in a ray when the both of the gates are not masked (below threshold).
+    # When one of the gates is masked the interpolation changes to a nearest
+    # neighbor interpolation. Nearest neighbor is also performed at the end
+    # points until the new range bin would be centered beyond half of the range
+    # spacing of the original range.
+    #
+    # Nearest neighbor interpolation is performed when linear_interp is False,
+    # this is equivalent to repeating each gate four times in each ray.
+    #
+    # No transformation of the raw data is performed prior to interpolation, so
+    # reflectivity will be interpolated in dB units, velocity in m/s, etc,
+    # this may not be the best method for interpolation.
+    #
+    # This method was adapted from Radx
+    cdef int ray_num, i, interp_ngates
+    cdef float gate_val, next_val, delta
+
+    interp_ngates = 2 * moment_ngates - 1 # number of gates interpolated
+
+    for ray_num in range(start, end+1):
+
+        # repeat each gate value 4 times
+        for i in range(moment_ngates):
+            gate_val = data[ray_num, i]
+            if i == moment_ngates - 1:
+                scratch_ray[i*2 + 0] = gate_val
+            else:
+                scratch_ray[i*2 + 0] = gate_val
+                scratch_ray[i*2 + 1] = gate_val
+
+        if linear_interp:
+            # linear interpolate
+            for i in range(1, interp_ngates - 2, 2):
+                gate_val = scratch_ray[i]
+                next_val = scratch_ray[i+2]
+                if gate_val == fill_value or next_val == fill_value:
+                    continue
+                delta = (next_val - gate_val) / 2.
+                scratch_ray[i+0] = gate_val + delta * 0.5
+                scratch_ray[i+1] = gate_val + delta * 1.5
+
+        for i in range(interp_ngates):
+            data[ray_num, i] = scratch_ray[i]

xradar/io/__init__.py

@@ -13,7 +13,24 @@
 .. automodule:: xradar.io.export
 
 """
-from .backends import *  # noqa
+from .backends import (
+    CfRadial1BackendEntrypoint,  # noqa
+    FurunoBackendEntrypoint,  # noqa
+    GamicBackendEntrypoint,  # noqa
+    IrisBackendEntrypoint,  # noqa
+    NexradLevel2BackendEntrypoint,  # noqa
+    OdimBackendEntrypoint,  # noqa
+    RainbowBackendEntrypoint,  # noqa
+    open_cfradial1_datatree,  # noqa
+    open_furuno_datatree,  # noqa
+    open_gamic_datatree,  # noqa
+    open_iris_datatree,  # noqa
+    open_nexradlevel2_datatree,  # noqa
+    open_odim_datatree,  # noqa
+    open_rainbow_datatree,  # noqa
+)
+
+# noqa
 from .export import *  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]

xradar/io/backends/__init__.py

@@ -15,6 +15,7 @@
 .. automodule:: xradar.io.backends.furuno
 .. automodule:: xradar.io.backends.rainbow
 .. automodule:: xradar.io.backends.iris
+.. automodule:: xradar.io.backends.nexrad_level2
 
 """
 
@@ -24,5 +25,7 @@
 from .iris import *  # noqa
 from .odim import *  # noqa
 from .rainbow import *  # noqa
-
-__all__ = [s for s in dir() if not s.startswith("_")]
+from .nexrad_level2 import (
+    NexradLevel2BackendEntrypoint,  # noqa
+    open_nexradlevel2_datatree,  # noqa
+)