# MetPy Case Study ## MetPy Case Study

This is a tutorial on building a case study map for Dynamic 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:¶

• Vertical Vorticity
• Advection of Temperature and Vorticity
• Horizontal Divergence
• Wind Speed
In :
from datetime import datetime

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from netCDF4 import Dataset, num2date
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
from metpy.plots import StationPlot
from metpy.units import units


### Case Study Data¶

There are a number of different sites that you can utilize to access past model output analyses and even forecasts. The most robust collection is housed at the National Center for Environmental Information (NCEI, formerly NCDC) on a THREDDS server. The general website to begin your search is

https://www.ncdc.noaa.gov/data-access

this link contains links to many different data sources (some of which we will come back to later in this tutorial). But for now, lets investigate what model output is avaiable

https://www.ncdc.noaa.gov/data-access/model-data/model-datasets

The gridded model output that are available

Reanalysis

• Climate Forecast System Reanalysis (CFSR)
• CFSR provides a global reanalysis (a best estimate of the observed state of the atmosphere) of past weather from January 1979 through March 2011 at a horizontal resolution of 0.5°.
• North American Regional Reanalysis (NARR)
• NARR is a regional reanalysis of North America containing temperatures, winds, moisture, soil data, and dozens of other parameters at 32km horizontal resolution.
• Reanalysis-1 / Reanalysis-2 (R1/R2)
• Reanalysis-1 / Reanalysis-2 are two global reanalyses of atmospheric data spanning 1948/1979 to present at a 2.5° horizontal resolution.

Numerical Weather Prediction

• Climate Forecast System (CFS)
• CFS provides a global reanalysis, a global reforecast of past weather, and an operational, seasonal forecast of weather out to nine months.
• Global Data Assimilation System (GDAS)
• GDAS is the set of assimilation data, both input and output, in various formats for the Global Forecast System model.
• Global Ensemble Forecast System (GEFS)
• GEFS is a global-coverage weather forecast model made up of 21 separate forecasts, or ensemble members, used to quantify the amount of uncertainty in a forecast. GEFS produces output four times a day with weather forecasts going out to 16 days.
• Global Forecast System (GFS)
• The GFS model is a coupled weather forecast model, composed of four separate models which work together to provide an accurate picture of weather conditions. GFS covers the entire globe down to a horizontal resolution of 28km.
• North American Mesoscale (NAM)
• NAM is a regional weather forecast model covering North America down to a horizontal resolution of 12km. Dozens of weather parameters are available from the NAM grids, from temperature and precipitation to lightning and turbulent kinetic energy.
• Rapid Refresh (RAP)
• RAP is a regional weather forecast model of North America, with separate sub-grids (with different horizontal resolutions) within the overall North America domain. RAP produces forecasts every hour with forecast lengths going out 18 hours. RAP replaced the Rapid Update Cycle (RUC) model on May 1, 2012.
• Navy Operational Global Atmospheric Prediction System (NOGAPS)
• NOGAPS analysis data are available in six-hourly increments on regularly spaced latitude-longitude grids at 1-degree and one-half-degree resolutions. Vertical resolution varies from 18 to 28 pressure levels, 34 sea level depths, the surface, and other various levels.

Ocean Models

• Hybrid Coordinate Ocean Model (HYCOM), Global
• The Navy implementation of HYCOM is the successor to Global NCOM. This site hosts regions covering U.S. coastal waters as well as a global surface model.
• Navy Coastal Ocean Model (NCOM), Global
• Global NCOM was run by the Naval Oceanographic Office (NAVOCEANO) as the Navy’s operational global ocean-prediction system prior to its replacement by the Global HYCOM system in 2013. This site hosts regions covering U.S., European, West Pacific, and Australian coastal waters as well as a global surface model.
• Navy Coastal Ocean Model (NCOM), Regional
• The Regional NCOM is a high-resolution version of NCOM for specific areas. NCEI serves the Americas Seas, U.S. East, and Alaska regions of NCOM.
• Naval Research Laboratory Adaptive Ecosystem Climatology (AEC)
• The Naval Research Laboratory AEC combines an ocean model with Earth observations to provide a synoptic view of the typical (climatic) state of the ocean for every day of the year. This dataset covers the Gulf of Mexico and nearby areas.
• National Centers for Environmental Prediction (NCEP) Real Time Ocean Forecast System (RTOFS)–Atlantic
• RTOFS–Atlantic is a data-assimilating nowcast-forecast system operated by NCEP. This dataset covers the Gulf of Mexico and most of the northern and central Atlantic.

Climate Prediction

• CM2 Global Coupled Climate Models (CM2.X)
• CM2.X consists of two climate models to model the changes in climate over the past century and into the 21st century.
• Coupled Model Intercomparison Project Phase 5 (CMIP5) (link is external)
• The U.N. Intergovernmental Panel on Climate Change (IPCC) coordinates global analysis of climate models under the Climate Model Intercomparison Project (CMIP). CMIP5 is in its fifth iteration. Data are available through the Program for Climate Model Diagnosis and Intercomparison (PCMDI) website.

Derived / Other Model Data

• Service Records Retention System (SRRS)
• SRRS is a store of weather observations, summaries, forecasts, warnings, and advisories generated by the National Weather Service for public use.
• The NOMADS Ensemble Probability Tool allows a user to query the Global Ensemble Forecast System (GEFS) to determine the probability that a set of forecast conditions will occur at a given location using all of the 21 separate GEFS ensemble members.
• National Digital Forecast Database (NDFD)
• NDFD are gridded forecasts created from weather data collected by National Weather Service field offices and processed through the National Centers for Environmental Prediction. NDFD data are available by WMO header or by date range.
• National Digital Guidance Database (NDGD)
• NDGD consists of forecasts, observations, model probabilities, climatological normals, and other digital data that complement the National Digital Forecast Database.

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

In :
# Case Study Date
year = 1993
month = 3
day = 13
hour = 0

dt = datetime(year, month, day, hour)

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)

ncss = ds.subset()
query = ncss.query().lonlat_box(north=60, south=18, east=300, west=225)
query.all_times().variables('Geopotential_height_isobaric', 'Temperature_isobaric',
'u-component_of_wind_isobaric',
data = ncss.get_data(query)

In :
# Back up in case of bad internet connection.
# Uncomment the following line to read local netCDF file of NARR data
# data = Dataset('../../data/NARR_19930313_0000.nc','r')


Let's see what dimensions are in the file:

In :
data.dimensions

Out:
{'time1': <class 'netCDF4._netCDF4.Dimension'>: name = 'time1', size = 1,
'isobaric1': <class 'netCDF4._netCDF4.Dimension'>: name = 'isobaric1', size = 29,
'y': <class 'netCDF4._netCDF4.Dimension'>: name = 'y', size = 119,
'x': <class 'netCDF4._netCDF4.Dimension'>: name = 'x', size = 268}

### Pulling Data for Calculation/Plotting¶

The object that we get from Siphon is netCDF-like, so we can pull data using familiar calls for all of the variables that are desired for calculations and plotting purposes.

NOTE: Due to the curvilinear nature of the NARR grid, there is a need to smooth the data that we import for calculation and plotting purposes. For more information about why, please see the following link: http://www.atmos.albany.edu/facstaff/rmctc/narr/

Additionally, we want to attach units to our values for use in MetPy calculations later and it will also allow for easy conversion to other units.

EXERCISE: Replace the 0's in the template below with your code:
• Use the gaussian_filter function to smooth the Temperature_isobaric, Geopotential_height_isobaric, u-component_of_wind_isobaric, and v-component_of_wind_isobaric variables from the netCDF object with a sigma value of 1.
• Assign the units of kelvin, meter, m/s, and m/s resectively.
• Extract the lat, lon, and isobaric1 variables.
In :
# Extract data and assign units
tmpk = gaussian_filter(data.variables['Temperature_isobaric'], sigma=1.0) * units.K
hght = 0
uwnd = 0
vwnd = 0

# Extract coordinate data for plotting
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lev = 0

SOLUTION
In :
# %load solutions/extract.py

# Cell content replaced by load magic replacement.
# Extract data and assign units
tmpk = gaussian_filter(data.variables['Temperature_isobaric'],
sigma=1.0) * units.K
hght = gaussian_filter(data.variables['Geopotential_height_isobaric'],
sigma=1.0) * units.meter
uwnd = gaussian_filter(data.variables['u-component_of_wind_isobaric'], sigma=1.0) * units('m/s')
vwnd = gaussian_filter(data.variables['v-component_of_wind_isobaric'], sigma=1.0) * units('m/s')

# Extract coordinate data for plotting
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lev = data.variables['isobaric1'][:]


Next we need to extract the time variable. It's not in very useful units, but the num2date function can be used to easily create regular datetime objects.

In :
time = data.variables['time1']
print(time.units)
vtime = num2date(time, units=time.units)
print(vtime)

Hour since 1993-03-13T00:00:00Z
1993-03-13 00:00:00


Finally, we need to calculate the spacing of the grid in distance units instead of degrees using the MetPy helper function lat_lon_grid_spacing.

In :
# Calcualte dx and dy for calculations
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)


### Finding Pressure Level Data¶

A robust way to parse the data for a certain pressure level is to find the index value using the np.where function. Since the NARR pressure data ('levels') is in hPa, then we'll want to search that array for our pressure levels 850, 500, and 300 hPa.

EXERCISE: Replace the 0's in the template below with your code:
• Find the index of the 850 hPa, 500 hPa, and 300 hPa levels.
• Extract the heights, temperature, u, and v winds at those levels.
In :
# Specify 850 hPa data
ilev850 = np.where(lev==850)
hght_850 = hght[ilev850]
tmpk_850 = 0
uwnd_850 = 0
vwnd_850 = 0

In :
# Specify 500 hPa data
ilev500 = 0
hght_500 = 0
uwnd_500 = 0
vwnd_500 = 0

In :
# Specify 300 hPa data
ilev300 = 0
hght_300 = 0
uwnd_300 = 0
vwnd_300 = 0

SOLUTION
In :
# %load solutions/get_850_500_300.py

# Cell content replaced by load magic replacement.
# Specify 850 hPa data
ilev850 = np.where(lev == 850)
hght_850 = hght[ilev850]
tmpk_850 = tmpk[ilev850]
uwnd_850 = uwnd[ilev850]
vwnd_850 = vwnd[ilev850]

# Specify 500 hPa data
ilev500 = np.where(lev == 500)
hght_500 = hght[ilev500]
uwnd_500 = uwnd[ilev500]
vwnd_500 = vwnd[ilev500]

# Specify 300 hPa data
ilev300 = np.where(lev == 300)
hght_300 = hght[ilev300]
uwnd_300 = uwnd[ilev300]
vwnd_300 = vwnd[ilev300]


### Using MetPy to Calculate Atmospheric Dynamic Quantities¶

MetPy has a large and growing list of functions to calculate many different atmospheric quantities. Here we want to use some classic functions to calculate wind speed, advection, planetary vorticity, relative vorticity, and divergence.

• Wind Speed: mpcalc.wind_speed()
• Advection: mpcalc.advection()
• Planetary Vorticity: mpcalc.coriolis_parameter()
• Relative Vorticity: mpcalc.vorticity()
• Divergence: mpcalc.divergence()

Note: For the above, MetPy Calculation module is imported in the following manner import metpy.calc as mpcalc.

A classic QG forcing term is 850-hPa temperature advection. MetPy has a function for advection

advection(scalar quantity, [advecting vector components], (grid spacing components))

So for temperature advection our scalar quantity would be the tempertaure, the advecting vector components would be our u and v components of the wind, and the grid spacing would be our dx and dy we computed in an earier cell.

EXERCISE:
• Uncomment and fill out the advection calculation below.
In :
# Temperature Advection

SOLUTION
In :
# %load solutions/advection_850.py

# Cell content replaced by load magic replacement.
(dx, dy), dim_order='yx').to('degC/s')


#### Vorticity Calculations¶

There are a couple of different vorticities that we are interested in for various calculations, planetary vorticity, relative vorticity, and absolute vorticity. Currently MetPy has two of the three as functions within the calc module.

Planetary Vorticity (Coriolis Parameter)

coriolis_parameter(latitude in radians)

Note: You must can convert your array of latitudes to radians...NumPy give a great function np.deg2rad() or have units attached to your latitudes in order for MetPy to convert them for you! Always check your output to make sure that your code is producing what you think it is producing.

Relative Vorticity

When atmospheric scientists talk about relative vorticity, we are really refering to the relative vorticity that is occuring about the vertical axis (the k-hat component). So in MetPy the function is

vorticity(uwind, vwind, dx, dy)

Absolute Vorticity

Currently there is no specific function for Absolute Vorticity, but this is easy for us to calculate from the previous two calculations because we just need to add them together!

ABS Vort = Rel. Vort + Coriolis Parameter

Here having units are great, becase we won't be able to add things together that don't have the same units! Its a nice safety check just in case you entered something wrong in another part of the calculation, you'll get a units error.

EXERCISE:
• Fill in the function calls below to complete the vorticity calculations.
In :
# Vorticity and Absolute Vorticity Calculations

# Planetary Vorticity
# f = mpcalc.coriolis_parameter(-- Fill in here --).to('1/s')

# Relative Vorticity
# vor_500 = mpcalc.vorticity(-- Fill in here --)

# Abosolute Vorticity
# avor_500 = vor_500 + f

SOLUTION
In :
# %load solutions/vort.py

# Cell content replaced by load magic replacement.
# Vorticity and Absolute Vorticity Calculations

# Planetary Vorticity

# Relative Vorticity
vor_500 = mpcalc.vorticity(uwnd_500, vwnd_500, dx, dy,
dim_order='yx')

# Abosolute Vorticity
avor_500 = vor_500 + f


We use the same MetPy function for temperature advection for our vorticity advection, we just have to change the scalar quantity (what is being advected) and have appropriate vector quantities for the level our scalar is from. So for vorticity advections well want our wind components from 500 hPa.

In :
# Vorticity Advection



#### Divergence and Stretching Vorticity¶

If we want to analyze another component of the vorticity tendency equation other than advection, we might want to assess the stretching forticity term.

-(Abs. Vort.)*(Divergence)

We already have absolute vorticity calculated, so now we need to calculate the divergence of the level, which MetPy has a function

divergence(uwnd, vwnd, dx, dy)

This function computes the horizontal divergence.

In :
# Stretching Vorticity
div_500 = mpcalc.divergence(uwnd_500, vwnd_500, dx, dy, dim_order='yx')

stretch_vort = -1 * avor_500 * div_500


#### Wind Speed, Geostrophic and Ageostrophic Wind¶

Wind Speed

Calculating wind speed is not a difficult calculation, but MetPy offers a function to calculate it easily keeping units so that it is easy to convert units for plotting purposes.

wind_speed(uwnd, vwnd)

Geostrophic Wind

The geostrophic wind can be computed from a given height gradient and coriolis parameter

geostrophic_wind(heights, coriolis parameter, dx, dy)

This function will return the two geostrophic wind components in a tuple. On the left hand side you'll be able to put two variables to save them off separately, if desired.

Ageostrophic Wind

Currently, there is not a function in MetPy for calculating the ageostrophic wind, however, it is again a simple arithmatic operation to get it from the total wind (which comes from our data input) and out calculated geostrophic wind from above.

Ageo Wind = Total Wind - Geo Wind

In :
# Divergence 300 hPa, Ageostrophic Wind
wspd_300 = mpcalc.wind_speed(uwnd_300, vwnd_300).to('kts')

div_300 = mpcalc.divergence(uwnd_300, vwnd_300, dx, dy, dim_order='yx')
ugeo_300, vgeo_300 = mpcalc.geostrophic_wind(hght_300, f, dx, dy, dim_order='yx')

uageo_300 = uwnd_300 - ugeo_300
vageo_300 = vwnd_300 - vgeo_300


### Maps and Projections¶

In :
# 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 :
def create_map_background():
fig=plt.figure(figsize=(14, 12))
ax=plt.subplot(111, projection=plotproj)
ax.set_extent([-125, -73, 25, 50],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
return fig, ax


• Add one contour (Temperature in Celsius with a dotted linestyle
EXERCISE:
• Add one contour (Temperature in Celsius with a dotted linestyle
In :
fig, ax = create_map_background()

# Contour 1 - Temperature, dotted

# Contour 2
clev850 = np.arange(0, 4000, 30)
cs = ax.contour(lon, lat, hght_850, clev850, colors='k',
linewidths=1.0, linestyles='solid', transform=dataproj)

plt.clabel(cs, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
rightside_up=True, use_clabeltext=True)

# Filled contours - Temperature advection
contours = [-3, -2.2, -2, -1.5, -1, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]

# Vector
ax.barbs(lon, lat, uwnd_850.to('kts').m, vwnd_850.to('kts').m,
regrid_shape=15, transform=dataproj)

# Titles
plt.title('850-hPa Geopotential Heights, Temperature (C), \
Temp Adv (C/h), and Wind Barbs (kts)', loc='left')
plt.title(f'VALID: {vtime}', loc='right')

plt.tight_layout()
plt.show()

/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/ipykernel_launcher.py:27: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations. SOLUTION
In :
# %load solutions/temp_adv_map_850.py

# Cell content replaced by load magic replacement.
fig, ax = create_map_background()

# Contour 1 - Temperature, dotted
cs2 = ax.contour(lon, lat, tmpk_850.to('degC'), range(-50, 50, 2),
colors='grey', linestyles='dotted', transform=dataproj)

plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
rightside_up=True, use_clabeltext=True)

# Contour 2
clev850 = np.arange(0, 4000, 30)
cs = ax.contour(lon, lat, hght_850, clev850, colors='k',
linewidths=1.0, linestyles='solid', transform=dataproj)

plt.clabel(cs, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
rightside_up=True, use_clabeltext=True)

# Filled contours - Temperature advection
contours = [-3, -2.2, -2, -1.5, -1, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]
cf = ax.contourf(lon, lat, tmpc_adv_850*3600, contours,
cmap='bwr', extend='both', transform=dataproj)
extendrect=True, ticks=contours)

# Vector
ax.barbs(lon, lat, uwnd_850.to('kts').m, vwnd_850.to('kts').m,
regrid_shape=15, transform=dataproj)

# Titles
plt.title('850-hPa Geopotential Heights, Temperature (C), \
Temp Adv (C/h), and Wind Barbs (kts)', loc='left')
plt.title(f'VALID: {vtime}', loc='right')

plt.tight_layout()
plt.show()

/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/ipykernel_launcher.py:37: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations. ### 500-hPa Absolute Vorticity¶

EXERCISE:
• Add code for plotting vorticity as filled contours with given levels and colors.
In :
fig, ax = create_map_background()

# Contour 1
clev500 = np.arange(0, 7000, 60)
cs = ax.contour(lon, lat, hght_500, clev500, colors='k',
linewidths=1.0, linestyles='solid', transform=dataproj)
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=4,
fmt='%i', rightside_up=True, use_clabeltext=True)

# Filled contours
# Set contour intervals for Absolute Vorticity
clevavor500 = [-4, -3, -2, -1, 0, 7, 10, 13, 16, 19,
22, 25, 28, 31, 34, 37, 40, 43, 46]

# Set colorfill colors for absolute vorticity
# purple negative
# yellow to orange positive
colorsavor500 = ('#660066', '#660099', '#6600CC', '#6600FF',
'#FFFFFF', '#ffE800', '#ffD800', '#ffC800',
'#ffB800', '#ffA800', '#ff9800', '#ff8800',
'#ff7800', '#ff6800', '#ff5800', '#ff5000',
'#ff4000', '#ff3000')

# Vector
ax.barbs(lon, lat, uwnd_500.to('kts').m, vwnd_500.to('kts').m,
regrid_shape=15, transform=dataproj)

# Titles
plt.title('500-hPa Geopotential Heights, Absolute Vorticity \
(1/s), and Wind Barbs (kts)', loc='left')
plt.title(f'VALID: {vtime}', loc='right')

plt.tight_layout()
plt.show()

/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/ipykernel_launcher.py:37: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations.