FIX: properly read CfRadial1 n_points files (#190)

* FIX: properly read CfRadial1 n_points files

* add history.md entry

* only strip dimensions which are available

* restructure coordinate assignments

* add test, bump open-radar-data

Author: kmuehlbauer

ci/notebooktests.yml

@@ -16,7 +16,7 @@ dependencies:
   - netCDF4
   - notebook
   - numpy
-  - open-radar-data>=0.1.0
+  - open-radar-data>=0.3.0
   - pip
   - pyproj
   - pytest

ci/unittests.yml

@@ -13,7 +13,7 @@ dependencies:
   - lat_lon_parser
   - netCDF4
   - numpy
-  - open-radar-data>=0.1.0
+  - open-radar-data>=0.3.0
   - pip
   - pyproj
   - pytest

docs/history.md

@@ -3,6 +3,7 @@
 ## Development Version
 
 * MNT: minimize CI ({pull}`192`) by [@kmuehlbauer](https://github.com/kmuehlbauer).
+* FIX: properly read CfRadial1 n_points files ({issue}`188`) by [@aladinor](https://github.com/aladinor), ({pull}`190`) by [@kmuehlbauer](https://github.com/kmuehlbauer).
 
 ## 0.6.0 (2024-08-05)
 

tests/conftest.py

@@ -15,6 +15,11 @@ def cfradial1_file(tmp_path_factory):
     return DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")
 
 
+@pytest.fixture(scope="session")
+def cfradial1n_file(tmp_path_factory):
+    return DATASETS.fetch("DES_VOL_RAW_20240522_1600.nc")
+
+
 @pytest.fixture(scope="session")
 def odim_file():
     return DATASETS.fetch("71_20181220_060628.pvol.h5")

tests/io/test_io.py

@@ -911,6 +911,59 @@ def test_cfradfial2_roundtrip(cfradial1_file, first_dim):
         xr.testing.assert_equal(dtree1[d1].ds, dtree2[d2].ds)
 
 
+def test_cfradial_n_points_file(cfradial1n_file):
+    dtree = open_cfradial1_datatree(
+        cfradial1n_file, first_dim="auto", site_coords=False
+    )
+    attrs = dtree.attrs
+
+    # root_attrs
+    assert attrs["Conventions"] == "CF-1.7"
+    assert attrs["version"] == "CF-Radial-1.4"
+    assert attrs["title"] == "VOL_A"
+    assert attrs["instrument_name"] == "Desio_Radar"
+    assert attrs["platform_is_mobile"] == "false"
+
+    # root vars
+    rvars = dtree.data_vars
+    assert rvars["volume_number"] == 1
+    assert rvars["platform_type"] == b"fixed"
+    assert rvars["instrument_type"] == b"radar"
+    assert rvars["time_coverage_start"] == b"2024-05-22T16:00:47Z"
+    assert rvars["time_coverage_end"] == b"2024-05-22T16:03:20Z"
+    np.testing.assert_almost_equal(rvars["latitude"].values, np.array(45.6272661))
+    np.testing.assert_almost_equal(rvars["longitude"].values, np.array(9.1963181))
+    np.testing.assert_almost_equal(rvars["altitude"].values, np.array(241.0))
+
+    # iterate over subgroups and check some values
+    moments = ["ZDR", "RHOHV", "KDP", "DBZ", "VEL", "PHIDP"]
+    elevations = [0.7, 1.3, 3.0, 5.0, 7.0, 10.0, 15.0, 25.0]
+    azimuths = [360] * 8
+    ranges = [416] * 5 + [383, 257, 157]
+    for grp in dtree.groups:
+        # only iterate sweep groups
+        if "sweep" not in grp:
+            continue
+        ds = dtree[grp].ds
+        i = int(ds.sweep_number.values)
+        assert i == int(grp[7:])
+        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",
+            "range",
+        }
+        assert np.round(ds.sweep_fixed_angle.values.item(), 1) == elevations[i]
+        assert ds.sweep_mode == "azimuth_surveillance"
+
+
 @pytest.mark.parametrize("sweep", ["sweep_0", 0, [0, 1], ["sweep_0", "sweep_1"]])
 @pytest.mark.parametrize(
     "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True

xradar/io/backends/cfradial1.py

@@ -34,7 +34,7 @@
 
 import numpy as np
 from datatree import DataTree
-from xarray import open_dataset
+from xarray import merge, open_dataset
 from xarray.backends import NetCDF4DataStore
 from xarray.backends.common import BackendEntrypoint
 from xarray.backends.store import StoreBackendEntrypoint
@@ -128,13 +128,20 @@ def _get_sweep_groups(
     ray_start_index = root.get("ray_start_index", False)
 
     # strip variables and attributes
-    anc_dims = list(set(root.dims) ^ {"time", "range", "sweep"})
+    anc_dims = set(root.dims) ^ {"time", "range", "sweep", "n_points"}
+    anc_dims &= set(root.dims)
+
     root = root.drop_dims(anc_dims)
 
     root = root.rename({"fixed_angle": "sweep_fixed_angle"})
 
     # conform to cfradial2 standard
     data = conform_cfradial2_sweep_group(root, optional, "time")
+    data_vars = {
+        k
+        for k, v in data.data_vars.items()
+        if any(d in v.dims for d in ["range", "n_points"])
+    }
 
     # which sweeps to load
     # sweep is assumed a list of strings with elements like "sweep_0"
@@ -172,23 +179,16 @@ def _get_sweep_groups(
             rslice = slice(0, current_ray_n_gates[0].values.astype(int))
             ds = ds.isel(range=rslice)
             ds = ds.isel(n_points=nslice)
-            ds = ds.stack(n_points=[dim0, "range"])
-            ds = ds.unstack("n_points")
-            # fix elevation/time additional range dimension in coordinate
-            ds = ds.assign_coords({"elevation": ds.elevation.isel(range=0, drop=True)})
-
-        # handling first dimension
-        # for CfRadial1 first dimension is time
-        if first_dim == "auto":
-            ds = ds.swap_dims({"time": dim0})
-            ds = ds.sortby(dim0)
-
-        # reassign azimuth/elevation coordinates
-        ds = ds.assign_coords({"azimuth": ds.azimuth})
-        ds = ds.assign_coords({"elevation": ds.elevation})
+            ds_vars = ds[data_vars]
+            ds_vars = merge([ds_vars, ds[[dim0, "range"]]])
+            ds_vars = ds_vars.stack(n_points=[dim0, "range"])
+            ds_vars = ds_vars.unstack("n_points")
+            ds = ds.drop_vars(ds_vars.data_vars)
+            ds = merge([ds, ds_vars])
 
         # assign site_coords
         if site_coords:
+
             ds = ds.assign_coords(
                 {
                     "latitude": root.latitude,
@@ -197,6 +197,20 @@ def _get_sweep_groups(
                 }
             )
 
+        # handling first dimension
+        # for CfRadial1 first dimension is time
+        if first_dim == "auto":
+            if "time" in ds.dims:
+                ds = ds.swap_dims({"time": dim0})
+            ds = ds.sortby(dim0)
+        else:
+            if "time" not in ds.dims:
+                ds = ds.swap_dims({dim0: "time"})
+            ds = ds.sortby("time")
+
+        # reassign azimuth/elevation coordinates
+        ds = ds.set_coords(["azimuth", "elevation"])
+
         sweep_groups[sw] = ds
 
     return sweep_groups

xradar/model.py

@@ -995,7 +995,12 @@ def determine_cfradial2_sweep_variables(obj, optional, dim0):
     keep_vars |= required_sweep_metadata_vars
     # all moment fields
     # todo: strip off non-conforming
-    keep_vars |= {k for k, v in obj.data_vars.items() if "range" in v.dims}
+    # this also handles cfradial1 n_points layout
+    keep_vars |= {
+        k
+        for k, v in obj.data_vars.items()
+        if any(d in v.dims for d in ["range", "n_points"])
+    }
     # optional variables
     if optional:
         keep_vars |= {k for k, v in obj.data_vars.items() if dim0 in v.dims}