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": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['CHE',\n",
       " 'CRL',\n",
       " 'EAX',\n",
       " 'HSI',\n",
       " 'KDSM',\n",
       " 'KFOE',\n",
       " 'KFRM',\n",
       " 'KFSD',\n",
       " 'KGRI',\n",
       " 'KLNK',\n",
       " 'KMCI',\n",
       " 'KMCW',\n",
       " 'KMHE',\n",
       " 'KMHK',\n",
       " 'KMKC',\n",
       " 'KOFK',\n",
       " 'KOMA',\n",
       " 'KRSL',\n",
       " 'KSLN',\n",
       " 'KSTJ',\n",
       " 'KSUX',\n",
       " 'KTOP',\n",
       " 'KYKN',\n",
       " 'OAX',\n",
       " 'P#8',\n",
       " 'P#9',\n",
       " 'P#A',\n",
       " 'P#G',\n",
       " 'P#I',\n",
       " 'RDD',\n",
       " 'WSC']"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "\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\")\n",
    "locations = DataAccessLayer.getAvailableLocationNames(request)\n",
    "locations.sort()\n",
    "list(locations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request.setLocationNames(\"KMCI\")\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": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parms    = ['uComp', 'cldCvr', 'temperature', 'vComp', 'pressure', 'omega', 'specHum']\n",
      "site     = KMCI\n",
      "geom     = POINT (-94.72000122070312 39.31999969482422)\n",
      "datetime = 1970-01-18 04:45:50.400000 (0)\n",
      "reftime  = Jan 18 70 04:45:50 GMT\n",
      "fcstHour = 0\n",
      "period   = (Jan 18 70 04:45:50 , Jan 18 70 04:45:50 )\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": [
    "## Create data arrays and calculate dewpoint from spec. humidity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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": 37,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "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": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mj/miniconda2/lib/python2.7/site-packages/MetPy-0.3.0+34.gcf954c5-py2.7.egg/metpy/calc/thermo.py:371: RuntimeWarning: divide by zero encountered in log\n",
      "  val = np.log(e / sat_pressure_0c)\n",
      "/Users/mj/miniconda2/lib/python2.7/site-packages/pint/quantity.py:1236: RuntimeWarning: divide by zero encountered in log\n",
      "  out = uf(*mobjs)\n",
      "/Users/mj/miniconda2/lib/python2.7/site-packages/pint/quantity.py:693: RuntimeWarning: invalid value encountered in true_divide\n",
      "  magnitude = magnitude_op(self._magnitude, other_magnitude)\n"
     ]
    }
   ],
   "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": 39,
   "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": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mj/miniconda2/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in log\n",
      "/Users/mj/miniconda2/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": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mj/miniconda2/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCQAAANoCAYAAAAYqzuQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8lOW9///XlYUdAsoWRNlFwQ1BKrUtSBUrIK4VRBEV\n0HqO3/bYo+fYnp9F29oqVo+ttWhBDlVBpFVkjYiyiEJFgUpYZF9CCBCWkITsM9fvj/sOmSQzWSCZ\n+07yfj4e80hm7u0z173M3J+5FmOtRUREREREREQkmmK8DkBEREREREREGh4lJEREREREREQk6pSQ\nEBEREREREZGoU0JCRERERERERKJOCQkRERERERERiTolJEREREREREQk6pSQEBEREREREannjDHP\nGGOs13GEMtb6Kh4RERERERERaQBUQ0JEREREREREok4JCRERERERERGJOiUkRERERERERCTqlJAQ\nERERERERkahTQkJEREREREREok4JCRERERERERGJOiUkRERERERERCTqlJAQERERERERkahTQkJE\nREREREREok4JCRERERERERGJOiUkRERERERERCTqlJAQERERERERkahTQkJEREREREREok4JCRER\nERERERGJOiUkRERERERERCTqlJAQERERERERkahTQkJEREREREREok4JCRERERERERGJOiUkRERE\nRERERCTqlJAQERERERERkahTQkJEREREREREok4JCRERERERERGJOiUkRERERERERCTqlJAQERER\nERERkahTQkJERERERETE54xjjDEm1utYaooSEiIiIiIiIiL+Nxl4F2jtdSA1xVhrvY5BRERERERE\nRCIwxowFZgCxQCNbT27kVUNCRERERERExKeMMd8DpgGNgRP1JRkBSkiIiIiIiIiI+JIxpiewGGjm\nvnTYw3BqnBISIiIiIiIiIj5jjDkfWAm0CHk5xZtoaocSEiIiIiIiIiI+YoxpDHwMtKP0ffsud3on\nY0xLL2KrSUpIiIiIiIiIiPiEMcYAs4FLgUYhkwqA/e7/qcCsKIdW4+K8DkBEREREREREzpgE3Eb5\nCgT5OImIYjujFlEtUQ0JEREREREREf9YBLwB5LiPYkFKJyTWRTOo2qCEhIiIiIiIiIhPWGsPWWv/\nDegIPAUcArKBBOBUyKx1PiFh6tEQpiIiIiIiIiL1ijEmBrgJWFJmUoyt4zf0qiEhIiIiIiIi4lPW\n2qC1Nslaa4A7QiYFjTGPeRVXTVANCREREREREZE6xBjTFfgC6OS+tBC4x1p72quYzoZqSIiIiIiI\niIjUIdbafdbaC4DGwEzgFiDbGJNnjOnraXDVoISEiIiIiIiISB1krS2w1j7oNue4DydBsdkYY40x\nD3gbXeXUZENERERERESknjDG9MYZgaOV+9K7wIPW2nzvogpPCQkRERERERGResYY0xR4C7jLfek4\ncK21dpd3UZWmhISIiIiIiIhIPWaMeRh4A8Bt3uELSkiIiIiIiIiISNSpU0sRERERERERiTolJERE\nREREREQk6pSQEBERqSZjzO+NMT/zOo7qMMZ8YYy50us4RMRfjDH/MMb8KOR5B2PMVmNMIy/jEqlt\nxphHjDH/W8V5HzPGPF/bMTVESkiIiEiNMsaMMcZ8aYzJNsYcMcb80xjzaMj0mcaYfGNMVsjjx+60\nCcaYbcaYTGPMYWPMYmNMiwjbWWmMmRDyfIgx5oQx5m73edDdfmzIPPHGmKPGmGCZdd1kjPnM3e5R\nd923uNMeMMasDpm3HTAOeD1ku0H3fWQaY74NHffbGNPYTWDsN8bkGGN2GGOeiPReQtb3Wpl5PjfG\njDfG/CKk3HKNMUUhz5MjlNUtwClr7TfhphfPY4zZ7K7nC2PMpRHm+9SNL+x3CGPM98vs2yx3/tvd\n6eONMV8bY04ZY1KMMS+U2UfnGWPmucfPPmPMPWXW/0O3jE8bY5YbYy4qM/0FY8wx9/F8mWm/McYk\nG2MKjTGTy0wb4ZbxSWNMmjFmWuixZ4z5g7vvMt1jdFyZ5a8yxqx34/rahCR/jDGvlymPPGNMZlXe\nszGmj7u+E8aYDHfffC9k+vXGmBXutL1lYmpnjHnXGJPqTv/cGDMw3H5z5x/jlu0pt/w+MMZ0Cpl+\nqVvmGcaYncaY285hXY+57yvPGPN/YZaf6G4jyxiTZIxJDJnW2hjzN+Oc30fC7MvvGmPWufvqG2PM\ndWWm/49xzsdTbvm0LDP9BmPMBnd/pBj3+uROO+vzxBjT2Riz0Bhz3D3GXi0+9o0x1xpjlrnTjhpj\n5hpjOoZZbyP3+EsJea3C/VzZsQ28APy2+Im19giwAng4ZB2XGWOWGmPSTZnrpzs94rFhjLm3zPF/\n2i2bfiHzRDxvI5TvWZ3nEdZV4fzGmIeNMbvc4+WrssdTmXkrvH6FzPcrtwyGViG+cvvcfX2fcT5T\nisv1owrWEfo5VfwYFzK9sTFmhvse04wxj1cWl7vcDHe93au6LlPBtdKd/ri73CljzJsmJDFWWfma\nSj4fyszbCPgfYEoVY5sG3Guc7wBSk6y1euihhx566FEjD+A/gcPAHUBz97WrgHeAePf5/wG/DrPs\nYHfZK93nbXBu/FtE2NYK4CH3/2HASWBUyPQgsA0YGfLaKOBbIBDy2l3AKeAhoKX72g+Av7r/PwCs\nDpn/SeCNkOdDgJSQ57cChcAl7vMFwD+BPjg/BHwH2AH8McJ7GQJkuTF1CZlnNXB/mTIYD3xWhf2y\nGLingum93O19143xKWAnEFtmvnuBVUAAiKniMTEYyASaus9/AlwHxAGdgK+B/w6Z/1330cydLwPo\n405r6z6/E2iE80Vybciyj7j7t5P72AI8EjL9fuBHwIfAr8rEeY97HDUBWgNLgKkh058BLnb/Hwic\nAAa5zxsB+4GfAfHA/wP24R7zYcrk/4DpVXzPCUA3wLiP/wccDln2Gne/TAL2ltlON+A/gA7uspOA\ndNxzM0xcFwLt3f+b45y3c9zncTjH7X+467oeyAZ6VXdd7mu345wrfwH+r8yyQ4AjwKVuef4FWFmm\n/N5z91UXYBfwgDvtPJxh7e5047zX3VetQ86ZbcAFblwfAjND1t3H3fZNOOdCG6B7TZwnwAdu7I3c\nfbIJ+H/utB+5MbcAmgJvAklhyvV/3HUfqMJ+blaVY9udZwfQP+T5d4HkkOcXAw/iXEODZZat7rEx\nHthZ1fM2zPJnfZ5HWF9F14WrcK7H/UKuX0dxBwYIs66I53LIPD3cfX8QGFqF+Mrtc/f1vVVZPuSc\nSqlg+u/dbSQAlwBpwE2VrPN7wEr3OO9elXVRybUS57w7jHPut8b5bPx9VcqXSj4fwsT/Y2BpyPNK\nr+PAX4H/rEqZ61H1h+cB6KGHHnroUT8e7pePbOD2SuaLlJB4AphXje2tACYAI3GSET8qMz0I/BKY\nG/LaP9zXgu5zAxyo6AsG5RMSnwJjQ56X+6KH84X1DuCHQC5wQZnpA4EiSm50yiYkUoA/AjNClgmX\nkCgVW4T4GwE5QKcK5nkMWBTy3LjLDA15LQHYjpNQCVL1hMT/AW9WMP1xYIH7f3MgH+gZMv1vxV9I\ncX6x/TxkWjM3zuJEwRpgYsj0B8N9IQXeBiZXEvftwKYKps8HHnf/HwYcLDN9P2G+0LvvMRP4flXe\nc5ll44B/BzaGmXYDZRISEeI+hXtzVcl8Ldw4/td9fhmQVWaepYQ5lytbV5lpv6F8QuIPwJ9Dnie6\nx1w393k6MCBk+i9wE3M414MtZda3PeT8+gfwRMi0QTjnaBP3+Wzg2do4T9zXfxTyfArweoRtXQ1k\nlnmtG7AV5+Y54s1lZfs53LGNc6P1q5DnccBp4MIy8/WkfEKiWscGzvXu6ZDnVTpvqzs/VTjPK5sf\nGAt8GfK8ubtfO4RZvkrnMpAE3EwVEgoV7XN3+R9W8b0NqeiYAVKBG0KePwu8W8H8ccAG4HK3PLpX\nZV1EvlYOc/+fDfw2ZNr
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107320210>"
      ]
     },
     "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}