python-awips/examples/notebooks/Map_Resources_and_Topography.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The python-awips package provides access to the entire AWIPS Maps Database for use in Python GIS applications.  Map objects are returned as <a href=\"http://toblerity.org/shapely/manual.html\">Shapely</a> geometries (*Polygon*, *Point*, *MultiLineString*, etc.) and can be easily plotted by Matplotlib, Cartopy, MetPy, and other packages. \n",
    "\n",
    "Each map database table has a geometry field called `the_geom`, which can be used to spatially select map resources for any column of type geometry,\n",
    "\n",
    "## Notes\n",
    "\n",
    "\n",
    "* This notebook requires: **python-awips, numpy, matplotplib, cartopy, shapely**\n",
    "* Use datatype **maps** and **addIdentifier('table', &lt;postgres maps schema&gt;)** to define the map table:\n",
    "        DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "        request = DataAccessLayer.newDataRequest('maps')\n",
    "        request.addIdentifier('table', 'mapdata.county')\n",
    "* Use **request.setLocationNames()** and **request.addIdentifier()** to spatially filter a map resource. In the example below, WFO ID **BOU** (Boulder, Colorado) is used to query counties within the BOU county watch area (CWA) \n",
    "    \n",
    "        request.addIdentifier('geomField', 'the_geom')\n",
    "        request.addIdentifier('inLocation', 'true')\n",
    "        request.addIdentifier('locationField', 'cwa')\n",
    "        request.setLocationNames('BOU')\n",
    "        request.addIdentifier('cwa', 'BOU')\n",
    "        \n",
    "See the <a href=\"http://unidata.github.io/awips2/python/maps-database/#mapdatacwa\">Maps Database Reference Page</a> for available database tables, column names, and types. \n",
    "  \n",
    " > Note the geometry definition of `the_geom` for each data type, which can be **Point**, **MultiPolygon**, or **MultiLineString**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from awips.dataaccess import DataAccessLayer\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import numpy as np\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "from cartopy.feature import ShapelyFeature,NaturalEarthFeature\n",
    "from shapely.geometry import Polygon\n",
    "from shapely.ops import cascaded_union\n",
    "\n",
    "# Standard map plot\n",
    "def make_map(bbox, projection=ccrs.PlateCarree()):\n",
    "    fig, ax = plt.subplots(figsize=(12,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.top_labels = gl.right_labels = False\n",
    "    gl.xformatter = LONGITUDE_FORMATTER\n",
    "    gl.yformatter = LATITUDE_FORMATTER\n",
    "    return fig, ax\n",
    "\n",
    "# Server, Data Request Type, and Database Table\n",
    "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest('maps')\n",
    "request.addIdentifier('table', 'mapdata.county')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Request County Boundaries for a WFO\n",
    "\n",
    "* Use **request.setParameters()** to define fields to be returned by the request."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 23 county MultiPolygons\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAI3CAYAAACCmn2HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAB4JklEQVR4nO39d5xc5X33/7+uM2V7kVYdISSQhEBCNIEQVQIEMr2IYhtjXIITO45zJ05if3+/b+47fyR2Yqc5vuNgYwwYg8GYIjpCsDQBKghJqPfedrW9TTnX94+ZHe1Ku9KudmfOmZn382GsnZkz53zmnCnvueY612WstYiIiIiISHeO1wWIiIiIiPiRgrKIiIiISA8UlEVEREREeqCgLCIiIiLSAwVlEREREZEeBL0uoDfz5s2zNTU1XpeRUyKRCOFw2Osy8p6Og3/oWHjjs88+Ix6Pc+655xIMBnUcfELHwR90HDJv+fLlb1hr5/V0m2+Dck1NDcuWLfO6jJxSXV3N7NmzvS4j7+k4+IeOhTeqqqo4fPgwb731FsOGDdNx8AkdB3/Qccg8Y8yw3m5T1wsREckoYwwAGsdfRPxOQVlERDJKQVlEsoWCsoiIZJSCsohkCwVlERHJKAVlEckWCsoiIpJRCsoiki0UlEVEJKMUlEUkWygoi4hIRikoi0i2UFAWEZGM6gzKIiJ+p6AsIiKeUIuyiPidgrKIiGSUul6ISLZQUBYRkYxSUBaRbKGgLCIiGaWgLCLZQkFZREQySkFZRLKFgrKIiGSUgrKIZAsFZRERySgFZRHJFgrKIiKSUQrKIpItFJRFRCSjFJRFJFsoKIuISEYpKItItlBQFhGRjFJQFpFsoaAsIiIZpaAsItlCQVlERDJKQVlEsoWCsoiIZJSCsohkCwVlERHJKAVlEckWCsoiIpJRnUFZRMTvFJRFRMQTalEWEb9TUBYRkYxS1wsRyRYKyiIiklEKyiKSLRSURUQkoxSURSRbKCiLiEhGKSiLSLZQUBYRkYxSUBaRbKGgLCIiGaWgLCLZQkFZREQySkFZRLKFgrKIiGSUgrKIZAsFZRERySjNzCci2UJBWUREPKEWZRHxOwVlERHJKHW9EJFsoaAsIiIZpaAsItlCQVlERDJKQVlEsoWCsoiIZJSCsohkCwVlERHJKAVlEckWCsoiIpJRCsoiki0UlEVEJKMUlEUkWygoi4hIRmnCERHJFgrKIiLiCbUoi4jf9TkoG2MCxpgVxpiXk5d/YoxZb4xZZYx53hhT2WXZnxhjlhljrkpeHm+MscaY73ZZ5ufGmAcG76GIiEg2UNcLEckW/WlR/h6wrsvlhcA0a+10YCPwQwBjzJTk7VcC3+my/EHge8aY8MmXKyIi2U5BWUSyRZ+CsjFmLHAj8HDnddbaN621seTFj4Gxyb8DgAtYoGtHtEPAIuCrA6xZRESymIKyiGSLvrYo/wfwtyQCcE++DrwGYK1dAxQDHwC/OGq5HwN/bYwJ9LtSERHJCQrKIpItzIneqIwxNwE3WGu/bYyZDXzfWntTl9v/f8AM4A7by8qMMeOBl62104wxj5PotjETWGatfbSn+0yfPt3+7Gc/6/cDkt41NzdTWlrqdRl5T8fBP3QsvLFx40aampqYPHkyZWVlOg4+oePgDzoOmTdnzpzl1toZPd0W7MP9LwNuMcbcABQC5caYJ6y19xljvgrcBFzTW0juwT8BzwLvHW+hcDjM7Nmz+7hK6Yvq6mrtUx/QcfAPHQtv/MM//APV1dUsWrSI2bNn6zj4hI6DP+g4+MsJu15Ya39orR1rrR0P3Au8nQzJ84C/A26x1rb2dYPW2vXAWhIBW0RE8oy6XohIthjIOMo/B8qAhcaYz4wx/9OP+/4jR07+ExGRPKKgLCLZoi9dL1KstdVAdfLvif2433ZgWpfLK9FkJyIieUlBWUSyhcKqiIhklIKyiGQLBWUREckoBWURyRYKyiIiklGdQVlExO8UlEVExBNqURYRv1NQFhGRjFLXCxHJFgrKIiKSUQrKIpItFJRFRCSjFJRFJFsoKIuISEYpKItItlBQFhGRjFJQFpFsoaAsIiIZpaAsItlCQVlERDJKQVlEsoWCsoiIZJSCsohkCwVlERHJKAVlEckWCsoiIpJRCsoiki0UlEVEJKMUlEUkWygoi4hIRikoi0i2UFAWEZGMUlAWkWyhoCwiIhmloCwi2UJBWUREMqozKIuI+J2CsoiIeEItyiLidwrKIiKSUep6ISLZQkFZREQySkFZRLKFgrKIiGSUgrKIZAsFZRERySgFZRHJFgrKIiKSUQrKIpItFJRFRCSjFJRFJFsoKIuISEYpKItItlBQFhGRjFJQFpFsoaAsIiIZpaAsItlCQVlERDJKQVlEsoWCsoiIZJSCsohkCwVlERHJKAVlEckWCsoiIpJRCsoiki0UlEVEJKMUlEUkWygoi4hIRnUGZRERv1NQFhERT6hFWUT8TkFZREQySl0vRCRbKCiLiEhGKSiLSLZQUBYRkYxSUBaRbKGgLCIiGaWgLCLZQkFZREQySkFZRLKFgrKIiGSUgrKIZAsFZRERySgFZRHJFgrKIiKSUQrKIpItFJRFRCSjFJRFJFsoKIuISEYpKItItlBQFhGRjFJQFpFsoaAsIiIZ1RmURUT8TkFZREQ8oRZlEfE7BWUREckodb0QkWyhoCwiIhmloCwi2UJBWUREMkpBWUSyhYKyiIhklIKyiGQLBWUREckoBWURyRYKyiIiklEKyiKSLRSURUQkoxSURSRbKCiLiEhGKSiLSLZQUBYRkYzSzHwiki0UlEVExBNqURYRv1NQFhGRjFLXCxHJFgrKIiKSUQrKIpItFJRFRCSjFJRFJFsoKIuISEYpKItItlBQFhGRjFJQFpFsoaAsIiIZpaAsItlCQVlERDJKQVlEsoWCsoiIZJQmHBGRbKGgLCIinlCLsoj4nYKyiIhklLpeiEi2UFAWEZGMUlAWkWyhoCwiIhmloCwi2UJBWUREMkpBWUSyhYKyiIhklIKyiGQLBWUREckoBWURyRYKyiIiklEKyiKSLRSURUQkoxSURSRbKCiLiEhGKSiLSLZQUBYRkYxSUBaRbKGgLCIiGaWgLCLZQkFZREQySkFZRLKFgrKIiGSUgrKIZAsFZRERyajOoCwi4ncKyiIi4gm1KIuI3ykoi4hIRqnrhYhkCwVlERHJKAVlEckWCsoiIpJRCsoiki0UlEVEJKMUlEUkWygoi4hIRikoi0i2UFAWEZGMUlAWkWzR56BsjAkYY1YYY15OXr7LGLPGGOMaY2YctexPjDHLjDFXJS+PN8ZYY8x3uyzzc2PMA4P0OEREJEsoKItItuhPi/L3gHVdLn8O3AG813UhY8yU5J9XAt/pctNB4HvGmPBJ1CkiIjlCQVlEskWwLwsZY8YCNwL/CPwVgLV2XfK2oxcPAC5gga43HgI+BL4K/GogRYuISPbq/NzYunUrzz//PMYYnn/++bRuK11yYZbBrl9cFixYkNZt9EfnF6l07OOTWedg12GtTa2zsLCQK664gnBYbYl+06egDPwH8LdA2YkWtNauMcYUAx8Af3PUzT8GXjPGPHKi9Sxfvjwn3oD85Kc//Slz5szp131mXnkVt37pvj4v77ouANa6WAsk3+hs8jps4s0h8QZoOW6DkrXYLn8DYMwxzwvrurjJ/+LxOG48juu6WOsSjUaJJf+LRCK48TjWWqLR6JHlXDdVU2odyetdN7H8uNMnsmv7VuKxWOIxJGtPPE5L1wdiAeva5F8ABuOY1LfGL981n7//px8d+3CT/zrGSa03sVrbp5a3RB0932YcQyAQwHV7XqBrfRiDwYABYxyM4cjlxJXH3t8kH2fyPt1vTHzAHH3fI/c59pgaY3CcY3/w6nzedHKT+8U5yfeKW6+/jp/89y9SqywuLaGouITagwePWbazxng8TigcoqCgkFCXDzXHcXACARzHSew3x+AkH4dxnMTtToBwQTj13HGcQOLxO537wUntl56Y1DE5cjxs8vmReB5321tYS6/rMwaM46TWYXrZpsV2fSofWc4
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define a WFO ID for location\n",
    "# tie this ID to the mapdata.county column \"cwa\" for filtering\n",
    "request.setLocationNames('BOU')\n",
    "request.addIdentifier('cwa', 'BOU')\n",
    "\n",
    "# enable location filtering (inLocation)\n",
    "# locationField is tied to the above cwa definition (BOU)\n",
    "request.addIdentifier('geomField', 'the_geom')\n",
    "request.addIdentifier('inLocation', 'true')\n",
    "request.addIdentifier('locationField', 'cwa')\n",
    "\n",
    "# This is essentially the same as \"'\"select count(*) from mapdata.cwa where cwa='BOU';\" (=1)\n",
    "\n",
    "# Get response and create dict of county geometries\n",
    "response = DataAccessLayer.getGeometryData(request, [])\n",
    "counties = np.array([])\n",
    "for ob in response:\n",
    "    counties = np.append(counties,ob.getGeometry())\n",
    "print(\"Using \" + str(len(counties)) + \" county MultiPolygons\")\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "# All WFO counties merged to a single Polygon\n",
    "merged_counties = cascaded_union(counties)\n",
    "envelope = merged_counties.buffer(2)\n",
    "boundaries=[merged_counties]\n",
    "\n",
    "# Get bounds of this merged Polygon to use as buffered map extent\n",
    "bounds = merged_counties.bounds\n",
    "bbox=[bounds[0]-1,bounds[2]+1,bounds[1]-1.5,bounds[3]+1.5]\n",
    "\n",
    "\n",
    "fig, ax = make_map(bbox=bbox)\n",
    "# Plot political/state boundaries handled by Cartopy\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='-', edgecolor='black')\n",
    "ax.add_feature(states, linestyle='-', edgecolor='black',linewidth=2)\n",
    "\n",
    "# Plot CWA counties\n",
    "for i, geom in enumerate(counties):\n",
    "    cbounds = Polygon(geom)\n",
    "    intersection = cbounds.intersection\n",
    "    geoms = (intersection(geom)\n",
    "         for geom in counties\n",
    "         if cbounds.intersects(geom))\n",
    "    shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), \n",
    "                        facecolor='none', linestyle=\"-\",edgecolor='#86989B')\n",
    "    ax.add_feature(shape_feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a merged CWA with cascaded_union"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAI3CAYAAACCmn2HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAB/h0lEQVR4nO3dd5wb1b028OfMjMr2Xe+u1zauYBuDjTGmGjDYBIPpvdcEAgkkITch7b65uTe5N4FcSCMkuRDTu0PHVGNYmsG94d572d5VRnPeP0arlbZ4d72rOSPp+X4+YM1oNPpJs5IeHZ05R0gpQUREREREiTTVBRARERERuRGDMhERERFRJxiUiYiIiIg6waBMRERERNQJBmUiIiIiok4YqgvoysyZM2VlZaXqMtJKKBSC1+tVXUbG43FwDx4LNZYvX45IJIJjjz0WhmHwOLgEj4M78Dg4b8mSJe9LKWd2dp1rg3JlZSUWL16suoy0Ul5ejmnTpqkuI+PxOLgHj4UaxcXFqK6uxocffoiSkhIeB5fgcXAHHgfnCSFKurqOXS+IiMhRQggAAMfxJyK3Y1AmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiRzEoE1GqYFAmIiJHtQZlIiK3Y1AmIiIl2KJMRG7HoExERI5i1wsiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkqNagTETkdgzKRESkBFuUicjtGJSJiMhR7HpBRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWZ+YgoVTAoExGREmxRJiK3Y1AmIiJHsesFEaUKBmUiInIUgzIRpQoGZSIichSDMhGlCgZlIiJyFIMyEaUKBmUiInIUgzIRpQoGZSIichSDMhGlCgZlIiJyFIMyEaUKBmUiInIUJxwholTBoExEREqwRZmI3K7HQVkIoQshlgkh5kSXHxBCrBNCrBRCvCaEKIzb9gEhxGIhxJnR5ZFCCCmE+H7cNg8LIW7tv4dCRESpgF0viChV9KZF+R4Aa+OW5wKYIKWcCGADgF8AgBBiXPT6MwDcHbf9AQD3CCG8h14uERGlOgZlIkoVPQrKQoihAC4AMKt1nZTyAymlGV38CsDQ6GUdgAVAAojviFYBYB6AW/pYMxERpTAGZSJKFT1tUf4zgJ/CDsCd+RaAdwFASrkaQDaAzwH8o9129wP4sRBC73WlRESUFhiUiShViO7eqIQQFwI4X0p5lxBiGoB7pZQXxl3//wCcAOBy2cXOhBAjAcyRUk4QQjwNu9vGyQAWSymf7Ow2EydOlA899FCvHxB1rbGxEbm5uarLyHg8Du7BY6HGhg0b0NDQgLFjxyIvL4/HwSV4HNyBx8F506dPXyKlPKGz64we3P40ABcLIc4H4AeQL4R4Vkp5oxDiFgAXAvhGVyG5E78D8DKATw+2kdfrxbRp03q4S+qJ8vJyPqcuwOPgHjwWavz6179GeXk55s2bh2nTpvE4uASPgzvwOLhLt10vpJS/kFIOlVKOBHAtgI+iIXkmgJ8BuFhK2dzTO5RSrgOwBnbAJiKiDMOuF0SUKvoyjvLDAPIAzBVCLBdC/F8vbvtbtJ38R0REGYRBmYhSRU+6XsRIKcsBlEcvj+7F7bYBmBC3vAKc7ISIKCMxKBNRqmBYJSIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIke1BmUiIrdjUCYiIiXYokxEbsegTEREjmLXCyJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROSo1qBMROR2DMpERKQEW5SJyO0YlImIyFHsekFEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgcxaBMRKmCQZmIiBzFoExEqYJBmYiIHMWgTESpgkGZiIgc1RqUiYjcjkGZiIiUYIsyEbkdgzIRETmKXS+IKFUwKBMRkaMYlIkoVTAoExGRoxiUiShVMCgTEZGjGJSJKFUwKBMRkaMYlIkoVTAoExGRoxiUiShVMCgTEZGjGJSJKFUwKBMRkaMYlIkoVTAoExGRoxiUiShVMCgTEZGjGJSJKFUwKBMRkaMYlIkoVTAoExGRo1qDMhGR2zEoExGREmxRJiK3Y1AmIiJHsesFEaUKBmUiInIUgzIRpQoGZSIichSDMhGlCgZlIiJyFIMyEaUKBmUiInIUgzIRpQoGZSIichSDMhGlCgZlIiJyFIMyEaUKBmUiInIUgzIRpQoGZSIichRn5iOiVMGgTERESrBFmYjcjkGZiIgcxa4XRJQqGJSJiMhRDMpElCoYlImIyFEMykSUKhiUiYjIUQzKRJQqGJSJiMhRDMpElCoYlImIyFEMykSUKhiUiYjIUQzKRJQqGJSJiMhRnHCEiFIFgzIRESnBFmUicjsGZSIichS7XhBRqmBQJiIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiRzEoE1GqYFAmIiJHMSgTUapgUCYiIkcxKBNRqmBQJiIiR7UGZSIit2NQJiIiJdiiTERux6BMRESOYtcLIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFQzKRETkKAZlIkoVDMpEROQoBmUiShUMykRE5CgGZSJKFT0OykIIXQixTAgxJ7p8lRBitRDCEkKc0G7bB4QQi4UQZ0aXRwohpBDi+3HbPCyEuLWfHgcREaUIBmUiShW9aVG+B8DauOWvAVwO4NP4jYQQ46IXzwBwd9xVBwDcI4TwHkKdRESUJhiUqSeklJCWpbqMpOHff2owerKREGIogAsA/BbAjwBASrk2el37zXUAFgAJIP7KCgBfALgFwD/7UjQ5I1x7AJvvvxRNGxcAPX2zEiL6nxb92xBx66PropchBAQOMkNX67b2gn1ZWh1r0XQITYfQDQjDA6F77HVCg/D6oXmz7P982dHrNGgeP4ThsWvSjWhtGqDrELrH3pemQ2gGoGkQmg4pZVv90XqE0ABNQ8Kfeuv61tqlhJQWEH1TNMVw7Hvtgc4fLwBYEbue1v0mPA8He7q0LreTVgTSDENoeufXx9UHSPtyQt2ty9HLHXeQeJv49V3dtsN9JBQMaUU6e5CJfzNa9Lv+IX6YhgtPwM7H/w042N9hXE0AIHQPrFALIi0NkKEW+yopASsCGQnbdVsR+3FZFmTEbFtnmbACTW2vgYgZffyWHQjin+/OSog9z1bb8xb7O9TbvR+3vmY6319rza376HAMYrsRCX/LrdtJSMCy7L8prX2bS8f7FrondryPbajH+zMBz74nsfC62Wg+5z+x4NHLo7frrAgNQCd/Jx3uVcY91Oh7hWitR098fXS2LyEghf1YhEz8m5Ka/T4By2x7D4j+bQsZf+zQ7n6ifzdSQkRC9u0BSMMPqXuQ+IBl4nMWPdaJoh+treuFhoTjKy0kvF8m3FTa+wQgRfSVFHd/zRf
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot CWA envelope\n",
    "for i, geom in enumerate(boundaries):\n",
    "    gbounds = Polygon(geom)\n",
    "    intersection = gbounds.intersection\n",
    "    geoms = (intersection(geom)\n",
    "         for geom in boundaries\n",
    "         if gbounds.intersects(geom))\n",
    "    shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), \n",
    "                        facecolor='none', linestyle=\"-\",linewidth=3.,edgecolor='#cc5000')\n",
    "    ax.add_feature(shape_feature)\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WFO boundary spatial filter for interstates\n",
    "\n",
    "Using the previously-defined **envelope=merged_counties.buffer(2)** in **newDataRequest()** to request geometries which fall inside the buffered boundary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 225 interstate MultiLineStrings\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAI3CAYAAACCmn2HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAADAbUlEQVR4nOzdd3xb5dXA8d+92t4zzt5775AFYQRCCHuFTYFCW9rSUkrH2z1poW9pS9uXPUrZe4RAEjCQhCRkh+y9471tzXvfPx55JU5sx7KvZJ3v56PEkq6ujvbRo/OcRzNNEyGEEEIIIURjutUBCCGEEEIIEY0kURZCCCGEEKIJkigLIYQQQgjRBEmUhRBCCCGEaIIkykIIIYQQQjTBbnUAJzNnzhyzsLDQ6jA6Fb/fj9PptDqMuCePg1VMKFkPaaNBswHyWFhl/fr1hEIhxowZg91ul8chSsjjEB3kceh4a9as+dA0zTlNnRe1iXJhYSGrV6+2OoxOJTc3l1mzZlkdRtyTx8Eixetg2Xy4eF3dSfJYWCMzM5Pi4mIWL15MVlaWPA5RQh6H6CCPQ8fTNC3rZOdJ6YUQIj6kDIZAOZRttTqSuKdpGgDSx18IEe0kURZCxAd7IqSNgqr9VkcS9yRRFkLECkmUhRDxofoIlKyDpH5WRxL3JFEWQsQKSZSFEPGhcDlkToGUIVZHEvckURZCxApJlIUQ8cGRAsEqq6MQSKIshIgdkigLIeLD0Y8ge7rVUQgkURZCxA5JlIUQ8eHAK9BnvtVRCCRRFkLEDkmUhRCdn79EHVKHWx2JoD5RFkKIaCeJshCi8yv9ClJHgCZvedFERpSFENFOPjWEEJ1f2RaVKIuoIKUXQohYIYmyEKLzqz4Eib2tjkKESaIshIgVkigLITo/TQdJyqKGJMpCiFghibIQovOzucDwWR2FCJNEWQgRKyRRFkJ0froTDL/VUYgwSZSFELFCEmUhRBzQpPQiikiiLISIFZIoCyE6P8MPNqfVUYgwSZSFELFCEmUhROcXqgGbx+ooRJgkykKIWCGJshCi8zOCoNmtjkKESaIshIgVkigLITo/MwC6w+ooRJgkykKIWCGJshCi86vYCQm9rI5ChEmiLISIFZIoCyE6v9RRcGyx1VGIMEmUhRCxQhJlIUTnN/R7cOgtqNpvdSQCSZSFELFDEmUhROfnyoSBX4fNf7Q6EoEkykKI2CGJshAiPgy9D/a/DBV7rY4k7tUmykIIEe0kURZCxAd3FmRNh6VXWx2JCJMRZSFEtJNEWQgRP0b8CErWQ+EqqyOJa1J6IYSIFZIoCyHiR5eZMOJnsOPvVkcS1yRRFkLECkmUhRDxZdj34ehCqNxjdSRxSxJlIUSskERZCBFfnKkw4E7Y/ojVkcQtSZSFELFCEmUhRPzpeg6UrLU6irglibIQIlZIoiyEiD9JA6Byt9VRxC1JlIUQsUISZSFE/EnoBd58wLA6krgkibIQIlZIoiyEiD+6HTzdwQhYHUlckkRZCBErJFEWQsQf0wR/CWg2qyOJS7IynxAiVkiiLISIP2YQDL9KmIVlZERZCBHtJFEWQsQf3QHDfgjVB62OJC5J6YUQIlZIoiyEiE/DfgiBMghWWR1J3JFEWQgRKyRRFkLEJ0cS2NxQ9KXVkcQdSZSFELFCEmUhRPxyZcH6H0Gg3OpI4ookykKIWCGJshAifrmyIWMCLDlHkuUOJImyECJWSKIshIhvE/8Jnh6w70WrI4kbkigLIWKFJMpCiPimaZAyWE3sEx1CEmUhRKyQRFkIITS76q0sOoQsOCKEiBWSKAshhKaDaVgdRdyREWUhRLRrcaKsaZpN07R1mqa9Fz7+oKZp2zRN26hp2puapqU12PZBTdNWa5p2Vvh4X03TTE3TvtNgm0c0Tbs1cjdFCCFOk2YDM2R1FHFDSi+EELGiNSPK9wBbGxxfBIw0TXM0sAP4CYCmaUPD558J3N1g+3zgHk3TnKcfrhBCtANXNlQdsDqKuCGJshAiVrQoUdY0rSdwEfBE7WmmaX5kmnVFfSuAnuG/bYABmEDDQrQCYAlwSxtjFkKIyOp5KRx8XVbp6yCSKAshYkVLR5QfBu5HJcBNuQ34AMA0zc1AArAU+Pdx2z0A/EDTNFurIxVCiPaS2BsSekHFbqsjiQuSKAshYoXW3BuVpmnzgLmmaX5L07RZwH2mac5rcP7/ABOBK8yT7EzTtL7Ae6ZpjtQ07TlU2cYUYLVpms80dZnRo0ebf//731t9g8TJVVZWkpSUZHUYcU8eh+jR6LGo2AHuruBIsTaoOLBjxw4qKioYPHgwycnJ8pqIEvI4RAd5HDre2WefvcY0zYlNnWdvweWnA5domjYXcAMpmqY9b5rmjZqm3QLMA849WZLchD8ArwGfnWojp9PJrFmzWrhL0RK5ublyn0YBeRyiR6PHYvMKqNkEE/9haUzx4Ne//jW5ubksWbKEWbNmyWsiSsjjEB3kcYguzZZemKb5E9M0e5qm2ReYD3wcTpLnAD8CLjFNs7qlV2ia5jZgCyrBFkKI6JA+Bsq3Wx1FXJDSCyFErGhLH+VHgGRgkaZp6zVN+79WXPb31E/+E0II6yX0gmrpfNERJFEWQsSKlpRe1DFNMxfIDf89sBWX2weMbHB8A7LYiRAimiT1Vy3iQl6wua2OplOTRFkIESskWRVCCAB7AmROgcPvWR1JpyeJshAiVkiiLIQQtfrfCrseA0ng2pUkykKIWCGJshBC1Op9DQTKYNXXwQg2v704LbWJshBCRDtJlIUQopbdA+csgar9kHshFK6yOqJOTUaUhRDRThJlIYRoyJEEsxZA93mw7Br4+AKoPmx1VJ2KlF4IIWKFJMpCCHE83QFD74GLd0KXmfDhZChYbnVUnYYkykKIWCGJshBCnIzugJE/g8mPwWeXgK/I6og6BUmUhRCxQhJlIYRoTo+LoMfFsOnX0hEjAiRRFkLECkmUhRCiJcb+CQpXwOdXgL/E6mhimiTKQohYIYmyEEK0hLsLzP4cEnrDBxOkI0YbSKIshIgVkigLIURL2Vww8W8w/iH4dB5s+5uUYpwGSZSFELFCEmUhhGitXlfABStg3/Pw2WUQKLc6opgiibIQIlZIoiyEEKcjqT/MXgbuHFhyLviKrY4oZkiiLISIFZIoCyHE6bI5YfKj0OVM+PhcCFZbHVFMkERZCBErJFGOJNOEij2QvxS8BW3fX8UuqDnW9v0IIdqPpsG4h8DTHfa/bHU0MUESZSFErLBbHUBM8xWpdlGFX6j/y7dBzWFAA0cqpAyG7nPVIX086DZ1uZAXfIVQfQT2vwRV+2DwtyFlCOhOKNsCm38HpV+pbR0pkDYK3NngzABXpvq5t9dV4Ey18h4QQoBKlrOmQuUuqyOJCZIoCyFihSTKLWWEoOyrcGK8XP1fcxQyJ6sPyKH3QtJASOgBjkQIeqHoCziyAFbcClX7VfIcKAfDD64slex2ORMSesKGn0D1QQj51PHBd8NZ76uVwSp2qiTcVwj+YvV/8Rq1+EHWNHAkqZZV3nzAhJRhkDYCUkeq5FoI0b5KN8O+/8LIX1odSUyQRFkIESskUT4ZXzEUrWyQGK8ETzeVFGdNg2H3Qcrw+lHi49ndkHO2Oox7UCXI/jI1AmxPViNQLZUyWB2OV7wGyndAsAIq90HyQNB0ldAfeFmNSNuckBpOmn1nQsFydVxGooVom2A1bP6j+kWpZB2MfQD6XGt1VDFBEmUhRKyQRNkIgr9UjeYWrQwnx1+osojMiZB5Bgz+LkybCu6s078eR4o6RFLGBHU4GdOEmiNQtlklzXmVsOZ7UL5FlXBkTAgn8+dC6vDWJe8doXYUv3gtHHwdkgbA8B9BQnerIxMC9jwD+Z/AiJ9B5iRVEiVaRIu29xohhDiJzp0o1xwLj/ZsAG9efemCv1iNGPtLIFStSiI83VXimHUGDPmeGoE92WhxrNA0VQqS0AO6nQ95uTBrFZgGVO5VXwryPoZtf1X3Q0IfsCeAzVP/v80NuhvsHvW/ZoNQjdo+VKNqqA0/JA+B7nMgcwroLXhamaZK4Cv3qjZbKUN
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n",
    "request.addIdentifier('table', 'mapdata.interstate')\n",
    "request.addIdentifier('geomField', 'the_geom')\n",
    "request.setParameters('name')\n",
    "interstates = DataAccessLayer.getGeometryData(request, [])\n",
    "print(\"Using \" + str(len(interstates)) + \" interstate MultiLineStrings\")\n",
    "\n",
    "# Plot interstates\n",
    "for ob in interstates:\n",
    "    shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), \n",
    "                        facecolor='none', linestyle=\"-\",edgecolor='orange')\n",
    "    ax.add_feature(shape_feature)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nearby cities\n",
    "\n",
    "Request the city table and filter by population and progressive disclosure level:\n",
    "\n",
    "**Warning**: the `prog_disc` field is not entirely understood and values appear to change significantly depending on WFO site.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Queried 1203 total cities\n",
      "Plotting 57 cities\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAKCCAYAAAA9RiPYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd1xV9f/A8de597IREHFrigNNQVFx5QA1R2qaKzUtzZzlKMvmr7K+tm1PG6aVqzQtR2maoDhy4saNEwcqKOMC997z++PAlctQQeAC9/18PK5yzj3jfQeX+z6fz+f9UVRVRQghhBBCCCHE3dHZOwAhhBBCCCGEKAskuRJCCCGEEEKIQiDJlRBCCCGEEEIUAkmuhBBCCCGEEKIQSHIlhBBCCCGEEIVAkishhBBCCCGEKASSXAkhhBBCCCFEIZDkSgghhBBCCCEKgSRXQgghhBBCCFEIJLkSQgghhBBCiEIgyZUQQgghhBBCFAJJroQQQgghhBCiEEhyJYQQQgghhBCFQJIrIYQQQgghhCgEklwJIYQQQgghRCGQ5EoIIYQQQgghCoEkV0IIIYQQQghRCCS5EkIIIYQQQohCIMmVEEIIIYQQQhQCSa6EEEIIIYQQohBIciWEEEIIIYQQhUCSKyGEEEIIIYQoBJJcCSGEEEIIIUQhkORKCCGEEEIIIQqBJFdCCCGEEEIIUQgkuRJCCCGEEEKIQiDJlRBCCCGEEEIUAkmuhBBCCCGEEKIQSHIlhBBCCCGEEIVAkishhBBCCCGEKASSXAkhhBBCCCFEIZDkSgghhBBCCCEKgSRXQgghhBBCCFEIJLkSQgghhBBCiEIgyZUQQgghhBBCFAJJroQQQgghhBCiEEhyJYQQQgghhBCFQJIrIYQQQgghhCgEklwJIYQQQgghRCGQ5EoIIYQQQgghCoEkV0IIIYQQQghRCCS5EkIIIYQQQohCIMmVEEIIIYQQQhQCSa6EEEIIIYQQohBIciWEEEIIIYQQhUCSKyGEEEIIIYQoBJJcCSGEEEIIIUQhkORKCCGEEEIIIQqBJFdCCCGEEEIIUQgkuRJCCCGEEEKIQiDJlRBCCCGEEEIUAkmuhBBCCCGEEKIQSHIlhBBCCCGEEIVAkishhBBCCCGEKASSXAkhhBBCCCFEIZDkSgghhBBCCCEKgSRXQgghhBBCCFEIJLkSQgghhBBCiEIgyZUQQgghhBBCFAJJroQQQgghhBCiEEhyJYQQQgghhBCFQJIrIYQQQgghhCgEklwJIYQQQgghRCGQ5EoIIYQQQgghCoHB3gGUdT169FDj4uLsHUaZk5aWhrOzs73DcHjyOthJ2jVIjYNy9bVFeR3swmg0cuDAAVxcXAgMDJTXoQSR16JkkNeh+O3cuXO1qqo97B2HI5PkqojFxcWxY8cOe4dR5oSHhxMWFmbvMByevA52cuxbOLcCQv8E5HWwl+joaO69915q1arFjh075HUoQeS1KBnkdSh+iqL42TsGRyfdAoUQorSpPQzOrwKLyd6RODRFUQBQVdXOkQghhCgpJLkSQohSRwc6JzCn2DsQhybJlRBCiOwkuRJCiNLm8MdQrTc4lbN3JA5NkishhBDZSXIlhBCljSkZfALtHYXDk+RKCCFEdpJcCSFEaZN2FQwe9o7C4UlyJYQQIjtJroQQojQxp8GZJVC9j70jcXiSXAkhhMhOkishhChNEg6Aix94Bdg7EocnyZUQQojsJLkSQojSJH4f+DSxdxSCm8mVEEIIkUmSKyGEKE0SDoC3FLMoSaTlSgghRCZJroQQojS5fgi877V3FALpFiiEECInSa6EEKI0MSWCk5e9oxBIciWEECInSa6EEKI0UQxgMdk7CoEkV0IIIXKS5EoIIUoTnQtYUu0dhUCSKyGEEDlJciWEEKWJzklarkoISa6EEEJkJ8mVEEKUKhZQ5KO7JJDkSgghRHbyF1oIIUoTUwro3ewdhUCSKyGEEDlJciWEEKWJJQ10zvaOQiDJlRBCiJwkuRJCiNJENYOit3cUgpvJlRBCCJFJkishhChN0hPA4GHvKEQW0nIlhBAikyRXQghRmiSfBRdfe0chkG6BQgghcpLkSgghShP/EXDoI3tHIZDkSgghRE6SXAkhRGnS6Hk4NR9SYu0dicOT5EoIIUR2klwJIURp4lYF/B+Dg+/bOxKHJ8mVEEKI7CS5EkKI0ubeaXB+NcQfsHckDk2SKyGEENlJciWEEKWNe3VtvqutI+0diUOT5EoIIUR2klwJIURp1G4epMbBjeP2jsRhSXIlhBAiO0muhBCiNPJrDbWHwcF37R2Jw5LkSgghRHaSXAkhRGnV4Gk4swQs6faOxCFJclW2LV26FEVRiI6OBiAmJobAwEAAwsPD6d27NwB//vkn774rFzmEEBpJroQQorRy9YMa/SD9mr0jcUiSXJVtCxYsoH379ixcuPCW2/Xp04cXX3yxmKISQpR0klwJIURp5hMI5lR7R+GQMpMrUfYkJiayadMmfvjhh9smV3PmzGHixIkALF++nNatW9OsWTPuv/9+Ll68CMD06dMZNWoUYWFh1KlTh88++8y6/08//USTJk1o2rQpjz76KACXL19mwIABtGzZkpYtW7Jp06YieqRCiMImyZUQQpRmnnXBIslVYbpy5QrBwcEEBwdTpUoVqlevbl1OS0vLsf2tWq5ee+011q5dW6A44uPj+eqrrwq0r7g7y5Yto0ePHgQEBODr68uuXbvuaL/27duzdetWdu/ezZAhQ3j//Zvz0UVHR7N69Wq2bdvGG2+8QXp6Op9++ikjRoxg1qxZ7Nmzh08//RSAKVOm8Mwzz7B9+3aWLFnC6NGj+eSTT0hOTrYez2AwEB8ff0dxTZ8+nZkzZwIwcuRIFi9efIfPxE15vZezdpEUQoDB3gEIIYS4C+Xqgfm0vaMoUypUqEBUVBSgfSn19PTkueeey7Fd9m6BZrPZ5n6z2cybb75Z4Dgyk6snn3yywMcQBbNgwQKefvppAIYMGcKCBQt46qmnbrvf2bNnGTx4MLGxsaSlpeHv72+9r1evXri4uODi4kKlSpW4ePEic+bMoUaNGqxevZq2bdvi6+sLwNq1azl48KB13+vXr/Pxxx8zfPhw3N3dAXB1dcXHx6fwHvRt3M17WQhHIi1XQghRmjl5gWqxdxRl3rp162jWrBlBQUGMGjWK1NRUa3KVnp5O+/btiYiIoHbt2rz55pu0b9+e3377zaaVoHbt2rz++us0b96coKAga6GEy5cv07VrV5o3b864ceOoVasWcXFxvPjiixw/fpzg4GCmTZuGqqpMmzaNwMBAgoKCWLRoEaC1HISFhTFw4EAaNmzIsGHDZBzYXbhy5Qr//vsvo0ePpnbt2nzwwQcsWrTojp7TSZMmMXHiRPbt28esWbMwGo3W+1xcXKw/6/V6EhISOHnyJL1797Z2Pcx8LRMSEkhJSaFx48bs3r2bF154gdjYWDp16kSnTp0ASE5OJi4uDoCPPvqIwMBAatWqRa1atazdEj///HOaNGnC119/zYIFCwA4cuQIzz77LK1bt6ZSpUrUq1eP4OBgxo0bh9lsxmw2M3LkSOv77OOPPwZsW7z+/vtvGjZsSPv27fn999+tjyspKYlRo0bRsmVLmjVrRmRk5N28FEKUSpJcCSFEaZZ0GnRO9o6iTDMajYwcOZJFixaxb98+TCYTX3/9tU3LVWRkJJ07dwa0FoXIyEiGDBmS41h+fn7s2rWLCRMmWLtpvfHGG3Tu3Jldu3bRr18/Tp/WWiLfffdd6tatS1RUFB988AG///47UVFR7Nmzh7Vr1zJt2jRiY2MB2L17N5988gkHDx7kxIkTMkbnLixevJjHHnuMU6dOERMTw5kzZ/D39+fs2bO33TchIYHq1asDMHfu3Ftuu2bNGjp37sy6devw8vJi165dXL9+nd27d/PAAw8watQo62vZsWNHqlWrxvr161m/fr3NcXbu3MmPP/7If//9R2RkJB4eHsyePZsOHTrw+uuv8++//zJhwgT69u3LF198wZkzZ2jRogWzZ8+mUaNGqKpKVFQUer2eefP
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n",
    "request.addIdentifier('table', 'mapdata.city')\n",
    "request.addIdentifier('geomField', 'the_geom')\n",
    "request.setParameters('name','population','prog_disc')\n",
    "cities = DataAccessLayer.getGeometryData(request, [])\n",
    "print(\"Queried \" + str(len(cities)) + \" total cities\")\n",
    "\n",
    "citylist = []\n",
    "cityname = []\n",
    "# For BOU, progressive disclosure values above 50 and pop above 5000 looks good\n",
    "for ob in cities:\n",
    "    if ob.getString(\"population\"):\n",
    "        if ob.getNumber(\"prog_disc\") > 50:\n",
    "            if int(ob.getString(\"population\")) > 5000:\n",
    "                citylist.append(ob.getGeometry())\n",
    "                cityname.append(ob.getString(\"name\"))\n",
    "print(\"Plotting \" + str(len(cityname)) + \" cities\")\n",
    "\n",
    "# Plot city markers\n",
    "ax.scatter([point.x for point in citylist],\n",
    "       [point.y for point in citylist],\n",
    "       transform=ccrs.PlateCarree(),marker=\"+\",facecolor='black')\n",
    "# Plot city names\n",
    "for i, txt in enumerate(cityname):\n",
    "    ax.annotate(txt, (citylist[i].x,citylist[i].y),\n",
    "                xytext=(3,3), textcoords=\"offset points\")\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lakes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 208 lake MultiPolygons\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAKCCAYAAAA9RiPYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzdd1xV9f/A8de5g8veiOBA3AMUwT3Rcptm6tfKSjOzMtP2+ParrG87sz1samVl2dQsTRO3KSrurThRQWWPy733/P44cAVBBb1wGe/n43H1nv0+d3He57MUVVURQgghhBBCCHFtdM4OQAghhBBCCCFqAkmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHECSKyGEEEIIIYRwAEmuhBBCCCGEEMIBJLkSQgghhBBCCAeQ5EoIIYQQQgghHMDg7ABquoEDB6opKSnODqNaMtts7M1Is0+bMOBpNOCi1+EJmEwmdqadRwVc0JNvVVF1NlC09eua3AhydXVK7FVJrtWKooBJp3f4vs1mMy4uLg7fr7gC83nISwGvZtpkwftwNDObNGsemPUoOpUG3m74GOX9qSi5ubns3LkTk8lERESEfB+qEHkvKp5FVdEBOkW55DryPlS+TZs2LVZVdaCz46jNJLmqYCkpKcTHxzs7jBonLi6O2NhYDmdmEGByxdtoZOHJY5itNjYln+fbE/shX0e4Tx1m9+6Kh6F2ftQ3nzvLLXGrUMxGLG75fNurJ+39Ahy2/8L3QVSyA5/AiYXQ+3cAli5fzttmhbPmZHLeiyF43C76Nw/grehOmPSOT6qFZs+ePbRq1YqwsDDi4+Pl+1CFyHtRsX4+doRHEjYAcHDoqEsmWPI+VD5FUQKdHUNtVzuvOEWNEe7pZX8+NLQBADc1CON4TiaL5roTP+ggx7KzaOnt46wQnWr8mjWcfDeazFUN8Ox3mPGs449+fajv7uHs0MS1aDQWNk4GmwV0BlRVZac5hbyEYHReZvJtNjafPSeJVQVTCi4oVVV1ciRCVK4hofU5lpmNQdFx6XKr0qWazaTk5XIoM4Pr64ZetuRLiOpIkitRY8w7nMgrW3eSpeRh0VlxrR/MDcFhPLp+M+6KCydys3gsqiU3Nmjo7FArxYGMdDLII/TuHVjv34JqyiddUTmanSXJVbWnA50RrDmg88Ko0/F1l17czkr69AGdJZC3ukY7O8gaT5KryqWqKjvTUjmWnUUDdw8ifP2cHVKtZdLrmdaqVbm325WWypCVf9unl/YZSJMiN0mFqAkkuRI1xpM7NpK1IZTzX0RhSXGj0XPrWXD6CLa9QaT82hD3zieZ63q01iRXK0+fBuDhTk3oHBhEpK8fmfn5+Ej99+pv71sQOhSMFy5KegQFs7zvIBp5eDoxsNpFkquKlWo2M2T5UlLMuZixAmBK88Ry3BuaJRPp58svfWOdG2QtsTX1HEuTktCpClNbtUJ/laVNTTy9+G/rdthsKh0DAiWxEjWSJFeixpgR1ZFH2Yia6kbKR9EceaUjxgYtMO/zJ3jqZtxij3Lc6u7sMCvNhKbNGNkwrFgyJYlVDWHJBt+IErMlsapcklxVrPR8M+kWM2aseOCCWbVi8czBWE8h97wrWR5WZ4dYa3y6+yB/pCQCMKFpU9afTebf02dxM+h5uE3rMidbJr2eu5s0r8BIhXA+Sa5EjWFQtJEFUj5uD4CaY8S8T+u8IfXvMHT10zjd/BxzDx1ibOPGTouzMkkyVUOZz4Fn7fgMV2WSXFWshh6ebB9yo31aVVUyLRaO52SxKy2VDv7Sbr+yvNkpmkdyWtLA3QOrqnJv/Fr7sntaNMNbeiUVwk6SK1EjHM3K5OUtOzn9YndQS95BazB5Oz7hOZy1qUzfuaXWJFeiBrKa4dhPcP0qZ0dS60lyVbkURcHLaKSV0ZdW3r7ODqdWMen19g6kDMBvPa8j32ajsaeXJFZCXESSK1GtZVss9F2yhNPWLLIWNSVnQ0gpa6ko/tm837kjBp2OA+kZlR6nEA6TthNMgeAtVWucTZIrUVu19fV3dgglLD11klyrFS+jkUgfP/xNJmeHJGopSa5EtaYCKNqFTd4pdyitU1gFVKMVF52eaP8AqUoiqrfU7eDb1tlRCC4kV0II58q2WLh74xr7dGuXQP4Y0MeJEYnaTOfsAIS4Fh4GA9cH1yP/iA/pC5vZ59d7NJ5GM+O0CRXUTCPHs7OcE6QQjpS2E3xKdmYhnEdKroRwLneDgaV9BjKuQXPq4cNdLZs4OyRRi0lyJaq9fKwYw9Iw1s0EwNgoFWOvw9A0mUa//4hbpySsgZl0Dazj5EiFcID03eBT/vFlhONJtUAhqo4mnl5Mj2rH6hv6c1NY7RhyRVRNklyJam/9yXOcndUeS7I7Ot9cGj6UQJ8gre2VurcO5n3+KFYdx3Ok5ErUAJZMMHo7OwqBJFdCCCFKkjZXoto7qqYSdOd26kzchlVvpb7Bh+XJyeSsaQApHjScs5BAnTuhbrVnjCtRgykGsFmcHYVAkishhBAlSXIlqr2fe/RFQaGplxeuOj06RWHe0cM8e+A4dE5i7YChBJlcnR2mEI6hM4Etz9lRCCS5EpVnb3oaryfs5nyemfe6d6Ceu9wsFKKqkuRKVHvt/QLsz789dJi3tu8hjRxyDJ60NHlLYiVqFp3R4SVXm86l4KY3kGHJp5N/oPSCV0aSXInK8taO3fyTdgyAEzlZklwJUYVJciVqlP9t30rmEU8M4Zn4NszlzW49nB2SEA5mA8WxzWVHrVkOgM6i5+GIVkxs3ByTXu/QY9REklyJijb/WCKLEpMY1aQBwxrVp6G7BxG+fs4OSwhxGZJciRojJS8Xi6qitxh4uW0MNzcMlzvwouax5IDezaG7HBUazvwN57DVT2PGnh0k5+WSmmmlQ1AAtzUJd+ixahJJrkRF+2j3Pg7lpbF803Fmtu0siZUQ1YD0FihqjA/37sWqWDE0PkfPoGBJrETNZDODzsWhu3TFCPXT7NNzDh/gt+TDPLMrnvvXbCTLIh1olEaSK1HRlvXvz4A69TEf9ebNbXucHY4QogwkuRI1Qo7FwpdH9nH+h5boLHrS8s3ODkmIiqFaQXFslb0XotvyfbdYeviEAhCu82N4Ha3EatG5RNYkn3bo8WoKuYEjKkMbX19cGqbTuU7AlVcWQjidVAsUNYKrXs89Ya35avQBshUz/9uyg697dceok/sHoobJTwO
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n",
    "request.addIdentifier('table', 'mapdata.lake')\n",
    "request.addIdentifier('geomField', 'the_geom')\n",
    "request.setParameters('name')\n",
    "\n",
    "# Get lake geometries\n",
    "response = DataAccessLayer.getGeometryData(request, [])\n",
    "lakes = np.array([])\n",
    "for ob in response:\n",
    "    lakes = np.append(lakes,ob.getGeometry())\n",
    "print(\"Using \" + str(len(lakes)) + \" lake MultiPolygons\")\n",
    "\n",
    "# Plot lakes\n",
    "for i, geom in enumerate(lakes):\n",
    "    cbounds = Polygon(geom)\n",
    "    intersection = cbounds.intersection\n",
    "    geoms = (intersection(geom)\n",
    "         for geom in lakes\n",
    "         if cbounds.intersects(geom))\n",
    "    shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), \n",
    "                        facecolor='blue', linestyle=\"-\",edgecolor='#20B2AA')\n",
    "    ax.add_feature(shape_feature)\n",
    "fig\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Major Rivers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 1400 river MultiLineStrings\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAKCCAYAAAA9RiPYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOzddXhUZ/bA8e8dj7snRAgBQgLBgmuLFbYspdSFsjVaSn3b7lZof+3W3bfeboW6UsO9aJBACAQCBEKI2ySj9/fHkBSKwyQTOZ/nyZPMXDujuee+73teRVVVhBBCCCGEEEKcHY2nAxBCCCGEEEKItkCSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDeQ5EoIIYQQQggh3ECSKyGEEEIIIYRwA0muhBBCCCGEEMINJLkSQgghhBBCCDfQeTqAtm7s2LFqSUlJk+zbqaoAaBTlrPZTa7ej12gwaDRU22z46fXuCO+Uj+2t07G3toYIL28sDgcWp4Mwo4k6h4MdNVWN6xrR4avXYdBq8AWMRiPZleWogAEtNoeKqnHCoacj1GAiysvrlOIwO+x4a9vmx6He4UBRwKjRun3fVqsVg8Hg9v2Kk7CWg6UE/Dq5bh56HfbUmKl0WMCuQQHi/L0I0Mvr01Tq6+vJzs7GaDSSlpZ2Sp8Hu6pyoM5MrLcPTlU96+9vcWyt4btpj7mGcKMXJq37v5ubg11V0XDic5DW8Dq0NWvXrv1VVdWxno6jPWubZ5MtSElJCWvWrGmSfX+2eyd7zLXcmNwF/1NMiCqsVqrtNpaXHOSiuATya2vIrqzAR6djeHgk169eziu9++NUVbx07n972J1OSiwWlpUU8ffYeK5auZj/9R+KcponGAsXLmT48OHsqqkmxGjCX6/nx/17sTqcrC0u55N928GmkBgQwfvDBuBzgsdyoK6O+zau4b99B6HXtK3G3HVlpVy6cAmKVY/dy8YnQ4fQMyjEbftveB1EM9vxX9j3Iwz7HoC5CxbwglWh1FpM3es9Cb9iK4MS/Xmz76Am+RwLl5ycHLp27Up8fDxr1qw55c+Dqqq8viOHYIORS+KTmj7Qdqgpv5tsTifLiovoHxrOf/O2MT25C9/v20OQwcjIiCie3rqJyxM6Eu3lfdS2FoeDmetW8mj33tTYbMT7+LbKBPvrvbu5M2sVAHkTLjzuY5D/Ec1PUZRQT8fQ3rWtM8l25pL4JHoHh/CfLRtOuN78okK+2psPwBUrFlFnt2O228mtruLJrZuYEBPHiIgoFEXhrcxBGDQarlq5mN21NRysr8PicLgl3k937+T7fXtRUam229AqCh8PGHbaidXhEn39GhPLCdFxXBAXz986xJDuF0zN3CTWmA+w11zbuL56qLWvQYXVypryEt7rN6TNJVYAU5ctY/9Lvdh2+XkUvN6dqYtXUHDY8yFaqYTLYf8ccNoB1/s621qMo9gLrb8Fu8bO6tJSdtXWeDjQtq3hu+uv3yunst0lHZI4NzKaz3bvbIrQxBk4WF+H3ekEXAkUQFF9HUV1ddicTj7Oz2NlyUHya2tYcPAAJq0WP50ep6rSKyiEUKMRgGQ/fyJNXnyyeyfziwoBWFhUyLWrllJtt3FrSjcC9AYSff1aZWIFMD46ltuS07irU3dO9xFUWK3sqK7it8J9jT1whGhL5JJmKzcgJJyR4VFsq6rkm4Ld3NM1nWdyNnNDcheKLfWgqqT6B7Lb7DrJ+mbIOeg1GlL8A1xXT/sMOGqfiqLwTr/B+OsNPLBxHeOj4+gfGnZWcdbZ7Zwf04H82mqivLyZmtjprPZ3LLN35fP4hmyqlXqcGiemWD0Twjswc8VqAjVebKgvItHPh1+Gj2H2nl1U2qxcHt+RrZUVjIqIxthKu2Ycz47qKqqxEH3dZhw3r0c12qhSVPaYa4n19vF0eOKsaECjB0cdaPzQazR81H8oV65czPgZpRw0+zI5KY7UgEBPB9qmnWlyBRBsNFJsqcfssFNhtVJUX0dn/wB3h9gqlVjqWV1awtioGPJqqjFptcR6+/B74T78dHryzTU4VZVL4pNOKTlxqCoLiwo5JzKaH/ftJdLLi95BIew215Lg48vsPbv4e0wHntuWzZjIGMKMJu7ftI6vB4/krbxcErx9uDg+iQ0VZaT6B5Ls68e/u/UA4Jok1/+yRF+/xuNNio0HYHBoODV2OyWWer7bt4ebkrsSYjASajQ1wbPWvIxaLbd27Xra222prGD84t8bb88dMZaOhz13QrQFkly1ct6Huvz46fWcH9MBq9NJnLcPvjod84vK0CgK58d0IPLQ2KPDW2dO1GLkf2icxj2p6egVDS/lbsGk0RDj7cN5UbG8sn0rNyR3odJmxajR4KPTY3E4GuM5nKqqjF/8O18NHkm3gCB3Pvwj3Lt5NbWroil/NwNHrY74Jxbx48Ei1BoDxa+m4DNYIXm4BYBzI6JRUfHR6bi7a3qTxeRJi4uKALgjsyP9QsNIDwyixmYjQPq/t37bnofoCaD/86RkcFgEC0aOI8bLm82V5cfskiTc62ySK4Awo4lpSSmsLDnIspKDdDzUktFaWzNOxcqSYgINBsKNJu7buIbX+gzkXxvWMDgsgkiTN78eKGBGp1TmHthP3+BQpixZSLXDigPXc2yo9MFe6oUzvpS3c3fw8znnYnE6j+gan3+oxXZ/nZkqmw2z3c7TOZsZERFFkMGAr05PTnUlL2zL5tXeAyiqr6PKbuOJHn0A1+v5zeCRKIrC/YeSKICnMvo2/m04hdeog49v49/P9+p3Vs+bp22oKGNuYSEaVWFm165oz/A92tHXj3+l9sDpVOkbEiqJlWiTJLlqI6K9vBtPphr68P/90NWzs+Grc/3DivP2obNfAL46HZU2G746PQaNhqe3buLCuAT2mGvZVFHOfanduWnNct7OHMx9G1xjzZ7M6MsPQ0edcNyTOzyT0Ze7WI1a4UX14ljsPhaq/pdK1eddiZi5Hq9B+yjD9RyFHOq+0ZZNS+7E5A7xRyRTkli1EXYzBKYddXeCjy/3b1xLncPBsz0zPRBY+3K2yVWD/qHh9A8N59LlC/lHUgrnRkaf8rZO1dUaneDji6qqZ9XNuilYHA50Gg03rl7OnV3SUFHJq6mio68fNyZ3QasozErriUGrZa+5lqsTOxFoMPBsr0z21NZgV504UPHBgFV14PCtx6BqqC/zRmvUMqewgE0V5fyzazq3rlvJS73
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n",
    "request.addIdentifier('table', 'mapdata.majorrivers')\n",
    "request.addIdentifier('geomField', 'the_geom')\n",
    "request.setParameters('pname')\n",
    "rivers = DataAccessLayer.getGeometryData(request, [])\n",
    "print(\"Using \" + str(len(rivers)) + \" river MultiLineStrings\")\n",
    "\n",
    "# Plot rivers\n",
    "for ob in rivers:\n",
    "    shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), \n",
    "                        facecolor='none', linestyle=\":\",edgecolor='#20B2AA')\n",
    "    ax.add_feature(shape_feature)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topography\n",
    "\n",
    "Spatial envelopes are required for topo requests, which can become slow to download and render for large (CONUS) maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[<awips.dataaccess.PyGridData.PyGridData object at 0x7ffd0f33c040>]\n",
      "Number of grid records: 1\n",
      "Sample grid data shape:\n",
      "(778, 1058)\n",
      "\n",
      "Sample grid data:\n",
      "[[1694. 1693. 1688. ...  757.  761.  762.]\n",
      " [1701. 1701. 1701. ...  758.  760.  762.]\n",
      " [1703. 1703. 1703. ...  760.  761.  762.]\n",
      " ...\n",
      " [1767. 1741. 1706. ...  769.  762.  768.]\n",
      " [1767. 1746. 1716. ...  775.  765.  761.]\n",
      " [1781. 1753. 1730. ...  766.  762.  759.]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy.ma as ma\n",
    "request = DataAccessLayer.newDataRequest(\"topo\")\n",
    "request.addIdentifier(\"group\", \"/\")\n",
    "request.addIdentifier(\"dataset\", \"full\")\n",
    "request.setEnvelope(envelope)\n",
    "gridData = DataAccessLayer.getGridData(request)\n",
    "print(gridData)\n",
    "print(\"Number of grid records: \" + str(len(gridData)))\n",
    "print(\"Sample grid data shape:\\n\" + str(gridData[0].getRawData().shape) + \"\\n\")\n",
    "print(\"Sample grid data:\\n\" + str(gridData[0].getRawData()) + \"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "623.0\n",
      "4328.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAKZCAYAAAD5854AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOy9d5gcV5n9/7m3qjpPntaMsiznIDnjbDAmB5NMXDILLDktsKQFlpxhgSVnWGDBC6xhSYuzDTbOOduyJGs0PXk6d9W9vz8qdHV19cxIMgK+vzrP0890V7j3Vpjuc98673mF1poECRIkSJAgQYIECRLsH8i/9gASJEiQIEGCBAkSJPj/ExICniBBggQJEiRIkCDBfkRCwBMkSJAgQYIECRIk2I9ICHiCBAkSJEiQIEGCBPsRCQFPkCBBggQJEiRIkGA/IiHgCRIkSJAgQYIECRLsRyQEPEGCBAkSJEiQIEGC/YiEgCdIkCBBggQJEiRIsB+REPAECRIkSJAgQYIECfYjEgKeIEGCBAkSJEiQIMF+RELAEyRIkCBBggQJEiTYj0gIeIIECRIkSJAgQYIE+xEJAU+QIEGCBAkSJEiQYD8iIeAJEiRIkCBBggQJEuxHJAQ8QYIECRIkSJAgQYL9iISAJ0iQIEGCBAkSJEiwH5EQ8AQJEiRIkCBBggQJ9iMSAp4gQYIECRIkSJAgwX5EQsATJEiQIEGCBAkSJNiPSAh4ggQJEiRIkCBBggT7EQkBT5AgQYIECRIkSJBgPyIh4AkSJEiQIEGCBAkS7EckBDxBggQJEiRIkCBBgv2IhIAnSJAgQYIECRIkSLAfkRDwBAkSJEiQIEGCBAn2IxICniBBggQJEiRIkCDBfkRCwBMkSJAgQYIECRIk2I9ICHiCBAkSJEiQIEGCBPsRCQFPkCBBggQJEiRIkGA/IiHgCRIkSJAgQYIECRLsRyQEPEGCBAkSJEiQIEGC/YiEgCdIkCBBggQJEiRIsB+REPAECRIkSJAgQYIECfYjEgKeIEGCBAkSJEiQIMF+RELAEyRIkCBBggQJEiTYj0gIeIIECRIkSJAgQYIE+xEJAU+QIEGCBAkSJEiQYD8iIeAJEiRIkCBBggQJEuxHJAQ8QYIECRIkSJAgQYL9iISAJ0iQIEGCBAkSJEiwH5EQ8AQJEiRIkCBBggQJ9iMSAp4gQYIECRIkSJAgwX5EQsATJEiQIEGCBAkSJNiPSAh4ggQJEiRIkCBBggT7EQkBT5AgQYIECRIkSJBgPyIh4AkSJEiQIEGCBAkS7EckBDxBggQJEiRIkCBBgv2IhIAnSJAgQYIECRIkSLAfkRDwBAkSJEiQIEGCBAn2IxICniBBggQJEiRIkCDBfkRCwBMkSJAgQYIECRIk2I9ICHiCBAkSJEiQIEGCBPsRCQFPkCBBggQJEiRIkGA/IiHgCRIkSJAgQYIECRLsRyQEPEGCBAkSJEiQIEGC/YiEgCdIkCBBggQJEiRIsB9h/rUH8PeGRz32UXp6arpjmYPq+GwsM6/xt1daI4Xo2q+lbQCkAKWVt71AIjGEgUTS0i0ayiEjLUxh9OyrpVuARqGQaAA0IDBQaDSQEWkAbO0ABO01lI0hBJP1BgOWhQYWWi2G0gYLTYeZpt3Rl4XBkCHAMpisNxBKoKUO1kstWFfIkDEkjlY9x11XLbQGU0oMJA4KA0lL2whAColCkRIWtnYwhERrUCgcrTCExBBuHwCGaF8Pgejqz9EOhjcWQ0iUt5/AbRMUhvevorxrJ71rpVGI0PW2tcYMXVPptdu5nSI699XYHe2AotlwSKctNAqtHQQCrRwQJmgNAoQwQKv2ldUKIQVaCNASvGMXSJRuuu+Feyxaa6QwUN51F2gcWggMNALtHacOjcrRChsnOP7oelMYtIL2wNYKyxuDf98bSExhYGun46+/H5H2wvC3sbzl7X3a59xfZ9R2AYp6erV7RlstSrZN2bFRFRMj57AunyVlSPz/Cnd/GbRte/dCVlpdY+vs3+3X/+ze3+3r2VAtUqE2rK7j6vwO8cfhL7eEQUt1bwNgyfZ447bz1+8Lpqem2PnAdgaGh1mzYT3YNlYqFbttS2nmmg2yhknaMLBk9/+cu117nOExtpSOfFZLtKFjl7ttimXWS7cvIWjp7u0sIWip6D7dbYaX+e+j4w2v7zqGZgMrle45zvAYrL24lNFj6NVWUzlIBHVHYUpJxujuzG/LDp2DrBl/bfYFLbWXx9rjXO7NOWw5ve8dAFuBGWrLMpY/Dy1HYxkitm1baUwpgvHt63UP2o18rZpGe5kZ+hq68fprf6u1ftze95RgT5EQ8D3E5FSJX13xa/qMQrBs0SkDdCxbDttqM7zvhjt4z1Fb2VTIs+AsAjBo9LG7tZu6qlFT8zRVA1MM49j9bFu0OXZ4iPvLFY4bHuLdN17Pcw9YS9NRLDhVThxeBcCQ2R/047c7b09TdUqkZAuNRV3lKRhD9Bl5Bs2+YPuG47CtUuXgvgJPvOhivnPacWyvVBlMWQynU+xulQBYk1rFnYszpESGvMyxIZ/DlJLrL7uaY04/gf/duYuD+grcW57l0gfn2FVucGNtFhs498BVPGHDAEWryJztnrtBs0DZqVBvWbz7xht501HrGLHyAPQbBQwB806ZYbMPA0HZqVAw8sw7ZWzvx7OqqvQbeQpGHkdrFp1KsNxtx21vwLtO806Zgrds0Ohj0SlTVRUMAUo3yMksOZnHoUpa5EnJPHVVpqXdzz4y0n3/w2338G833o1E8Jz1a3nr1k3BOltXsWSu5/3gaHeMpsihvPeXX3Ijp525FaVnQdUx9ACtRhUhsxhmDp+jKAVaVdG6jtYNREqgZSbSg/vZpg64/N3RaQyRDW0zR0ZmKNuTNDBxdIqW1hSMERwNZVWmrjQLToWW1iw6Ne9+G2bQLOCESMysd11HrT7m7DKmEAx517F9zN74vc/DZh8lewGAEbOfydY8Rat9LwPsbrrrV1n9TLTce3vcGmCiWWHccu/jiYY7rg03voVmdi0TB7+ZMSvH//3hIt6nqxhXriZ/0i6eunYdLzhkEwibUmuRotXX0de4NeC215oP2nY/+/26yyZbC6yy+plsuWObtBcxEBQt9z4rtcreeerrbLdZifTX3UfsOu/4gnXpbNfy8VQ+6MNfvy/4z29+m7e++nU84gmP41+/+Fmmr7mJI085uef2/73jAZ68ZgPTrYY3xniCGT4H46k8Ew1/+/aYb1qcpZiK3s9umxONRuzxdZwL//zUWu1lWatj2Ximx/hqylvfPQGbqPZgtsugVHcoZgxvHAa3XHkJR550Zqhdb4yhr4tgWXZ5kre7XuU799/F89cfxR2VXazPjFJzbK6d386LNh2GozVGKEigtRuI+cWD99KwDR5R3Mh4Lp7x9TrmuO0najp2zL2Wr2Sb87bfw6fvvo5zRg7m7VuP7trvlqsu5siHPXzJcZdqmmJWMJ7r7n+y3MlWS9VOolz09ilVdfA+un0xJxjLr4wx7160Y5ePFdz7Y3doPP6yPUFpod3+qr7u/ScXnWDd6sHU6B53kGCfkBDwPYQfTV10ygHhXop4++Tch7/txuwwnznxKH5w7wOcOTaCKQTzrRaHDUBLa4bNUb509wNsGTYZNHL8Zsckhw/0I4B7yy7p/MSxxwNw3cwsC40aC60WX7lzG286fDOmlAyZ/ZSdCnP2DH1GDkQKRR8CSMsUphDsbpVYdCoBEd9ZrfGD++/n/Vu38O3TjsWSgnymStbIM2MvoHSDPiNH2alwwuBGZu0FQFNWFVBuxHPOXuTUMfc4V+WGOHKoj9/tepBb7pijUcrwO2ORE1Yb3DJfZX2mj5pSzKcrXLJrgTNWFTn3gCIjVh5DuOTbJ9IpISk7FSwh6Dfz1FQFR9cRQEZmyUq3z7K3fUvX6DfyjJhjHdeg32xfL5+EzzkuafLJtyUIyLlPFJuqgkSQFvmAWPuYaS7wgRvvZteHTkGYiu+/448c1J/luZs3L3dLAWCIHI6uYusqpmj/8ppC4CjQOo0QWQxTI2QOraogc/hBxGajipQNUrk8Wma
      "text/plain": [
       "<Figure size 864x864 with 2 Axes>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid=gridData[0]\n",
    "topo=ma.masked_invalid(grid.getRawData()) \n",
    "lons, lats = grid.getLatLonCoords()\n",
    "print(topo.min()) # minimum elevation in our domain (meters)\n",
    "print(topo.max()) # maximum elevation in our domain (meters)\n",
    "\n",
    "# Plot topography\n",
    "cs = ax.contourf(lons, lats, topo, 80, cmap=plt.get_cmap('terrain'),alpha=0.1, extend='both')\n",
    "cbar = fig.colorbar(cs, shrink=0.5, orientation='horizontal')\n",
    "cbar.set_label(\"topography height in meters\")\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}