Automagic rotation of wind barbs to north-relative grid

[sgdecker]

Feature Request: I am trying to replicate GEMPAK's gdprof to a reasonable extent.

My GEMPAK settings:

 GPOINT   = yyr
 GDATTIM  = f0
 GVCORD   = PRES
 GFUNC    = TMPC
 GVECT    = WND
 GDFILE   = $MODEL/nam/19060312_nam212.gem
 LINE     = 3
 MARKER   = 0
 BORDER   = 10
 PTYPE    = skewt
 SCALE    = 0
 XAXIS    = -20/20/5
 YAXIS    = 1000/200/100/1;1;1
 WIND     = bk1
 REFVEC   =  
 WINPOS   = 1
 FILTER   = YES
 TITLE    = 1
 PANEL    = 0
 CLEAR    = YES
 TEXT     = 1
 DEVICE   = gf|yyr_sound.gif|800;1000
 OUTPUT   = T
 THTALN   = 0
 THTELN   = 0
 MIXRLN   = 0

The plot:

My MetPy attempt:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

from siphon.catalog import TDSCatalog
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units

cat_url = ('https://thredds-test.unidata.ucar.edu/thredds/catalog/grib/'
           'NCEP/NAM/CONUS_40km/conduit/'
           'NAM_CONUS_40km_conduit_20190603_1200.grib2/catalog.xml')
cat = TDSCatalog(cat_url)

ds = cat.datasets[-1]
nam = ds.remote_access(use_xarray=True).metpy.parse_cf()

# Grab data for F0
T = nam['Temperature_isobaric'][0,:]
u = nam['u-component_of_wind_isobaric'][0,:]
v = nam['v-component_of_wind_isobaric'][0,:]
p, = T.metpy.coordinates('vertical')

nam_crs = T.metpy.cartopy_crs

# Desired station: Goose Bay (53.319 N, 60.426 W)
stat_lat = 53.319
stat_lon = -60.426

stat_x, stat_y = nam_crs.transform_point(stat_lon, stat_lat,
                                         ccrs.PlateCarree())

# Compute the profile at station
T_stat = T.interp(x=stat_x, y=stat_y)
u_stat = u.interp(x=stat_x, y=stat_y)
v_stat = v.interp(x=stat_x, y=stat_y)

# Clean up the data so it plots nicely
T_stat = (T_stat.values * units('degK')).to('degC')
u_stat = (u_stat.values * units('m/s')).to('kt')
v_stat = (v_stat.values * units('m/s')).to('kt')
p = (p.values * units('Pa')).to('hPa')

# Make a basic Skew-T
fig = plt.figure(figsize=(8,10))
skew = SkewT(fig, rotation=45)
skew.plot(p, T_stat, 'r')
skew.plot_barbs(p[::2], u_stat[::2], v_stat[::2])
skew.ax.set_ylim(1000, 200)
skew.ax.set_xlim(-20, 20)
plt.show()

The plot:

For this issue, I am focusing on the wind profile, but as side issues, is my code as clean as currently possible, and is there a way not to have to figure out manually the lat and lon for Goose Bay?

Looking at the wind barbs, you'll see they differ. GEMPAK automagically plots them in north-relative coordinates. The MetPy sounding is in grid-relative coordinates, which is misleading in a sounding (or cross section for that matter). I tried to find functionality to do the conversion, but couldn't find any, and I do see VECR and VECN are blank in the GEMPAK Conversion Guide. I think ideally this would make use of the crs as well as a table of station IDs to be automatic.

• Python 3.7.3
• MetPy 0.10.0

(Edited so my GPOINT shows up)

[kgoebber]

@sgdecker Yes, the NAM data is grid relative, which is not easy information to find! I recently added an example of 500 hPa model output to the Python Gallery (http://unidata.github.io/python-gallery/examples/500hPa_Absolute_Vorticity_winds.html#sphx-glr-examples-500hpa-absolute-vorticity-winds-py) that addresses this and has an example function for converting to Earth-relative u and v winds. I haven't created a MetPy PR for this yet, but may, so any feedback would be great.

[kgoebber]

def earth_relative_wind_components(da, ugrd, vgrd):
    """Calculate the north-relative components of the wind from the grid-relative
    components using Cartopy transform_vectors.

    Parameters
    ----------
        ds : Xarray DataArray
            Either of u- or v-component of the wind
        ugrd : (M, N) Quantity
            grid relative u-component of the wind
        vgrd : (M, N) Quantity
            grid relative v-component of the wind

    Returns
    -------
        unr, vnr : tuple of array-like Quantity
            The north-relative wind components in the X (East-West) and Y (North-South)
            directions, respectively.
    """
    if 'crs' not in list(da.coords):
        raise ValueError('No CRS in coordinate, be sure to use the MetPy accessor parse_cf()')

    data_crs = da.metpy.cartopy_crs

    x = da.x.values
    y = da.y.values

    xx, yy = np.meshgrid(x, y)

    ut, vt = ccrs.PlateCarree().transform_vectors(data_crs, xx, yy, ugrd.m, vgrd.m)
    uer = ut * ugrd.units
    ver = vt * vgrd.units

    return uer, ver

[dopplershift]

What's the purpose of the da argument in the function above? Otherwise, it seems like a nice example of what we want future API calls to look like, at least as input. The only other question would be: Should we go ahead and to the work to return a DataArray?

[kgoebber]

Upon closer inspection, I guess it doesn't need the DataArray since all of the requisite information needed would all be contained within either the ugrd or vgrd object - since they are DataArray's themselves. It would likely work well to return as a DataArray as it would be easy to fit back into a Dataset if needed.

[dopplershift]

That was my thinking as well.