Plot an 850 hPa map with calculating advection using MetPy.
Beyond just plotting 850-hPa level data, this uses calculations from metpy.calc to find
the temperature advection. Currently, this needs an extra helper function to calculate
the distance between lat/lon grid points.
Imports
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 xrCreate NCSS object to access the NetcdfSubset¶
Data from NCEI GFS 0.5 deg Analysis Archive
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/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.lat.values
lon = ds.lon.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)---------------------------------------------------------------------------
HTTPError Traceback (most recent call last)
Cell In[2], line 6
3 # Assemble our URL to the THREDDS Data Server catalog,
4 # and access our desired dataset within via NCSS
5 base_url = 'https://www.ncei.noaa.gov/thredds/model-gfs-g4-anl-files-old/'
----> 6 cat = TDSCatalog(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml')
7 ncss = cat.datasets[f'gfsanl_4_{dt:%Y%m%d}_{dt:%H}00_000.grb2'].subset()
9 # Create NCSS query for our desired time, region, and data variables
File ~/micromamba/envs/metpy-cookbook/lib/python3.13/site-packages/siphon/catalog.py:289, in TDSCatalog.__init__(self, catalog_url)
287 # get catalog.xml file
288 resp = self.session.get(catalog_url)
--> 289 resp.raise_for_status()
291 # top level server url
292 self.catalog_url = resp.url
File ~/micromamba/envs/metpy-cookbook/lib/python3.13/site-packages/requests/models.py:1026, in Response.raise_for_status(self)
1021 http_error_msg = (
1022 f"{self.status_code} Server Error: {reason} for url: {self.url}"
1023 )
1025 if http_error_msg:
-> 1026 raise HTTPError(http_error_msg, response=self)
HTTPError: 404 Client Error: for url: https://www.ncei.noaa.gov/thredds/model-gfs-g4-anl-files-old/201704/20170405/catalog.xmlBegin data calculations¶
# 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)Begin map creation¶
# 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);