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

Problem with units display on SkewT #1335

Closed Soudaaa closed 4 years ago

Soudaaa commented 4 years ago

So, I'm trying to adaptate a SkewT example script to plot UWYO soundings with some severe weather parameters included, which is running fine.

The thing that is bothering me, is the fact that the parameters units are not abbreviated when plotting the sounding, which consumes a lot of space and ends up getting out of the bounds of the image.

Although not a major issue, I tried to fix this with no success and I'm not sure if this somekind of bug or I'm not seeing what I'm doing wrong (probably the latter), so I've included the code and the output.

import posixpath
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import matplotlib.gridspec as gridspec
import numpy as np
from metpy.plots import add_metpy_logo, add_timestamp, SkewT, Hodograph
from metpy.units import units
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from siphon.simplewebservice.wyoming import WyomingUpperAir

def plot_skewt(df):
    # We will pull the data out of the example dataset into individual variables
    # and assign units.
    height_AGL = (df['height']-df['elevation']).values * units.m
    z = df['height'].values * units.m
    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.
    ax = plt.figure(figsize=(9,9))

    # Grid for plots
    gs = gridspec.GridSpec(3, 3)
    skew = SkewT(ax, rotation=30, subplot=gs[:,:2])

    # 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[::2], u[::2], v[::2], y_clip_radius=0.01)
    skew.ax.set_ylim(1000, 100)
    skew.ax.set_xlim(-40, 40)

    # Calculate LCL height and plot as black dot
    lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
    lfc_pressure, lfc_temperature = mpcalc.lfc(p, T, Td)
    lcl_hgt = np.round(mpcalc.pressure_to_height_std(lcl_pressure), decimals=3).to(units.meter)
        lfc_hgt = np.round(mpcalc.pressure_to_height_std(lfc_pressure), decimals=3).to(units.meter)

    skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
    skew.plot(lfc_pressure, lfc_temperature, 'ko', markerfacecolor='blue')

    sb_cape, sb_cin = mpcalc.surface_based_cape_cin(p, T, Td)
    ml_cape, ml_cin = mpcalc.mixed_layer_cape_cin(p, T, Td)
    mu_cape, mu_cin = mpcalc.most_unstable_cape_cin(p, T, Td)
    lr_700_500 = np.round(-1 * np.divide(T[20]-T[14], (z[20]-z[14]).to(units.kilometer)),2)
    sbcape = np.round(sb_cape, 1)
    sbcin = np.round(sb_cin, 1)
    mlcape = np.round(ml_cape, 1)
    mlcin = np.round(ml_cin, 1)
    mucape = np.round(mu_cape, 1)

    u_shear01, v_shear01 = mpcalc.bulk_shear(p, u.to(units.meter/units.second), v.to(units.meter/units.second), depth = 1000 * units.meter)
    shear01 = np.round((np.sqrt(u_shear01**2 + v_shear01**2)), 1)
    u_shear06, v_shear06 = mpcalc.bulk_shear(p, u.to(units.meter/units.second), v.to(units.meter/units.second), depth = 6000 * units.meter)
    shear06 = np.round((np.sqrt(u_shear06**2 + v_shear06**2)), 1)
    rmover, lmover, mean = mpcalc.bunkers_storm_motion(p, u.to(units.meter/units.second), v.to(units.meter/units.second), z)
    srh_01_pos, srh_01_neg, srh_01_tot = mpcalc.storm_relative_helicity(u.to(units.meter/units.second), v.to(units.meter/units.second), z, depth = 1000 * units.meter, bottom = height_AGL[0], storm_u = lmover[0], storm_v = lmover[1])
    srh_01 = np.round(srh_01_neg, 1)
    srh_03_pos, srh_03_neg, srh_03_tot = mpcalc.storm_relative_helicity(u.to(units.meter/units.second), v.to(units.meter/units.second), z, depth = 3000 * units.meter, bottom = height_AGL[0], storm_u = lmover[0], storm_v = lmover[1])
    srh_03 = np.round(srh_03_neg, 1)

    plt.figtext( 0.65, 0.60, 'LCL Height:')
    plt.figtext( 0.8, 0.60, '{}'.format(lcl_hgt))
    plt.figtext( 0.65, 0.58, 'LFC Height:')
    plt.figtext( 0.8, 0.58, '{}'.format(lfc_hgt))
    plt.figtext( 0.65, 0.56, 'MLLR:')
    plt.figtext( 0.8, 0.56, '{}'.format(lr_700_500))
    plt.figtext( 0.65, 0.54, 'SBCAPE:')
    plt.figtext( 0.8, 0.54, '{0} J/kg'.format(sbcape))
    plt.figtext( 0.65, 0.52, 'SBCIN:')
    plt.figtext( 0.8, 0.52, '{}'.format(sbcin))
    plt.figtext( 0.65, 0.50, 'MLCAPE:')
    plt.figtext( 0.8, 0.50, '{}'.format(mlcape))
    plt.figtext( 0.65, 0.48, 'MLCIN:')
    plt.figtext( 0.8, 0.48, '{}'.format(mlcin))
        plt.figtext( 0.65, 0.46, 'MUCAPE:')
    plt.figtext( 0.8, 0.46, '{}'.format(mucape))
    plt.figtext( 0.65, 0.44, 'Shear 0-1 km:')
    plt.figtext( 0.8, 0.44, '{}'.format(shear01))
    plt.figtext( 0.65, 0.42, 'Shear 0-6 km:')
    plt.figtext( 0.8, 0.42, '{}'.format(shear06))
    plt.figtext( 0.65, 0.40, 'SRH 0-1 km:')
    plt.figtext( 0.8, 0.40, '{}'.format(srh_01))
    plt.figtext( 0.65, 0.38, 'SRH 0-3 km:')
    plt.figtext( 0.8, 0.38, '{}'.format(srh_03))

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

    # An example of a slanted line at constant T -- in this case the 0
    # isotherm
    skew.ax.axvline(0, color='c', linestyle='--', linewidth=1)
    skew.ax.axvline(-20, color='c', linestyle='--', linewidth=1)

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

    # Create a hodograph
    ax1 = ax.add_subplot(gs[0,-1])
    h = Hodograph(ax1, component_range=80.)
    h.add_grid(increment=20)
    wind_speed = np.sqrt(u**2 + v**2)
    h.plot(u, v)

    return skew

def make_name(site, time):
    return '{site}_{dt:%Y%m%d_%H%M}.png'.format(site=site, dt=time)
if __name__ == '__main__':
    import argparse
    from datetime import datetime, timedelta
    import tempfile
        # Set up argument parsing for the script. Provides one argument for the site, and another
    # that controls whether the plot should be shown or saved as an image.
    parser = argparse.ArgumentParser(description='Download sounding data and plot.')
    parser.add_argument('-s', '--site', help='Site to obtain data for', type=str,
                        default='DNR')
    parser.add_argument('--show', help='Whether to show the plot rather than save to disk',
                        action='store_true')
    parser.add_argument('-d', '--date', help='Date and time to request data for in YYYYMMDDHH.'
                        ' Defaults to most recent 00/12 hour.', type=str)
    parser.add_argument('-g', '--gdrive', help='Google Drive upload path', type=str)
    parser.add_argument('-f', '--filename', help='Image filename', type=str)
    args = parser.parse_args()

    if args.date:
        request_time = datetime.strptime(args.date, '%Y%m%d%H')
    else:
        # Figure out the most recent sounding, 00 or 12. Subtracting two hours
        # helps ensure that we choose a time with data available.
        now = datetime.utcnow() - timedelta(hours=2)
        request_time = now.replace(hour=(now.hour // 12) * 12, minute=0, second=0)

    # Request the data and plot
    df = WyomingUpperAir.request_data(request_time, args.site)
    skewt = plot_skewt(df)

    # Add the timestamp for the data to the plot
    plt.title('Valid Time: {}'.format(request_time), loc='right', x=1.3, y=1)
    plt.title('{} Sounding'.format(args.site), fontsize=14, fontweight='bold', loc='left', x=-2.7, y=1)

    if args.show:
        plt.show()
    else:
        fname = args.filename if args.filename else make_name(args.site, request_time)
        if args.gdrive:
            uploader = DriveUploader()
            with tempfile.NamedTemporaryFile(suffix='.png') as f:
                skewt.ax.figure.savefig(f.name)
                uploader.upload_to(f.name, posixpath.join(args.gdrive, fname))
        else:
            skewt.ax.figure.savefig(make_name(args.site, request_time))

Figure_1

eliteuser26 commented 4 years ago

I just did the same research about converting to abbreviated units. As you might know, Pint module is part of the Metpy code as a dependency. This is the module that attach the units to the values and uses the long format by default. As such, I searched in the Pint documentation where it gives examples in converting from long format to abbreviated format for units.

There are 2 ways to do this with the print function depending on the Python version. You need to use a Pint print formatter to abbreviate the units. As an example:

Python 3.6+:
plt.figtext( 0.8, 0.54, f'{sbcape:~P}')

Older formatting method in Python:
plt.figtext( 0.8, 0.54, '{:~P}'.format(sbcape)) # removed J/kg in the string

I was debating if I was going to put the severity indexes inside or outside the Skew-T graphic. Because I display on a web page I decided to put it outside in a web table. I was looking for a way to abbreviate units which I did find. So it is not a bug after all.

N.B. Please be aware that this will convert to a string. It won't retain the Pint format.

Soudaaa commented 4 years ago

I just did the same research about converting to abbreviated units. As you might know, Pint module is part of the Metpy code as a dependency. This is the module that attach the units to the values and uses the long format by default. As such, I searched in the Pint documentation where it gives examples in converting from long format to abbreviated format for units.

There are 2 ways to do this with the print function depending on the Python version. You need to use a Pint print formatter to abbreviate the units. As an example:

Python 3.6+:
plt.figtext( 0.8, 0.54, f'{sbcape:~P}')

Older formatting method in Python:
plt.figtext( 0.8, 0.54, '{:~P}'.format(sbcape)) # removed J/kg in the string

I was debating if I was going to put the severity indexes inside or outside the Skew-T graphic. Because I display on a web page I decided to put it outside in a web table. I was looking for a way to abbreviate units which I did find. So it is not a bug after all.

N.B. Please be aware that this will convert to a string. It won't retain the Pint format.

Thanks for your answer.

Problem solved.