python-awips/examples/notebooks/Watch_and_Warning_Polygons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Request for Data Type Warning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "433 records\n"
     ]
    }
   ],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "\n",
    "DataAccessLayer.changeEDEXHost('edex-el7.westus.cloudapp.azure.com')\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype(\"warning\")\n",
    "request.setParameters('act','countyheader', 'endtime','etn','floodbegin', \n",
    "                      'floodcrest', 'floodend','floodrecordstatus', \n",
    "                      'floodseverity', 'id', 'immediatecause', 'inserttime', \n",
    "                      'issuetime', 'loc', 'locationid', 'motdir', 'motspd', \n",
    "                      'overviewtext', 'phen', 'phensig', 'pil', 'productclass',\n",
    "                      'purgetime', 'rawmessage', 'seg', 'segtext', 'sig', \n",
    "                      'starttime', 'ugczones', 'vtecstr', 'wmoid', 'xxxid')\n",
    "response = DataAccessLayer.getGeometryData(request)\n",
    "print(str(len(response)) + \" records\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a dictionary of lists for each parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 433 geometries.\n",
      "['act', 'countyheader', 'endtime', 'etn', 'floodbegin', 'floodcrest', 'floodend', 'floodrecordstatus', 'floodseverity', 'id', 'immediatecause', 'inserttime', 'issuetime', 'loc', 'locationid', 'motdir', 'motspd', 'overviewtext', 'phen', 'phensig', 'pil', 'productclass', 'purgetime', 'rawmessage', 'seg', 'segtext', 'sig', 'starttime', 'ugczones', 'vtecstr', 'wmoid', 'xxxid']\n"
     ]
    }
   ],
   "source": [
    "parameters = {}\n",
    "for x in request.getParameters():\n",
    "    parameters[x] = np.array([])\n",
    "geometries=np.array([])\n",
    "siteid=np.array([])\n",
    "period=np.array([])\n",
    "reftime=np.array([])\n",
    "\n",
    "for ob in response:\n",
    "    for parm in parameters:\n",
    "        parameters[parm] = np.append(parameters[parm],ob.getString(parm))\n",
    "    geometries = np.append(geometries,ob.getGeometry())\n",
    "    siteid = np.append(siteid,str(ob.getLocationName()))\n",
    "    period = np.append(period,ob.getDataTime().getValidPeriod())\n",
    "    reftime = np.append(reftime,ob.getDataTime().getRefTime())\n",
    "\n",
    "print(\"Using \" + str(len(geometries)) + \" geometries.\")\n",
    "print(request.getParameters())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Print out basic info for each warning**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FL.W\n",
      "KMKX: WGUS83 KMKX 042104\n",
      "issued 2017-02-04 14:04:00\n",
      "start 2017-02-04 14:04:00\n",
      "end 2017-02-07 23:00:00\n",
      "purge 2017-02-05 02:04:00\n",
      "POLYGON ((-89.7 43.51, -89.66 43.65, -89.37 43.59, -89.39 43.39, -89.65000000000001 43.44, -89.56999999999999 43.5, -89.7 43.51))\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i, val in enumerate(geometries):\n",
    "    issue = datetime.fromtimestamp(\n",
    "        float(int(parameters['issuetime'][i])/1000)\n",
    "    ).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    purge = datetime.fromtimestamp(\n",
    "        float(int(parameters['purgetime'][i])/1000)\n",
    "    ).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    start = datetime.fromtimestamp(\n",
    "        float(int(parameters['starttime'][i])/1000)\n",
    "    ).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    end = datetime.fromtimestamp(\n",
    "        float(int(parameters['endtime'][i])/1000)\n",
    "    ).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    \n",
    "if parameters['phensig'][i]:\n",
    "    print(parameters['phensig'][i])  \n",
    "print(siteid[i] +\": \"+ str(parameters['wmoid'][i]) \\\n",
    "      + \"\\nissued \" + issue \\\n",
    "      + \"\\nstart \" + start \\\n",
    "      + \"\\nend \" + end \\\n",
    "      + \"\\npurge \" + purge \\\n",
    "      + \"\\n\" + str(geometries[i]) + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Warning Polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIcAAAIfCAYAAADqlLyPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8Tff/wPHXDULESEKNNghirxgxE0JsIYIatUe1+Cql\nSq3WqFGj9qitoZUUsSmJCEJDao+IxggSZBJC1v390Sa/IDv35q738/HIg3vvuee8z733vO/nvO/n\n8zkKpVKJEEIIIYQQQgghhDBMRpoOQAghhBBCCCGEEEJojhSHhBBCCCGEEEIIIQyYFIeEEEIIIYQQ\nQgghDJgUh4QQQgghhBBCCCEMmBSHhBBCCCGEEEIIIQxYfk0H8D6FQiGXTxNCCCGEEEIIIYRQMaVS\nqUjrfq0rDgEolVIfErmnUPz7mZfPk1AX+YwJVZHPkmGQ91mIvBcdHU2lSpX4/fffadeuHQBhYWEs\nXbqUvXv38vLlSzp16kTTpk1p2rQpNWrUwMgoe4MrXr58yWeffcbr169xd3fHwsJCHbuSZaNGjeLa\ntWts27aNypUrazQWIXRR8vc1pP2d/f73uS59v6fet/fJsDKhUd7e3mpdf+oDNqMDQQihvdSdJ7SB\nUqnUiQaFyJ7k757UjUZ5n9XDEPKEyJnixYvj5ubG8OHDcXJy4ty5c5QsWZJ58+Zx69YtTpw4Qd26\ndTl58iTOzs6Ym5vTrl07pk+fzsGDB3n+/Hmm2yhatCgeHh7Url0bZ2dnjR/nq1atokePHjRp0oQD\nBw5oNBZtITlCZEfyMZzZsaxv3+9SHBJ6Sx8OUKHd9OWLQAghhNBnjo6O3Llzh65du9KrVy+mTZtG\nQkICANWrV2fs2LG4urpy9+5d7t69y/jx41EoFKxYsYIqVapgbW3NyJEj2bNnD9HR0WluI1++fCxd\nupTY2FhcXV3zcvfSjGXChAksWrSIxYsXazQWIXRVRm385HOA1Mtoa0eE7HSSUGjbiY1CoVAGBwdj\naWmp6VCEEEKIPKFL3ZGFEEKXPXv2jP79+2NqaoqbmxvGxsYZLp+UlJTSw+jIkSOcPXuW+vXr065d\nO5o3b06tWrUoXbp0Sh5ft24dvr6+bN++PS92J0OOjo4MHDiQIUOGaDoUIXRKRu2ytB7LbBiaJqU1\nBC69OYe0sji0YMECJk+erOlQhB7S5gNX6B45oReqIp8l/SLvpxDaLS4ujr59+5KYmIi7u3umBaLU\nXr9+jY+PD56enpw/f56AgABiYmKwsrKiVKlSXL16lbVr19KnTx817kHWlC9fHk9PT6pUqaLpUITQ\nKam/x3V5fqG0ZFQc0sphZQ0bNtR0CCKPyPhfIURmDCFPyBBFIXLHEPKEUB1jY2N27dpFvnz56NWr\nF2/fvs3ycwsXLkzHjh1ZtGgRp0+f5tmzZzx79gw3NzemT5/O5cuXtaIwBFC3bl2uXr2q6TC0guQI\nkR2G2i7TyuJQ27ZtNR2C0FNpjQ8VIqfksySESIvkBiG0X4ECBdi1axeFCxemRYsWBAUF5XhdRYoU\noXbt2rRt25by5curMMrccXZ2Zvbs2dy/f1/ToQih01J/p6c115C2zjcEejDnkLbFJPSPrncHFELo\nF8lJQgihGUqlkpUrVzJ37ly2bdtGp06dNB2SyiiVSn7++WcWLFjAhg0bcHZ21nRIQugdbW/DpTNE\nTnfmHNK2mIT+SW/saOr7hMiMtn8ZCN0hnyX9Iu+nELrH19cXFxcXli5dSv/+/TUdjkqdP3+eTz/9\nlFmzZjFs2DBNhyOE0JDMikNaOaxMGA5Njv+VRrsQusEQ5gmQYUhC5I4h5AmhXs2bN8fLy4spU6aw\nevVqTYejUk2bNuXw4cNMnTqVpKQkTYejEZIjhKHLytAyKQ4Jg5TcrS71r7upT860feyo0A5yQi+E\nSIvkBiF0U61atTh9+jTTp08nJCRE0+GoVJ06dbCwsODSpUuaDkUIvZJ83qjN545Z7dGslcUhfUvG\nIn0ODg6aDkEIoeUMJU9oe8NCCG1mKHlCqN+bN28oWLAgwcHBmg5F5aysrHj69Kmmw9AIyRFCXfTp\nxyCtLA7VqlVLLxOy0C7yy67ILTmZF0KkRXKDELolISEBDw8Pevfujb29PfPmzaNx48aaDkvlKleu\nzGeffcbFixc1HYoQeiWteWy1SVbPe7WyOBQZGcn8+fPlxN0AyPhfIURmDCVPSMFaiJwzlDwhVC8m\nJoYCBQrg4uKCu7s7ERERTJkyhXr16rF48WJev36t6RBVpmDBgkRHR3PlyhVNh5LnJEcIdUnvB6Gc\n/FCkrh+XsrrO/CrfsgrUrl2btWvXEhoaiqurK4ULF9Z0SMLApD5Be/9KZnIVGpFMPgNCiLRIbhBC\nd5iamrJ06VI++eQTKlWqRLly5QAICAhgyZIleHh44OHhQcmSJTUcae5s2bKFQ4cOERYWRokSJTQd\njhB6Rx3f/ao871QoFJmuRysvZf/mzRsGDBjAH3/8QYECBQgLC6NYsWKaDk0YKLnsvRAiL0jhWQgh\ntEtSUhLTpk3jl19+YdiwYXz++edUqVJFY0NHlEol165d46+//sLa2hoHB4csxXLw4EGGDRvGyZMn\nqVWrVh5EKoRILXVhJr1zy7Taf6pqG6Zez3+xpJk4tLI4lNy1ftOmTXz++ecMGTKELVu2aDo0IT4g\nJ3OGTd5/oUryedIf8l4KoV8ePnzIzz//jJubG8bGxnTs2JFhw4Zha2ur9m1HR0dz4sQJjhw5wtGj\nRylYsCBlypTB19eX58+fZ9qjac+ePXz55ZccPHhQL+dREkLbZTYKJatthux2UEhvvRkVh7RyziH4\nN+gRI0bw/PlzFi1apOlwhJrI+F8hRGYMJU/InENC5Jyh5AmhGeXLl+fnn3/m0aNHHDx4ECsrK3r1\n6kWbNm04duyY2nJ3REQEZmZmTJs2jTp16uDl5UVAQADx8fEsX748w8LQ7du3cXZ2ZsKECRw+fNjg\nC0OSI4SmJOeH1P+mzhna1P7T2uJQspIlS+r8GF+hv94/mLV1hnqhHtqUzIUQ2kNygxD6SaFQUKtW\nLSZPnszdu3cZPnw448ePp0ePHoSGhqp8eyYmJkyYMIGIiAhWrVrF8OHDsba2pnTp0owePTrd53l6\nemJnZ4e9vT23b9+mUaNGKo9NCJF1qmgTJLctsnpltJy0RbR2WJkQukbmIhJC5IYMRRJCCN3z9u1b\nZs+ezcaNG1myZAn9+/dX+Y+FCQkJBAYG8vz5cywsLKhVq1a621AqlVhbW7N27Vrat2+v0jiEENoh\nN+edOjnnkBC6KK1Z4OWET3/JeytUSQrM+kNygxCGx9/fn6FDh1KhQgXWr1/Pxx9/rJE4jh49yqhR\nowgKCpIe7ULosYyuPpbZJNc6N+eQMAz6Nv5XTgSEUD19yxPpeb+7sBAi6wwlTwjt1bBhQy5evEiD\nBg1o2LAhfn5+eR7D/Pnz6d+/P1u2bJHC0HskRwh9kzy5dfKfKuRXyVpUbMSIESgUCqKjo4mKikqp\nwBsZSS1L6Ib3q7Vysqef5H0VQqRFcoMQhsnY2JhZs2bRsGFDunfvjr+/P2XLls2z7V+/fp1u3bph\nZ2eXZ9sUQmhOqkvTp/lYdmnlsLL377OysiIwMJD8+bWyliXEBzK6dGBa9wshBEiOEEIIfTF79mxO\nnjyJp6enyn/gDg8Px9zc/IP1RkVF0adPHwB27dqFmZmZSrcrhNB+mbUldW5YWXh4OElJSSnd6+/d\nuyeFIaFTZGiIYVBlN04hhP6Q3CCEmDZtGm/evOGXX35RyfqePn3Kzz//TP369SlXrhw9evT4YBkz\nMzMOHTpEtWrVaNasGXfv3lXJtoUQ2kuVbQ6tLA5ZWFhIo8pAGNr43/cvPyifcyEyZ0h5QgrLQuSM\nIeUJoRvy5cvHpk2bmDFjBsHBwblal0KhoEyZMly9epUlS5YQHh7OxYsXmTRpEr///jsREREpy+bP\nn58VK1Ywbtw47OzsuHL
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106633c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cartopy.feature as cfeat\n",
    "from cartopy.feature import ShapelyFeature,NaturalEarthFeature\n",
    "from shapely.geometry import Polygon\n",
    "from awips.tables import vtec\n",
    "\n",
    "def warning_color(phensig):\n",
    "    if vtec[phensig]['color']:\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",
    "    ax.coastlines(resolution='50m')\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\n",
    "\n",
    "%matplotlib inline\n",
    "bbox=[-127,-64,24,49]\n",
    "fig, ax = make_map(bbox=bbox)\n",
    "for i, geom in enumerate(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",
    "    # For now do not plot SPS\n",
    "    if (parameters['phensig'][i]):\n",
    "        color = warning_color(parameters['phensig'][i].astype('U13'))\n",
    "        shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), \n",
    "                        facecolor='none', edgecolor=color)\n",
    "        ax.add_feature(shape_feature)\n",
    "\n",
    "\n",
    "political_boundaries = NaturalEarthFeature(category='cultural',\n",
    "                               name='admin_0_boundary_lines_land',\n",
    "                               scale='50m', facecolor='none')\n",
    "states = NaturalEarthFeature(category='cultural',\n",
    "                               name='admin_1_states_provinces_lines',\n",
    "                               scale='50m', facecolor='none')\n",
    "ax.add_feature(political_boundaries, linestyle='-')\n",
    "ax.add_feature(states, linestyle=':')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Warning Polygons Plotted Over Level 3 Radar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed 1 grid.\n"
     ]
    }
   ],
   "source": [
    "from awips import ThriftClient, RadarCommon\n",
    "from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange\n",
    "from dynamicserialize.dstypes.com.raytheon.uf.common.dataplugin.radar.request import GetRadarDataRecordRequest\n",
    "from datetime import timedelta\n",
    "from numpy import ma\n",
    "from metpy.plots import ctables\n",
    "\n",
    "site = 'khnx'\n",
    "request = DataAccessLayer.newDataRequest()\n",
    "request.setDatatype('radar')\n",
    "request.setLocationNames(site)\n",
    "\n",
    "# Get latest time for site\n",
    "datatimes = DataAccessLayer.getAvailableTimes(request)\n",
    "dateTimeStr = str(datatimes[-1])\n",
    "buffer = 60\n",
    "dateTime = datetime.strptime(dateTimeStr, '%Y-%m-%d %H:%M:%S')\n",
    "beginRange = dateTime - timedelta(0, buffer)\n",
    "endRange = dateTime + timedelta(0, buffer)\n",
    "timerange = TimeRange(beginRange, endRange)\n",
    "\n",
    "# GetRadarDataRecordRequest to query site with timerange\n",
    "client = ThriftClient.ThriftClient('edex-cloud.unidata.ucar.edu')\n",
    "request = GetRadarDataRecordRequest()\n",
    "request.setTimeRange(timerange)\n",
    "request.setRadarId(site)\n",
    "\n",
    "nexrad = {}\n",
    "nexrad[\"N0Q\"] = {\n",
    "    'id': 94,\n",
    "    'unit':'dBZ',\n",
    "    'name':'0.5 deg Base Reflectivity',\n",
    "    'ctable': ['NWSStormClearReflectivity',-20., 0.5],\n",
    "    'res': 1000.,\n",
    "    'elev': '0.5'\n",
    "}\n",
    "\n",
    "grids = []\n",
    "for code in nexrad:\n",
    "    request.setProductCode(nexrad[code]['id'])\n",
    "    request.setPrimaryElevationAngle(nexrad[code]['elev'])\n",
    "    response = client.sendRequest(request)\n",
    "\n",
    "    if response.getData():\n",
    "        for record in response.getData():\n",
    "            # Get record hdf5 data\n",
    "            idra = record.getHdf5Data()\n",
    "            rdat,azdat,depVals,threshVals = RadarCommon.get_hdf5_data(idra)\n",
    "            dim = rdat.getDimension()\n",
    "            lat,lon = float(record.getLatitude()),float(record.getLongitude())\n",
    "            radials,rangeGates = rdat.getSizes()\n",
    "\n",
    "            # Convert raw byte to pixel value\n",
    "            rawValue=np.array(rdat.getByteData())\n",
    "            array = []\n",
    "            for rec in rawValue:\n",
    "                if rec<0:\n",
    "                    rec+=256\n",
    "                array.append(rec)\n",
    "\n",
    "            if azdat:\n",
    "                azVals = azdat.getFloatData()\n",
    "                az = np.array(RadarCommon.encode_radial(azVals))\n",
    "                dattyp = RadarCommon.get_data_type(azdat)\n",
    "                az = np.append(az,az[-1])\n",
    "\n",
    "            header = RadarCommon.get_header(record, format, rangeGates, radials, azdat, 'description')\n",
    "            rng = np.linspace(0, rangeGates, rangeGates + 1)\n",
    "\n",
    "            # Convert az/range to a lat/lon\n",
    "            from pyproj import Geod\n",
    "            g = Geod(ellps='clrk66')\n",
    "            center_lat = np.ones([len(az),len(rng)])*lat\n",
    "            center_lon = np.ones([len(az),len(rng)])*lon\n",
    "            az2D = np.ones_like(center_lat)*az[:,None]\n",
    "            rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None])*nexrad[code]['res']\n",
    "            lons,lats,back=g.fwd(center_lon,center_lat,az2D,rng2D)\n",
    "            bbox = [lons.min(), lons.max(), lats.min(), lats.max()]\n",
    "\n",
    "            # Create 2d array\n",
    "            multiArray = np.reshape(array, (-1, rangeGates))\n",
    "            data = ma.array(multiArray)\n",
    "\n",
    "            # threshVals[0:2] contains halfwords 31,32,33 (min value, increment, num levels)\n",
    "            data = ma.array(threshVals[0]/10. + (multiArray)*threshVals[1]/10.)\n",
    "\n",
    "            if nexrad[code]['unit'] == 'kts':\n",
    "                data[data<-63] = ma.masked\n",
    "                data *= 1.94384 # Convert to knots\n",
    "            else:\n",
    "                data[data<=((threshVals[0]/10.)+threshVals[1]/10.)] = ma.masked\n",
    "\n",
    "            # Save our requested grids so we can render them multiple times\n",
    "            product = {\n",
    "                \"code\": code,\n",
    "                \"bbox\": bbox,\n",
    "                \"lats\": lats,\n",
    "                \"lons\": lons,\n",
    "                \"data\": data\n",
    "            }\n",
    "            grids.append(product)\n",
    "\n",
    "print(\"Processed \"+str(len(grids))+\" grid.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAJ7CAYAAADOeXaMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8U+X3wPHP7aIthdJC2XvvJXuJbBEQBCyyp8pS+QIC\nPwQFUQRxW0FAluy9CsgsQxkioowyyxYKpQW6V+7vjzTpSjrTZvS8X6+8zE3uOD2m7enDuc+jqKqK\nEEIIIYQQInPszB2AEEIIIYQQ1kgKaSGEEEIIIbJACmkhhBBCCCGyQAppIYQQQgghskAKaSGEEEII\nIbJACmkhhBBCCCGywMHcAWSFoigyZ58QQgghhMgVqqoqhl63ykIaQOa/zh4/Pz/atm1r7jAskuTG\nMMmLYZIXw6w1L5cvX2b9+vXMnj07x65hrbnJaZIXwyQv5qcoBmtoQFo7hBBCCL2aNWsmK6KDgoLM\nGI0QwtIp1jiyqyiKao1xCyGEsB5RUVE0btyYEydO4O7ubu5whBBmoiiK0dYOKaSFEEIII2JjY3F0\ndDR3GEIIM0qrkJbWjjzKz8/P3CFYLMmNYZIXwyQvhtlKXpIW0WPHjuX48ePZPqet5MbUJC+GSV4s\nm9XebCiEEELkpvfff5/SpUubOwwhhAWR1g4hhBAik/7++2/u3btHjx49zB2KECKHSWuHEEIIYULx\n8fFoNBpzhyGEMDMppPMo6bkyTnJjmOTFMMmLYbael0aNGtGzZ09Au67BpUuXMnysrecmqyQvhkle\nLJsU0kIIIUQ23L9/n0mTJskItRB5kPRICyGEECYUFxeHg4Pcyy+ErZAeaSGEECIXqKpKmzZtCAgI\nMHcoQohcIIV0HiU9V8ZJbgyTvBgmeTEsr+ZFURS2bdtGxYoVje6TV3OTHsmLYZIXyyaFtBBCCGFC\nxYoV0z9ftmwZ+/fvN2M0QoicJD3SQgghRA75888/cXd3p2rVquYORQiRRWn1SEshLYQQQuSC58+f\n888//9CmTRtzhyKEyAS52VCkIj1XxkluDJO8GCZ5MUzyktrNmzc5ePCg5MYIyYthkhfLJoW0EEII\nkQsaNmzI7Nmz9dvPnj0zYzRCCFOw2daO4OBgVq9ezZo1a8ifPz81atSgevXq1KhRgxo1alCyZEkU\nxeAovRBCCJGjwsLCaNq0KX/++Seurq7mDkcIkYY80yOtqionTpxg8eLF7Nq1i9dee42hQ4cC4O/v\nj7+/P1euXMHf35+IiAiaNGnChx9+SMeOHaWoFkIIkauio6PJly+fucMQQqTD5nukg4KC+Prrr6lZ\nsyZvv/02DRs25MaNG6xZs4aOHTvSsWNH3nvvPRYuXMiRI0d49OgRd+7cYfjw4UyYMIFmzZrh6+uL\nOf+ouHbtGhs3bsTPzw9/f3+Cg4NzNB7puTJOcmOY5MUwyYthkhfjdLlJWkSPHj2aP//800wRWQb5\nzBgmebFsVltIq6qKn58f/fv3p3Llypw/f54lS5Zw+fJlJkyYQJEiRdI83sPDg/79+3PhwgUmT57M\ntGnTaNSoETt27MjVgvrcuXP07duXVq1asX79embOnEnPnj2pWLEibm5uTJ06lRcvXuRaPEIIIXLf\n6NGjqVOnjrnDEEJkktW2dlStWhVHR0fefvttBg4ciKenZ7bOqdFo2LlzJ7Nnz0aj0TBjxgx69eqF\nnZ3p/9ZQVZVjx47x+eefc+nSJSZNmsSoUaPInz9/sv0ePHjARx99xL59+5g1axbDhw/HwcHB5PEI\nIYSwHGfPnuX58+e0b9/e3KEIIbDRHunRo0fj4+Nj8t5mVVXx9fVl9uzZREREMGPGDPr06YO9vb3J\nzv35558TFBTElClTGDhwYLo9cufOnWPChAkEBwfz9ddf07Fjx2zHIoQQwjIdP36ckJAQevToYe5Q\nhBDYaI90uXLlcuQGQUVR6NatG6dPn2bBggV8++231K5dm9WrVxMXF5elc8bFxbFu3Trq1avHjBkz\n+OCDD/D392fEiBEZutGkYcOG+Pn5MXv2bEaPHs1rr72Gv79/lmLRkZ4r4yQ3hkleDJO8GCZ5MS69\n3LRu3VpfRKuqyvXr13MhKvOTz4xhkhfLZrV9Atlt5UiPoih06dKFzp07c+jQIWbNmsWsWbOYPn06\nAwYMwNHRMd1zREdHs3LlSubNm0epUqWYN28eXbp0ydIfAIqi0KtXL7p27YqPjw+tW7emUKFCNG/e\nnGbNmtG8eXPq1q2bbuuHRqMhIiKCkJAQbt26RUREBOHh4YSHhxMREYFGo8HLy4uiRYtStGhR8ufP\nbzDe4OBg/Uwoly9fJiAggDp16tClSxeaNm0qLShCCGECN27c4P3338fX11dmlxLCAllta8fmzZvp\n3bt3rl1TVVWOHj3K7NmzuX37Nu+//z516tShQoUKlClTJlnhGBoays8//8w333xD/fr1mTZtGq1a\ntTJpPBqNBn9/f06ePMmpU6c4efIkd+7coWHDhuTPn19fGCf9b3h4ONHR0bi4uJA/f35cXV2T/VfX\nox0UFMTjx48JDAwE0BfVRYsWJSIigsuXLxMREUGNGjWoWbMmNWrUoEKFCpw7d459+/Zx+/Zt2rdv\nT+fOnencuTNly5Y16dcuhBB5iaqq+iJao9HkyL07QgjjbLJH+siRI7Rt29Ys1z9x4gTLli3j5s2b\n3Lp1i8DAQEqWLEmFChUoUaIE+/fvp3379kydOpX69evnWlzPnj3j7NmzxMbGJiuOkxbMLi4umRrV\nCA8P5/Hjxzx58oTAwECcnZ2pUaMGpUqV0n2wiIqKwsXFRX/Mo0eP2L9/P/v27ePAgQN4eXnpR/fb\ntGmTbF8hhBAZEx8fT4sWLdi2bRslS5Y0dzhC5Bk2WUj/888/1K1b19yhABATE8Pdu3e5desWd+/e\npU2bNlSpUsXcYaXJz88vy3+IqKrKxYsXWbduHevXr+fevXsULlyYSpUqpXpUqFCBv/76i7Vr17J9\n+3Y0Gg0vv/wyAwYMYNCgQab9okwkO7mxZZIXwyQvhklejMtObh48eECpUqVMG5CFkM+MYZIX80ur\nkLbaRtac7pHODCcnJ8qXL4+7uzvlypXj6dOn3Lhxg+DgYIKDgwkJCdE/f/78OQ4ODvrR4pStFRnZ\ndnFxMcs/7d24cYP169ezbt06QkND6devH1u2bKFevXr8999/3Lx5U//YvXs3N2/eJCAgADs7OwoX\nLkz9+vUpXLgwnp6ePHv2LNfjF0IIa5e0iP7ll1+oVq2ayVsHhRAZZ7Uj0uHh4bi6uubqdZ8+fcrW\nrVs5ePAgQUFB+gI5JCSE8PBwChUqhKenJx4eHnh6ehp87u7uTlxcnL5nOWUPc0a2o6KicHZ2znIh\nnt62s7Ozvv3j3r17bNy4kXXr1nH//n369u3LW2+9RbNmzaRPTwghzMjPz48yZcpQqVIlc4cihE2z\nydYOjUaTK3cwv3jxgu3bt7NhwwZOnDhB586d6datGyVLlkxWKBcoUCDXCkuNRkNkZKTRQjsrxXnS\n7ZiYGH1xHRcXR69evejXrx9t27aV2TiEEMICBQcHc/PmTRo3bmzuUISwOTbZ2pGTRXR4eDi+vr6s\nX7+eQ4cO0bZtWwYOHMiGDRtwc3PLsetmlJ2dXbJZNrIirZ6ruLg4faHu6emJk5NTlq9jjaQfzTDJ\ni2GSF8MkL8blRG78/f3Zv3+/VRfS8pkxTPJi2ay2kDa16Oho9u3bx/r169m7dy/NmjXD29ubX375\nBQ8PD3OHl6scHBwoUKAABQoUMHcoQphMz3d/1T/fvsgyb3QVIqtatmxJy5Yt9duhoaHyM1yIXGC1\nrR2mijsiIoJPPvmEpUuXUrduXfr160fv3r3x8vIyyfmFEOaTtHhOSgppYctCQkJo1aoV586dy9Dq\nuUKItNlkj7Qp4j569CgjRoygSZMmzJ8/n9KlS5sgOiGEJZGRaJEXRUZGypz9QphIWoV0npx2ITQ0\nlLFjx9K/f3+++uor1q5
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1084e12d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cartopy.feature as cfeat\n",
    "from cartopy.feature import ShapelyFeature\n",
    "from shapely.geometry import Polygon\n",
    "\n",
    "for rec in grids:\n",
    "    code = rec[\"code\"]\n",
    "    bbox = rec[\"bbox\"]\n",
    "    lats = rec[\"lats\"]\n",
    "    lons = rec[\"lons\"]\n",
    "    data = rec[\"data\"]\n",
    "    # Create figure\n",
    "    %matplotlib inline\n",
    "    fig, ax = make_map(bbox=bbox)\n",
    "    # Colortable filename, beginning value, increment\n",
    "    ctable = nexrad[code]['ctable'][0]\n",
    "    beg = nexrad[code]['ctable'][1]\n",
    "    inc = nexrad[code]['ctable'][2]\n",
    "\n",
    "    norm, cmap = ctables.registry.get_with_steps(ctable, beg, inc)\n",
    "    cs = ax.pcolormesh(lons, lats, data, norm=norm, cmap=cmap)\n",
    "    ax.set_aspect('equal', 'datalim')\n",
    "\n",
    "    cbar = plt.colorbar(cs, extend='both', shrink=0.75, orientation='horizontal')\n",
    "    cbar.set_label(site.upper()+\" \"+ str(nexrad[code]['res']/1000.) +\"km \" \\\n",
    "                   +nexrad[code]['name']+\" (\"+code+\") \" \\\n",
    "                   +nexrad[code]['unit']+\" \" \\\n",
    "                   +str(record.getDataTime()))\n",
    "\n",
    "    # Zoom to within +- some value deg of center\n",
    "    ax.set_xlim(lon-4., lon+4.)\n",
    "    ax.set_ylim(lat-4., lat+4.)\n",
    "    for i, geom in enumerate(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",
    "        # For now do not plot SPS\n",
    "        if (parameters['phensig'][i]):\n",
    "            color = warning_color(parameters['phensig'][i].astype('U13'))\n",
    "            shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(),\n",
    "                            linewidth=2,\n",
    "                            facecolor='none', edgecolor=color)\n",
    "            ax.add_feature(shape_feature)\n",
    "\n",
    "\n",
    "    political_boundaries = NaturalEarthFeature(category='cultural',\n",
    "                                   name='admin_0_boundary_lines_land',\n",
    "                                   scale='50m', facecolor='none')\n",
    "    states = NaturalEarthFeature(category='cultural',\n",
    "                                   name='admin_1_states_provinces_lines',\n",
    "                                   scale='50m', facecolor='none')\n",
    "    ax.add_feature(political_boundaries, linestyle='-')\n",
    "    ax.add_feature(states, linestyle=':')\n",
    "    plt.show()"
   ]
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}