Source code for radarx.vis.maxcappi

"""
Plot Max-CAPPI

This module provides a function to plot a Maximum Constant Altitude Plan Position Indicator (Max-CAPPI)
from radar data using an xarray dataset. The function includes options for adding map features, range rings,
color bars, and customized visual settings.

Author: Syed Hamid Ali (@syedhamidali)

.. autosummary::
   :nosignatures:
   :toctree: generated/

   {}
"""

__all__ = ["plot_maxcappi"]
__doc__ = __doc__.format("\n   ".join(__all__))

import os
import cmweather  # noqa
import cartopy.crs as ccrs
import cartopy.feature as feat
import matplotlib.pyplot as plt
import numpy as np
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from matplotlib.ticker import NullFormatter

# warnings.filterwarnings("ignore")


[docs] def plot_maxcappi( ds, data_var, 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, ): """ Plot a Maximum Constant Altitude Plan Position Indicator (Max-CAPPI) using an xarray Dataset. Parameters ---------- ds : xarray.Dataset Xarray Dataset containing the gridded radar data to be plotted. data_var : str The radar data_var to be plotted (e.g., "REF", "VEL", "WIDTH"). cmap : str or matplotlib colormap, optional Colormap to use for the plot. Default is "NWSRef". vmin : float, optional Minimum value for the color scaling. Default is set to the minimum value of the data if not provided. vmax : float, optional Maximum value for the color scaling. Default is set to the maximum value of the data if not provided. title : str, optional Title of the plot. If None, the title is set to "Max-{data_var}". lat_lines : array-like, optional Latitude lines to be included in the plot. Default is calculated based on dataset coordinates. lon_lines : array-like, optional Longitude lines to be included in the plot. Default is calculated based on dataset coordinates. add_map : bool, optional Whether to include a map background in the plot. Default is True. projection : cartopy.crs.Projection, optional The map projection for the plot. Default is automatically determined based on dataset coordinates. colorbar : bool, optional Whether to include a colorbar in the plot. Default is True. range_rings : bool, optional Whether to include range rings at 50 km intervals. Default is False. dpi : int, optional DPI (dots per inch) for the plot. Default is 100. savedir : str, optional Directory where the plot will be saved. If None, the plot is not saved. show_figure : bool, optional Whether to display the plot. Default is True. add_slogan : bool, optional Whether to add a slogan like "Powered by Radarx" to the plot. Default is False. **kwargs : dict, optional Additional keyword arguments to pass to matplotlib's `pcolormesh` function. Returns ------- None This function does not return any value. It generates and optionally displays or saves a plot. Notes ----- - The function extracts the maximum value across the altitude (z) dimension to create the Max-CAPPI. - It supports customizations such as map projections, color scales, and range rings. - If the radar_name attribute in the dataset is a byte string, it will be decoded and limited to 4 characters. - If add_map is True, map features and latitude/longitude lines are included. - The plot can be saved to a specified directory in PNG format. Author: Syed Hamid Ali (@syedhamidali) """ # Define default latitude and longitude lines if not provided if lon_lines is None: lon_lines = np.arange(int(ds.lon.min().values), int(ds.lon.max().values) + 1) if lat_lines is None: lat_lines = np.arange(int(ds.lat.min().values), int(ds.lat.max().values) + 1) # Configure plotting parameters plt.rcParams.update( { "font.weight": "bold", "axes.labelweight": "bold", "xtick.direction": "in", "ytick.direction": "in", "xtick.major.size": 10, "ytick.major.size": 10, "xtick.minor.size": 7, "ytick.minor.size": 7, "font.size": 14, "axes.linewidth": 2, "ytick.labelsize": 12, "xtick.labelsize": 12, } ) # Extract maximum values across dimensions max_c = ds[data_var].max(dim="z") max_x = ds[data_var].max(dim="y") max_y = ds[data_var].max(dim="x").T # Extract coordinate and dimension information trgx, trgy, trgz = ds["x"].values, ds["y"].values, ds["z"].values max_height = int(np.floor(trgz.max()) / 1e3) sideticks = np.arange(max_height / 4, max_height + 1, max_height / 4).astype(int) # Set colormap and value ranges cmap = cmap or "ChaseSpectral" vmin = vmin if vmin is not None else ds[data_var].min().values vmax = vmax if vmax is not None else ds[data_var].max().values title = title or f"Max-{data_var.upper()[:3]}" projection = _get_projection(ds) # FIG fig = plt.figure(figsize=[10.3, 10]) left, bottom, width, height = 0.1, 0.1, 0.6, 0.2 ax_xy = plt.axes((left, bottom, width, width), projection=projection) ax_x = plt.axes((left, bottom + width, width, height)) ax_y = plt.axes((left + width, bottom, height, width)) ax_cnr = plt.axes((left + width, bottom + width, left + left, height)) if colorbar: # pragma: no cover ax_cb = plt.axes((left - 0.015 + width + height + 0.02, bottom, 0.02, width)) # Set axis label formatters ax_x.xaxis.set_major_formatter(NullFormatter()) ax_y.yaxis.set_major_formatter(NullFormatter()) ax_cnr.yaxis.set_major_formatter(NullFormatter()) ax_cnr.xaxis.set_major_formatter(NullFormatter()) ax_x.set_ylabel("Height (km)", size=13) ax_y.set_xlabel("Height (km)", size=13) # Draw CAPPI plt.sca(ax_xy) xy = ax_xy.pcolormesh(trgx, trgy, max_c, cmap=cmap, vmin=vmin, vmax=vmax, **kwargs) # Add map features if add_map: # pragma: no cover _add_map_features(ax_xy, lat_lines, lon_lines) ax_xy.minorticks_on() if range_rings: # pragma: no cover _plot_range_rings(ax_xy, trgx.max()) ax_xy.set_xlim(trgx.min(), trgx.max()) ax_xy.set_ylim(trgx.min(), trgx.max()) # Draw colorbar if colorbar: # pragma: no cover cb = plt.colorbar(xy, cax=ax_cb) units = ds[data_var].attrs.get("units", "") # Provide a default cb.set_label(units, size=15) background_color = ax_xy.get_facecolor() color = "k" if sum(background_color[:3]) / 3 > 0.5 else "w" plt.sca(ax_x) plt.pcolormesh(trgx / 1e3, trgz / 1e3, max_x, cmap=cmap, vmin=vmin, vmax=vmax) plt.yticks(sideticks) ax_x.set_xlim(trgx.min() / 1e3, trgx.max() / 1e3) ax_x.grid(axis="y", lw=0.5, color=color, alpha=0.5, ls=":") ax_x.minorticks_on() plt.sca(ax_y) plt.pcolormesh(trgz / 1e3, trgy / 1e3, max_y, cmap=cmap, vmin=vmin, vmax=vmax) ax_y.set_xticks(sideticks) ax_y.set_ylim(trgx.min() / 1e3, trgx.max() / 1e3) ax_y.grid(axis="x", lw=0.5, color=color, alpha=0.5, ls=":") ax_y.minorticks_on() plt.sca(ax_cnr) plt.tick_params( axis="both", # changes apply to both axes which="both", # both major and minor ticks are affected bottom=False, # ticks along the bottom edge are off top=False, # ticks along the top edge are off left=False, right=False, labelbottom=False, ) # Title and metadata _add_plot_metadata(ax_cnr, ds, title, trgx, trgz) # pragma: no cover ax_xy.set_aspect("auto") if add_slogan: # pragma: no cover _add_slogan(fig) _save_and_display_plot( fig, ds, title, savedir, dpi, show_figure ) # pragma: no cover # plt.rcParams.update(original_rc_params) plt.rcdefaults() return fig
def _add_slogan(fig): """Add Slogan""" text = "Plot by Radarx" + " | " + "Powered by Xradar" fig.text( 0.5, 0.05, text, fontsize=10, fontname="Courier New", ha="center", # bbox=dict(facecolor='none', boxstyle='round,pad=0.5') ) def _save_and_display_plot(fig, ds, title, savedir, dpi, show_figure): """ Save the plot to file or display it based on user preference. """ if savedir: # pragma: no cover time_str = ds["time"].dt.strftime("%Y%m%d%H%M%S").values.item() radar_name = ds.attrs.get("instrument_name", "Radar") filepath = f"{savedir}{os.sep}{title}_{radar_name}_{time_str}.png" plt.savefig(filepath, dpi=dpi, bbox_inches="tight") print(f"Figure saved as {filepath}") plt.show() if show_figure else plt.close() # pragma: no cover def _add_plot_metadata(ax_cnr, ds, title, trgx, trgz): # pragma: no cover """ Add metadata such as radar name, date, time, and title to the plot. """ radar_name = ds.attrs.get("radar_name", "").split(",")[0][:4] date_str = ds["time"].dt.strftime("%d %b, %Y UTC").values.item() time_str = ds["time"].dt.strftime("%H:%M:%S Z").values.item() max_range, max_height = int(trgx.max() / 1e3), int(trgz.max() / 1e3) ax_cnr.text(0.5, 0.9, radar_name, size=13, weight="bold", ha="center", va="center") ax_cnr.text(0.5, 0.76, title, size=13, weight="bold", ha="center", va="center") ax_cnr.text( 0.5, 0.63, f"Max Range: {max_range} km", size=11, ha="center", va="center" ) ax_cnr.text( 0.5, 0.47, f"Max Height: {max_height} km", size=11, ha="center", va="center" ) ax_cnr.text(0.5, 0.28, time_str, weight="bold", size=16, ha="center", va="center") ax_cnr.text(0.5, 0.13, date_str, size=13.5, ha="center", va="center") def _plot_range_rings(ax_xy, max_range): # pragma: no cover """ Plots range rings at 50 km intervals. Parameters ---------- ax_xy : matplotlib.axes.Axes The axis on which to plot the range rings. max_range : float The maximum range for the range rings. Returns ------- None """ background_color = ax_xy.get_facecolor() # pragma: no cover color = "k" if sum(background_color[:3]) / 3 > 0.5 else "w" # pragma: no cover for i, r in enumerate( np.arange(5e4, np.floor(max_range) + 1, 5e4) ): # pragma: no cover label = f"Ring Dist. {int(r/1e3)} km" if i == 0 else None # pragma: no cover ax_xy.plot( r * np.cos(np.arange(0, 360) * np.pi / 180), r * np.sin(np.arange(0, 360) * np.pi / 180), color=color, ls="--", linewidth=0.4, alpha=0.3, label=label, ) # pragma: no cover ax_xy.legend(loc="upper right", prop={"weight": "normal", "size": 8}) def _get_projection(ds): # pragma: no cover """ Determine the central latitude and longitude from a dataset and return the corresponding projection. Parameters ---------- ds : xarray.Dataset The dataset from which to extract latitude and longitude information. Returns ------- projection : cartopy.crs.Projection A Cartopy projection object centered on the extracted or calculated latitude and longitude. """ def _get_coord_or_attr(ds, coord_name, attr_name): # pragma: no cover """Helper function to get a coordinate or attribute, or calculate median if available. """ if coord_name in ds: # pragma: no cover return ( ds[coord_name].values.item() if ds[coord_name].values.ndim == 0 # pragma: no cover else ds[coord_name].values[0] # pragma: no cover ) if f"origin_{coord_name}" in ds.coords: # pragma: no cover return ds.coords[f"origin_{coord_name}"].median().item() if f"radar_{coord_name}" in ds.coords: # pragma: no cover return ds.coords[f"radar_{coord_name}"].median().item() return ds.attrs.get(attr_name, None) lat_0 = _get_coord_or_attr(ds, "latitude", "origin_latitude") or _get_coord_or_attr( ds, "radar_latitude", "origin_latitude" ) lon_0 = _get_coord_or_attr( ds, "longitude", "origin_longitude" ) or _get_coord_or_attr(ds, "radar_longitude", "origin_longitude") if lat_0 is None or lon_0 is None: # pragma: no cover lat_0 = ds.lat.mean().item() lon_0 = ds.lon.mean().item() projection = ccrs.LambertAzimuthalEqualArea(lon_0, lat_0) return projection def _add_map_features(ax, lat_lines, lon_lines): # pragma: no cover """ Adds map features and gridlines to the plot. Parameters ---------- ax : matplotlib.axes.Axes The axis on which to add map features and gridlines. lat_lines : array-like Latitude lines for the gridlines. lon_lines : array-like Longitude lines for the gridlines. Returns ------- None """ background_color = ax.get_facecolor() color = "k" if sum(background_color[:3]) / 3 > 0.5 else "w" # Labeling gridlines depending on the projection if isinstance(ax.projection, (ccrs.PlateCarree, ccrs.Mercator)): # pragma: no cover gl = ax.gridlines( xlocs=lon_lines, ylocs=lat_lines, linewidth=1, alpha=0.5, linestyle="--", draw_labels=True, ) gl.top_labels = False gl.right_labels = False gl.xlines = False gl.ylines = False gl.xlabel_style = {"color": color} gl.ylabel_style = {"color": color} ax.add_feature(feat.COASTLINE, alpha=0.8, lw=1, ec=color) ax.add_feature(feat.BORDERS, alpha=0.7, lw=0.7, ls="--", ec=color) ax.add_feature( feat.STATES.with_scale("10m"), alpha=0.6, lw=0.5, ls=":", ec=color ) elif isinstance( ax.projection, (ccrs.LambertConformal, ccrs.LambertAzimuthalEqualArea), # pragma: no cover ): ax.figure.canvas.draw() gl = ax.gridlines( crs=ccrs.PlateCarree(), xlocs=lon_lines, ylocs=lat_lines, linewidth=1, alpha=0.5, linestyle="--", draw_labels=False, ) gl.xlines = False gl.ylines = False gl.xlabel_style = {"color": color} gl.ylabel_style = {"color": color} ax.add_feature(feat.COASTLINE, alpha=0.8, lw=1, ec=color) ax.add_feature(feat.BORDERS, alpha=0.7, lw=0.7, ls="--", ec=color) ax.add_feature( feat.STATES.with_scale("10m"), alpha=0.6, lw=0.5, ls=":", ec=color ) # Label the end-points of the gridlines using custom tick makers ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) ax.yaxis.set_major_formatter(LATITUDE_FORMATTER) import os from contextlib import redirect_stdout # Suppress output during import with open(os.devnull, "w") as fnull, redirect_stdout(fnull): from pyart.graph.gridmapdisplay import lambert_xticks, lambert_yticks lambert_xticks(ax, lon_lines) # pragma: no cover lambert_yticks(ax, lat_lines) # pragma: no cover else: # pragma: no cover ax.gridlines(xlocs=lon_lines, ylocs=lat_lines) # pragma: no cover