3-geopandas_advanced.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GeoPandas: Advanced topics\n",
    "[Emilio Mayorga, University of Washington](https://github.com/emiliom/). 2018-9-9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We covered the basics of GeoPandas in the previous episode and notebook. Here, we'll extend that introduction to illustrate additional aspects of GeoPandas and its interactions with other Python libraries, covering fancier mapping, analysis (unitary and binary spatial operators), raster zonal stats + GeoPandas. \n",
    "\n",
    "**Here are the main sections in this episode / notebook:**\n",
    "- Read HydroBASINS for Western Washington from a PostGIS / PostgreSQL relational database on the cloud\n",
    "- Dissolve into larger watersheds, and reproject\n",
    "- Plot `choropleth` map based on calculated watershed areas\n",
    "- Choropleth map as an interactive map with folium\n",
    "- Spatial join, `sjoin`, of polygons on points\n",
    "- rasterstats: \"zonal\" statistics from polygons on rasters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Set up packages and data file path\n",
    "We'll use these throughout the rest of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import (absolute_import, division, print_function)\n",
    "import os\n",
    "import json\n",
    "import psycopg2\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from shapely.geometry import Point\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from geopandas import GeoSeries, GeoDataFrame\n",
    "\n",
    "data_pth = \"./data\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.__version__, pd.__version__, gpd.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(data_pth, \"db.json\")) as f:\n",
    "    db_conn_dict = json.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fix for pypro module\n",
    "envName = 'my_env'\n",
    "\n",
    "#Fix the issue with pyproj\n",
    "import os\n",
    "localappPth = os.environ.get('localappdata', '')\n",
    "os.environ[\"PROJ_LIB\"] = \"{}\\\\ESRI\\\\conda\\\\envs\\\\{}\\\\Library\\\\share\".format(localappPth,envName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Read HydroBASINS North America dataset, extracting Western Washington\n",
    "Read [HydroBASINS](http://hydrosheds.org/page/hydrobasins) \"all-levels\" (lev00) hierarchical watersheds dataset for North America and the Caribbean (`hybas_na_lev00_v1c`), from Amazon Cloud PostgreSQL/PostGIS database. **Watersheds in the dataset are at the finest (highest resolution) \"Pfastetter\" hierarchical level, level 12.** HydroBASINS dataset technical documentation is [here](http://hydrosheds.org/images/inpages/HydroBASINS_TechDoc_v1c.pdf).\n",
    "\n",
    "`read_postgis` is called as before, except now we'll apply a SQL filter (server side) to the PostGIS dataset to select only the Pfastetter level-4 watershed with code 7831: `WHERE pfaf_4 = 7831`. This is **most of Western Washington.** Watershed polygons will still be read at their original level 12 resolution.\n",
    "\n",
    "For a more in-depth look at interacting with spatial relational databases, see the eScience Institute tutorial [Introduction to SQL and Geospatial Data Processing](https://uwescience.github.io/SQL-geospatial-tutorial/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn = psycopg2.connect(**db_conn_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww = gpd.read_postgis(\n",
    "    \"SELECT * FROM hybas_na_lev00_v1c WHERE pfaf_4 = 7831\", conn, \n",
    "    geom_col='polygongeom',\n",
    "    coerce_float=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(hydrobas_ww)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "413 polygon features returned. Let's examine the attributes available, using the first feature as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww.iloc[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot a categorical map with coloring based on the aggregating column `pfaf_7`.** Watershed *boundaries* are at the high-resolution Pfastetter level 12.   \n",
    "Note: pick a color map (`cmap`) appropriate for your data. [Get to know the matplotlib color maps.](https://matplotlib.org/tutorials/colors/colormaps.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww.plot(column='pfaf_7', cmap='tab20', categorical=True, figsize=(14, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rename the GeoDataFrame geometry column from `polygongeom` to `geometry` to avoid issues with other packages\n",
    "Unfortunately, `folium` choropleth and `rasterstats` (demonstrated below) require the geometry column to be named \"geometry\". So, we'll rename it here first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww = hydrobas_ww.rename(columns={'polygongeom': 'geometry'})\n",
    "hydrobas_ww._geometry_column_name = 'geometry'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Dissolve into larger watersheds, and reproject"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dissolve source polygons into larger watersheds based on attribute values\n",
    "Apply GeoDataFrame [dissolve](http://geopandas.org/aggregation_with_dissolve.html) aggregation method (implemented from lower level `shapely` operators) on level-7 Pfastetter codes (`pfaf_7`) shown in the plot above. Aggregate attributes, retaining only `pfaf_7` and `pfaf_6` (plus `geometry`, of course). This operation results in only 17 polygons, from the original 413."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = ['pfaf_6', 'pfaf_7', 'geometry']\n",
    "hydrobas_ww_p7 = hydrobas_ww[cols].dissolve(by='pfaf_7', aggfunc='first', as_index=False)\n",
    "len(hydrobas_ww_p7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's examine some of the features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results. Looks like the previous plot, except the polygon boundaries are now the pfaf_7 watersheds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7.plot(column='pfaf_7', cmap='tab20', categorical=True, edgecolor='white',\n",
    "                    figsize=(14, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_NOTE/WATCH:_**   \n",
    "**Beware that `dissolve` may fail if there are \"invalid\" geometries.**\n",
    "This code is based on a GeoDataFrame examined in the previous, intro notebook. The 6 geometries/points reported are invalid (and are reported by the `is_valid()` method). This dissolve statement does work, though.\n",
    "```python\n",
    "seas_grp = seas[['oceans', 'geometry']]\n",
    "seas_oceans_diss = seas_grp[seas_grp.geometry.is_valid].dissolve(by='oceans')\n",
    "\n",
    "Ring Self-intersection at or near point 10.407218181818182 54.821390909090908\n",
    "Self-intersection at or near point -79.365827272727287 76.296645454545455\n",
    "Ring Self-intersection at or near point 10.979445510225332 54.380555030408686\n",
    "Ring Self-intersection at or near point 133.61550925464189 -4.3005540903175188\n",
    "Ring Self-intersection at or near point 121.91067196634913 -5.0593090510592447\n",
    "Ring Self-intersection at or near point 115.29553592754269 -7.0082630551828515\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reproject (transform) to WA State Plane South, epsg:2927\n",
    "Partly so we can calculate polygon areas in linear units, not geodetic degrees. But also because that's the projection used by most state and local governments in Washington.\n",
    "- http://epsg.io/2927\n",
    "- http://spatialreference.org/ref/epsg/2927/\n",
    "- [Report from http://www.epsg-registry.org](http://www.epsg-registry.org/report.htm?type=selection&entity=urn:ogc:def:crs:EPSG::2927&reportDetail=short&style=urn:uuid:report-style:default-with-code&style_name=OGP%20Default%20With%20Code&title=EPSG:2927)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No need to go to a web site to learn more about what `epsg:2927` is. Use `pyepsg`, which issues queries to http://epsg.io web services."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyepsg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the epsg code from the string returned by `crs['init']`, then query epsg `2927`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7_epsg_str = hydrobas_ww_p7.crs['init'].split(':')[1]\n",
    "pyepsg.get(hydrobas_ww_p7_epsg_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyepsg.get('2927')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Apply the crs transformation (reprojection)** using `to_crs` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7.crs = {'init' :'epsg:4326'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7_wasp = hydrobas_ww_p7.to_crs(epsg=2927)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot the reprojected map.** Note that, being in a planar project (not geodetic), the shape looks different compared to the previous map. More \"normal\". And the axes are now in `feet` relative to some origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7_wasp.plot(column='pfaf_7', cmap='tab20', categorical=True, edgecolor='white',\n",
    "                         figsize=(14, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Plot `choropleth` map based on calculated watershed areas\n",
    "As the projection is in `feet`, auto-calculated polygon areas will be in feet<sup>2</sup>. So let's convert to miles<sup>2</sup> first (why not!). We'll add a new column to the GeoDataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7_wasp['area_mi2'] = hydrobas_ww_p7_wasp.geometry.area / 27878400\n",
    "hydrobas_ww_p7_wasp.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_NOTE/FUN:_**  \n",
    "Now you could get the area of a pfaf_6 watershed via simple Pandas DataFrame `groupby` aggregation (sum)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the choloropleth, using `area_mi2`.\n",
    "\n",
    "The \"fisher_jenks\" value segmentation `scheme` (using 7 segments, k=7) used is one of the available `pysal.esda.mapclassify.Map_Classifier` classifiers from the powerful [PySAL package](http://pysal.org/) (Python Spatial Analysis Library); GeoPandas can use these classifiers if PySAL is installed, as it is here. To get the list of classifiers, use:\n",
    "```python\n",
    "import pysal\n",
    "print(pysal.esda.mapclassify.Map_Classifier.__doc__)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, figsize=(10, 6))\n",
    "ax.set_title('Watersheds by area ($mi^2$)')\n",
    "hydrobas_ww_p7_wasp.plot(column='area_mi2', scheme='fisher_jenks', k=7, \n",
    "                         cmap=plt.cm.Blues, legend=True, ax=ax)\n",
    "ax.set_axis_off()\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_NOTE/FUN:_**  \n",
    "Let's stop for a bit to explore on your own, hack with your neighbors, ask questions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Choropleth map as an interactive map with folium\n",
    "[Folium](https://github.com/python-visualization/folium) is very cool, specially for use in Jupyter notebooks; or to export into stand-alone HTML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folium.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`m.choropleth` internally splits the geometry from the other attributes in `hydrobas_ww_p7_wasp`, and rejoins them based on the key `pfaf_7`. `key_on` uses an attribute reference derived from GeoJSON representations; this is awkward, and hopefully will be simplified in future folium implementations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrobas_ww_p7_wasp.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47.8, -122.5], zoom_start=7, tiles=\"cartodbpositron\")\n",
    "\n",
    "m.choropleth(\n",
    "    geo_data=hydrobas_ww_p7_wasp,\n",
    "    data=hydrobas_ww_p7_wasp,\n",
    "    columns=['pfaf_7', 'area_mi2'],\n",
    "    key_on='feature.properties.pfaf_7',\n",
    "    legend_name='Area (sq mi)', \n",
    "    fill_color='YlGn',\n",
    "    fill_opacity=0.4,\n",
    "    highlight=True)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This map is interactive, so play with it (zoom and pan). There is a lot more to explore in Folium! This is just a teaser."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Spatial join, `sjoin`, of polygons on points\n",
    "We'll use an old, local snapshot of NANOOS coastal and marine monitoring stations in the Pacific NW, from the [NANOOS Visualization System (NVS) Data Explorer](http://nvs.nanoos.org/Explorer). While many stations are moorings on marine waters, some are inshore or in tidal shores and will overlap the watershed boundaries. The point file is in the [GeoPackage](http://www.geopackage.org) format, an OGC format implemented in SQLite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanoosstations_gdf = gpd.read_file(os.path.join(data_pth, \"nanoos_nvs.gpkg\"))\n",
    "len(nanoosstations_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanoosstations_gdf.iloc[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points are on the coasts of the Pacific NW (BC, WA, OR) and out in the open ocean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanoosstations_gdf.plot(markersize=15);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Apply \"inner\" spatial join with the `sjoin` operator**. An inner join will retain only overlapping features. Then plot as a map overlay on top of `hydrobas_ww_p7`, categorizing (coloring) each point by the `pfaf_6` watershed it's in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nanoossta_hydrobas_ww_gdf = gpd.sjoin(nanoosstations_gdf, hydrobas_ww_p7, how=\"inner\")\n",
    "len(nanoossta_hydrobas_ww_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, figsize=(10, 6))\n",
    "ax.set_axis_off()\n",
    "plt.axis('equal')\n",
    "hydrobas_ww_p7.plot(ax=ax, cmap='Greys_r', linewidth=0.5, edgecolor='red')\n",
    "nanoossta_hydrobas_ww_gdf.plot(column='pfaf_6', markersize=30, \n",
    "                               categorical=True, legend=True, ax=ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. rasterstats: \"zonal\" statistics from polygons on rasters\n",
    "We'll end by mixing features from a GeoDataFrame with a raster, applying zonal statistics using the cool and light weight [rasterstats](https://github.com/perrygeo/python-rasterstats) package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monthly Juy long-term climatology precipitation. The original monthly time series data are from the [PRISM Climate Group](http://prism.oregonstate.edu); the monthly climatology and Pacific NW clip were created by your truly and Don Setiawan for the [BiGCZ project](http://bigcz.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppt_july_tif_pth = os.path.join(data_pth, 'prism_precipitation_july_climatology.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rasterstas` uses [rasterio](https://mapbox.github.io/rasterio) to read rasters (and `fiona` to read vector datasets), so we'll first do a quick exploration of rasterio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import rasterio.plot as rioplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterio.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppt_july = rasterio.open(ppt_july_tif_pth)\n",
    "ppt_july"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the metadata read from the raster file (we can confirm CRS is epsg:4326), then plot the raster. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppt_july.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rioplot.show(ppt_july, with_bounds=True, cmap=plt.cm.Blues);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply rasterstas `zonal_stats`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply `zonal_stats` from `rasterstats` package. Can pass a `GeoDataFrame` directly (instead of the file path to a GIS file) because it implements our old friend, the `__geo_interface__` method. For the raster, we pass its file path.\n",
    "\n",
    "`zonal_stats` returns a geojson with the original properties plus the zonal statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterstats as rs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zonal_ppt_gjson = rs.zonal_stats(hydrobas_ww_p7, ppt_july_tif_pth, prefix='pptjuly_',\n",
    "                                 geojson_out=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(zonal_ppt_gjson), len(zonal_ppt_gjson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zonal_ppt_gdf = GeoDataFrame.from_features(zonal_ppt_gjson)\n",
    "zonal_ppt_gdf.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### And finally, a choropleth map of July precipitation by watershed! With a good bit of plot tweaking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, figsize=(8, 6))\n",
    "ax.set_title('Mean July precipitation ($mm/month$) by watershed')\n",
    "zonal_ppt_gdf.plot(ax=ax, column='pptjuly_mean', scheme='Equal_Interval', k=5, \n",
    "                   cmap=plt.cm.Blues, linewidth=1, edgecolor='black'), \n",
    "#                   legend=True, legend_kwds={'loc': 'upper left'})\n",
    "ax.set_facecolor(\"lightgray\")\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Other resources, tools, and overlap with other tutorials\n",
    "- Advanced vector geospatial analytics\n",
    "  - Go over all the analytics that have been wrapped into GeoPandas, including [Geometric Manipulations](http://geopandas.org/geometric_manipulations.html), [Set-Operations with Overlay](http://geopandas.org/set_operations.html), [Aggregation with Dissolve](http://geopandas.org/aggregation_with_dissolve.html) and [Merging Data](http://geopandas.org/aggregation_with_dissolve.html)\n",
    "  - Use [PySAL](http://pysal.org), the Python Spatial Analysis Library! It's in the conda environment (version 1.14.4). This is a powerful, multi-faceted package. But watch out for the current transition from \"version 1\" to \"version 2\".\n",
    "- Visualizations\n",
    "  - [GeoPlot](https://residentmario.github.io/geoplot/index.html) makes it easy to generate a variety of useful plots from GeoDataFrames. It's available in the conda environment. [Here is a gallery of GeoPlot plots from GeoPandas](http://geopandas.org/gallery/plotting_with_geoplot.html)\n",
    "  - More spatial visualization options are coming in the [Visualization tutorial](https://geohackweek.github.io/visualization/). Stay tuned, or look it up now. The mapping packages it'll cover include [cartopy](https://scitools.org.uk/cartopy/docs/latest/), [GeoViews](http://geoviews.org/) and [Folium](http://python-visualization.github.io/folium/) (we only covered a small subset of Folium capabilities, just to give you a taste)\n",
    "  - See the [Python mapping libraries tour](https://github.com/oceanhackweek/ohw2018_tutorials/tree/master/day3/geospatial_and_mapping_tools) from OceanHackWeek 2018 last month.\n",
    "- Overlap with raster processing\n",
    "  - We illustrated `rasterstats` and [rasterio](https://rasterio.readthedocs.io/en/latest/). `rasterio` will be a pretty important component of your raster handling and manipulation toolbox. And it interacts with the GeoJSON-like objects we've examined; for example, see its [features module](https://rasterio.readthedocs.io/en/latest/topics/features.html).\n",
    "  - [regionmask](https://regionmask.readthedocs.io/) works nicely with GeoDataFrames to support gridded operations, including ones in `xarray` that you'll see in the [nDarrays tutorial](https://geohackweek.github.io/nDarrays/)"
   ]
  }
 ],
 "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": 2
}