import math
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import warnings
from .. import retrieval
from matplotlib.axes import Axes
from datatree import DataTree
try:
from cartopy.mpl.geoaxes import GeoAxes
CARTOPY_AVAILABLE = True
except ImportError:
CARTOPY_AVAILABLE = False
GeoAxes._pcolormesh_patched = Axes.pcolormesh
[docs]
def plot_horiz_xsection_quiver(
Grids,
ax=None,
background_field="reflectivity",
level=1,
cmap="ChaseSpectral",
vmin=0,
vmax=70,
u_vel_contours=None,
v_vel_contours=None,
w_vel_contours=None,
wind_vel_contours=None,
u_field="u",
v_field="v",
w_field="w",
show_lobes=True,
title_flag=True,
axes_labels_flag=True,
colorbar_flag=True,
colorbar_contour_flag=False,
bg_grid_no=0,
scale=3,
quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0,
contour_alpha=0.7,
quiverkey_len=5.0,
quiverkey_loc="best",
quiver_width=0.01,
):
"""
This procedure plots a horizontal cross section of winds from wind fields
generated by PyDDA using quivers. The length of the quivers varies
with horizontal wind speed.
Parameters
----------
Grids: list or DataTree
List of Py-DDA Grids to visualize
ax: matplotlib axis handle
The axis handle to place the plot on. Set to None to plot on the
current axis.
background_field: str
The name of the background field to plot the quivers on.
level: int
The number of the vertical level to plot the cross section through.
cmap: str or matplotlib colormap
The name of the matplotlib colormap to use for the background field.
vmin: float
The minimum bound to use for plotting the background field. None will
automatically detect the background field minimum.
vmax: float
The maximum bound to use for plotting the background field. None will
automatically detect the background field maximum.
u_vel_contours: 1-D array
The contours to use for plotting contours of u. Set to None to not
display such contours.
v_vel_contours: 1-D array
The contours to use for plotting contours of v. Set to None to not
display such contours.
w_vel_contours: 1-D array
The contours to use for plotting contours of w. Set to None to not
display such contours.
wind_vel_contours: 1-D array
The contours to use for plotting contours of horizontal wind speed.
Set to None to not display such contours
u_field: str
Name of zonal wind (u) field in Grids.
v_field: str
Name of meridional wind (v) field in Grids.
w_field: str
Name of vertical wind (w) field in Grids.
show_lobes: bool
If True, the dual doppler lobes from each pair of radars will be shown.
title_flag: bool
If True, PyDDA will generate a title for the plot.
axes_labels_flag: bool
If True, PyDDA will generate axes labels for the plot
colorbar_flag: bool
If True, PyDDA will generate a colorbar for the plot background field.
colorbar_contour_flag: bool
If True, PyDDA will generate a colorbar for the contours.
bg_grid_no: int
Number of grid in Grids to take background field from.
Set to -1 to use maximum value from all grids.
quiver_spacing_x_km: float
Spacing in km between quivers in x axis.
quiver_spacing_y_km: float
Spacing in km between quivers in y axis.
contour_alpha: float
Alpha (transparency) of velocity contours. 0 = transparent, 1 = opaque
quiverkey_len: float
Length to use for the quiver key in m/s.
quiverkey_loc: str
Location of quiverkey. One of:
'best'
'top_left'
'top'
'top_right'
'bottom_left'
'bottom'
'bottom_right'
'left'
'right'
'top_left_outside'
'top_right_outside'
'bottom_left_outside'
'bottom_right_outside'
'best' will put the quiver key in the corner with the fewest amount of
valid data points while keeping the quiver key inside the plot.
The rest of the options will put the quiver key in that
particular part of the plot.
quiver_width: float
The width of the lines for the quiver. Use this to specify
the thickness of the quiver lines. Units are in fraction of plot
width.
Returns
-------
ax: Matplotlib axis handle
The matplotlib axis handle associated with the plot.
"""
if isinstance(Grids, DataTree):
child_list = list(Grids.children.keys())
grid_list = []
rad_names = []
for child in child_list:
if "radar" in child:
grid_list.append(Grids[child].to_dataset())
rad_names.append(child)
bca_min = math.radians(Grids[u_field].attrs["min_bca"])
bca_max = math.radians(Grids[u_field].attrs["max_bca"])
else:
grid_list = Grids
bca_min = math.radians(grid_list[0][u_field].attrs["min_bca"])
bca_max = math.radians(grid_list[0][u_field].attrs["max_bca"])
grid_bg = grid_list[bg_grid_no][background_field].values.squeeze()
grid_bg = np.ma.masked_invalid(grid_bg)
if not CARTOPY_AVAILABLE:
raise ModuleNotFoundError(
"Cartopy needs to be installed in order to use plotting module!"
)
if vmin is None:
vmin = grid_bg.min()
if vmax is None:
vmax = grid_bg.max()
grid_h = grid_list[0]["point_altitude"].values / 1e3
grid_x = grid_list[0]["point_x"].values / 1e3
grid_y = grid_list[0]["point_y"].values / 1e3
dx = np.diff(grid_x, axis=2)[0, 0, 0]
dy = np.diff(grid_y, axis=1)[0, 0, 0]
if isinstance(Grids, DataTree):
u = Grids[u_field].values.squeeze()
v = Grids[v_field].values.squeeze()
w = Grids[w_field].values.squeeze()
else:
u = grid_list[0][u_field].values.squeeze()
v = grid_list[0][v_field].values.squeeze()
w = grid_list[0][w_field].values.squeeze()
qloc_x, qloc_y = _parse_quiverkey_string(
quiverkey_loc, grid_h[level], grid_x[level], grid_y[level], grid_bg[level]
)
if ax is None:
ax = plt.gca()
the_mesh = ax.pcolormesh(
grid_x[level, :, :],
grid_y[level, :, :],
grid_bg[level, :, :],
cmap=cmap,
vmin=vmin,
vmax=vmax,
)
np.ma.sqrt(u**2 + v**2)
quiver_density_x = int((1 / dx) * quiver_spacing_x_km)
quiver_density_y = int((1 / dy) * quiver_spacing_y_km)
q = ax.quiver(
grid_x[level, ::quiver_density_y, ::quiver_density_x],
grid_y[level, ::quiver_density_y, ::quiver_density_x],
u[level, ::quiver_density_y, ::quiver_density_x],
v[level, ::quiver_density_y, ::quiver_density_x],
color="k",
scale=25.0 * quiverkey_len,
scale_units="width",
width=quiver_width,
)
quiver_font = {
"family": "sans-serif",
"style": "normal",
"variant": "normal",
"weight": "bold",
"size": "medium",
}
ax.quiverkey(
q,
qloc_x,
qloc_y,
quiverkey_len,
label=(str(quiverkey_len) + " m/s"),
fontproperties=quiver_font,
)
if colorbar_flag is True:
cp = grid_list[bg_grid_no][background_field].attrs["long_name"]
cp.replace(" ", "_")
cp = cp + " [" + grid_list[bg_grid_no][background_field].attrs["units"]
cp = cp + "]"
plt.colorbar(the_mesh, ax=ax, label=(cp))
if u_vel_contours is not None:
u_filled = np.ma.filled(u[level, :, :], fill_value=np.nan)
cs = ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
u_filled,
levels=u_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="U [m/s]")
if v_vel_contours is not None:
v_filled = np.ma.filled(v[level, :, :], fill_value=np.nan)
cs = ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
v_filled,
levels=u_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="V [m/s]")
if w_vel_contours is not None:
w_filled = np.ma.filled(w[level, :, :], fill_value=np.nan)
cs = ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
w_filled,
levels=w_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="W [m/s]")
if wind_vel_contours is not None:
vel = np.ma.sqrt(u[level, :, :] ** 2 + v[level, :, :] ** 2)
# vel = vel.filled(fill_value=np.nan)
cs = ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
vel,
levels=wind_vel_contours,
linewidths=2,
)
cs.set_clim([np.min(wind_vel_contours), np.max(wind_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="|V| [m/s]")
if show_lobes is True:
for i in range(len(grid_list)):
for j in range(len(grid_list)):
if i != j:
bca = retrieval.get_bca(
grid_list[i],
grid_list[j],
)
ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
bca,
levels=[bca_min, bca_max],
color="k",
)
if axes_labels_flag is True:
ax.set_xlabel(("X [km]"))
ax.set_ylabel(("Y [km]"))
if title_flag is True:
ax.set_title(("PyDDA retreived winds @" + str(grid_h[level, 0, 0]) + " km"))
ax.set_xlim([grid_x.min(), grid_x.max()])
ax.set_ylim([grid_y.min(), grid_y.max()])
return ax
[docs]
def plot_horiz_xsection_quiver_map(
Grids,
ax=None,
background_field="reflectivity",
level=1,
cmap="ChaseSpectral",
vmin=0,
vmax=70,
u_vel_contours=None,
v_vel_contours=None,
w_vel_contours=None,
wind_vel_contours=None,
u_field="u",
v_field="v",
w_field="w",
show_lobes=True,
title_flag=True,
axes_labels_flag=True,
colorbar_flag=True,
colorbar_contour_flag=False,
bg_grid_no=0,
contour_alpha=0.7,
coastlines=True,
quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0,
gridlines=True,
quiverkey_len=5.0,
quiverkey_loc="best",
quiver_width=0.01,
):
"""
This procedure plots a horizontal cross section of winds from wind fields
generated by PyDDA using quivers onto a geographical map. The length of
the quivers varies with wind speed.
Parameters
----------
Grids: list or DataTree
List of Py-DDA Grids to visualize
ax: matplotlib axis handle (with cartopy ccrs)
The axis handle to place the plot on. Set to None to create a new map.
Note: the axis needs to be in a PlateCarree() projection. Support for
other projections is planned in the future.
background_field: str
The name of the background field to plot the quivers on.
level: int
The number of the vertical level to plot the cross section through.
cmap: str or matplotlib colormap
The name of the matplotlib colormap to use for the background field.
vmin: float
The minimum bound to use for plotting the background field. None will
automatically detect the background field minimum.
vmax: float
The maximum bound to use for plotting the background field. None will
automatically detect the background field maximum.
u_vel_contours: 1-D array
The contours to use for plotting contours of u. Set to None to not
display such contours.
v_vel_contours: 1-D array
The contours to use for plotting contours of v. Set to None to not
display such contours.
w_vel_contours: 1-D array
The contours to use for plotting contours of w. Set to None to not
display such contours.
u_field: str
Name of zonal wind (u) field in Grids.
v_field: str
Name of meridional wind (v) field in Grids.
w_field: str
Name of vertical wind (w) field in Grids.
show_lobes: bool
If True, the dual doppler lobes from each pair of radars will be shown.
title_flag: bool
If True, PyDDA will generate a title for the plot.
axes_labels_flag: bool
If True, PyDDA will generate axes labels for the plot.
colorbar_flag: bool
If True, PyDDA will generate a colorbar for the plot background field.
colorbar_contour_flag: bool
If True, PyDDA will generate a colorbar for the contours.
bg_grid_no: int
Number of grid in Grids to take background field from.
Set to -1 to use maximum value from all grids.
contour_alpha: float
Alpha (transparency) of velocity contours. 0 = transparent, 1 = opaque
coastlines: bool
Set to true to display coastlines.
quiver_spacing_x_km: float
Spacing in km between quivers in x axis.
quiver_spacing_y_km: float
Spacing in km between quivers in y axis.
gridlines: bool
Set to true to show grid lines.
quiverkey_len: float
Length to use for the quiver key in m/s.
quiverkey_loc: str
Location of quiverkey. One of:
'best'
'top_left'
'top'
'top_right'
'bottom_left'
'bottom'
'bottom_right'
'left'
'right'
'top_left_outside'
'top_right_outside'
'bottom_left_outside'
'bottom_right_outside'
'best' will put the quiver key in the corner with the fewest amount of
valid data points while keeping the quiver key inside the plot.
The rest of the options will put the quiver key in that
particular part of the plot.
quiver_width: float
The width of the lines for the quiver given as a fraction
relative to the plot width. Use this to specify the thickness
of the quiver lines.
Returns
-------
ax: matplotlib axis
Axis handle to output axis
"""
if isinstance(Grids, DataTree):
child_list = list(Grids.children.keys())
grid_list = []
rad_names = []
for child in child_list:
if "radar" in child:
grid_list.append(Grids[child].to_dataset())
rad_names.append(child)
bca_min = math.radians(Grids[u_field].attrs["min_bca"])
bca_max = math.radians(Grids[u_field].attrs["max_bca"])
else:
grid_list = Grids
bca_min = math.radians(grid_list[0][u_field].attrs["min_bca"])
bca_max = math.radians(grid_list[0][u_field].attrs["max_bca"])
if not CARTOPY_AVAILABLE:
raise ModuleNotFoundError(
"Cartopy needs to be installed in order to use plotting module!"
)
if bg_grid_no > -1:
grid_bg = grid_list[bg_grid_no][background_field].values.squeeze()
else:
grid_array = np.ma.stack(
[x[background_field].values.squeeze() for x in grid_list]
)
grid_bg = grid_array.max(axis=0)
grid_bg = np.ma.masked_invalid(grid_bg)
if vmin is None:
vmin = grid_bg.min()
if vmax is None:
vmax = grid_bg.max()
grid_h = grid_list[0]["point_altitude"].values / 1e3
grid_x = grid_list[0]["point_x"].values / 1e3
grid_y = grid_list[0]["point_y"].values / 1e3
grid_lat = grid_list[0].point_latitude.values[level]
grid_lon = grid_list[0].point_longitude.values[level]
qloc_x, qloc_y = _parse_quiverkey_string(
quiverkey_loc, grid_h[level], grid_x[level], grid_y[level], grid_bg[level]
)
dx = np.diff(grid_x, axis=2)[0, 0, 0]
dy = np.diff(grid_y, axis=1)[0, 0, 0]
if isinstance(Grids, DataTree):
u = Grids[u_field].values.squeeze()
v = Grids[v_field].values.squeeze()
w = Grids[w_field].values.squeeze()
else:
if np.ma.isMaskedArray(grid_list[0][u_field].values.squeeze()):
u = grid_list[0][u_field].values.squeeze().filled(fill_value=np.nan)
else:
u = grid_list[0][u_field].values.squeeze()
if np.ma.isMaskedArray(grid_list[0][v_field].values):
v = grid_list[0][v_field].values.squeeze().filled(fill_value=np.nan)
else:
v = grid_list[0][v_field].values.squeeze()
if np.ma.isMaskedArray(grid_list[0][u_field].values):
w = grid_list[0][w_field].values.squeeze().filled(fill_value=np.nan)
else:
w = grid_list[0][w_field].values.squeeze()
transform = ccrs.PlateCarree()
if ax is None:
ax = plt.axes(projection=transform)
the_mesh = ax.pcolormesh(
grid_lon[:, :],
grid_lat[:, :],
grid_bg[level, :, :],
cmap=cmap,
transform=transform,
zorder=0,
vmin=vmin,
vmax=vmax,
)
np.ma.sqrt(u**2 + v**2)
quiver_density_x = int((1 / dx) * quiver_spacing_x_km)
quiver_density_y = int((1 / dy) * quiver_spacing_y_km)
q = ax.quiver(
grid_lon[::quiver_density_y, ::quiver_density_x],
grid_lat[::quiver_density_y, ::quiver_density_x],
u[level, ::quiver_density_y, ::quiver_density_x],
v[level, ::quiver_density_y, ::quiver_density_x],
transform=transform,
width=quiver_width,
scale=25.0 * quiverkey_len,
)
quiver_font = {
"family": "sans-serif",
"style": "normal",
"variant": "normal",
"weight": "bold",
"size": "medium",
}
ax.quiverkey(
q,
qloc_x,
qloc_y,
quiverkey_len,
label=(str(quiverkey_len) + " m/s"),
fontproperties=quiver_font,
)
if colorbar_flag is True:
cp = grid_list[bg_grid_no][background_field].attrs["long_name"]
cp.replace(" ", "_")
cp = cp + " [" + grid_list[bg_grid_no][background_field].attrs["units"]
cp = cp + "]"
plt.colorbar(the_mesh, ax=ax, label=(cp))
if u_vel_contours is not None:
u_filled = np.ma.masked_where(
u[level, :, :] < np.min(u_vel_contours), u[level, :, :]
)
try:
cs = ax.contour(
grid_lon[:, :],
grid_lat[:, :],
u_filled,
levels=u_vel_contours,
linewidths=2,
zorder=2,
extend="both",
)
cs.set_clim([np.min(u_vel_contours), np.max(u_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_over(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(
cs, ax=ax, label="U [m/s]", extend="both", spacing="proportional"
)
except ValueError:
warnings.warn(
(
"Cartopy does not support blank contour plots, "
+ "contour color map not drawn!"
),
RuntimeWarning,
)
if v_vel_contours is not None:
v_filled = np.ma.masked_where(
v[level, :, :] < np.min(v_vel_contours), v[level, :, :]
)
try:
cs = ax.contour(
grid_lon[:, :],
grid_lat[:, :],
v_filled,
levels=u_vel_contours,
linewidths=2,
zorder=2,
extend="both",
)
cs.set_clim([np.min(v_vel_contours), np.max(v_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_over(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(
cs, ax=ax, label="V [m/s]", extend="both", spacing="proportional"
)
except ValueError:
warnings.warn(
(
"Cartopy does not support blank contour plots, "
+ "contour color map not drawn!"
),
RuntimeWarning,
)
if w_vel_contours is not None:
w_filled = np.ma.masked_where(
w[level, :, :] < np.min(w_vel_contours), w[level, :, :]
)
try:
cs = ax.contour(
grid_lon[::, ::],
grid_lat[::, ::],
w_filled,
levels=w_vel_contours,
linewidths=2,
zorder=2,
extend="both",
)
cs.set_clim([np.min(w_vel_contours), np.max(w_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_over(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(
cs,
ax=ax,
label="W [m/s]",
extend="both",
spacing="proportional",
ticks=w_vel_contours,
)
except ValueError:
warnings.warn(
(
"Cartopy does not support color maps on blank "
+ "contour plots, contour color map not drawn!"
),
RuntimeWarning,
)
if wind_vel_contours is not None:
vel = np.ma.sqrt(u[level, :, :] ** 2 + v[level, :, :] ** 2)
vel = vel.filled(fill_value=np.nan)
try:
cs = ax.contour(
grid_x[level, :, :],
grid_y[level, :, :],
vel,
levels=wind_vel_contours,
linewidths=2,
)
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(
cs,
ax=ax,
label="|V\ [m/s]",
extend="both",
spacing="proportional",
ticks=w_vel_contours,
)
except ValueError:
warnings.warn(
(
"Cartopy does not support color maps on blank "
+ "contour plots, contour color map not drawn!"
),
RuntimeWarning,
)
if show_lobes is True:
for i in range(len(grid_list)):
for j in range(len(grid_list)):
if i != j:
bca = retrieval.get_bca(grid_list[i], grid_list[j])
ax.contour(
grid_lon[:, :],
grid_lat[:, :],
bca,
levels=[bca_min, bca_max],
color="k",
zorder=1,
)
if axes_labels_flag is True:
ax.set_xlabel(("Latitude [$^{\circ}$]"))
ax.set_ylabel(("Longitude [$^{\circ}$]"))
if title_flag is True:
ax.set_title(("PyDDA retreived winds @" + str(grid_h[level, 0, 0]) + " km"))
if coastlines is True:
ax.coastlines(resolution="10m")
if gridlines is True:
ax.gridlines()
ax.set_extent([grid_lon.min(), grid_lon.max(), grid_lat.min(), grid_lat.max()])
num_tenths = int(round((grid_lon.max() - grid_lon.min()) * 10) + 1)
the_ticks_x = np.round(np.linspace(grid_lon.min(), grid_lon.max(), num_tenths), 1)
num_tenths = int(round((grid_lat.max() - grid_lat.min()) * 10) + 1)
the_ticks_y = np.round(np.linspace(grid_lat.min(), grid_lat.max(), num_tenths), 1)
ax.set_xticks(the_ticks_x)
ax.set_yticks(the_ticks_y)
return ax
[docs]
def plot_xz_xsection_quiver(
Grids,
ax=None,
background_field="reflectivity",
level=1,
cmap="ChaseSpectral",
vmin=0,
vmax=70,
u_vel_contours=None,
v_vel_contours=None,
w_vel_contours=None,
wind_vel_contours=None,
u_field="u",
v_field="v",
w_field="w",
title_flag=True,
axes_labels_flag=True,
colorbar_flag=True,
colorbar_contour_flag=False,
bg_grid_no=0,
quiver_spacing_x_km=10.0,
quiver_spacing_z_km=1.0,
contour_alpha=0.7,
quiverkey_len=5.0,
quiverkey_loc="best",
quiver_width=0.01,
):
"""
This procedure plots a cross section of winds from wind fields
generated by PyDDA in the X-Z plane using quivers.
The length of the quivers varies with wind speed.
Parameters
----------
Grids: list or DataTree
List of Py-DDA Grids to visualize
ax: matplotlib axis handle
The axis handle to place the plot on. Set to None to plot on the
current axis.
background_field: str
The name of the background field to plot the quivers on.
level: int
The number of the Y level to plot the cross section through.
cmap: str or matplotlib colormap
The name of the matplotlib colormap to use for the background field.
vmin: float
The minimum bound to use for plotting the background field. None will
automatically detect the background field minimum.
vmax: float
The maximum bound to use for plotting the background field. None will
automatically detect the background field maximum.
u_vel_contours: 1-D array
The contours to use for plotting contours of u. Set to None to not
display such contours.
v_vel_contours: 1-D array
The contours to use for plotting contours of v. Set to None to not
display such contours.
w_vel_contours: 1-D array
The contours to use for plotting contours of w. Set to None to not
display such contours.
u_field: str
Name of zonal wind (u) field in Grids.
v_field: str
Name of meridional wind (v) field in Grids.
w_field: str
Name of vertical wind (w) field in Grids.
show_lobes: bool
If True, the dual doppler lobes from each pair of radars will be shown.
title_flag: bool
If True, PyDDA will generate a title for the plot.
axes_labels_flag: bool
If True, PyDDA will generate axes labels for the plot
colorbar_flag: bool
If True, PyDDA will generate a colorbar for the plot background field.
colorbar_contour_flag: bool
If True, PyDDA will generate a colorbar for the contours.
bg_grid_no: int
Number of grid in Grids to take background field from.
Set to -1 to use maximum value from all grids.
quiver_spacing_x_km: float
Spacing in km between quivers in x axis.
quiver_spacing_z_km: float
Spacing in km between quivers in z axis.
contour_alpha: float
Alpha (transparency) of velocity contours. 0 = transparent, 1 = opaque
quiverkey_len: float
Length to use for the quiver key in m/s.
quiverkey_loc: str
Location of quiverkey. One of:
'best'
'top_left'
'top'
'top_right'
'bottom_left'
'bottom'
'bottom_right'
'left'
'right'
'top_left_outside'
'top_right_outside'
'bottom_left_outside'
'bottom_right_outside'
'best' will put the quiver key in the corner with the fewest amount of
valid data points while keeping the quiver key inside the plot.
The rest of the options will put the quiver key in that
particular part of the plot.
quiver_width: float
The width of the lines for the quiver. Use this to specify
the thickness of the quiver lines. Units are fraction of the plot
width.
Returns
-------
ax: matplotlib axis
Axis handle to output axis
"""
if isinstance(Grids, DataTree):
child_list = list(Grids.children.keys())
grid_list = []
rad_names = []
for child in child_list:
if "radar" in child:
grid_list.append(Grids[child].to_dataset())
rad_names.append(child)
else:
grid_list = Grids
grid_bg = grid_list[bg_grid_no][background_field].values.squeeze()
grid_bg = np.ma.masked_invalid(grid_bg)
if vmin is None:
vmin = grid_bg.min()
if vmax is None:
vmax = grid_bg.max()
grid_h = grid_list[0]["point_altitude"].values / 1e3
grid_x = grid_list[0]["point_x"].values / 1e3
grid_y = grid_list[0]["point_y"].values / 1e3
dx = np.diff(grid_x, axis=2)[0, 0, 0]
dz = np.diff(grid_y, axis=1)[0, 0, 0]
if isinstance(Grids, DataTree):
u = Grids[u_field].values.squeeze()
v = Grids[v_field].values.squeeze()
w = Grids[w_field].values.squeeze()
else:
u = grid_list[0][u_field].values.squeeze()
v = grid_list[0][v_field].values.squeeze()
w = grid_list[0][w_field].values.squeeze()
qloc_x, qloc_y = _parse_quiverkey_string(
quiverkey_loc,
grid_h[:, level, :],
grid_x[:, level, :],
grid_y[:, level, :],
grid_bg[:, level, :],
xsection="xz",
)
if ax is None:
ax = plt.gca()
the_mesh = ax.pcolormesh(
grid_x[:, level, :],
grid_h[:, level, :],
grid_bg[:, level, :],
cmap=cmap,
vmin=vmin,
vmax=vmax,
)
np.ma.sqrt(u**2 + w**2)
quiver_density_x = int((1 / dx) * quiver_spacing_x_km)
quiver_density_z = int((1 / dz) * quiver_spacing_z_km)
q = ax.quiver(
grid_x[::quiver_density_z, level, ::quiver_density_x],
grid_h[::quiver_density_z, level, ::quiver_density_x],
u[::quiver_density_z, level, ::quiver_density_x],
w[::quiver_density_z, level, ::quiver_density_x],
color="k",
scale=25.0 * quiverkey_len,
scale_units="width",
width=quiver_width,
)
quiver_font = {
"family": "sans-serif",
"style": "normal",
"variant": "normal",
"weight": "bold",
"size": "medium",
}
ax.quiverkey(
q,
qloc_x,
qloc_y,
quiverkey_len,
label=(str(quiverkey_len) + " m/s"),
fontproperties=quiver_font,
)
if colorbar_flag is True:
cp = grid_list[bg_grid_no][background_field].attrs["long_name"]
cp.replace(" ", "_")
cp = cp + " [" + grid_list[bg_grid_no][background_field].attrs["units"]
cp = cp + "]"
plt.colorbar(the_mesh, ax=ax, label=(cp))
if u_vel_contours is not None:
u_filled = np.ma.filled(u[:, level, :], fill_value=np.nan)
cs = ax.contour(
grid_x[:, level, :],
grid_h[:, level, :],
u_filled,
levels=u_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="U [m/s]")
if v_vel_contours is not None:
v_filled = np.ma.filled(w[:, level, :], fill_value=np.nan)
cs = ax.contour(
grid_x[:, level, :],
grid_h[:, level, :],
v_filled,
levels=v_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="V [m/s]")
if w_vel_contours is not None:
w_filled = np.ma.filled(w[:, level, :], fill_value=np.nan)
cs = ax.contour(
grid_x[:, level, :],
grid_h[:, level, :],
w_filled,
levels=w_vel_contours,
linewidths=2,
)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="W [m/s]")
if wind_vel_contours is not None:
vel = np.ma.sqrt(u[:, level, :] ** 2 + v[:, level, :] ** 2)
vel = vel.filled(fill_value=np.nan)
cs = ax.contour(
grid_x[:, level, :],
grid_h[:, level, :],
vel,
levels=wind_vel_contours,
linewidths=2,
)
cs.set_clim([np.min(wind_vel_contours), np.max(wind_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="|V| [m/s]")
if axes_labels_flag is True:
ax.set_xlabel(("X [km]"))
ax.set_ylabel(("Z [km]"))
if title_flag is True:
if grid_y[0, level, 0] > 0:
ax.set_title(
(
"PyDDA retreived winds @"
+ str(grid_y[0, level, 0])
+ " km north of origin."
)
)
else:
ax.set_title(
(
"PyDDA retreived winds @"
+ str(-grid_y[0, level, 0])
+ " km south of origin."
)
)
ax.set_xlim([grid_x.min(), grid_x.max()])
ax.set_ylim([grid_h.min(), grid_h.max()])
return ax
[docs]
def plot_yz_xsection_quiver(
Grids,
ax=None,
background_field="reflectivity",
level=1,
cmap="ChaseSpectral",
vmin=0,
vmax=70,
u_vel_contours=None,
v_vel_contours=None,
w_vel_contours=None,
wind_vel_contours=None,
u_field="u",
v_field="v",
w_field="w",
title_flag=True,
axes_labels_flag=True,
colorbar_flag=True,
colorbar_contour_flag=False,
bg_grid_no=0,
quiver_spacing_y_km=10.0,
quiver_spacing_z_km=1.0,
contour_alpha=0.7,
quiverkey_len=5.0,
quiverkey_loc="best",
quiver_width=0.01,
):
"""
This procedure plots a cross section of winds from wind fields
generated by PyDDA in the Y-Z plane using quivers.
The length of the quivers varies with wind speed.
Parameters
----------
Grids: list or DataTree
List of Py-DDA Grids to visualize
ax: matplotlib axis handle
The axis handle to place the plot on. Set to None to plot on the
current axis.
background_field: str
The name of the background field to plot the quivers on.
level: int
The number of the X level to plot the cross section through.
cmap: str or matplotlib colormap
The name of the matplotlib colormap to use for the background field.
vmin: float
The minimum bound to use for plotting the background field. None will
automatically detect the background field minimum.
vmax: float
The maximum bound to use for plotting the background field. None will
automatically detect the background field maximum.
u_vel_contours: 1-D array
The contours to use for plotting contours of u. Set to None to not
display such contours.
v_vel_contours: 1-D array
The contours to use for plotting contours of v. Set to None to not
display such contours.
w_vel_contours: 1-D array
The contours to use for plotting contours of w. Set to None to not
display such contours.
u_field: str
Name of zonal wind (u) field in Grids.
v_field: str
Name of meridional wind (v) field in Grids.
w_field: str
Name of vertical wind (w) field in Grids.
show_lobes: bool
If True, the dual doppler lobes from each pair of radars will be shown.
title_flag: bool
If True, PyDDA will generate a title for the plot.
axes_labels_flag: bool
If True, PyDDA will generate axes labels for the plot.
colorbar_flag: bool
If True, PyDDA will generate a colorbar for the plot background field.
colorbar_contour_flag: bool
If True, PyDDA will generate a colorbar for the contours.
bg_grid_no: int
Number of grid in Grids to take background field from.
Set to -1 to use maximum value from all grids.
quiver_spacing_y_km: float
Spacing in km between quivers in y axis.
quiver_spacing_z_km: float
Spacing in km between quivers in z axis.
contour_alpha: float
Alpha (transparency) of velocity contours. 0 = transparent, 1 = opaque
quiverkey_loc: str
Location of quiverkey. One of:
'best'
'top_left'
'top'
'top_right'
'bottom_left'
'bottom'
'bottom_right'
'left'
'right'
'top_left_outside'
'top_right_outside'
'bottom_left_outside'
'bottom_right_outside'
'best' will put the quiver key in the corner with the fewest amount of
valid data points while keeping the quiver key inside the plot.
The rest of the options will put the quiver key in that
particular part of the plot.
quiver_width: float
The width of each quiver for the plot, given as a fraction relative
to the plot width.
Returns
-------
ax: matplotlib axis
Axis handle to output axis
"""
if isinstance(Grids, DataTree):
child_list = list(Grids.children.keys())
grid_list = []
rad_names = []
for child in child_list:
if "radar" in child:
grid_list.append(Grids[child].to_dataset())
rad_names.append(child)
else:
grid_list = Grids
grid_bg = grid_list[bg_grid_no][background_field].values.squeeze()
grid_bg = np.ma.masked_invalid(grid_bg)
if vmin is None:
vmin = grid_bg.min()
if vmax is None:
vmax = grid_bg.max()
grid_h = grid_list[0]["point_altitude"].values / 1e3
grid_x = grid_list[0]["point_x"].values / 1e3
grid_y = grid_list[0]["point_y"].values / 1e3
dx = np.diff(grid_x, axis=2)[0, 0, 0]
dz = np.diff(grid_y, axis=1)[0, 0, 0]
u = grid_list[0][u_field].values.squeeze()
v = grid_list[0][v_field].values.squeeze()
w = grid_list[0][w_field].values.squeeze()
qloc_x, qloc_y = _parse_quiverkey_string(
quiverkey_loc,
grid_h[:, :, level],
grid_x[:, :, level],
grid_y[:, :, level],
grid_bg[:, :, level],
xsection="yz",
)
if ax is None:
ax = plt.gca()
the_mesh = ax.pcolormesh(
grid_y[:, :, level],
grid_h[:, :, level],
grid_bg[:, :, level],
cmap=cmap,
vmin=vmin,
vmax=vmax,
)
np.ma.sqrt(v**2 + w**2)
quiver_density_y = int((1 / dx) * quiver_spacing_y_km)
quiver_density_z = int((1 / dz) * quiver_spacing_z_km)
q = ax.quiver(
grid_y[::quiver_density_z, ::quiver_density_y, level],
grid_h[::quiver_density_z, ::quiver_density_y, level],
v[::quiver_density_z, ::quiver_density_y, level],
w[::quiver_density_z, ::quiver_density_y, level],
color="k",
cmap="coolwarm",
scale=25.0 * quiverkey_len,
scale_units="width",
width=quiver_width,
)
quiver_font = {
"family": "sans-serif",
"style": "normal",
"variant": "normal",
"weight": "bold",
"size": "medium",
}
ax.quiverkey(
q,
qloc_x,
qloc_y,
quiverkey_len,
label=(str(quiverkey_len) + " m/s"),
fontproperties=quiver_font,
)
if colorbar_flag is True:
cp = grid_list[bg_grid_no][background_field].attrs["long_name"]
cp.replace(" ", "_")
cp = cp + " [" + grid_list[bg_grid_no][background_field].attrs["units"]
cp = cp + "]"
plt.colorbar(the_mesh, ax=ax, label=(cp))
if u_vel_contours is not None:
u_filled = np.ma.filled(u[:, :, level], fill_value=np.nan)
cs = ax.contour(
grid_y[:, :, level],
grid_h[:, :, level],
u_filled,
levels=u_vel_contours,
linewidths=2,
)
plt.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="U [m/s]")
if v_vel_contours is not None:
v_filled = np.ma.filled(v[:, :, level], fill_value=np.nan)
cs = ax.contour(
grid_y[:, :, level],
grid_h[:, :, level],
v_filled,
levels=w_vel_contours,
linewidths=2,
)
plt.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="V [m/s]")
if w_vel_contours is not None:
w_filled = np.ma.filled(w[:, :, level], fill_value=np.nan)
cs = ax.contour(
grid_y[:, :, level],
grid_h[:, :, level],
w_filled,
levels=w_vel_contours,
linewidths=2,
)
plt.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="W [m/s]")
if wind_vel_contours is not None:
vel = np.ma.sqrt(u[:, :, level] ** 2 + v[:, :, level] ** 2)
vel = vel.filled(fill_value=np.nan)
cs = ax.contourf(
grid_y[:, :, level],
grid_h[:, :, level],
vel,
levels=wind_vel_contours,
linewidths=2,
)
cs.set_clim([np.min(wind_vel_contours), np.max(wind_vel_contours)])
cs.cmap.set_under(color="white", alpha=0)
cs.cmap.set_bad(color="white", alpha=0)
ax.clabel(cs)
if colorbar_contour_flag is True:
plt.colorbar(cs, ax=ax, label="|V| [m/s]")
if axes_labels_flag is True:
ax.set_xlabel(("Y [km]"))
ax.set_ylabel(("Z [km]"))
if title_flag is True:
if grid_x[0, 0, level] > 0:
ax.set_title(
(
"PyDDA retreived winds @"
+ str(grid_x[0, level, 0])
+ " km east of origin."
)
)
else:
ax.set_title(
(
"PyDDA retreived winds @"
+ str(-grid_x[0, level, 0])
+ " km west of origin."
)
)
ax.set_xlim([grid_y.min(), grid_y.max()])
ax.set_ylim([grid_h.min(), grid_h.max()])
return ax
def _parse_quiverkey_string(qloc, grid_z, grid_y, grid_x, grid_bg, xsection="xy"):
"""
This is a private method for parsing the quiverkey location string.
This is only used internally by PyDDA.
"""
if qloc == "best":
# Get top left corner of grid
if xsection == "xy":
top_y = grid_y.max() - (grid_y.max() - grid_y.min()) / 10
right_x = grid_x.max() - (grid_x.max() - grid_x.min()) / 10
bottom_y = grid_y.min() + (grid_y.max() - grid_y.min()) / 10
left_x = grid_x.min() + (grid_x.max() - grid_x.min()) / 10
mid_y1 = grid_y.max() - 4 * (grid_y.max() - grid_y.min()) / 10
mid_x1 = grid_x.max() - 4 * (grid_x.max() - grid_x.min()) / 10
mid_y2 = grid_y.max() - 6 * (grid_y.max() - grid_y.min()) / 10
mid_x2 = grid_x.max() - 6 * (grid_x.max() - grid_x.min()) / 10
top_right = len(
np.where(
grid_bg[np.logical_and(grid_y > top_y, grid_x > right_x)].mask
is True
)[0]
)
top_left = len(
np.where(
grid_bg[np.logical_and(grid_y > top_y, grid_x < left_x)].mask
is True
)[0]
)
bot_left = len(
np.where(
grid_bg[np.logical_and(grid_y < bottom_y, grid_x < left_x)].mask
is True
)[0]
)
bot_right = len(
np.where(
grid_bg[np.logical_and(grid_y > bottom_y, grid_x > right_x)].mask
is True
)[0]
)
top_inds = np.logical_and.reduce(
(grid_y > top_y, grid_x > mid_x2, grid_x < mid_x1)
)
bottom_inds = np.logical_and.reduce(
(grid_y < bottom_y, grid_x > mid_x2, grid_x < mid_x1)
)
top = len(np.where(grid_bg[top_inds].mask is True)[0])
bottom = len(np.where(grid_bg[bottom_inds].mask is True)[0])
left_inds = np.logical_and.reduce(
(grid_x < left_x, grid_y > mid_y2, grid_y < mid_y1)
)
left = len(np.where(grid_bg[left_inds].mask is True)[0])
right_inds = np.logical_and.reduce(
(grid_x > right_x, grid_y > mid_y2, grid_y < mid_y1)
)
right = len(np.where(grid_bg[right_inds].mask is True)[0])
elif xsection == "xz":
top_z = grid_z.max() - (grid_z.max() - grid_y.min()) / 10
right_x = grid_x.max() - (grid_x.max() - grid_x.min()) / 10
bottom_z = grid_z.min() + (grid_z.max() - grid_z.min()) / 10
left_x = grid_x.min() + (grid_x.max() - grid_x.min()) / 10
mid_z1 = grid_z.max() - 4 * (grid_z.max() - grid_z.min()) / 10
mid_x1 = grid_x.max() - 4 * (grid_x.max() - grid_x.min()) / 10
mid_z2 = grid_z.max() - 6 * (grid_z.max() - grid_z.min()) / 10
mid_x2 = grid_x.max() - 6 * (grid_x.max() - grid_x.min()) / 10
top_right = len(
np.where(
grid_bg[np.logical_and(grid_z > top_z, grid_x > right_x)].mask
is True
)[0]
)
top_left = len(
np.where(
grid_bg[np.logical_and(grid_z > top_z, grid_x < left_x)].mask
is True
)[0]
)
bot_left = len(
np.where(
grid_bg[np.logical_and(grid_z < bottom_z, grid_x < left_x)].mask
is True
)[0]
)
bot_right = len(
np.where(
grid_bg[np.logical_and(grid_z > bottom_z, grid_x > right_x)].mask
is True
)[0]
)
top_inds = np.logical_and.reduce(
(grid_z > top_z, grid_x > mid_x2, grid_x < mid_x1)
)
bottom_inds = np.logical_and.reduce(
(grid_z < bottom_z, grid_x > mid_x2, grid_x < mid_x1)
)
top = len(np.where(grid_bg[top_inds].mask is True)[0])
bottom = len(np.where(grid_bg[bottom_inds].mask is True)[0])
left_inds = np.logical_and.reduce(
(grid_x < left_x, grid_z > mid_z2, grid_z < mid_z1)
)
left = len(np.where(grid_bg[left_inds].mask is True)[0])
right_inds = np.logical_and.reduce(
(grid_x > right_x, grid_z > mid_z2, grid_z < mid_z1)
)
right = len(np.where(grid_bg[right_inds].mask is True)[0])
elif xsection == "yz":
top_z = grid_z.max() - (grid_z.max() - grid_y.min()) / 10
right_y = grid_y.max() - (grid_y.max() - grid_y.min()) / 10
bottom_z = grid_z.min() + (grid_z.max() - grid_z.min()) / 10
left_y = grid_y.min() + (grid_y.max() - grid_y.min()) / 10
mid_z1 = grid_z.max() - 4 * (grid_z.max() - grid_z.min()) / 10
mid_y1 = grid_y.max() - 4 * (grid_y.max() - grid_y.min()) / 10
mid_z2 = grid_z.max() - 6 * (grid_z.max() - grid_z.min()) / 10
mid_y2 = grid_y.max() - 6 * (grid_y.max() - grid_y.min()) / 10
top_right = len(
np.where(
grid_bg[np.logical_and(grid_z > top_z, grid_y > right_y)].mask
is True
)[0]
)
top_left = len(
np.where(
grid_bg[np.logical_and(grid_z > top_z, grid_y < left_y)].mask
is True
)[0]
)
bot_left = len(
np.where(
grid_bg[np.logical_and(grid_z < bottom_z, grid_y < left_y)].mask
is True
)[0]
)
bot_right = len(
np.where(
grid_bg[np.logical_and(grid_z > bottom_z, grid_y > right_y)].mask
is True
)[0]
)
top_inds = np.logical_and.reduce(
(grid_z > top_z, grid_y > mid_y2, grid_x < mid_y1)
)
bottom_inds = np.logical_and.reduce(
(grid_z < bottom_z, grid_y > mid_y2, grid_y < mid_y1)
)
top = len(np.where(grid_bg[top_inds].mask is True)[0])
bottom = len(np.where(grid_bg[bottom_inds].mask is True)[0])
left_inds = np.logical_and.reduce(
(grid_y < left_y, grid_z > mid_z2, grid_z < mid_z1)
)
left = len(np.where(grid_bg[left_inds].mask is True)[0])
right_inds = np.logical_and.reduce(
(grid_y > right_y, grid_z > mid_z2, grid_z < mid_z1)
)
right = len(np.where(grid_bg[right_inds].mask is True)[0])
loc_array = np.array(
[top_right, top_left, bot_left, bot_right, top, bottom, left, right]
)
if loc_array.max() == top_right:
qloc_x = 0.9
qloc_y = 0.93
elif loc_array.max() == top_left:
qloc_x = 0.07
qloc_y = 0.9
elif loc_array.max() == bot_left:
qloc_x = 0.07
qloc_y = 0.05
elif loc_array.max() == bot_right:
qloc_x = 0.9
qloc_y = 0.05
elif loc_array.max() == top:
qloc_x = 0.45
qloc_y = 0.93
elif loc_array.max() == bottom:
qloc_x = 0.45
qloc_y = 0.05
elif loc_array.max() == left:
qloc_x = 0.07
qloc_y = 0.45
elif loc_array.max() == right:
qloc_x = 0.9
qloc_y = 0.45
elif qloc == "top_left":
qloc_x = 0.07
qloc_y = 0.93
elif qloc == "top_right":
qloc_x = 0.9
qloc_y = 0.93
elif qloc == "bottom_left":
qloc_x = 0.07
qloc_y = 0.05
elif qloc == "bottom_right":
qloc_x = 0.9
qloc_y = 0.05
elif qloc == "top":
qloc_x = 0.45
qloc_y = 0.93
elif qloc == "bottom":
qloc_x = 0.45
qloc_y = 0.05
elif qloc == "left":
qloc_x = 0.07
qloc_y = 0.45
elif qloc == "right":
qloc_x = 0.9
qloc_y = 0.45
elif qloc == "top_left_outside":
qloc_x = 0.07
qloc_y = 1.01
elif qloc == "top_right_outside":
qloc_x = 0.9
qloc_y = 1.01
elif qloc == "bottom_left_outside":
qloc_x = 0.07
qloc_y = -0.09
elif qloc == "bottom_right_outside":
qloc_x = 0.9
qloc_y = -0.09
else:
raise ValueError("Invalid quiver key location!")
return qloc_x, qloc_y