Gridding RHI and PPI sweeps of single radar objects separately into xarray.Dataset objects

pyart/map/__init__.py

@@ -11,5 +11,6 @@
 from .grid_mapper import example_roi_func_dist_beam  # noqa
 from .grid_mapper import grid_from_radars  # noqa
 from .grid_mapper import map_to_grid  # noqa
+from .grid_mapper import grid_ppi_sweeps, grid_rhi_sweeps  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]

pyart/map/_gate_to_grid_map.pyx

@@ -134,10 +134,11 @@ cdef class DistBeamRoI(RoIFunction):
             z_offset = self.offsets[i, 0]
             y_offset = self.offsets[i, 1]
             x_offset = self.offsets[i, 2]
-            roi = (self.h_factor[0] * ((z - z_offset) / 20.0) +
-                   sqrt((self.h_factor[1] * (y - y_offset))**2 +
+            roi = (sqrt((self.h_factor[0] * (z - z_offset))**2 +
+                        (self.h_factor[1] * (y - y_offset))**2 +
                         (self.h_factor[2] * (x - x_offset))**2) *
-                   self.beam_factor)
+                   self.beam_factor
+            )
             if roi < self.min_radius:
                 roi = self.min_radius
             if roi < min_roi:

pyart/map/gates_to_grid.py

@@ -35,9 +35,9 @@ def map_gates_to_grid(
     constant_roi=None,
     z_factor=0.05,
     xy_factor=0.02,
-    min_radius=500.0,
+    min_radius=250.0,
     h_factor=(1.0, 1.0, 1.0),
-    nb=1.5,
+    nb=1.0,
     bsp=1.0,
     dist_factor=(1.0, 1.0, 1.0),
     **kwargs
@@ -95,6 +95,19 @@ def map_gates_to_grid(
     if len(radars) == 0:
         raise ValueError("Length of radars tuple cannot be zero")
 
+    # set min_radius depending on whether processing ARM radars
+    try:
+        if "platform_id" in radars[0].metadata.keys():
+            if np.any(
+                [
+                    x in radars[0].metadata["platform_id"].lower()
+                    for x in ["sacr", "sapr"]
+                ]
+            ):
+                min_radius = 100.0
+    except AttributeError:
+        pass
+
     skip_transform = False
     if len(radars) == 1 and grid_origin_alt is None and grid_origin is None:
         skip_transform = True

pyart/map/grid_mapper.py

@@ -8,6 +8,7 @@
 import netCDF4
 import numpy as np
 import scipy.spatial
+import xarray as xr
 
 from ..config import get_fillvalue, get_metadata
 from ..core.grid import Grid
@@ -26,7 +27,7 @@ def grid_from_radars(
     grid_limits,
     gridding_algo="map_gates_to_grid",
     copy_field_dtypes=True,
-    **kwargs
+    **kwargs,
 ):
     """
     Map one or more radars to a Cartesian grid returning a Grid object.
@@ -292,11 +293,11 @@ def map_to_grid(
     constant_roi=None,
     z_factor=0.05,
     xy_factor=0.02,
-    min_radius=500.0,
-    h_factor=1.0,
-    nb=1.5,
+    min_radius=250.0,
+    h_factor=(1.0, 1.0, 1.0),
+    nb=1.0,
     bsp=1.0,
-    **kwargs
+    **kwargs,
 ):
     """
     Map one or more radars to a Cartesian grid.
@@ -394,8 +395,12 @@ def map_to_grid(
     h_factor, nb, bsp, min_radius : float
         Radius of influence parameters for the built in 'dist_beam' function.
         The parameter correspond to the height scaling, virtual beam width,
-        virtual beam spacing, and minimum radius of influence. These
-        parameters are only used when `roi_func` is 'dist_mean'.
+        virtual beam spacing, and minimum radius of influence.
+        NOTE: the default `min_radius` value is smaller for ARM SACR and SAPR
+        radars (those radars are operated at range resolution of 100 m or
+        higher).
+        to reflect their higher resolution relative to precipitation radars..
+        These parameters are only used when `roi_func` is 'dist_mean'.
     copy_field_data : bool
         True to copy the data within the radar fields for faster gridding,
         the dtype for all fields in the grid will be float64. False will not
@@ -436,6 +441,19 @@ def map_to_grid(
     if len(radars) == 0:
         raise ValueError("Length of radars tuple cannot be zero")
 
+    # set min_radius depending on whether processing ARM radars
+    try:
+        if "platform_id" in radars[0].metadata.keys():
+            if np.any(
+                [
+                    x in radars[0].metadata["platform_id"].lower()
+                    for x in ["sacr", "sapr"]
+                ]
+            ):
+                min_radius = 100.0
+    except AttributeError:
+        pass
+
     skip_transform = False
     if len(radars) == 1 and grid_origin_alt is None and grid_origin is None:
         skip_transform = True
@@ -866,3 +884,231 @@ def roi(zg, yg, xg):
         return min(r)
 
     return roi
+
+
+def grid_ppi_sweeps(
+    radar,
+    target_sweeps=None,
+    grid_size=801,
+    grid_limits="auto",
+    max_z=12000.0,
+    el_rounding_frac=0.25,
+    add_grid_altitude=True,
+    **kwargs,
+):
+    """
+    Separately grid PPI sweeps to an X-Y plane considering only horizontal distances
+    in grid RoI and weighting function.
+    Gridding is performed using the `grid_from_radars` method, which receives any
+    additional input parameters.
+    Note that `h_factor` and `dist_factor` should not be included in kwargs (required
+    for valid gridding results)
+
+    Parameters
+    ----------
+    radar : Radar
+        Radar volume containing PPI sweeps.
+    target_sweeps : int or list
+        sweeps to grid. Using all sweeps in `radar` if None.
+    grid_size: int or 2-tuple
+        grid dimension size. Using sizes for the X-Y plane if tuple.
+        This input parameter is ignored if `grid_shape` is given
+        explicitly via kwargs.
+    grid_limits: 3-tuple with 2-tuple elements or 'auto'
+        if 'auto' using the maximum horizontal range rounded up to the nearest kilometer
+        and limiting vertically up to `max_z`.
+    max_z: float
+        maximum height to consider in gridding (only used if `grid_size` is 'auto')
+    el_rounding_frac: float
+        A fraction for rounding the elevation elements. This variables is also used to
+        represent the sweep for altitude estimation.
+    add_grid_altitude: bool
+        adding a sweep-dependent altitude estimate corresponding to the X-Y plane if True.
+        This output field is useful considering the slanted PPI scans.
+
+    Returns
+    -------
+    radar_ds : xarray.Dataset
+        Radar data gridded to the X-Y plane with a third dimension
+        representing the different sweep elevations.
+
+    """
+    if target_sweeps is None:
+        target_sweeps = radar.sweep_number["data"]
+    elif isinstance(target_sweeps, int):
+        target_sweeps = [target_sweeps]
+
+    # Set grid shape
+    if "grid_shape" not in kwargs.keys():
+        if isinstance(grid_size, int):
+            grid_shape = (1, grid_size, grid_size)
+        elif isinstance(grid_size, tuple):
+            grid_shape = (1, *grid_size)
+        else:
+            raise TypeError("'grid_shape' must be of type int or tuple")
+
+    # Set grid limits in 'auto' option
+    if isinstance(grid_limits, str):
+        if grid_limits == "auto":
+            max_xy = np.max(
+                [np.max(radar.get_gate_x_y_z(sweep=sw)[0]) for sw in target_sweeps]
+            )
+            max_xy = np.ceil(max_xy / 1e3) * 1e3
+            grid_limits = ((0.0, max_z), (-max_xy, max_xy), (-max_xy, max_xy))
+        else:
+            raise ValueError(f"Unknown 'grid_limits' processing string {grid_limits}")
+
+    # Calling the gridding method
+    radar_ds = None
+    for sweep in target_sweeps:
+        radar_sw = radar.extract_sweeps([sweep])
+        sweep_grid = grid_from_radars(
+            (radar_sw,),
+            grid_shape=grid_shape,
+            grid_limits=grid_limits,
+            h_factor=(0.0, 1.0, 1.0),
+            dist_factor=(0.0, 1.0, 1.0),
+            **kwargs,
+        )
+
+        # Convert to xarray.Dataset and finalize
+        el_round = (
+            np.mean(radar_sw.elevation["data"]) / el_rounding_frac
+        ).round() * el_rounding_frac
+        radar_ds_tmp = sweep_grid.to_xarray().squeeze()
+        if add_grid_altitude:
+            alt_est = (radar_ds_tmp["x"] ** 2 + radar_ds_tmp["y"] ** 2) ** 0.5 * np.tan(
+                el_round * np.pi / 180.0
+            )
+        radar_ds_tmp["altitude_est"] = xr.DataArray(
+            alt_est,
+            coords=radar_ds_tmp.coords,
+            dims=radar_ds_tmp.dims,
+            attrs={"long_name": "Estimated altitude in PPI scan", "units": "m"},
+        )
+        radar_ds_tmp = radar_ds_tmp.expand_dims(elevation=[el_round])
+        radar_ds_tmp["elevation"].attrs = {
+            "long_name": "Elevation angle",
+            "units": "deg",
+        }
+        if radar_ds is None:
+            radar_ds = radar_ds_tmp
+        else:
+            radar_ds = xr.concat((radar_ds, radar_ds_tmp), dim="elevation")
+
+    return radar_ds
+
+
+def grid_rhi_sweeps(
+    radar,
+    target_sweeps=None,
+    grid_size=801,
+    grid_limits="auto",
+    max_z=12000.0,
+    az_rounding_frac=0.25,
+    **kwargs,
+):
+    """
+    Separately grid RHI sweeps to a Y-Z plane considering only cross-sectional distances
+    in grid RoI and weighting function.
+    Gridding is performed using the `grid_from_radars` method, which receives any
+    additional input parameters.
+    Note that `h_factor` and `dist_factor` should not be included in kwargs (required
+    for valid gridding results)
+
+    Parameters
+    ----------
+    radar : Radar
+        Radar volume containing PPI sweeps.
+    target_sweeps : int or list
+        sweeps to grid. Using all sweeps in `radar` if None.
+    grid_size: int or 2-tuple
+        grid dimension size. Using sizes for the Y-Z plane if tuple.
+        This input parameter is ignored if `grid_shape` is given
+        explicitly via kwargs.
+    max_z: float
+        maximum height in grid (only used if `grid_size` is 'auto').
+    grid_limits: 3-tuple with 2-tuple elements or 'auto'
+        if 'auto' using the maximum horizontal range and limiting vertically up to 12 km.
+    az_rounding_frac: float
+        A fraction for rounding the azimuth elements.
+
+    Returns
+    -------
+    radar_ds : xarray.Dataset
+        Radar data gridded to the Y-Z plane with a third dimension
+        representing the different sweep azimuths.
+
+    """
+    if target_sweeps is None:
+        target_sweeps = radar.sweep_number["data"]
+    elif isinstance(target_sweeps, int):
+        target_sweeps = [target_sweeps]
+
+    # Set grid shape
+    if "grid_shape" not in kwargs.keys():
+        if isinstance(grid_size, int):
+            grid_shape = (grid_size, grid_size, 1)
+        elif isinstance(grid_size, tuple):
+            grid_shape = (*grid_size, 1)
+        else:
+            raise TypeError("'grid_shape' must be of type int or tuple")
+
+    # Set grid limits in 'auto' option
+    if isinstance(grid_limits, str):
+        if grid_limits == "auto":
+            max_xy = np.max(
+                [np.max(radar.get_gate_x_y_z(sweep=sw)[0]) for sw in target_sweeps]
+            )
+            max_xy = np.ceil(max_xy / 1e3) * 1e3
+            grid_limits = ((0.0, max_z), (-max_xy, max_xy), (-max_xy, max_xy))
+        else:
+            raise ValueError(f"Unknown 'grid_limits' processing string {grid_limits}")
+
+    # Calling the gridding method
+    radar_ds = None
+    for sweep in target_sweeps:
+        radar_sw = radar.extract_sweeps([sweep])
+        if (
+            np.max(radar_sw.azimuth["data"]) - np.min(radar_sw.azimuth["data"]) > 180.0
+        ):  # 0 or 180 deg sweep
+            if np.any(radar_sw.azimuth["data"] > 180.0):
+                diff_center = 180.0  # 0 to 360 deg
+            else:
+                diff_center = 0.0  # -180 to 180
+            az_round = (
+                np.abs(
+                    np.mean(
+                        radar_sw.azimuth["data"]
+                        - 360.0 * (radar_sw.azimuth["data"] > diff_center).astype(int)
+                    )
+                    / az_rounding_frac
+                ).round()
+                * az_rounding_frac
+            )
+        else:
+            az_round = (
+                np.mean(radar_sw.azimuth["data"]) / az_rounding_frac
+            ).round() * az_rounding_frac
+        radar_sw.azimuth[
+            "data"
+        ] -= az_round  # centering azimuth values to 0 to maximize grid range
+        sweep_grid = grid_from_radars(
+            (radar_sw,),
+            grid_shape=grid_shape,
+            grid_limits=grid_limits,
+            h_factor=(1.0, 1.0, 0.0),
+            dist_factor=(1.0, 1.0, 0.0),
+            **kwargs,
+        )
+
+        # Convert to xarray.Dataset and finalize
+        radar_ds_tmp = sweep_grid.to_xarray().squeeze()
+        radar_ds_tmp = radar_ds_tmp.expand_dims(azimuth=[az_round])
+        radar_ds_tmp["azimuth"].attrs = {"long_name": "Azimuth angle", "units": "deg"}
+        if radar_ds is None:
+            radar_ds = radar_ds_tmp
+        else:
+            radar_ds = xr.concat((radar_ds, radar_ds_tmp), dim="azimuth")
+
+    return radar_ds

pyart/testing/sample_objects.py

@@ -149,6 +149,23 @@ def make_target_radar():
     return radar
 
 
+def make_target_rhi_radar():
+    """
+    Return an RHI radar with a target like reflectivity field.
+    """
+    radar = radar = make_empty_rhi_radar(50, 180, 1)
+    fields = {"reflectivity": get_metadata("reflectivity")}
+    fdata = np.zeros((180, 50), dtype="float32")
+    fdata[:, 0:10] = 0.0
+    fdata[:, 10:20] = 10.0
+    fdata[:, 20:30] = 20.0
+    fdata[:, 30:40] = 30.0
+    fdata[:, 40:50] = 40.0
+    fields["reflectivity"]["data"] = fdata
+    radar.fields = fields
+    return radar
+
+
 def make_velocity_aliased_radar(alias=True):
     """
     Return a PPI radar with a target like reflectivity field.

tests/map/test_gates_to_grid.py

@@ -136,11 +136,11 @@ def test_map_to_grid_tiny_grid():
     grids = pyart.map.map_gates_to_grid(
         (radar,),
         grid_shape=(1, 1, 1),
-        grid_limits=((-400.0, 400.0), (-900.0, 900.0), (-900, 900)),
+        grid_limits=((0.0, 300.0), (-800.0, 800.0), (-800, 800)),
         fields=["reflectivity"],
     )
     assert grids["reflectivity"].shape == (1, 1, 1)
-    assert abs(np.round(grids["reflectivity"][0]) - 40.0) < 0.01
+    assert abs(np.round(grids["reflectivity"][0]) - 40.0) < 5.0
 
 
 def test_grid_from_radars_gates_to_grid():