python-awips/examples/notebooks/Upper_Air_BUFR_Soundings.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following script takes you through the steps of retrieving an Upper Air vertical profile from an AWIPS EDEX server and plotting a Skew-T/Log-P chart with Matplotlib and MetPy.\n",
    "\n",
    "The **bufrua** plugin returns separate objects for parameters at **mandatory levels** and at **significant temperature levels**.  For the Skew-T/Log-P plot, significant temperature levels are used to plot the pressure, temperature, and dewpoint lines, while mandatory levels are used to plot the wind profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+YAAANsCAYAAADWS7KZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XecVNX5x/HPw9IREKQIUiwQOyKoqEnskTR7NzbUaKyJ\nxpLEGGNM7N3oT7Gghtg1tqjR2KNiAaSpqHQp0vsubDm/P86d3WGc7XPvnZn7fb9e89o7M3fueZ49\ns+WZe8655pxDREREREREROLRIu4ARERERERERJJMhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiI\niIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRI\nhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iI\niIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhI\njFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmI\niIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiI\niMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmIiIiIiIhIjFSYi4iIiIiIiMRIhbmISEjM7E0zq6rn\nVmlmf0p7TQ8zW2Jmh6Y9dkUdrz072KeTmS00s6NzGP/JDYi/ysym13OcB81sRq7iqqOdKjO7Pcvj\nVwbP3Rnc7x/cv7Ce483MyLPczOaY2d1m1i1j33r7qJ62WprZWWb2vpmtMLO1ZvaZmV1jZl0b+70I\njnmcmf26gfu2MLPfmdkbZrbYzNab2Swzu93MutTz2u3MbF2Q75AGtvdXM3vBzL4JXvdAPfuPMLMP\nzWx18L0Zb2aHNKCdn5nZQ2Y2Mcipspb9djGzB8xsqpmVmdkqM3vHzA5qSD5pxzkg6MM1ZrbIzEaZ\nWfcs+7UM3jMzgvY+N7Nz87itQWZ2v5l9HXz/15rZl0Gb+2bsm/pZqDSzzbMcq72ZrUzvd2vC78o6\nYn3XzG7NeGwbM3sm+D6tD953d5vZJvUca3TQ9vMZj+9dT6x31RdncJzzgv4oM7PpZvYnM2vZkNdG\nxczeMrM3Qjx+SfC77pKw2hApFHn1wy8iUmTOAjql3f85cBlwCjA17fFv0ravAr5yzj2bcSwHDAdW\nZjw+A8A5t9LMrgduNLNnnXPrmx8+LwK7Zzw2BngSuCntsXX1HMcFt8iZL9TPAf7mnLu8kS93wP+A\n3wIGtAKG4vtoCLBblv1r7aM6YmwHvAzsCdwD/AUoBfYALgaON7MDnHNfNTL+44HtgdsasG87/Hvz\nUeAuYD4wCJ/rvma2i3PuO/1sZi2AB4CFQO9GxPYbYALwHHBqXTua2d3AScDNwO+AqiC21g1o5zBg\nGDAe/z6t7YODY4Jj3h3s2xo4E3jOzE5yzo2uryEz2xt4CXgBuBzoAVwP/Df4/pWn7f5/wC+APwKf\n4N83t5nZRs65a/OsrTOBO4AvgFuBKfj3+rb499h/zWyAcy7zfb4KGAFckfH4Ufj/P9N/RzXld2W2\nWI8DBgOHpD3WE3gfWIT/XfANsDPwN2BX/M90tmP9LDjOiixPj+W7vxsBzgZOBJ6pK87g+Jfhf9av\nBl4LYvkb/ufoV/W9PkKh/u52zqU+cBllZg865xaG2Z5IXnPO6aabbrrpFsENOBmoBIbU8nwffPFw\nVMbjVwSv61rP8TsBq4EzQ8yhCri9ka8ZBUyP4PtbHRtQAjwcfN8uyNivf7DvhfUcbwbwfJbHLwuO\nO6CxfVRLO/cErz0yy3MDgGXARMAaedwXGvp9x4+g65Ll8Z8F36vja3ndRcAsfMFT63u7nrZXAQ/U\n8tyhQftH5OD9cQdQWctz3Wp5/CPgywYe/yNgEtAi7bE9gvjPTHtsu+B7dUmW98FqYON8aQv4PlAB\n/AtoWcs+RwCbpt1P/SzcA8zMsv87wOh6+r3O35V1xPsZcGfGY6cHx9oz4/ELgsd3ynKcTsAc4Ne1\n/R6opf2vG/IzB3QF1gJ3ZTz+++D7vU1z3++5ugFvAm+E3EaL4PfINXHnq5tucd40lF1EJH+chj+L\n9Hx9O2bjnFsJ/Bs4I5dBNYaZnWJmXwRDM6eY2YkxxNAGf8bqWOA059wtOW5idfC12WeSgrN5I4BX\nnHNPZT7vnPsauA5/5vvQjNceHwxlXhV8vyeZ2enBc2/ii+rUsP0qq2UYd9BOlXNuWZanPg2+9s0S\n+0DgSvxZwtWZz+fIr4EZzrmnQzo+AM65xbU8NYEsuWcys97ALsDDzrmqtON+AHyJP3Ofktp+MOMw\no4D2wI/zpS3gD/hC8UznXEW2HZxzTzvnFmR56gGgn5n9KC327+GL/TqnLjRFMIpgG+CRjKdS7/tV\nGY+n7mfL62b8mfXvTI2po/39gC1pWG4/AdqQvV9akPazbn4q0Coz29rM/mN+Osc8M/t98Pz3g+H7\na8xstpn9soHxtjKzP6YNpV9ofjpHt1y8Npg6MLOW148xs7Gp+8H7+DHg1GAUjkgi6c0vIpI/DgLe\nc1mGDAdaBvPxUrdsv8NfBwabWZ/wwszOzE7B/1P6Cf4fzz/ih2LvF2EYHfHDwg8EjnbOPdjM41na\n97uNmf0AXyyOcc5Ny7J/Q/oo3b74Yb3P1bHPs/ih9OkFzl/wZx2/xA8N/hF+GPamwS5nAe8BC/BD\nuXfHn1FtrB/hP4CYkuW5+/FnEv/dhOPWy8xK8HGPN7MLzM/5rzCzaWb22zDazNL+3sDkjMdT84vT\n5zvvgP8+TcpyqInB8ynbA4vcd4fsTkw7VuRtZQreu/sAn2R5fUN8BbzLhlMVTsWfRQ9jzvLB+LPQ\nH2Q8/gz+TPaNZjYg+DneA7gUeMk5t8F728wOAE4ATnfONebDt9PwHwI82IB9tw++btCHwQcci9mw\nXxx+Gs0T+GlEP8b/vvibmd2IH5nwd+Cn+JEU9wT51crMDP8B8EXAfcAB+N9r+wBvBh9uNve1o4C+\nwQcW6a/fBj8N6P6MQ78OdMNP6RFJJM0xFxHJA+YX/BlE7WdoDF9kpfsG6Jfx2ORg312pZz5mLgX/\nrP0VeNc5d0La42OAaXw39rCcjP9H9gz33Xn6TfEzoDzjsXHA4Vn2bWgfpeuHj7eueegz0vbF/IJa\nvwfucc6dlbbfu6kN59wXZrYcWOec+7iOY9fKzDYDrgE+ds69mPHcufhh0tm+D7nSDX9W8QD8nOCL\ngHn4odM3mNnGrvHrBjTGlcBW+IIvncOfZU0fgZBaRGxpluMsTXs+te939nPOrTWz9Rn7RtlWpm74\ntQdmZT4RFO2WdrzaRmM8APyfmW2MX3vhRPyc9zAMAb5IH0UQxLbC/AJ1L+I/yEr5F35UTTUz6wCM\nBG7ILNjrYmad8aMTXnXONeT37ib4n83SLM9l9iH4wvx3zrmXg/bew89/vwDY3jn3RfD4ePz79Vi+\n+wFFumPwH17+1Dn3n7Q8JuI/LDgFX/A357Uv4deeGAGkfxAzAj9l69GM46b/7fpfHbGLFC0V5iIi\n+aE7/ndybWemHLA/Gy4slm2Bt9TrG7MQVy5sHbR5TfqDzrn5ZvYuMDCiON4BdgSuMLO3ajmr3Rjv\n4hcqM/wos+/h59D+z8z2dM4tStu3oX3UXAcGsdwdwrEBML8S+0v4nI7JeG5z/IJV59cxBDwXUqMN\nOgJ7OedSZ3nfD4ZzX2hmV9dS3DRLMCXgD/gCbYMRAc65d2jYwnPNFmVbjTQW2Cl1x8wucs7dnGW/\nJ/EfNv4CX+D3oGFnlJuiN/5DwA2Y2ab43wsr8B/qzMV/qHQV8IqZDU8bpn8d/mf2qka2fQL+Q6T7\nmhZ6vSqAV1J3nHPOzL4
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f50bd1e2950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "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",
    "import numpy as np\n",
    "import math\n",
    "from metpy.calc import get_wind_speed, get_wind_components, lcl, dry_lapse, parcel_profile\n",
    "from metpy.plots import SkewT, Hodograph\n",
    "from metpy.units import units, concatenate\n",
    "\n",
    "# Set host\n",
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "\n",
    "# Set data type\n",
    "request.setDatatype(\"bufrua\")\n",
    "availableLocs = DataAccessLayer.getAvailableLocationNames(request)\n",
    "availableLocs.sort()\n",
    "\n",
    "# Set Mandatory and Significant Temperature level parameters\n",
    "MAN_PARAMS = set(['prMan', 'htMan', 'tpMan', 'tdMan', 'wdMan', 'wsMan'])\n",
    "SIGT_PARAMS = set(['prSigT', 'tpSigT', 'tdSigT'])\n",
    "request.setParameters(\"wmoStaNum\", \"validTime\", \"rptType\", \"staElev\", \"numMand\",\n",
    "                      \"numSigT\", \"numSigW\", \"numTrop\", \"numMwnd\", \"staName\")\n",
    "request.getParameters().extend(MAN_PARAMS)\n",
    "request.getParameters().extend(SIGT_PARAMS)\n",
    "\n",
    "# Set station ID (not name)\n",
    "request.setLocationNames(\"72562\") #KLBF\n",
    "\n",
    "# Get all times\n",
    "datatimes = DataAccessLayer.getAvailableTimes(request)\n",
    "\n",
    "# Get most recent record\n",
    "response = DataAccessLayer.getGeometryData(request,times=datatimes[-1].validPeriod)\n",
    "\n",
    "# Initialize data arrays\n",
    "tdMan,tpMan,prMan,wdMan,wsMan = np.array([]),np.array([]),np.array([]),np.array([]),np.array([])\n",
    "prSig,tpSig,tdSig = np.array([]),np.array([]),np.array([])\n",
    "manGeos = []\n",
    "sigtGeos = []\n",
    "\n",
    "# Build  arrays\n",
    "for ob in response:\n",
    "    if set(ob.getParameters()) & MAN_PARAMS:\n",
    "        manGeos.append(ob)\n",
    "        prMan = np.append(prMan,ob.getNumber(\"prMan\"))\n",
    "        tpMan = np.append(tpMan,ob.getNumber(\"tpMan\"))\n",
    "        tdMan = np.append(tdMan,ob.getNumber(\"tdMan\"))\n",
    "        wdMan = np.append(wdMan,ob.getNumber(\"wdMan\"))\n",
    "        wsMan = np.append(wsMan,ob.getNumber(\"wsMan\"))\n",
    "        continue\n",
    "    if set(ob.getParameters()) & SIGT_PARAMS:\n",
    "        sigtGeos.append(ob)\n",
    "        prSig = np.append(prSig,ob.getNumber(\"prSigT\"))\n",
    "        tpSig = np.append(tpSig,ob.getNumber(\"tpSigT\"))\n",
    "        tdSig = np.append(tdSig,ob.getNumber(\"tdSigT\"))\n",
    "        continue\n",
    "\n",
    "# Sort mandatory levels (but not sigT levels) because of the 1000.MB interpolation inclusion\n",
    "ps = prMan.argsort()[::-1]\n",
    "wpres = prMan[ps]\n",
    "direc = wdMan[ps]\n",
    "spd   = wsMan[ps]\n",
    "tman = tpMan[ps]\n",
    "dman = tdMan[ps]\n",
    "\n",
    "# Flag missing data\n",
    "prSig[prSig <= -9999] = np.nan\n",
    "tpSig[tpSig <= -9999] = np.nan\n",
    "tdSig[tdSig <= -9999] = np.nan\n",
    "wpres[wpres <= -9999] = np.nan\n",
    "tman[tman <= -9999] = np.nan\n",
    "dman[dman <= -9999] = np.nan\n",
    "direc[direc <= -9999] = np.nan\n",
    "spd[spd <= -9999] = np.nan\n",
    "\n",
    "# assign units\n",
    "p = (prSig/100) * units.mbar\n",
    "T = (tpSig-273.15) * units.degC\n",
    "Td = (tdSig-273.15) * units.degC\n",
    "wpres = (wpres/100) * units.mbar\n",
    "tman = tman * units.degC\n",
    "dman = dman * units.degC\n",
    "u,v = get_wind_components(spd, np.deg2rad(direc))\n",
    "\n",
    "# Create SkewT/LogP\n",
    "plt.rcParams['figure.figsize'] = (8, 10)\n",
    "skew = SkewT()\n",
    "skew.plot(p, T, 'r', linewidth=2)\n",
    "skew.plot(p, Td, 'g', linewidth=2)\n",
    "skew.plot_barbs(wpres, u, v)\n",
    "skew.ax.set_ylim(1000, 100)\n",
    "skew.ax.set_xlim(-30, 30)\n",
    "\n",
    "title_string = \" T(F)   Td   \" \n",
    "title_string +=  \" \" + str(ob.getString(\"staName\"))\n",
    "title_string += \" \" + str(ob.getDataTime().getRefTime())\n",
    "title_string += \" (\" + str(ob.getNumber(\"staElev\")) + \"m elev)\"\n",
    "title_string += \"\\n\" + str(round(T[0].to('degF').item(),1))\n",
    "title_string += \"  \" + str(round(Td[0].to('degF').item(),1))\n",
    "plt.title(title_string, loc='left')\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",
    "# Calculate full parcel profile and add to plot as black line\n",
    "prof = parcel_profile(p, T[0], Td[0]).to('degC')\n",
    "skew.plot(p, prof, 'k', linewidth=2)\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, '30%', '30%', loc=3)\n",
    "h = Hodograph(ax_hod, component_range=max(wsMan))\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
}