Hovmoller_Diagram

Hovmoller Diagram Example

By: Kevin Goebbert

Northern Hemispheric v-wind component over the mid-latitudes in a Hovmoller diagram. This diagram can be used to illustrate upper-level wave and energy propogation (e.g., downstream baroclinic development)

In [1]:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import numpy as np
import xarray as xr

Get the data

Using NCEP/NCAR reanalysis data via xarray remote access using the OPeNDAP protocol.

Set the time range, parameter, and level to desired values

In [2]:
# Create time slice from dates
start_time = '2011-01-20'
end_time = '2011-02-06'

# Select NCEP/NCAR parameter and level
param = 'vwnd'
level = 250

# Remote get dataset using OPeNDAP method via xarray
ds = xr.open_dataset('http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/'
                     'ncep.reanalysis/pressure/{}.{}.nc'.format(param, start_time[:4]))

# Create slice variables subset domain
time_slice = slice(start_time, end_time)
lat_slice = slice(60, 40)
lon_slice = slice(0, 360)

# Get data, selecting time, level, lat/lon slice
data = ds[param].sel(time=time_slice,
                     level=level,
                     lat=lat_slice,
                     lon=lon_slice)

# Compute weights and take weighted average over latitude dimension
weights = np.cos(np.deg2rad(data.lat.values))
avg_data = (data * weights[None, :, None]).sum(dim='lat') / np.sum(weights)

# Get times and make array of datetime objects
vtimes = data.time.values.astype('datetime64[ms]').astype('O')

# Specify longitude values for chosen domain
lons = data.lon.values

Make the Hovmoller Plot

Pretty simple to use common matplotlib/cartopy to create the diagram. Cartopy is used to create a geographic reference map to highlight the area being averaged as well as the visual reference for longitude.

In [3]:
# Start figure
fig = plt.figure(figsize=(10, 13))

# Use gridspec to help size elements of plot; small top plot and big bottom plot
gs = gridspec.GridSpec(nrows=2, ncols=1, height_ratios=[1, 6], hspace=0.03)

# Tick labels
x_tick_labels = [u'0\N{DEGREE SIGN}E', u'90\N{DEGREE SIGN}E',
                 u'180\N{DEGREE SIGN}E', u'90\N{DEGREE SIGN}W',
                 u'0\N{DEGREE SIGN}E']

# Top plot for geographic reference (makes small map)
ax1 = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree(central_longitude=180))
ax1.set_extent([0, 357.5, 35, 65], ccrs.PlateCarree(central_longitude=180))
ax1.set_yticks([40, 60])
ax1.set_yticklabels([u'40\N{DEGREE SIGN}N', u'60\N{DEGREE SIGN}N'])
ax1.set_xticks([-180, -90, 0, 90, 180])
ax1.set_xticklabels(x_tick_labels)
ax1.grid(linestyle='dotted', linewidth=2)

# Add geopolitical boundaries for map reference
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES.with_scale('50m'), color='black', linewidths=0.5)

# Set some titles
plt.title('Hovmoller Diagram', loc='left')
plt.title('NCEP/NCAR Reanalysis', loc='right')

# Bottom plot for Hovmoller diagram
ax2 = fig.add_subplot(gs[1, 0])
ax2.invert_yaxis()  # Reverse the time order to do oldest first

# Plot of chosen variable averaged over latitude and slightly smoothed
clevs = np.arange(-50, 51, 5)
cf = ax2.contourf(lons, vtimes, mpcalc.smooth_n_point(
    avg_data, 9, 2), clevs, cmap=plt.cm.bwr, extend='both')
cs = ax2.contour(lons, vtimes, mpcalc.smooth_n_point(
    avg_data, 9, 2), clevs, colors='k', linewidths=1)
cbar = plt.colorbar(cf, orientation='horizontal', pad=0.04, aspect=50, extendrect=True)
cbar.set_label('m $s^{-1}$')

# Make some ticks and tick labels
ax2.set_xticks([0, 90, 180, 270, 357.5])
ax2.set_xticklabels(x_tick_labels)
ax2.set_yticks(vtimes[4::8])
ax2.set_yticklabels(vtimes[4::8])

# Set some titles
plt.title('250-hPa V-wind', loc='left', fontsize=10)
plt.title('Time Range: {0:%Y%m%d %HZ} - {1:%Y%m%d %HZ}'.format(vtimes[0], vtimes[-1]),
          loc='right', fontsize=10)
Out[3]:
Text(1.0, 1.0, 'Time Range: 20110120 00Z - 20110206 18Z')