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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAJ4CAYAAAA5uar5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd4m9XZ/z/He8R7Jh5ZZBL2DoGwXl7GW6CFl1lCoAVKoMyWPQKkBCiBFiijzEAaxo/RQkvhZYVCWE0TZkjixLa8l7xt2Zal8/vjkYLseMjWIz3D53NdvmxLj85zJPsr3c997vt7hJQShUKhUCgUCj2IMnoCCoVCoVAo7IMKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLBQKhUKhUOiGCiwUCoVCoVDohgosFAqFQqFQ6IYKLHwIIY4VQvw1iOP2FEJ8Gok5+c63VAjxSaTOp1DYASHExUKIPxg9j0gihJBCiN1CHCNPCPGDECJer3kpJh62DyyEEMVCiM6ALymE6Ar4/TDfoXcBdwc8bvBxrQBSym+AViHET8I4l7AghFguhFgT8HuBEGKLEOJBobFOCPHLIR43zTdX/zzrhRCPCCFiA44pF0K4Bj2/KcPM40YhRJnvmCohxEther6ZQojXfa+xQwhx9qD7c4QQa4UQrUKIFiHEX0YY6zIhxAYhRK8Q4tkh7j/d94bcIYTYLIQ4ZYSxjhRCfCiEaBNClA+6L1cI8YIQosZ3/3ohxEGjPM+zfc+vSwjxVyFEZrCvQZjHmiWEeFEI0SiEaBdClAghHhJCFPruP8L3f/XaoMft5bt93Xg0I4SIA24Gfh9w21FCiI2+eZQKIS4aZs7PiEEf0IPO3ymE8AghHhrheV8lhKjz/f2eFhb6kJZS1gMfAkO+PgpFMNg+sJBSVkgpJ/m/fDfvFXDbx0KIA4A0KeXngx4eeFx6wO1/AS4Ox1zG/gzHhxBiKvAv4A0p5eVSShnEw9J9894DOAS4dND9Pwl8flLKmiHOex5wLnCMb6z9gfdDejLD8yegD8gDzgEeFULsHnD/a0AdMBXIBe4bYawaYAXw9OA7hBAFwBrgaiAV+C2wVgiRO8xYXb5xfjvEfZOAfwP7AZnAauAfQohJQxyL7/k8jvaa5gHdwCMBh4z2GoRrrN2AL9Bet32klKnAocAOYFHAoY3AQiFEVsBt5wHbYNyaORnYIqWs9s0lFnjd99zSgDOA+4UQew2a8yJg5uDBBp0/D3AB/2+Y5/3fwPXA0cA0YAZw+1DHmphxvb8pFDuRUk6oL0ACuw267VbgydGOC7ivAO3NJT4Mc8kC3gDagS+BO4FPdHruy9E+AGcCDuDOQfevA345xOOm+eYaE3DbvcCfA34vRwsWRpvDw8AfRrg/DXgKqAWq0T7Mo333LQXWAw8BbcAW4OhhxklG+xCcHXDb88Ddvp+P9c05eoyv4Qrg2UG3HQQ0DLqtEThklLGOAcqDOGc7sN8w990FrA34fabveaeM9hqEeaw1wJujPK8jgCrgMeBS323RvttuBdYFo5khjnkauDng9zzf45ICbvs3cFbA7zHAJmDPkc6BFvSUAmKY+9cCdwX8fjRQN2j+y4ASoANN3zOBz3x/55eBuBGe2wXAD0AL8A4wdajXBohHC5QrgHrfa5zou+8H4H8GPfcmYN+A37sDx1Zf6mssX7bPWATJHsDWYA+W2pWQG5gThrn8CegBJqO9iVyg8/gz0DIVj0spbxnPAL4ljv8GBmd4guFzYIkQ4rdCiP2FENGD7l8N9AO7AfugBQCByzMHob2xZwO3Aa8FpusDmA14pJTbAm77GvBfYR+M9jdfLYRwCiH+LYRYPI7nA7AB+EEIcZIQItq3DNILfDPO8XYihNgbiAO2+34v9i3dFPsO2R3teQEgpdyBLwBglNdAz7GG4Bjg1SCf5nPAEt/P/w18j5bpGC8D9Cy19P4LwPm+v88haFmqwNqlq4B/SW2pcyTOA56TUg6X4RvwGvp+zhuUkTkOLSN1MHAt8Ge0DFARsAA4a6iBff9XNwI/A3KAj33PayjuQfub7Y2mpQK0YA3fYwLP8d9Ak5RyI4CUsh/t/21ARkehCBYVWGiko109DGaj7423VQjx4KD7OnyP0w3fh+ypwK1Syi4p5XdoH7R6sgDt6nM8dQ1NQqs1qUZL578y6P6/BrxeQxbCSinXAL9GezP7CGgQQlwPWuEYcDxwpe/5NwAPAGcGDNGAlvFwSylfQvsAOXGIU01Cy2oE0oZ29Q1QiBa0fAjkA6uAvwkhskd7EYZ4Th60D8e1aAHFWuBiKWXXWMcKRAiRipYVuF1K2eY7V4WUMl1KWeE7bKTnOeJroOdYQ5CNtszkfy6X+f4vOoUQTwQeKKX8FMgUQsxBCzCeG2bMYBlKzy+gfbD2on0g3ySlrPTNrQgt9X8rI+ALwBYzsiYHv07+nwNfp3uklO1Syu+B74D/k1KW+v7G/0QLqIfiYmCllPIH34f/XcDevmXNwHkK4ELgKills5Syw3esX0drgZOEEEm+38/23RaI7u9viomDCiw0Whj6DXJf3xtvupTy8kH3pQCtOs8jBy0NWRlwm0Pnc7yBlir+YPAbUhBkS63WJAltSeLtQfefEvB6DVu8KKX8i5TyGLQ3rl8Bd/jWpqcCsUCtP0BBWxcPrFWoHnS16ACGKhLtRKt3CCSVHz9wXGjLEE/5gpQX0V73Q4d99sMghDgGbWnoCLTswmLgSV+2YVwIIRKBN4HPpZQrRzh0pOc52msQzrGcaFk3AKSUD/v+d/6A9jcezPPAZcCRaPUQoTBAz0KIuWiB9BK0v8/uwLVCCH9A+gfgDn/wNgJL0JYly0Y4ZvDr5P858HWqD/jZNcTvQ9bToOnjjwHaaAYEWjYikBw0jf4n4Ni3fbcjpdyOthzyE19wcRK7BhbheH9TTBBUYKHxDVraMCh8SwFxjGH5JEga0ZYBigJuKx7m2HEjpbwa+DtacDH4TSmYx7uAZ4FDxnOFHzCOW0r5/9Be/wVoH+y9+AIY31eqlDIw3V7guyLzU8zQafNtQIwQYlbAbXuhpdnxnTOYgtVg2Bstjb5BSumVUv4brXDxmPEM5usi+CtaZmi0IrrvCUhZCyFmoK2vb2P01yCcY72PlrIPlufRag/eklJ2j+FxQzFYzwuArVLKd3x/n63AP9CyY6DVQfze18nhz7J8JnbtelnC6BnEAa+h7+d6KaVzPE9kEJVombD0gK9EX8YnkCa0AGX3gOPS5I/Fr/DjcsjJwGZfsAGAECIGbfkkcElHoQgaFVhovIV2lRksRwAfSCl79ZyEL6X+GrBcCJEkhJiPtqYbDi4DPgDe9y1B+IkRQiQEfO1yden74DsXLdU9pjdMoflynCiESBFCRAkhjke7gvxCSlkL/B+wSgiR6rt/5qDah1zgciFErBDif4F5aH+/AfiWIV5Dy4YkCyEORXsTfd53yOtAhhDiPN+6+2loV37rh5l3jBAiAa24MNr32sT47v43cJg/QyGE2Ac4jGFqLHzPKwHtyl34xorz3ReLtsTkApZIKb2jvKR/QbvyPEwIkQzcAbwmpewI4jUI51jLfa/J/f7g1ReEzhvqYF8WYDFw0yjPNxgG63kTMEtoLadCCDET+B9+/OCcjRYA7O37AvgJAZkTIcRCtP+PIbtBAngO+IUQYr4QIgOt7fXZ0J7OTh4DbhC+ThwhRJpPAwPw/c88ATwgfJ1JQmst/++Aw15EWwq8hF2zFQeiZfP0zpYqJgpGVo4a8cUwFd9oHw4HjXac775/ACeFYy5o6cq/E8aukIDfo9DeCL9FWxNf55tT4NcafuwK6fR9taLVRxwQMFY5wXWF/Aztw7vF9xy/BZYG3J8GPIrWGdC
      "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(\"149.165.156.89\")\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 parameter\n",
    "#[b'staName', \n",
    "# b'numSigT',\n",
    "# b'numSigW', \n",
    "# b'numMand', \n",
    "# b'numTrop', \n",
    "# b'wmoStaNum', \n",
    "# b'validTime', \n",
    "# b'staElev', \n",
    "# b'rptType', \n",
    "# b'numMwnd']\n",
    "\n",
    "\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.getString(b\"prMan\"))\n",
    "        tpMan = np.append(tpMan,ob.getString(b\"tpMan\"))\n",
    "        tdMan = np.append(tdMan,ob.getString(b\"tdMan\"))\n",
    "        wdMan = np.append(wdMan,ob.getString(b\"wdMan\"))\n",
    "        wsMan = np.append(wsMan,ob.getString(b\"wsMan\"))\n",
    "        continue\n",
    "    if set(parm_array) & SIGT_PARAMS:\n",
    "        sigtGeos.append(ob)\n",
    "        prSig = np.append(prSig,ob.getString(b\"prSigT\"))\n",
    "        tpSig = np.append(tpSig,ob.getString(b\"tpSigT\"))\n",
    "        tdSig = np.append(tdSig,ob.getString(b\"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",
    "#u,v = wind_components(spd, direc)\n",
    "\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 +=  \" \" + ob.getString(b\"staName\").decode('UTF-8')\n",
    "title_string += \" \" + str(ob.getDataTime().getRefTime())\n",
    "title_string += \" (\" + str(ob.getString(b\"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
}