From 934098044b789bd5224c26a9a92d7634576089a5 Mon Sep 17 00:00:00 2001 From: syedhamidali Date: Mon, 2 Sep 2024 03:52:08 -0400 Subject: [PATCH 1/5] ADD: New function to create a CAPPI product --- pyart/retrieve/__init__.py | 1 + pyart/retrieve/cappi.py | 199 +++++++++++++++++++++++++++++++++++ tests/retrieve/test_cappi.py | 47 +++++++++ 3 files changed, 247 insertions(+) create mode 100644 pyart/retrieve/cappi.py create mode 100644 tests/retrieve/test_cappi.py 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..7cff039a60 --- /dev/null +++ b/pyart/retrieve/cappi.py @@ -0,0 +1,199 @@ +""" +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 + fields_data = {} + for field in fields: + # Reshape original data to 3D + original_data = radar.fields[field]["data"] + dim0 = original_data[radar.get_slice(0)].shape + data_3d = original_data.reshape((radar.nsweeps, dim0[0], dim0[1])) + + # Sort azimuth for all sweeps + for ns in range(radar.nsweeps): + sorted_idx = np.argsort(radar.azimuth["data"][radar.get_slice(ns)]) + data_3d[ns] = data_3d[ns, sorted_idx, :] + + # Generate cartesian coordinates from radar spherical coordinates + azimuths = np.linspace(0, 359, dim0[0]) + elevation_angles = radar.fixed_angle["data"] + 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, + ) 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" From af54fe6b59ed630b20928092730d8e479f99abd0 Mon Sep 17 00:00:00 2001 From: syedhamidali Date: Mon, 2 Sep 2024 05:09:14 -0400 Subject: [PATCH 2/5] FIX: Fixed some minor errors --- pyart/retrieve/cappi.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/pyart/retrieve/cappi.py b/pyart/retrieve/cappi.py index 7cff039a60..a9b2d088a9 100644 --- a/pyart/retrieve/cappi.py +++ b/pyart/retrieve/cappi.py @@ -106,22 +106,18 @@ def create_cappi( sweep_data for i, sweep_data in enumerate(data_stack) if nyquist_mask[i] ] - # Generate CAPPI for each field + # Generate CAPPI for each field using data_stack fields_data = {} for field in fields: - # Reshape original data to 3D - original_data = radar.fields[field]["data"] - dim0 = original_data[radar.get_slice(0)].shape - data_3d = original_data.reshape((radar.nsweeps, dim0[0], dim0[1])) + # 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 - for ns in range(radar.nsweeps): - sorted_idx = np.argsort(radar.azimuth["data"][radar.get_slice(ns)]) - data_3d[ns] = data_3d[ns, sorted_idx, :] - - # Generate cartesian coordinates from radar spherical coordinates + dim0 = data_3d.shape[1:] azimuths = np.linspace(0, 359, dim0[0]) - elevation_angles = radar.fixed_angle["data"] + elevation_angles = radar.fixed_angle["data"][: data_3d.shape[0]] ranges = radar.range["data"] theta = (450 - azimuths) % 360 From b241595d629935e5c45b65c9b95a30ff634ac4f8 Mon Sep 17 00:00:00 2001 From: syedhamidali Date: Tue, 3 Sep 2024 13:52:50 -0400 Subject: [PATCH 3/5] FIX: Citation --- pyart/retrieve/cappi.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyart/retrieve/cappi.py b/pyart/retrieve/cappi.py index a9b2d088a9..8eaaee06ee 100644 --- a/pyart/retrieve/cappi.py +++ b/pyart/retrieve/cappi.py @@ -21,7 +21,7 @@ def create_cappi( """ Create a Constant Altitude Plan Position Indicator (CAPPI) from radar data. - Parameters: + Parameters ----------- radar : Radar Py-ART Radar object containing the radar data. @@ -37,16 +37,16 @@ def create_cappi( nyquist_vector_idx : int, optional Index for the Nyquist velocity vector if `same_nyquist` is True (default is 0). - Returns: + Returns -------- Radar A Py-ART Radar object containing the CAPPI at the specified height. - Notes: + 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 + Atlas, David, ed. Radar in Meteorology: Battan Memorial and 40th Anniversary Radar Meteorology Conference. Springer, 2015. Author: Hamid Ali Syed (@syedhamidali) """ From bd6eba5a681ab9cafe3e0e330c3f79cafbf0f4c6 Mon Sep 17 00:00:00 2001 From: syedhamidali Date: Tue, 3 Sep 2024 14:06:13 -0400 Subject: [PATCH 4/5] FIX: Added AMS Glossary citation --- pyart/retrieve/cappi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyart/retrieve/cappi.py b/pyart/retrieve/cappi.py index 8eaaee06ee..7099b78bb1 100644 --- a/pyart/retrieve/cappi.py +++ b/pyart/retrieve/cappi.py @@ -46,9 +46,9 @@ def create_cappi( -------- CAPPI (Constant Altitude Plan Position Indicator) is a radar visualization technique that provides a horizontal view of meteorological data at a fixed altitude. - Atlas, David, ed. Radar in Meteorology: Battan Memorial and 40th Anniversary Radar Meteorology Conference. Springer, 2015. + Reference : https://glossary.ametsoc.org/wiki/Cappi - Author: Hamid Ali Syed (@syedhamidali) + Author : Hamid Ali Syed (@syedhamidali) """ if fields is None: From 62c77fd39cd373039b9e2b25eb1c3ad446edfc91 Mon Sep 17 00:00:00 2001 From: syedhamidali Date: Tue, 3 Sep 2024 16:45:34 -0400 Subject: [PATCH 5/5] ADD: Added example of PPI vs CAPPI --- examples/plotting/plot_cappi.py | 54 +++++++++++++++++++++++ pyart/retrieve/cappi.py | 77 ++++++++++++++++++++++++--------- 2 files changed, 111 insertions(+), 20 deletions(-) create mode 100644 examples/plotting/plot_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/cappi.py b/pyart/retrieve/cappi.py index 7099b78bb1..8d9ff9ac34 100644 --- a/pyart/retrieve/cappi.py +++ b/pyart/retrieve/cappi.py @@ -15,6 +15,7 @@ def create_cappi( fields=None, height=2000, gatefilter=None, + vel_field="velocity", same_nyquist=True, nyquist_vector_idx=0, ): @@ -22,33 +23,41 @@ def create_cappi( 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). + 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). + 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). + 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). + 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). + 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 + Reference: https://glossary.ametsoc.org/wiki/Cappi - Author : Hamid Ali Syed (@syedhamidali) + Author + ------ + Hamid Ali Syed (@syedhamidali) """ if fields is None: @@ -64,7 +73,15 @@ def create_cappi( # 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)) + 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: @@ -76,6 +93,7 @@ def create_cappi( 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) @@ -109,7 +127,6 @@ def create_cappi( # 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 ) @@ -134,31 +151,51 @@ def create_cappi( ] ).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) + # 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 = np.ma.masked_array(CAPPI, fill_value=1e20) - CAPPI.set_fill_value(-32768) + 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": -9999.0, + "_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 + # Since we are using the whole volume scan, report mean time try: dtime = to_datetime( num2date(radar.time["data"], radar.time["units"]).astype(str),