diff --git a/docs/history.md b/docs/history.md index d03e6815..f670a742 100644 --- a/docs/history.md +++ b/docs/history.md @@ -4,13 +4,15 @@ This is the first version which uses datatree directly from xarray. Thus, xarray is pinned to version >= 2024.10.0. +* FIX: Convert volumes to_cfradial1 containing sweeps with different range and azimuth shapes, raise for different range bin sizes ({issue}`233`) by [@syedhamidali](https://github.com/syedhamidali), ({pull}`234`) by [@kmuehlbauer](https://github.com/kmuehlbauer). +* FIX: Correctly handle 8bit/16bit, big-endian/little-endian in nexrad reader (PHI and ZDR) ({issue}`230`) by [@syedhamidali](https://github.com/syedhamidali), ({pull}`231`) by [@kmuehlbauer](https://github.com/kmuehlbauer). +* ENH: Refactoring all xradar backends to use `from_dict` datatree constructor. Test for `_get_required_root`, `_get_subgroup`, and `_get_radar_calibration` were also added ({pull}`221`) by [@aladinor](https://github.com/aladinor) * ENH: Added pytests to the missing functions in the `test_xradar` and `test_iris` in order to increase codecov in ({pull}`228`) by [@syedhamidali](https://github.com/syedhamidali). * ENH: Updated Readme ({pull}`226`) by [@syedhamidali](https://github.com/syedhamidali). * ADD: Added new module `transform` for transforming CF1 data to CF2 and vice versa ({pull}`224`) by [@syedhamidali](https://github.com/syedhamidali). * Use DataTree from xarray and add xarray nightly run ({pull}`213`, {pull}`214`, {pull}`215`, {pull}`218`) by [@kmuehlbauer](https://github.com/kmuehlbauer). * ADD: Added new accessor `map_over_sweeps` for volume operations on DataTrees and a matching decorator ({pull}`203`) by [@syedhamidali](https://github.com/syedhamidali). - ## 0.7.0 (2024-10-26) This is the last version which uses datatree from xarray-contrib/datatree. Thus, xarray is pinned to version 2024.9.0. diff --git a/tests/io/test_io.py b/tests/io/test_io.py index 3e295f42..1b7c871c 100644 --- a/tests/io/test_io.py +++ b/tests/io/test_io.py @@ -109,7 +109,7 @@ def test_open_odim_datatree_sweep(odim_file, sweep): lswp = len([sweep]) else: lswp = len(sweep) - assert len(dtree.groups[1:]) == lswp + assert len(dtree.match("sweep_*")) == lswp def test_open_odim_datatree(odim_file): @@ -164,7 +164,7 @@ def test_open_odim_datatree(odim_file): 200, 200, ] - for i, grp in enumerate(dtree.groups[1:]): + for i, grp in enumerate(dtree.match("sweep_*")): ds = dtree[grp].ds assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]} assert set(ds.data_vars) & ( @@ -183,7 +183,7 @@ def test_open_odim_datatree(odim_file): "range", } assert np.round(ds.elevation.mean().values.item(), 1) == elevations[i] - assert ds.sweep_number.values == int(grp[7:]) + assert ds.sweep_number.values == int(grp[6:]) @pytest.mark.parametrize("first_dim", ["auto", "time"]) @@ -258,7 +258,7 @@ def test_open_gamic_datatree_sweep(gamic_file, sweep): lswp = len([sweep]) else: lswp = len(sweep) - assert len(dtree.groups[1:]) == lswp + assert len(dtree.match("sweep_*")) == lswp def test_open_gamic_datatree(gamic_file): @@ -319,7 +319,7 @@ def test_open_gamic_datatree(gamic_file): 1000, 1000, ] - for i, grp in enumerate(dtree.groups[1:]): + for i, grp in enumerate(dtree.match("sweep_*")): ds = dtree[grp].ds assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]} assert set(ds.data_vars) & ( @@ -545,7 +545,7 @@ def test_open_rainbow_datatree(rainbow_file): ] azimuths = [361] * 14 ranges = [400] * 14 - for i, grp in enumerate(dtree.groups[1:]): + for i, grp in enumerate(dtree.match("sweep_*")): ds = dtree[grp].ds assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]} assert set(ds.data_vars) & ( @@ -641,28 +641,27 @@ def test_open_iris_datatree(iris0_file): azimuths = [360] * 10 ranges = [664] * 10 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 + for grp in dtree.match("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 def test_open_iris0_dataset(iris0_file): @@ -879,36 +878,35 @@ def test_open_datamet_datatree(datamet_file): 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 + for grp in dtree.match("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 + assert len(dtree.groups) == 5 # Try to read list of sweeps dtree = open_datamet_datatree(datamet_file, sweep=[1, 2]) - assert len(dtree.groups) == 3 + assert len(dtree.groups) == 6 @pytest.mark.parametrize("first_dim", ["time", "auto"]) @@ -993,6 +991,7 @@ def test_cfradial_n_points_file(cfradial1n_file): assert ds.sweep_mode == "azimuth_surveillance" +@pytest.mark.run(order=1) @pytest.mark.parametrize("sweep", ["sweep_0", 0, [0, 1], ["sweep_0", "sweep_1"]]) @pytest.mark.parametrize( "nexradlevel2_files", ["nexradlevel2_gzfile", "nexradlevel2_bzfile"], indirect=True @@ -1003,7 +1002,7 @@ def test_open_nexradlevel2_datatree_sweep(nexradlevel2_files, sweep): lswp = len([sweep]) else: lswp = len(sweep) - assert len(dtree.groups[1:]) == lswp + assert len(dtree.match("sweep*")) == lswp @pytest.mark.parametrize( @@ -1080,8 +1079,8 @@ def test_open_nexradlevel2_datatree(nexradlevel2_files): 308, 232, ] - assert len(dtree.groups[1:]) == 11 - for i, grp in enumerate(dtree.groups[1:]): + assert len(dtree.groups[1:]) == 14 + for i, grp in enumerate(dtree.match("sweep_*")): print(i) ds = dtree[grp].ds assert dict(ds.sizes) == {"azimuth": azimuths[i], "range": ranges[i]} @@ -1101,4 +1100,4 @@ def test_open_nexradlevel2_datatree(nexradlevel2_files): "range", } assert np.round(ds.elevation.mean().values.item(), 1) == elevations[i] - assert ds.sweep_number.values == int(grp[7:]) + assert ds.sweep_number.values == int(grp[6:]) diff --git a/tests/test_util.py b/tests/test_util.py index d33e470d..ed0a80e6 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -11,6 +11,11 @@ import xradar as xd from xradar import io, model, util +from xradar.io.backends.common import ( + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, +) @pytest.fixture( @@ -424,3 +429,57 @@ def dummy_function(ds, refl="none"): sweep_0.dummy_field.values ) # Convert NaNs to zero for comparison assert np.all(np.isclose(non_nan_values, 0)) + + +def test_get_required_root_dataset(): + + filename = DATASETS.fetch("cor-main131125105503.RAW2049") + sweeps = [f"sweep_{i}" for i in range(10)] + ls_ds = [xr.open_dataset(filename, engine="iris", group=sweep) for sweep in sweeps] + root = _get_required_root_dataset(ls_ds, optional=True) + elevations = [ + 0.5, + 1.0, + 2.0, + 3.0, + 5.0, + 7.0, + 10.0, + 15.0, + 20.0, + 30.0, + ] + assert len(root.variables) == 10 + assert root.variables["time_coverage_start"] == "2013-11-25T10:55:04Z" + assert root.variables["time_coverage_end"] == "2013-11-25T10:59:24Z" + np.testing.assert_equal( + root.variables["sweep_fixed_angle"].values, np.array(elevations) + ) + assert len(list(root.attrs.keys())) == 10 + assert root.attrs["instrument_name"] == "Corozal, Radar" + assert root.attrs["scan_name"] == "SURV_HV_300 " + assert root.attrs["comment"] == "AEROCIVIL OPERATIONAL DUAL POLE SCAN" + + +def test_get_radar_calibration(): + filename = DATASETS.fetch("DWD-Vol-2_99999_20180601054047_00.h5") + sweeps = [f"sweep_{i}" for i in range(10)] + ls_ds = [xr.open_dataset(filename, engine="gamic", group=sweep) for sweep in sweeps] + subgroup = _get_radar_calibration(ls_ds, model.radar_calibration_subgroup) + assert len(subgroup.variables) == 6 + assert subgroup["noise_power_h"] == "-3.8298" + assert subgroup["rx_loss_h"] == "3" + assert subgroup["ant_gain_v"] == "43" + assert subgroup["ant_gain_h"] == "43" + + +def test_get_subgroup(): + filename = DATASETS.fetch("71_20181220_060628.pvol.h5") + sweeps = [f"sweep_{i}" for i in range(10)] + ls_ds = [xr.open_dataset(filename, engine="odim", group=sweep) for sweep in sweeps] + subgroup = _get_subgroup(ls_ds, model.radar_parameters_subgroup) + assert len(subgroup.variables) == 3 + assert list(subgroup.variables) == ["longitude", "latitude", "altitude"] + np.testing.assert_almost_equal(subgroup.longitude.values.item(), 151.20899963378906) + np.testing.assert_almost_equal(subgroup.latitude.values.item(), -33.700801849365234) + assert isinstance(subgroup.altitude.values.item(), float) diff --git a/tests/transform/test_cfradial.py b/tests/transform/test_cfradial.py index 5767c5c2..f846b55f 100644 --- a/tests/transform/test_cfradial.py +++ b/tests/transform/test_cfradial.py @@ -2,39 +2,71 @@ # Copyright (c) 2024, openradar developers. # Distributed under the MIT License. See LICENSE for more info. +import pytest import xarray as xr -from open_radar_data import DATASETS +from xarray import MergeError import xradar as xd -def test_to_cfradial1(): +def test_to_cfradial1(cfradial1_file): """Test the conversion from DataTree to CfRadial1 format.""" - file = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc") - dtree = xd.io.open_cfradial1_datatree(file) + with xd.io.open_cfradial1_datatree(cfradial1_file) as dtree: - # Call the conversion function - ds_cf1 = xd.transform.to_cfradial1(dtree) + # Call the conversion function + ds_cf1 = xd.transform.to_cfradial1(dtree) - # Verify key attributes and data structures in the resulting dataset - assert isinstance(ds_cf1, xr.Dataset), "Output is not a valid xarray Dataset" - assert "Conventions" in ds_cf1.attrs and ds_cf1.attrs["Conventions"] == "Cf/Radial" - assert "sweep_mode" in ds_cf1.variables, "Missing sweep_mode in converted dataset" - assert ds_cf1.attrs["version"] == "1.2", "Incorrect CfRadial version" + # Verify key attributes and data structures in the resulting dataset + assert isinstance(ds_cf1, xr.Dataset), "Output is not a valid xarray Dataset" + assert ( + "Conventions" in ds_cf1.attrs and ds_cf1.attrs["Conventions"] == "Cf/Radial" + ) + assert ( + "sweep_mode" in ds_cf1.variables + ), "Missing sweep_mode in converted dataset" + assert ds_cf1.attrs["version"] == "1.2", "Incorrect CfRadial version" -def test_to_cfradial2(): +def test_to_cfradial2(cfradial1_file): """Test the conversion from CfRadial1 to CfRadial2 DataTree format.""" - file = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc") - dtree = xd.io.open_cfradial1_datatree(file) + with xd.io.open_cfradial1_datatree(cfradial1_file) as dtree: - # Convert to CfRadial1 dataset first - ds_cf1 = xd.transform.to_cfradial1(dtree) + # Convert to CfRadial1 dataset first + ds_cf1 = xd.transform.to_cfradial1(dtree) - # Call the conversion back to CfRadial2 - dtree_cf2 = xd.transform.to_cfradial2(ds_cf1) + # Call the conversion back to CfRadial2 + dtree_cf2 = xd.transform.to_cfradial2(ds_cf1) - # Verify key attributes and data structures in the resulting datatree - assert isinstance(dtree_cf2, xr.DataTree), "Output is not a valid DataTree" - assert "radar_parameters" in dtree_cf2, "Missing radar_parameters in DataTree" - assert dtree_cf2.attrs == ds_cf1.attrs, "Attributes mismatch between formats" + # Verify key attributes and data structures in the resulting datatree + assert isinstance(dtree_cf2, xr.DataTree), "Output is not a valid DataTree" + assert "radar_parameters" in dtree_cf2, "Missing radar_parameters in DataTree" + assert dtree_cf2.attrs == ds_cf1.attrs, "Attributes mismatch between formats" + + +def test_to_cfradial1_with_different_range_shapes(nexradlevel2_bzfile): + with xd.io.open_nexradlevel2_datatree(nexradlevel2_bzfile) as dtree: + ds_cf1 = xd.transform.to_cfradial1(dtree) + # Verify key attributes and data structures in the resulting dataset + assert isinstance(ds_cf1, xr.Dataset), "Output is not a valid xarray Dataset" + assert ( + "Conventions" in ds_cf1.attrs and ds_cf1.attrs["Conventions"] == "Cf/Radial" + ) + assert ( + "sweep_mode" in ds_cf1.variables + ), "Missing sweep_mode in converted dataset" + assert ds_cf1.attrs["version"] == "1.2", "Incorrect CfRadial version" + assert ds_cf1.sizes.mapping == {"time": 5400, "range": 1832, "sweep": 11} + + # Call the conversion back to CfRadial2 + dtree_cf2 = xd.transform.to_cfradial2(ds_cf1) + # Verify key attributes and data structures in the resulting datatree + assert isinstance(dtree_cf2, xr.DataTree), "Output is not a valid DataTree" + # todo: this needs to be fixed in nexrad level2reader + # assert "radar_parameters" in dtree_cf2, "Missing radar_parameters in DataTree" + assert dtree_cf2.attrs == ds_cf1.attrs, "Attributes mismatch between formats" + + +def test_to_cfradial1_error_with_different_range_bin_sizes(gamic_file): + with xd.io.open_gamic_datatree(gamic_file) as dtree: + with pytest.raises(MergeError): + xd.transform.to_cfradial1(dtree) diff --git a/xradar/io/backends/cfradial1.py b/xradar/io/backends/cfradial1.py index 7888e80b..7a8a7f51 100644 --- a/xradar/io/backends/cfradial1.py +++ b/xradar/io/backends/cfradial1.py @@ -33,7 +33,7 @@ __doc__ = __doc__.format("\n ".join(__all__)) import numpy as np -from xarray import DataTree, merge, open_dataset +from xarray import Dataset, DataTree, merge, open_dataset from xarray.backends import NetCDF4DataStore from xarray.backends.common import BackendEntrypoint from xarray.backends.store import StoreBackendEntrypoint @@ -49,7 +49,7 @@ required_global_attrs, required_root_vars, ) -from .common import _maybe_decode +from .common import _attach_sweep_groups, _maybe_decode def _get_required_root_dataset(ds, optional=True): @@ -89,7 +89,6 @@ def _get_required_root_dataset(ds, optional=True): root.sweep_group_name.encoding["dtype"] = root.sweep_group_name.dtype # remove cf standard name root.sweep_group_name.attrs = [] - return root @@ -297,6 +296,8 @@ def _get_radar_calibration(ds): subgroup = subgroup.rename_vars(calib_vars) subgroup.attrs = {} return subgroup + else: + return Dataset() def open_cfradial1_datatree(filename_or_obj, **kwargs): @@ -354,24 +355,22 @@ def open_cfradial1_datatree(filename_or_obj, **kwargs): "/georeferencing_correction": _get_subgroup( ds, georeferencing_correction_subgroup ), + "/radar_calibration": _get_radar_calibration(ds), } # radar_calibration (connected with calib-dimension) - calib = _get_radar_calibration(ds) - if calib: - dtree["/radar_calibration"] = calib - - sweep_child = list( - _get_sweep_groups( - ds, - sweep=sweep, - first_dim=first_dim, - optional=optional, - site_coords=site_coords, - ).values() + dtree = _attach_sweep_groups( + dtree, + list( + _get_sweep_groups( + ds, + sweep=sweep, + first_dim=first_dim, + optional=optional, + site_coords=site_coords, + ).values() + ), ) - for i, sw in enumerate(sweep_child): - dtree[f"sweep_{i}"] = sw return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/common.py b/xradar/io/backends/common.py index f055c168..aa91a961 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -20,6 +20,13 @@ import numpy as np import xarray as xr +from ...model import ( + optional_root_attrs, + optional_root_vars, + required_global_attrs, + required_root_vars, +) + def _maybe_decode(attr): try: @@ -56,7 +63,8 @@ def _fix_angle(da): def _attach_sweep_groups(dtree, sweeps): """Attach sweep groups to DataTree.""" for i, sw in enumerate(sweeps): - dtree[f"sweep_{i}"] = xr.DataTree(sw) + # remove attributes only from Dataset's not DataArrays + dtree[f"sweep_{i}"] = xr.DataTree(sw.drop_attrs(deep=False)) return dtree @@ -219,6 +227,80 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): return data +def _get_required_root_dataset(ls_ds, optional=True): + """Extract Root Dataset.""" + # keep only defined mandatory and defined optional variables per default + # by checking in all nodes + data_var = {x for xs in [sweep.variables.keys() for sweep in ls_ds] for x in xs} + remove_root = set(data_var) ^ set(required_root_vars) + if optional: + remove_root ^= set(optional_root_vars) + remove_root ^= {"sweep_number", "fixed_angle"} + remove_root &= data_var + root = [sweep.drop_vars(remove_root) for sweep in ls_ds] + root_vars = {x for xs in [sweep.variables.keys() for sweep in root] for x in xs} + # rename variables + # todo: find a more easy method not iterating over all variables + for k in root_vars: + rename = optional_root_vars.get(k, None) + if rename: + root = [sweep.rename_vars({k: rename}) for sweep in root] + + ds_vars = [sweep[root_vars] for sweep in ls_ds] + _vars = xr.concat(ds_vars, dim="sweep").reset_coords() + + # Creating the root group using _assign_root function + ls = ls_ds.copy() + ls.insert(0, xr.Dataset()) + root = _assign_root(ls) + + # merging both the created and the variables within each dataset + root = xr.merge([root, _vars], compat="override") + + attrs = root.attrs.keys() + remove_attrs = set(attrs) ^ set(required_global_attrs) + if optional: + remove_attrs ^= set(optional_root_attrs) + for k in remove_attrs: + root.attrs.pop(k, None) + # Renaming variable + if "sweep_number" in data_var and "sweep_group_name" not in data_var: + root = root.rename_vars({"sweep_number": "sweep_group_name"}) + elif "sweep_group_name" in data_var: + root["sweep_group_name"].values = np.array( + [f"sweep_{i}" for i in range(len(root["sweep_group_name"].values))] + ) + return root + + +def _get_subgroup(ls_ds: list[xr.Dataset], subdict): + """Get iris-sigmet root metadata group. + Variables are fetched from the provided Dataset according to the subdict dictionary. + """ + meta_vars = subdict + data_vars = {x for xs in [ds.variables.keys() for ds in ls_ds] for x in xs} + extract_vars = set(data_vars) & set(meta_vars) + subgroup = xr.merge([ds[extract_vars] for ds in ls_ds]) + for k in subgroup.data_vars: + rename = meta_vars[k] + if rename: + subgroup = subgroup.rename_vars({k: rename}) + subgroup.attrs = {} + return subgroup + + +def _get_radar_calibration(ls_ds: list[xr.Dataset], subdict: dict) -> xr.Dataset: + """Get radar calibration root metadata group.""" + meta_vars = subdict + data_vars = {x for xs in [ds.attrs for ds in ls_ds] for x in xs} + extract_vars = set(data_vars) & set(meta_vars) + if extract_vars: + var_dict = {var: ls_ds[0].attrs[var] for var in extract_vars} + return xr.Dataset({key: xr.DataArray(value) for key, value in var_dict.items()}) + else: + return xr.Dataset() + + # IRIS Data Types and corresponding python struct format characters # 4.2 Scalar Definitions, Page 23 # https://docs.python.org/3/library/struct.html#format-characters diff --git a/xradar/io/backends/datamet.py b/xradar/io/backends/datamet.py index e5cdd59b..a64ffadd 100644 --- a/xradar/io/backends/datamet.py +++ b/xradar/io/backends/datamet.py @@ -31,6 +31,7 @@ import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -40,6 +41,7 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -48,9 +50,16 @@ get_range_attrs, get_time_attrs, moment_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, sweep_vars_mapping, ) -from .common import _assign_root, _attach_sweep_groups +from .common import ( + _attach_sweep_groups, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, +) #: mapping from DataMet names to CfRadial2/ODIM datamet_mapping = { @@ -361,7 +370,12 @@ def get_variables(self): ) def get_attrs(self): - return FrozenDict() + attributes = { + "scan_name": self.root.scan_metadata["scan_type"], + "instrument_name": self.root.scan_metadata["origin"], + "source": "Datamet", + } + return FrozenDict(attributes) class DataMetBackendEntrypoint(BackendEntrypoint): @@ -476,7 +490,7 @@ def open_datamet_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = kwargs.pop("optional", True) kwargs["backend_kwargs"] = backend_kwargs sweep = kwargs.pop("sweep", None) @@ -499,16 +513,20 @@ def open_datamet_datatree(filename_or_obj, **kwargs): f"sweep_{i}" for i in range(0, dmet.scan_metadata["elevation_number"]) ] - ds = [ + ls_ds: list[xr.Dataset] = [ xr.open_dataset( filename_or_obj, group=swp, engine=DataMetBackendEntrypoint, **kwargs ) for swp in sweeps ] - ds.insert(0, xr.Dataset()) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/furuno.py b/xradar/io/backends/furuno.py index 0ddfd7af..8e34aead 100644 --- a/xradar/io/backends/furuno.py +++ b/xradar/io/backends/furuno.py @@ -47,6 +47,7 @@ import lat_lon_parser import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -56,6 +57,7 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -64,6 +66,8 @@ get_range_attrs, get_time_attrs, moment_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, sweep_vars_mapping, ) from .common import ( @@ -72,10 +76,12 @@ UINT1, UINT2, UINT4, - _assign_root, _attach_sweep_groups, _calculate_angle_res, _get_fmt_string, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, _unpack_dictionary, ) @@ -686,8 +692,13 @@ def get_variables(self): ) def get_attrs(self): - # attributes = {"fixed_angle": float(self.ds.fixed_angle)} - return FrozenDict() + attributes = {"source": "Furuno", "version": self.root.header["format_version"]} + return FrozenDict(attributes) + + def get_calibration_parameters(self): + vars = [var for var in self.root.header if var in radar_calibration_subgroup] + calibration = {var: self.root.header[var] for var in vars} + return FrozenDict(calibration) class FurunoBackendEntrypoint(BackendEntrypoint): @@ -767,7 +778,7 @@ def open_dataset( "altitude": ds.altitude, } ) - + ds.attrs.update(store.get_calibration_parameters()) return ds @@ -806,14 +817,18 @@ def open_furuno_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("optional", True) kwargs["backend_kwargs"] = backend_kwargs - ds = [xr.open_dataset(filename_or_obj, engine="furuno", **kwargs)] - - ds.insert(0, xr.Dataset()) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + ls_ds = [xr.open_dataset(filename_or_obj, engine="furuno", **kwargs)] + + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/gamic.py b/xradar/io/backends/gamic.py index 642ca322..f48a06bc 100644 --- a/xradar/io/backends/gamic.py +++ b/xradar/io/backends/gamic.py @@ -39,6 +39,7 @@ import h5netcdf import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import ( AbstractDataStore, BackendEntrypoint, @@ -53,13 +54,25 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_azimuth_attrs, get_elevation_attrs, get_time_attrs, moment_attrs, + optional_root_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, + required_global_attrs, sweep_vars_mapping, ) -from .common import _assign_root, _attach_sweep_groups, _fix_angle, _get_h5group_names +from .common import ( + _attach_sweep_groups, + _fix_angle, + _get_h5group_names, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, +) from .odim import H5NetCDFArrayWrapper, _get_h5netcdf_encoding, _H5NetCDFMetadata HDF5_LOCK = SerializableLock() @@ -334,7 +347,29 @@ def get_variables(self): ) def get_attrs(self): - return FrozenDict() + _attributes = { + attrs: self.root.grp.attrs[attrs] + for attrs in (dict(self.root.grp.attrs)) + if attrs in required_global_attrs | optional_root_attrs + } + _attributes.update( + { + attrs: self.root.what.attrs[attrs] + for attrs in (dict(self.root.what)) + if attrs in required_global_attrs | optional_root_attrs + } + ) + _attributes["source"] = "gamic" + return FrozenDict(_attributes) + + def get_calibration_parameters(self): + calib_vars = [ + var + for var in dict(self.root.how).keys() + if var in radar_calibration_subgroup + ] + calibration = {var: self.root.how[var] for var in calib_vars} + return FrozenDict(calibration) class GamicBackendEntrypoint(BackendEntrypoint): @@ -445,7 +480,7 @@ def open_dataset( "altitude": ds.altitude, } ) - + ds.attrs.update(store.get_calibration_parameters()) return ds @@ -484,7 +519,7 @@ def open_gamic_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("Optional", True) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -501,14 +536,18 @@ def open_gamic_datatree(filename_or_obj, **kwargs): else: sweeps = _get_h5group_names(filename_or_obj, "gamic") - ds = [ + ls_ds: list[xr.Dataset] = [ xr.open_dataset(filename_or_obj, group=swp, engine="gamic", **kwargs) for swp in sweeps ] - ds.insert(0, xr.open_dataset(filename_or_obj, group="/")) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/hpl.py b/xradar/io/backends/hpl.py index 79849d64..4197a00c 100644 --- a/xradar/io/backends/hpl.py +++ b/xradar/io/backends/hpl.py @@ -37,6 +37,7 @@ import numpy as np import pandas as pd import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -44,14 +45,22 @@ from xarray.core.utils import FrozenDict from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, get_latitude_attrs, get_longitude_attrs, get_time_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, +) +from .common import ( + _attach_sweep_groups, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, ) -from .common import _assign_root, _attach_sweep_groups variable_attr_dict = {} variable_attr_dict["intensity"] = { @@ -571,6 +580,11 @@ def open_dataset( return ds +def _get_h5group_names(filename_or_obj): + store = HplStore.open(filename_or_obj) + return [f"sweep_{i}" for i in store.root.data["sweep_number"]] + + def open_hpl_datatree(filename_or_obj, **kwargs): """Open Halo Photonics processed Doppler lidar dataset as :py:class:`xarray.DataTree`. @@ -606,7 +620,7 @@ def open_hpl_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("optional", None) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -621,19 +635,22 @@ def open_hpl_datatree(filename_or_obj, **kwargs): else: sweeps.extend(sweep) else: - sweeps = ["sweep_0"] + sweeps = _get_h5group_names(filename_or_obj) - ds = [ + ls_ds: list[xr.Dataset] = [ xr.open_dataset(filename_or_obj, group=swp, engine="hpl", **kwargs) for swp in sweeps ] - ds.insert(0, xr.Dataset()) # open_dataset(filename_or_obj, group="/")) - - # create datatree root node with required data - root = _assign_root(ds) - root["fixed_angle"] = ("sweep", [x["sweep_fixed_angle"].values for x in ds[1:]]) - root["sweep_group_name"] = ("sweep", [x["sweep_group_name"].values for x in ds[1:]]) - dtree = xr.DataTree(dataset=root, name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional).rename( + {"sweep_fixed_angle": "fixed_angle"} + ), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/iris.py b/xradar/io/backends/iris.py index 46a1d5b9..59f0d126 100644 --- a/xradar/io/backends/iris.py +++ b/xradar/io/backends/iris.py @@ -43,6 +43,7 @@ import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -52,6 +53,7 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -59,14 +61,16 @@ get_longitude_attrs, get_range_attrs, moment_attrs, - optional_root_attrs, - optional_root_vars, + radar_calibration_subgroup, radar_parameters_subgroup, - required_global_attrs, - required_root_vars, sweep_vars_mapping, ) -from .common import _assign_root, _attach_sweep_groups +from .common import ( + _attach_sweep_groups, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, +) #: mapping from IRIS names to CfRadial2/ODIM iris_mapping = { @@ -3948,6 +3952,16 @@ def get_attrs(self): attributes.update( {"elevation_lower_limit": ll, "elevation_upper_limit": ul} ) + attributes["source"] = "Sigmet" + attributes["scan_name"] = self.root.product_hdr["product_configuration"][ + "task_name" + ] + attributes["instrument_name"] = self.root.ingest_header["ingest_configuration"][ + "site_name" + ].strip() + attributes["comment"] = self.root.ingest_header["task_configuration"][ + "task_end_info" + ]["task_description"] return FrozenDict(attributes) @@ -3959,65 +3973,6 @@ def _get_iris_group_names(filename): return keys -def _get_required_root_dataset(ls_ds, optional=True): - """Extract Root Dataset.""" - # keep only defined mandatory and defined optional variables per default - data_var = {x for xs in [sweep.variables.keys() for sweep in ls_ds] for x in xs} - remove_root = set(data_var) ^ set(required_root_vars) - if optional: - remove_root ^= set(optional_root_vars) - remove_root ^= { - "fixed_angle", - "sweep_group_name", - "sweep_fixed_angle", - } - remove_root &= data_var - root = [sweep.drop_vars(remove_root) for sweep in ls_ds] - root_vars = {x for xs in [sweep.variables.keys() for sweep in root] for x in xs} - # rename variables - # todo: find a more easy method not iterating over all variables - for k in root_vars: - rename = optional_root_vars.get(k, None) - if rename: - root = [sweep.rename_vars({k: rename}) for sweep in root] - - root_vars = {x for xs in [sweep.variables.keys() for sweep in root] for x in xs} - ds_vars = [sweep[root_vars] for sweep in ls_ds] - - root = xr.concat(ds_vars, dim="sweep").reset_coords() - # keep only defined mandatory and defined optional attributes per default - attrs = root.attrs.keys() - remove_attrs = set(attrs) ^ set(required_global_attrs) - if optional: - remove_attrs ^= set(optional_root_attrs) - for k in remove_attrs: - root.attrs.pop(k, None) - # creating a copy of the dataset list for using the _assing_root function. - # and get the variabes/attributes for the root dataset - ls = ls_ds.copy() - ls.insert(0, xr.Dataset()) - dtree = xr.DataTree(dataset=_assign_root(ls), name="root") - root = root.assign(dtree.variables) - root.attrs = dtree.attrs - return root - - -def _get_subgroup(ls_ds: list[xr.Dataset], subdict): - """Get iris-sigmet root metadata group. - Variables are fetched from the provided Dataset according to the subdict dictionary. - """ - meta_vars = subdict - data_vars = {x for xs in [ds.variables.keys() for ds in ls_ds] for x in xs} - extract_vars = set(data_vars) & set(meta_vars) - subgroup = xr.concat([ds[extract_vars] for ds in ls_ds], "sweep") - for k in subgroup.data_vars: - rename = meta_vars[k] - if rename: - subgroup = subgroup.rename_vars({k: rename}) - subgroup.attrs = {} - return subgroup - - class IrisBackendEntrypoint(BackendEntrypoint): """Xarray BackendEntrypoint for IRIS/Sigmet data.""" @@ -4138,7 +4093,7 @@ def open_iris_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = kwargs.pop("optional", True) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -4159,13 +4114,13 @@ def open_iris_datatree(filename_or_obj, **kwargs): xr.open_dataset(filename_or_obj, group=swp, engine="iris", **kwargs) for swp in sweeps ] - # get the datatree root - root = _get_required_root_dataset(ls_ds) - # create datatree root node with required data - dtree = xr.DataTree(dataset=root, name="root") - # get radar_parameters group - subgroup = _get_subgroup(ls_ds, radar_parameters_subgroup) - # attach radar_parameter group - dtree["radar_parameters"] = xr.DataTree(subgroup) - # return Datatree attaching the sweep child nodes - return _attach_sweep_groups(dtree, ls_ds) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/metek.py b/xradar/io/backends/metek.py index dda81965..09ade801 100644 --- a/xradar/io/backends/metek.py +++ b/xradar/io/backends/metek.py @@ -25,6 +25,7 @@ import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -32,12 +33,21 @@ from xarray.core.utils import FrozenDict from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, get_latitude_attrs, get_longitude_attrs, get_time_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, +) +from .common import ( + _attach_sweep_groups, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, ) __all__ = [ @@ -236,6 +246,8 @@ def open(self, filename_or_obj): temp_number = np.zeros((self.n_gates, 64)) spec_var = "" for file_line in self._fp: + if isinstance(file_line, bytes): + file_line = file_line.decode("utf-8") if file_line[:3] == "MRR": if num_times > 0: self._data[spec_var].append(temp_spectra) @@ -651,7 +663,7 @@ def open_metek_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("optional", True) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -668,14 +680,17 @@ def open_metek_datatree(filename_or_obj, **kwargs): else: sweeps = ["sweep_0"] - dtree = {"/": xr.Dataset()} - dtree.update( - { - swp: xr.open_dataset( - filename_or_obj, group=swp, engine="metek", **kwargs - ).copy() - for swp in sweeps - } - ) - - return xr.DataTree.from_dict(dtree) + ls_ds: list[xr.Dataset] = [ + xr.open_dataset(filename_or_obj, group=swp, engine="metek", **kwargs) + for swp in sweeps + ].copy() + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 394b6e83..06290b73 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -41,6 +41,7 @@ import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -52,8 +53,11 @@ from xradar.io.backends.common import ( _assign_root, _attach_sweep_groups, + _get_radar_calibration, + _get_subgroup, ) from xradar.model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -62,6 +66,8 @@ get_range_attrs, get_time_attrs, moment_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, sweep_vars_mapping, ) @@ -594,16 +600,15 @@ def get_data(self, sweep_number, moment=None): ngates = moments[name]["ngates"] word_size = moments[name]["word_size"] data_offset = moments[name]["data_offset"] - ws = {8: 1, 16: 2} - width = ws[word_size] + width = {8: 1, 16: 2}[word_size] data = [] self.rh.pos += data_offset - data.append(self._rh.read(ngates, width=width).view(f"uint{word_size}")) + data.append(self._rh.read(ngates, width=width).view(f">u{width}")) while self.init_next_record() and self.record_number <= stop: if self.record_number in intermediate_records: continue self.rh.pos += data_offset - data.append(self._rh.read(ngates, width=width).view(f"uint{word_size}")) + data.append(self._rh.read(ngates, width=width).view(f">u{width}")) moments[name].update(data=data) def get_data_header(self): @@ -1241,9 +1246,9 @@ def __init__(self, datastore, name, var): - len(datastore.ds["intermediate_records"]) ) nbins = max([v["ngates"] for k, v in datastore.ds["sweep_data"].items()]) - self.dtype = np.dtype("uint8") - if name == "PHI": - self.dtype = np.dtype("uint16") + word_size = datastore.ds["sweep_data"][name]["word_size"] + width = {8: 1, 16: 2}[word_size] + self.dtype = np.dtype(f">u{width}") self.shape = (nrays, nbins) def _getitem(self, key): @@ -1253,8 +1258,14 @@ def _getitem(self, key): except KeyError: self.datastore.root.get_data(self.group, self.name) data = self.datastore.ds["sweep_data"][self.name]["data"] - if self.name == "PHI": + # see 3.2.4.17.6 Table XVII-I Data Moment Characteristics and Conversion for Data Names + word_size = self.datastore.ds["sweep_data"][self.name]["word_size"] + if self.name == "PHI" and word_size == 16: + # 10 bit mask, but only for 2 byte data x = np.uint16(0x3FF) + elif self.name == "ZDR" and word_size == 16: + # 11 bit mask, but only for 2 byte data + x = np.uint16(0x7FF) else: x = np.uint8(0xFF) if len(data[0]) < self.shape[1]: @@ -1537,18 +1548,22 @@ def open_nexradlevel2_datatree(filename_or_obj, **kwargs): engine = NexradLevel2BackendEntrypoint # todo: only open once! Needs new xarray built in datatree! - ds = [] + ls_ds: list[xr.Dataset] = [] for swp in sweeps: try: dsx = xr.open_dataset(filename_or_obj, group=swp, engine=engine, **kwargs) except IndexError: break else: - ds.append(dsx) - - ds.insert(0, xr.Dataset()) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + ls_ds.append(dsx) + ls_ds.insert(0, xr.Dataset()) + dtree: dict = { + "/": _assign_root(ls_ds), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds[1:]) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/odim.py b/xradar/io/backends/odim.py index 6d85a2c0..69528f57 100644 --- a/xradar/io/backends/odim.py +++ b/xradar/io/backends/odim.py @@ -39,6 +39,7 @@ import h5netcdf import numpy as np import xarray as xr +from xarray import DataTree from xarray.backends.common import ( AbstractDataStore, BackendArray, @@ -54,6 +55,7 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -62,13 +64,17 @@ get_range_attrs, get_time_attrs, moment_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, sweep_vars_mapping, ) from .common import ( - _assign_root, _attach_sweep_groups, _fix_angle, _get_h5group_names, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, _maybe_decode, ) @@ -742,7 +748,8 @@ def get_variables(self): ) def get_attrs(self): - return FrozenDict() + attributes = {"Conventions": "ODIM_H5/V2_2"} + return FrozenDict(attributes) class OdimBackendEntrypoint(BackendEntrypoint): @@ -890,7 +897,7 @@ def open_odim_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("optional", True) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -907,15 +914,18 @@ def open_odim_datatree(filename_or_obj, **kwargs): else: sweeps = _get_h5group_names(filename_or_obj, "odim") - ds = [ + ls_ds: list[xr.Dataset] = [ xr.open_dataset(filename_or_obj, group=swp, engine="odim", **kwargs) for swp in sweeps ] - # todo: apply CfRadial2 group structure below - ds.insert(0, xr.open_dataset(filename_or_obj, group="/")) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/backends/rainbow.py b/xradar/io/backends/rainbow.py index 3623914a..8067e73f 100644 --- a/xradar/io/backends/rainbow.py +++ b/xradar/io/backends/rainbow.py @@ -39,6 +39,7 @@ import numpy as np import xarray as xr import xmltodict +from xarray import DataTree from xarray.backends.common import AbstractDataStore, BackendArray, BackendEntrypoint from xarray.backends.file_manager import CachingFileManager from xarray.backends.store import StoreBackendEntrypoint @@ -48,6 +49,7 @@ from ... import util from ...model import ( + georeferencing_correction_subgroup, get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -56,10 +58,16 @@ get_range_attrs, get_time_attrs, moment_attrs, + radar_calibration_subgroup, + radar_parameters_subgroup, sweep_vars_mapping, ) -from .common import _attach_sweep_groups -from .odim import _assign_root +from .common import ( + _attach_sweep_groups, + _get_radar_calibration, + _get_required_root_dataset, + _get_subgroup, +) #: mapping of rainbow moment names to CfRadial2/ODIM names rainbow_mapping = { @@ -904,7 +912,7 @@ def open_rainbow_datatree(filename_or_obj, **kwargs): """ # handle kwargs, extract first_dim backend_kwargs = kwargs.pop("backend_kwargs", {}) - # first_dim = backend_kwargs.pop("first_dim", None) + optional = backend_kwargs.pop("optional", True) sweep = kwargs.pop("sweep", None) sweeps = [] kwargs["backend_kwargs"] = backend_kwargs @@ -921,14 +929,18 @@ def open_rainbow_datatree(filename_or_obj, **kwargs): else: sweeps = _get_rainbow_group_names(filename_or_obj) - ds = [ + ls_ds: list[xr.Dataset] = [ xr.open_dataset(filename_or_obj, group=swp, engine="rainbow", **kwargs) for swp in sweeps ] - ds.insert(0, xr.Dataset()) # open_dataset(filename_or_obj, group="/")) - - # create datatree root node with required data - dtree = xr.DataTree(dataset=_assign_root(ds), name="root") - # return datatree with attached sweep child nodes - return _attach_sweep_groups(dtree, ds[1:]) + dtree: dict = { + "/": _get_required_root_dataset(ls_ds, optional=optional), + "/radar_parameters": _get_subgroup(ls_ds, radar_parameters_subgroup), + "/georeferencing_correction": _get_subgroup( + ls_ds, georeferencing_correction_subgroup + ), + "/radar_calibration": _get_radar_calibration(ls_ds, radar_calibration_subgroup), + } + dtree = _attach_sweep_groups(dtree, ls_ds) + return DataTree.from_dict(dtree) diff --git a/xradar/io/export/cfradial1.py b/xradar/io/export/cfradial1.py index 64d6f32c..8e5cfb65 100644 --- a/xradar/io/export/cfradial1.py +++ b/xradar/io/export/cfradial1.py @@ -57,9 +57,8 @@ def _calib_mapper(calib_params): attrs=data_array.attrs, ) radar_calib_renamed = xr.Dataset(new_data_vars) - dummy_ds = radar_calib_renamed.rename_vars({"r_calib": "fake_coord"}) - del dummy_ds["fake_coord"] - return dummy_ds + radar_calib_renamed = radar_calib_renamed.drop_vars("r_calib", errors="ignore") + return radar_calib_renamed def _main_info_mapper(dtree): @@ -135,12 +134,15 @@ def _variable_mapper(dtree, dim0=None): # Convert to a dataset and append to the list sweep_datasets.append(data) - result_dataset = xr.concat( + # need to use combine_by_coords to correctly test for + # incompatible attrs on DataArray's + result_dataset = xr.combine_by_coords( sweep_datasets, - dim="time", + data_vars="all", compat="no_conflicts", - join="right", - combine_attrs="drop_conflicts", + join="outer", + coords="minimal", + combine_attrs="no_conflicts", ) drop_variables = [ @@ -304,11 +306,11 @@ def to_cfradial1(dtree=None, filename=None, calibs=True): # Add additional parameters if they exist in dtree if "radar_parameters" in dtree: - radar_params = dtree["radar_parameters"].to_dataset() + radar_params = dtree["radar_parameters"].to_dataset().reset_coords() dataset.update(radar_params) if "georeferencing_correction" in dtree: - radar_georef = dtree["georeferencing_correction"].to_dataset() + radar_georef = dtree["georeferencing_correction"].to_dataset().reset_coords() dataset.update(radar_georef) # Ensure that the data type of sweep_mode and similar variables matches diff --git a/xradar/io/export/odim.py b/xradar/io/export/odim.py index a0a47d3d..549fadff 100644 --- a/xradar/io/export/odim.py +++ b/xradar/io/export/odim.py @@ -185,7 +185,7 @@ def to_odim( h5_how = h5.create_group("how") _write_odim(how, h5_how) - grps = dtree.groups[1:] + grps = list(dtree.match("sweep_*")) # what group, object, version, date, time, source, mandatory # p. 10 f @@ -216,7 +216,7 @@ def to_odim( # datasets ds_list = [f"dataset{i + 1}" for i in range(len(grps))] for idx in range(len(ds_list)): - ds = dtree[grps[idx]].ds + ds = dtree[grps[idx]].to_dataset() dim0 = "elevation" if ds.sweep_mode == "rhi" else "azimuth" # datasetN group diff --git a/xradar/model.py b/xradar/model.py index c925afc8..47a0c2fb 100644 --- a/xradar/model.py +++ b/xradar/model.py @@ -228,6 +228,8 @@ ("radar_antenna_gain_v", None), ("radar_beam_width_h", None), ("radar_beam_width_v", None), + ("half_power_beam_width_h", "radar_beam_width_h"), + ("half_power_beam_width_v", "radar_beam_width_v"), ("radar_receiver_bandwidth", None), # cfradial2.1 ("radar_rx_bandwidth", "radar_receiver_bandwidth"), # cfradial1.X ] @@ -241,8 +243,12 @@ ("pulse_width", None), ("antenna_gain_h", None), ("antenna_gain_v", None), + ("ant_gain_h", "antenna_gain_h"), # gamic + ("ant_gain_v", "antenna_gain_v"), # gamic ("xmit_power_h", None), ("xmit_power_v", None), + ("tx_power_h", "xmit_power_h"), + ("tx_power_v", "xmit_power_v"), ("two_way_waveguide_loss_h", None), ("two_way_waveguide_loss_v", None), ("two_way_radome_loss_h", None), @@ -250,6 +256,8 @@ ("receiver_mismatch_loss", None), ("receiver_mismatch_loss_h", None), ("receiver_mismatch_loss_v", None), + ("rx_loss_h", "receiver_mismatch_loss_h"), + ("rx_loss_v", "receiver_mismatch_loss_v"), ("radar_constant_h", None), ("radar_constant_v", None), ("probert_jones_correction", None), @@ -276,6 +284,10 @@ ("sun_power_vx", None), ("noise_source_power_h", None), ("noise_source_power_v", None), + ("noise_power_short_pulse_h", "noise_source_power_h"), + ("noise_power_short_pulse_v", "noise_source_power_v"), + ("noise_power_h", "noise_source_power_h"), # Gamic + ("noise_power_v", "noise_source_power_v"), ("power_measure_loss_h", None), ("power_measure_loss_v", None), ("coupler_forward_loss_h", None), diff --git a/xradar/transform/cfradial.py b/xradar/transform/cfradial.py index c88914d1..d1cde26a 100644 --- a/xradar/transform/cfradial.py +++ b/xradar/transform/cfradial.py @@ -106,11 +106,11 @@ def to_cfradial1(dtree=None, filename=None, calibs=True): # Add additional parameters if they exist in dtree if "radar_parameters" in dtree: - radar_params = dtree["radar_parameters"].to_dataset() + radar_params = dtree["radar_parameters"].to_dataset().reset_coords() dataset.update(radar_params) if "georeferencing_correction" in dtree: - radar_georef = dtree["georeferencing_correction"].to_dataset() + radar_georef = dtree["georeferencing_correction"].to_dataset().reset_coords() dataset.update(radar_georef) # Ensure that the data type of sweep_mode and similar variables matches