From f43533ccc84ec51c24d6783ae606d948ae064a26 Mon Sep 17 00:00:00 2001 From: Max Grover Date: Wed, 4 Dec 2024 08:17:01 -0600 Subject: [PATCH] DOC: Add cfradial2 read example (#1700) * FIX: Fix the gallery build issue with new sphinx * ENH: Improve compatibility of xradar integration with cfradial2 * ADD: Add cfradial2 example * FIX: Fix change of nbsphinx * FIX: Fix the io portion of the cfradial2 examples * FIX: Fix the io portion of the cfradial2 examples --- examples/io/plot_read_cfradial2.py | 28 ++++++++++++++++ pyart/xradar/accessor.py | 54 ++++++++++++++++++++++++++++-- tests/xradar/test_accessor.py | 21 ++++++++++++ 3 files changed, 101 insertions(+), 2 deletions(-) create mode 100644 examples/io/plot_read_cfradial2.py diff --git a/examples/io/plot_read_cfradial2.py b/examples/io/plot_read_cfradial2.py new file mode 100644 index 0000000000..1113cb0848 --- /dev/null +++ b/examples/io/plot_read_cfradial2.py @@ -0,0 +1,28 @@ +""" +========================================================== +Read and Plot Cfradial2/FM301 data Using Xradar and Py-ART +========================================================== + +An example which uses xradar and Py-ART to read and plot Cfradial2/FM301 data. + +""" + +# Author: Max Grover (mgrover@anl.gov) +# License: BSD 3 clause + + +import xarray as xr +from open_radar_data import DATASETS + +import pyart + +# Locate the test data and read in using xradar +filename = DATASETS.fetch("cfrad2.20080604_002217_000_SPOL_v36_SUR.nc") +tree = xr.open_datatree(filename) + +# Give the tree Py-ART radar methods +radar = tree.pyart.to_radar() + +# Plot the Reflectivity Field (corrected_reflectivity_horizontal) +display = pyart.graph.RadarMapDisplay(radar) +display.plot_ppi("DBZ", cmap="pyart_ChaseSpectral", vmin=-20, vmax=70) diff --git a/pyart/xradar/accessor.py b/pyart/xradar/accessor.py index 09c084fc9b..34ad3ab1f5 100644 --- a/pyart/xradar/accessor.py +++ b/pyart/xradar/accessor.py @@ -25,7 +25,7 @@ from xarray import DataArray, Dataset, concat from xarray.core import utils from xradar.accessors import XradarAccessor -from xradar.util import get_sweep_keys +from xradar.util import apply_to_sweeps, get_sweep_keys from ..config import get_metadata from ..core.transforms import ( @@ -283,7 +283,23 @@ def to_xarray(self): class Xradar: def __init__(self, xradar, default_sweep="sweep_0", scan_type=None): - self.xradar = xradar + # Make sure that first dimension is azimuth + self.xradar = apply_to_sweeps(xradar, ensure_dim) + # Run through the sanity check for latitude/longitude/altitude + for coord in ["latitude", "longitude", "altitude"]: + if coord not in self.xradar: + raise ValueError( + f"{coord} not included in xradar object, cannot georeference" + ) + + # Ensure that lat/lon/alt info is applied across the sweeps + self.xradar = apply_to_sweeps( + self.xradar, + ensure_georeference_variables, + latitude=self.xradar["latitude"], + longitude=self.xradar["longitude"], + altitude=self.xradar["altitude"], + ) self.scan_type = scan_type or "ppi" self.sweep_group_names = get_sweep_keys(self.xradar) self.nsweeps = len(self.sweep_group_names) @@ -626,6 +642,12 @@ def _combine_sweeps(self): # Add in number of gates variable cleaned["ngates"] = ("gates", np.arange(len(cleaned.gates))) + # Ensure latitude/longitude/altitude are length 1 + if cleaned["latitude"].values.shape != (): + cleaned["latitude"] = cleaned.latitude.isel(gates=0) + cleaned["longitude"] = cleaned.longitude.isel(gates=0) + cleaned["altitude"] = cleaned.altitude.isel(gates=0) + # Return the non-missing times, ensuring valid data is returned return cleaned @@ -837,3 +859,31 @@ def to_radar(self, scan_type=None) -> DataTree: """ dt = self.xarray_obj return Xradar(dt, scan_type=scan_type) + + +def ensure_dim(ds, dim="azimuth"): + """ + Ensure the first dimension is a certain coordinate (ex. azimuth) + """ + core_dims = ds.dims + if dim not in core_dims: + if "time" in core_dims: + ds = ds.swap_dims({"time": dim}) + else: + return ValueError( + "Poorly formatted data: time/azimuth not included as core dimension." + ) + return ds + + +def ensure_georeference_variables(ds, latitude, longitude, altitude): + """ + Ensure georeference variables are included in the sweep information + """ + if "latitude" not in ds: + ds["latitude"] = latitude + if "longitude" not in ds: + ds["longitude"] = longitude + if "altitude" not in ds: + ds["altitude"] = altitude + return ds diff --git a/tests/xradar/test_accessor.py b/tests/xradar/test_accessor.py index 1ec9a2af5d..001dbb17c0 100644 --- a/tests/xradar/test_accessor.py +++ b/tests/xradar/test_accessor.py @@ -7,6 +7,7 @@ import pyart filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc") +cfradial2_file = DATASETS.fetch("cfrad2.20080604_002217_000_SPOL_v36_SUR.nc") def test_get_field(filename=filename): @@ -167,3 +168,23 @@ def test_to_xradar_object(): assert reflectivity.shape == (480, 996) assert radar.nsweeps == 9 assert radar.ngates == 996 + + +def test_cfradial2_data(): + """ + Test swapping the dimensions when the first is time + """ + dtree = xr.open_datatree( + cfradial2_file, + ) + + # Before the transformation, time should be in the dimensions + assert "time" in dtree["sweep_0"].dims + + # After the transformation, azimuth should be moved from coord to dim + radar = dtree.pyart.to_radar() + assert "azimuth" in radar["sweep_0"].dims + + # Ensure georeferencing works as expected + x, y, z = radar.get_gate_x_y_z(0) + assert x.shape == (480, 996)