Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ADD: Create CAPPI from Radar #1640

Merged
merged 6 commits into from
Sep 3, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 54 additions & 0 deletions examples/plotting/plot_cappi.py
Original file line number Diff line number Diff line change
@@ -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 ([email protected])
# 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()
1 change: 1 addition & 0 deletions pyart/retrieve/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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("_")]
232 changes: 232 additions & 0 deletions pyart/retrieve/cappi.py
Original file line number Diff line number Diff line change
@@ -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.
syedhamidali marked this conversation as resolved.
Show resolved Hide resolved
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,
)
47 changes: 47 additions & 0 deletions tests/retrieve/test_cappi.py
Original file line number Diff line number Diff line change
@@ -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"