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

Improve SkewT layout examples #2460

Closed
dopplershift opened this issue Apr 27, 2022 · 11 comments · Fixed by #3161
Closed

Improve SkewT layout examples #2460

dopplershift opened this issue Apr 27, 2022 · 11 comments · Fixed by #3161
Labels
Area: Examples Affects examples good first issue Straightforward issues suitable for new and inexperienced contributors to the project Type: Bug Something is not working like it should Type: Enhancement Enhancement to existing functionality
Milestone

Comments

@dopplershift
Copy link
Member

  • Need to fix up the Gridspec example since it currently looks like crap:

misaligned plots

That can be fixed by using fig = plt.figure(figsize=(12, 9)) (makes more room for the right sized skewT)

  • Also, add some version of an example that sets the adjustable to 'datalim' and also uses some manual boxes with rect. This example does that:
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import pandas as pd

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, Hodograph, SkewT
from metpy.units import units

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)

p = df['pressure'].values * units.hPa
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)

# 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')

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
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)

# Change to adjust data limits and give it a semblance of what we want
skew.ax.set_adjustable('datalim')
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-20, 30)

# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()


# Create a hodograph
ax = plt.axes((0.7, 0.75, 0.2, 0.2))
h = Hodograph(ax, component_range=60.)
h.add_grid(increment=20)
h.plot(u, v)

gives:
image

but it would be nice to come up with something different, maybe multiple skew-T's or adding some other panels. Potential inspiration:

image

image

@dopplershift dopplershift added Type: Bug Something is not working like it should Area: Examples Affects examples Type: Enhancement Enhancement to existing functionality labels Apr 27, 2022
@dopplershift dopplershift added this to the May 2022 milestone Apr 27, 2022
@dcamron dcamron added the good first issue Straightforward issues suitable for new and inexperienced contributors to the project label Apr 27, 2022
@eliteuser26
Copy link
Contributor

It would be nice to view the Skewt temperature vs dewpoint with parcel profile on one side with a table of Metpy severe calculated indices on the other side in the same graphic like we see on other websites. Also have the ability to calculate severe indices based on current temperature and dewpoint not just morning values and use it to display current parcel profile. This is something that I do right now on my website. Just an idea.

@dopplershift dopplershift modified the milestones: May 2022, July 2022 Jun 1, 2022
@dopplershift dopplershift modified the milestones: 1.4.1, April 2023 Mar 16, 2023
@kylejgillett
Copy link
Contributor

Hey @dopplershift!

I'd like to contribute to this issue. I have several versions of sounding plots that have a more complex layout. Attached is an example. Is a super-simple version of something like this what you are looking for?

07122023-02z_BattleCreekMi_HRRR_

@dopplershift
Copy link
Member Author

@kylejgillett That would be great! I think the trick will be figuring out the balance of putting in a lot of the cool elements without making the example too complex to serve as an "example". But yes, that's definitely along the lines of what I'm looking for.

@kylejgillett
Copy link
Contributor

@dopplershift, here is a mock idea of a more advanced/complex sounding viewer. It's based on the original example notebook and uses a few matplotlib, numpy, and metpy functions to increase readability and make it look a little cooler.

What do you think of this? Too much, too little? Any ideas?

metpy_advSkewT_test1

@dopplershift
Copy link
Member Author

I think that looks nice. The real decider on complexity, though, will be what code it takes to get that (though I'm guessing it's not too bad)

@kylejgillett
Copy link
Contributor

Essentially I took the original "Skew-T with Complex Layout" notebook and added this to it:

This layout isn't bad, but we could add a few simple tricks to greatly increase the readability and complexity of our Skew-T/Hodograph layout

Lets try another Skew-T:
 
'''
    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')

# Change to adjust data limits and give it a semblance of what we want
skew.ax.set_adjustable('datalim')
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-20, 30)

# Set some better labels than the default to increase readability
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
fig.set_facecolor('#ffffff')         
skew.ax.set_facecolor('#ffffff')     

# Here we can use some basic math and python functionality to make a cool
# shaded isotherm pattern. 
x1 = np.linspace(-100, 40, 8)                                                          
x2 = np.linspace(-90, 50, 8)                                                         
y = [1100, 50]                                                                      
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. 
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'
# the wind barbs for a cleaner output with increased readability. 
# Something like this would work.
interval = np.logspace(2, 3, 40) * units.hPa
idx = mpcalc.resample_nn_1d(p, interval) 
skew.plot_barbs(pressure=p[idx], u=u[idx], v=v[idx])

# 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
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
# ~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
# 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')

# Shade areas of CAPE and CIN
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))
h = Hodograph(hodo_ax, component_range=80.)
# Add two seperate 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
h.ax.set_box_aspect(1) 
h.ax.set_yticklabels([])
h.ax.set_xticklabels([])
h.ax.set_xticks([])
h.ax.set_yticks([])
h.ax.set_xlabel(' ')
h.ax.set_ylabel(' ')

# 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))
for i in range(10,120,10):
    h.ax.annotate(str(i),(i,0),xytext=(0,2),textcoords='offset pixels',clip_on=True,fontsize=10,weight='bold',alpha=0.3,zorder=0)
for i in range(10,120,10):
    h.ax.annotate(str(i),(0,i),xytext=(0,2),textcoords='offset pixels',clip_on=True,fontsize=10,weight='bold',alpha=0.3,zorder=0)

# plot the hodograph itself, using plot_colormapped, colored
# by height 
h.plot_colormapped(u, v, c=z, linewidth=6, label="0-12km WIND")
# compute Bunkers storm motion so we can plot it on the hodograph! 
RM, LM, MW = mpcalc.bunkers_storm_motion(p, u, v, z)
h.ax.text((RM[0].m +0.5), (RM[1].m -0.5), 'RM', weight='bold', ha='left', fontsize=13, alpha=0.6) 
h.ax.text((LM[0].m +0.5), (LM[1].m -0.5), 'LM', weight='bold', ha='left', fontsize=13, alpha=0.6) 
h.ax.text((MW[0].m +0.5), (MW[1].m -0.5), 'MW', weight='bold', ha='left', fontsize=13, alpha=0.6) 
h.ax.arrow(0,0,RM[0].m-0.3, RM[1].m-0.3, linewidth=2, color='black', 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' 
# fucntionality 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,
                                  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
# metpy's calculations 
# here are some classic severe parameters!
kindex = mpcalc.k_index(p, T, Td)
total_totals = mpcalc.total_totals_index(p, T, Td)

# mixed layer parcel properties!
ml_t, ml_td = mpcalc.mixed_layer(p, T, Td, depth=50 * units.hPa)
ml_p, _, _ = mpcalc.mixed_parcel(p, T, Td, depth=50 * units.hPa)
mlcape, mlcin = mpcalc.mixed_layer_cape_cin(p, T, prof, depth=50 * units.hPa)

# most unstable parcel properties!
mu_p, mu_t, mu_td, _ = mpcalc.most_unstable_parcel(p, T, Td, depth=50 * units.hPa)
mucape, mucin = mpcalc.most_unstable_cape_cin(p, T, Td, depth=50 * units.hPa)

# Estimate height of LCL in meters from hydrostatic thickness (for sig_tor)
new_p = np.append(p[p > lcl_pressure], lcl_pressure)
new_t = np.append(T[p > lcl_pressure], lcl_temperature)
lcl_height = mpcalc.thickness_hydrostatic(new_p, new_t)

# 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,
                                                    storm_u=u_storm, storm_v=v_storm)
*_, total_helicity3 = mpcalc.storm_relative_helicity(z, u, v, depth=3 * units.km,
                                                    storm_u=u_storm, storm_v=v_storm)
*_, total_helicity6 = mpcalc.storm_relative_helicity(z, u, v, depth=6 * units.km,
                                                    storm_u=u_storm, storm_v=v_storm)

# Copmute Bulk Shear components and then magnitude
ubshr1, vbshr1 = mpcalc.bulk_shear(p, u, v, height=z, depth=1 * units.km)
bshear1 = mpcalc.wind_speed(ubshr1, vbshr1)
ubshr3, vbshr3 = mpcalc.bulk_shear(p, u, v, height=z, depth=3 * units.km)
bshear3 = mpcalc.wind_speed(ubshr3, vbshr3)
ubshr6, vbshr6 = mpcalc.bulk_shear(p, u, v, height=z, depth=6 * units.km)
bshear6 = mpcalc.wind_speed(ubshr6, vbshr6)

# Use all computed pieces to calculate the Significant Tornado parameter
sig_tor = mpcalc.significant_tornado(sbcape, lcl_height,
                                     total_helicity3, bshear3).to_base_units()

# 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
# these values right on the plot, in the box we made
# first lets plot some thermodynamic parameters
plt.figtext( 0.53, 0.32,  f'SBCAPE: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.66, 0.32,  f'{int(sbcape.m)} J/kg', weight='bold', fontsize=15, color='orangered', ha='right')
plt.figtext( 0.53, 0.29,  f'SBCIN: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.66, 0.29,  f'{int(sbcin.m)} J/kg', weight='bold', fontsize=15, color='lightblue', ha='right')

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

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

plt.figtext( 0.53, 0.08,  f'TT-INDEX: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.66, 0.08,  f'{int(total_totals.m)} Δ°C', weight='bold', fontsize=15, color='orangered', ha='right')
plt.figtext( 0.53, 0.05,  f'K-INDEX: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.66, 0.05,  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,  f'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}', weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext( 0.68, 0.29,  f'0-1km SHEAR: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.83, 0.29,  f'{int(bshear1.m)} kts', weight='bold', fontsize=15, color='blue', ha='right')

plt.figtext( 0.68, 0.24,  f'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}', weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext( 0.68, 0.21,  f'0-3km SHEAR: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.83, 0.21,  f'{int(bshear3.m)} kts', weight='bold', fontsize=15, color='blue', ha='right')

plt.figtext( 0.68, 0.16,  f'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}', weight='bold', fontsize=15, color='navy', ha='right')
plt.figtext( 0.68, 0.13,  f'0-6km SHEAR: ', weight='bold', fontsize=15, color='black', ha='left')
plt.figtext( 0.83, 0.13,  f'{int(bshear6.m)} kts', weight='bold', fontsize=15, color='blue', ha='right')

plt.figtext( 0.68, 0.08,  f'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, color='orangered', ha='right')
plt.figtext( 0.68, 0.05,  f'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, color='orangered', ha='right')

# 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,  f'OUN | MAY 4TH 1999 - 00Z VERTICAL PROFILE', weight='bold', fontsize=20, ha='center')
# Show the plot
plt.show()

@dopplershift
Copy link
Member Author

I think that complexity is fine. I would replace the triple-quote strings with e.g.:

###########################################
# Step 1: Create the Skew-T object and modify it to create a nice, clean plot

Those blocks get converted to text in the notebook/html blocks in the rendered examples.

@kylejgillett
Copy link
Contributor

@dopplershift, ok, good idea!

Should I open a pull request for this issue for the .py: https://github.com/Unidata/MetPy/blob/main/examples/plots/Skew-T_Layout.py?

@dopplershift
Copy link
Member Author

Since your example relies on some fixed boxes, I think it would be good to leave "Skew-T Layout" examples as is using GridSpec and add yours as a new example.

@kylejgillett
Copy link
Contributor

Okay, nice. Is there anything I need to do to officially contribute this as a new example?

@dopplershift
Copy link
Member Author

If you add it in examples/ alongside the existing Advanced_Sounding.py example, submit a Pull Request with the new file, that should be enough. The rest of the machinery should take care of the rest.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Examples Affects examples good first issue Straightforward issues suitable for new and inexperienced contributors to the project Type: Bug Something is not working like it should Type: Enhancement Enhancement to existing functionality
Projects
Status: Done
Development

Successfully merging a pull request may close this issue.

4 participants