python-awips/examples/notebooks/METAR_Station_Plot_with_MetPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise creates a METAR plot for North America using AWIPS METAR observations (datatype *obs*) and MetPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange\n",
    "from datetime import datetime, timedelta\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt\n",
    "from metpy.calc import wind_components\n",
    "from metpy.plots import StationPlot, StationPlotLayout\n",
    "from metpy.units import units\n",
    "import warnings\n",
    "%matplotlib inline\n",
    "warnings.filterwarnings(\"ignore\",category =RuntimeWarning)\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",
    "# Pull out these specific stations (prepend K for AWIPS identifiers)\n",
    "selected = ['PDX', 'OKC', 'ICT', 'GLD', 'MEM', 'BOS', 'MIA', 'MOB', 'ABQ', 'PHX', 'TTF',\n",
    "            'ORD', 'BIL', 'BIS', 'CPR', 'LAX', 'ATL', 'MSP', 'SLC', 'DFW', 'NYC', 'PHL',\n",
    "            'PIT', 'IND', 'OLY', 'SYR', 'LEX', 'CHS', 'TLH', 'HOU', 'GJT', 'LBB', 'LSV',\n",
    "            'GRB', 'CLT', 'LNK', 'DSM', 'BOI', 'FSD', 'RAP', 'RIC', 'JAN', 'HSV', 'CRW',\n",
    "            'SAT', 'BUY', '0CO', 'ZPC', 'VIH', 'BDG', 'MLF', 'ELY', 'WMC', 'OTH', 'CAR',\n",
    "            'LMT', 'RDM', 'PDT', 'SEA', 'UIL', 'EPH', 'PUW', 'COE', 'MLP', 'PIH', 'IDA', \n",
    "            'MSO', 'ACV', 'HLN', 'BIL', 'OLF', 'RUT', 'PSM', 'JAX', 'TPA', 'SHV', 'MSY',\n",
    "            'ELP', 'RNO', 'FAT', 'SFO', 'NYL', 'BRO', 'MRF', 'DRT', 'FAR', 'BDE', 'DLH',\n",
    "            'HOT', 'LBF', 'FLG', 'CLE', 'UNV']\n",
    "selected = ['K{0}'.format(id) for id in selected]\n",
    "data_arr = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDEX Request\n",
    "edexServer = \"edex-cloud.unidata.ucar.edu\"\n",
    "DataAccessLayer.changeEDEXHost(edexServer)\n",
    "request = DataAccessLayer.newDataRequest(\"obs\")\n",
    "availableProducts = DataAccessLayer.getAvailableParameters(request)\n",
    "\n",
    "single_value_params = [\"timeObs\", \"stationName\", \"longitude\", \"latitude\", \n",
    "                       \"temperature\", \"dewpoint\", \"windDir\",\n",
    "                       \"windSpeed\", \"seaLevelPress\"]\n",
    "multi_value_params = [\"presWeather\", \"skyCover\", \"skyLayerBase\"]\n",
    "pres_weather, sky_cov, sky_layer_base = [],[],[]\n",
    "params = single_value_params + multi_value_params\n",
    "obs = dict({params: [] for params in params})\n",
    "\n",
    "request.setParameters(*(params))\n",
    "request.setLocationNames(*(selected))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use the Python-AWIPS class **TimeRange** to prepare a beginning and end time span for requesting observations (the last hour):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time range\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)\n",
    "\n",
    "response = DataAccessLayer.getGeometryData(request,timerange)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "station_names = []\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",
    "        # If we already have a record for this stationName, skip\n",
    "        if ob.getString('stationName') not in station_names:\n",
    "            station_names.append(ob.getString('stationName'))\n",
    "            for param in single_value_params:            \n",
    "                if param in avail_params:\n",
    "                    if param == 'timeObs':\n",
    "                        obs[param].append(datetime.fromtimestamp(ob.getNumber(param)/1000.0))\n",
    "                    else:\n",
    "                        try:\n",
    "                            obs[param].append(ob.getNumber(param))\n",
    "                        except TypeError:\n",
    "                            obs[param].append(ob.getString(param))\n",
    "                else:\n",
    "                    obs[param].append(None)\n",
    "        \n",
    "            obs['presWeather'].append(pres_weather);\n",
    "            obs['skyCover'].append(sky_cov);\n",
    "            obs['skyLayerBase'].append(sky_layer_base);\n",
    "            pres_weather = []\n",
    "            sky_cov = []\n",
    "            sky_layer_base = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next grab the simple variables out of the data we have (attaching correct units), and\n",
    "put them into a dictionary that we will hand the plotting function later:\n",
    "\n",
    "- Get wind components from speed and direction\n",
    "- Convert cloud fraction values to integer codes [0 - 8]\n",
    "- Map METAR weather codes to WMO codes for weather symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = dict()\n",
    "data['stid'] = np.array(obs['stationName'])\n",
    "data['latitude']  = np.array(obs['latitude'])\n",
    "data['longitude'] = np.array(obs['longitude'])\n",
    "data['air_temperature'] = np.array(obs['temperature'], dtype=float)* units.degC\n",
    "data['dew_point_temperature'] = np.array(obs['dewpoint'], dtype=float)* units.degC\n",
    "data['air_pressure_at_sea_level'] = np.array(obs['seaLevelPress'])* units('mbar')\n",
    "\n",
    "direction = np.array(obs['windDir'])\n",
    "direction[direction == -9999.0] = 'nan'\n",
    "\n",
    "u, v = wind_components(np.array(obs['windSpeed']) * units('knots'),\n",
    "                           direction * units.degree)\n",
    "data['eastward_wind'], data['northward_wind'] = u, v\n",
    "data['cloud_coverage'] = [int(get_cloud_cover(x)*8) for x in obs['skyCover']]\n",
    "data['present_weather'] = obs['presWeather']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['K0CO', 'KHOT', 'KSHV', 'KIND', 'KBDE', 'KPSM', 'KORD', 'KDFW', 'KPHL', 'KTTF', 'KBDG', 'KOLY', 'KNYC', 'KABQ', 'KLEX', 'KDRT', 'KELP', 'KRUT', 'KRIC', 'KPIT', 'KMSP', 'KHSV', 'KUNV', 'KSAT', 'KCLE', 'KPHX', 'KMIA', 'KBOI', 'KBRO', 'KLAX', 'KLBB', 'KMSO', 'KPDX', 'KTLH', 'KUIL', 'KTPA', 'KVIH', 'KBIL', 'KMLF', 'KCPR', 'KATL', 'KBIS', 'KCLT', 'KOKC', 'KRAP', 'KACV', 'KEPH', 'KELY', 'KFAR', 'KFAT', 'KMSY', 'KOLF', 'KPDT', 'KLMT', 'KHLN', 'KHOU', 'KICT', 'KIDA', 'KPIH', 'KPUW', 'KGJT', 'KGLD', 'KGRB', 'KLBF', 'KMLP', 'KBOS', 'KSYR', 'KDLH', 'KCOE', 'KOTH', 'KCRW', 'KSEA', 'KCAR', 'KDSM', 'KJAN', 'KSLC', 'KBUY', 'KLNK', 'KMEM', 'KNYL', 'KRDM', 'KCHS', 'KFSD', 'KJAX', 'KMOB', 'KRNO', 'KSFO', 'KWMC', 'KFLG', 'KLSV']\n"
     ]
    }
   ],
   "source": [
    "print(obs['stationName'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MetPy Surface Obs Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2wAAAJBCAYAAAA6KTYNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XdcVfX/wPHXh8sG2aIsBRX3QMWBe5TmTHPlLitLy/HL1NRvw2yXaZm5cuRIceTee++JGxVQFBERkT3uPb8/zgXZXFAU9fN8PHg84HzOZ5x7L3De57OEoihIkiRJkiRJkiRJxY/R826AJEmSJEmSJEmSlDMZsEmSJEmSJEmSJBVTMmCTJEmSJEmSJEkqpmTAJkmSJEmSJEmSVEzJgE2SJEmSJEmSJKmYkgGbJEmSJEmSJElSMSUDNkmSJEmSJEmSpGJKBmySJEkSAEKIl25jTiHEAiHEO0+prK+FEIufRlkFqHOPEOL9Iii3hRAi9GmXqy/7ghCixbOu90UlhHhfCLHnebdDkqTiSwZskvQKE0KYCSHmCiFChBAxQojTQoh2Wc5pLYS4LISIF0LsFkKUzZDWUwhxSJ+2J4fyWwkhTgkhHgkhbgghBufTnklCiAAhRKoQ4usc0vvo2xonhFgjhHDIo6zqQoitQoj7OQUi+hvhRCFErP7ryhOUtVgIEaa/zqt53WALITyFEEqGemOFEF9kSF8ghEjOkq7Jo7y83h8zIcQ8fbvuCiE+za2c/OiDFUUIMTzL8ZH641/rf24hhNBlaX+sEMJPfyOf9rM2y+s/PkOZLfRljsnntQsWQnxe2GuSioaiKNUURdnzpOUIId4RQhx4Ck2SJEl6ocmATZJebcbALaA5YAt8ASwXQngCCCGcgP/0xx2AE4B/hvwPgKnAj1kLFkKYAKuBWfqyewG/CSFq5dGea8AYYGMO5VXTl9UfKAXEA3/lUVYKsBx4L49zPlEUxVr/VekJyvoB8FQUxQboDHwrhKibR3kAdhnqnpQl7ecMadaKomhzKsCA9+drwBsoC7QExggh3sinXXm5CgzMcmyA/nhGd7K031pRlMP6G3lrRVGsgf1kfv2/z5B/IOpnK2tdaez0ZXQHvhBCvP4E1yRJhSZU8l5KkqQiJf/ISNIrTFGUOEVRvlYUJVhRFJ2iKBuAICAt2HgLuKAoygpFURJRA4BaQojK+vw7FEVZDtzJoXgHwAZYpKiOA5eAqnm05x9FUTYDMTkk9wXWK4qyT1GUWNQg5S0hRIlcyrqiKMpc4EK+L0Q+8itLUZQLiqIkpf2o/yr/pPUaIM/3BzWYmqQoSpSiKJeAOcA7T1DfccBSHzynBdEW+uNPhRDCEjUQ+xjwFkL45nauoignUN8Tnyeor6G+l/ihEOKsyDCUTwjhJYTYK9Te5+2AkyF5hRAOQohQIUQn/c/WQohrQogBebTjTSHEGX1v6PWcAmshhJEQ4n/6XuZ7QoiFQghbfVq2oYb6HsjX9N9b6Htvo4QQF4F6ebQlv7K+FkIs19cfo+859c3l3DzrFUJ8rr/eGCHERSFEV/3xKsBMwE/fm/pQf7yDUEcCPBJC3BI59MRnKDvTENYMPbTG+p8dhBDzhRB39O1boz9uL4TYIISI0B/fIIRwz1DOHiHEd0KIg6gPjsrlUHdVIcQOIcQDofaAd8uQVlJf5iMhxBHAK0NaBZGlF18IcUA8pWG9kiS9mGTAJklSOiFEKaAijwOTasDZtHRFUeKA6/rjeVIUJRxYCrwrhNAIIfxQe3oKO8Qpa1uuA8n69hbWD0Id5nhQ5DLnxlBCiL+EEPHAZSAM2JQh7aEQokmWLCH6m/r5Qu0py2io/kbvZMYbvRzKyvX9EULYA64Z0/Xf5/ve5WMRaiAIag/YwicsL6tuQCywAtiaoa5shBANgeqoPbMFJoRwQ+3N/Rb1AcNnwCohREn9Kf8CJ1EDtUlk6PHLK6+iKA+AQcAcIYQzMAU4oyhKjq+VEKI+6us4GrADmgHBOZz6jv6rJWqQYA38aeDlfoX6EKE80Jbcey8N1RlYhtredXm0I796rwNNUXvhJwKLhRAu+gcMHwGH9T2wdvrz41A/E3ZAB2CIEKJLIa9hEWCJ+juR9j6Bem80H/XvVRkgIYfr6w8MBkoAIRkThPoQaTvqe+qM+rBpthAirRd/BupDqdL6MgYVsv2SJL0iZMAmSRKQPoRxCfCPoiiX9Yetgegsp0aj3qQYYinwJZCEOgRugqIotwrZxCdtS1ZjUW963YDZwHohRKF7xRRFGapvS1PUYYpJGdLsFEVJC1Tvo/YylEXtySyB+rqn+QN1GKMzai/iAiFE41zKyus1sc7wc9a0J7EY6K3/vLyt/zkrV31gmfHLysDyBwL++mGg/2aoK6P7QogE4DDqsNg1hbsU+gGbFEXZpO9h3o46rLS9EKIM6vv0haIoSYqi7APWG5IXQFGUbahB507UwOLDPNrxHjBPUZTt+rJuZ/gdzKgv8JuiKDf0vczjgLfTeozy0RP4TlGUB/rfwT8MyJOXA/pr16IGPrkNdc6zXn3v8B39dfsDgUD93CpVFGWPoigB+vPPof6NaV7QxgshXIB2wEf6HugURVH26uuIVBRllaIo8YqixADf5VDHAn3PeqqiKClZ0joDVxVFWahPP4n6Ge2u/yx3Qf1cxeuvYVFB2y9J0qtFBmySJCHUORiLUHusPsmQFIs6rDEjG3Iespi1zMqo86kGAKaoT7HHCCE66NMzLkDR1IBm5toWIUTfDGVtNqAsFEU5qihKjP5m/B/gIPqb7cJSFEWrD6bcgSG5nBOrKMoJ/Y1cOOrr3UYIYaNPP6W/YUxVFGUTajD3Vi5V5vX+xGb4OWtaoSmKchO1R+t7IDCXAPyOPrDM+BWXX9lCCA/U3qO0AHYtYI4a8GTkhBqQfga0ALIGdIYqC/TIGFgCTQAX1N7JqCztDjEwb5rZqD2A8xVFidRfY5kMn9W098gDtacpP65Z2hCCOg+1lIF5M75XIbmdaKC7Gb6PB8xzCRzzrFcIMUCoQ0HTXsPqZBl6muX8BkJdXCdCCBGN2guX6/l58AAeKIoSlUMdlkKIWUIdevoI2AfYicyL/+T14Kks0DjLZ6MX6mejFKDh6b4XkiS95GTAJkmvOCGEAOai3kh0y/K0+AIZnpzre0nKY9i8sOrAFUVRtuqfhl9BHULWDtJXkktbcGK/AeVlbUs5wAz1SfaSDGW1y7WEvCmAKGTerIwxfA5b2nyV3OrOq125vj/6G9EwMvd81OIpzOlDHeo1iqc/HLI/6v+l9UKIu8AN1IAt27BIfXA8GUgEhhayvluocywzBpZWiqL8iPra2WfpGSxjYF70N/ezUF+jIUKICvp238zwWbXOUJYhn5c7qMFAxvakAuGoQwUt0xL09ZfMcG4YapCS07VklV9ZBZFrvUJd0XQO6kMLR/2wx/M8/rzntM3Ev6hDMD0URbFFneeW2+9HputAHYKY5hbgIISwI7tRQCWggaIuJNQsrckZzslrC4xbwM4snw1rRVE+QX2vdOT+XsRB+lzOnNotSdIrSAZskiTNAKoAnRRFSciSthqoLoToJoQwRx3eeC5tuJZ+bpo5aoBiJIQwzzB87TTqohGthKo80JHMc6oyEUKY6MszAoz15aU91V4CdBJCNNXfRH8D/KcfspRTWUJflqn+Z3MhhJn+ezshRFv9MWMhRF/Um7KthSjLWQjxtlAXltAIIdoCvYFduZTVQAhRSagLSDiiDhHboyhKtD69u74sIyFEG9Shd+tyecnyfH9Qg4X/CXURhcrAB8CCXMoqCH+gDerKmU/TANR5TD4ZvroBHfSvVU5+RO25NS9EfYtRP1Nt0z7LQl1ww11RlBDUIY4ThRCmQp032MmQvPr0tG0KBgG/AgtF7tszzEWd69la/767iccLx2S0FPg/oS6GYo3ay+mvKEoq6kqd5kJdlMME+B/qA400y4Fx+s+COzAsj9clv7IKIq96rVADnwgAIcS7qA960oQD7kII0wzHSqD2jCUKde5fnzzqPgM00/dq2qI
      "text/plain": [
       "<Figure size 1440x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,\n",
    "                             standard_parallels=[35])\n",
    "\n",
    "# Change the DPI of the figure\n",
    "plt.rcParams['savefig.dpi'] = 255\n",
    "\n",
    "# Winds, temps, dewpoint, station id\n",
    "custom_layout = StationPlotLayout()\n",
    "custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')\n",
    "custom_layout.add_value('NW', 'air_temperature', fmt='.0f', units='degF', color='darkred')\n",
    "custom_layout.add_value('SW', 'dew_point_temperature', fmt='.0f', units='degF', color='darkgreen')\n",
    "custom_layout.add_value('E', 'precipitation', fmt='0.1f', units='inch', color='blue')\n",
    "\n",
    "# Create the figure\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=proj)\n",
    "\n",
    "# Add various map elements\n",
    "ax.add_feature(cfeature.LAND)\n",
    "ax.add_feature(cfeature.OCEAN)\n",
    "ax.add_feature(cfeature.LAKES)\n",
    "ax.add_feature(cfeature.COASTLINE)\n",
    "ax.add_feature(cfeature.STATES)\n",
    "ax.add_feature(cfeature.BORDERS, linewidth=2)\n",
    "\n",
    "# Set plot bounds\n",
    "ax.set_extent((-118, -73, 23, 50))\n",
    "ax.set_title(str(ob.getDataTime()) + \" | METAR | \" + edexServer)\n",
    "\n",
    "stationplot = StationPlot(ax, data['longitude'], data['latitude'], clip_on=True,\n",
    "                          transform=ccrs.PlateCarree(), fontsize=10)\n",
    "stationplot.plot_text((2, 0), data['stid'])\n",
    "custom_layout.plot(stationplot, data)\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
}