Plotting a sounding with indices and a hodograph

Out:

##############
   INDICES
##############
0-1 km SRH: 355 m2/s2
0-1 km Shear: 31 kts
0-3 km SRH: 584 m2/s2
0-6 km Shear: 55 kts
EBWD: 59 kts
Eff. SRH: 527 m2/s2
K-index: 37
MLCAPE: 4205 J/kg
MLCIN: -59 J/kg
MLEL: 13558 m AGL
MLLCL: 1002 m AGL
MLLFC: 2393 m AGL
MLLI: -11 C
MUCAPE: 5769 J/kg
MUCIN: 0 J/kg
MUEL: 13882 m AGL
MULCL: 512 m AGL
MULFC: 612 m AGL
MULI: -13 C
PWV: 1.53 inch
SBCAPE: 5769 J/kg
SBCIN: 0 J/kg
SBEL: 13882 m AGL
SBLCL: 512 m AGL
SBLFC: 612 m AGL
SBLI: -13 C
SCP: 60.9
SHIP: 4.9
STP(cin): 13.9
STP(fix): 13.7

#############
 SARS OUTPUT
#############
Supercell
-----------
NO QUALITY MATCHES
Total Loose Matches: 31
# of Loose Matches that met Criteria: 24
SVR Probability: 0.7741935483870968

Hail
-----------
NO QUALITY MATCHES
Total Loose Matches: 22.0
# of Loose Matches that met Criteria: 20.0
SVR Probability: 0.9090909090909091

('SHARPpy Skew-T image output at:', '20140616.1900_OAX.png')

import warnings # Silence the warnings from SHARPpy
warnings.filterwarnings("ignore")
import sharppy.plot.skew as skew
from matplotlib.ticker import ScalarFormatter, MultipleLocator
from matplotlib.collections import LineCollection
import matplotlib.transforms as transforms
import matplotlib.pyplot as plt
from datetime import datetime
import numpy as np
from matplotlib import gridspec
from sharppy.sharptab import winds, utils, params, thermo, interp, profile
from sharppy.io.spc_decoder import SPCDecoder

def decode(filename):

    dec = SPCDecoder(filename)

    if dec is None:
        raise IOError("Could not figure out the format of '%s'!" % filename)

    # Returns the set of profiles from the file that are from the "Profile" class.
    profs = dec.getProfiles()
    stn_id = dec.getStnId()

    for k in list(profs._profs.keys()):
        all_prof = profs._profs[k]
        dates = profs._dates
        for i in range(len(all_prof)):
            prof = all_prof[i]
            new_prof = profile.create_profile(pres=prof.pres, hght=prof.hght, tmpc=prof.tmpc, dwpc=prof.dwpc, wspd=prof.wspd, \
                                              wdir=prof.wdir, strictQC=False, profile='convective', date=dates[i])
            return new_prof, dates[i], stn_id

FILENAME = 'data/14061619.OAX'

prof, time, location = decode(FILENAME)
# Bounds of the pressure axis
pb_plot=1050
pt_plot=100
dp_plot=10
plevs_plot = np.arange(pb_plot,pt_plot-1,-dp_plot)
# Open up the text file with the data in columns (e.g. the sample OAX file distributed with SHARPpy)
title = time.strftime('%Y%m%d/%H%M') + ' ' + location + '   (Observed)'

# Set up the figure in matplotlib.
fig = plt.figure(figsize=(9, 8))
gs = gridspec.GridSpec(4,4, width_ratios=[1,5,1,1])
ax = plt.subplot(gs[0:3, 0:2], projection='skewx')
skew.draw_title(ax, title)
ax.grid(True)
plt.grid(True)

# Plot the background variables
presvals = np.arange(1000, 0, -10)

ax.semilogy(prof.tmpc[~prof.tmpc.mask], prof.pres[~prof.tmpc.mask], 'r', lw=2)
ax.semilogy(prof.dwpc[~prof.dwpc.mask], prof.pres[~prof.dwpc.mask], 'g', lw=2)
ax.semilogy(prof.vtmp[~prof.dwpc.mask], prof.pres[~prof.dwpc.mask], 'r--')
ax.semilogy(prof.wetbulb[~prof.dwpc.mask], prof.pres[~prof.dwpc.mask], 'c-')

# Plot the parcel trace, but this may fail.  If it does so, inform the user.
try:
    ax.semilogy(prof.mupcl.ttrace, prof.mupcl.ptrace, 'k--')
except:
    print("Couldn't plot parcel traces...")

# Highlight the 0 C and -20 C isotherms.
l = ax.axvline(0, color='b', ls='--')
l = ax.axvline(-20, color='b', ls='--')

# Disables the log-formatting that comes with semilogy
ax.yaxis.set_major_formatter(ScalarFormatter())
ax.set_yticks(np.linspace(100,1000,10))
ax.set_ylim(1050,100)

# Plot the hodograph data.
inset_axes = skew.draw_hodo_inset(ax, prof)
skew.plotHodo(inset_axes, prof.hght, prof.u, prof.v, color='r')
#inset_axes.text(srwind[0], srwind[1], 'RM', color='r', fontsize=8)
#inset_axes.text(srwind[2], srwind[3], 'LM', color='b', fontsize=8)

# Draw the wind barbs axis and everything that comes with it.
ax.xaxis.set_major_locator(MultipleLocator(10))
ax.set_xlim(-50,50)
ax2 = plt.subplot(gs[0:3,2])
ax3 = plt.subplot(gs[3,0:3])
skew.plot_wind_axes(ax2)
skew.plot_wind_barbs(ax2, prof.pres, prof.u, prof.v)
srwind = params.bunkers_storm_motion(prof)
gs.update(left=0.05, bottom=0.05, top=0.95, right=1, wspace=0.025)

# Calculate indices to be shown.  More indices can be calculated here using the tutorial and reading the params module.
p1km = interp.pres(prof, interp.to_msl(prof, 1000.))
p6km = interp.pres(prof, interp.to_msl(prof, 6000.))
sfc = prof.pres[prof.sfc]
sfc_1km_shear = winds.wind_shear(prof, pbot=sfc, ptop=p1km)
sfc_6km_shear = winds.wind_shear(prof, pbot=sfc, ptop=p6km)
srh3km = winds.helicity(prof, 0, 3000., stu = srwind[0], stv = srwind[1])
srh1km = winds.helicity(prof, 0, 1000., stu = srwind[0], stv = srwind[1])
scp = params.scp(prof.mupcl.bplus, prof.right_esrh[0], prof.ebwspd)
stp_cin = params.stp_cin(prof.mlpcl.bplus, prof.right_esrh[0], prof.ebwspd, prof.mlpcl.lclhght, prof.mlpcl.bminus)
stp_fixed = params.stp_fixed(prof.sfcpcl.bplus, prof.sfcpcl.lclhght, srh1km[0], utils.comp2vec(prof.sfc_6km_shear[0], prof.sfc_6km_shear[1])[1])
ship = params.ship(prof)

# A routine to perform the correct formatting when writing the indices out to the figure.
def fmt(value, fmt='int'):
    if fmt == 'int':
        try:
            val = int(value)
        except:
            val = str("M")
    else:
        try:
            val = round(value,1)
        except:
            val = "M"
    return val

# Setting a dictionary that is a collection of all of the indices we'll be showing on the figure.
# the dictionary includes the index name, the actual value, and the units.
indices = {'SBCAPE': [fmt(prof.sfcpcl.bplus), 'J/kg'],\
           'SBCIN': [fmt(prof.sfcpcl.bminus), 'J/kg'],\
           'SBLCL': [fmt(prof.sfcpcl.lclhght), 'm AGL'],\
           'SBLFC': [fmt(prof.sfcpcl.lfchght), 'm AGL'],\
           'SBEL': [fmt(prof.sfcpcl.elhght), 'm AGL'],\
           'SBLI': [fmt(prof.sfcpcl.li5), 'C'],\
           'MLCAPE': [fmt(prof.mlpcl.bplus), 'J/kg'],\
           'MLCIN': [fmt(prof.mlpcl.bminus), 'J/kg'],\
           'MLLCL': [fmt(prof.mlpcl.lclhght), 'm AGL'],\
           'MLLFC': [fmt(prof.mlpcl.lfchght), 'm AGL'],\
           'MLEL': [fmt(prof.mlpcl.elhght), 'm AGL'],\
           'MLLI': [fmt(prof.mlpcl.li5), 'C'],\
           'MUCAPE': [fmt(prof.mupcl.bplus), 'J/kg'],\
           'MUCIN': [fmt(prof.mupcl.bminus), 'J/kg'],\
           'MULCL': [fmt(prof.mupcl.lclhght), 'm AGL'],\
           'MULFC': [fmt(prof.mupcl.lfchght), 'm AGL'],\
           'MUEL': [fmt(prof.mupcl.elhght), 'm AGL'],\
           'MULI': [fmt(prof.mupcl.li5), 'C'],\
           '0-1 km SRH': [fmt(srh1km[0]), 'm2/s2'],\
           '0-1 km Shear': [fmt(utils.comp2vec(sfc_1km_shear[0], sfc_1km_shear[1])[1]), 'kts'],\
           '0-3 km SRH': [fmt(srh3km[0]), 'm2/s2'],\
           '0-6 km Shear': [fmt(utils.comp2vec(sfc_6km_shear[0], sfc_6km_shear[1])[1]), 'kts'],\
           'Eff. SRH': [fmt(prof.right_esrh[0]), 'm2/s2'],\
           'EBWD': [fmt(prof.ebwspd), 'kts'],\
           'PWV': [round(prof.pwat, 2), 'inch'],\
           'K-index': [fmt(params.k_index(prof)), ''],\
           'STP(fix)': [fmt(stp_fixed, 'flt'), ''],\
           'SHIP': [fmt(ship, 'flt'), ''],\
           'SCP': [fmt(scp, 'flt'), ''],\
           'STP(cin)': [fmt(stp_cin, 'flt'), '']}

# List the indices within the indices dictionary on the side of the plot.
trans = transforms.blended_transform_factory(ax.transAxes,ax.transData)

# Write out all of the indices to the figure.
print("##############")
print("   INDICES    ")
print("##############")
string = ''
keys = np.sort(list(indices.keys()))
x = 0
counter = 0
for key in keys:
    string = string + key + ': ' + str(indices[key][0]) + ' ' + indices[key][1] + '\n'
    print((key + ": " + str(indices[key][0]) + ' ' + indices[key][1]))
    if counter < 7:
        counter += 1
        continue
    else:
        counter = 0
        ax3.text(x, 1, string, verticalalignment='top', transform=ax3.transAxes, fontsize=11)
        string = ''
        x += 0.3
ax3.text(x, 1, string, verticalalignment='top', transform=ax3.transAxes, fontsize=11)
ax3.set_axis_off()

# Show SARS matches (edited for Keith Sherburn)
try:
    supercell_matches = prof.supercell_matches
    hail_matches = prof.matches
except:
    supercell_matches = prof.right_supercell_matches
    hail_matches = prof.right_matches

print()
print("#############")
print(" SARS OUTPUT ")
print("#############")
for mtype, matches in zip(['Supercell', 'Hail'], [supercell_matches, hail_matches]):
    print(mtype)
    print('-----------')
    if len(matches[0]) == 0:
        print("NO QUALITY MATCHES")
    for i in range(len(matches[0])):
        print(matches[0][i] + ' ' + matches[1][i])
    print("Total Loose Matches:", matches[2])
    print("# of Loose Matches that met Criteria:", matches[3])
    print("SVR Probability:", matches[4])
    print()

# Finalize the image formatting and alignments, and save the image to the file.
gs.tight_layout(fig)
fn = time.strftime('%Y%m%d.%H%M') + '_' + location + '.png'
fn = fn.replace('/', '')
print(("SHARPpy Skew-T image output at:", fn))
plt.savefig(fn, bbox_inches='tight', dpi=180)