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

@@ -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 ([email protected])
+# 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)

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

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)