ADD: Create CAPPI from Radar

pyart/retrieve/__init__.py

@@ -32,5 +32,6 @@
 from .spectra_calculations import dealias_spectra, spectra_moments  # noqa
 from .vad import vad_browning, vad_michelson  # noqa
 from .cfad import create_cfad  # noqa
+from .cappi import create_cappi  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]

pyart/retrieve/cappi.py

@@ -0,0 +1,195 @@
+"""
+Constant Altitude Plan Position Indicator
+"""
+
+import numpy as np
+from netCDF4 import num2date
+from pandas import to_datetime
+from scipy.interpolate import RectBivariateSpline
+
+from pyart.core import Radar
+
+
+def create_cappi(
+    radar,
+    fields=None,
+    height=2000,
+    gatefilter=None,
+    same_nyquist=True,
+    nyquist_vector_idx=0,
+):
+    """
+    Create a Constant Altitude Plan Position Indicator (CAPPI) from radar data.
+
+    Parameters:
+    -----------
+    radar : Radar
+        Py-ART Radar object containing the radar data.
+    fields : list of str, optional
+        List of radar fields to be used for creating the CAPPI.
+        If None, all available fields will be used (default is None).
+    height : float, optional
+        The altitude at which to create the CAPPI (default is 2000 meters).
+    gatefilter : GateFilter, optional
+        A GateFilter object to apply masking/filtering to the radar data (default is None).
+    same_nyquist : bool, optional
+        Whether to only stack sweeps with the same Nyquist velocity (default is True).
+    nyquist_vector_idx : int, optional
+        Index for the Nyquist velocity vector if `same_nyquist` is True (default is 0).
+
+    Returns:
+    --------
+    Radar
+        A Py-ART Radar object containing the CAPPI at the specified height.
+
+    Notes:
+    --------
+    CAPPI (Constant Altitude Plan Position Indicator) is a radar visualization
+    technique that provides a horizontal view of meteorological data at a fixed altitude.
+    For more information, see: https://en.wikipedia.org/wiki/Constant_altitude_plan_position_indicator
+
+    Author: Hamid Ali Syed (@syedhamidali)
+    """
+
+    if fields is None:
+        fields = list(radar.fields.keys())
+
+    # Initialize the first sweep as the reference
+    first_sweep = 0
+
+    # Initialize containers for the stacked data and nyquist velocities
+    data_stack = []
+    nyquist_stack = []
+
+    # Process each sweep individually
+    for sweep in range(radar.nsweeps):
+        sweep_slice = radar.get_slice(sweep)
+        nyquist = np.round(radar.get_nyquist_vel(sweep=sweep))
+        sweep_data = {}
+
+        for field in fields:
+            data = radar.get_field(sweep, field)
+
+            # Apply gatefilter if provided
+            if gatefilter is not None:
+                data = np.ma.masked_array(
+                    data, gatefilter.gate_excluded[sweep_slice, :]
+                )
+            time = radar.time["data"][sweep_slice]
+            # Extract and sort azimuth angles
+            azimuth = radar.azimuth["data"][sweep_slice]
+            azimuth_sorted_idx = np.argsort(azimuth)
+            azimuth = azimuth[azimuth_sorted_idx]
+            data = data[azimuth_sorted_idx]
+
+            # Store initial lat/lon for reordering
+            if sweep == first_sweep:
+                azimuth_final = azimuth
+                time_final = time
+            else:
+                # Interpolate data for consistent azimuth ordering across sweeps
+                interpolator = RectBivariateSpline(azimuth, radar.range["data"], data)
+                data = interpolator(azimuth_final, radar.range["data"])
+
+            sweep_data[field] = data[np.newaxis, :, :]
+
+        data_stack.append(sweep_data)
+        nyquist_stack.append(nyquist)
+
+    nyquist_stack = np.array(nyquist_stack)
+
+    # Filter for sweeps with similar Nyquist velocities
+    if same_nyquist:
+        nyquist_range = nyquist_stack[nyquist_vector_idx]
+        nyquist_mask = np.abs(nyquist_stack - nyquist_range) <= 1
+        data_stack = [
+            sweep_data for i, sweep_data in enumerate(data_stack) if nyquist_mask[i]
+        ]
+
+    # Generate CAPPI for each field using data_stack
+    fields_data = {}
+    for field in fields:
+        # Collect the 3D grid for this field from the stacked data
+        data_3d = np.concatenate(
+            [sweep_data[field] for sweep_data in data_stack], axis=0
+        )
+
+        # Sort azimuth for all sweeps
+        dim0 = data_3d.shape[1:]
+        azimuths = np.linspace(0, 359, dim0[0])
+        elevation_angles = radar.fixed_angle["data"][: data_3d.shape[0]]
+        ranges = radar.range["data"]
+
+        theta = (450 - azimuths) % 360
+        THETA, PHI, R = np.meshgrid(theta, elevation_angles, ranges)
+        Z = R * np.sin(PHI * np.pi / 180)
+
+        # Extract the data slice corresponding to the requested height
+        height_idx = np.argmin(np.abs(Z - height), axis=0)
+        CAPPI = np.array(
+            [
+                data_3d[height_idx[j, i], j, i]
+                for j in range(dim0[0])
+                for i in range(dim0[1])
+            ]
+        ).reshape(dim0)
+
+        # Apply valid_min and valid_max masking if they exist
+        field_attrs = radar.fields[field].get("attrs", {})
+        valid_min = field_attrs.get("valid_min", None)
+        valid_max = field_attrs.get("valid_max", None)
+        if valid_min is not None:
+            CAPPI = np.ma.masked_less(CAPPI, valid_min)
+        if valid_max is not None:
+            CAPPI = np.ma.masked_greater(CAPPI, valid_max)
+
+        # Convert to masked array with the specified fill value
+        CAPPI = np.ma.masked_array(CAPPI, fill_value=1e20)
+        CAPPI.set_fill_value(-32768)
+
+        fields_data[field] = {
+            "data": CAPPI,
+            "units": radar.fields[field]["units"],
+            "long_name": f"CAPPI {field} at {height} meters",
+            "comment": f"CAPPI {field} calculated at a height of {height} meters",
+            "_FillValue": -9999.0,
+        }
+
+    # Set the elevation to zeros for CAPPI
+    elevation_final = np.zeros(dim0[0], dtype="float32")
+
+    # since we are using the whole volume scan, report mean time
+    try:
+        dtime = to_datetime(
+            num2date(radar.time["data"], radar.time["units"]).astype(str),
+            format="ISO8601",
+        )
+    except ValueError:
+        dtime = to_datetime(
+            num2date(radar.time["data"], radar.time["units"]).astype(str)
+        )
+    dtime = dtime.mean()
+
+    time = radar.time.copy()
+    time["data"] = time_final
+    time["mean"] = dtime
+
+    # Create the Radar object with the new CAPPI data
+    return Radar(
+        time=radar.time.copy(),
+        _range=radar.range.copy(),
+        fields=fields_data,
+        metadata=radar.metadata.copy(),
+        scan_type=radar.scan_type,
+        latitude=radar.latitude.copy(),
+        longitude=radar.longitude.copy(),
+        altitude=radar.altitude.copy(),
+        sweep_number=radar.sweep_number.copy(),
+        sweep_mode=radar.sweep_mode.copy(),
+        fixed_angle=radar.fixed_angle.copy(),
+        sweep_start_ray_index=radar.sweep_start_ray_index.copy(),
+        sweep_end_ray_index=radar.sweep_end_ray_index.copy(),
+        azimuth=radar.azimuth.copy(),
+        elevation={"data": elevation_final},
+        instrument_parameters=radar.instrument_parameters,
+    )

tests/retrieve/test_cappi.py

@@ -0,0 +1,47 @@
+import numpy as np
+from open_radar_data import DATASETS
+
+import pyart
+from pyart.retrieve import create_cappi
+
+
+def test_create_cappi():
+    # Load radar data
+    file = DATASETS.fetch("RAW_NA_000_125_20080411190016")
+    radar = pyart.io.read(file)
+
+    # Create CAPPI at 10000 meters for the 'reflectivity' field
+    cappi = create_cappi(radar, fields=["reflectivity"], height=10000)
+
+    # Retrieve the 'reflectivity' field from the generated CAPPI
+    reflectivity_cappi = cappi.fields["reflectivity"]
+
+    # Test 1: Check the shape of the reflectivity CAPPI data
+    expected_shape = (360, 992)  # As per the sample data provided
+    assert (
+        reflectivity_cappi["data"].shape == expected_shape
+    ), "Shape mismatch in CAPPI data"
+
+    # Test 2: Check the units of the reflectivity CAPPI
+    assert (
+        reflectivity_cappi["units"] == "dBZ"
+    ), "Incorrect units for CAPPI reflectivity"
+
+    # Test 3: Check that the elevation data is correctly set to zero
+    assert np.all(
+        cappi.elevation["data"] == 0
+    ), "Elevation data should be all zeros in CAPPI"
+
+    # Test 4: Verify the fill value
+    assert (
+        reflectivity_cappi["_FillValue"] == -9999.0
+    ), "Incorrect fill value in CAPPI reflectivity"
+
+    # Test 5: Check the long name and comment
+    assert (
+        reflectivity_cappi["long_name"] == "CAPPI reflectivity at 10000 meters"
+    ), "Incorrect long name"
+    assert (
+        reflectivity_cappi["comment"]
+        == "CAPPI reflectivity calculated at a height of 10000 meters"
+    ), "Incorrect comment"