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": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Apr 10 16 12:52:00 , Apr 10 16 12:52:00 )\n",
      "425\n",
      "85\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAADFcAAAfrCAYAAADTfaYrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAnNwAAJzcBVAXW8gAAIABJREFUeJzs3Xd0VOXaxuF7Jh1IgRCK9CCIoCBVSCgJSBUsgIICKkVQ\nioKAR7GhoIIUKUpRuhSBIFU4CpqgCSAgIE2q9BJaKAnpM98fnG+HUaROsifJ71rrrDXvk8m893Ay\niWvv/ezHYrfbBQAAAAAAAAAAAAAAAAAAAAAAkFtZzQ4AAAAAAAAAAAAAAAAAAAAAAABgJporAAAA\nAAAAAAAAAAAAAAAAAABArkZzBQAAAAAAAAAAAAAAAAAAAAAAyNVorgAAAAAAAAAAAAAAAAAAAAAA\nALkazRUAAAAAAAAAAAAAAAAAAAAAACBXo7kCAAAAAAAAAAAAAAAAAAAAAADkajRXAAAAAAAAAAAA\nAAAAAAAAAACAXI3mCgAAAAAAAAAAAAAAAAAAAAAAkKvRXAEAAAAAAAAAAAAAAAAAAAAAAHI1misA\nAAAAAAAAAAAAAAAAAAAAAECuRnMFAAAAAAAAAAAAAAAAAAAAAADI1WiuAAAAAAAAAAAAAAAAAAAA\nAAAAuRrNFQAAAAAAAAAAAAAAAAAAAAAAIFejuQIAAAAAAAAAAAAAAAAAAAAAAORqNFcAAAAAAAAA\nAAAAAAAAAAAAAIBcjeYKAAAAAAAAAAAAAAAAAAAAAACQq9FcAQAAAAAAAAAAAAAAAAAAAAAAcjWa\nKwAAAAAAAAAAAAAAAAAAAAAAQK5GcwUAAAAAAAAAAAAAAAAAAAAAAMjVaK4AAAAAAAAAAAAAAAAA\nAAAAAAC5Gs0VAAAAAAAAAAAAAAAAAAAAAAAgV6O5AgAAAAAAAAAAAAAAAAAAAAAA5Grumb2BxWKZ\nJOnhzN4HAAAAAAAAAAAAAAAAAAAAAADkPHa7PTSz98j05gpda6wIyYJ9AAAAAAAAAAAAAAAAAAAA\nAAAA7pjV7AAAAAAAAAAAAAAAAAAAAAAAAAD/xmKx2C0WS0xm7pEVkysMXj55Vax0hazcMsdJS03W\nqWMHlJqc5FAvEBikssHFZbFYTErmGuIuxWv/3n3GutB9ZZTPv4BpeWKP/6WEK3HGOp+fvypWKGta\nHgAAAAAAAAB3zy7p0lWL0tLNTnJ3rFapQF672TEAAAAAAAAAAJksPiFRe/fuV3pamlFzc3NX0VLl\n5enlY2KyzJOakqQTh/fKlp7xnr3z5FWlB8vJzY15BNnZxo1bsmyvLG2uKFa6ggaO+i4rt8yREq9e\n0fTPXteOjT8ZtQvnz6p4mfJaOHesChX0NzGduU7FxqlihVrGumrdFnrqpTdNy5Nw5ZI+7tVUcedO\nSZLiL19So8db683X2pmWCQAAAAAAAMDdS0yRftzuocTU7Hmjm0aVUhXkR4MFAAAAAAAAAORUa9ft\nUsfnOjs0VuQPuk99P5mrQsXKmJgs81w4c0IjB7Z1aKwodX9Frf7vNwoK9DMxGZyhadN2WdZgQRtO\nNuSTx1evvPe1Grfp4VDfvjlGDcKf1badh0xKZr6ihfMroEAhY33yyF4T00h5ff310sAxDhNFRn4y\nVL9vP2hiKgAAAAAAAAB3y8dTqlshTW6W7NmgcOQcpwUAAAAAAAAAIKf6709b9NwznRR/Oc6oBd1X\nWv0/W5hjGysux53V2EEdFHf2pFErWiJYK5ZOo7ECd4yzKNmU1c1NrbsO0gv9Rsrd3dOonz7+lx5v\n3lbfrVhnYjpzlSpb3nhsdnOFJJV/uLaaPNPTWKcmJ6lbt/5KTEoxMRUAAAAAAACAuxWYz66aZdPN\njnFXjp63ymYzOwUAAAAAAAAAwNm+W7FOL3Z4SYlXrxi1oqXKq/9nCxRYuLiJyTJPwpVLGvduJ505\nmXFj+sBCxbR0yQwVvy/QxGTIrmiuyObqNH5GfT+dq3z+Gb8ArsZf1ssvdtEno+aYmMw85cpnNFec\nP31MSVfjTUxzTcsOfVWqXGVjfXj/LvV7a6yJiQAAAAAAAADci9JBNj1YLPs1WKSkWXTqkuXWTwQA\nAAAAAAAAZBvfLPhJPTq/rJTkRKNW8v6H9cbwBfIvUNjEZJkn6Wq8vnj/RZ049KdR8y8QpMWLZ6pc\ncFETkyE7o7kiByhbqabeGrNMxUpXMGo2W7pGDB2sF14erOTkVPPCmaBSpQcc1qeO7jcpSQZ3D091\neXOcPL18jNr8mV/pu+/Xm5gKAAAAAAAAwL2oXCJdxfJnvzEQR85yagAAAAAAAAAAcopJ01eo76u9\nlJaWYtTKVqp57ebtfvlNTJZ5UlOSNGnIyzq8d6tRy+Prr28XzNDDFUuZmAzZHWdQcojAwsU1YOQi\nVa7d2KG+PGKOGj/eTafPXDQpWdar8nA5h/XJI3tNSuKoULEyevaVwQ61/q8P1MnYOHMCAQAAAAAA\nALgnFotUu1ya/H2yV4PFiTirUrPf0A0AAAAAAAAAwN+M+mKhBvV/QzZbxkHfB6vWU58hs+ST18/E\nZJknPS1VX3/aS3v/WGfUvHzy6ps501S7enkTkyEnoLkiB/HOk0893v1KTZ/p6VDf8fs6hTV8Rlu2\nHzQpWdaqVrmsw/rkkX0mJfmnkCbtVDW0ubG+eD5WXXu8J5vNbmIqAAAAAAAAAHfLw02qVyFNnu7Z\n5xhfus2iExc4PQAAAAAAAAAA2dmHw2dq6HuDZLdnHJ+uUruJXh08VV7eeUxMlnls6emaObq/dvy2\nxqh5eHrpq2mT1bBeZROTIafg7EkOY7Va9VTn/+il/p/L3d3TqMeeOKyWLZ5RxLJoE9NlDX+/PCpU\ntKSxPnnYNSZXSJLFYlGH14YpILCIUduw9geN/jLCxFQAAAAAAAAA7kU+b6lu+TRZLNmnweLIOU4P\nAAAAAAAAAEB2ZLPZ9eb7EzVm2FCHes2wJ/XyoAny8PAyKVnmstvtmjfhXW2KWmrUrG7uGjvhCz3R\n7FETkyEn4exJDvVoo9bqN/xb+QYUNGqJCVfUvXNXDflsVo6flFD6/oyxPiePuE5zhSTl9Q3QSwM+\nl8ViMWqfDn5P9Rq9qKEjZ2v/oVNGPS7uop59tpsmTpyuP//c79BdCAAAAAAAAMB1FPK3q3qZ9Fs/\n0UWcvmhRUqrZKQAAAAAAAAAAd8Jms6v3gFH6evxoh3rdZs/ppf6fy83dw6Rkmctut2vxtE8VvWqu\nUbNYLPp01Cg91ybMvGDIEll5/TTNFTlY8IPV9daYZSoeXNGo2W02jf50iDp1+0DJyTn3zNkDFR4w\nHl+OO6v4SxdMTPNPD1QJUeO2rxhrmy1dO7es06iPP1StavVVI6St3hr8lWbNWabVq9dq0KBPFBLS\nQpUq1VPPnm9q4cJlOnv2vInvAAAAAAAAAMDf3V/YpnKFs0eDhV0WHWV6BQAAAAAAAABkG+npNr30\nykeaN32yQ73RU131fJ9PZXVzMylZ5vvv/C+0epHj+37no4/V/cUWJiVCVjpz9mKW7eWeZTvBFAUK\nFVP/ERGaMbKf/lj/g1FfuXieGh0+pIXzxqlo4fwmJswclSqWc1ifPLJX5SvXMSnNjbXq+IaO7N+u\nvdti/vG1g3/+oYN//iFdN91Ckk6ditW8eYs1b95iSdJDD1VQeHhdhYeHqnbtGvLx8c6S7AAAAAAA\nAABurGrpdF1KtOjMZddvXDhyzqryRW1mxwAAAAAAAAAA3EJKapo6dH5Xa75f5FBv8fzratmhnyx/\nu940J/l56XQtmzXSofbawEHq3/sZkxIhK6Wn23Tq1Kks28/1z+7gnnn75FX3dyapWbteDvVdWzco\nLLytNm/bb1KyzPNI5fIO65NH9pqU5N+5e3jqtSHfqOt/vlD1ei3l5Z3nn0+6xRibnTv3aPz4KWrd\nurOCg2vo6adf1LhxX2vHjt2y2TgpCgAAAAAAAGQ1q1UKLZ+mfF5ZN6L6bp2Ptyo+yewUAAAAAAAA\nAICbSUxKUev2b/yjsaJ
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe54eff6150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from awips.dataaccess import DataAccessLayer\n",
    "\n",
    "from metpy.calc import get_wind_components\n",
    "from metpy.cbook import get_test_data\n",
    "from metpy.plots import StationPlot, StationPlotLayout, simple_layout\n",
    "from metpy.units import units\n",
    "\n",
    "# Initialize\n",
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "\n",
    "data,latitude,longitude,stationName,temperature,dewpoint,seaLevelPress,windDir,windSpeed = [],[],[],[],[],[],[],[],[]\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"obs\")\n",
    "\n",
    "\n",
    "#\n",
    "# we need to set one station to query latest time.  this is hack-y and should be fixed\n",
    "# because when you DON'T set a location name, you tend to get a single observation\n",
    "# that came in a second ago, so your \"latest data for the last time for all stations\"\n",
    "# data array consists of one village in Peru and time-matching is suspect right now.\n",
    "#\n",
    "# So here take a known US station (OKC) and hope/assume that a lot of other stations \n",
    "# are also reporting (and that this is a 00/20/40 ob). \n",
    "#\n",
    "request.setLocationNames(\"KOKC\")\n",
    "datatimes = DataAccessLayer.getAvailableTimes(request)\n",
    "\n",
    "# Get most recent time for location\n",
    "time = datatimes[0].validPeriod\n",
    "\n",
    "# \"presWeather\",\"skyCover\",\"skyLayerBase\"\n",
    "# are multi-dimensional(??) and returned seperately (not sure why yet)... deal with those later\n",
    "request.setParameters(\"presWeather\",\"skyCover\", \"skyLayerBase\",\"stationName\",\"temperature\",\"dewpoint\",\"windDir\",\"windSpeed\",\n",
    "                      \"seaLevelPress\",\"longitude\",\"latitude\")\n",
    "request.setLocationNames()\n",
    "response = DataAccessLayer.getGeometryData(request,times=time)\n",
    "print time\n",
    "PRES_PARAMS = set([\"presWeather\"])\n",
    "SKY_PARAMS = set([\"skyCover\", \"skyLayerBase\"])\n",
    "# Build ordered arrays\n",
    "wx,cvr,bas=[],[],[]\n",
    "for ob in response:\n",
    "    #print ob.getParameters()\n",
    "    if set(ob.getParameters()) & PRES_PARAMS :\n",
    "        wx.append(ob.getString(\"presWeather\"))\n",
    "        continue\n",
    "    if set(ob.getParameters()) & SKY_PARAMS :\n",
    "        cvr.append(ob.getString(\"skyCover\"))\n",
    "        bas.append(ob.getNumber(\"skyLayerBase\"))\n",
    "        continue\n",
    "    latitude.append(float(ob.getString(\"latitude\")))\n",
    "    longitude.append(float(ob.getString(\"longitude\")))\n",
    "    #stationName.append(ob.getString(\"stationName\"))\n",
    "    temperature.append(float(ob.getString(\"temperature\")))\n",
    "    dewpoint.append(float(ob.getString(\"dewpoint\")))\n",
    "    seaLevelPress.append(float(ob.getString(\"seaLevelPress\")))\n",
    "    windDir.append(float(ob.getString(\"windDir\")))\n",
    "    windSpeed.append(float(ob.getString(\"windSpeed\")))\n",
    "    \n",
    "    \n",
    "print len(wx)\n",
    "print len(temperature)\n",
    "\n",
    "\n",
    "# Convert\n",
    "data = dict()\n",
    "data['latitude']  = np.array(latitude)\n",
    "data['longitude'] = np.array(longitude)\n",
    "data['air_temperature'] = np.array(temperature)* units.degC\n",
    "data['dew_point_temperature'] = np.array(dewpoint)* units.degC\n",
    "#data['air_pressure_at_sea_level'] = np.array(seaLevelPress)* units('mbar')\n",
    "u, v = get_wind_components(np.array(windSpeed) * units('knots'),\n",
    "                           np.array(windDir) * units.degree)\n",
    "data['eastward_wind'], data['northward_wind'] = u, v\n",
    "\n",
    "# Convert the fraction value into a code of 0-8, which can be used to pull out\n",
    "# the appropriate symbol\n",
    "#data['cloud_coverage'] = (8 * data_arr['cloud_fraction']).astype(int)\n",
    "\n",
    "# Map weather strings to WMO codes, which we can use to convert to symbols\n",
    "# Only use the first symbol if there are multiple\n",
    "#wx_text = make_string_list(data_arr['weather'])\n",
    "#wx_codes = {'':0, 'HZ':5, 'BR':10, '-DZ':51, 'DZ':53, '+DZ':55,\n",
    "#            '-RA':61, 'RA':63, '+RA':65, '-SN':71, 'SN':73, '+SN':75}\n",
    "#data['present_weather'] = [wx_codes[s.split()[0] if ' ' in s else s] for s in wx]\n",
    "\n",
    "# Set up the map projection\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as feat\n",
    "from matplotlib import rcParams\n",
    "rcParams['savefig.dpi'] = 255\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, 10))\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((-118, -73, 23, 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 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}