diff --git a/docs/history.md b/docs/history.md index a464149e..da3b489d 100644 --- a/docs/history.md +++ b/docs/history.md @@ -22,6 +22,7 @@ * TST: every sweep of every engine passes ``xradar.util.is_sweep``, and ``is_sweep(strict=True)`` for all engines except HPL and Metek; the ``Open-Any-Format`` notebook shows the check ({issue}`537`, {pull}`542`) by [@syedhamidali](https://github.com/syedhamidali) * FIX: the calibration values of Furuno SCNX files are scaled (the header gives antenna gain, transmit power, radar constant and noise power in tenths) and have ``units``; the transmit power is ``xmit_power_h/v`` in dBm and is written to CfRadial1 files ({issue}`538`, {pull}`543`) by [@syedhamidali](https://github.com/syedhamidali) * FIX: the CfRadial1 reader gave each sweep the rays of another sweep when the sweeps of the file were not in time order (e.g. written by ``to_cfradial1`` from ``sweep=[5, 0]``); the sweeps are sliced in file order and sorted by time within each sweep ({issue}`550`, {pull}`551`) by [@syedhamidali](https://github.com/syedhamidali) +* ENH: ``.xradar.reindex_angle`` and ``.xradar.reindex_range`` for Dataset (a sweep) and DataTree (every sweep), the same arguments as ``xradar.util.reindex_angle`` and ``reindex_range`` ({issue}`539`, {pull}`544`) by [@syedhamidali](https://github.com/syedhamidali) ## 0.13.0 (2026-10-21) diff --git a/docs/notebooks/angle_reindexing.md b/docs/notebooks/angle_reindexing.md index d6d61b2e..18df5658 100644 --- a/docs/notebooks/angle_reindexing.md +++ b/docs/notebooks/angle_reindexing.md @@ -145,3 +145,59 @@ plt.gca().grid() ```{code-cell} ds_out2.DBTH.sortby("time").plot(y="time") ``` + +## Reindexing after reading: the accessor methods + +`reindex_coord` reindexes while reading. Data which is already open can be reindexed with the `.xradar` accessor, on a sweep (Dataset) and on a whole volume (DataTree). The methods take the arguments of {func}`xradar.util.reindex_angle` and {func}`xradar.util.reindex_range`. + +### Angle of a sweep + +The same as `fix_angle` and `ipol_time` above, in one call: + +```{code-cell} +params = xd.util.extract_angle_parameters(ds_in) +ds_acc = ds_in.xradar.reindex_angle( + start_angle=params["start_angle"], + stop_angle=params["stop_angle"], + angle_res=params["angle_res"], + direction=params["direction"], +) +xr.testing.assert_identical(ds_acc, ds_out2) +ds_acc.sizes +``` + +### Range of all sweeps of a volume + +The sweeps of this volume have different gate spacing, 100 m in the first six sweeps and 150 m in the others, and different lengths. CfRadial1 has one `range` axis, so the export asks for one gate geometry: + +```{code-cell} +import numpy as np + +dtree = xd.open_datatree(filename, engine="gamic") +for name in ("sweep_0", "sweep_5", "sweep_6", "sweep_9"): + rng = dtree[name].dataset.range.values + print(f"{name}: {rng.size:5d} gates, first {rng[0]:5.0f} m, spacing {np.median(np.diff(rng)):4.0f} m") +``` + +```{code-cell} +try: + dtree.xradar.to_cfradial1() +except ValueError as err: + print(err) +``` + +`reindex_range` on the DataTree puts every sweep on one grid. Each gate of the grid takes the nearest gate of the sweep (up to half a grid spacing), so this is a choice of the resolution: here 300 m, which both spacings fit, from 150 m to 36150 m. Gates beyond the last gate of the grid are not kept. + +```{code-cell} +common = dtree.xradar.reindex_range(start_range=150.0, stop_range=36150.0, range_res=300.0) +ranges = [common[name].dataset.range.values for name in ("sweep_0", "sweep_9")] +assert np.array_equal(*ranges) +print(f"{ranges[0].size} gates from {ranges[0][0]:.0f} m to {ranges[0][-1]:.0f} m") +``` + +Now the volume can be exported, and the sweeps can be stacked (see the Multi-Volume-Concatenation notebook). + +```{code-cell} +cf1 = common.xradar.to_cfradial1() +dict(cf1.sizes) +``` diff --git a/tests/test_accessors.py b/tests/test_accessors.py index 70abd2fa..8688f188 100644 --- a/tests/test_accessors.py +++ b/tests/test_accessors.py @@ -272,3 +272,79 @@ def test_create_xradar_dataarray_accessor(): ) xr.testing.assert_equal(da.xradar_test.double(), da * 2) assert float(da.xradar_test.total()) == 6.0 + + +def _jittered_sweep(azimuth, range_start, ngates=6, res=250.0): + # a small PPI whose angles and ranges differ slightly between sweeps + rng = (range_start + np.arange(ngates) * res).astype("float32") + data = np.arange(4 * ngates, dtype="float32").reshape(4, ngates) + return xr.Dataset( + {"DBZH": (("azimuth", "range"), data, {"_FillValue": np.float32(np.nan)})}, + coords={ + "azimuth": ("azimuth", np.array(azimuth, dtype="float32")), + "elevation": ("azimuth", np.full(4, 0.5, dtype="float32")), + "time": ("azimuth", np.datetime64("2024-01-01", "s") + np.arange(4)), + "range": ("range", rng, {"units": "meters"}), + }, + ) + + +def test_accessor_reindex_dataset(): + ds = _jittered_sweep([44.0, 136.0, 224.0, 316.0], 125.4) + # the same as the functions + angle = dict(start_angle=0, stop_angle=360, angle_res=90.0, direction=1) + rng = dict(start_range=125.0, stop_range=1375.0, range_res=250.0) + xr.testing.assert_identical( + ds.xradar.reindex_angle(**angle), xd.util.reindex_angle(ds, **angle) + ) + xr.testing.assert_identical( + ds.xradar.reindex_range(**rng), xd.util.reindex_range(ds, **rng) + ) + out = ds.xradar.reindex_angle(**angle).xradar.reindex_range(**rng) + np.testing.assert_array_equal(out.azimuth, [45.0, 135.0, 225.0, 315.0]) + np.testing.assert_array_equal(out.range, 125.0 + np.arange(6) * 250.0) + + +def test_accessor_reindex_datatree(): + sweeps = { + "/sweep_0": _jittered_sweep([44.0, 136.0, 224.0, 316.0], 125.4), + "/sweep_1": _jittered_sweep([46.0, 134.0, 226.0, 314.0], 124.6), + } + dtree = xr.DataTree.from_dict({"/": xr.Dataset(attrs={"title": "t"}), **sweeps}) + angle = dict(start_angle=0, stop_angle=360, angle_res=90.0, direction=1) + rng = dict(start_range=125.0, stop_range=1375.0, range_res=250.0) + out = dtree.xradar.reindex_angle(**angle).xradar.reindex_range(**rng) + # every sweep is on the same grid, so that the sweeps can be stacked + for name in ("sweep_0", "sweep_1"): + np.testing.assert_array_equal(out[name].azimuth, [45.0, 135.0, 225.0, 315.0]) + np.testing.assert_array_equal(out[name].range, 125.0 + np.arange(6) * 250.0) + stacked = xr.concat( + [out["sweep_0"].to_dataset(), out["sweep_1"].to_dataset()], dim="volume_time" + ) + assert stacked.sizes == {"volume_time": 2, "azimuth": 4, "range": 6} + assert not np.isnan(stacked.DBZH).any() + # the other nodes are not touched, the input is not changed + assert out.attrs["title"] == "t" + assert dtree["sweep_0"].azimuth[0] == 44.0 + # same result as the function on a sweep + expected = xd.util.reindex_range(sweeps["/sweep_0"], **rng) + xr.testing.assert_equal( + out["sweep_0"].to_dataset().xradar.reindex_range(**rng).DBZH, + expected.pipe(xd.util.reindex_angle, **angle).DBZH, + ) + + +def test_accessor_reindex_angle_as_reindex_coord(): + # duplicate rays (sweep_4, sweep_7) and the same result as reading with + # reindex_coord + fname = DATASETS.fetch("DWD-Vol-2_99999_20180601054047_00.h5") + angle = dict(start_angle=0, stop_angle=360, angle_res=1.0, direction=1) + out = xd.open_datatree(fname, engine="gamic").xradar.reindex_angle(**angle) + read = xd.open_datatree(fname, engine="gamic", reindex_coord={"angle": angle}) + for name in ("sweep_0", "sweep_4", "sweep_7"): + xr.testing.assert_identical( + out[name].to_dataset(), read[name].to_dataset(inherit="all_coords") + ) + # missing rays get a time + ds = read["sweep_0"].to_dataset().isel(azimuth=np.r_[0:100, 110:360]) + assert not np.isnat(ds.xradar.reindex_angle(**angle).time).any() diff --git a/xradar/accessors.py b/xradar/accessors.py index 76aa11b2..f0b4dbe5 100644 --- a/xradar/accessors.py +++ b/xradar/accessors.py @@ -29,7 +29,22 @@ from .georeference import add_crs, add_crs_tree, get_crs, get_x_y_z, get_x_y_z_tree from .transform import to_cfradial1, to_cfradial2 -from .util import is_sweep, map_over_sweeps +from .util import ( + apply_to_sweeps, + ipol_time, + is_sweep, + map_over_sweeps, + reindex_angle, + reindex_range, + remove_duplicate_rays, +) + + +def _fix_angle(ds, **kwargs): + # as reindex_coord while reading: remove duplicate rays, fill the time of + # missing rays + ds = reindex_angle(remove_duplicate_rays(ds), **kwargs) + return ipol_time(ds, **kwargs) def accessor_constructor(self, xarray_obj): @@ -182,6 +197,25 @@ def is_sweep(self, strict=False): """ return is_sweep(self.xarray_obj, strict=strict) + def reindex_angle(self, **kwargs): + """Reindex the sweep along its first angle onto a regular grid. + + See :func:`xradar.util.reindex_angle` for the keyword arguments + (``start_angle``, ``stop_angle``, ``angle_res``, ``direction``, ``method`` + and ``tolerance``). As ``reindex_coord`` while reading, duplicate rays are + removed first and the time of missing rays is interpolated. + """ + return _fix_angle(self.xarray_obj, **kwargs) + + def reindex_range(self, **kwargs): + """Reindex the sweep along range onto a regular grid. + + See :func:`xradar.util.reindex_range` for the keyword arguments + (``start_range``, ``stop_range``, ``range_res``, ``method`` and + ``tolerance``). + """ + return reindex_range(self.xarray_obj, **kwargs) + def to_cfradial2_datatree(self): """Convert a CfRadial1 Dataset to CfRadial2 DataTree.""" return to_cfradial2(self.xarray_obj) @@ -245,6 +279,27 @@ def is_sweep(self, strict=False): """ return is_sweep(self.xarray_obj, strict=strict) + def reindex_angle(self, **kwargs): + """Reindex every sweep along its first angle onto a regular grid. + + See :func:`xradar.util.reindex_angle` for the keyword arguments. The + angle grid is the same for all sweeps. As ``reindex_coord`` while reading, + duplicate rays are removed first and the time of missing rays is + interpolated. + """ + return apply_to_sweeps(self.xarray_obj, _fix_angle, **kwargs) + + def reindex_range(self, **kwargs): + """Reindex every sweep along range onto a regular grid. + + See :func:`xradar.util.reindex_range` for the keyword arguments. Without + arguments each sweep keeps its own range grid and only gets a regular one + (the first gate, the last gate and the median spacing of the sweep); give + ``start_range``, ``stop_range`` and ``range_res`` to bring all sweeps onto + the same grid. + """ + return apply_to_sweeps(self.xarray_obj, reindex_range, **kwargs) + def map_over_sweeps(self, func, *args, **kwargs): """ Apply a function across all sweep nodes in the DataTree using the xradar accessor.