Plotting 500‑hPa Geopotential Height, Wind Speed, and Wind Barbs with Python, xarray, and Cartopy

This example demonstrates how to read NAM model output from a NetCDF file dated 31 Oct 2016 12:00 UTC using xarray , and visualize the 500‑hPa geopotential height, wind speed, and wind barbs with Cartopy and Matplotlib .

After opening the dataset, latitude and longitude arrays are extracted, and the 500‑hPa fields for geopotential height, u‑ and v‑wind components are selected and smoothed with a Gaussian filter.

MetPy is used to compute wind speed in knots, and a datetime object is built from the dataset time coordinate for labeling.

A Lambert Conformal projection is defined for the map, and the data are plotted: a color‑filled contour of wind speed, contour lines of geopotential height, and wind barbs (regridded) are added; titles show variable names and the valid time.

The figure layout is adjusted and the plot is saved to x.png . If the shapely library causes installation errors, reinstall it with pip uninstall shapely followed by pip install shapely --no-binary shapely .

<code># https://unidata.github.io/python-training/gallery/500hpa_hght_winds/
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter
import xarray as xr
ds = xr.open_dataset('/Users/Downloads/NAM_20161031_1200.nc')
lats = ds.lat.data
lons = ds.lon.data
hght = ds['Geopotential_height_isobaric']
uwnd = ds['u-component_of_wind_isobaric']
vwnd = ds['v-component_of_wind_isobaric']
hght_500 = gaussian_filter(hght.sel(isobaric=500).data[0], sigma=3.0)
uwnd_500 = gaussian_filter(uwnd.sel(isobaric=500).data[0], sigma=3.0) * units('m/s')
vwnd_500 = gaussian_filter(vwnd.sel(isobaric=500).data[0], sigma=3.0) * units('m/s')
sped_500 = mpcalc.wind_speed(uwnd_500, vwnd_500).to('kt')
vtime = datetime.strptime(str(ds.time.data[0].astype('datetime64[ms]')), '%Y-%m-%dT%H:%M:%S.%f')
mapcrs = ccrs.LambertConformal(central_longitude=-100, central_latitude=35, standard_parallels=(30, 60))
datacrs = ccrs.PlateCarree()
fig = plt.figure(1, figsize=(14, 12))
ax = plt.subplot(111, projection=mapcrs)
ax.set_extent([-130, -72, 20, 55], ccrs.PlateCarree())
clevs_500_sped = np.arange(30, 150, 20)
cf = ax.contourf(lons, lats, sped_500, clevs_500_sped, cmap=plt.cm.BuPu, transform=datacrs)
plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50)
clevs_500_hght = np.arange(0, 8000, 60)
cs = ax.contour(lons, lats, hght_500, clevs_500_hght, colors='black', transform=datacrs)
plt.clabel(cs, fmt='%d')
ax.barbs(lons, lats, uwnd_500.to('kt').m, vwnd_500.to('kt').m, pivot='middle', color='black', regrid_shape=20, transform=datacrs)
plt.title('500-hPa NAM Geopotential Heights (m), Wind Speed (kt), and Wind Barbs (kt)', loc='left')
plt.title('Valid Time: {}'.format(vtime), loc='right')
plt.subplots_adjust(bottom=0, top=1)
plt.savefig("x.png", bbox_inches='tight')
</code>