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": 1,
   "metadata": {
    "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 wind_speed, wind_direction, thermo, vapor_pressure\n",
    "from metpy.plots import SkewT, Hodograph\n",
    "from metpy.units import units, concatenate\n",
    "\n",
    "DataAccessLayer.changeEDEXHost(\"149.165.156.89\")\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": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[b'CHE',\n",
       " b'CRL',\n",
       " b'EAX',\n",
       " b'HSI',\n",
       " b'KDSM',\n",
       " b'KFOE',\n",
       " b'KFRM',\n",
       " b'KFSD',\n",
       " b'KGRI',\n",
       " b'KLNK',\n",
       " b'KMCI',\n",
       " b'KMCW',\n",
       " b'KMHE',\n",
       " b'KMHK',\n",
       " b'KMKC',\n",
       " b'KOFK',\n",
       " b'KOMA',\n",
       " b'KRSL',\n",
       " b'KSLN',\n",
       " b'KSTJ',\n",
       " b'KSUX',\n",
       " b'KTOP',\n",
       " b'KYKN',\n",
       " b'OAX',\n",
       " b'P#8',\n",
       " b'P#9',\n",
       " b'P#A',\n",
       " b'P#G',\n",
       " b'P#I',\n",
       " b'RDD',\n",
       " b'WSC']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "locations = DataAccessLayer.getAvailableLocationNames(request)\n",
    "locations.sort()\n",
    "list(locations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parms    = [b'temperature', b'pressure', b'vComp', b'uComp', b'cldCvr', b'specHum', b'omega']\n",
      "site     = b'KFRM'\n",
      "geom     = POINT (-94.41999816894531 43.65000152587891)\n",
      "datetime = 2018-09-06 06:00:00\n",
      "reftime  = Sep 06 18 06:00:00 GMT\n",
      "fcstHour = 0\n",
      "period   = (Sep 06 18 06:00:00 , Sep 06 18 06:00:00 )\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.getString(b\"temperature\"))\n",
    "    prs = np.append(prs,ob.getString(b\"pressure\"))\n",
    "    sh = np.append(sh,ob.getString(b\"specHum\"))\n",
    "    uc = np.append(uc,ob.getString(b\"uComp\"))\n",
    "    vc = np.append(vc,ob.getString(b\"vComp\"))\n",
    "    om = np.append(om,ob.getString(b\"omega\"))\n",
    "    cld = np.append(cld,ob.getString(b\"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()))"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "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 = wind_speed(u, v) * units.knots\n",
    "dir = wind_direction(u, v) * units.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/metpy/calc/thermo.py:690: RuntimeWarning: divide by zero encountered in log\n",
      "  val = np.log(e / sat_pressure_0c)\n",
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/pint/quantity.py:1403: RuntimeWarning: divide by zero encountered in log\n",
      "  out = uf(*mobjs)\n",
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/pint/quantity.py:802: RuntimeWarning: invalid value encountered in true_divide\n",
      "  magnitude = magnitude_op(new_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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/metpy/calc/thermo.py:690: RuntimeWarning: divide by zero encountered in log\n",
      "  val = np.log(e / sat_pressure_0c)\n",
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/pint/quantity.py:1403: RuntimeWarning: divide by zero encountered in log\n",
      "  out = uf(*mobjs)\n",
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/pint/quantity.py:802: RuntimeWarning: invalid value encountered in true_divide\n",
      "  magnitude = magnitude_op(new_self._magnitude, other._magnitude)\n"
     ]
    }
   ],
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in log\n",
      "  \n",
      "/Users/mjames/miniconda2/envs/python3-awips/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvcAAANjCAYAAAAu/zVaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd8VHW6+PHPNwUS0iG0kJDQpCodRQGliij2gqiIXXfX3+7d7r3rvdvu6t57d+9et6jIqqhYEAuuYkdUdAVRsaBJICEJCem9l5nv74/vmZwhJiRAkjPleb9e80pmzsyZZ86cmfnOM895vkprjRBCCCGEEML/hTgdgBBCCCGEEKJ3yOBeCCGEEEKIACGDeyGEEEIIIQKEDO6FEEIIIYQIEDK4F0IIIYQQIkDI4F4IIYQQQogAIYN7IYQQQggRUJRSw5VSWim1zDqvlVKPef3/rV7wSqkwpdQspdSdSqmnlVL5nut2c0rt78d3LGFOByCEEEIIIURv0loXK6UA3gQUsBW4Tin1KpANjO1sgN+JQmAX8IF1+lxr3do3UfcOJZNYCSGEEEKIQKOUOgP4ZzdX24k9cP+n1rqqr+PqazK4F0IIIYQQAUkpdQUQiRm8Z+sgGPjK4F4IIYQQQogAIQfUCiGEEEIIESBkcC+EEEIIIUSAkMG9EEIIIYQQAUIG90IIIYQQQgQIGdwLIYQQQggRIGRwL4QQQgghRICQwb0QQgghhBABQgb3QgghhBBCBAgZ3AshhBBCCBEgZHAvhBBCCCFEgJDBvRBCCCGEEAFCBvdCCCGEEEIECBncCyGEEEIIESBkcC+EEEIIIUSAkMG9EEIIIYQQAUIG90IIIYQQQgQIGdwLIYQQQggRIGRwL4QQQgghRICQwb0QQgghhBABQgb3QgghhBBCBAgZ3AshhBBCCBEgZHAvhBBCCCFEgJDBvRBCCCGEEAFCBvdCCCGEEEIECBncCyGEEEIIESBkcC+EEEIIIUSAkMG9EEIIIYQQAUIG90IIIYQQQgQIGdwLIYQQQggRIGRwL4QQQgghRICQwb0QQgghhAhoSinldAz9RQb3QgghhBDC7yljsPV/qlJqu/X/dMCtlBrgaID9RAb3QgghhBAiEJwPlFuDeDdwnlJqEfCFtfwZxyLrR0pr7XQMQgghhBBCnBSr9MYNvKi1vkQppQG01kop9Wfge0Co1trtZJx9TTL3QgghhBDC72mTsf4rcLFSKgQ4DUApNQP4vnW1BxwKr99I5l4IIYQQQgQEa1DvAjZorW/rkL1/BrgSCNEBPACWzL0QQgghhAgIVsnNFuBWq0znLACl1HjgWutqv3UovH4hmXshhBBCCBEwlFLhQAvwW6313Vb2XmutQ5RSO4DFWuuTao2plBoBnAcsAO7SWpecdOC9RAb3QgghhBAioCil3gHOscpxVgGvAElANVAP/IvW+k/drCMROBdYZZ3iu7jqqVrrr3ot+JMkg3shhBBCCBFQlFKDMIP472ut77Oy9yVa6+FKqS+BadbAPw5Yjj2AH36M1X4KbMd8UfhYa+3q20dxYmRwLwSglJoCPKa1nuN0LD2llLoNmKy1/oHTsQghfINS6jTgAa31mV6X/RHI1FoHfJcQIbwppfYDU4Bo4O/AVUAJMOwYN/sKM4DfDnyotW7t6zh7ndZaTnL61glYA+zGfOstsf7/DvYXwkcx9Wx1XqerrGULgA8xP31VAB8Ac7u4n18CT3idHwWkA/cBCtgJNHW4n/nWdbUVXx1QAPwR07/Ws66d1nWmd7jPF63Lz/G67Dlgjdf5HKDRWncx8AgQ7bX8AmCPdf/lwGYg2Wv5emBXh/UVA1Fel91sxTi6w+Pzflx1wMJOttsA4DAw6hjP4Shgm/Uc5AO3d3G96637vLkH+8UE6/nwfs5GAi8BR6z1pB1PHMBqzJtpnbXfTPFaNhD4X2vdlcDfgHCv5WmYN+BKoAj4CxDmtXwDkIHpe7y+k8czFngZqAXKgP86jnV3GXeH+9hhbRfv23rvX3XAGx1eexmY108JsAmI9Vr+PWAv0Aw8ehyv6UesOMZ7XfYEUAjUAJnd7QPdbK+dHP1azehw27VALmbffhEY7LVsMPCCtSwXWHsct+1ye1jPoebo19fdXsv/BzhgPZ50YF2H23d8LW7ssG8+gHldVwD/wOv12OE+6zDdO/7stfxK4Bvrvr8GLj6O/f6Y+wBmv13d4TV6GBjQw33lfGAXUIXZ9x8CYjrE97C13xQBP+zw3rQVs48f9T7bk+3WRTz/Yt1PtXW/A72WnYl5L67FTFS0oJt1pQHvAA3Wc76sp/t4H6/rhF8Dfbkufz8BCdZ+2PF0AHgaUzMf4XScvf64nQ5ATr53An5kvfFeDsRgBtkzMQPYgdZ1HsUcqNLxtrHWB8LVQCgQCawATuvivn6JNVAEUoEsvj1g6HTAgddABRiPGeDf0uG2GcAfvC4bYn1IlGB96GA++Cq8X+CYD6Zl1v+jMIO4e63zl2M+1K6xHt8IzAdODpBgXWc93x7clwP/6nXZzcDOYz2uYzxHVwBvdnOdd4A/AeHAdOsxLu5wnQTMh9JXXW3nDtd/A3ifowf3wzFf/ObT+eC+yzgwXxZqMF8Iw4C7gINYA2HgP6z7GwwMBT4CfuW17u3WvhhhPQ9fAv/Pa/l3gaWYgdD6DnENsPa3HwJR1jpO68m6u4vbax3XAO/R+eB+WRfbOAVItP6Pxrzu7vNafilwMXA/PRzcW3F64vAe3E/Ffk1Pwrw2Znexju62186u9iHrfmqBRdZjehJ42mv5U5iZI6OtWKuBqT28bZfbA3twH9ZFXL+yHncIcDpmIH1mT16LwE+BzzH7fwTwOPB8F9eNwgzwF3m9p7RgBhYKM5huAIb1cL8/5j5g7Xcvd7jsTeDyHu4va4GVwCDMe8SrmF8DPMvvseJLACZb+81Kr/3kB9bzWMi3B/c93m7W9c/FfB5Nte5vJ/Z78WDMoPkKzOfNtdZzmHCM9f0TkwiKBC7DfF4N7ck+3sfrOuHXQF+uKxBOQCJeyblgODkegJx86wTEYb7NX9bN9R6l88H9HKDqOO7vl5js4ThMJuE3HZbvpAeDe+v8FuCvHW7775hscah12fcwH4j52IP7dcBbHdadg9fgC/hvTAZGWXH+tMP1QzAD5F9b59fz7cH9zzED23jrspMZ3D8M/OIYy6Ot9Qz1umwD8HiH6z2AGZh3uZ29rrvG2sa/xGtw77U8jA6D++7isJ6PVzpsx0ZgqXV+L3CF1/K1wGGv898Aqzo8Tw92Etsuvj24vxV4/xiPt8t1dxe312spEziD4xjcd/I8PgZs72TZb+nB4N56Xj7DTOZyrMHqRMxg7Moulne3vbrch4DfAU96nR+HGdzGYAY+LcApXssfxx7AdXnb7rYH3QzuO4nzJeBHXuePtb3u5+hExPl0+LXCa9n1QDb2L5+nY2p/va9Tiv2r5DH3++72AcyXh0aOznD/G/BIT7ZDJ+u7FPjS63wBsMLr/G/oZICI1/vsiWw3a/mTwO+8zi8Fiqz/LwD2d7h+JnBTF+s6BfNrh/evEO9j/ZrY3T7eh+vqlddAb69LTv57kj73oqP5mJ9Nt53g7TMBl1Jqk1LqPKVUQg9uMxaTVXxQa333idypUmoSsBCTPfV2BPOT9wrr/DrMYMnbqZgMf1frTsEcZPMZZgA0GnjW+zra9NV9DnNQTlf2YgZAPz7GdXrqmDFjvoR4//X8P639jFLzMF/Guq3DVUrFAr/G/KpzPLqLQ3WyrLvlydYBUAD/B6xRSg1SSo3CZEJf62FsZwA5SqlXlVJlSqmdSqlTvZYfa93dxQ3mQ/R+TFazM5uVUqVKqTeUUtO9FyilFiilqjEZtsswv3ycqH8B3tNaf9HZQqXU35RSnrKCQswvFp3pbnsB3GMt+0ApdY7X5VMx2VoAtNZZWAMQ6+TSWmd6Xf9z6zbd3bancpVS+UqpR6zuF9+ilIoE5gL7Oyx6TylVpJR6XimV5nX534GzlFJJ1oF712Ay3J25HnNMj7bO7wW+UUpdqJQ
      "text/plain": [
       "<Figure size 864x1008 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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( forecastModel + \" \" \\\n",
    "          + ob.getLocationName().decode('UTF-8') \\\n",
    "          + \"(\"+ str(ob.getGeometry()) + \")\" \\\n",
    "          + \", \" + str(ob.getDataTime())\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=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 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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}