850 hPa Temperature Advection

from datetime import datetime

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import numpy as np
from metpy.units import units
from siphon.catalog import TDSCatalog
from xarray.backends import NetCDF4DataStore
import xarray as xr

dt = datetime(2017, 4, 5, 12)

# Assemble our URL to the THREDDS Data Server catalog,
# and access our desired dataset within via NCSS
base_url = 'https://www.ncei.noaa.gov/thredds/catalog/model-gfs-g4-anl-files-old/'
cat = TDSCatalog(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml')
ncss = cat.datasets[f'gfsanl_4_{dt:%Y%m%d}_{dt:%H}00_000.grb2'].subset()

# Create NCSS query for our desired time, region, and data variables
query = ncss.query()

query.time(dt)
query.lonlat_box(north=65, south=15, east=310, west=220)
query.accept('netcdf')
query.variables('Geopotential_height_isobaric',
                'Temperature_isobaric',
                'u-component_of_wind_isobaric',
                'v-component_of_wind_isobaric')

# Obtain the queried data
data = ncss.get_data(query)

# Make into an xarray Dataset object
ds = xr.open_dataset(NetCDF4DataStore(data)).metpy.parse_cf()

# Pull out variables you want to use
level = 850 * units.hPa
hght_850 = ds.Geopotential_height_isobaric.metpy.sel(
    vertical=level).squeeze()
temp_850 = ds.Temperature_isobaric.metpy.sel(
    vertical=level).squeeze()
u_wind_850 = ds['u-component_of_wind_isobaric'].metpy.sel(
    vertical=level).squeeze()
v_wind_850 = ds['v-component_of_wind_isobaric'].metpy.sel(
    vertical=level).squeeze()
time = hght_850.metpy.time
lat = ds.latitude.values
lon = ds.longitude.values

# Convert number of hours since the reference time into an actual date
vtime = time.values.astype('datetime64[ms]').astype('O')

# Combine 1D latitude and longitudes into a 2D grid of locations
lon_2d, lat_2d = np.meshgrid(lon, lat)

# Calculate temperature advection using metpy function
adv = mpcalc.advection(temp_850, u_wind_850, v_wind_850)

# Smooth heights and advection a little
# Be sure to only put in a 2D lat/lon or Y/X array for smoothing
Z_850 = mpcalc.smooth_gaussian(hght_850, 2)
adv = mpcalc.smooth_gaussian(adv, 2)

# Set Projection of Data
datacrs = ccrs.PlateCarree()

# Set Projection of Plot
plotcrs = ccrs.LambertConformal(central_latitude=45,
                                central_longitude=-100, standard_parallels=[30, 60])

# Create new figure
fig = plt.figure(figsize=(14, 12))

# Add the map and set the extent
ax = plt.subplot(111, projection=plotcrs)
plt.title(f'850mb Temperature Advection for {vtime:%d %B %Y %H:%MZ}',
          fontsize=16)
ax.set_extent([235., 290., 20., 55.])

# Add state/country boundaries to plot
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.BORDERS)

# Plot Height Contours
clev850 = np.arange(900, 3000, 30)
cs = ax.contour(lon_2d, lat_2d, Z_850, clev850, colors='black', linewidths=1.5,
                linestyles='solid', transform=datacrs)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
           rightside_up=True, use_clabeltext=True)

# Plot Temperature Contours
clevtemp850 = np.arange(-20, 20, 2)
cs2 = ax.contour(lon_2d, lat_2d, temp_850.metpy.convert_units('degC'),
                 clevtemp850, colors='grey', linewidths=1.25,
                 linestyles='dashed', transform=datacrs)
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
           rightside_up=True, use_clabeltext=True)

# Plot Colorfill of Temperature Advection
cint = np.arange(-8, 9)
cf = ax.contourf(lon_2d, lat_2d, 3*adv.metpy.convert_units('delta_degC/hour'),
                 cint[cint != 0],
                 extend='both', cmap='bwr', transform=datacrs)
cb = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50,
                  extendrect=True, ticks=cint)
cb.set_label(r'$^{o}C/3h$', size='large')

wind_slice = (slice(None, None, 10), slice(None, None, 10))
# Plot Wind Barbs
ax.barbs(lon_2d[wind_slice], lat_2d[wind_slice],
         u_wind_850.metpy.convert_units('kt').values[wind_slice],
         v_wind_850.metpy.convert_units('kt').values[wind_slice],
         length=6, pivot='middle', transform=datacrs);