ENH: accessor methods reindex_angle and reindex_range

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)
 

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)
+```

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()

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.