This is a tutorial demonstrates common analyses for Synoptic Meteorology courses with use of Unidata tools, specifically MetPy and Siphon. In this tutorial we will cover accessing, calculating, and plotting model output.
Let's investigate The Storm of the Century, although it would easy to change which case you wanted (please feel free to do so).
Reanalysis Output: NARR 00 UTC 13 March 1993
In [ ]:
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from scipy.ndimage import gaussian_filter
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import metpy.constants as mpconstants
from metpy.units import units
import xarray as xr
Lets investigate what specific NARR output is available to work with from NCEI.
We specifically want to look for data that has "TDS" data access, since that is short for a THREDDS server data access point. There are a total of four different GFS datasets that we could potentially use.
Choosing our data source
Let's go ahead and use the NARR Analysis data to investigate the past case we identified (The Storm of the Century).
And we will use a python package called Siphon to read this data through the NetCDFSubset (NetCDFServer) link.
First we can set out date using the datetime module
In [ ]:
# Case Study Date
year = 1993
month = 3
day = 13
hour = 0
dt = datetime(year, month, day, hour)
Next, we set up access to request subsets of data from the model. This uses the NetCDF Subset Service (NCSS) to make requests from the GRIB collection and get results in netCDF format.
In [ ]:
# Read NARR Data from THREDDS server
base_url = 'https://www.ncei.noaa.gov/thredds/catalog/narr-a-files/'
# Programmatically generate the URL to the day of data we want
cat = TDSCatalog(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml')
# Have Siphon find the appropriate dataset
ds = cat.datasets.filter_time_nearest(dt)
# Download data using the NetCDF Subset Service
ncss = ds.subset()
query = ncss.query().lonlat_box(north=60, south=18, east=300, west=225)
query.time(dt).variables('Geopotential_height_isobaric',
'Temperature_isobaric',
'u-component_of_wind_isobaric',
'v-component_of_wind_isobaric').add_lonlat().accept('netcdf')
data = ncss.get_data(query)
In [ ]:
# Open data with xarray, and parse it with MetPy
ds = xr.open_dataset(xr.backends.NetCDF4DataStore(data)).metpy.parse_cf()
ds
In [ ]:
# Back up in case of bad internet connection.
# Uncomment the following line to read local netCDF file of NARR data
# ds = xr.open_dataset('../../data/NARR_19930313_0000.nc').metpy.parse_cf()
In [ ]:
# Save time of data to datetime format
time_format = '%Y-%m-%dT%H:%M:%S.%f'
vtime = datetime.strptime(ds.Temperature_isobaric.metpy.time.dt.strftime(time_format)[0],
time_format)
# Grab lat/lon values from file as unit arrays
lats = ds.lat.metpy.unit_array
lons = ds.lon.metpy.unit_array
# Calculate distance between grid points
# will need for computations later
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
# Grabbing data for specific variable contained in file (as a unit array)
# 700 hPa Geopotential Heights
hght_700 = ds.Geopotential_height_isobaric.metpy.sel(vertical=700 * units.hPa,
time=vtime).metpy.unit_array
# Equivalent form needed if there is a dash in name of variable
# (e.g., 'u-component_of_wind_isobaric')
# hght_700 = ds['Geopotential_height_isobaric'].metpy.sel(vertical=700 * units.hPa, time=vtime).metpy.unit_array
# 700 hPa Temperature
tmpk_700 = ds.Temperature_isobaric.metpy.sel(vertical=700 * units.hPa,
time=vtime).metpy.unit_array
# 700 hPa u-component_of_wind
uwnd_700 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=700 * units.hPa,
time=vtime).metpy.unit_array
# 700 hPa v-component_of_wind
vwnd_700 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=700 * units.hPa,
time=vtime).metpy.unit_array
Write the code to access the remaining necessary pieces of data from our file to calculate the QG Omega forcing terms valid at 700 hPa.
Data variables desired:
In [ ]:
# 500 hPa Geopotential Height
# 500 hPa u-component_of_wind
# 500 hPa v-component_of_wind
# 900 hPa u-component_of_wind
# 900 hPa v-component_of_wind
In [ ]:
# %load solutions/QG_data.py
Here is the QG Omega equation from Bluesetein (1992; Eq. 5.6.11) with the two primary forcing terms on the right hand side of this equation.
$$\left(\nabla_p ^2 + \frac{f^2}{\sigma}\frac{\partial ^2}{\partial p^2}\right)\omega = \frac{f_o}{\sigma}\frac{\partial}{\partial p}\left[\vec{V_g} \cdot \nabla_p \left(\zeta_g + f \right)\right] + \frac{R}{\sigma p} \nabla_p ^2 \left[\vec{V_g} \cdot \nabla_p T \right]$$We want to write code that will calculate the differential vorticity advection term (the first term on the r.h.s.) and the laplacian of the temperature advection. We will compute these terms so that they are valid at 700 hPa. Need to set constants for static stability, f0, and Rd.
In [ ]:
# Set constant values that will be needed in computations
# Set default static stability value
sigma = 2.0e-6 * units('m^2 Pa^-2 s^-2')
# Set f-plane at typical synoptic f0 value
f0 = 1e-4 * units('s^-1')
# Use dry gas constant from MetPy constants
Rd = mpconstants.Rd
In [ ]:
# Smooth Heights
# For calculation purposes we want to smooth our variables
# a little to get to the "synoptic values" from higher
# resolution datasets
# Number of repetitions of smoothing function
n_reps = 50
# Apply the 9-point smoother
hght_700s = mpcalc.smooth_n_point(hght_700, 9, n_reps)
hght_500s = mpcalc.smooth_n_point(hght_500, 9, n_reps)
tmpk_700s = mpcalc.smooth_n_point(tmpk_700, 9, n_reps)
tmpc_700s = tmpk_700s.to('degC')
uwnd_700s = mpcalc.smooth_n_point(uwnd_700, 9, n_reps)
vwnd_700s = mpcalc.smooth_n_point(vwnd_700, 9, n_reps)
uwnd_500s = mpcalc.smooth_n_point(uwnd_500, 9, n_reps)
vwnd_500s = mpcalc.smooth_n_point(vwnd_500, 9, n_reps)
uwnd_900s = mpcalc.smooth_n_point(uwnd_900, 9, n_reps)
vwnd_900s = mpcalc.smooth_n_point(vwnd_900, 9, n_reps)
In [ ]:
# Absolute Vorticity Calculation
avor_900 = mpcalc.absolute_vorticity(uwnd_900s, vwnd_900s, dx, dy,
lats, dim_order='yx')
avor_500 = mpcalc.absolute_vorticity(uwnd_500s, vwnd_500s, dx, dy,
lats, dim_order='yx')
# Advection of Absolute Vorticity
vortadv_900 = mpcalc.advection(avor_900, (uwnd_900s, vwnd_900s), (dx, dy),
dim_order='yx').to_base_units()
vortadv_500 = mpcalc.advection(avor_500, (uwnd_500s, vwnd_500s), (dx, dy),
dim_order='yx').to_base_units()
# Differential Vorticity Advection between two levels
diff_avor = ((vortadv_900 - vortadv_500)/(400 * units.hPa)).to_base_units()
# Calculation of final differential vorticity advection term
term_A = (-f0 / sigma * diff_avor).to_base_units()
print(term_A.units)
In [ ]:
# Temperature Advection
# Laplacian of Temperature Advection
# Calculation of final Laplacian of Temperature Advection term
In [ ]:
# %load solutions/term_B_calc.py
Upper-left Panel: 700-hPa Geopotential Heights, Temperature, and Winds
Upper-right Panel: 500-hPa Geopotential Heights, Absolute Vorticity, and Winds
Lower-left Panel: Term B (Laplacian of Temperature Advection)
Lower-right Panel: Term A (Laplacian of differential Vorticity Advection)
In [ ]:
# Set some contour intervals for various parameters
# CINT 500 hPa Heights
clev_hght_500 = np.arange(0, 7000, 60)
# CINT 700 hPa Heights
clev_hght_700 = np.arange(0, 7000, 30)
# CINT 700 hPa Temps
clev_tmpc_700 = np.arange(-40, 40, 5)
# CINT Omega terms
clev_omega = np.arange(-20, 21, 2)
In [ ]:
# Set some projections for our data (Plate Carree)
# and output maps (Lambert Conformal)
# Data projection; NARR Data is Earth Relative
dataproj = ccrs.PlateCarree()
# Plot projection
# The look you want for the view, LambertConformal for mid-latitude view
plotproj = ccrs.LambertConformal(central_longitude=-100.,
central_latitude=40.,
standard_parallels=[30, 60])
In [ ]:
# Set figure size
fig=plt.figure(1, figsize=(24.5,17.))
# Upper-Left Panel
ax=plt.subplot(221, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)
# Contour #1
cs = ax.contour(lons, lats, hght_700, clev_hght_700,colors='k',
linewidths=1.5, linestyles='solid', transform=dataproj)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=3, fmt='%i',
rightside_up=True, use_clabeltext=True)
# Contour #2
cs2 = ax.contour(lons, lats, tmpc_700s, clev_tmpc_700, colors='grey',
linewidths=1.0, linestyles='dotted', transform=dataproj)
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=3, fmt='%d',
rightside_up=True, use_clabeltext=True)
# Colorfill
cf = ax.contourf(lons, lats, tadv_700*10**4, np.arange(-10,10.1,0.5),
cmap=plt.cm.bwr, extend='both', transform=dataproj)
plt.colorbar(cf, orientation='horizontal', pad=0.0, aspect=50, extendrect=True)
# Vector
ax.barbs(lons.m, lats.m, uwnd_700.to('kts').m, vwnd_700.to('kts').m,
regrid_shape=15, transform=dataproj)
# Titles
plt.title('700-hPa Geopotential Heights (m), Temperature (C),\n'
'Winds (kts), and Temp Adv. ($*10^4$ C/s)',loc='left')
plt.title('VALID: %s' %(vtime),loc='right')
# Upper-Right Panel
ax=plt.subplot(222, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)
# Contour #1
clev500 = np.arange(0,7000,60)
cs = ax.contour(lons, lats, hght_500, clev500, colors='k',
linewidths=1.5, linestyles='solid', transform=dataproj)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=3, fmt='%i',
rightside_up=True, use_clabeltext=True)
# Contour #2
cs2 = ax.contour(lons, lats, avor_500*10**5, np.arange(-40, 50, 3),colors='grey',
linewidths=1.0, linestyles='dotted', transform=dataproj)
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=3, fmt='%d',
rightside_up=True, use_clabeltext=True)
# Colorfill
cf = ax.contourf(lons, lats, vortadv_500*10**8, np.arange(-2, 2.2, 0.2),
cmap=plt.cm.BrBG, extend='both', transform=dataproj)
plt.colorbar(cf, orientation='horizontal', pad=0.0, aspect=50, extendrect=True)
# Vector
ax.barbs(lons.m, lats.m, uwnd_500.to('kts').m, vwnd_500.to('kts').m,
regrid_shape=15, transform=dataproj)
# Titles
plt.title('500-hPa Geopotential Heights (m), Winds (kt), and\n'
'Absolute Vorticity Advection ($*10^{8}$ 1/s^2)',loc='left')
plt.title('VALID: %s' %(vtime),loc='right')
# Lower-Left Panel
ax=plt.subplot(223, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)
# Contour #1
cs = ax.contour(lons, lats, hght_700s, clev_hght_700, colors='k',
linewidths=1.5, linestyles='solid', transform=dataproj)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=3, fmt='%i',
rightside_up=True, use_clabeltext=True)
# Contour #2
cs2 = ax.contour(lons, lats, tmpc_700s, clev_tmpc_700, colors='grey',
linewidths=1.0, transform=dataproj)
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=3, fmt='%d',
rightside_up=True, use_clabeltext=True)
# Colorfill
cf = ax.contourf(lons, lats, term_B*10**12, clev_omega,
cmap=plt.cm.RdYlBu_r, extend='both', transform=dataproj)
plt.colorbar(cf, orientation='horizontal', pad=0.0, aspect=50, extendrect=True)
# Vector
ax.barbs(lons.m, lats.m, uwnd_700s.to('kts').m, vwnd_700s.to('kts').m,
regrid_shape=15, transform=dataproj)
# Titles
plt.title('700-hPa Geopotential Heights (m), Winds (kt), and\n'
'Term B QG Omega ($*10^{12}$ kg m$^{-3}$ s$^{-3}$)',loc='left')
plt.title('VALID: %s' %(vtime),loc='right')
# # Lower-Right Panel
ax=plt.subplot(224, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)
# Contour #1
cs = ax.contour(lons, lats, hght_500s, clev500, colors='k',
linewidths=1.5, linestyles='solid', transform=dataproj)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=3, fmt='%i',
rightside_up=True, use_clabeltext=True)
# Contour #2
cs2 = ax.contour(lons, lats, avor_500*10**5, np.arange(-40, 50, 3), colors='grey',
linewidths=1.0, linestyles='dotted', transform=dataproj)
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=3, fmt='%d',
rightside_up=True, use_clabeltext=True)
# Colorfill
cf = ax.contourf(lons, lats, term_A*10**12, clev_omega,
cmap=plt.cm.RdYlBu_r, extend='both', transform=dataproj)
plt.colorbar(cf, orientation='horizontal', pad=0.0, aspect=50, extendrect=True)
# Vector
ax.barbs(lons.m, lats.m, uwnd_500s.to('kt').m, vwnd_500s.to('kt').m,
regrid_shape=15, transform=dataproj)
# Titles
plt.title('500-hPa Geopotential Heights (m), Winds (kt), and\n'
'Term B QG Omega ($*10^{12}$ kg m$^{-3}$ s$^{-3}$)',loc='left')
plt.title('VALID: %s' %(vtime),loc='right')
plt.show()
In [ ]:
In [ ]:
# %load solutions/qg_omega_total_fig.py