ENH: Max CAPPI (#1633)

* ENH: xarray grid compatibility

* ENH: xarray grid output

* adding back xarray availablity check

* pre-commit run

* fix minor error in grid_to_xarary

* modified coords in function to_xarray()

* modified coords in function to_xarray() + pre-commit

* Add plot_max_cappi example and update related modules

* Update tests and fix Example

* drop a few  args

* Fix grid.to_xarray()

* FIX: Fixed some tests and to_xarray() function

* FIX: Fixed to_xarray() for more than one radar

---------

Co-authored-by: Zach Sherman <[email protected]>

examples/plotting/plot_max_cappi.py

@@ -0,0 +1,111 @@
+"""
+=============
+Plot Max-CAPPI
+=============
+
+This is an example of how to plot a Max-CAPPI
+within a Py-ART grid display object.
+
+"""
+
+print(__doc__)
+
+# Author: Hamid Ali Syed ([email protected])
+# License: BSD 3 clause
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pyart
+from pyart.testing import get_test_data
+
+#########################################
+# ** MAX-CAPPI Display
+#
+
+# Define and Read in the test data
+grid_file = get_test_data("20110520100000_nexrad_grid.nc")
+grid = pyart.io.read_grid(grid_file)
+
+
+# Create a grid display
+gdisplay = pyart.graph.GridMapDisplay(grid)
+gdisplay.plot_maxcappi(field="REF", range_rings=True, add_slogan=True)
+
+
+#########################################
+# ** Second Example
+#
+# Let's read in a cfradial file and create a grid.
+
+
+import logging
+from datetime import datetime
+
+import fsspec
+import pytz
+
+
+def download_nexrad(timezone, date, site, local_date=False):
+    """Download NEXRAD radar data from an S3 bucket."""
+    try:
+        utc_date = (
+            pytz.timezone(timezone).localize(date).astimezone(pytz.utc)
+            if local_date
+            else date
+        )
+        logging.info(f"Time: {utc_date}")
+
+        fs = fsspec.filesystem("s3", anon=True)
+        nexrad_path = utc_date.strftime(
+            f"s3://noaa-nexrad-level2/%Y/%m/%d/{site}/{site}%Y%m%d_%H*"
+        )
+        files = sorted(fs.glob(nexrad_path))
+
+        return [file for file in files if not file.endswith("_MDM")]
+    except Exception as e:
+        logging.error("Error in processing: %s", e)
+        return []
+
+
+# Load NEXRAD data from S3 Bucket
+site = "PHWA"
+timezone = "UTC"
+date = datetime(2024, 8, 25, 8, 29)
+
+# Correctly passing the site and timezone
+file = download_nexrad(timezone, date, site, local_date=False)[0]
+
+
+# Read the data using nexrad_archive reader
+radar = pyart.io.read_nexrad_archive("s3://" + file)
+
+# Create a 3D grid
+# Mask out last 10 gates of each ray, this removes the "ring" around the radar.
+radar.fields["reflectivity"]["data"][:, -10:] = np.ma.masked
+
+# Exclude masked gates from the gridding
+gatefilter = pyart.filters.GateFilter(radar)
+gatefilter.exclude_transition()
+gatefilter.exclude_masked("reflectivity")
+gatefilter.exclude_outside("reflectivity", 10, 80)
+
+# Perform Cartesian mapping, limit to the reflectivity field.
+max_range = np.ceil(radar.range["data"].max())
+if max_range / 1e3 > 250:
+    max_range = 250 * 1e3
+
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    gatefilters=(gatefilter,),
+    grid_shape=(30, 441, 441),
+    grid_limits=((0, 10000), (-max_range, max_range), (-max_range, max_range)),
+    fields=["reflectivity"],
+)
+
+# Create a grid display
+gdisplay = pyart.graph.GridMapDisplay(grid)
+with plt.style.context("dark_background"):
+    gdisplay.plot_maxcappi(
+        field="reflectivity", cmap="pyart_HomeyerRainbow", add_slogan=True
+    )

pyart/core/grid.py

@@ -320,68 +320,136 @@ def to_xarray(self):
         x, y, z : dict, 1D
             Distance from the grid origin for each Cartesian coordinate axis
             in a one dimensional array.
-
         """
-
         if not _XARRAY_AVAILABLE:
             raise MissingOptionalDependency(
-                "Xarray is required to use Grid.to_xarray but is not " + "installed!"
+                "Xarray is required to use Grid.to_xarray but is not installed!"
             )
 
+        def _process_radar_name(radar_name):
+            """Process radar_name to handle different formats."""
+            if radar_name.dtype.kind == "S" and radar_name.ndim > 1:
+                # Join each row of bytes into a single byte string
+                return np.array(
+                    [b"".join(row) for row in radar_name],
+                    dtype=f"|S{radar_name.shape[1]}",
+                )
+            return radar_name
+
         lon, lat = self.get_point_longitude_latitude()
         z = self.z["data"]
         y = self.y["data"]
         x = self.x["data"]
 
-        time = np.array([num2date(self.time["data"][0], self.time["units"])])
+        time = np.array([num2date(self.time["data"][0], units=self.time["units"])])
 
         ds = xarray.Dataset()
-        for field in list(self.fields.keys()):
-            field_data = self.fields[field]["data"]
+        for field, field_info in self.fields.items():
+            field_data = field_info["data"]
             data = xarray.DataArray(
                 np.ma.expand_dims(field_data, 0),
                 dims=("time", "z", "y", "x"),
                 coords={
-                    "time": (["time"], time),
-                    "z": (["z"], z),
+                    "time": time,
+                    "z": z,
                     "lat": (["y", "x"], lat),
                     "lon": (["y", "x"], lon),
-                    "y": (["y"], y),
-                    "x": (["x"], x),
+                    "y": y,
+                    "x": x,
                 },
             )
-            for meta in list(self.fields[field].keys()):
-                if meta != "data":
-                    data.attrs.update({meta: self.fields[field][meta]})
-
+            data.attrs.update({k: v for k, v in field_info.items() if k != "data"})
             ds[field] = data
-            ds.lon.attrs = [
-                ("long_name", "longitude of grid cell center"),
-                ("units", "degree_E"),
-                ("standard_name", "Longitude"),
-            ]
-            ds.lat.attrs = [
-                ("long_name", "latitude of grid cell center"),
-                ("units", "degree_N"),
-                ("standard_name", "Latitude"),
-            ]
-
-            ds.z.attrs = get_metadata("z")
-            ds.y.attrs = get_metadata("y")
-            ds.x.attrs = get_metadata("x")
-
-            ds.z.encoding["_FillValue"] = None
-            ds.lat.encoding["_FillValue"] = None
-            ds.lon.encoding["_FillValue"] = None
-
-            # Grab original radar(s) name and number of radars used to make grid
-            ds.attrs["nradar"] = self.nradar
-            ds.attrs["radar_name"] = self.radar_name
-
-            # Grab all metadata
-            ds.attrs.update(self.metadata)
-
-            ds.close()
+
+        ds.lon.attrs = {
+            "long_name": "longitude of grid cell center",
+            "units": "degree_E",
+            "standard_name": "Longitude",
+        }
+        ds.lat.attrs = {
+            "long_name": "latitude of grid cell center",
+            "units": "degree_N",
+            "standard_name": "Latitude",
+        }
+
+        for attr in [ds.z, ds.lat, ds.lon]:
+            attr.encoding["_FillValue"] = None
+
+        from ..io.grid_io import _make_coordinatesystem_dict
+
+        ds.coords["ProjectionCoordinateSystem"] = xarray.DataArray(
+            data=np.array(1, dtype="int32"),
+            attrs=_make_coordinatesystem_dict(self),
+        )
+
+        projection = self.projection.copy()
+        if "_include_lon_0_lat_0" in projection:
+            projection["_include_lon_0_lat_0"] = str(
+                projection["_include_lon_0_lat_0"]
+            ).lower()
+        ds.coords["projection"] = xarray.DataArray(
+            data=np.array(1, dtype="int32"),
+            attrs=projection,
+        )
+
+        # Handle origin and radar attributes with appropriate dimensions
+        for attr_name in [
+            "origin_latitude",
+            "origin_longitude",
+            "origin_altitude",
+        ]:
+            if hasattr(self, attr_name):
+                attr_data = getattr(self, attr_name)
+                if attr_data is not None:
+                    attr_value = np.ma.expand_dims(attr_data["data"][0], 0)
+                    ds.coords[attr_name] = xarray.DataArray(
+                        attr_value, dims=("time",), attrs=get_metadata(attr_name)
+                    )
+
+        # Radar-specific attributes that should have the nradar dimension
+        for attr_name in [
+            "radar_altitude",
+            "radar_latitude",
+            "radar_longitude",
+            "radar_time",
+        ]:
+            if hasattr(self, attr_name):
+                attr_data = getattr(self, attr_name)
+                if attr_data is not None:
+                    ds.coords[attr_name] = xarray.DataArray(
+                        attr_data["data"],
+                        dims=("nradar",),
+                        attrs=get_metadata(attr_name),
+                    )
+
+        if "radar_time" in ds.variables:
+            ds.radar_time.attrs.pop("calendar")
+
+        # Handle radar_name and ensure it has the correct dimension
+        if self.radar_name is not None:
+            radar_name = _process_radar_name(self.radar_name["data"])
+            ds.coords["radar_name"] = xarray.DataArray(
+                radar_name, dims=("nradar",), attrs=get_metadata("radar_name")
+            )
+        else:
+            radar_name = np.array(["ExampleRadar"], dtype="U")
+            ds.coords["radar_name"] = xarray.DataArray(
+                radar_name, dims=("nradar",), attrs=get_metadata("radar_name")
+            )
+
+        # Add radar_name to attributes
+        ds.attrs["radar_name"] = (
+            radar_name.tolist() if radar_name.size > 1 else radar_name.item()
+        )
+        ds.attrs["nradar"] = radar_name.size
+        ds.attrs.update(self.metadata)
+        for key in ds.attrs:
+            try:
+                ds.attrs[key] = ds.attrs[key].decode("utf-8")
+            except AttributeError:
+                pass
+
+        ds.close()
         return ds
 
     def add_field(self, field_name, field_dict, replace_existing=False):

pyart/graph/__init__.py

@@ -62,5 +62,6 @@
 from .radardisplay_airborne import AirborneRadarDisplay  # noqa
 from .radarmapdisplay import RadarMapDisplay  # noqa
 from .radarmapdisplay_basemap import RadarMapDisplayBasemap  # noqa
+from .max_cappi import plot_maxcappi  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]