diff --git a/tests/io/test_odim.py b/tests/io/test_odim.py index 77886c9f..56338245 100644 --- a/tests/io/test_odim.py +++ b/tests/io/test_odim.py @@ -57,6 +57,186 @@ def test_get_azimuth_where(nrays): assert udiff[0] == 360.0 / nrays +def test_odim_azimuth_nominal(odim_file): + import xarray as xr + + ds = xr.open_dataset( + odim_file, engine="odim", angle_spec="nominal", group="sweep_0" + ) + + azimuth = ds.azimuth.values + udiff = np.unique(np.diff(azimuth)) + assert len(azimuth) == 360 + assert len(udiff) == 1 + assert udiff[0] == 1.0 + + +def test_odim_azimuth_auto(odim_file): + import xarray as xr + + ds = xr.open_dataset(odim_file, engine="odim", angle_spec="auto", group="sweep_0") + + azimuth = ds.azimuth.values + udiff = np.unique(np.diff(azimuth)) + assert len(azimuth) == 360 + assert len(udiff) == 1 + assert udiff[0] == 1.0 + + +def test_odim_azimuth_reindex_deprecation(odim_file): + import xarray as xr + + with pytest.warns( + DeprecationWarning, match="The 'reindex_angle' kwarg is deprecated" + ): + ds = xr.open_dataset( + odim_file, engine="odim", reindex_angle=False, group="sweep_0" + ) + azimuth = ds.azimuth.values + udiff = np.unique(np.diff(azimuth)) + assert len(azimuth) == 360 + assert len(udiff) == 1 + assert udiff[0] == 1.0 + + +def test_odim_azimuth_raw(odim_file): + import xarray as xr + + ds = xr.open_dataset(odim_file, engine="odim", group="sweep_0", angle_spec="raw") + azimuth = ds.azimuth.values + assert len(azimuth) == 360 + # raw should still give uniform spacing for this file + udiff = np.unique(np.diff(azimuth)) + assert len(udiff) == 1 + assert udiff[0] == 1.0 + + +@pytest.mark.parametrize("angle_spec", ["nominal", "auto", "raw"]) +def test_odim_azimuth_nonstandard_nrays(tmp_path, angle_spec): + import h5netcdf + + nrays = 361 + ascale = 360.0 / nrays + + startazA = np.arange(0, 360, ascale, dtype=np.float32) + stopazA = np.arange(ascale, 360 + ascale, ascale, dtype=np.float32) + + filepath = tmp_path / "test_nonstandard.h5" + with h5netcdf.File(filepath, "w") as f: + f.attrs["Conventions"] = "ODIM_H5/V2_2" + wg = f.create_group("where") + wg.attrs["lon"] = 0.0 + wg.attrs["lat"] = 0.0 + wg.attrs["height"] = 0.0 + dg = f.create_group("dataset1") + dw = dg.create_group("where") + dw.attrs["nrays"] = nrays + dw.attrs["nbins"] = 100 + dw.attrs["rstart"] = 0.0 + dw.attrs["rscale"] = 100.0 + dw.attrs["elangle"] = 0.5 + dw.attrs["a1gate"] = 0 + dh = dg.create_group("how") + dh.attrs["startazA"] = startazA + dh.attrs["stopazA"] = stopazA + dh.attrs["startazT"] = np.zeros(nrays, dtype=np.float64) + dh.attrs["stopazT"] = np.ones(nrays, dtype=np.float64) + dwhat = dg.create_group("what") + dwhat.attrs["quantity"] = "DBZH" + dwhat.attrs["startdate"] = "20000101" + dwhat.attrs["starttime"] = "000000" + dwhat.attrs["enddate"] = "20000101" + dwhat.attrs["endtime"] = "000030" + dg.create_group("data1") + + import xarray as xr + + if angle_spec == "nominal": + context = pytest.raises(ValueError, match="Unexpected number of rays") + elif angle_spec == "auto": + context = pytest.warns(RuntimeWarning, match="Unexpected number of rays") + wanted = np.arange(0.5, 360, 1.0) + else: + context = nullcontext() + # should match per-ray midpoint + wanted = (startazA + np.where(stopazA < startazA, stopazA + 360, stopazA)) / 2 + wanted[wanted >= 360] -= 360 + + with context: + ds = xr.open_dataset( + filepath, engine="odim", angle_spec=angle_spec, group="sweep_0" + ) + + if angle_spec != "nominal": + azimuth = ds.azimuth.values + assert len(azimuth) == 360 if angle_spec == "auto" else nrays + np.testing.assert_array_almost_equal(azimuth, wanted, decimal=4) + + +@pytest.mark.parametrize("angle_spec", ["nominal", "auto", "raw"]) +def test_odim_azimuth_fallback_nonstandard_nrays_with_duplicate(tmp_path, angle_spec): + import h5netcdf + + nrays = 361 + startazA = create_startazA(nrays) + stopazA = create_stopazA(nrays) + + filepath = tmp_path / "test_nonstandard_dup.h5" + with h5netcdf.File(filepath, "w") as f: + f.attrs["Conventions"] = "ODIM_H5/V2_2" + wg = f.create_group("where") + wg.attrs["lon"] = 0.0 + wg.attrs["lat"] = 0.0 + wg.attrs["height"] = 0.0 + dg = f.create_group("dataset1") + dw = dg.create_group("where") + dw.attrs["nrays"] = nrays + dw.attrs["nbins"] = 100 + dw.attrs["rstart"] = 0.0 + dw.attrs["rscale"] = 100.0 + dw.attrs["elangle"] = 0.5 + dw.attrs["a1gate"] = 0 + dh = dg.create_group("how") + dh.attrs["startazA"] = startazA + dh.attrs["stopazA"] = stopazA + dh.attrs["startazT"] = np.zeros(nrays, dtype=np.float64) + dh.attrs["stopazT"] = np.ones(nrays, dtype=np.float64) + dwhat = dg.create_group("what") + dwhat.attrs["quantity"] = "DBZH" + dwhat.attrs["startdate"] = "20000101" + dwhat.attrs["starttime"] = "000000" + dwhat.attrs["enddate"] = "20000101" + dwhat.attrs["endtime"] = "000030" + dg.create_group("data1") + + import xarray as xr + + if angle_spec == "nominal": + context = pytest.raises(ValueError, match="Unexpected number of rays") + elif angle_spec == "auto": + context = pytest.warns(RuntimeWarning, match="Unexpected number of rays") + wanted = np.arange(0.5, 360, 1.0) + else: + context = nullcontext() + # should match per-ray midpoint + wanted = (startazA + np.where(stopazA < startazA, stopazA + 360, stopazA)) / 2 + wanted[wanted >= 360] -= 360 + + with context: + ds = xr.open_dataset( + filepath, engine="odim", angle_spec=angle_spec, group="sweep_0" + ) + + if angle_spec != "nominal": + azimuth = ds.azimuth.values + assert len(azimuth) == 360 if angle_spec == "auto" else nrays + # assert len(azimuth) == nrays + # wanted = (startazA + np.where(stopazA < startazA, stopazA + 360, stopazA)) / 2 + # wanted[wanted >= 360] -= 360 + print(azimuth) + np.testing.assert_array_almost_equal(azimuth, wanted, decimal=4) + + @pytest.mark.parametrize( "ang", [("az_angle", "elevation"), ("az_angle", "elevation"), ("elangle", "azimuth")], diff --git a/xradar/io/backends/odim.py b/xradar/io/backends/odim.py index f0825ec2..f1e12f42 100644 --- a/xradar/io/backends/odim.py +++ b/xradar/io/backends/odim.py @@ -100,8 +100,6 @@ def _get_azimuth_how(how): startaz = how["startazA"] stopaz = how.get("stopazA", False) if stopaz is False: - # stopazA missing - # create from startazA stopaz = np.roll(startaz, -1) stopaz[-1] += 360 zero_index = np.where(stopaz < startaz) @@ -116,6 +114,26 @@ def _get_azimuth_where(where): return np.arange(res / 2.0, 360.0, res, dtype="float32") +def _get_azimuth_nominal(nrays, how): + + if nrays not in [180, 360, 450, 720, 900]: + raise ValueError(f"xradar: Unexpected number of rays ({nrays})") + + ascale = 360 / nrays + try: + astart = how["startazA"][0] + astart = min(astart, astart - 360, key=abs) + astart = np.round(astart / (ascale / 2)) * ascale / 2 + except (KeyError, TypeError, AttributeError): + try: + astart = how["astart"] + except KeyError: + warnings.warn("xradar: No startazA or astart found, using 0 as default.") + astart = 0 + azimuth = np.arange(astart + ascale / 2, 360, ascale) + return azimuth + + def _get_fixed_dim_and_angle(where): dim = "elevation" @@ -407,12 +425,20 @@ class _OdimH5NetCDFMetadata(_H5NetCDFMetadata): h5netcdf filehandle. group : str odim group to acquire + angle_spec : str + ``"nominal"`` for uniform sweep from first ``startazA``, + ``"auto"`` try ``nominal`` first, fallback to ``raw``, + ``"raw"`` for per-ray midpoint from ``startazA``/``stopazA``. Returns ------- object : metadata object """ + def __init__(self, fileobj, group, angle_spec="nominal"): + super().__init__(fileobj, group) + self._angle_spec = angle_spec + @property def ds_what(self): return self._get_dset_what() @@ -461,10 +487,25 @@ def a1gate(self): @property def _azimuth(self): - try: - azimuth = _get_azimuth_how(self.how) - except (AttributeError, KeyError, TypeError): - azimuth = _get_azimuth_where(self.where) + nrays = self.where["nrays"] + if self._angle_spec in ["nominal", "auto"]: + try: + azimuth = _get_azimuth_nominal(nrays, self.how) + except ValueError as exc: + if self._angle_spec == "nominal": + raise + warnings.warn( + "xradar: Failed to extract nominal angles. " + f"Original error: {exc}", + category=RuntimeWarning, + stacklevel=2, + ) + self._angle_spec = "raw" + if self._angle_spec == "raw": + try: + azimuth = _get_azimuth_how(self.how) + except (AttributeError, KeyError, TypeError): + azimuth = _get_azimuth_where(self.where) return Variable((self.dim0,), azimuth, get_azimuth_attrs()) @property @@ -607,6 +648,7 @@ def __init__( store, group=None, lock=False, + angle_spec="nominal", ): if not isinstance(store, OdimStore): raise TypeError( @@ -619,11 +661,14 @@ def __init__( self._filename = store.filename self.is_remote = is_remote_uri(self._filename) self.lock = ensure_lock(lock) + self._angle_spec = angle_spec @property def root(self): with self._manager.acquire_context(False) as root: - return _OdimH5NetCDFMetadata(root, self._group.lstrip("/")) + return _OdimH5NetCDFMetadata( + root, self._group.lstrip("/"), angle_spec=self._angle_spec + ) def _acquire(self, needs_lock=True): with self._manager.acquire_context(needs_lock) as root: @@ -663,7 +708,7 @@ def get_variables(self): class OdimStore(AbstractDataStore): """Store for reading ODIM dataset groups via h5netcdf.""" - def __init__(self, manager, group=None, lock=False): + def __init__(self, manager, group=None, lock=False, angle_spec="nominal"): if isinstance(manager, (h5netcdf.File, h5netcdf.Group)): if group is None: root, group = find_root_and_group(manager) @@ -683,6 +728,7 @@ def __init__(self, manager, group=None, lock=False): self.lock = ensure_lock(lock) self._substore = None self._need_time_recalc = False + self._angle_spec = angle_spec @classmethod def open( @@ -695,6 +741,7 @@ def open( invalid_netcdf=None, phony_dims=None, decode_vlen_strings=True, + angle_spec="nominal", ): if isinstance(filename, bytes): raise ValueError( @@ -717,7 +764,7 @@ def open( lock = False manager = CachingFileManager(h5netcdf.File, filename, mode=mode, kwargs=kwargs) - return cls(manager, group=group, lock=lock) + return cls(manager, group=group, lock=lock, angle_spec=angle_spec) @property def filename(self): @@ -741,6 +788,7 @@ def substore(self): self, group=group, lock=self.lock, + angle_spec=self._angle_spec, ) for group in subgroups ] @@ -780,8 +828,15 @@ class OdimBackendEntrypoint(BackendEntrypoint): first dimension will be either ``azimuth`` or ``elevation`` depending on type of sweep. Defaults to ``auto``. reindex_angle : bool or dict + Deprecated, Use angle_spec instead. See below. Defaults to False, no reindexing. Given dict should contain the kwargs to reindex_angle. Only invoked if `decode_coord=True`. + angle_spec : "nominal", "auto", "raw" or dict + ``"nominal"`` (future default) for uniform sweep from first ``startazA``, + ``"auto"`` try ``nominal`` first, fallback to ``raw``, automatic reindexing, + ``"raw"`` (deprecated default) for per-ray midpoint from ``startazA``/``stopazA``, + ``reindex_dict`` use ``raw`` with reindexing. + Only invoked if `decode_coord=True`. fix_second_angle : bool If True, fixes erroneous second angle data. Defaults to ``False``. site_as_coords : bool @@ -810,10 +865,35 @@ def open_dataset( phony_dims="access", decode_vlen_strings=True, first_dim="auto", - reindex_angle=False, + reindex_angle=None, + angle_spec=None, fix_second_angle=False, site_as_coords=True, ): + # kwarg change deprecation + # reindex_angle -> angle_spec + # old kwarg requested + if reindex_angle is not None and angle_spec is None: + warnings.warn( + "The 'reindex_angle' kwarg is deprecated and will be removed in a future release. " + "Use 'angle_spec' instead: angle_spec=\"raw\" replaces reindex_angle=False, " + "'angle_spec=\"nominal\"' returns nominal angles, angle_spec=reindex_dict " + "replaces reindex_angle=reindex_dict.", + DeprecationWarning, + ) + if reindex_angle is False: + angle_spec = "raw" + else: + angle_spec = reindex_angle + if angle_spec is None: + warnings.warn( + "In a future release nominal angles will be returned by default. " + "To silence this warning use: 'angle_spec=\"raw\"' to return raw measured " + "angles or 'angle_spec=\"nominal\"' to return nominal angles.", + DeprecationWarning, + ) + angle_spec = "raw" + if isinstance(filename_or_obj, io.IOBase): filename_or_obj.seek(0) @@ -824,6 +904,7 @@ def open_dataset( invalid_netcdf=invalid_netcdf, phony_dims=phony_dims, decode_vlen_strings=decode_vlen_strings, + angle_spec=angle_spec, ) store_entrypoint = StoreBackendEntrypoint() @@ -849,10 +930,15 @@ def open_dataset( ds.encoding["engine"] = "odim" # handle duplicates and reindex - if decode_coords and reindex_angle is not False: - ds = ds.pipe(util.remove_duplicate_rays) - ds = ds.pipe(util.reindex_angle, **reindex_angle) - ds = ds.pipe(util.ipol_time, **reindex_angle) + if decode_coords: + if angle_spec == "auto": + angle_spec = util.extract_angle_parameters(ds) + allowed_keys = {"start_angle", "stop_angle", "angle_res", "direction"} + angle_spec = {k: v for k, v in angle_spec.items() if k in allowed_keys} + if isinstance(angle_spec, dict): + ds = ds.pipe(util.remove_duplicate_rays) + ds = ds.pipe(util.reindex_angle, **angle_spec) + ds = ds.pipe(util.ipol_time, **angle_spec) # handling first dimension dim0 = "elevation" if ds.sweep_mode.load() == "rhi" else "azimuth" @@ -898,8 +984,15 @@ def open_odim_datatree(filename_or_obj, **kwargs): first dimension will be either ``azimuth`` or ``elevation`` depending on type of sweep. Defaults to ``auto``. reindex_angle : bool or dict + Deprecated, Use angle_spec instead. See below. Defaults to False, no reindexing. Given dict should contain the kwargs to reindex_angle. Only invoked if `decode_coord=True`. + angle_spec : "nominal", "auto", "raw" or dict + ``"nominal"`` (future default) for uniform sweep from first ``startazA``, + ``"auto"`` try ``nominal`` first, fallback to ``raw``, automatic reindexing, + ``"raw"`` (deprecated default) for per-ray midpoint from ``startazA``/``stopazA``, + ``reindex_dict`` use ``raw`` with reindexing. + Only invoked if `decode_coord=True`. fix_second_angle : bool If True, fixes erroneous second angle data. Defaults to ``False``. site_as_coords : bool