Plot Vorticity at Pole from GFS Data

gfs_vort.py

from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from metpy.units import units
import metpy.calc as mpcalc

# Grab a recent GFS run on Grid 3
base_url = 'https://www.ncei.noaa.gov/thredds/dodsC/gfs-003-files/'
dt = datetime(2018, 7, 23, 0)
data = xr.open_dataset('{}{dt:%Y%m}/{dt:%Y%m%d}/gfs_3_{dt:%Y%m%d}_'
                       '{dt:%H}00_000.grb2'.format(base_url, dt=dt),
                       decode_times=True)

# Compute grid spacing in lat/lon space
x = data.lon.values.astype('float64')
y = data.lat.values.astype('float64')

lons = np.deg2rad(x)
lats = np.deg2rad(y)

# Compute lat/lon deltas
dlon = np.diff(lons)
dlon2d = np.ones((lats.size, dlon.size)) * dlon

dlat = np.diff(lats)
dlat2d = np.ones((dlat.size, lons.size)) * np.expand_dims(dlat, axis=1)

# Convert deltas to physical distances
earth_radius = data.LatLon_Projection.earth_radius * units.m

dx = earth_radius * dlon2d * np.expand_dims(np.cos(lats), axis=1)
dy = earth_radius * dlat2d

# Get the valid times from the file
vtimes = []
for t in range(data.time.size):
    vtimes.append(datetime.utcfromtimestamp(data.time[t].data.astype('O') / 1e9))

idx_500 = 18  # Selecting isobaric=500 returns the data from index 3; bug in xarray?
u = data['u-component_of_wind_isobaric'][0,idx_500,:]
v = data['v-component_of_wind_isobaric'][0,idx_500,:]

vort = mpcalc.vorticity(u, v, dx, dy, dim_order='yx')

# Set up data and plot projections
datacrs = ccrs.PlateCarree()
plotcrs = ccrs.Stereographic(central_latitude=90., central_longitude=-90.)

# Plot map
fig = plt.figure(figsize=(17., 11.))
ax = plt.axes(projection=plotcrs)
ax.coastlines('50m', edgecolor='black')
ax.add_feature(cfeature.STATES, linewidth=0.5)
ax.set_extent([-179, 180, 60, 90], ccrs.PlateCarree())

# Add vorticity contours
clev500 = np.arange(-40,40,4)
cs = ax.contour(x, y, vort*1e5, clev500, transform=datacrs)
tl = plt.clabel(cs, fontsize=12, colors='k', inline=1, inline_spacing=8,
               fmt='%i', rightside_up=True, use_clabeltext=True)

# Finish the plot
plt.title('MetPy Vorticity, GFS Grid 3', loc='left')
plt.tight_layout()
plt.show()