Advanced MetPy: Quasi-Geostrophic Analysis

Objectives¶

Download NARR output from TDS
Calculate QG-Omega Forcing Terms
Create a four-panel plot of QG Forcings

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

Data from Reanalysis on pressure surfaces:¶

Geopotential Heights
Temperature
u-wind component
v-wind component

Calculations:¶

Laplacian of Temperature Advection
Differential Vorticity Advection
Wind Speed

In [1]:

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 metpy.constants as mpconstants
import numpy as np
import xarray as xr

from metpy.units import units
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS

Downloading NARR Output¶

Lets investigate what specific NARR output is available to work with from NCEI.

https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-regional-reanalysis-narr

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).

https://www.ncei.noaa.gov/thredds/catalog/narr-a-files/199303/19930313/catalog.html?dataset=narr-a-files/199303/19930313/narr-a_221_19930313_0000_000.grb

And we will use a python package called Siphon to read this data through the NetCDFSubset (NetCDFServer) link.

https://www.ncei.noaa.gov/thredds/ncss/grid/narr-a-files/199303/19930313/narr-a_221_19930313_0000_000.grb/dataset.html

First we can set out date using the datetime module

In [2]:

# 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 [3]:

# Read NARR Data from THREDDS server
base_url = 'https://www.ncei.noaa.gov/thredds/catalog/model-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
tds_ds = cat.datasets.filter_time_nearest(dt)

# Interface with the data through the NetCDF Subset Service (NCSS) 
ncss = tds_ds.subset()

# Create an NCSS query with our desired specifications
query = ncss.query()
query.lonlat_box(north=60, south=18, east=300, west=225)
query.time(dt)
query.add_lonlat()
query.accept('netcdf')
query.variables('Geopotential_height_isobaric',
                'Temperature_isobaric',
                'u-component_of_wind_isobaric',
                'v-component_of_wind_isobaric')

# Use the query to obtain our NetCDF data
data = ncss.get_data(query)

In [4]:

# Open data with xarray, and parse it with MetPy
ds = xr.open_dataset(xr.backends.NetCDF4DataStore(data)).metpy.parse_cf()
ds

Out[4]:

<xarray.Dataset>
Dimensions:                       (isobaric1: 29, time1: 1, x: 268, y: 119)
Coordinates:
  * time1                         (time1) datetime64[ns] 1993-03-13
  * isobaric1                     (isobaric1) float32 100.0 125.0 ... 1000.0
  * y                             (y) float32 -3116548.0 ... 714085.94
  * x                             (x) float32 -3324470.8 ... 5343150.5
    crs                           object Projection: lambert_conformal_conic
Data variables:
    u-component_of_wind_isobaric  (time1, isobaric1, y, x) float32 ...
    LambertConformal_Projection   int32 0
    lat                           (y, x) float64 17.87 17.96 ... 34.85 34.64
    lon                           (y, x) float64 -135.0 -134.8 ... -43.13 -42.91
    Geopotential_height_isobaric  (time1, isobaric1, y, x) float32 ...
    v-component_of_wind_isobaric  (time1, isobaric1, y, x) float32 ...
    Temperature_isobaric          (time1, isobaric1, y, x) float32 ...
Attributes:
    Originating_or_generating_Center:     US National Weather Service, Nation...
    Originating_or_generating_Subcenter:  North American Regional Reanalysis ...
    GRIB_table_version:                   0,131
    Generating_process_or_model:          North American Regional Reanalysis ...
    Conventions:                          CF-1.6
    history:                              Read using CDM IOSP GribCollection v3
    featureType:                          GRID
    History:                              Translated to CF-1.0 Conventions by...
    geospatial_lat_min:                   10.753308882144761
    geospatial_lat_max:                   46.8308828962289
    geospatial_lon_min:                   -153.88242040519995
    geospatial_lon_max:                   -42.666108129242815

In [5]:

# 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()

Subset Pressure Levels¶

Using xarray gives great funtionality for selecting pieces of your dataset to use within your script/program. MetPy also includes helpers for unit- and coordinate-aware selection and getting unit arrays from xarray DataArrays.

In [6]:

# This is the time we're using
vtime = ds['Temperature_isobaric'].metpy.time[0]

# 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)

# 700 hPa Temperature
tmpk_700 = ds['Temperature_isobaric'].metpy.sel(vertical=700 * units.hPa,
                                                time=vtime)

# 700 hPa u-component_of_wind
uwnd_700 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=700 * units.hPa,
                                                        time=vtime)

# 700 hPa v-component_of_wind
vwnd_700 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=700 * units.hPa,
                                                        time=vtime)

EXERCISE: 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:

hght_500: 500-hPa Geopotential_height_isobaric
uwnd_500: 500-hPa u-component_of_wind_isobaric
vwnd_500: 500-hPa v-component_of_wind_isobaric
uwnd_900: 900-hPa u-component_of_wind_isobaric
vwnd_900: 900-hPa v-component_of_wind_isobaric

In [7]:

# 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

SOLUTION

In [8]:

# %load solutions/QG_data.py

# Cell content replaced by load magic replacement.
# Remaining variables needed to compute QG Omega forcing terms
hght_500 = ds.Geopotential_height_isobaric.metpy.sel(vertical=500 * units.hPa,
                                                     time=vtime)
uwnd_500 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=500 * units.hPa,
                                                        time=vtime)
vwnd_500 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=500 * units.hPa,
                                                        time=vtime)
uwnd_900 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=900 * units.hPa,
                                                        time=vtime)
vwnd_900 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=900 * units.hPa,
                                                        time=vtime)

QG Omega Forcing Terms¶

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 [9]:

# 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 [10]:

# 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)#.metpy.unit_array
hght_500s = mpcalc.smooth_n_point(hght_500, 9, n_reps)#.metpy.unit_array

tmpk_700s = mpcalc.smooth_n_point(tmpk_700, 9, n_reps)#.metpy.unit_array
tmpc_700s = tmpk_700s.to('degC')

uwnd_700s = mpcalc.smooth_n_point(uwnd_700, 9, n_reps)#.metpy.unit_array
vwnd_700s = mpcalc.smooth_n_point(vwnd_700, 9, n_reps)#.metpy.unit_array

uwnd_500s = mpcalc.smooth_n_point(uwnd_500, 9, n_reps)#.metpy.unit_array
vwnd_500s = mpcalc.smooth_n_point(vwnd_500, 9, n_reps)#.metpy.unit_array

uwnd_900s = mpcalc.smooth_n_point(uwnd_900, 9, n_reps)#.metpy.unit_array
vwnd_900s = mpcalc.smooth_n_point(vwnd_900, 9, n_reps)#.metpy.unit_array

Compute Term A - Differential Vorticity Advection¶

Need to compute:

absolute vorticity at two levels (e.g., 500 and 900 hPa)
absolute vorticity advection at same two levels
centered finite-difference between two levels (e.g., valid at 700 hPa)
apply constants to calculate value of full term

In [11]:

# Absolute Vorticity Calculation
avor_900 = mpcalc.absolute_vorticity(uwnd_900s, vwnd_900s, dx, dy, lats)
avor_500 = mpcalc.absolute_vorticity(uwnd_500s, vwnd_500s, dx, dy, lats)

# Advection of Absolute Vorticity
vortadv_900 = mpcalc.advection(avor_900, (uwnd_900s, vwnd_900s), (dx, dy)).to_base_units()
vortadv_500 = mpcalc.advection(avor_500, (uwnd_500s, vwnd_500s), (dx, dy)).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)

kilogram / meter ** 3 / second ** 3

EXERCISE: Compute Term B - Laplacian of Temperature Advection Need to compute:

Temperature advection at 700 hPa (tadv_700)
Laplacian of Temp Adv. at 700 hPa (lap_tadv_700)
Final term B with appropriate constants (term_B)

For information on how to calculate a Laplacian using MetPy, see the documentation on this function.

In [12]:

# Temperature Advection


# Laplacian of Temperature Advection


# Calculation of final Laplacian of Temperature Advection term

SOLUTION

In [13]:

# %load solutions/term_B_calc.py

# Cell content replaced by load magic replacement.
# 700-hPa Temperature Advection
tadv_700 = mpcalc.advection(tmpk_700s, (uwnd_700s, vwnd_700s), (dx, dy)).to_base_units()
# Laplacian of Temperature Advection
lap_tadv_700 = mpcalc.laplacian(tadv_700, deltas=(dy, dx))

# Final term B calculation with constants
term_B = (-Rd / (sigma * (700 * units.hPa)) * lap_tadv_700).to_base_units()
print(term_B.units)

kilogram / meter ** 3 / second ** 3

Four Panel Plot¶

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 [14]:

# 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 [15]:

# 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])

Start 4-panel Figure¶

In [16]:

# Set figure size
fig=plt.figure(1, figsize=(24.5,17.))

# Format the valid time
vtime_str = str(vtime.dt.strftime('%Y-%m-%d %H%MZ').values)

# Upper-Left Panel
ax=plt.subplot(221, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
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_700s.to('kts').m, vwnd_700s.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(f'VALID: {vtime_str}', loc='right')



# Upper-Right Panel
ax=plt.subplot(222, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
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_500s.to('kts').m, vwnd_500s.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(f'VALID: {vtime_str}', loc='right')



# Lower-Left Panel
ax=plt.subplot(223, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
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(f'VALID: {vtime_str}', loc='right')



# # Lower-Right Panel
ax=plt.subplot(224, projection=plotproj)
ax.set_extent([-125., -73, 25., 50.],ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
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 A QG Omega ($*10^{12}$ kg m$^{-3}$ s$^{-3}$)',loc='left')
plt.title(f'VALID: {vtime_str}', loc='right')

plt.show()

EXERCISE:

Plot the combined QG Omega forcing terms (term_A + term_B) in a single panel.
BONUS: Compute a difference map of Term A and Term B and plot.

In [17]:

# YOUR CODE GOES HERE

SOLUTION

In [18]:

# %load solutions/qg_omega_total_fig.py

# Cell content replaced by load magic replacement.
fig=plt.figure(1, figsize=(15.,12.))

# Upper-Left Panel
ax=plt.subplot(111,projection=plotproj)
ax.set_extent([-125.,-73,25.,50.],ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
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, 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+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, Temperature (C),\n'
          'Winds (kt), and QG Omega Forcings ($*10^{12}$ kg m$^{-3}$ s$^{-3}$)',loc='left')
plt.title('VALID: ' + vtime_str, loc='right')

plt.show()