From 2451dc9d6720f15b254f39eeeef407a96e45bcfc Mon Sep 17 00:00:00 2001
From: Hamid Ali Syed <35923822+syedhamidali@users.noreply.github.com>
Date: Tue, 3 Sep 2024 18:03:09 -0400
Subject: [PATCH] ADD: Create CAPPI from Radar (#1640)

* ADD: New function to create a CAPPI product

* FIX: Fixed some minor errors

* FIX: Citation

* FIX: Added  AMS Glossary citation

* ADD: Added example of PPI vs CAPPI

---------

Co-authored-by: Max Grover <maxi2312@gmail.com>
---
 examples/plotting/plot_cappi.py |  54 ++++++++
 pyart/retrieve/__init__.py      |   1 +
 pyart/retrieve/cappi.py         | 232 ++++++++++++++++++++++++++++++++
 tests/retrieve/test_cappi.py    |  47 +++++++
 4 files changed, 334 insertions(+)
 create mode 100644 examples/plotting/plot_cappi.py
 create mode 100644 pyart/retrieve/cappi.py
 create mode 100644 tests/retrieve/test_cappi.py

diff --git a/examples/plotting/plot_cappi.py b/examples/plotting/plot_cappi.py
new file mode 100644
index 0000000000..7a0b9a17f5
--- /dev/null
+++ b/examples/plotting/plot_cappi.py
@@ -0,0 +1,54 @@
+"""
+====================
+Compare PPI vs CAPPI
+====================
+
+This example demonstrates how to create and compare PPI (Plan Position Indicator)
+and CAPPI (Constant Altitude Plan Position Indicator) plots using radar data.
+
+In this example, we load sample radar data, create a CAPPI at 2,000 meters
+for the 'reflectivity' field, and then plot
+both the PPI and CAPPI for comparison.
+
+"""
+
+print(__doc__)
+
+# Author: Hamid Ali Syed (syed44@purdue.edu)
+# License: BSD 3 clause
+
+import matplotlib.pyplot as plt
+from open_radar_data import DATASETS
+
+import pyart
+
+# Load the sample radar data
+file = DATASETS.fetch("RAW_NA_000_125_20080411190016")
+radar = pyart.io.read(file)
+
+# Apply gate filtering to exclude unwanted data
+gatefilter = pyart.filters.GateFilter(radar)
+gatefilter.exclude_transition()
+
+# Create CAPPI at 2,000 meters for the 'reflectivity' and 'differential_reflectivity' fields
+cappi = pyart.retrieve.create_cappi(
+    radar, fields=["reflectivity"], height=2000, gatefilter=gatefilter
+)
+
+# Create RadarMapDisplay objects for both PPI and CAPPI
+radar_display = pyart.graph.RadarMapDisplay(radar)
+cappi_display = pyart.graph.RadarMapDisplay(cappi)
+
+# Plotting the PPI and CAPPI for comparison
+fig, ax = plt.subplots(1, 2, figsize=(13, 5))
+
+# Plot PPI for 'reflectivity' field
+radar_display.plot_ppi("reflectivity", vmin=-10, vmax=60, ax=ax[0])
+ax[0].set_title("PPI Reflectivity")
+
+# Plot CAPPI for 'reflectivity' field
+cappi_display.plot_ppi("reflectivity", vmin=-10, vmax=60, ax=ax[1])
+ax[1].set_title("CAPPI Reflectivity at 2000 meters")
+
+# Show the plots
+plt.show()
diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py
index 563ced9a2c..217bbddcb2 100644
--- a/pyart/retrieve/__init__.py
+++ b/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("_")]
diff --git a/pyart/retrieve/cappi.py b/pyart/retrieve/cappi.py
new file mode 100644
index 0000000000..8d9ff9ac34
--- /dev/null
+++ b/pyart/retrieve/cappi.py
@@ -0,0 +1,232 @@
+"""
+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,
+    vel_field="velocity",
+    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.
+    vel_field : str, optional
+        The name of the velocity field to be used for determining the Nyquist velocity.
+        Default is 'velocity'.
+    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.
+    Reference: https://glossary.ametsoc.org/wiki/Cappi
+
+    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)
+        try:
+            nyquist = radar.get_nyquist_vel(sweep=sweep)
+            nyquist = np.round(nyquist)
+        except LookupError:
+            print(
+                f"Nyquist velocity unavailable for sweep {sweep}. Estimating using maximum velocity."
+            )
+            nyquist = radar.fields[vel_field]["data"][sweep_slice].max()
+
+        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:
+        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)
+
+        # Retrieve units and handle case where units might be missing
+        units = radar.fields[field].get("units", "").lower()
+
+        # Determine valid_min and valid_max based on units
+        if units == "dbz":
+            valid_min, valid_max = -10, 80
+        elif units in ["m/s", "meters per second"]:
+            valid_min, valid_max = -100, 100
+        elif units == "db":
+            valid_min, valid_max = -7.9, 7.9
+        else:
+            # If units are not found or don't match known types, set default values or skip masking
+            valid_min, valid_max = None, None
+
+        # If valid_min or valid_max are still None, set them to conservative defaults or skip
+        if valid_min is None:
+            print(f"Warning: valid_min not set for {field}, using default of -1e10")
+            valid_min = -1e10  # Conservative default
+        if valid_max is None:
+            print(f"Warning: valid_max not set for {field}, using default of 1e10")
+            valid_max = 1e10  # Conservative default
+
+        # Apply valid_min and valid_max masking
+        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.set_fill_value(radar.fields[field].get("_FillValue", np.nan))
+        CAPPI = np.ma.masked_invalid(CAPPI)
+        CAPPI = np.ma.masked_outside(CAPPI, valid_min, valid_max)
+
+        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": radar.fields[field].get("_FillValue", np.nan),
+        }
+
+    # 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,
+    )
diff --git a/tests/retrieve/test_cappi.py b/tests/retrieve/test_cappi.py
new file mode 100644
index 0000000000..3a56537f1b
--- /dev/null
+++ b/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"