python-awips/examples/notebooks/Watch_and_Warning_Polygons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example uses matplotlib, cartopy, shapely, and python-awips to plot watch and warning polygons requested from a real-time AWIPS EDEX server.\n",
    "\n",
    "First, set up our imports and define functions to be used later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "from awips.tables import vtec\n",
    "from datetime import datetime\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "from cartopy.feature import ShapelyFeature,NaturalEarthFeature\n",
    "from shapely.geometry import MultiPolygon,Polygon\n",
    "\n",
    "def warning_color(phensig):\n",
    "    return vtec[phensig]['color']\n",
    "\n",
    "def make_map(bbox, projection=ccrs.PlateCarree()):\n",
    "    fig, ax = plt.subplots(figsize=(20,12),\n",
    "            subplot_kw=dict(projection=projection))\n",
    "    ax.set_extent(bbox)\n",
    "    gl = ax.gridlines(draw_labels=True)\n",
    "    gl.xlabels_top = gl.ylabels_right = False\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a request for the \"warning\" data type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 82 records\n",
      "{'phensig': array([], dtype=float64)}\n"
     ]
    }
   ],
   "source": [
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"warning\")\n",
    "request.setParameters('phensig')\n",
    "times = DataAccessLayer.getAvailableTimes(request)\n",
    "\n",
    "# Get records for last 50 available times\n",
    "response = DataAccessLayer.getGeometryData(request, times[-50:-1])\n",
    "print(\"Using \" + str(len(response)) + \" records\")\n",
    "\n",
    "# Each record will have a numpy array the length of the number of \"parameters\"\n",
    "# Default is 1 (request.setParameters('phensig'))\n",
    "parameters = {}\n",
    "for x in request.getParameters():\n",
    "    parameters[x] = np.array([])\n",
    "print(parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now loop through each record and plot it as either Polygon or MultiPolygon, with appropriate colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIwAAAIjCAYAAABhzICuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XdYFNfbxvHv7KIgqChSlCgW7F2wd1FREbuCvffYe+81lsQYu2jsBXvvYuyaxK4x9i6IgiAKUub9w8RffFWkLMwCz+e6uNTZmTP3JrvDzDPnnFFUVUUIIYQQQgghhBBCiH/ptA4ghBBCCCGEEEIIIYyLFIyEEEIIIYQQQgghxCekYCSEEEIIIYQQQgghPiEFIyGEEEIIIYQQQgjxCSkYCSGEEEIIIYQQQohPSMFICCGEEEIIIYQQQnxCCkZCCCGEEEIIIYQQ4hNSMBJCCCGEEEIIIYQQn5CCkRBCCCGEEEIIIYT4hInWAb7EwSGr+ujRE61jCCGEEEIIIYQQQiQnD1RVzRGTFRVVVRM4S+wpiqIGBNzSOoZIJH/8cQdnZ8cEaXv27IU8ePCIOXMmJ0j7QgDkz1+eX39dTtmy+bSOIpKwCRNmYW6ehkGDemqy/4Q8Fgu4c+c+Q4dO4MmTZ8yYMY6KFctoHSnZkc+wiK3IyEjc3FrSsmUL+vYdpHUcfHx8qFq1qkHbvHTpT+rWdefyZR90OhlcIhJWcjsO16nTgnbtPGjevNFnr4WGhrF372FWrNiApWU6ChTIy/Ll61i3bhFOTkU1SBtzGTPmQVVVJSbrGmUPIyEMpX17T8qWrUOaNGbY2tpQunQJOUkXQhglOztr7tx5oHUMYWBv377jxx8XsmzZOvr370a3bm1JlSqV1rGEEIBer+fnnydTt25LGjRoRI4cyedC919FihTHwsKcP/+8TMmSxbWOI0SScffuA+7de0DTpvW++PqoUVO5dOkqzZo1oHXrppibp6Fo0UK0bNmds2f3YWmZPpETf5uqqjx+/CxW20iZWSRrVlYZ2blzNZkz2/L6dRA9egzhl1+8tI4lhBCfsbW1wc/vhdYxhAEdPHiMcuXcuHPnAceP76BXr05SLBLCyOTLl5vu3dvRtWsXoqKitI5jcDqdjsaNG7J9+z6towiRpFy6dI2SJYtjYvJ5HxtVVdm79zALFsyga9c2mJunAcDMLDWhoWFGdyxRVZXx42eSI4czrq7NYrWtFIxEspcvX2769evG+PFD2L17DTNmzOP5cz+tYwkhxCdsba15/lwKRslJr17DmDZtFMuW/YS9fWat4wghvqJPny48ffqMlSuXaR0lQXh4NGf79n0Y41QkQhirFy9ekjmz7Rdfe/DgEZGRkTg65vhk+dSpc5g/fzoZM2ZIhIQxd+nSNTZt2snp03u4ceNkrLaVgpFIURwcstKvXzdcXBpz754M/RBCGA87OxtevPDXOoYwoIiISEqXLqF1DCHEN6ROnZq5c6cwZMgwnj9/qnUcgytWzAlT09RcuHBF6yhCJBn29nYcO3aKQ4d+IygomHfvQomKiuLq1b+YOvVnqlQpj6L8bxqghw+f8Ndft6lQobSGqb8sMPA1jo454nTzSuYwEilO//7dCAsLY/bshcydO1XrOCIZ+O8vCyHiytbWGj8/KRglJ3I3X4iko0SJInh6NqR37554e2/TOo5B/Tssbdu2vUY/Ga8QxqJu3ZoEBLxm9uwFXLx4lYiISFKlMsHaOhPu7jXp27frx3VVVaVly+6MGNHXKOcuCgl5i4WFeZy2lR5GIkW6ePEqWbPaax1DCCE+SpcuLZGRkbx5E6J1FGFAUlAWIukYPrwvf/xxgS1bNmodxeCaNWvOjh37pZAtRAwpikKbNs3Ys2cdT59ewdf3GpcuHeXixSNMnjwCW1vrj+uqqsqtW3do1665hom/LiTkLWnSmMVpWykYiRSpbVsPvLzWcOvWXa2jCCEE8OHExNbWWoalJSNyYSZE0mJunoY5cybTp09fAgNfaR3HoJycSqKqKteu3dQ6ihBJkqIoWFtn+uKNIJ1Oh4NDVh49eqJBsm9LnTo1YWFhQOzPTaRgJFIkN7cafP99R/r0GUlQUPDH5ffuPZAhIUIIzdja2uDrK8eg5EJVVelhJEQSU6lSWWrUqMLAgf21jmJQN29e582bEL77TibgFyIhZM+elfv3H2kd44ucnIpw8uR5goKCGTx4fKy2lYKRSLF69erEgwePOH/+IgC3bt3FxaUJpUvXYuDAsVy5cl3uDgshEpWdnTV+fvKkNCGE0NL48UPYt28/hw/v1zqKwYwZM4oePdob3dObhEgusmfPxoMHxlkwcnDIStWq5Rk/fiZeXmtita0UjESKpdfriYyMZNeuAzRo0BYXl8YMHvw9587tJ2NGS9q370PBghUZNmwiUVFRWscVQqQAtrY2PH8uBaPkQnoYCZE0WVqmZ8aM8XTt2o23b5P+vHK//36O3347Sffu7bSOIkSyFBUVhb//S4KC3mgd5at69+7M/v1H2blzday2k4KRSNHWr19M9uxZ6dWrE5cuHaVnzw7Y2lozatQA/vjjEFOnjmLRopW8fh2kdVQhRApgZydzGCU3UjASImlyc6tOsWKFGDVqmNZR4m3kyOEMGNCdtGkttI4iRLI0ZcpP+Pr606VLa62jfFXx4oUJDHxN4cL5Y7WdFIxEilaiRBH69etGzZpVsLLK+MlrFy9eZfDgcaxc+Yt03xVCJAobG2uZRy0ZkWHNQiRt06aNZvXqdZw5c1LrKHHm43OYGzdu0r69cT69SYjk4N69h7Rt2wxLy/R07NiPwoUrs2PHPqM6D/D3f0Xq1KnJkMEyVttJwUiIr9DpFFRVJTj4De/fv9c6jhAiBbCzs5GCUbIjPYyESKpsba2ZMGEYXbp0TpLnglFRUYwcOYJhw3pjamqqdRwhkq3q1SuxadNO3rwJ4eBBHyZPHs7MmfOpXr0Je/YcilEbqqpy9uyf7N9/lMuXr33yYCZDuHPnPo6O2WO9nYlBUwiRjBQtWghvby/GjJnOuHEzaNmyCbVqVUNRFCIjIwgNfc+FC1do0qQuOXPG/ssnhBD/n62t9DBKXoznzqIQIm48PRuwadNOhg0byIwZP6HX67WOFGO7d2/H3/8lHh4NtI4iRLLWpIk7CxeuoFq1RpQqVZwGDepQr14t9u07Qp8+I8iVKzv58+eJto0lS1bx889LyZ8/D8+e+XL//kOsrDJStqwzlSuXo2HDOqRLlzZO+aKioti9+yAFC+aL9bZSMBIiGiVKFGHnztXcunWXFSs2MGbMNHQ6PXq9DhMTE6ysMrB16x4OH96MmZncuRFCxI+trQ2+vjLpdXIiUxgJkbQpisLPP0+mffveVKpUjsWLl1C4cDGtY31TZGQkI0eOYuTIfpiYyCWfEAnJ1NSUnTvXsHLlBlxcKgGg0+lwc6vB48fP+P77YcyfP518+XIDcPz4Gc6e/ZMcObJRvXolrl//mxkz5nHwoDc5cjgAH3oc3b59j7Nn/2DfvqOMGTOdVq2aMGhQz1gNKzty5DgjR07B3NycNWsWxPq9ydFDiBjIkycXkyYN/2y5qqq0b9+bESMmM2PG2CR110kIYXxsbT9Mei1P10oejGjqAiFEPNjbZ2bv3vV4ea2lSpVq9OzZlVGjxhv1MK/161djYqKnXr1aWkcRIkWwtExH796dP1veqVNLgoPfUL9+G9KnT0eRIgXx8TlJixaNuHDhCr17Dyd9+nQsWDDjY7EIPhSr8+TJRZ48uWjduhlPnz5n5sz5lCrlyvTpY2jcuO43M0VERNCt2yB++mkSbm414nRuKXMYCREPiqLw00+TuHXrLg0btuPJk2daRxJCJGFmZqaYmZkRGPha6yjCQKTwJ0TyoNfr6dq1DceObef8+XOUKFGU48d9tI71ReHh4YwbN57RowfIMUgIjen1egYO7MGNGydZuvRHateuxubNy5g8eQRr1izg6dMr3Lx5mho1Kkfbjr19ZmbPnsDmzcsYM2Y6v/zi9c19nzhxlmzZvqNu3ZpxPhZIwUi
      "text/plain": [
       "<Figure size 1440x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "bbox=[-127,-64,24,49]\n",
    "fig, ax = make_map(bbox=bbox)\n",
    "\n",
    "siteids=np.array([])\n",
    "periods=np.array([])\n",
    "reftimes=np.array([])\n",
    "\n",
    "for ob in response:\n",
    "    \n",
    "    poly = ob.getGeometry()\n",
    "    site = ob.getLocationName()\n",
    "    pd   = ob.getDataTime().getValidPeriod()\n",
    "    ref  = ob.getDataTime().getRefTime()\n",
    "    \n",
    "    # do not plot if phensig is blank (SPS)\n",
    "    if ob.getString(b'phensig'):\n",
    "        \n",
    "        phensigString = ob.getString(b'phensig')\n",
    "        \n",
    "        siteids = np.append(siteids,site)\n",
    "        periods = np.append(periods,pd)\n",
    "        reftimes = np.append(reftimes,ref)\n",
    "\n",
    "        for parm in parameters:\n",
    "            byte_str = str.encode(parm)\n",
    "            parameters[parm] = np.append(parameters[parm],ob.getString(byte_str))\n",
    "\n",
    "        if poly.geom_type == 'MultiPolygon':\n",
    "            geometries = np.array([])\n",
    "            geometries = np.append(geometries,MultiPolygon(poly))\n",
    "            geom_count = \", \" + str(len(geometries)) +\" geometries\"\n",
    "        else:\n",
    "            geometries = np.array([])\n",
    "            geometries = np.append(geometries,Polygon(poly))\n",
    "            geom_count=\"\"\n",
    "\n",
    "        for geom in geometries:\n",
    "            bounds = Polygon(geom)\n",
    "            intersection = bounds.intersection\n",
    "            geoms = (intersection(geom)\n",
    "                 for geom in geometries\n",
    "                 if bounds.intersects(geom))\n",
    "    \n",
    "        #print(vtec[phensigString]['hdln'] \n",
    "        #  + \" (\" + phensigString + \") issued at \" + str(ref)\n",
    "        #  + \" (\"+str(poly.geom_type) + geom_count + \")\")\n",
    "        \n",
    "        color = warning_color(phensigString.decode('UTF-8'))\n",
    "        shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), \n",
    "                        facecolor=color, edgecolor=color)\n",
    "        ax.add_feature(shape_feature)\n",
    "            \n",
    "states_provinces = cfeature.NaturalEarthFeature(\n",
    "        category='cultural',\n",
    "        name='admin_1_states_provinces_lines',\n",
    "        scale='50m',\n",
    "        facecolor='none')\n",
    "political_boundaries = cfeature.NaturalEarthFeature(category='cultural',\n",
    "                               name='admin_0_boundary_lines_land',\n",
    "                               scale='50m', facecolor='none')\n",
    "ax.add_feature(cfeature.LAND)\n",
    "ax.add_feature(cfeature.COASTLINE)\n",
    "ax.add_feature(states_provinces, edgecolor='black')\n",
    "ax.add_feature(political_boundaries, edgecolor='black')\n",
    "\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
}