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": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAJ4CAYAAADWe8HHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXd8W+XZ//++vVe84hmP7ISwAqUJu4RRRhnhW1paSKFQHjaFsFoeZihlP4wSNiWUEaAtv5aRDnaAMkPZM05sy5a35SXbsrXu3x9HCsI4tiwd6Yzc79dLL8vS0X0ujY90neu+r88RUkoUCoVCoVAo9CLF6AAUCoVCoVDYC5VcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVcKBQKhUKh0BWVXCgUCoVCodAVlVwoFAqFQqHQFZVchBBCHCyEeDqK7XYWQryVjJhC+ztJCPGfZO1PobAyQojThRC3Gx1HMhFCSCHEvDjHKBdCfCmEyNQrLsW2je2TCyFErRBiMOIihRBDEf/vG9r0OuCGiMeN3a4PQEr5CdAnhDgygbEkBCHEKiHEYxH/VwkhvhJC3CE01gsh/mecx80KxRqOs0MIcbcQIj1im0YhhGfM85uxlTgKhRD3CCHahRDDQohPhRAnR/kcJk22hBCZQog1QoiB0D4umOC5DAohrphgrHOEEO8LIUaFEH8a5/5jQ1/KbiHEF0KIoycYa38hxKtCiH4hROOY+8qEEE8IIVpD978phNh9kud5vBDCEfoMPS2EKI64r1gI8ffQfQ4hxPFJHGu+EOJJIURX6D2oE0KsFkJUh+5fFnoP/jbmcYtDt6+PRStCiAzgcuDmiNsOEEJ8EIqjXghx2lZifkiM+ZEes/9BIURACLF6gud9fujz1h/6/Fnmh1pK2QG8Coz7+igUU8X2yYWUsklKmRe+hG5eHHHbG0KIJUCBlPKdMQ+P3K4w4va1wOmJiGXqzzA2hBAzgdeBZ6WU50opZRQPKwzFvROwJ3D2mPuPjHx+UsrWcfabAbwEzAyNUQBcDNwQmQTEySpgfmgf+wO/EUIcOt5zCV2umWCsVuD3wJqxdwghqoDHgAuAfLTn8bgQomwrYw2Fxrl4nPvygA3AbkAx8DDwDyFE3jjbIoTYAbgPOAEoB4aBuyM2uQvwhu5bAdwTekyix5oHvIv2uu0qpcwH9gY2A/tEbNoF7CWEmB5x2y+BjRCzVpYDX0kpW0KxpAN/Dz23AuBnwK1CiMVjYt4HmDt2sDH7Lwc8wF+38rwPAS4BDgRmAXOAq8fb1sTE9L2mUIyLlHKbugASmDfmtiuBP062XcR9VWhfNJkJiGU68CwwALwHXAP8R6fnvgrtx3Au4ACuGXP/euB/xnncrFCsaRG33QTcH/F/I3BQFDGcAnQCuWNu/xkwCOSH/q8B/ob2I+QC7gQWASNAILRt31b20QIcHPH/NcCTW3suUb52vwf+NOa23YHOMbd1AXtOMtZBQGMU+xwAdtvKfdcBj0f8PxctAZgG5IauL4i4/1HghiSM9Rjw3CTPaxngBO4Fzg7dlhq67UpgfTRaGWebNcDlEf+Xhx6XE3HbBuC4iP/TgA+BnSfaB1riUw+Irdz/OHBdxP8HAu1j4j8LqAPcoc/kXODt0Pv8FyBjguf2K+BLoBd4Hpg53msDZAL/BzQBHaHXODt035fAEWOeezfwvYj/hyPHVhd1ifVi+8pFlOwEfB3txlI7MvIBCxMQy11oP6CVaF8ov9J5/DloFYv7pJRbnQ6YiNB0xyHA2EpPNPwQ+JeUcmjM7f8fkAXsKYRIBdahJUCz0JK5J6WUXwJnAG/L71aTwrEVATOAjyNu/hgYe6TtEEI4Q+XwkhieB8D7wJdCiKOEEKmhKZFR4JMYx9uCEGIXIAPYFPq/VgjRJ4SoDW2yAxHPUUq5mVASELoEpJQbI4bc8hroOdY4HIT2XkbDI8CJoeuHAJ+jVTxi5Vs6llqp/wng5ND7sydaNStyWu184HWpTXdOxC+BR6SUW6vwfes1DF0vH1OZORStMrUH8BvgfrRKUA2wI3DceAOHPleXAj8GSoE3Qs9rPG5Ee892AeahaefK0H1PjNnHIUC3lPIDACmlH+3z9q3KjkIRCyq50ChEO5oYywehL+E+IcQdY+5zhx6nG6Ef1WOAK6WUQ1LKz9DK43qyI9rR6J9jeGy30NaetKCV+J8ac//TEa/X1hbHlgBtY28MfbF1h+5fipYgXBx6HUaklNEuag2X0PsjbutHOwontI8laD8yu4VuXxvl2GNjDqD9QD6OllQ8Dpw+TuI0JYQQ+WjVgaullP2hfTVJKQullE2hzfL49nOEb57nRPfpOtY4lADtEc/lnNDnYVAI8UDkhlLKt4BiIcRCtCTjka2MGS3j6fgJtB/XUbQf5cuklM2h2GrQpgGuZAJCSdh+TKzFsa9T+Hrk63SjlHJASvk58BnwgpSyPvQe/wvYdStjnw5cL6X8MqST64BdQlObkXEK4FTgfCllj5TSHdr256FNHgeOEkLkhP4/PnRbJLp/rym2TVRyodHL+F+W3wt9CRdKKc8dc980oE/nOErRSpPNEbc5dN7Hs2jl41fGfjlFQUmoWpADvAn8e8z9R0e8Xltb2NiNVpX5FkKINLQfpm60IzlH6It0qgyG/uZH3JZP6EdHSjkopXxfSukPHdmeAxwc+kGfEkKIg9Cmh5ahVRn2A/4YqjrEhBAiG3gOeEdKef0Emw7y7ecI3zzPie5L9FguIt5fKeWdoc/M7UD6ONs/ivYe7I+2PiIevqVjIcR2aEn0iWjvzw5o628OD21yO/C7cAI3ASeiTU02TLDN2NcpfD3ydeqIuO4Z5/9x19egJcJ/CCfuQA8g0KoSkZSiafO/Edv+O3Q7UspNaFMjR4YSjKP4bnKRiO81xTaISi40PkErJUZFaFoggylMpURJF+BH+3ENU7uVbWNGSnkB2rTDK6FFiVN9vAf4E9oUxlSnFF4CDhNC5I65/Ri0o8t30JKr2lDC8Z3dTxJbL1plJLK0uxit5D7uQ0J/xSRxj8cuaCX196WUQSnlBrTFjAfFMBah7oKn0SpDky2s+5yI5yiEmIM2374xdEkTQsyP2H6i10DPsV5GK99Hy6NoaxH+KaUcnsLjxmOsjncEvpZSPh96f74G/gEcFrr/QODmUIdHuNry9jjdMCcyeQXxW69h6HqHlNIVyxMZQzNaRaww4pIdqvxE0o2WpOwQsV2B/GZBLHwzNbIc+CKUcABbEvx5fHt6R6GICZVcaPwT7agzWpYBr0gpR/UMIlRm/xuwSgiRI4TYHm2uNxGcA7wCvCyEKI+4PU0IkRVx+c7RZuhH8AS08vdUvzwfRVu491ehtYWmh1ba3wGsCh1FvoeWINwghMgNxbF36PEdQHWo62RrPAJcLoQoCh29noqWDCGE2F0IsVAIkRKaD78DbQHhuEevQog0IUQW2oLD1FAs4aRnA7BvuFIhhNgV2JetrLkI7TML7QhehMbKCN2XjjbN5AFOlFIGJ3wVtamcI4UQ+4YStd8Bf5NSukPTMn8Dfhd6/fZG+zF5NAljrQq9JreGE9dQArpovI1D1YD9gMsmeb7RMFbHHwLzhdaOKoQQc4Ej+ObHcwFaErBL6AJwJBEVFCHEXmgVgnG7RCJ4BDhFCLF9aN3P5YQ+czpwL/C/ItShI4QoEEL8dOxGoc/MA8BtItSxJLR280MiNnsSOBg4k+9WLZaiLTTWu1qq2BYxekVpsi9sZUU42g/F7pNtF7rvH8BRiYgFrYS5jgR2i0T8n4L2pfgp2pTE+lBMkZfH+KbDYjB06QNeA5ZEjNVIFN0ioW2L0doDO9B+TD9nTJcKWsXmabTkpRu4I3R7Ruj170FbjDbe+Jl
      "text/plain": [
       "<Figure size 576x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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 wind_speed, 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",
    "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",
    "locations = DataAccessLayer.getAvailableLocationNames(request)\n",
    "locations.sort()\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",
    "    parm_array = [x.decode('utf-8') for x in ob.getParameters()]\n",
    "    if set(parm_array) & 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(parm_array) & 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 = 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(-60, 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",
    "lcl_pressure, lcl_temperature = lcl(p[0], T[0], Td[0])\n",
    "skew.plot(lcl_pressure, lcl_temperature, '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 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
}