"""
Skew-T Plotting Library
"""
import numpy as np
from scipy.integrate import odeint
import matplotlib
import pylab
import matplotlib.transforms as transforms
from matplotlib.scale import LogScale
from matplotlib.patches import Circle, Polygon, Path
from matplotlib.ticker import NullLocator
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from mpl_toolkits.axes_grid1.mpl_axes import Axes as InsetAxes
import copy
from pycaps.derive import moist_lapse
from pycaps.derive import temp_from_vapr
from pycaps.derive import sb_parcel
from pycaps.derive import compute_lcl
def _malr(temp, pres):
return moist_lapse(t=temp, p=pres)
[docs]class SkewTPlotter(object):
"""
Plot Skew-T log-p diagrams
Attributes:
p_lines (list): A list of the pressure lines (in hPa) to plot.
p_color (str): The color of the pressure lines
T_plot_min (float): Minimum isotherm on the plot in degrees C.
T_step (float): The temperature difference (in degrees C) between adjacent isotherms.
T_color_below_0 (str): The color of isotherms below 0 C
T_color_0 (str): The color of the 0 C isotherm
T_color_above_0 (str): The color of isotherms above 0 C
th_min (float): Minimum dry adiabat on the plot in degrees C (references temperature at 1000 mb).
th_max (float): Maximum dry adiabat on the plot in degrees C (references temperature at 1000 mb).
th_step (float): Temperature difference (in degrees C) between adjacent dry adiabats.
th_color (str): The color of the dry adiabats.
the_min (float): Minimum moist adiabat on the plot in degrees C (references temperature at 1000 mb).
the_max (float): Maximum moist adiabat on the plot in degrees C (references temeprature at 1000 mb).
the_step (float): Temperature difference (in degrees C) between adjacent moist adiabats.
the_color (float): The color of the moist adiabats.
w_lines (list): A list of mixing ratio lines (in g/kg) to plot.
w_min_p (float): The starting pressure (in hPa) for mixing ratio lines.
w_max_p (float): The ending pressure (in hPa) for mixing ratio lines.
w_color (str): The color of the mixing ratio lines.
"""
def __init__(self):
"""
Constructor for the SkewTPlotter class
"""
self._initialized = False
self._stlp_data_transform, self._stlp_xlabel_transform, self._stlp_ylabel_transform = None, None, None
self._p_min, self._p_max = 100, 1050
self.p_lines = [1000, 925, 850, 700, 500, 300, 200, 100]
self.p_color = 'k'
self.T_plot_min, self._T_min, self._T_max, self.T_step = -110, -30, 50, 10
self.T_color_below_0, self.T_color_0, self.T_color_above_0 = 'c', 'b', '#880000'
self.th_min, self.th_max, self.th_step = -30, 200, 10
self.th_color = '#802a2a'
self.the_min, self.the_max, self.the_step = -20, 50, 5
self.the_color = '#2e8b57'
self.w_lines = [1, 2, 3, 5, 8, 12, 16, 20]
self.w_min_p, self.w_max_p = 600, 1000
self.w_color = '#7fff00'
self._initialized = True
def _build_transform(self):
"""
Builds the skew-T log-p transform used in plotting. [private]
"""
current_axes = pylab.gca()
data_figure_trans = current_axes.transData + current_axes.transAxes.inverted()
pylab.xlim((self.T_min, self.T_max))
pylab.ylim((self.p_min, self.p_max))
identity_matrix = np.zeros((3, 3))
for idx in range(3): identity_matrix[idx, idx] = 1
# Create the affine matrix for the skew transform. This only works in data coordinates. We'll fix that later ...
skew_matrix = np.copy(identity_matrix)
skew_matrix[0, 1] = np.tan(45 * np.pi / 180)
skew_transform = transforms.Affine2D(skew_matrix)
# Create the logarithmic transform in the y.
log_p_transform = transforms.blended_transform_factory(transforms.Affine2D(),
LogScale(current_axes.yaxis, basey=10).get_transform())
# The log transform shrinks everything to log(p) space, so define a scale factor to blow it back up to a
# reasonable size. Also, need to scale everything to -30 to 50 C in T and 1000 to 100 hPa in p to get the
# skewing to work properly.
p_bnd_pre_trans = log_p_transform.transform(np.array([[0, 100], [0, 1000]]))[:, 1]
pre_p_scale_factor = (self.p_max - self.p_min) / (p_bnd_pre_trans[1] - p_bnd_pre_trans[0])
pre_T_scale_factor = (self.T_max - self.T_min) / (50. - -30.)
# This will scale it to the user's requested size.
p_bnd_post_trans = log_p_transform.transform(np.array([[0, self.p_min], [0, self.p_max]]))[:, 1]
post_p_scale_factor = (self.p_max - self.p_min) / (p_bnd_post_trans[1] - p_bnd_post_trans[0])
# Define the affine transform for the flip and another for the scale back to reasonable coordinates after the
# log transform.
flip_transform = transforms.Affine2D.identity().scale(1, -1)
preskew_scale_transform = transforms.Affine2D().translate(30, p_bnd_pre_trans[1])\
.scale(pre_T_scale_factor, pre_p_scale_factor).translate(self.T_min, self.p_min)
postskew_scale_transform = preskew_scale_transform.inverted().translate(0, p_bnd_post_trans[1])\
.scale(1, post_p_scale_factor).translate(0, self.p_min)
# Define a transform that basically does everything but the skew so we can figure out where the 1000 mb level is
# and skew around that line.
prelim_data_transform = log_p_transform + flip_transform + preskew_scale_transform + data_figure_trans
marker = prelim_data_transform.transform(np.array([[self.T_min, 1000]]))[0, 1]
# Add a translation to that marker point into the data-figure transform matrix.
data_figure_trans += transforms.Affine2D().translate(0, -marker)
# Define our skew transform in figure coordinates.
figure_skew_transform = data_figure_trans + skew_transform + data_figure_trans.inverted()
# Create our skew-T log-p transform matrix. It does the log-p transform first, then the flip, then the scale,
# then the skew.
self._stlp_data_transform = log_p_transform + flip_transform + preskew_scale_transform + \
figure_skew_transform + postskew_scale_transform + current_axes.transData
# Create a blended transform where the y axis is the log-p, but the x axis is the axes transform (for adding
# pressure labels and wind barbs).
self._stlp_xlabel_transform = \
transforms.blended_transform_factory(self._stlp_data_transform, current_axes.transAxes)
self._stlp_ylabel_transform = \
transforms.blended_transform_factory(current_axes.transAxes, self._stlp_data_transform)
return
@property
def trans(self):
"""
The skew-T log-p transform as a matplotlib transform object.
"""
if self._stlp_data_transform is None:
self._build_transform()
return self._stlp_data_transform
@property
def p_min(self):
"""
float: The minimum pressure of the plot in hPa
"""
return self._p_min
@p_min.setter
def p_min(self, p_min):
"""
Setter for the plot's minimum pressure attribute
Args:
p_min (float): The new minimum pressure in hPa
"""
self._p_min = p_min
if self._initialized:
self._build_transform()
@property
def p_max(self):
"""
float: The maximum pressure of the plot in hPa
"""
return self._p_max
@p_max.setter
def p_max(self, p_max):
"""
Setter for the plot's maximum pressure attribute
Args:
p_max (float): The new maximum pressure in hPa
"""
self._pmax = p_max
if self._initialized:
self._build_transform()
@property
def T_min(self):
"""
float: The minimum temperature of the plot (at 1000 hPa) in degrees C.
"""
return self._T_min
@T_min.setter
def T_min(self, T_min):
"""
Setter for the minimum temperature of the plot (at 1000 hPa).
Args:
T_min: The new minimum temperature at 1000 hPa for the plot in degrees C.
"""
self._T_min = T_min
if self._initialized:
self._build_transform()
@property
def T_max(self):
"""
float: The maximum temperature of the plot (at 1000 hPa) in degrees C.
"""
return self._T_max
@T_max.setter
def T_max(self, T_max):
"""
Setter for the maximum temperature of the plot (at 1000 hPa).
Args:
T_max: The new maximum temperature of the plot at 1000 hPa in degrees C.
"""
self._T_max = T_max
if self._initialized:
self._build_transform()
[docs] def plot_skewt_background(self, staff=True, clip=False):
"""
Plots the background for the Skew-T diagram. Must be called first to initialize the plotting.
Args:
staff (bool): If True, draw a line for the wind staff. Default is True.
clip (bool): If True, clip the top right corner of the plot. Default is False.
"""
current_axes = pylab.gca()
current_figure = current_axes.figure
current_axes.xaxis.set_major_locator(NullLocator())
current_axes.xaxis.set_minor_locator(NullLocator())
current_axes.yaxis.set_major_locator(NullLocator())
current_axes.yaxis.set_minor_locator(NullLocator())
if self._stlp_data_transform is None:
self._build_transform()
_stlp_axes_transform = self._stlp_data_transform + current_axes.transAxes.inverted()
linekw = {
'transform':self._stlp_data_transform,
'zorder':999,
}
plinekw = linekw
if clip:
for s in ['top', 'bottom', 'left', 'right']:
current_axes.spines[s].set_visible(False)
(_dummy, clip_p_1), (clip_T_2, clip_p_2) = _stlp_axes_transform.transform(np.array([[0, 600], [0, 400]]))
clip_region = np.array([[0., 0.], [1., 0.], [1., clip_p_1], [clip_T_2, clip_p_2], [clip_T_2, 1.], [0., 1.]])
clip_neg = np.array([[1., 1.], [1., clip_p_1], [clip_T_2, clip_p_2], [clip_T_2, 1.]])
clip_patch = Polygon(clip_region, ec='k', fc='none', lw=1.5, transform=current_axes.transAxes, zorder=1001)
clip_neg_patch = Polygon(clip_neg, ec='w', fc='w', transform=current_axes.transAxes, zorder=1000)
plinekw = copy.copy(linekw)
linekw['clip_path'] = (Path(clip_region), current_axes.transAxes)
linekw['clip_on'] = True
# Draw the isobars.
for p_line in self.p_lines:
pylab.plot([self.T_plot_min, self.T_max], [p_line, p_line], color='k', lw=0.75, **plinekw)
ax_x, ax_y = \
(self._stlp_ylabel_transform + current_axes.transAxes.inverted()).transform(np.array([[0, p_line]]))[0]
if 0 <= ax_y <= 1:
pylab.text(-0.01, p_line, "%d" % p_line,
va='center', ha='right', transform=self._stlp_ylabel_transform, zorder=999)
# Draw the isotherms.
for T_line in range(self.T_plot_min, self.T_max + self.T_step, self.T_step):
if T_line < 0:
weight = 0.75
color = self.T_color_below_0
elif T_line > 0:
color = self.T_color_above_0
weight = 0.75
else:
weight = 1
color = self.T_color_0
pylab.plot([T_line, T_line], [self.p_min, self.p_max], color=color, lw=weight, **linekw)
txt = "%d" if T_line < 0 else " %d"
offset = (self.p_max - self.p_min) / (1050. - 100.) * 25
if self.T_min <= T_line <= self.T_max:
pylab.text(T_line, self.p_max + offset, txt % T_line,
ha='center', va='top', transform=self._stlp_data_transform, zorder=999)
# Draw the dry adiabats.
for th_line in range(self.th_min, self.th_max + self.th_step, self.th_step):
p_th = np.arange(self.p_min, self.p_max + 10., 10.)
t_th = (th_line + 273.15) * (p_th / 1000.) ** (2. / 7.) - 273.15
pylab.plot(t_th, p_th, color=self.th_color, lw=0.75, **linekw)
for w_line in self.w_lines:
p_w = np.arange(self.w_min_p, self.w_max_p + 10., 10.)
td_w = temp_from_vapr(qv=w_line, p=p_w)
pylab.plot(td_w, p_w, color=self.w_color, lw=0.75, linestyle='--', **linekw)
ax_x, ax_y = _stlp_axes_transform.transform(np.array([[td_w[0], p_w[0]]]))[0]
if 0 < ax_x < 1 and 0 < ax_y < 1:
offset = (self.p_max - self.p_min) / (1050. - 100.) * 15
p_lbl = p_w[0] - offset
td_lbl = temp_from_vapr(qv=w_line, p=p_lbl)
pylab.text(td_lbl, p_lbl, "%d" % w_line,
color=self.w_color, ha='center', transform=self._stlp_data_transform, zorder=999,
fontsize=(matplotlib.rcParams['font.size'] * 0.6))
# Draw the moist adiabats.
for the_line in xrange(self.the_min, self.the_max + self.the_step, self.the_step):
p_the_above = np.arange(self.p_min, 1010., 10.)
p_the_below = np.arange(1000., self.p_max + 10., 10.)
t_the_above = odeint(_malr, the_line + 273.15, 100 * p_the_above[::-1])[::-1, 0] - 273.15
t_the_below = odeint(_malr, the_line + 273.15, 100 * p_the_below)[:, 0] - 273.15
p_the = np.concatenate((p_the_above, p_the_below[1:]))
t_the = np.concatenate((t_the_above, t_the_below[1:]))
pylab.plot(t_the, p_the, color=self.the_color, lw=0.75, linestyle=':', **linekw)
if staff:
# Draw wind staff
t_base = self.T_max - 5
pylab.axvline(x=t_base, color='k', lw=1., zorder=1002)
if clip:
pylab.gca().add_patch(clip_patch)
pylab.gca().add_patch(clip_neg_patch)
pylab.xlim((self.T_min, self.T_max))
pylab.ylim((self.p_min, self.p_max))
pylab.xlabel("Temperature ($^{\circ}$C)", labelpad=15)
pylab.ylabel("Pressure (mb)", labelpad=35)
# Make the resolution that savefig() uses be the same as the figure resolution (weird figures result if we don't
# do this).
matplotlib.rcParams['savefig.dpi'] = pylab.gcf().dpi
return
[docs] def plot_profile(self, t_snd, p_snd, **kwargs):
"""
Plot a temperature or dewpoint profile on the skew-T
Args:
t_snd (np.array): A 1-dimensional array containing the profile temperature in degrees C.
p_snd (np.array): A 1-dimensional array containing the profile pressure in hPa.
**kwargs: Any keywords that can be passed to matplotlib.pyplot.plot() (see
http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot; note that if you pass in a transform,
it will be overwritten).
"""
kwargs['transform'] = self._stlp_data_transform
line = pylab.plot(t_snd, p_snd, **kwargs)
pylab.xlim((self.T_min, self.T_max))
pylab.ylim((self.p_min, self.p_max))
return line
[docs] def plot_winds(self, u_snd, v_snd, p_snd, **kwargs):
"""
Plot a wind profile on the right side of the skew-T
Args:
u_snd (np.array): A 1-dimensional array containing the u wind profile
v_snd (np.array): A 1-dimensional array containing the v wind profile
p_snd (np.array): A 1-dimenionsal array containing the pressure profile
**kwargs: Any keywords that can be passed to matplotlib.pyplot.barbs() (see
http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.barbs; note that if you pass in clip_on or
z_order, it will be overwritten)
"""
t_base = self.T_max - 5
stride = 1
idx_min = np.argmin(np.abs(self.p_max - p_snd))
idx_max = np.argmin(np.abs(self.p_min - p_snd))
if p_snd[idx_min] > self.p_max:
idx_min += 1
if p_snd[idx_max] > self.p_min:
idx_max += 1
thin_data = (slice(idx_min, idx_max, stride))
points = np.vstack((t_base * np.ones(p_snd.shape), p_snd)).T
x, base_y = self._stlp_data_transform.transform(np.array([[t_base, 1000]]))[0]
trans_points = self._stlp_data_transform.transform(points)
trans_points[:, 0] = x
axes_points = pylab.gca().transData.inverted().transform(trans_points)
kwargs['clip_on'] = False
kwargs['zorder'] = 1002
plot_xs, plot_ys = zip(*axes_points)
pylab.barbs(plot_xs[thin_data], plot_ys[thin_data], u_snd[thin_data], v_snd[thin_data], **kwargs)
return
def _plot_hodograph_background(self, max_ring, step_ring=10):
pylab.gca().axes.get_xaxis().set_visible(False)
pylab.gca().axes.get_yaxis().set_visible(False)
for w_spd in xrange(step_ring, max_ring + step_ring, step_ring):
ring = Circle((0, 0), w_spd, edgecolor='#444444', facecolor='none', linestyle='dashed', zorder=999)
pylab.gca().add_patch(ring)
if not (w_spd % (step_ring * 2)):
pylab.text(w_spd, -2, "%d" % w_spd,
color='#444444', ha='center', va='top', clip_on=True, clip_box=pylab.gca().get_position(),
zorder=999)
pylab.axvline(x=0, color='k', zorder=999)
pylab.axhline(y=0, color='k', zorder=999)
return
[docs] def plot_hodograph(self, u_snd, v_snd, p_snd, clip_p=200., bounds=(-20, -20, 40, 40), **kwargs):
"""
Plot a hodograph inset in the upper-right corner
Args:
u_snd (np.array): A 1-dimensional array containing the u wind component
v_snd (np.array): A 1-dimensional array containing the v wind component
p_snd (np.array): A 1-dimensional array containing the pressure in hPa.
clip_p (float): The pressure in hPa at which to cut off the line on the hodograph. The default is 200 hPa.
bounds (tuple): The bounds for the hodograph as a tuple (u_min, v_min, u_max, v_max). The default is (-20, -20,
40, 40)
**kwargs: Any keywords that can be passed to matplotlib.pyplot.plot() (see
http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot).
"""
u_min, v_min, u_max, v_max = bounds
if not any(isinstance(a, InsetAxes) for a in pylab.gcf().axes):
par_ax = pylab.gca()
ax_width_inches, ax_height_inches = \
[f * a for f, a in zip(pylab.gcf().get_size_inches(), par_ax.get_position().size)]
ins_ax = inset_axes(pylab.gca(), 0.4 * ax_width_inches, 0.4 * ax_width_inches, loc=1)
pylab.sca(ins_ax)
max_speed = np.hypot(max(abs(u_min), abs(u_max)), max(abs(v_min), abs(v_max)))
self._plot_hodograph_background(int(max_speed))
else:
ins_ax = [a for a in pylab.gcf().axes if isinstance(a, InsetAxes)][0]
par_ax = [a for a in pylab.gcf().axes if not isinstance(a, InsetAxes)][0]
pylab.sca(ins_ax)
pylab.plot(u_snd[p_snd >= clip_p], v_snd[p_snd >= clip_p], **kwargs)
pylab.xlim((u_min, u_max))
pylab.ylim((v_min, v_max))
pylab.sca(par_ax)
return
[docs]def compute_parcel_trace(t_snd, td_snd, p_snd, prc_func=sb_parcel):
"""
Compute the parcel temperature trace
Args:
t_snd (np.array): A 1-dimensional float array containing the environmental temperature
td_snd (np.array): A 1-dimensional float array containing the environmental dewpoint
p_snd (np.array): A 1-dimensional float array containing the pressure
prc_func (function): A function that returns the initial parcel temperature, pressure, and dewpoint (in that
order) given the sounding temperature, pressure, and dewpoint (in that order). Some pre-defined functions
are available in derived.parcel.
Returns:
A 1-dimensional array containing the parcel temperature trace.
"""
t_prc, p_prc, td_prc = prc_func(t_snd, p_snd, td_snd)
ub_p = np.argmin(np.abs(p_snd - p_prc))
if p_snd[ub_p] < p_prc:
ub_p -= 1
lcl_temp, lcl_pres = compute_lcl(t_prc, p_prc, sfc_td=td_prc)
lcl_idx = np.argmin(np.abs(lcl_pres - p_snd))
if p_snd[lcl_idx] > lcl_pres:
lcl_idx += 1
dry_adiabatic_pres = np.append(p_snd[ub_p:lcl_idx], lcl_pres)
moist_adiabatic_pres = np.append([lcl_pres], p_snd[lcl_idx:])
dry_adiabatic_temp = (t_prc + 273.15) * (dry_adiabatic_pres / p_prc) ** (2. / 7.) - 273.15
moist_adiabatic_temp = odeint(_malr, lcl_temp + 273.15, 100 * moist_adiabatic_pres)[:, 0] - 273.15
return np.append(dry_adiabatic_temp[:-1], moist_adiabatic_temp), \
np.append(dry_adiabatic_pres[:-1], moist_adiabatic_pres)
[docs]def plot_sounding(u, v, p, t, td, file_name, hodograph=False, hodo_bounds=(-20, -20, 40, 40)):
"""
Plot a single sounding, optionally with a hodograph, and save to a file
Args:
u (np.array): A 1-dimensional array containing the u wind component
v (np.array): A 1-dimensional array containing the v wind component
p (np.array): A 1-dimensional array containing the pressure in hPa
t (np.array): A 1-dimensional array containing the temperature in degrees C
td (np.array): A 1-dimensional array containing the dewpoint in degrees C
file_name (str): The name of the file in which to save the image
hodograph (bool): If True, plot a hodograph in the upper-right corner
hodo_bounds (tuple): Bounds on the hodograph as a tuple (u_min, v_min, u_max, v_max), if one is plotted. The
default is (-20, -20, 40, 40)
"""
pylab.figure(figsize=(8, 10), dpi=150)
stp = SkewTPlotter()
stp.plot_skewt_background()
stp.plot_profile(t, p, color='r', lw=1.5)
stp.plot_profile(td, p, color='g', lw=1.5)
stp.plot_winds(u, v, p)
if hodograph:
stp.plot_hodograph(u, v, p, bounds=hodo_bounds)
pylab.savefig(file_name)
pylab.close()
return