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": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAALtCAYAAABHKTIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4FNX6gN+TSgiE3hICCV3hIqKiiIj16s+GFVEsWFApIioqCLYLqCBWFAtSVBDswEWwoChiuYAFKUIoqUAglBBIQur8/pgNLjFly8zuzsz3Ps8+yU45873zbbLfnDNFaZqGIAiCIAiCYC3Cgh2AIAiCIAiC4D1SxAmCIAiCIFgQKeIEQRAEQRAsiBRxgiAIgiAIFkSKOEEQBEEQBAsiRZwgCIIgCIIFkSJOEARBEATBgkgRJwiCIAiCYEGkiBMEQRAEQbAgUsQJgiAIgiBYECniBEEQBEEQLIgUcYIgCIIgCBZEijhBEARBEAQLIkWcIAiCIAiCBZEiThAEQRAEwYJIEScIgiAIgmBBpIgTBEEQBEGwIFLECYIgCIIgWBAp4gRBEARBECyIFHGCIAiCIAgWRIo4QRAEQRAECyJFnCAIgiAIggWRIk4QBEEQBMGCSBEnCIIgCIJgQaSIEwRBEARBsCBSxAmCIAiCIFgQKeIEQRAEQRAsiBRxLpRS/1ZKLfRgue5KqZ8CEZNre4OVUqsCtT1BsDKe/h3bCaVUmlLqAj/biFZKbVZKNTcqLkEQzMf2RZxSqo1S6ojbS1NK5bu97+ta9GngWbf1Ki+XC6Bp2p9ArlLqchNjMQWl1JNKqblu7xNc/7hfUTrfKaXurGK9JFesFXHuUUpNV0pFui2TppQqrOQXX00cDZVSryulspVSBUqp9Uqp2zx0qLWodX0hzVJK5bm28UCl+Xcqpba5Yvyiujhdyw5QSv3kivO7Kuafp5T6zbWtHUqpu2poq5tS6kul1D6llFZFzDOVUulKqcNKqd+VUv9Xi+f5rvwVKKVWKKXaeroPTG6rlVJqhlJql2sf71BKzVFKdXHNr/g8/VZpvaZKqWKlVJrrvftnqbzS52tQNZuv/HfcQyn1g1LqkFIqSyn1eDUxP+GK6QK3aRsrxVCqlPpvDd43uvKXr5RaqJRqXNN+CiU0TSsCZgGPBDsWQRA8x/ZFnKZpGZqm1at4uSaf5DbtB6XUaUADTdN+qbS6+3IN3abPA+42IxbvDX3D9SW9ElisadpITdO02tYBGrri/hfQGxheaf7l7n6apu2qYrtRwHKgrauNBsBDwLO1FQde8CTQ0bWNc4GHlVIXu7bfD/2Lvj/QGEgF5tfQ1gHgJdwKAzeXSOAz4E2Xx/XAC0qpk6ppqwT4ELijinkRQCbQz9XWY8CHSqmkqhpSSjUFPnUt1xhYC3zgtsiTVLMPTG6rCfATUBfoC9QHegLfAxdWWjxWKdXN7f2N6PkAoNLfSgbHf77mVbHtqv6O30f/nDdG37dDlVJXVFqvPXAtsNt9uqZpXd22X98Vw0fVeHdF/xzcDLQACoDpVS0bwrwP3KqUig52IIIgeIimaY56ARrQodK0x4G3a1vObV4CUAhEmxBLE2AxkAesBiYAqwxyfxKYC7QH0oEJleZ/B9xZxXpJrlgj3KZNAd5ye58GXOBBDHcAe4HYStOvB44Aca73ieiFRQ6wH3gVOAE4CpS5ls2tZhs7gX+7vZ8ALHD9PhV4zW1evMutfS1x3wl8V2laC9e6dd2mrQFuqKWtDvqfXq376k/gmmrm3QX85PY+1vWZ7FLbPjC5rYnAOiCsBq+Kz9N44Dm36WuBcUBaFevU+vmi6r/jAuBEt/cfAWMrLbMMuKSmbaAXgEcqf27d5j8NvO/2vj1QDNR3i/8hV07zgZmuz88y4DD6gU2jGtwuA/4ActGL5O5V7Rv0A/MxwHb0v5sPgcaueV8AIyq1uw642u39VqBfbZ9NeclLXqHxsn1PnIf8C9ji6cKapu1E71XpbEIsr6EXKq2A210vI2mH3jPxpqZpj/nSgGv48SKgcs+lJ1wILNM0Lb/S9E+AOkBvpVQ4sAS90ExCL5oXaJr2F3AP8LP2z97RitgaoRdm69wmrwO6VizieuH2HsC9R8gjNE3bg96Ld5tSKlwp1Ru9t8rvcxiVUi2ATsBGt2m5SqmzXG+74ubo2p/bga4e7AND26rEBcBnmqaVe6A5Fxjo2ncnoPd2/c+D9aqjqr/jl4BblFKRSqnO6L2/yytmKqWuA4o1TVtaS9u3Ah9X8bmtoPI+3I5exHVyW+Ya9M9/J+By9ALuUaApevE1sqqGlVI90Yc670Y/yHsTWFxNj9lI4Er0ojMeOIj+PwX0nrYb3No9Ef3z+rnb+n8B1fUkC4IQYkgRp9MQ/Wi4Mr+5vuxylVKvVJp32LWeYbiKl2uAxzVNy9c0bQPwjpHbQC9WYjl+uMxT9in93MCd6L0JH1eav9Btf1V3cnlTKg1bAWiaVgrsc83vhf4F9JBrPxzVNM3TwqhimPqQ27RD6AUCwFJggNIvUIlB773R0If/fGG+q40i4AdgnKZpmT62BRwbpp0HvKNp2uaK6ZqmNXTbD/U43hH+9qxtHxjaViWaAtluLle4Pg+HlVJfVVo2C73ougC9SHq3mjY9paq/4yXoQ6WFwGZgpqZpa1yx1UPvQRtVU6NKqbquNubUsFhN+7CCaZqm7XEdBP4A/E/TtN81/Xy0z4CTq2l7CPpB1/80TSvTNO0d9M/bGVUsezf6ZzDL1e6TwLVKqQjXNnq4ne84CPjUtVwFhv9fEwTBPKSI0zlI1V9KPV1fdg01Tat8lFwffWjDSJrx97lRFaQbvI3F6Ef137qfvO4hTV29X3WBH9GHZ9y50m1/XVlNG/vQexmPw/Ul09Q1PxFIdxV23nLE9TPObVocri93TdO+AZ5A7/lLRx+KOoxeUHiF60T9D4BbgCj03piHlVKX+hB3RZthwHvovTgjalj0CMc7wt+eNe4Dk9vaj1t+NU1b7PrM3I++jyrzLjAYvYdobhXzveG4v2PXhQVfAP9B7+VNBC5SSg1zLfIU8J6maamVG6rE1ejnRn5fwzI17cMK9rj9XljF+3pUTVvgQbcDpFx0l6ouyGkLfOa23F/opx+00DTtMHqv20DXsgPRDxbcMeP/miAIJiFFnM6fHD/sUSOu4cQovBiC9ZAcoBT9H3QFbQzeBpqmPYDeQ/GtUirBh/UL0XslertOiveG5cD/KaViK02/Br134Rf0IraNq7D7x+Zrie0gek+f+5DQSbgNS2qa9pqmaR01TWuOXsxFABu89AC9V3OLpmlfappWrmnaFvQvyRqvKq0OpZTi73OlrtE0raSGxTfi5ujan+2BjZ7sAxPb+ga40lWMesInwKXADk3T/D1gqfx33A4o0zTtXU3TSjVNywIWoJ//BnA+MNJ1xW02+t/dh0qpyldo3gq8q2laTZ+9yvuwHRANpPhlpJMJTHI7QGqoaVpdTdOquiAnE/i/SsvWcfX+gd5zfINr6D8GWFFp/RM4fuhcEIQQRoo4naXo55B4yjnAt5WGIfxG07Qy9JP5n1RK1XWds3KrkdtwYwTwLfCN6/yrCiKUUnXcXpGVV3Sdi3Mz+rDZfi+3+x56r9dHSr/VRKRS6iLgFeBJTdMOoV/QsRv9itVYVxx9XOvvAVq7rnKtjneB8UqpRq7esiG4hsJcbXVTOm2At4CXXcXKP3Cdr1UHvdALq7RPfgc6Kv02I8p1leNlVPMl6FqmDq4eKVdb7uc1vY7+JXq5q1Cuic+Abkqpa1xtPg786Tb8Wu0+MLmtF4BGwHtKqfYu5/pAj6oWdp1jdh76hSP+UvnvOAV9t9+olApTSrVEv4CmIj/noxfiPVyvXejDkRXnkKGUao1+RW5tpzXMAy5XSvV1FcH/QR+qrK7H0htmAPcopU537c9YpdSlrv1amTeASRW97EqpZkqp/m7zl6L31v0H+MD93EXXAV1jfDvXVRCEYBDsKysC/aKaq07Rryo8vbblXPM+B64wIxb0IdUlmHh1qtv7MPQv6PX
      "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 = 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'] = (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
}