python-awips/examples/notebooks/Model_Sounding_Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EDEX modelsounding plugin creates 64-level vertical profiles from GFS and ETA (NAM) BUFR products distirubted over NOAAport. Paramters which are requestable are **pressure**, **temperature**, **specHum**, **uComp**, **vComp**, **omega**, **cldCvr**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"modelsounding\")\n",
    "forecastModel = \"ETA\"\n",
    "request.addIdentifier(\"reportType\", forecastModel)\n",
    "request.setParameters(\"pressure\",\"temperature\",\"specHum\",\"uComp\",\"vComp\",\"omega\",\"cldCvr\")\n",
    "\n",
    "request.setLocationNames(\"KDSM\")\n",
    "cycles = DataAccessLayer.getAvailableTimes(request, True)\n",
    "times = DataAccessLayer.getAvailableTimes(request)\n",
    "try:\n",
    "    fcstRun = DataAccessLayer.getForecastCycle(cycles[-1], times)\n",
    "    list(fcstRun)\n",
    "    response = DataAccessLayer.getGeometryData(request,[fcstRun[0]])\n",
    "except:\n",
    "    print('No times available')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parms    = ['uComp', 'cldCvr', 'temperature', 'vComp', 'pressure', 'omega', 'specHum']\n",
      "site     = KDSM\n",
      "geom     = POINT (-93.68000030517578 41.52000045776367)\n",
      "datetime = 1970-01-18 02:22:12 (0)\n",
      "reftime  = Jan 18 70 02:22:12 GMT\n",
      "fcstHour = 0\n",
      "period   = (Jan 18 70 02:22:12 , Jan 18 70 02:22:12 )\n"
     ]
    }
   ],
   "source": [
    "# initialize arrays\n",
    "tmp,prs,sh = np.array([]),np.array([]),np.array([])\n",
    "uc,vc,om,cld = np.array([]),np.array([]),np.array([]),np.array([])\n",
    "\n",
    "for ob in response:\n",
    "    tmp = np.append(tmp,ob.getNumber(\"temperature\"))\n",
    "    prs = np.append(prs,ob.getNumber(\"pressure\"))\n",
    "    sh = np.append(sh,ob.getNumber(\"specHum\"))\n",
    "    uc = np.append(uc,ob.getNumber(\"uComp\"))\n",
    "    vc = np.append(vc,ob.getNumber(\"vComp\"))\n",
    "    om = np.append(om,ob.getNumber(\"omega\"))\n",
    "    cld = np.append(cld,ob.getNumber(\"cldCvr\"))\n",
    "\n",
    "print(\"parms    = \" + str(ob.getParameters()))\n",
    "print(\"site     = \" + str(ob.getLocationName()))\n",
    "print(\"geom     = \" + str(ob.getGeometry()))\n",
    "print(\"datetime = \" + str(ob.getDataTime()))\n",
    "print(\"reftime  = \" + str(ob.getDataTime().getRefTime()))\n",
    "print(\"fcstHour = \" + str(ob.getDataTime().getFcstTime()))\n",
    "print(\"period   = \" + str(ob.getDataTime().getValidPeriod()))\n",
    "sounding_title = forecastModel + \" \"  + str(ob.getLocationName())  + \"(\"+ str(ob.getGeometry())+\")\" + str(ob.getDataTime())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create data arrays and calculate dewpoint from spec. humidity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.tri as mtri\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "from math import exp, log\n",
    "import numpy as np\n",
    "\n",
    "from metpy.calc import get_wind_components, lcl, dry_lapse, parcel_profile, dewpoint\n",
    "from metpy.calc import get_wind_speed,get_wind_dir, thermo, vapor_pressure\n",
    "from metpy.plots import SkewT, Hodograph\n",
    "from metpy.units import units, concatenate\n",
    "\n",
    "# we can use units.* here...\n",
    "t = (tmp-273.15) * units.degC\n",
    "p = prs/100 * units.mbar\n",
    "\n",
    "u,v = uc*1.94384,vc*1.94384 # m/s to knots\n",
    "spd = get_wind_speed(u, v) * units.knots\n",
    "dir = get_wind_dir(u, v) * units.deg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dewpoint from Specific Humidity\n",
    "\n",
    "Because the modelsounding plugin does not return dewpoint values, we must calculate the profile ourselves.  Here are three examples of dewpoint calculated from specific humidity, including a manual calculation following NCEP AWIPS/NSHARP. \n",
    "\n",
    "### 1) metpy calculated mixing ratio and vapor pressure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rmix = (sh/(1-sh)) *1000 * units('g/kg')\n",
    "e = vapor_pressure(p, rmix)\n",
    "td = dewpoint(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) metpy calculated assuming spec. humidity = mixing ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "td2 = dewpoint(vapor_pressure(p, sh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) NCEP AWIPS soundingrequest plugin\n",
    "based on GEMPAK/NSHARP, from https://github.com/Unidata/awips2-ncep/blob/unidata_16.2.2/edex/gov.noaa.nws.ncep.edex.plugin.soundingrequest/src/gov/noaa/nws/ncep/edex/plugin/soundingrequest/handler/MergeSounding.java#L1783"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/awips2/python/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in log\n",
      "/awips2/python/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: invalid value encountered in divide\n"
     ]
    }
   ],
   "source": [
    "# new arrays\n",
    "ntmp = tmp\n",
    "\n",
    "# where p=pressure(pa), T=temp(C), T0=reference temp(273.16)\n",
    "rh = 0.263*prs*sh / (np.exp(17.67*ntmp/(ntmp+273.15-29.65)))\n",
    "vaps =  6.112 * np.exp((17.67 * ntmp) / (ntmp + 243.5))\n",
    "vapr = rh * vaps / 100\n",
    "dwpc = np.array(243.5 * (np.log(6.112) - np.log(vapr)) / (np.log(vapr) - np.log(6.112) - 17.67)) * units.degC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot with MetPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABb4AAAS7CAYAAABeoax/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmYFNXVx/HfYQYYFmVVUUQRFEURFBT3JagY44JxQVEU\no6+i0ZiYuG8xLkk0cY+oGE2iIO57VFxwjyuLKIILiKCigqDINsDMff+4t6Ho6Z7pGWaqerq/n+fp\nZ6arbledPl1V3X3q9i1zzgkAAAAAAAAAgELRJOkAAAAAAAAAAACoTxS+AQAAAAAAAAAFhcI3AAAA\nAAAAAKCgUPgGAAAAAAAAABQUCt8AAAAAAAAAgIJC4RsAAAAAAAAAUFAofAMAAAAAAAAACgqFbwAA\nAAAAAABAQaHwDQAAAAAAAAAoKBS+AQAAAAAAACBPmZkzM1fX+8WKwjcAAAAAAAAAoKBQ+AYAAAAA\nAAAAFBQK3wAAAAAAAACQ32YmHUBjQ+EbAAAAAAAAAPKQmXUI/04N9y3cn55MRI0HhW8AAAAAAAAA\nyE9bhr/Twt/1ovcjhfBP0h7HxS2TDgAAAAAAAAAAkNFW4e+0LPc3iN6PFMJT84tWadIBAAAAAAAA\nAAAy6hT+3m5mO0j6INxPFbZ7hr9Tw9/OafeLFoVvAAAAAAAAAMhPf5N0sKSdJZ0cmd4y/E0VvtN7\nhBd94ZuhTgAAAAAAAAAgDznnVjjndnHOmaRTIrNuNDMn6ZZwP9tQKEWLwjcAAAAAAAAA5Dnn3B2h\nAN5Z0vi02W+a2U6i8L2KOVf0F/gEAAAAAAAAgEbHzC6SdGWGWR2cc/PjjiefUPgGAAAAAAAAgEbM\nzDpLekrSdpHJyyXt5px7L5moksVQJwAAAAAAAADQiDnnvnLObS9f7z09TG4m6V0zc2Z2o5k1TS7C\n+NHjGwAAAAAAAAAKjJl1kfRfSdtGJi+VtIdzLn2M8IJDj28AAAAAAAAAKDDOudnOud7yNeDfhMkt\nJL0XeoFfZ2alyUXYsOjxDQAAAAAAAABFwMw2kfSMpK0jkxfJ9wKflExUDYMe3wAAAAAAAABQBJxz\ns5xz28jXhX8XJreWNDH0Av97ofQCp8c3AAAAAAAAABQpM+sq3wt8q8jkhZL2dM69n0RM9YEe3wAA\nAAAAAABQpJxzM51zPeVrxX8Ik9eVNCn0At8zuejqjsI3AAAAAAAAABQ5513nnDNJ3SR9Gmbtk2BY\ndcZQJwAAAAAAAACAgkKPbwAAAAAAAABAQaHwDQAAAAAAAAAoKBS+AQAAAAAAAAAFhcI3AAAAAAAA\nAKCgUPgGAAAAAAAAABQUCt8AAAAAAAAAgIJC4RsAAAAAAAAAUFAofAMAAAAAAAAACgqFbwAAAAAA\nAABAQaHwDQAAAAAAAAAoKBS+AQAAAAAAAAAFhcI3AAAAAAAAAKCgUPgGAAAAAAAAABQUCt8AAAAA\nAAAAgIJC4RsAAAAAAAAAUFAofAMAAAAAAAAACgqFbwAAAAAAAABAQaHwDQAAAAAAAAAoKBS+AQAA\nAAAAAAAFhcI3AAAAAAAAAKCgUPgGAAAAAAAAABQUCt8AAAAAAAAAgIJC4RsAAAAAAAAAUFAofAMA\nAAAAAAAACgqFbwAAAAAAAABAQaHwDQAAAAAAAAAoKBS+AQAAAAAAAAAFhcI3AAAAAAAAAKCgUPgG\nAAAAAAAAABQUCt8AAAAAAAAAgIJC4RsAAAAAAAAAUFAofAMAAAAAAABAHjOz6Wb2QNJxNCbmnEs6\nBgAAAAAAAABAFmbmJMk5Z+F+W0mLnXMrEg0sj9HjGwAAAAAAAADyX3nk/wWSvk4qkMaAwjcAAAAA\nAAAA5CkzKwv/TkqbVRJpc6OZ9Y8vqvxH4RsAAAAAAAAA8tfW4e8kSTKzdaP3gzMl3R9nUPmOwjcA\nAAAAAAAA5K/twt9Uobt39L6ZdQz3Pw73zcycmZ0WX4j5h8I3AAAAAAAAAOSv9MJ3+v0+afc7h7/7\nNXBceY3CNwAAAAAAAADkr1Sh+4O0+9kK4en3ixKFbwAAAAAAAADIX9tJknNucbif6uE9LTpfFL7X\nQOEbAAAAAAAAAPLXOmn3U4Xw5dH7kj5Nu1/UhW9zziUdAwAAAAAAAAAgAzNzkuScsxzvfyapu6Qm\nroiLvxS+AQAAAAAAACBP1aHwvcb9YsVQJwAAAAAAAACAgkLhGwAAAAAAAADy27d1faCZtTCz0voM\npjEouicMAAAAAAAAAI2BmXUI/25gZiWSKsP9L2qxmCWpxdVbYI0APb4BAABqycyGmVlllluFme1p\nZn+spk30Ni7D8h8P866vRUwl4THXZZh3ZZh3U7jfPS2G5WY2z8zeNrNrzWyrLOvY3sweNrNvzGyF\nmf1oZhPN7DYzaxlpd09Y7g9mVpZhOZtF1n1hhvmnmtkX0cea2ZdpMf9kZm+a2TEZHl9qZmeE+T+a\n2RIz+8jMrjKzdhnav25mE9KmpdZ3c4b2+4R5h0TyXtOtwsx2zZTXyHI7mdl8MzuounaR9jub2XMh\nFwvN7EUz2zlDu6vN7P3QZoWZfW1m/zGzTXJZT1hGRzO7ycxmhu1lrpk9a2brprU7xMxeCzGVh7yf\na2ZVvneY2cDwGi0xs+/M7E5b/cUu2q6pmf3JzD43s6VhmadlibO7mT0Wtr2fzGysmfXO0vZYM5sU\nlvll2PZbprXpG5bxZXjei83sXTM7Kcsydwivw0/htXzIzDZNa5Ntm6kws9+nte1lZrea2f/MbFFo\nV2U7imyT2W435dC2wsz6puWyumU+kRbD1mb2qPljyXIzm21mI8ysfaZcVcf8MaDSzOZnmLeHmf3T\nzMab2bLQbqNaLPv1HJ/PQDO7L2x3y81sgZk9bWa7ZFluja99pO1vzWxaiH+GmV1svpCR3m59M7vb\n/P62OMS+d5Zl5rQ/pT3m55HXPn1fviJLnhamtTuphu3k95G2r1XTriI9XjNrZf796xPzx5MFZvaS\nmW0WaVPbfbTUzM42s8nm9/2FIW/9c3juqdthkbZ3mdmTGdazeXgtdsgwL6f9OrRd18z+EZ7fsrDd\n/MHMLK1dbXNb6+0lQ2w1bsdm1sXMbjCzl80flystw/t2DevJeZ83s43NfwaZY/797hszG2VmXXJc\nV07vN2Z2uvn3wG/Ddve1mf3LzDrnuJ7+oX1q214YcnRgfa8rLONKM3vKzL4K+RuZpV19rKuFmU03\nszPSprc2/zni65DbCWZ2RIbHj7MMn72gpZH/V2p1EXuSJNnqz1kfpz2u6C/sSI9vAACAunGSTlDV\nD5iS9JGkzyQ9E5m2oaRHJd0oaUxkenoRYQNJB8h/oD3WzM5xzq2sS4Dhi/EtkoZLutw5d1lak+sk\nPSDfGaKdpL6SfiXpN2Z2rnPuhsiydpD0qqTxkn4n38OkraQdJB0Z/l8SWXa5pOZh3j1p6z1R0o+S\n1k2bLjNrI+lKSec555ZFZjlJr0g6V76nygbh/1Fm1sI5d2d4fCv5vO8k6TZJl4VYdpV0tnxO93HO\nTU9bdjoXbieb2fXOuRkZ5ss5V2FVi81/CuvbV2v2qpmSYT1RV0n6yDn3VA3tFNb5sqRxkg6T/1x/\ngaRxZranc+69SPMWku4K618iaWtJf5S0j5lt45z7sYZ1dZb0enjsZZI+ld9e9pbUNNLuQPlt/BlJ\nR8h/STtE0l8lrSfpnEjbfST9N7S/SFInSddI2t7M+qdt8yMlDQ7tJsjvH/8ws1bOub9Hlrl+iPNb\nScdLWhEe86qZ9Yu+5mY2TNK/JN0q6beStgzr31JS9MRDW/n9+U5JX0pqLekYSXeYWQfn3DWRZW4j\n6SVJ70g6XFIr+W35NTPr45xbkJba+yTdkDYtvedW/xDPREkvpsUW9bakKic9JP0mxPtI2nQnv/+8\nljZ9auT/2VmWeYSkP0SXaWYbSvqfpDmSTpP0lfzx5Cr5Y0T/qovJLBSorg7LaJWhyX6Sfiafkx8l\n7ZXrsgMn6RNJx2nN/TP
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fce90d2b0d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (12, 14)\n",
    "\n",
    "# Create a skewT plot\n",
    "skew = SkewT()\n",
    "\n",
    "# Plot the data\n",
    "skew.plot(p, t, 'r', linewidth=2)\n",
    "skew.plot(p, td, 'b', linewidth=2)\n",
    "skew.plot(p, td2, 'y')\n",
    "skew.plot(p, dwpc, 'g', linewidth=2)\n",
    "\n",
    "skew.plot_barbs(p, u, v)\n",
    "skew.ax.set_ylim(1000, 100)\n",
    "skew.ax.set_xlim(-40, 60)\n",
    "\n",
    "plt.title(sounding_title)\n",
    "\n",
    "# Calculate LCL height and plot as black dot\n",
    "l = lcl(p[0], t[0], td[0])\n",
    "lcl_temp = dry_lapse(concatenate((p[0], l)), t[0])[-1].to('degC')\n",
    "skew.plot(l, lcl_temp, 'ko', markerfacecolor='black')\n",
    "\n",
    "\n",
    "# An example of a slanted line at constant T -- in this case the 0 isotherm\n",
    "l = skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)\n",
    "\n",
    "# Draw hodograph\n",
    "ax_hod = inset_axes(skew.ax, '40%', '40%', loc=2)\n",
    "h = Hodograph(ax_hod, component_range=get_wind_speed(u, v).max())\n",
    "h.add_grid(increment=20)\n",
    "h.plot_colormapped(u, v, spd)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}