Add DataMet backend for Italian Radar data (#175)

* ADD: datamet backend for italian C-band data

* ADD: add unit test for datamet backend

* ENH: support for .gz files in DataMet reader

* [ENH] black style formatting

* remove version.py

* make DataMetBackend discoverable

* fix tests

* fix entrypoint

* ENH: add tests for datamet xarray backend

* FIX: remove version.py

* ENH: added more tests for datamet format

* ENH: black fix

* FIX: fix E721 in type comparison (ruff)

* DOC: edited history.md

* FIX: changed xradar.io.backends.DataMetBackendEntrypoint to 'datamet'

* FIX: black formatting

---------

Co-authored-by: Kai Mühlbauer <[email protected]>
Co-authored-by: Kai Mühlbauer <[email protected]>

docs/history.md

@@ -1,6 +1,7 @@
 # History
 
 ## Development Version (unreleased)
+* ADD: DataMet reader ({pull}`175`) by [@wolfidan](https://github.com/wolfidan).
 * FIX: Fix use of ruff, CI and numpy2 ({pull}`177`) by [@mgrover1](https://github.com/mgrover1) and [@kmuehlbauer](https://github.com/kmuehlbauer).
 * ADD: Add Alfonso to citation doc ({pull}`169`) by [@mgrover1](https://github.com/mgrover1).
 * ENH: Adding global variables and attributes to iris datatree ({pull}`166`) by [@aladinor](https://github.com/aladinor).

pyproject.toml

@@ -40,6 +40,7 @@ iris = "xradar.io.backends:IrisBackendEntrypoint"
 odim = "xradar.io.backends:OdimBackendEntrypoint"
 rainbow = "xradar.io.backends:RainbowBackendEntrypoint"
 nexradlevel2 = "xradar.io.backends:NexradLevel2BackendEntrypoint"
+datamet = "xradar.io.backends:DataMetBackendEntrypoint"
 
 
 [build-system]

tests/conftest.py

@@ -25,6 +25,11 @@ def odim_file2():
     return DATASETS.fetch("T_PAGZ35_C_ENMI_20170421090837.hdf")
 
 
+@pytest.fixture(scope="session")
+def datamet_file():
+    return DATASETS.fetch("H-000-VOL-ILMONTE-201907100700.tar.gz")
+
+
 @pytest.fixture(scope="session")
 def furuno_scn_file():
     return DATASETS.fetch("0080_20210730_160000_01_02.scn.gz")

tests/io/test_datamet.py

@@ -0,0 +1,58 @@
+import tarfile
+
+import pytest
+
+from xradar.io.backends import datamet
+from xradar.util import _get_data_file
+
+
+@pytest.fixture
+def data(datamet_file):
+    with _get_data_file(datamet_file, "file") as datametfile:
+        print(datametfile)
+        data = datamet.DataMetFile(datametfile)
+        assert data.filename == datametfile
+        print(data.scan_metadata)
+    data.get_sweep(0)
+    return data
+
+
+def test_file_types(data):
+    assert isinstance(data, datamet.DataMetFile)
+    assert isinstance(data.tarfile, tarfile.TarFile)
+
+
+def test_basic_content(data):
+    assert data.moments == ["UZ", "CZ", "V", "W", "ZDR", "PHIDP", "RHOHV", "KDP"]
+    assert len(data.data[0]) == len(data.moments)
+    assert data.first_dimension == "azimuth"
+    assert data.scan_metadata["origin"] == "ILMONTE"
+    assert data.scan_metadata["orig_lat"] == 41.9394
+    assert data.scan_metadata["orig_lon"] == 14.6208
+    assert data.scan_metadata["orig_alt"] == 710
+
+
+def test_moment_metadata(data):
+    mom_metadata = data.get_mom_metadata("UZ", 0)
+    print(mom_metadata)
+    assert mom_metadata["Rangeoff"] == 0.0
+    assert mom_metadata["Eloff"] == 16.05
+    assert mom_metadata["nlines"] == 360
+    assert mom_metadata["ncols"] == 493
+
+
+@pytest.mark.parametrize(
+    "moment, expected_value",
+    [
+        ("UZ", -3.5),
+        ("CZ", -3.5),
+        ("V", 2.3344999999999985),
+        ("W", 16.0),
+        ("ZDR", 0.6859999999999999),
+        ("PHIDP", 94.06648),
+        ("RHOHV", 1.9243000000000001),
+        ("KDP", 0.5190000000000001),
+    ],
+)
+def test_moment_data(data, moment, expected_value):
+    assert data.data[0][moment][(4, 107)] == expected_value

tests/io/test_io.py

@@ -17,6 +17,7 @@
 import xradar.io
 from xradar.io import (
     open_cfradial1_datatree,
+    open_datamet_datatree,
     open_gamic_datatree,
     open_iris_datatree,
     open_nexradlevel2_datatree,
@@ -793,6 +794,99 @@ def test_write_odim_source(rainbow_file2):
     assert ds["what"].attrs["source"].decode("utf-8") == "NOD:bewid,WMO:06477"
 
 
+def test_open_datamet_dataset(datamet_file):
+    # open first sweep group
+    ds = xr.open_dataset(
+        datamet_file,
+        group="sweep_0",
+        engine="datamet",
+    )
+    assert dict(ds.sizes) == {"azimuth": 360, "range": 493}
+    assert set(ds.data_vars) & (
+        sweep_dataset_vars | non_standard_sweep_dataset_vars
+    ) == {"DBTH", "DBZH", "KDP", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"}
+    assert ds.sweep_number == 0
+
+    # open last sweep group
+    ds = xr.open_dataset(
+        datamet_file,
+        group="sweep_10",
+        engine="datamet",
+    )
+    assert dict(ds.sizes) == {"azimuth": 360, "range": 1332}
+    assert set(ds.data_vars) & (
+        sweep_dataset_vars | non_standard_sweep_dataset_vars
+    ) == {"DBTH", "DBZH", "KDP", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"}
+    assert ds.sweep_number == 10
+
+
+def test_open_datamet_dataset_reindex(datamet_file):
+    # open first sweep group
+    reindex_angle = dict(start_angle=0, stop_angle=360, angle_res=2.0, direction=1)
+    ds = xr.open_dataset(
+        datamet_file,
+        group="sweep_10",
+        engine="datamet",
+        decode_coords=True,
+        reindex_angle=reindex_angle,
+    )
+    assert dict(ds.sizes) == {"azimuth": 180, "range": 1332}
+
+
+def test_open_datamet_datatree(datamet_file):
+    dtree = open_datamet_datatree(datamet_file)
+
+    # root_attrs
+    attrs = dtree.attrs
+    assert attrs["Conventions"] == "None"
+
+    # root vars
+    rvars = dtree.data_vars
+    assert rvars["volume_number"] == 0
+    assert rvars["platform_type"] == "fixed"
+    assert rvars["instrument_type"] == "radar"
+    assert rvars["time_coverage_start"] == "2019-07-10T07:00:00Z"
+    assert rvars["time_coverage_end"] == "2019-07-10T07:00:00Z"
+    np.testing.assert_almost_equal(rvars["latitude"].values, np.array(41.9394))
+    np.testing.assert_almost_equal(rvars["longitude"].values, np.array(14.6208))
+    np.testing.assert_almost_equal(rvars["altitude"].values, np.array(710))
+
+    # iterate over subgroups and check some values
+    moments = ["DBTH", "DBZH", "KDP", "PHIDP", "RHOHV", "VRADH", "WRADH", "ZDR"]
+    elevations = [16.1, 13.9, 11.0, 9.0, 7.0, 5.5, 4.5, 3.5, 2.5, 1.5, 0.5]
+    azimuths = [360] * 11
+    ranges = [493, 493, 493, 664, 832, 832, 1000, 1000, 1332, 1332, 1332]
+    i = 0
+    for grp in dtree.groups:
+        if grp.startswith("/sweep_"):
+            ds = dtree[grp].ds
+            assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]}
+            assert set(ds.data_vars) & (
+                sweep_dataset_vars | non_standard_sweep_dataset_vars
+            ) == set(moments)
+            assert set(ds.data_vars) & (required_sweep_metadata_vars) == set(
+                required_sweep_metadata_vars ^ {"azimuth", "elevation"}
+            )
+            assert set(ds.coords) == {
+                "azimuth",
+                "elevation",
+                "time",
+                "latitude",
+                "longitude",
+                "altitude",
+                "range",
+            }
+            assert np.round(ds.elevation.mean().values.item(), 1) == elevations[i]
+            assert ds.sweep_number == i
+            i += 1
+    # Try to reed single sweep
+    dtree = open_datamet_datatree(datamet_file, sweep=1)
+    assert len(dtree.groups) == 2
+    # Try to read list of sweeps
+    dtree = open_datamet_datatree(datamet_file, sweep=[1, 2])
+    assert len(dtree.groups) == 3
+
+
 @pytest.mark.parametrize("first_dim", ["time", "auto"])
 def test_cfradfial2_roundtrip(cfradial1_file, first_dim):
     dtree0 = open_cfradial1_datatree(cfradial1_file, first_dim=first_dim)

xradar/io/backends/__init__.py

@@ -16,6 +16,7 @@
 .. automodule:: xradar.io.backends.rainbow
 .. automodule:: xradar.io.backends.iris
 .. automodule:: xradar.io.backends.nexrad_level2
+.. automodule:: xradar.io.backends.datamet
 
 """
 
@@ -26,5 +27,6 @@
 from .odim import *  # noqa
 from .rainbow import *  # noqa
 from .nexrad_level2 import *  # noqa
+from .datamet import *  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]

xradar/io/backends/common.py

@@ -25,7 +25,7 @@
 def _maybe_decode(attr):
     try:
         # Decode the xr.DataArray differently than a byte string
-        if type(attr) == xr.core.dataarray.DataArray:
+        if type(attr) is xr.core.dataarray.DataArray:
             decoded_attr = attr.astype(str).str.rstrip()
         else:
             decoded_attr = attr.decode()