Smoothing Contours

from datetime import datetime, UTC

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

# Specify our date/time of product desired
dt = datetime.now(UTC)

# Construct the URL for our THREDDS Data Server Catalog,
# and access our desired dataset within via NCSS
cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/'
                 'NCEP/NAM/CONUS_12km/latest.xml')
ncss = cat.datasets[0].subset()

# Create our NCSS query with desired specifications
query = ncss.query()
query.time(dt)
query.add_lonlat()
query.variables('Geopotential_height_isobaric',
                'u-component_of_wind_isobaric',
                'v-component_of_wind_isobaric')

# Obtain the data we've queried for
data = ncss.get_data(query)

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

# Get dimension names to pull appropriate variables
# dtime = ds.Geopotential_height_isobaric.dims[0]
# dlev = ds.Geopotential_height_isobaric.dims[1]
# dlat = ds.Geopotential_height_isobaric.dims[2]
# dlon = ds.Geopotential_height_isobaric.dims[3]

# Get lat and lon data, as well as time data and metadata
lats = ds.lat
lons = ds.lon

# Need 2D lat/lons for plotting, do so if necessary
if lats.ndim < 2:
    lons, lats = np.meshgrid(lons, lats)

# Determine the level of 500 hPa
lev_500 = 500 * units.hPa

# Create more useable times for output
times = ds.Geopotential_height_isobaric.metpy.time.squeeze()
vtimes = times.values.astype('datetime64[ms]').astype('O')

# Pull out the 500 hPa Heights
hght_500 = ds.Geopotential_height_isobaric.metpy.sel(
    vertical=lev_500).squeeze()
uwnd_500 = ds['u-component_of_wind_isobaric'].metpy.sel(
    vertical=lev_500).squeeze()
vwnd_500 = ds['v-component_of_wind_isobaric'].metpy.sel(
    vertical=lev_500).squeeze()

# Calculate the magnitude of the wind speed in kts
sped = mpcalc.wind_speed(uwnd_500, vwnd_500).metpy.convert_units('knots')

plotcrs = ccrs.LambertConformal(central_longitude=-100.0,
                                central_latitude=45.0)
datacrs = ccrs.PlateCarree(central_longitude=0.)

# Smooth the 500-hPa geopotential height field
# Be sure to only smooth the 2D field
Z_500 = mpcalc.smooth_gaussian(hght_500, 50)

# Start plot with new figure and axis
fig = plt.figure(figsize=(17., 11.))
ax = plt.subplot(1, 1, 1, projection=plotcrs)

# Add some titles to make the plot readable by someone else
plt.title('500-hPa Geo Heights (m; black), Smoothed 500-hPa Geo. '
          'Heights (m; red)', loc='left')
plt.title(f'VALID: {vtimes}', loc='right')

# Set GAREA and add map features
ax.set_extent([-125., -67., 22., 52.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='black', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)

# Set the CINT
clev500 = np.arange(5100, 6000, 60)

# Plot smoothed 500-hPa contours
cs2 = ax.contour(lons, lats, Z_500, clev500, colors='black',
                 linewidths=3, linestyles='solid', transform=datacrs)
c2 = plt.clabel(cs2, fontsize=12, colors='black', inline=1, inline_spacing=8,
                fmt='%i', rightside_up=True, use_clabeltext=True)

# Contour the 500 hPa heights with labels
cs = ax.contour(lons, lats, hght_500, clev500, colors='red',
                linewidths=2.5, linestyles='solid', alpha=0.6,
                transform=datacrs)
cl = plt.clabel(cs, fontsize=12, colors='red', inline=1, inline_spacing=8,
                fmt='%i', rightside_up=True, use_clabeltext=True)