Repository: ARM-DOE/pyart (language: python)

Issue:
Constant Altitude Plan Position Indicator (CAPPI)
It would be convenient to create a griddisplay class in line with the radardisplay class to be able to plot CAPPIs in km units.

In addition it may be convenient to be able to plot gridmapdisplay and radarmapdisplay in km units.



Discussion / hints:
@meteoswiss-mdr I agree completely with the recommendation. Thanks for submitting this.  Adding labels.

So Just checking I understand @meteoswiss-mdr : x and y axes in km instead of lat/lon? 

Yes, that is what I meant

Sent from my iPhone

On 13 Jul 2017, at 18:23, Scott <notifications@github.com<mailto:notifications@github.com>> wrote:


So Just checking I understand @meteoswiss-mdr<https://github.com/meteoswiss-mdr> : x and y axes in km instead of lat/lon?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<https://github.com/ARM-DOE/pyart/issues/597#issuecomment-315129417>, or mute the thread<https://github.com/notifications/unsubscribe-auth/APYptep7uo8CbGSqCF6DlM-0yoGiL-eyks5sNkR-gaJpZM4KAKPN>.

I recently needed to do a CAPPI, so I have coded up a notebook that do just that:
https://github.com/vlouf/radar_cappi

It wouldn't be a lot of work to add this to Py-ART, so I can do this if you want. 

Strictly speaking, no interpolation is made, just looking for the nearest neighbour at each grid point for a given altitude.
My own code for this can be found here:

https://meteoswiss-mdr.github.io/pyrad/_modules/pyrad/graph/plots_vol.html#plot_cappi


I guess this feature hasn't been added to pyart yet? I need some help on this.
@mustafaalicelik Seems this one went under the radar, no pun intended, I'll see if I can get the pyrad implementation provided by meteoswiss into pyart.
Any improvement?
Can’t you just make a single level grid?

-sent from a mobile device-

> On Feb 14, 2021, at 6:53 AM, malicelik <notifications@github.com> wrote:
> 
> ﻿
> Any improvement?
> 
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub, or unsubscribe.

@mustafaalicelik I'm hoping to have a try at this sometime tomorrow or early next week. 
Hi @hamidrizvi ,

Frankly I don't remember well, it was quite some time ago.
1/ It's to convert the azimuth angle (North is 0 degree and going clockwise) to mathematical angle (0 degree is x-axis and going counter-clockwise).
2/ Altitude is in meter I believe, you should use 2000, not 2.
3/ Nearest neighbour.
Has this CAPPI product already been added to pyart?
@gabrielcassol Hello! At the moment, it has not, still on our list for an addition before the next big release.
@gabrielcassol I would suggest having a look at this, https://github.com/syedhamidali/CAPPI-NETCDF/blob/main/simple_cappi.ipynb

You can make the data gridded, then plot it at whatever height you want. Here I have plotted Max-CAPPI in this example.
@syedhamidali - is there any interest in contributing this to the main Py-ART repo? If not, we can draft something from our side and ping you for review.
@mgrover1 Sounds good! I’ll get started on a script for plotting max_cappi with any projection. If that aligns with what you’re thinking, awesome. If not, just let me know what you’d like me to focus on instead.

Fix patch:
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,
+    )


Test patch:
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"