Max cappi

examples/plotting/plot_max_cappi.py

@@ -0,0 +1,112 @@
+"""
+=============
+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

@@ -333,55 +333,129 @@ def to_xarray(self):
         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()):
+
+            for meta, value in field_info.items():
                 if meta != "data":
-                    data.attrs.update({meta: self.fields[field][meta]})
+                    data.attrs.update({meta: value})
 
             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"),
+        ]
+
+        ds.z.attrs = get_metadata("z")
+        ds.y.attrs = get_metadata("y")
+        ds.x.attrs = get_metadata("x")
+
+        for attr in [ds.z, ds.lat, ds.lon]:
+            attr.encoding["_FillValue"] = None
+
+        # Delayed import
+        from ..io.grid_io import _make_coordinatesystem_dict
+
+        ds.coords["ProjectionCoordinateSystem"] = xarray.DataArray(
+            data=np.array(1, dtype="int32"),
+            attrs=_make_coordinatesystem_dict(self),
+        )
+
+        # write the projection dictionary as a scalar
+        projection = self.projection.copy()
+        # NetCDF does not support boolean attribute, covert to string
+        if "_include_lon_0_lat_0" in projection:
+            include = projection["_include_lon_0_lat_0"]
+            projection["_include_lon_0_lat_0"] = ["false", "true"][include]
+        ds.coords["projection"] = xarray.DataArray(
+            data=np.array(1, dtype="int32"),
+            dims=None,
+            attrs=projection,
+        )
+
+        for attr_name in [
+            "origin_latitude",
+            "origin_longitude",
+            "origin_altitude",
+            "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:
+                    if attr_name in [
+                        "origin_latitude",
+                        "origin_longitude",
+                        "origin_altitude",
+                    ]:
+                        # Adjusting the dims to 'time' for the origin attributes
+                        attr_value = np.ma.expand_dims(attr_data["data"][0], 0)
+                        dims = ("time",)
+                    else:
+                        if "radar_time" not in attr_name:
+                            attr_value = np.ma.expand_dims(attr_data["data"][0], 0)
+                        else:
+                            attr_value = [
+                                np.array(
+                                    num2date(
+                                        attr_data["data"][0],
+                                        units=attr_data["units"],
+                                    ),
+                                    dtype="datetime64[ns]",
+                                )
+                            ]
+                        dims = ("nradar",)
+
+                    ds.coords[attr_name] = xarray.DataArray(
+                        attr_value, dims=dims, attrs=get_metadata(attr_name)
+                    )
+
+        if "radar_time" in ds.variables:
+            ds.radar_time.attrs.pop("calendar")
+
+        if self.radar_name is not None:
+            radar_name = self.radar_name["data"]
+            ds["radar_name"] = xarray.DataArray(
+                np.array([b"".join(radar_name)]),
+                dims=("nradar"),
+                attrs=get_metadata("radar_name"),
+            )
+
+        ds.attrs = self.metadata
+        for key in ds.attrs:
+            try:
+                ds.attrs[key] = ds.attrs[key].decode("utf-8")
+            except AttributeError:
+                # If the attribute is not a byte string, just pass
+                pass
+
+        ds.close()
         return ds
 
     def add_field(self, field_name, field_dict, replace_existing=False):

pyart/graph/__init__.py

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

pyart/graph/gridmapdisplay.py

@@ -50,6 +50,8 @@
 except ImportError:
     _LAMBERT_GRIDLINES = False
 
+from . import max_cappi  # noqa
+
 
 class GridMapDisplay:
     """
@@ -120,7 +122,7 @@ def plot_grid(
         add_grid_lines=True,
         ticks=None,
         ticklabs=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Plot the grid using xarray and cartopy.
@@ -257,7 +259,7 @@ def plot_grid(
             vmin=vmin,
             vmax=vmax,
             add_colorbar=False,
-            **kwargs
+            **kwargs,
         )
 
         self.mappables.append(pm)
@@ -414,7 +416,7 @@ def plot_latitudinal_level(
         fig=None,
         ticks=None,
         ticklabs=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Plot a slice along a given latitude.
@@ -575,7 +577,7 @@ def plot_longitudinal_level(
         fig=None,
         ticks=None,
         ticklabs=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Plot a slice along a given longitude.
@@ -718,7 +720,7 @@ def plot_cross_section(
         fig=None,
         ticks=None,
         ticklabs=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Plot a cross section through a set of given points (latitude,
@@ -849,7 +851,7 @@ def plot_cross_section(
             add_colorbar=False,
             ax=ax,
             cmap=cmap,
-            **kwargs
+            **kwargs,
         )
 
         self.mappables.append(plot)
@@ -1110,6 +1112,46 @@ def cartopy_coastlines(self):
             category="physical", name="coastline", scale="10m", facecolor="none"
         )
 
+    def plot_maxcappi(
+        self,
+        field,
+        cmap=None,
+        vmin=None,
+        vmax=None,
+        title=None,
+        lat_lines=None,
+        lon_lines=None,
+        add_map=True,
+        projection=None,
+        colorbar=True,
+        range_rings=False,
+        dpi=100,
+        savedir=None,
+        show_figure=True,
+        add_slogan=False,
+        **kwargs,
+    ):
+        # Call the plot_maxcappi function from the max_cappi module or object
+        max_cappi.plot_maxcappi(
+            grid=self.grid,  # Assuming `self.grid` holds the Grid object in your class
+            field=field,
+            cmap=cmap,
+            vmin=vmin,
+            vmax=vmax,
+            title=title,
+            lat_lines=lat_lines,
+            lon_lines=lon_lines,
+            add_map=add_map,
+            projection=projection,
+            colorbar=colorbar,
+            range_rings=range_rings,
+            dpi=dpi,
+            savedir=savedir,
+            show_figure=show_figure,
+            add_slogan=add_slogan,
+            **kwargs,
+        )
+
 
 # These methods are a hack to allow gridlines when the projection is lambert
 # https://nbviewer.jupyter.org/gist/ajdawson/dd536f786741e987ae4e