python-awips/examples/notebooks/Forecast_Model_Vertical_Sounding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ModelSounding class allows us to create a vertical sounding through any available AWIPS model with isobaric levels.\n",
    "\n",
    "* A Shapely Point geometry is used to select longitude and latitude:\n",
    "        from shapely.geometry import Point\n",
    "        point = Point(-104.67,39.87)\n",
    "* Parameters `['T','DpT','uW','vW']` are requested for all isobaric levels available for the selected model.\n",
    "\n",
    "* There is a single-record query performed for `level = \"0.0FHAG\"` to determine the surface pressure level.\n",
    "* Pay attention to units when switching models.  This notebook was written for the NAM 40km AWIPS model where temperature and dewpoint are returned as Kelvin and wind components as m/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using NAM40 forecast time 2018-10-14 12:00:00\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from awips.dataaccess import DataAccessLayer, ModelSounding\n",
    "from awips import ThriftClient\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from metpy.plots import SkewT, Hodograph\n",
    "from metpy.units import units\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "from math import sqrt\n",
    "from datetime import datetime, timedelta\n",
    "from shapely.geometry import Point, Polygon\n",
    "import shapely.wkb\n",
    "import timeit\n",
    "model=\"NAM40\"\n",
    "parms = ['T','DpT','uW','vW']\n",
    "server = 'edex-cloud.unidata.ucar.edu'\n",
    "DataAccessLayer.changeEDEXHost(server)\n",
    "\n",
    "# note the order is LON,lat and not lat,LON\n",
    "point = Point(-104.67,39.87)\n",
    "\n",
    "inc = 0.005\n",
    "bbox=[point.y-inc, point.y+inc, point.x-inc, point.x+inc]\n",
    "polygon = Polygon([(bbox[0],bbox[2]),(bbox[0],bbox[3]), \n",
    "                   (bbox[1],bbox[3]),(bbox[1],bbox[2]),\n",
    "                   (bbox[0],bbox[2])])\n",
    "\n",
    "# Get latest forecast cycle run\n",
    "timeReq = DataAccessLayer.newDataRequest(\"grid\")\n",
    "timeReq.setLocationNames(model)\n",
    "cycles = DataAccessLayer.getAvailableTimes(timeReq, True)\n",
    "times = DataAccessLayer.getAvailableTimes(timeReq)\n",
    "fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)\n",
    "\n",
    "print(\"Using \" + model + \" forecast time \" + str(fcstRun[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found surface record at 832.6MB\n",
      "Using 32 levels between 832.6 and 50.0MB\n"
     ]
    }
   ],
   "source": [
    "p,t,d,u,v = [],[],[],[],[]\n",
    "use_parms = ['T','DpT','uW','vW','P']\n",
    "use_level = \"0.0FHAG\"\n",
    "sndObject = ModelSounding.getSounding(model, use_parms, \n",
    "                                      [use_level], point, timeRange=[fcstRun[0]])\n",
    "if len(sndObject) > 0:\n",
    "    for time in sndObject._dataDict:\n",
    "        p.append(float(sndObject._dataDict[time][use_level]['P']))\n",
    "        t.append(float(sndObject._dataDict[time][use_level]['T']))\n",
    "        d.append(float(sndObject._dataDict[time][use_level]['DpT']))\n",
    "        u.append(float(sndObject._dataDict[time][use_level]['uW']))\n",
    "        v.append(float(sndObject._dataDict[time][use_level]['vW']))\n",
    "    print(\"Found surface record at \" + \"%.1f\" % p[0] + \"MB\")\n",
    "else:\n",
    "    raise ValueError(\"sndObject returned empty for query [\" \n",
    "                    + ', '.join(str(x) for x in (model, use_parms, point, use_level)) +\"]\")\n",
    "    \n",
    "\n",
    "# Get isobaric levels with our requested parameters\n",
    "levelReq = DataAccessLayer.newDataRequest(\"grid\", envelope=point)\n",
    "levelReq.setLocationNames(model)\n",
    "levelReq.setParameters('T','DpT','uW','vW')\n",
    "availableLevels = DataAccessLayer.getAvailableLevels(levelReq)\n",
    "\n",
    "# Clean levels list of unit string (MB, FHAG, etc.)\n",
    "levels = []\n",
    "for lvl in availableLevels:\n",
    "    name=str(lvl)\n",
    "    if 'MB' in name and '_' not in name:\n",
    "        # If this level is above (less than in mb) our 0.0FHAG record\n",
    "        if float(name.replace('MB','')) < p[0]:\n",
    "            levels.append(lvl)\n",
    "\n",
    "# Get Sounding\n",
    "sndObject = ModelSounding.getSounding(model, parms, levels, point, \n",
    "                        timeRange=[fcstRun[0]])\n",
    "\n",
    "if not len(sndObject) > 0:\n",
    "    raise ValueError(\"sndObject returned empty for query [\" \n",
    "                    + ', '.join(str(x) for x in (model, parms, point, levels)) +\"]\")\n",
    "    \n",
    "for time in sndObject._dataDict:\n",
    "    for lvl in sndObject._dataDict[time].levels():\n",
    "        for parm in sndObject._dataDict[time][lvl].parameters():\n",
    "            if parm == \"T\":\n",
    "                t.append(float(sndObject._dataDict[time][lvl][parm]))\n",
    "            elif parm == \"DpT\":\n",
    "                d.append(float(sndObject._dataDict[time][lvl][parm]))\n",
    "            elif parm == 'uW':\n",
    "                u.append(float(sndObject._dataDict[time][lvl][parm]))\n",
    "            elif parm == 'vW':\n",
    "                v.append(float(sndObject._dataDict[time][lvl][parm]))\n",
    "            else:\n",
    "                print(\"WHAT IS THIS\")\n",
    "                print(sndObject._dataDict[time][lvl][parm])\n",
    "        # Pressure is our requested level rather than a returned parameter\n",
    "        p.append(float(lvl.replace('MB','')))\n",
    "\n",
    "# convert to numpy.array()\n",
    "p = np.array(p, dtype=float)\n",
    "t = (np.array(t, dtype=float) - 273.15) * units.degC\n",
    "d = (np.array(d, dtype=float) - 273.15) * units.degC\n",
    "u = (np.array(u, dtype=float) * units('m/s')).to('knots')\n",
    "v = (np.array(v, dtype=float) * units('m/s')).to('knots')\n",
    "w = np.sqrt(u**2 + v**2)\n",
    "\n",
    "print(\"Using \" + str(len(levels)) + \" levels between \" + \n",
    "      str(\"%.1f\" % max(p)) + \" and \" + str(\"%.1f\" % min(p)) + \"MB\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Skew-T/Log-P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAM2CAYAAAAuNu86AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8VNXd/98nE0JCEgh72HeDhJCE7IgoiwRZBEWRRQE3iq2t2j72eVqr7a+rfV59WlsrKlWLyi67LKKCFUEEEgi7YQ9LwpKEkD2TzJzfH3cGxzATJpLMvTdz3q/XvDK599xzvnPnzvd+7znnc75CSolCoVAoFAqFQqFofAL0NkChUCgUCoVCofAXVPCtUCgUCoVCoVD4CBV8KxQKhUKhUCgUPkIF3wqFQqFQKBQKhY9QwbdCoVAoFAqFQuEjVPCtUCgUCoVCoVD4CBV8KxRNHCHEaCHEGr3tqA9CiD8JIZ7T246miBCiuRDiGyFEB5dt9wkhluppl0KhUPgLKvhWKLxACHFGCHFJCBHqsu1JIcR/apUTQohTQogjbur4jxBCCiFia21f49h+t5tjtjr2Bbps6ymE+FwIUe4IokbdxPw/Aq+4HC+FEGVCiFIhxAUhxF+FEBaX/bOFEAcd9V8UQrwhhIhw2f8bIcTCWvUdFEIEuGz7vRBigRDiTkc7pY42pcv/pUKI7m4+c3tgJvCWpw8khHhGCJEhhKgSQixws3+k49yUO85VDzdl2gghrgghttd18oQQvYUQ64UQJUKIfCHE/7rsK631sgkhXvNQz3DHeSoSQhQIIVYLIbq47O8ihFgrhCgUQpwXQsy9iV0/FkKcFkIUO87FUJd9m2rZZRVCHASQUlYB7wL/7SwvpVwHDBRCDPLQVnMhxDtCiBzHedgnhLi3VhmP51wIMUUI8ZVj33/c1D9CCLHX8VlOCSHm3OSz/85xLmuEEL9xs3+6w9Yyx++rTR11DRRCbHZ8tx4TXwgh+gkhKl2v/frU5c05dFNfXee0uRDiXcc5uyiE+KkP63reUe6a47jmLvvq9E91HatQ+Asq+FYovCcQePYmZYYBHYDeQogkN/uPoQWWAAgh2gKpwJXaBYUQMxxt1mYJsA9oC7wIrBBawHoDDhtaSSm/rrUrVkoZBowEpgNPOcr/DPgz8ALQymFbD+BTIUSQh88M0BmYWnujlPJLKWWYo61ox+YI5zYp5Vk3dc0GNkopK+poLxf4PVoQ+R2EEO2AVcBLQBsgA1jmpo4/A0fraAPHZ/4U2ApEAl2B68GXy+cIAzoCFcCHHqo7AqRLKSPQztdx4A2X/QuB0456xgF/FEIM92BXCtoD1YNo39M7wGrheIiSUt5by7avatm1GJhVK/BZAngKegOBc8BdjvZeApYLIXo67LnZOS8EXsXlIdDlszQDVqM9bLUCHgb+Kmo9pNbiBPBzYIOb+qIddT2Kdi7LgXl11FUNLAeeqKMMwOvAnpuUqauuOs9hbbw4p78B+qH9PocDPxdCjPFBXenA/6D5jp5Ab+D/uRTx6J+8OFah8A+klOqlXup1kxdwBu2mUYgWPAI8CfynVrl3gUVoN7p/1tr3H+Bl4DxgcWx7Bi0AOw/c7VK2FVqgngpIINCx/TagCgh3KfslMNeD3S8Db9faJoG+Lv9/CPwTaAmUAlNqlQ8DLgOPO/7/DbCwVn3/jRZMOu38PbCgVj09XT9LHed6K/CIl9+Lu3bmAF+5/B+KFhT3d9mWBuwEHgO211H/HOBLL22ZBZwChBdlmwN/Ao64nGMJtHcpMx/4wMPxDwO7a31GCXRyU7YnYAN61dp+HLjL5f87gNP1+E0cACZ7e87r+M10dNjewmXbHmCaFzYsBH5Ta9sfgcUu//cBrK6/GQ919QWkh31T0YLq71z736cuT+fQw7Xn8ZwCF4DRLvt/Byz1QV2LgT+6/D8SuOh4X6d/qutY9VIvf3qpnm+Fwnsy0ALo/3K3UwjRAq0ncpHjNdVNb3EuWg/oaMf/M4H33VT3R7Sg/GKt7dHAKSllicu2/Xzbq1ybGCDbwz6EEAOAO9F6qoYAwWgPDteRUpYCm4B7PNXjOKYYrdf6VqnTZi+IRjsnAEgpy4CTju04eodfR3vw8TjNwEEqcMYxjSNfaFOHYjyUnQW8L6Wsa+pCdyFEEVrg81+AcwqLqPXX+X6gh6o2ARYhRIrj8zwOZHHj9QLaNfallPJ0re1Hgdha//cUQrT0ZL/L5+iIFmgddmyq85zXhZTyElpv6WNCCIsQIg2tB7bO6UB1UNuWk2jB923fpzLH+fgt8LPvaY+nemufQ4Q2Jck5fcjjORVCtEYbPdn/bY3f9QMNWVctot2U7egYxbuZf6rrWIXCb1DBt0JRP14GfuxhmscDaL0+nwDr0YaZx7kp9z4wUwgRhdaLvtN1pxAiEa0X0t3c4TDgWq1t14BwD/ZGACVutu8VQlwFPgLeBv4NtAPypZQ1bsrnOfZ7QqINab/cAHM4PdnsLTc7Rz8BdkkpM72oqytar+c/0AKUDcDa2g9VQpu7fhfwXl2VSSnPSm3aSTvgV8A3ju0lwA7gJSFEsBBiMDAZaOGhqhJgJVqAWgX8GpjjIfCfCSzwUEdErf+pte0GHNNEFgHvSSm/cWyu73VZmyVov60qtJ7SF6WU57w8tja3akttfge8cwv23ICHc4iUMkJK6XzoqOtzhLn8X3tfg9dVi9p1Od+Hu9lXu666jlUo/AYVfCsU9UBKeQgtsP4fN7tnAcullDVSE7WtcmyrzSpgBPBj4APXHUITLc4DnvUQBJeiTQ9xpSWeg9WruL+xDZZStpZS9pFS/kpKaQfygXbCRdzpQifHfo9IKTcCZ/E8b9hbvmOz+K54cIYXx3s8R0KIzmjB94te2lKBNi1lk5TSCvwFbS7r7bXKzXSUq9277BYpZSFaoL7W5XzPAHqhzQt+Ay04O++hiifRerujgSDgEWC94/Ndx9HzGQmscFNHOFBU639qbfsOjuvzA7Se5GdcdtX3unStsz/a/OOZjs8SjTbneJxj/2GX7//Om9VXly1CiBkudW3ywrY4YBTwNy/a9Yo6zmFt6jqnpS7/197n67qc70vc7KtdV13HKhR+gwq+FYr682s0gaLrShVd0QLqRxxK/otoU1DGOsRO15FSlqNNG3iaWsE32s0oEVjmqMMp8DrvCDwOo4k5XQPqWFyGrmtxAO+H23ei9Tw+4LpRaCu83Ats8aKOX6EFtp56bL3hOzbL74oHF3lx/GFcplM47O/j2J6M9iBxxHF+/w4kO74zi5u6DnDzqSmgBY519nq7IRBNnNsSQEqZI6UcL6VsL6VMQQvyd3s4Nhb4SEp5TEppl1J+jDY6MaRWuVnAKsfUodrcznenANwOnJFSFrtrUAgh0ISdHdHmKVe77K7rnN+MgUC2lHKz47Nko40w3AsgpYx2+f6/9KK+2rb0Rptjf0xKucilrjpXGnFwN9qc+bOO6+W/gMlCiL1eHHsDNzmHN/sc18+plPIq2vftOm2oLj/QaHU53l+SUhZwc/9U17EKhUeEEL30tqFB0XvSuXqplxleaILLUS7//wsowCEeA36BNmc2stbrFPBjR5n/AE863ncGhrrUdx7tRi9qHZ+EFvx1AYIcZb9G64ENBu5H66ls78HuwWhBh+u27wgua+37OXAJGAM0Qws8NgJ7geaOMr/hRsGlq4DzU8e5WVCr7p54J7j8KTD/JmUCHZ//T2gPMMF8K/ZsjzacPdmx/c/A1459zWud32eBXUCkh3ai0FbLGAVYgOfR5soGuZQZApRxc0HfA476Ahw2Lgf2uuy/Ha332dmTnV/H9zoLTZDb23HN3OOw01VUGuK4Nka4Ob6L4ztq7rLtl8C8Oux/03HthbnZ5/GcO/ZbHNvnAtsc75s59vVB6xEd4fgsfdBWM3mqDluaOepYjCa6DeZbEXM0mv7gTjRh4UI8iAcd5YXj+AGO6zOYb6/1FrWul7+gjSJ4+l481nWzc/g9zukrwBdAa6A/WgA9xgd1jUHTFgxwlN8KvOKy36N/utm
      "text/plain": [
       "<Figure size 864x1008 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = (12, 14)\n",
    "\n",
    "# Skew-T\n",
    "skew = SkewT(rotation=45)\n",
    "skew.plot(p, t, 'r', linewidth=2)\n",
    "skew.plot(p, d, 'g', linewidth=2)\n",
    "skew.plot_barbs(p, u, v)\n",
    "skew.plot_dry_adiabats()\n",
    "skew.plot_moist_adiabats()\n",
    "skew.plot_mixing_lines(linestyle=':')\n",
    "\n",
    "skew.ax.set_ylim(1000, np.min(p))\n",
    "skew.ax.set_xlim(-50, 40)\n",
    "\n",
    "# Title\n",
    "plt.title( model + \" (\" + str(point) + \") \" + str(time.getRefTime()))\n",
    "\n",
    "# Hodograph\n",
    "ax_hod = inset_axes(skew.ax, '40%', '40%', loc=2)\n",
    "h = Hodograph(ax_hod, component_range=max(w.magnitude))\n",
    "h.add_grid(increment=20)\n",
    "h.plot_colormapped(u, v, w)\n",
    "\n",
    "# Dotted line at 0C isotherm\n",
    "l = skew.ax.axvline(0, color='c', linestyle='-', linewidth=1)\n",
    "\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
}