A 250-hPa Hemispheric Map using Python

from datetime import datetime, UTC

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

# Latest GFS Dataset
cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/'
                 'NCEP/GFS/Global_0p5deg/latest.xml')
best_ds = list(cat.datasets.values())[0]
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])

# Set the time to current
now = datetime.now(UTC)

# Query for Latest GFS Run
gfsdata_hght = ncss.query().time(now).accept('netcdf4')
gfsdata_hght.variables('Geopotential_height_isobaric')

# Set the lat/lon box for the data you want to pull in.
# lonlat_box(north_lat,south_lat,east_lon,west_lon)
gfsdata_hght.lonlat_box(0, 359.9, 0, 90)

# Set desired level 50000 = 50000 Pa = 500 hPa
gfsdata_hght.vertical_level(25000)

# Actually getting the data
data_hght = ncss.get_data(gfsdata_hght)

# Make into an xarray Dataset object
ds_hght = xr.open_dataset(NetCDF4DataStore(data_hght))

# Query for Latest GFS Run
gfsdata_wind = ncss.query().time(now).accept('netcdf4')
gfsdata_wind.variables('u-component_of_wind_isobaric',
                       'v-component_of_wind_isobaric')

# Set the lat/lon box for the data you want to pull in.
# lonlat_box(north_lat,south_lat,east_lon,west_lon)
gfsdata_wind.lonlat_box(0, 359.9, 0, 90)

# Set desired level 50000 = 50000 Pa = 500 hPa
gfsdata_wind.vertical_level(25000)

# Actually getting the data
data_wind = ncss.get_data(gfsdata_wind)

# Make into an xarray Dataset object
ds_wind = xr.open_dataset(NetCDF4DataStore(data_wind))

lat = ds_hght.latitude.values
lon = ds_hght.longitude.values

# Converting times using the num2date function available through netCDF4
vtimes = ds_hght.Geopotential_height_isobaric.metpy.time.values.astype('datetime64[ms]').astype('O')

# Smooth the 250-hPa heights using a gaussian filter from scipy.ndimage
hgt = ds_hght.Geopotential_height_isobaric.squeeze()
hgt_250, lon = cutil.add_cyclic_point(hgt, coord=lon)

Z_250 = mpcalc.smooth_n_point(hgt_250, 9, 2)

# Calculate wind speed from u and v components, add cyclic point,
# and smooth slightly
u250 = ds_wind['u-component_of_wind_isobaric'].squeeze()
v250 = ds_wind['v-component_of_wind_isobaric'].squeeze()

wspd250 = mpcalc.wind_speed(u250, v250).metpy.convert_units('knots')
wspd250 = cutil.add_cyclic_point(wspd250)

smooth_wspd250 = mpcalc.smooth_n_point(wspd250, 9, 2)

datacrs = ccrs.PlateCarree()
plotcrs = ccrs.NorthPolarStereo(central_longitude=-100.0)

# Make a grid of lat/lon values to use for plotting.
lons, lats = np.meshgrid(lon, lat)

fig = plt.figure(1, figsize=(12., 15.))
ax = plt.subplot(111, projection=plotcrs)

# Set some titles for the plots
ax.set_title('250-hPa Geopotential Heights (m)', loc='left')
ax.set_title(f'VALID: {vtimes[0]}', loc='right')

# Set the extent of the image for the NH and add
#   ax.set_extent([west long, east long, south lat, north lat])
ax.set_extent([-180, 180, 10, 90], ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='black',
               linewidth=0.5)
ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=0.5)

# Add geopotential height contours every 120 m
clev250 = np.arange(9000, 12000, 120)
cs = ax.contour(lons, lats, Z_250, clev250, colors='k',
                linewidths=1.0, linestyles='solid', transform=datacrs)
plt.clabel(cs, fontsize=8, inline=1, inline_spacing=10, fmt='%i',
           rightside_up=True, use_clabeltext=True)

# Add colorfilled wind speed in knots every 20 kts
clevsped250 = np.arange(50, 200, 20)
cf = ax.contourf(lons, lats, smooth_wspd250, clevsped250, cmap=plt.cm.BuPu,
                 extend='max', transform=datacrs)
cbar = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50,
                    extendrect=True)