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": 46,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "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",
    "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",
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"modelsounding\")\n",
    "forecastModel = \"GFS\"\n",
    "request.addIdentifier(\"reportType\", forecastModel)\n",
    "request.setParameters(\"pressure\",\"temperature\",\"specHum\",\"uComp\",\"vComp\",\"omega\",\"cldCvr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Available Locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['KFRM', 'KFSD', 'KMHE', 'KYKN', 'OAX']"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "locations = DataAccessLayer.getAvailableLocationNames(request)\n",
    "locations.sort()\n",
    "list(locations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "request.setLocationNames(\"KFRM\")\n",
    "cycles = DataAccessLayer.getAvailableTimes(request, True)\n",
    "times = DataAccessLayer.getAvailableTimes(request)\n",
    "\n",
    "try:\n",
    "    fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)\n",
    "    list(fcstRun)\n",
    "    response = DataAccessLayer.getGeometryData(request,[fcstRun[0]])\n",
    "except:\n",
    "    print('No times available')\n",
    "    exit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Sounding Parameters\n",
    "\n",
    "Construct arrays for each parameter to plot (temperature, pressure, moisutre (spec. humidity), wind components, and cloud cover)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parms    = ['uComp', 'cldCvr', 'temperature', 'vComp', 'pressure', 'omega', 'specHum']\n",
      "site     = KFRM\n",
      "geom     = POINT (-94.41999816894531 43.65000152587891)\n",
      "datetime = 1970-01-18 13:45:50.400000 (0)\n",
      "reftime  = Jan 18 70 13:45:50 GMT\n",
      "fcstHour = 15\n",
      "period   = (Jan 18 70 13:46:05 , Jan 18 70 13:46:05 )\n"
     ]
    }
   ],
   "source": [
    "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": [
    "## Calculating 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": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "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": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "td2 = dewpoint(vapor_pressure(p, sh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3) NCEP AWIPS soundingrequest plugin**\n",
    "\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": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mjames/miniconda2/envs/python-awips/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in log\n",
      "  \n",
      "/Users/mjames/miniconda2/envs/python-awips/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in divide\n",
      "  \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": [
    "## MetPy SkewT and Hodograph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAANjCAYAAAAET8T1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXl8FFW2+L83CdlIAgSSEHYI+yIg\niKyKogiIgMsoO7jgKOjM+N7M2+bNjLO8dWbebI6OOjqoo4IyCKjjwo9VVgEVBNmTEBICBLJ2tk7S\n9/fHrU4XTWch6U4n4Xw/n/4kXXep07eqbp06de45SmuNIAiCIAiCIAhNR0iwBRAEQRAEQRCE6w1R\nwgVBEARBEAShiRElXBAEQRAEQRCaGFHCBUEQBEEQBKGJESVcEARBEARBEJoYUcIFQRAEQRAEoYkR\nJVwQBEEQBEFolSil5imlCpVS2vYpUkpNs8q1Ukpb/6+3viul1Czr/w+UUnle7b0/PRsimyjhgiAI\ngiAIQmtlCBDjtS0G+C+lVDTwawClVD4wyyp3Aeut/+8G2gMO4E/AJKCN1lrZPmcaIlhYQxoJgiAI\ngiAIQgtgNKBs351AGTACKLZtbwecAvYCbwGbtdZlgRRMlHBBEARBEAShtRIJ5AJpwBHgIHASOAuc\n1loXBUswJWnrBUEQBEEQBKFpEZ9wQRAEQRAEQWhiRAkXBEEQBEEQhCZGlHBBEARBEARBaGJECRcE\nQRAEQRCEJkaUcEEQBEEQBEFoYkQJFwRBEARBEIQmRpRwQRAEQRAEQWhiRAkXBEEQBEEQhCZGlHBB\nEARBEARBaGJECRcEQRAEQRCEJkaUcEEQBEEQBEFoYkQJFwRBEARBEIQmRpRwQRAEQRAEQWhiRAkX\nBEEQBEEQhCZGlHBBEARBEARBaGJECRcEQRAEQRCEJkaUcEEQBEEQBEFoYkQJFwRBEARBEIQmRpRw\nQRAEQRAEQWhiRAkXBEEQBEEQhCZGlHBBEARBEARBaGJECRcEQRAEQRCEJkaUcEEQBEEQBEFoYkQJ\nFwRBEARBEIQmRpRwQRAEQRAEQWhiRAkXBEEQBEEQhCZGlHBBEARBEARBaGJECRcEQRAEQRCEJkaU\ncEEQBEEQBEFoYkQJFwRBEARBEK5rlKGH9X+SUupioPcpSrggCIIgCIJwXaCUaqOUusH6f6BSqsQq\nGgKcUUolAA4gQSn13UDKIkq4IAiCIAiC0GpRSvVRSv3U+joBOKiU6gBkA1FKqUe01oet8i+11sXA\n18BvAymXKOGCIAiCIAhCa+Z24MdKqaFa663Wtt1a6wIgDXjF2rYY6KqU6ghMBFBKPREooUQJFwRB\nEARBEFozQ6y//2n9fQIYoJSKA24CUEot1Fq/YZXv01oXAseBF9ydKKXethR0v6C01v7qSxAEQRAE\nQRCaFUqpzcBtQCkwSmt9VCmlga+11jcopbKALlprpZR6DHgZ6GA1zwMe0Vr/xWpzVmvdwx9yiSVc\nEARBEARBaM30s/6GA/9h/f89YJhSKgYYAaCUelBr/WerfJfWOh84A7xqbVsKdLf8yRuNWMIFQRAE\nQRCEVolSKhQoA8KsTWXAcK31CcuyvU9rPUYpdQnoaFnDVwDPAXEYxf0SsFBr/abV5pTWut/Ve7s2\nxBIuCIIgCIIgtFZ6AeW2722AX1j//zNwk1IqGhgKoJS6V2v9R6t8m9b6MnAB+Ku17dtAX6VUD6XU\nfKXUy0qpxIYIJpZwQRAEQRAEoVWilJoOvA20s20uA4ZqrU9blu3PtNa3KKUcQFvLGv4M8H/AGGAB\nUFvM8GG2EIf1RizhgiAIgiAIQmulPxDhtS0M+Ln1/4+ASUqpOcB2AEsx/z+r/HOuVMD/AtyLpaxb\nn2tWwEGUcKEZoZQarJTaH2w5rgWl1LeVUgEN5i8IQstCKXWDUmqX17b/C2S8YUEINkqpqUqpdfWs\n+x2l1H8HWiaLG4BIr21hwDxL2XYr4+8B063/M4EfA8OBEJuyrbTWj2it12mtS2gsWmv5NNMPMBfY\nCxQDF63/l+NxI1oJODHpVd2fh6yyicAuoADIBXYCN9Wwn2eBv9q+dwWOAb8HFLAV8+rGvp9xVl1t\nyecAsjBPjqG2vrZadYZ77XOdtX2ybdvfgLm27+mYcEIOjD/WX4AYW/lMzBNqMXAZeBPoZitfCuzw\n6u8C5unVve0xS8YeXr/P/rscwCQf4xYOnAW61nIMuwLrrWOQCTxRQ70l1j4fq8d50c86HvZjlgxs\nAM5Z/fS6FjmAe4DD1m/dBQy2lUUAv7H6zgOeB9rYynsBf7fKzmMWs4TZyl/CxFp1AUt9/J4+wAdA\nEWbxy/9eQ981yu21j83WuNjb2s8vB/Cp17V3HHP9XAReA+Js5U8B+zF+hiuv4Zr+iyVHX9u2v2Ky\nthUCJ+o6B+oYr61cea0e92o7H7PSvxhzDcbbyuIxN6Fiq878a2hb43hYx1Bz5fX1I1v5r4CT1u85\nBiz2au99Lf7Z69z8E+a6zgXex3Y9eu3TAVQBf7CVPwgctfb9DTDnGs77Ws8BzHl7j9c1ehYIt237\nOSYrXyXwrFd7BfwQyLDOjVVe5+ARr99WCbxvKx8BHABKrL8j6jivpljjXwJsAXp6jdMuq2xrPc7z\nWutjEqd8Yf2uVODxWvoaCnyCOde1j/Jrun5s7X5inVt32Lat5Op7amgN7XtR+3kdgYmoUYiZu/6h\nnnL5mqt6WcekxDpGd3i1ecbaR4G1z4imaOtD9v3A2Pq0xyjFmUBifcalMR/buViMuV5PYe6XG4E1\nwEJs81lTfpp8h/Kp54GBf8TcWB4AYjET8kiMohlh1VkJ/MJH2zggH5gHhAJRwFTghhr29SyWQgf0\nBE5z9Y3d58SGTaEA+mIU8WVebY8Dv7Zt62hd9BexlHDMDSoXiLTVS3dftBgl8jDw39b3BzCT2wLr\n93W2JpB0oINVZylXK+GXgX+zbXsM3zeJKxSlGn77t4CNddTZgkl72wbzRJ0L3OZVp4M1QR2uaZy9\n6n8KfMaVSngS5gFtHL6V8BrlwCj1hZgHtzDgXzGTVJhV/hNrf/FAArAH+Kmt779b52KkdRy+Br5j\nK1+BucHvx0sJxzzInAb+AWhr9XFDffquS25bHwswrxh9KeE+bypAd6CT9X8M5rr7va38PmAOJonD\nynpe0xNtctiV8CF4rumBmGtjVA191DVeW2s6h6z9FAG3WL/pLWCVrfxtYLVVNhFzUx5Sz7Y1jgce\nZSWsBrl+av3uEOBmjMI7vj7XIvBPwEHM+R8JvAGsraFuW4yidIttTnFiLF8KuBtzo06s53lf6zlg\nnXcfeG3bCDxg+77E2v96rlbCl2Dmhe7WmK8HXqvhtymMMrvYdp6cwShYEcB3rO/hNbTvZB3vb1nj\n+Etgj638Doxi/WPqp4TXWB8zBxVgFrcpTKIUB16GGlv9AcCjwGx8K+H1vn5sbVIwc8k5rlbCr7qn\n1tBHXef1f1nnTwdgkCXXtDr6rGmu2o0xcEUB92Pu7wlW2V0YXWGIta+tWPfJQLb1IftNwEmvbbW2\nx8Ti/n59xrsxH8x9cSrG2BYS6P1dk2zBFkA+Pg6KWTxQDNxfRz2fEwYwGsi/hv09i7EmpGAm6p97\nlW+lHkq49f0d4I9ebX+MeeINtbY9hblxZeJRwhcD/8+r73SvCfKXGAugsuT8J6/6IRhF9mfW96Vc\nrYT/C0YBbW9ta4wS/irw77WUx1j92Cedl4A3vOr9CaNA1zjOtrpzrTF+FpsSbisPw0sJr0sO63h8\n6DWOpcAU6/t+4Fu28vmYZAXu70eBGV7H6UUfsu3gaiX8ccyCmJp+b4191yW37Vo6AYzlGpRwH8fx\ndeDvPsp+QT2UcOu4fIl5LVqbUjkAY9V7sIbyusarxnMIkynuLdv3FIwSGotRUJ1Af1v5G3geemts\nW9d4UIey4kPODcA/2r7
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10280d1d0>"
      ]
     },
     "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",
    "# 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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}