Merge pull request #76 from kgoebber/DT_example

adds a Dynamic Tropopause example

pages/gallery/dynamic_tropopause_example.ipynb

@@ -0,0 +1,197 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Dynamic Tropopause Calculation\n",
+    "\n",
+    "By: Kevin Goebbert\n",
+    "\n",
+    "This example uses MetPy calculation ability to determine the potential temperature on the dynamic tropopause (2 PVU surface), add the derived variables to the xarray dataset and plot using the MetPy declarative syntax."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from datetime import datetime, timedelta\n",
+    "\n",
+    "import metpy.calc as mpcalc\n",
+    "from metpy.interpolate import interpolate_to_isosurface\n",
+    "from metpy.plots.declarative import *\n",
+    "from metpy.units import units\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Get GFS Data\n",
+    "\n",
+    "Obtain and subset GFS data to cover the CONUS region"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "date = datetime.utcnow() - timedelta(days=1)\n",
+    "ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/'\n",
+    "                     f'Global_onedeg_ana/GFS_Global_onedeg_ana_{date:%Y%m%d}_1200.grib2').metpy.parse_cf()\n",
+    "\n",
+    "ds = ds.sel(lat=slice(80, -10), lon=slice(360-140, 360-40))\n",
+    "\n",
+    "vtime = ds.time.values[0].astype('datetime64[ms]').astype('O')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Compute Potential Temperature at 2 PVU\n",
+    "The following cell takes the necessary data from the GFS analysis, smooths and calculates needed variables to obtain the potential temperature on the 2 PVU surface (e.g., the dynamic tropopause), as well as interpolate the wind components to that level as well."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pressure = ds.isobaric.values * units.Pa\n",
+    "\n",
+    "lats = ds.lat.values\n",
+    "lons = ds.lon.values\n",
+    "\n",
+    "dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)\n",
+    "\n",
+    "uwind = mpcalc.smooth_n_point(ds['u-component_of_wind_isobaric'].squeeze(), 9, 2)\n",
+    "vwind = mpcalc.smooth_n_point(ds['v-component_of_wind_isobaric'].squeeze(), 9, 2)\n",
+    "\n",
+    "LL_avor = mpcalc.smooth_n_point(ds.Absolute_vorticity_isobaric.metpy.sel(vertical=slice(850 * units.hPa,\n",
+    "                                                                                        925 * units.hPa)).squeeze(), 9, 2)\n",
+    "avg_LL_avor = LL_avor.mean(axis=0)\n",
+    "\n",
+    "potemp = mpcalc.smooth_n_point(mpcalc.potential_temperature(pressure[:, None, None], ds.Temperature_isobaric.squeeze()), 9, 2)\n",
+    "avg_LL_rvor = avg_LL_avor - mpcalc.coriolis_parameter(lats[:, None] * units.degree)\n",
+    "\n",
+    "pvor = mpcalc.potential_vorticity_baroclinic(potemp, pressure[:, None, None], uwind, vwind,\n",
+    "                                             dx[None, :, :], dy[None, :, :], lats[None, :, None] * units.degrees)\n",
+    "\n",
+    "DT_potemp = interpolate_to_isosurface(pvor.m*1e6, potemp.m, 2, bottom_up_search=False)\n",
+    "DT_uwnd = (interpolate_to_isosurface(pvor.m*1e6, uwind.m, 2, bottom_up_search=False) * units('m/s')).to(units.knots)\n",
+    "DT_vwnd = (interpolate_to_isosurface(pvor.m*1e6, vwind.m, 2, bottom_up_search=False) * units('m/s')).to(units.knots)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Add Variables to Dataset\n",
+    "This next cell adds the variables calculated/derived above to the xarray dataset, which will make them available for plotting with the MetPy declarative syntax."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = ds.assign(dynamic_trop=(tuple(('lat', 'lon')), DT_potemp,\n",
+    "                           {'grid_mapping': ds['u-component_of_wind_isobaric'].grid_mapping,\n",
+    "                            'units': 'PVU'}))\n",
+    "ds = ds.assign(uwnd_DT=(tuple(('lat', 'lon')), DT_uwnd.m,\n",
+    "                           {'grid_mapping': ds['u-component_of_wind_isobaric'].grid_mapping,\n",
+    "                            'units': 'knots'}))\n",
+    "ds = ds.assign(vwnd_DT=(tuple(('lat', 'lon')), DT_vwnd.m,\n",
+    "                           {'grid_mapping': ds['u-component_of_wind_isobaric'].grid_mapping,\n",
+    "                            'units': 'knots'}))\n",
+    "ds = ds.assign(avg_LL_rel_vort=(tuple(('lat', 'lon')), avg_LL_rvor*1e4,\n",
+    "                           {'grid_mapping': ds['u-component_of_wind_isobaric'].grid_mapping,\n",
+    "                            'units': '1/s'}))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Create the Plot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cntr = ContourPlot()\n",
+    "cntr.data = ds\n",
+    "cntr.level = None\n",
+    "cntr.field = 'avg_LL_rel_vort'\n",
+    "cntr.clabels = True\n",
+    "cntr.contours = [0.5, 1.5, 2.5, 3.5, 4.5]\n",
+    "\n",
+    "cntr2 = FilledContourPlot()\n",
+    "cntr2.data = ds\n",
+    "cntr2.level = None\n",
+    "cntr2.field = 'dynamic_trop'\n",
+    "cntr2.contours = list(range(250, 420, 1))\n",
+    "cntr2.colormap = 'coolwarm'\n",
+    "cntr2.colorbar = 'horizontal'\n",
+    "\n",
+    "barbs = BarbPlot()\n",
+    "barbs.data = ds\n",
+    "barbs.field = ['uwnd_DT', 'vwnd_DT']\n",
+    "barbs.skip = (3, 3)\n",
+    "\n",
+    "panel = MapPanel()\n",
+    "panel.projection = 'lcc'\n",
+    "panel.area = 'us'\n",
+    "panel.layers = ['states', 'borders', 'coastline']\n",
+    "panel.title = ('Dynamic Tropopause Potential Temperature (K), Wind Barbs (kts), and LL Rel. Vort. (s$^{-1}$) at '\n",
+    "              f'{vtime}')\n",
+    "panel.plots = [cntr2, cntr, barbs]\n",
+    "\n",
+    "pc = PanelContainer()\n",
+    "pc.size = (18, 14)\n",
+    "pc.panels = [panel]\n",
+    "\n",
+    "pc.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}