Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create Advanced_Sounding_With_Complex_Layout.py #3204

Merged
merged 8 commits into from
Oct 10, 2023
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,18 @@
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause


"""
==================================
Sounding Plot with Complex Layout
==================================
==========================================
Advanced Sounding Plot with Complex Layout
==========================================

This example combines simple MetPy plotting functionality, `metpy.calc`
kylejgillett marked this conversation as resolved.
Show resolved Hide resolved
computation functionality, and a few basic tricks to create an
advanced sounding plotter with a clean layout & high readability.
computation functionality, and a few basic Matplotlib tricks to create
an advanced sounding plot with a complex layout & high readability.
"""

# First let's start with some simple imports
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
Expand All @@ -24,36 +26,29 @@
###########################################
# Upper air data can easily be obtained using the siphon package,
# but for this example we will use some of MetPy's sample data.

col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']

df = pd.read_fwf(get_test_data('may4_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)

# Drop any rows with all NaN values for T, Td, winds
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'),
how='all').reset_index(drop=True)

###########################################
# We will pull the data out of the example dataset into
# individual variables and assign units.

p = df['pressure'].values * units.hPa
z = df['height'].values * units.m
T = df['temperature'].values * units.degC
Td = df['dewpoint'].values * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)

###########################################
# Now lets make a Skew-T Log-P diagram using some simply
# Now let's make a Skew-T Log-P diagram using some simply
# MetPy functionality

# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
add_metpy_logo(fig, 430, 30, size='large')

add_metpy_logo(fig, 90, 80, size='small')
skew = SkewT(fig, rotation=45, rect=(0.1, 0.1, 0.55, 0.85))

# Plot the data using normal plotting functions, in this case using
Expand All @@ -77,22 +72,20 @@
h = Hodograph(ax, component_range=60.)
h.add_grid(increment=20)
h.plot(u, v)

###########################################
# This layout isn't bad, especially for how little code it required,
# This layout isn't bad, especially for how little code it requires,
# but we could add a few simple tricks to greatly increase the
# readability and complexity of our Skew-T/Hodograph layout. Lets
# readability and complexity of our Skew-T/Hodograph layout. Let's
# try another Skew-T with a few more advanced features:
kylejgillett marked this conversation as resolved.
Show resolved Hide resolved

###########################################
# STEP 1: CREATE THE SKEW-T OBJECT AND MODIFY IT TO CREATE A
# NICE, CLEAN PLOT

# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(18, 12))
skew = SkewT(fig, rotation=45, rect=(0, 0, 0.50, 0.90))
# add the metpy logo
add_metpy_logo(fig, 100, 80, size='small')
skew = SkewT(fig, rotation=45, rect=(0.05, 0.05, 0.50, 0.90))

# add the Metpy logo
add_metpy_logo(fig, 105, 85, size='small')

# Change to adjust data limits and give it a semblance of what we want
skew.ax.set_adjustable('datalim')
Expand All @@ -103,7 +96,7 @@
skew.ax.set_xlabel(str.upper(f'Temperature ({T.units:~P})'), weight='bold')
skew.ax.set_ylabel(str.upper(f'Pressure ({p.units:~P})'), weight='bold')

# Set the facecolor of the Skew Object and the Figure to white
# Set the facecolor of the skew-t object and the figure to white
fig.set_facecolor('#ffffff')
skew.ax.set_facecolor('#ffffff')

Expand All @@ -115,18 +108,16 @@
for i in range(0, 8):
skew.shade_area(y=y, x1=x1[i], x2=x2[i], color='gray', alpha=0.02, zorder=1)

###########################################
# STEP 2: PLOT DATA ON THE SKEW-T. TAKE A COUPLE EXTRA STEPS TO
# INCREASE READABILITY

# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
# set the linewidth to 4 for increased readability.
# We will also add the 'label' kew word argument for our legend.
# Set the linewidth to 4 for increased readability.
# We will also add the 'label' keyword argument for our legend.
skew.plot(p, T, 'r', lw=4, label='TEMPERATURE')
skew.plot(p, Td, 'g', lw=4, label='DEWPOINT')

# again we can use some simple python math functionality to 'resample'
# Again we can use some simple Python math functionality to 'resample'
# the wind barbs for a cleaner output with increased readability.
# Something like this would work.
interval = np.logspace(2, 3, 40) * units.hPa
Expand All @@ -135,20 +126,19 @@

# Add the relevant special lines native to the Skew-T Log-P diagram &
# provide basic adjustments to linewidth and alpha to increase readability
# first we add a matplotlib axvline to highlight the 0 degree isotherm
# first, we add a matplotlib axvline to highlight the 0-degree isotherm
skew.ax.axvline(0 * units.degC, linestyle='--', color='blue', alpha=0.3)
skew.plot_dry_adiabats(lw=1, alpha=0.3)
skew.plot_moist_adiabats(lw=1, alpha=0.3)
skew.plot_mixing_lines(lw=1, alpha=0.3)

# Calculate LCL height and plot as black dot. Because `p`'s first value is
# Calculate LCL height and plot as a black dot. Because `p`'s first value is
# ~1000 mb and its last value is ~250 mb, the `0` index is selected for
# `p`, `T`, and `Td` to lift the parcel from the surface. If `p` was inverted,
# i.e. start from low value, 250 mb, to a high value, 1000 mb, the `-1` index
# i.e. start from a low value, 250 mb, to a high value, 1000 mb, the `-1` index
# should be selected.
lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')

# Calculate full parcel profile and add to plot as black line
prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2, label='SB PARCEL PATH')
Expand All @@ -157,20 +147,20 @@
skew.shade_cin(p, T, prof, Td, alpha=0.2, label='SBCIN')
skew.shade_cape(p, T, prof, alpha=0.2, label='SBCAPE')

###########################################
# STEP 3: CREATE THE HODOGRAPH INSET. TAKE A FEW EXTRA STEPS TO
# INCREASE READABILITY

# Create a hodograph object: first we need to add an axis
# then we can create the metpy Hodograph
hodo_ax = plt.axes((0.43, 0.40, 0.5, 0.5))
# then we can create the Metpy Hodograph
hodo_ax = plt.axes((0.48, 0.45, 0.5, 0.5))
h = Hodograph(hodo_ax, component_range=80.)

# Add two separate grid increments for a cooler look. This also
# helps to increase readability
h.add_grid(increment=20, ls='-', lw=1.5, alpha=0.5)
h.add_grid(increment=10, ls='--', lw=1, alpha=0.2)

# The next few steps makes for a clean hodograph inset, removing the
# tick marks, tick labels and axis labels
# tick marks, tick labels, and axis labels
h.ax.set_box_aspect(1)
h.ax.set_yticklabels([])
h.ax.set_xticklabels([])
Expand All @@ -179,7 +169,7 @@
h.ax.set_xlabel(' ')
h.ax.set_ylabel(' ')

# Here we can add a simple python for loop that adds tick marks
# Here we can add a simple Python for loop that adds tick marks
# to the inside of the hodograph plot to increase readability!
plt.xticks(np.arange(0, 0, 1))
plt.yticks(np.arange(0, 0, 1))
Expand All @@ -205,22 +195,19 @@
alpha=0.2, label='Bunkers RM Vector',
length_includes_head=True, head_width=2)

###########################################
# STEP 4: ADD A FEW EXTRA ELEMENTS TO REALLY MAKE A NEAT PLOT


# First we want to actually add values of data to the plot for easy viewing
# to do this, lets first add a simple rectangle using matplotlib's 'patches'
# To do this, let's first add a simple rectangle using Matplotlib's 'patches'
# functionality to add some simple layout for plotting calculated parameters
# xloc yloc xsize ysize
fig.patches.extend([plt.Rectangle((0.513, 0.00), 0.334, 0.37,
fig.patches.extend([plt.Rectangle((0.563, 0.05), 0.334, 0.37,
edgecolor='black', facecolor='white',
linewidth=1, alpha=1, transform=fig.transFigure,
figure=fig)])

# now lets take a moment to calculate some simple severe-weather parameters using
# Now let's take a moment to calculate some simple severe-weather parameters using
# metpy's calculations
# here are some classic severe parameters!
# Here are some classic severe parameters!
kindex = mpcalc.k_index(p, T, Td)
total_totals = mpcalc.total_totals_index(p, T, Td)

Expand All @@ -240,7 +227,6 @@

# Compute Surface-based CAPE
sbcape, sbcin = mpcalc.surface_based_cape_cin(p, T, Td)

# Compute SRH
(u_storm, v_storm), *_ = mpcalc.bunkers_storm_motion(p, u, v, z)
*_, total_helicity1 = mpcalc.storm_relative_helicity(z, u, v, depth=1 * units.km,
Expand All @@ -265,89 +251,85 @@
# Perform the calculation of supercell composite if an effective layer exists
super_comp = mpcalc.supercell_composite(mucape, total_helicity3, bshear3)

# there is a lot we can do with this data operationally, so lets plot some of
# There is a lot we can do with this data operationally, so let's plot some of
# these values right on the plot, in the box we made
# first lets plot some thermodynamic parameters
plt.figtext(0.53, 0.32, 'SBCAPE: ', weight='bold', fontsize=15,
# First lets plot some thermodynamic parameters
plt.figtext(0.58, 0.37, 'SBCAPE: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.32, f'{int(sbcape.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.37, f'{int(sbcape.m)} J/kg', weight='bold',
fontsize=15, color='orangered', ha='right')
plt.figtext(0.53, 0.29, 'SBCIN: ', weight='bold',
plt.figtext(0.58, 0.34, 'SBCIN: ', weight='bold',
fontsize=15, color='black', ha='left')
plt.figtext(0.66, 0.29, f'{int(sbcin.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.34, f'{int(sbcin.m)} J/kg', weight='bold',
fontsize=15, color='lightblue', ha='right')

plt.figtext(0.53, 0.24, 'MLCAPE: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.29, 'MLCAPE: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.24, f'{int(mlcape.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.29, f'{int(mlcape.m)} J/kg', weight='bold',
fontsize=15, color='orangered', ha='right')
plt.figtext(0.53, 0.21, 'MLCIN: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.26, 'MLCIN: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.21, f'{int(mlcin.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.26, f'{int(mlcin.m)} J/kg', weight='bold',
fontsize=15, color='lightblue', ha='right')

plt.figtext(0.53, 0.16, 'MUCAPE: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.21, 'MUCAPE: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.16, f'{int(mucape.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.21, f'{int(mucape.m)} J/kg', weight='bold',
fontsize=15, color='orangered', ha='right')
plt.figtext(0.53, 0.13, 'MUCIN: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.18, 'MUCIN: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.13, f'{int(mucin.m)} J/kg', weight='bold',
plt.figtext(0.71, 0.18, f'{int(mucin.m)} J/kg', weight='bold',
fontsize=15, color='lightblue', ha='right')

plt.figtext(0.53, 0.08, 'TT-INDEX: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.13, 'TT-INDEX: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.08, f'{int(total_totals.m)} Δ°C', weight='bold',
plt.figtext(0.71, 0.13, f'{int(total_totals.m)} Δ°C', weight='bold',
fontsize=15, color='orangered', ha='right')
plt.figtext(0.53, 0.05, 'K-INDEX: ', weight='bold', fontsize=15,
plt.figtext(0.58, 0.10, 'K-INDEX: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.66, 0.05, f'{int(kindex.m)} °C', weight='bold',
plt.figtext(0.71, 0.10, f'{int(kindex.m)} °C', weight='bold',
fontsize=15, color='orangered', ha='right')

# now some kinematic parameters
met_per_sec = (units.m * units.m) / (units.sec * units.sec)
plt.figtext(0.68, 0.32, '0-1km SRH: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.37, '0-1km SRH: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.32, f'{int(total_helicity1.m)* met_per_sec:~P}',
plt.figtext(0.88, 0.37, f'{int(total_helicity1.m)* met_per_sec:~P}',
weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext(0.68, 0.29, '0-1km SHEAR: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.34, '0-1km SHEAR: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.29, f'{int(bshear1.m)} kts', weight='bold',
plt.figtext(0.88, 0.34, f'{int(bshear1.m)} kts', weight='bold',
fontsize=15, color='blue', ha='right')

plt.figtext(0.68, 0.24, '0-3km SRH: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.29, '0-3km SRH: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.24, f'{int(total_helicity3.m)* met_per_sec:~P}',
plt.figtext(0.88, 0.29, f'{int(total_helicity3.m)* met_per_sec:~P}',
weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext(0.68, 0.21, '0-3km SHEAR: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.26, '0-3km SHEAR: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.21, f'{int(bshear3.m)} kts', weight='bold',
plt.figtext(0.88, 0.26, f'{int(bshear3.m)} kts', weight='bold',
fontsize=15, color='blue', ha='right')

plt.figtext(0.68, 0.16, '0-6km SRH: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.21, '0-6km SRH: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.16, f'{int(total_helicity6.m)* met_per_sec:~P}',
plt.figtext(0.88, 0.21, f'{int(total_helicity6.m)* met_per_sec:~P}',
weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext(0.68, 0.13, '0-6km SHEAR: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.18, '0-6km SHEAR: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.13, f'{int(bshear6.m)} kts', weight='bold',
plt.figtext(0.88, 0.18, f'{int(bshear6.m)} kts', weight='bold',
fontsize=15, color='blue', ha='right')

plt.figtext(0.68, 0.08, 'SIG TORNADO: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.13, 'SIG TORNADO: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.08, f'{int(sig_tor.m)}', weight='bold', fontsize=15,
plt.figtext(0.88, 0.13, f'{int(sig_tor[0].m)}', weight='bold', fontsize=15,
color='orangered', ha='right')
plt.figtext(0.68, 0.05, 'SUPERCELL COMP: ', weight='bold', fontsize=15,
plt.figtext(0.73, 0.10, 'SUPERCELL COMP: ', weight='bold', fontsize=15,
color='black', ha='left')
plt.figtext(0.83, 0.05, f'{int(super_comp.m)}', weight='bold', fontsize=15,
plt.figtext(0.88, 0.10, f'{int(super_comp[0].m)}', weight='bold', fontsize=15,
color='orangered', ha='right')

# add legends to the skew and hodo
# Add legends to the skew and hodo
skewleg = skew.ax.legend(loc='upper left')
hodoleg = h.ax.legend(loc='upper left')

# add a plot title
plt.figtext(0.40, 0.92, 'OUN | MAY 4TH 1999 - 00Z VERTICAL PROFILE',
# add a quick plot title, this could be automated by
# declaring a station and datetime variable when using
# realtime observation data from Siphon.
plt.figtext(0.45, 0.97, 'OUN | MAY 4TH 1999 - 00Z VERTICAL PROFILE',
weight='bold', fontsize=20, ha='center')

# Show the plot
Expand Down