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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAALtCAYAAACo61k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4W+XZh+/XI96O42xnb0agjFBmaKFsSimkbCijkLTsAkkZLWV/QFmFUkbK3rSEPVoIo0DZAcKMncRJHMcZTjzieEvv98eREsXYjmSf/T73dfmyLB2d89x6ZOl33rOU1hpBEARBEAQhOKR5XYAgCIIgCIKQGhLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAU4QBEEQBCFgSIATBEEQBEEIGBLgBEEQBEEQAoYEOEEQBEEQhIAhAS6GUuoApdRzSUy3vVLqf27UFFveKUqp99xaniAEmWT/j8OEUmqJUmq/Xs4jSyn1vVJqkF11CYLgLKEPcEqpkUqphoQfrZTakPD31Nik1wHXJzyv43S1AFrr+UCtUuowB2txBKXUFUqpRxP+Hhb70L5dWbytlDq9k+eNjtUar3OVUurvSqnMhGmWKKWaOviVdFFHkVLqLqXUSqVUo1LqK6XUqUk6bDHQKqWOVkr9Lzbvtzt5vOPr/o9ezGtfpdQ8pVS9UmqxUmp6N/OarJT6t1KqWimlOzyWpZS6Tym1VCm1Xin1uVLq4C14/izWv0al1FtKqVEd5nd/rK6VSqkLXJzXUKXUbKXUitjru1gp9aBSaqvY4/H307wOzxuglGpVSi2J/Z34Xop2eH+d0MXiO/4f76CUelcpVaeUWq6UuryLmv8cq2m/hPu+6VBDu1LqxW68j4/1b4NS6jmlVHF3r5Of0Fq3APcDf/C6FkEQkiP0AU5rvUxrnR//id39o4T73lVK7QL01Vp/2OHpidMVJdz/GDDDiVpSN+wZsS/o/wIvaK3P1VrrLT0HKIrVvR2wO3BWh8cPS/TTWq/oZLl9gDeAUbF59AVmAtdvKRikwDrgNhK+yDsh8XX/QWhNZl6xAPsscA+WxzHALUqpH3UxrzbgaeA3nTyWAVQAP4nN60/A00qp0Z3NSCk1AJgTm64Y+BR4KmGSK4AJWK/zPsAspdRBLsyrP/A/IBeYChQAOwHvAPt3mDxPKTU54e/jgfL4Hx3+V5ax+fvrsU6W3dn/8eNY7/NirNf2d0qpX3R43jjgV0BV4v1a620Tll8Qq+GfXXhvi/U+OAkYDDQCf+9sWh/zOHCyUirL60IEQdgyoQ9wSXIw1hdMsrwN/MyJDzqlVH+l1Aux0Y6PgXEOLGMc1pfa41rrWak+X2u9Gngd2KYHiz8JGAkcpbUu11q3aa1fA84FrlJKFcZqHKGUmqOUWqOUWquU+ptSamvgbmB3lTAq2kl9b2itnwZ+ECBTZQvzKgYKgUe0xSfAd3TxumitF2it7wO+6eSxDVrrK7TWS7TWUa31S1hhZucuSjsS+EZr/U+tdTNWyPpRfJQL+DVwtda6Rmv9HTAbOMWFef0eqAdO0lovir0utVrrB7TWd3SY9hHg5IS/fw083MV8k6Gz/+PRwGNa64jWehHwHrBth2n+hjXy1NrNvPcGBgHPdPH4CcCLWuv/aq0bsMLwkUqpAtg4Qj1TKTU/NkJ3n1JqsFLq1diI6xtKqX5dLVwp9XOl1BdKqVpljQhv38V0aUqpi5VSi2L/N0/HRwKVUq8ppc7uMP2XSqkjAbTWy4EaYLduXgdBEHyCBDiL7YAFyU6sta7EGk2Z5EAtdwLNwFDgtNiPnYzFCm/3aK3/1JMZKGvT6IFAxxHLZNgfeFVrvaHD/c8A2VjhLB14CViK9QU8DHgyFh5+C3zQyahoqvw3tjlwTlejXFtCa70KeAI4VSmVrpTaHWuUqtf7LCqlBgMTSQh7sS/vvWJ/bgt8mVDLBmARsG0sCJQkPh67va0T8+rAfsCzWutoEpqPAsfGXrutsUa5PkrieV3R2f/xbcCvlVKZSqlJWKO+b8QfVEodBbRqrV/ZwrxPBv7Vyfs2TsfXcBFWIJyYMM00rPf/ROAw4FXgUmAA1mfxuZ3NWCm1E9bmzRlAf6yRvhe6WIE8F/gl1mhjCVYguzP22OPAcQnz3Qbr/fpywvO/A7oaQRYEwUdIgLMoAtZ3cv+82BddrVLq9g6PrY89zzZiwWUacHlsROZr4CE7lwFMBvLYfBNZslTHRr0qgQ3Avzo8/lzC69XVjuQD6LCpCkBr3Q5Uxx7/MdaXz8zY69CstbbzQI6fYAXDrbBG1l5SSmX0cF5PAJcDLcC7wGVa64reFBfbNPsY8JDW+vv4/VrrooTXIR+o6/DUOqwQlJ/wd8fHbJ9XBwYAKxNcfhF7P6xXSv2nw7TLsQLXflgBqTejb9D5//FLWJtHm4DvgftiI6UopfKx9pk7v7uZKqVyY/N4sJvJunsN49yhtV4VWwF8F/hIa/15bP+zZ4Edu5j3GVgrXB/FRhIfwnq/dTZSNgPrPbg8Nt8rgF/F3t/PAjuoTfs3ngDMiU0Xx/bPNUEQnEECnEUNnX8h7RT7oivSWndcOy4AOt2E1wsGsmlfqDhLbV7GC1hr828mfJAny4DYqFcu8D7wWofHf5nwev2yi3lUY40ubkbsC2ZA7PERwNJYqLOd2GauVq11LXAeMAbYOtX5xDYxPoW16a8P1ijMLKXUoT2tTSmVhrVpsRU4u5tJG7A23yZSiPUF3JDwd8fHnJ7XWhL6q7V+Ifae+T3Wa9SRh7E2xx6HNSLXGzb7P45tOnwNuAprdHcEcKBS6szYJFdibf4u7zijDhyJtS9kd7tZdPcaxlmVcLupk7/z6ZxRwIUJK0e1WC6dHSQ0Cng2YbrvgAgwWGu9Hmu07djYtMdirSgk4sTnmiAIDiABzmI+m2/q6JbYJsQ+pLDZNUnWAO1YH85xRtq8DLTWF2CNTLyplBrWg+c3YY1G7B7bAT4V3gAOVkrldbh/GtaowodYAXZkF6NiyRxskSoaUD143mRggdb637H91hZgfUF2e/RoVyilFHAf1k7w07TWbd1M/g0Jm7pir+c4rH3ZarBGORM3hf2ITva9c2Bec4FfxoJoMjwDHAos1lr3dmWl4//xWCCitX5Ya90e28frSeCQ2OM/A86NbUpfifV/97RSquORmCcDD2/hQJ+Or+FYIAso7ZWRRQVwbcLKUZHWOldr/UQX0x7cYdrs2KgfWCPGx8U29+cAb3V4/tZsvrlcEASfIgHO4hWszWrJ8lPgzQ6bHnqN1jqCdTTgFUqp3Ng+Kidv4Wk95WzgTWBubH+rOBlKqeyEn8yOT4zte3MS1qaytSku9xGsTWf/VNbpJDKVUgcCtwNXaK3rgI+xQsP1Sqm8WB17xp6/ChiurKNZOyW2T1U21mhmWqKHUmpbZZ1aIj22Ce1mrE3C36U6L+BzYIKyTiWilHVwyM/p4gswNk02sZGo2LwS92O6C+sL9LBYSO6OZ4HJSqlpsXleDsxP2OT6MPBHpVS/2EjhGXS9CdDOed0C9AMeUUqNizkXADt0NnFsn7J9ge6OBE6Wjv/HpVgv+/HK2rl/CNaRwvH+/AwrhO8Q+1mBtQkyvs8YSqnhWEfebmlXhseAw5RSU2MB+CqszZNdjVSmwmzgt0qpXWOvZ55S6tDY69qRu4Fr46PrSqmBSqnDEx5/BWuU7irgqcR9FWMrc8X0bN9WQRDcRmtt1A/WaMv4Tu7/BNh1S9PFHnsZ+IUTtWBtRn0J60i+j4Grgfdscr8CeDTh7zSsL+evsDZ
      "text/plain": [
       "<Figure size 720x864 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 = ob.getParameters()\n",
    "    if set(parm_array) & MAN_PARAMS:\n",
    "        manGeos.append(ob)\n",
    "        prMan = np.append(prMan,ob.getNumber(\"prMan\"))\n",
    "        tpMan, tpUnit = np.append(tpMan,ob.getNumber(\"tpMan\")), ob.getUnit(\"tpMan\")\n",
    "        tdMan, tdUnit = np.append(tdMan,ob.getNumber(\"tdMan\")), ob.getUnit(\"tdMan\")\n",
    "        wdMan = np.append(wdMan,ob.getNumber(\"wdMan\"))\n",
    "        wsMan, wsUnit = np.append(wsMan,ob.getNumber(\"wsMan\")), ob.getUnit(\"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",
    "wpres = (wpres/100) * units.mbar\n",
    "u,v = wind_components(spd * units.knots, np.deg2rad(direc))\n",
    "\n",
    "if tpUnit is 'K':\n",
    "    T = (tpSig-273.15) * units.degC\n",
    "    Td = (tdSig-273.15) * units.degC\n",
    "    tman = tman * units.degC\n",
    "    dman = dman * units.degC\n",
    "\n",
    "# Create SkewT/LogP\n",
    "plt.rcParams['figure.figsize'] = (10, 12)\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
}