python-awips/examples/notebooks/Surface_Obs_Plot_with_MetPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the MetPy example [\"Station Plot with Layout\"](http://metpy.readthedocs.org/en/latest/examples/generated/Station_Plot_with_Layout.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pprint\n",
    "from awips.dataaccess import DataAccessLayer\n",
    "from metpy.calc import get_wind_components\n",
    "from metpy.cbook import get_test_data\n",
    "from metpy.plots.wx_symbols import sky_cover, current_weather\n",
    "from metpy.plots import StationPlot, StationPlotLayout, simple_layout\n",
    "from metpy.units import units\n",
    "\n",
    "def get_cloud_cover(code):\n",
    "    if 'OVC' in code:\n",
    "        return 1.0\n",
    "    elif 'BKN' in code:\n",
    "        return 6.0/8.0\n",
    "    elif 'SCT' in code:\n",
    "        return 4.0/8.0\n",
    "    elif 'FEW' in code:\n",
    "        return 2.0/8.0\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "state_capital_wx_stations = {'Washington':'KOLM', 'Oregon':'KSLE', 'California':'KSAC',\n",
    "                             'Nevada':'KCXP', 'Idaho':'KBOI', 'Montana':'KHLN',\n",
    "                             'Utah':'KSLC', 'Arizona':'KDVT', 'New Mexico':'KSAF',\n",
    "                             'Colorado':'KBKF', 'Wyoming':'KCYS', 'North Dakota':'KBIS',\n",
    "                             'South Dakota':'KPIR', 'Nebraska':'KLNK', 'Kansas':'KTOP',\n",
    "                             'Oklahoma':'KPWA', 'Texas':'KATT', 'Louisiana':'KBTR',\n",
    "                             'Arkansas':'KLIT', 'Missouri':'KJEF', 'Iowa':'KDSM',\n",
    "                             'Minnesota':'KSTP', 'Wisconsin':'KMSN', 'Illinois':'KSPI',\n",
    "                             'Mississippi':'KHKS', 'Alabama':'KMGM', 'Nashville':'KBNA',\n",
    "                             'Kentucky':'KFFT', 'Indiana':'KIND', 'Michigan':'KLAN',\n",
    "                             'Ohio':'KCMH', 'Georgia':'KFTY', 'Florida':'KTLH',\n",
    "                             'South Carolina':'KCUB', 'North Carolina':'KRDU',\n",
    "                             'Virginia':'KRIC', 'West Virginia':'KCRW',\n",
    "                             'Pennsylvania':'KCXY', 'New York':'KALB', 'Vermont':'KMPV',\n",
    "                             'New Hampshire':'KCON', 'Maine':'KAUG', 'Massachusetts':'KBOS',\n",
    "                             'Rhode Island':'KPVD', 'Connecticut':'KHFD', 'New Jersey':'KTTN',\n",
    "                             'Delaware':'KDOV', 'Maryland':'KNAK'}\n",
    "single_value_params = [\"timeObs\", \"stationName\", \"longitude\", \"latitude\", \n",
    "                       \"temperature\", \"dewpoint\", \"windDir\",\n",
    "                       \"windSpeed\", \"seaLevelPress\"]\n",
    "multi_value_params = [\"presWeather\", \"skyCover\", \"skyLayerBase\"]\n",
    "all_params = single_value_params + multi_value_params\n",
    "obs_dict = dict({all_params: [] for all_params in all_params})\n",
    "pres_weather = []\n",
    "sky_cov = []\n",
    "sky_layer_base = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange\n",
    "from datetime import datetime, timedelta\n",
    "lastHourDateTime = datetime.utcnow() - timedelta(hours = 1)\n",
    "start = lastHourDateTime.strftime('%Y-%m-%d %H')\n",
    "beginRange = datetime.strptime( start + \":00:00\", \"%Y-%m-%d %H:%M:%S\")\n",
    "endRange = datetime.strptime( start + \":59:59\", \"%Y-%m-%d %H:%M:%S\")\n",
    "timerange = TimeRange(beginRange, endRange)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"obs\")\n",
    "request.setParameters(*(all_params))\n",
    "request.setLocationNames(*(state_capital_wx_stations.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "response = DataAccessLayer.getGeometryData(request,timerange)\n",
    "for ob in response:\n",
    "    avail_params = ob.getParameters()\n",
    "    if \"presWeather\" in avail_params:\n",
    "        pres_weather.append(ob.getString(\"presWeather\"))\n",
    "    elif \"skyCover\" in avail_params and \"skyLayerBase\" in avail_params:\n",
    "        sky_cov.append(ob.getString(\"skyCover\"))\n",
    "        sky_layer_base.append(ob.getNumber(\"skyLayerBase\"))\n",
    "    else:\n",
    "        for param in single_value_params:\n",
    "            if param in avail_params:\n",
    "                if param == 'timeObs':\n",
    "                    obs_dict[param].append(datetime.fromtimestamp(ob.getNumber(param)/1000.0))\n",
    "                else:\n",
    "                    try:\n",
    "                        obs_dict[param].append(ob.getNumber(param))\n",
    "                    except TypeError:\n",
    "                        obs_dict[param].append(ob.getString(param))\n",
    "            else:\n",
    "                obs_dict[param].append(None)\n",
    "\n",
    "        obs_dict['presWeather'].append(pres_weather);\n",
    "        obs_dict['skyCover'].append(sky_cov);\n",
    "        obs_dict['skyLayerBase'].append(sky_layer_base);\n",
    "        pres_weather = []\n",
    "        sky_cov = []\n",
    "        sky_layer_base = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use pandas to retrieve desired subsets of our observations.\n",
    "\n",
    "In this case, return the most recent observation for each station."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "df = pandas.DataFrame(data=obs_dict, columns=all_params)\n",
    "#sort rows with the newest first\n",
    "df = df.sort_values(by='timeObs', ascending=False)\n",
    "#group rows by station\n",
    "groups = df.groupby('stationName')\n",
    "#create a new DataFrame for the most recent values\n",
    "df_recent = pandas.DataFrame(columns=all_params)\n",
    "#retrieve the first entry for each group, which will\n",
    "#be the most recent observation\n",
    "for rid, station in groups:\n",
    "    row = station.head(1)\n",
    "    df_recent = pandas.concat([df_recent, row])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert DataFrame to something metpy-readable by \n",
    "attaching units and calculating derived values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = dict()\n",
    "data['stid'] = np.array(df_recent[\"stationName\"])\n",
    "data['latitude']  = np.array(df_recent['latitude'])\n",
    "data['longitude'] = np.array(df_recent['longitude'])\n",
    "data['air_temperature'] = np.array(df_recent['temperature'], dtype=float)* units.degC\n",
    "data['dew_point'] = np.array(df_recent['dewpoint'], dtype=float)* units.degC\n",
    "data['slp'] = np.array(df_recent['seaLevelPress'])* units('mbar')\n",
    "u, v = get_wind_components(np.array(df_recent['windSpeed']) * units('knots'),\n",
    "                           np.array(df_recent['windDir']) * units.degree)\n",
    "data['eastward_wind'], data['northward_wind'] = u, v\n",
    "data['cloud_frac'] = [int(get_cloud_cover(x)*8) for x in df_recent['skyCover']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10e95ec10>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABiIAAAPmCAYAAABpY4e+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4FMUfx/H3pNAh9C499F6kKU0RUPQnTRCQIgoqIFgo\nCopSBLGgIIIgVaUIKkUUC71I71JEikhvoZNAcvv7YzfxCCH1Lkf5vJ5nn5C92Znv7u2FZGbnO8ay\nLERERERERERERERERLzBz9cBiIiIiIiIiIiIiIjI3UsDESIiIiIiIiIiIiIi4jUaiBARERERERER\nEREREa/RQISIiIiIiIiIiIiIiHiNBiJERERERERERERERMRrNBAhIiIiIiIiIiIiIiJeo4EIERER\nERERERERERHxGg1EiIiIiIiIiIiIiIiI12ggQkREREREREREREREvEYDESIiIiIiSWCMWWqM2e7r\nOHzBOfclvo7jVowxDY0xW4wxV40xLmNMBl/HJLdmjCngvE/tPVzvQWPMJE/WKSIiIiIJo4EIERER\nES8yxnRwOtZcxpiaMbxujDH/Oq/P91IMuY0x7xhjysWzvHvMLmPMdWPMEWPMVGNMPm/E6C3GmBrG\nmAHGmKAEHtfYGLPQGHPa6cTeY4z5wBiT+RaHWB4I97ZkjCnp3D/5Y3jZ4jY9d2NMFuBb4DLwEtAW\nuOLlNssYY2Y7Hd9XjTGHjTG/GmO6RSv3pjHmf0loJ7b3JMmMMTmMMR8aY3YbYy4bYy4ZYzYYY/ol\n9LOUCDfcU8aYR40xAzxZp4iIiIgkPw1EiIiIiCSPq0DrGPbXBvIAYXivoyw38DYQr4EIN29hd952\nAX4CngaWG2PSeDY8r6oBDADi3XlqjPkQmAdkB4YBXYHfgW7AVmNMUS/EeTsriX3/xNTpXR94JHnD\nibcqQDrgLcuyJlmWNc2yrHBvNWaMqQFsAMoA47Dvm/GAC3g5WvE3gUQPRBD7e5IkxpgqwA7gRWAZ\n8ArwKrAZ6Is9uOMVlmUdBFIDX7vtfhT7MywiIiIid7AAXwcgIiIico/4GWhhjHnZsqwIt/2tgY1A\n1mSIwSSw/M+WZW1y/j3RGHMa6AM8CUzzaGTeF69zN8Y8jd3pOgNoY1lW5ODQRGPMZGAJMMsYUzHa\n++hzxphUQJhbzB5vIvoOb3bse0B25+t5T1VojEljWdatZlX0A0KAKpZlXYh2XPTPt0XCP48xhuSB\nOv6rzJiMwA/AdaCCZVl/ub08zhjTD3jOk21GZ1nWtZh2e7NNEREREfE+zYgQERERSR7TgSzYT5AD\nYIxJATQDvonpAGNMWmPMR07qplAnTcprMZSrb4xZaYwJMcZcdMoNcV6rA6xzik5yS7fULhHnsNL5\nel+09os76WjOOOlo1htjHo8hzozGmBFO2ppQ57ymOCl0IsukNMa8a4z52ylzyBjzvnOt3OtyGWNG\nGWOeNMbscMruMMY0cCvzDjDc+faA27nHll5qAHAW6By9Q9+yrPXA+9hPvDeP4fwqGWNWG2OuGGP2\nG2O6xFCmuzHmTyfdzVnnWj0drUweY8xEY8wJt/PqGK1MHedcWhpjBhtjjmCnIKp4q/fXGNPAee1R\n5/v8xpjPnbRTV5w0VN+6p/sxxnTgvyfgl7hdw1rO6zetEWGMyW6MmeDEf9XYazS0i1Ymci2A14wx\nnY0x+5xzXWeMqRytbE5jzCRjpzkKNcYcNcbMMbGkJTLGLAUmO9+ud9qa5PZ6C2PMRue8TxljvjLG\n5I5Wx2Tn81TIGPOTMeYCt/isOgoDf0YfhACwLOu0W70uIC3Q3u16TnJeS/J74pRpZIxZYeyUSheM\nMT8aY0rGEnukLtgzqF6NNggReR4nLct6z62d/xljFhg7dVuo87ntb4y54e9M5z7ZHtdnxO2+aOd8\nPxk7rZZxO0+XW/nXnfpOO3VuMMY0i+skjTGBxk7Ztte5R0871+vheFwjEREREUkEzYgQERERSR4H\ngT+w0xstdPY1wk4ZNBPo6V7YGGOw0wPVAb4EtgANgQ+MMXksy3rVKVcK+NF5/S3sFE/B2CmJAHZi\np3AZCHwBrHD2r07EORRwvh53i7MUsAr4FxiK3RneEphjjGlmWdYcp1w6p+3iwARgE5ANeBw7NdUZ\np/NyHlDTiXUXUBY7NUxRoEm0eB4AmgKjgUvY6W++M8bksyzrLPCdcy2exr6+kZ3Bp4mBMSbYaWeS\nZVmXbnENpgLvAo2x37dImYEFzr5vnGswxhhzzbKsyE7m54FPgVnACCAVdrqs+7EHqjDG5ADWABHA\nSOAUdmqaCcaYDJZlfRotnsj3fDiQEvv93g885cTqriX2IMsvzveVgerYs1sOAwWx0/EsNcaUtCzr\nKnZqnpHY13YI9nuC21e4MZ9/amApdqf8KOCAE8tkY0xGy7JGRoupNZAeGON83xv43hhTyG22xXfY\nqYhGYn+OcgAPYw+I/UPMBmMPFnV2rtEBYJ8TYwdgIvYAXV8gJ9ADqGmMqWBZlvsMigDneq0AXiP2\nNSYOAtWNMaUsy/ozlnLPYH+m12KncCIyNux0Uol9T3Y75/cM9iDMQuzrmdapY6Vzfre6ZgBPOOc4\nO5Yy7toDF4CPsD+DD2H/rMngtB3JAjIRx2ckBmOBXNgDuG1jeP1lYC7wFZAC+7M+yxjT2LKsn2KJ\n+x3s93489n0QhP15qICdhk1EREREPM2yLG3atGnTpk2bNm1e2oAO2DniK2I/2XseSOm89i3wu/Pv\ng8A8t+P+5xz3RrT6vsXupC7kfN/TKZc5lhgqO2XaJTDmetgpo/Jiz9w4CRwF0ruV/R17ECQwWh0r\ngT1u37/r1Pm/WNptC4QDNaLt7+wcW91tnwt73Y2CbvvKOPu7uu173dmXLx7nHXnNX46j3Hlgvdv3\nS53jerrtC8QebDkO+Dv75gDb4qj7S+wO6EzR9k/DTvsTee/UcdrcG7nPrewQ7MGJjG77UjjHj3fb\nlyqG9qs69bZ129fc2VcrhvJLgcVu3/dwyj7tti8Ae7DqApDO2VfAKXcSCHIr+7iz/zHn+4zO968m\n5bMX7X05AWwFUrjtf9Qp+47bvsnOviHxbO9h7JRG17EH+t7H7kAPiKHsRWBiDPuT9J5gr4kRAoyN\ntj+7s/+LOM7hLLApAdc4pnjHYA9KuF/fpcT+GQmIdl+0cyv3GeC6RfvR7/0AYBvOz1W3/Qfcrzf2\nz6x5cZ2fNm3atGnTpk2bNs9tSs0kIiIikny+xV6ItbExJj32U/W3WmvhUexO+ehPkH+EnRe+ofN9\niPP1yejpUDzgd+yO4kPYT/H/CzxoWdZFAGNMZqCu81qQMSZr5Ab8CgQbY3I5dTUDtliWNTeW9lpg\nP929J1pdkal/6kaPz7KsA5HfWJa1Hbuzu2Aizze98/ViHOUuYj/x7e469iyOyFgiv88OVHJ2hwD3\nRU89FMmZBdMMmA/4x3A9g7AHtNxNsSwrLNq+mdidvE3d9j3Cf7NvImMMdWs70NgpsvYB57CfDE+M\nR4FjlmVNd2sn8j5Oh704+w2xWjfOQIhM/xX5Hl4FrgF1jb1+QVJVxp6J87nlthaBZT89vxt4LIZj\nxsSw7yaWZf2OPZthHvZMnl7YsymOmBhSld2ijqS+J/Wx3+cZ0e4fF/aT/9E/Q9FlIO77/1bxpnfa\nWgmkAYpFKx6fz0iCuN/7xphM2ANXK7n5cxJdCFDaGFMkMe2KiIiISMJpIEJEREQkmVh2nvjfgTbY\nncR+3DoFSn7gqGVZl6Pt3+32Otgdy6uwn6Q/boyZ7uS/98Qiti9hP+XdHPgJKI+dRihSEexBkUHY\nAxbu2zvY6VgiFwwuDOyIo71goBR2OiL3uvY4dWWLVv5QDHWEYKeASYzIDtj0sZayX4/eWXvUstPm\nuNvrfC3gfH0f+0nxdcaYv4wxnxljariVz4bdidyFm6/nRG68npEORPsey7K2Yd8nLd12t8S+rosj\ndxhjUhtjBhpj/gVC+e+
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d383e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as feat\n",
    "from matplotlib import rcParams\n",
    "rcParams['savefig.dpi'] = 100\n",
    "proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,\n",
    "                             standard_parallels=[35])\n",
    "state_boundaries = feat.NaturalEarthFeature(category='cultural',\n",
    "                                            name='admin_1_states_provinces_lines',\n",
    "                                            scale='110m', facecolor='none')\n",
    "# Create the figure\n",
    "fig = plt.figure(figsize=(20, 15))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=proj)\n",
    "\n",
    "# Add map elements \n",
    "ax.add_feature(feat.LAND, zorder=-1)\n",
    "ax.add_feature(feat.OCEAN, zorder=-1)\n",
    "ax.add_feature(feat.LAKES, zorder=-1)\n",
    "ax.coastlines(resolution='110m', zorder=2, color='black')\n",
    "ax.add_feature(state_boundaries)\n",
    "ax.add_feature(feat.BORDERS, linewidth='2', edgecolor='black')\n",
    "ax.set_extent((-120, -70, 20, 50))\n",
    "\n",
    "# Start the station plot by specifying the axes to draw on, as well as the\n",
    "# lon/lat of the stations (with transform). We also set the fontsize to 12 pt.\n",
    "stationplot = StationPlot(ax, data['longitude'], data['latitude'],\n",
    "                          transform=ccrs.PlateCarree(), fontsize=12)\n",
    "\n",
    "# The layout knows where everything should go, and things are standardized using\n",
    "# the names of variables. So the layout pulls arrays out of `data` and plots them\n",
    "# using `stationplot`.\n",
    "simple_layout.plot(stationplot, data)\n",
    "\n",
    "# Plot the temperature and dew point to the upper and lower left, respectively, of\n",
    "# the center point. Each one uses a different color.\n",
    "stationplot.plot_parameter('NW', np.array(data['air_temperature']), color='red')\n",
    "stationplot.plot_parameter('SW', np.array(data['dew_point']), color='darkgreen')\n",
    "\n",
    "# A more complex example uses a custom formatter to control how the sea-level pressure\n",
    "# values are plotted. This uses the standard trailing 3-digits of the pressure value\n",
    "# in tenths of millibars.\n",
    "stationplot.plot_parameter('NE', np.array(data['slp']),\n",
    "                           formatter=lambda v: format(10 * v, '.0f')[-3:])\n",
    "\n",
    "# Plot the cloud cover symbols in the center location. This uses the codes made above and\n",
    "# uses the `sky_cover` mapper to convert these values to font codes for the\n",
    "# weather symbol font.\n",
    "stationplot.plot_symbol('C', data['cloud_frac'], sky_cover)\n",
    "\n",
    "# Also plot the actual text of the station id. Instead of cardinal directions,\n",
    "# plot further out by specifying a location of 2 increments in x and 0 in y.\n",
    "stationplot.plot_text((2, 0), np.array(obs_dict[\"stationName\"]))\n",
    "\n",
    "plt.title(\"Most Recent Observations for State Capitals\")"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}