Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.21k stars 408 forks source link

Improve SkewT layout examples #2460

Closed dopplershift closed 10 months ago

dopplershift commented 2 years ago

misaligned plots

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

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:
<img width="574" alt="image" src="https://user-images.githubusercontent.com/221526/165427282-8da0d00f-dd68-43aa-b6b0-c4d648693327.png">

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

![image](https://user-images.githubusercontent.com/221526/165427030-cc0de28d-eea7-4763-8615-a5edf1b62e23.png)

![image](https://user-images.githubusercontent.com/221526/165427524-38d4d4d2-b5ff-4f1f-bae7-9df2fa4aa448.png)
eliteuser26 commented 2 years ago

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.

kylejgillett commented 11 months ago

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 commented 11 months ago

@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 commented 11 months ago

@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 commented 11 months ago

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 commented 11 months ago

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 commented 11 months ago

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 commented 11 months ago

@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 commented 11 months ago

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 commented 11 months ago

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

dopplershift commented 10 months ago

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.