[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[python #XTF-300986]: Follow-up on Pint Quantity problem

Hi,

No problem. Thanks for sending it over. Working with this script 
was...illuminating. Not because there's anything fundamentally wrong there, it 
just helps to see how people tackle things, and more importantly be inspired to 
see what challenges exist. And there were *many* challenges as I tweaked your 
script to work like I wanted it to. So here's the version I arrived at:

    import matplotlib as mpl
    mpl.use('Agg')
    import matplotlib.pyplot as plt
    import numpy as np
    from datetime import datetime
    from siphon.catalog import TDSCatalog
    import metpy.calc as mpcalc
    from metpy.units import units
    from metpy.plots import SkewT

    # Use siphon to get a random weather model dataset
    dt = datetime(2020, 7, 19, 0)
    cat = 
TDSCatalog(f"https://www.ncei.noaa.gov/thredds/catalog/model-rap130/{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml";)
    ds = cat.datasets[f"rap_130_{dt:%Y%m%d_%H%M_000.grb2}"]
    data = ds.remote_access(use_xarray=True)
    data = data.metpy.parse_cf()

    # Get the data fields from the file
    tmpk = data['Temperature_isobaric'][0,:,150,150] # Give it any random y,x

    # Flip here and everything else will be consistent
    tmpk = np.flip(tmpk)

    # Pressure needs to go bottom to top--vertical always points to the 
vertical coordinate
    pres = tmpk.metpy.vertical

    # Converts XArray to Pint-wrapping numpy using units metadata
    # This *shouldn't* be needed but makes plotting work far more easily
    tmpk = tmpk.metpy.unit_array
    pres = pres.metpy.unit_array

    # Initialize the plot and overlay the data
    fig = plt.figure(figsize=(1600/96, 1600/96))
    skew = SkewT(fig, rotation=45)
    skew.plot(pres, tmpk, 'r', linewidth=3.) # This works fine

    # Calculate a parcel profile
    # Just use some dummy surface temps and dew points for the sake of this 
example
    tmpcsfc = 40. * units.degC
    tdsfc = 20. * units.degC
    parcel_prof = mpcalc.parcel_profile(pres, tmpcsfc, tdsfc)
    skew.shade_cape(pres, tmpk, parcel_prof, facecolor='#ffeb66', alpha=0.4)

    # Set the plot to use degC
    skew.ax.xaxis.set_units(units.degC)

Instead of using NCSS, I had it use OPeNDAP and xarray--which paves the way to 
avoiding a lot of hand-manipulation of units. Unfortunately, the happy path is 
not nearly as robust as I'd like.

If any of that is unclear, I'm happy to clarify. Thanks for some inspiration on 
how we can make things better!

Ryan


> Here's the example script I was talking about on Twitter. Thanks again for
> all you do!


Ticket Details
===================
Ticket ID: XTF-300986
Department: Support Python
Priority: Low
Status: Closed
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.