Isentropic Analysis

The MetPy function mpcalc.isentropic_interpolation allows for isentropic analysis from model analysis data in isobaric coordinates.

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, add_timestamp
from metpy.units import units

Getting the data

In this example, NARR reanalysis data for 18 UTC 04 April 1987 from the National Centers for Environmental Information (https://www.ncdc.noaa.gov/data-access/model-data) will be used.

data = xr.open_dataset(get_test_data('narr_example.nc', False))

print(list(data.variables))

Out:

['Temperature', 'time', 'isobaric', 'y', 'x', 'Lambert_Conformal', 'lat', 'lon', 'u_wind', 'v_wind', 'Geopotential_height', 'Specific_humidity']

We will reduce the dimensionality of the data as it is pulled in to remove an empty time dimension.

# Assign data to variable names
lat = data['lat']
lon = data['lon']
lev = data['isobaric']
times = data['time']

tmp = data['Temperature'][0]
uwnd = data['u_wind'][0]
vwnd = data['v_wind'][0]
spech = data['Specific_humidity'][0]

# pint doesn't understand gpm
data['Geopotential_height'].attrs['units'] = 'meter'
hgt = data['Geopotential_height'][0]

To properly interpolate to isentropic coordinates, the function must know the desired output isentropic levels. An array with these levels will be created below.

isentlevs = [296.] * units.kelvin

Conversion to Isentropic Coordinates

Once three dimensional data in isobaric coordinates has been pulled and the desired isentropic levels created, the conversion to isentropic coordinates can begin. Data will be passed to the function as below. The function requires that isentropic levels, isobaric levels, and temperature be input. Any additional inputs (in this case relative humidity, u, and v wind components) will be linearly interpolated to isentropic space.

isent_anal = mpcalc.isentropic_interpolation(isentlevs,
                                             lev,
                                             tmp,
                                             spech,
                                             uwnd,
                                             vwnd,
                                             hgt,
                                             tmpk_out=True)

The output is a list, so now we will separate the variables to different names before plotting.

isentprs, isenttmp, isentspech, isentu, isentv, isenthgt = isent_anal
isentu.ito('kt')
isentv.ito('kt')

A quick look at the shape of these variables will show that the data is now in isentropic coordinates, with the number of vertical levels as specified above.

print(isentprs.shape)
print(isentspech.shape)
print(isentu.shape)
print(isentv.shape)
print(isenttmp.shape)
print(isenthgt.shape)

Out:

(1, 118, 292)
(1, 118, 292)
(1, 118, 292)
(1, 118, 292)
(1, 118, 292)
(1, 118, 292)

Converting to Relative Humidity

The NARR only gives specific humidity on isobaric vertical levels, so relative humidity will have to be calculated after the interpolation to isentropic space.

isentrh = 100 * mpcalc.relative_humidity_from_specific_humidity(isentspech, isenttmp, isentprs)

Plotting the Isentropic Analysis

# Set up our projection
crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)

# Coordinates to limit map area
bounds = [(-122., -75., 25., 50.)]
# Choose a level to plot, in this case 296 K
level = 0

fig = plt.figure(figsize=(17., 12.))
add_metpy_logo(fig, 120, 245, size='large')
ax = fig.add_subplot(1, 1, 1, projection=crs)
ax.set_extent(*bounds, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)

# Plot the surface
clevisent = np.arange(0, 1000, 25)
cs = ax.contour(lon, lat, isentprs[level, :, :], clevisent,
                colors='k', linewidths=1.0, linestyles='solid', transform=ccrs.PlateCarree())
ax.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
          fmt='%i', rightside_up=True, use_clabeltext=True)

# Plot RH
cf = ax.contourf(lon, lat, isentrh[level, :, :], range(10, 106, 5),
                 cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', extend='max', aspect=65, shrink=0.5, pad=0.05,
                  extendrect='True')
cb.set_label('Relative Humidity', size='x-large')

# Plot wind barbs
ax.barbs(lon.values, lat.values, isentu[level, :, :].m, isentv[level, :, :].m, length=6,
         regrid_shape=20, transform=ccrs.PlateCarree())

# Make some titles
ax.set_title('{:.0f} K Isentropic Pressure (hPa), Wind (kt), Relative Humidity (percent)'
             .format(isentlevs[level].m), loc='left')
add_timestamp(ax, times[0].dt, y=0.02, high_contrast=True)
fig.tight_layout()

Montgomery Streamfunction

The Montgomery Streamfunction, \({\psi} = gdz + CpT\), is often desired because its gradient is proportional to the geostrophic wind in isentropic space. This can be easily calculated with mpcalc.montgomery_streamfunction.

# Calculate Montgomery Streamfunction and scale by 10^-2 for plotting
msf = mpcalc.montgomery_streamfunction(isenthgt, isenttmp) / 100.

# Choose a level to plot, in this case 296 K
level = 0

fig = plt.figure(figsize=(17., 12.))
add_metpy_logo(fig, 120, 250, size='large')
ax = plt.subplot(111, projection=crs)
ax.set_extent(*bounds, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=0.5)

# Plot the surface
clevmsf = np.arange(0, 4000, 5)
cs = ax.contour(lon, lat, msf[level, :, :], clevmsf,
                colors='k', linewidths=1.0, linestyles='solid', transform=ccrs.PlateCarree())
ax.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
          fmt='%i', rightside_up=True, use_clabeltext=True)

# Plot RH
cf = ax.contourf(lon, lat, isentrh[level, :, :], range(10, 106, 5),
                 cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', extend='max', aspect=65, shrink=0.5, pad=0.05,
                  extendrect='True')
cb.set_label('Relative Humidity', size='x-large')

# Plot wind barbs.
ax.barbs(lon.values, lat.values, isentu[level, :, :].m, isentv[level, :, :].m, length=6,
         regrid_shape=20, transform=ccrs.PlateCarree())

# Make some titles
ax.set_title('{:.0f} K Montgomery Streamfunction '.format(isentlevs[level].m) +
             r'($10^{-2} m^2 s^{-2}$), ' +
             'Wind (kt), Relative Humidity (percent)', loc='left')
add_timestamp(ax, times[0].dt, y=0.02, pretext='Valid: ', high_contrast=True)

fig.tight_layout()
plt.show()

Total running time of the script: ( 0 minutes 17.544 seconds)