{ "cells": [ { "cell_type": "markdown", "metadata": { "cell_marker": "\"\"\"" }, "source": [ "Hovmoller Diagram Example\n", "=========================\n", "\n", "By: Kevin Goebbert\n", "\n", "Northern Hemispheric v-wind component over the mid-latitudes in a\n", "Hovmoller diagram. This diagram can be used to illustrate upper-level\n", "wave and energy propogation (e.g., downstream baroclinic development)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "import matplotlib.gridspec as gridspec\n", "import matplotlib.pyplot as plt\n", "import metpy.calc as mpcalc\n", "import numpy as np\n", "import xarray as xr" ] }, { "cell_type": "markdown", "metadata": { "cell_marker": "######################################################################" }, "source": [ "Get the data\n", "------------\n", "\n", "Using NCEP/NCAR reanalysis data via xarray remote access using the\n", "OPeNDAP protocol.\n", "\n", "Set the time range, parameter, and level to desired values\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create time slice from dates\n", "start_time = '2011-01-20'\n", "end_time = '2011-02-06'\n", "\n", "# Select NCEP/NCAR parameter and level\n", "param = 'vwnd'\n", "level = 250\n", "\n", "# Remote get dataset using OPeNDAP method via xarray\n", "ds = xr.open_dataset('http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/'\n", " f'ncep.reanalysis/pressure/{param}.{start_time[:4]}.nc')\n", "\n", "# Create slice variables subset domain\n", "time_slice = slice(start_time, end_time)\n", "lat_slice = slice(60, 40)\n", "lon_slice = slice(0, 360)\n", "\n", "# Get data, selecting time, level, lat/lon slice\n", "data = ds[param].sel(time=time_slice,\n", " level=level,\n", " lat=lat_slice,\n", " lon=lon_slice)\n", "\n", "# Compute weights and take weighted average over latitude dimension\n", "weights = np.cos(np.deg2rad(data.lat.values))\n", "avg_data = (data * weights[None, :, None]).sum(dim='lat') / np.sum(weights)\n", "\n", "# Get times and make array of datetime objects\n", "vtimes = data.time.values.astype('datetime64[ms]').astype('O')\n", "\n", "# Specify longitude values for chosen domain\n", "lons = data.lon.values" ] }, { "cell_type": "markdown", "metadata": { "cell_marker": "######################################################################" }, "source": [ "Make the Hovmoller Plot\n", "-----------------------\n", "\n", "Pretty simple to use common matplotlib/cartopy to create the diagram.\n", "Cartopy is used to create a geographic reference map to highlight the\n", "area being averaged as well as the visual reference for longitude.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Start figure\n", "fig = plt.figure(figsize=(10, 13))\n", "\n", "# Use gridspec to help size elements of plot; small top plot\n", "# and big bottom plot\n", "gs = gridspec.GridSpec(nrows=2, ncols=1, height_ratios=[1, 6], hspace=0.03)\n", "\n", "# Tick labels\n", "x_tick_labels = [u'0\\N{DEGREE SIGN}E', u'90\\N{DEGREE SIGN}E',\n", " u'180\\N{DEGREE SIGN}E', u'90\\N{DEGREE SIGN}W',\n", " u'0\\N{DEGREE SIGN}E']\n", "\n", "# Top plot for geographic reference (makes small map)\n", "ax1 = fig.add_subplot(gs[0, 0],\n", " projection=ccrs.PlateCarree(central_longitude=180))\n", "ax1.set_extent([0, 357.5, 35, 65], ccrs.PlateCarree(central_longitude=180))\n", "ax1.set_yticks([40, 60])\n", "ax1.set_yticklabels([u'40\\N{DEGREE SIGN}N', u'60\\N{DEGREE SIGN}N'])\n", "ax1.set_xticks([-180, -90, 0, 90, 180])\n", "ax1.set_xticklabels(x_tick_labels)\n", "ax1.grid(linestyle='dotted', linewidth=2)\n", "\n", "# Add geopolitical boundaries for map reference\n", "ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))\n", "ax1.add_feature(cfeature.LAKES.with_scale('50m'), color='black',\n", " linewidths=0.5)\n", "\n", "# Set some titles\n", "plt.title('Hovmoller Diagram', loc='left')\n", "plt.title('NCEP/NCAR Reanalysis', loc='right')\n", "\n", "# Bottom plot for Hovmoller diagram\n", "ax2 = fig.add_subplot(gs[1, 0])\n", "ax2.invert_yaxis() # Reverse the time order to do oldest first\n", "\n", "# Plot of chosen variable averaged over latitude and slightly smoothed\n", "clevs = np.arange(-50, 51, 5)\n", "cf = ax2.contourf(lons, vtimes, mpcalc.smooth_n_point(\n", " avg_data, 9, 2), clevs, cmap=plt.cm.bwr, extend='both')\n", "cs = ax2.contour(lons, vtimes, mpcalc.smooth_n_point(\n", " avg_data, 9, 2), clevs, colors='k', linewidths=1)\n", "cbar = plt.colorbar(cf, orientation='horizontal', pad=0.04, aspect=50,\n", " extendrect=True)\n", "cbar.set_label('m $s^{-1}$')\n", "\n", "# Make some ticks and tick labels\n", "ax2.set_xticks([0, 90, 180, 270, 357.5])\n", "ax2.set_xticklabels(x_tick_labels)\n", "ax2.set_yticks(vtimes[4::8])\n", "ax2.set_yticklabels(vtimes[4::8])\n", "\n", "# Set some titles\n", "plt.title('250-hPa V-wind', loc='left', fontsize=10)\n", "plt.title(f'Time Range: {vtimes[0]:%Y%m%d %HZ} - {vtimes[-1]:%Y%m%d %HZ}',\n", " loc='right', fontsize=10)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.5" } }, "nbformat": 4, "nbformat_minor": 4 }