{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overlays\n",
    "\n",
    "Spatial overlays allow you to compare two GeoDataFrames containing polygon or multipolygon geometries \n",
    "and create a new GeoDataFrame with the new geometries representing the spatial combination *and*\n",
    "merged properties. This allows you to answer questions like\n",
    "\n",
    "> What are the demographics of the census tracts within 1000 ft of the highway?\n",
    "\n",
    "The basic idea is demonstrated by the graphic below but keep in mind that overlays operate at the dataframe level, \n",
    "not on individual geometries, and the properties from both are retained\n",
    "\n",
    "![illustration](http://docs.qgis.org/testing/en/_images/overlay_operations.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can load up two GeoDataFrames containing (multi)polygon geometries..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from shapely.geometry import Point\n",
    "from geopandas import GeoDataFrame, read_file\n",
    "from geopandas.tools import overlay\n",
    "from geodatasets import get_path\n",
    "\n",
    "# NYC Boros\n",
    "zippath = get_path(\"nybb\")\n",
    "polydf = read_file(zippath)\n",
    "\n",
    "# Generate some circles\n",
    "b = [int(x) for x in polydf.total_bounds]\n",
    "N = 10\n",
    "polydf2 = GeoDataFrame(\n",
    "    [\n",
    "        {\"geometry\": Point(x, y).buffer(10000), \"value1\": x + y, \"value2\": x - y}\n",
    "        for x, y in zip(\n",
    "            range(b[0], b[2], int((b[2] - b[0]) / N)),\n",
    "            range(b[1], b[3], int((b[3] - b[1]) / N)),\n",
    "        )\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first dataframe contains multipolygons of the NYC boros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polydf.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the second GeoDataFrame is a sequentially generated set of circles in the same geographic space. We'll plot these with a [different color palette](https://matplotlib.org/examples/color/colormaps_reference.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polydf2.plot(cmap=\"tab20b\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `geopandas.tools.overlay` function takes three arguments:\n",
    "\n",
    "* df1\n",
    "* df2\n",
    "* how\n",
    "\n",
    "Where `how` can be one of:\n",
    "\n",
    "    ['intersection',\n",
    "    'union',\n",
    "    'identity',\n",
    "    'symmetric_difference',\n",
    "    'difference']\n",
    "\n",
    "So let's identify the areas (and attributes) where both dataframes intersect using the `overlay` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf = polydf.overlay(polydf2, how=\"intersection\")\n",
    "newdf.plot(cmap=\"tab20b\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And take a look at the attributes; we see that the attributes from both of the original GeoDataFrames are retained. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polydf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polydf2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the other `how` operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf = polydf.overlay(polydf2, how=\"union\")\n",
    "newdf.plot(cmap=\"tab20b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf = polydf.overlay(polydf2, how=\"identity\")\n",
    "newdf.plot(cmap=\"tab20b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "newdf = polydf.overlay(polydf2, how=\"symmetric_difference\")\n",
    "newdf.plot(cmap=\"tab20b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf = polydf.overlay(polydf2, how=\"difference\")\n",
    "newdf.plot(cmap=\"tab20b\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}