Spatial Autocorrelation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Take-notice!\" data-toc-modified-id=\"Take-notice!-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Take notice!</a></span></li><li><span><a href=\"#Spatial-Autocorrelation\" data-toc-modified-id=\"Spatial-Autocorrelation-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Spatial Autocorrelation</a></span><ul class=\"toc-item\"><li><span><a href=\"#Methodology\" data-toc-modified-id=\"Methodology-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Methodology</a></span></li><li><span><a href=\"#Libraries-to-use\" data-toc-modified-id=\"Libraries-to-use-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Libraries to use</a></span></li><li><span><a href=\"#Libraries-to-use\" data-toc-modified-id=\"Libraries-to-use-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Libraries to use</a></span></li><li><span><a href=\"#Data-preparation\" data-toc-modified-id=\"Data-preparation-2.4\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>Data preparation</a></span></li><li><span><a href=\"#Block-Groups\" data-toc-modified-id=\"Block-Groups-2.5\"><span class=\"toc-item-num\">2.5&nbsp;&nbsp;</span>Block Groups</a></span></li><li><span><a href=\"#Data-cleanup\" data-toc-modified-id=\"Data-cleanup-2.6\"><span class=\"toc-item-num\">2.6&nbsp;&nbsp;</span>Data cleanup</a></span><ul class=\"toc-item\"><li><span><a href=\"#Subsetting-the-data\" data-toc-modified-id=\"Subsetting-the-data-2.6.1\"><span class=\"toc-item-num\">2.6.1&nbsp;&nbsp;</span>Subsetting the data</a></span></li></ul></li><li><span><a href=\"#Map-the-census-block-groups\" data-toc-modified-id=\"Map-the-census-block-groups-2.7\"><span class=\"toc-item-num\">2.7&nbsp;&nbsp;</span>Map the census block groups</a></span><ul class=\"toc-item\"><li><span><a href=\"#Using-Matplotlib-subplots\" data-toc-modified-id=\"Using-Matplotlib-subplots-2.7.1\"><span class=\"toc-item-num\">2.7.1&nbsp;&nbsp;</span>Using Matplotlib subplots</a></span></li></ul></li><li><span><a href=\"#Get-Arrest-Data-from-LA-Open-Data-Portal\" data-toc-modified-id=\"Get-Arrest-Data-from-LA-Open-Data-Portal-2.8\"><span class=\"toc-item-num\">2.8&nbsp;&nbsp;</span>Get Arrest Data from LA Open Data Portal</a></span><ul class=\"toc-item\"><li><span><a href=\"#Convert-data-to-a-geodataframe\" data-toc-modified-id=\"Convert-data-to-a-geodataframe-2.8.1\"><span class=\"toc-item-num\">2.8.1&nbsp;&nbsp;</span>Convert data to a geodataframe</a></span></li><li><span><a href=\"#The-0,0-conundrum\" data-toc-modified-id=\"The-0,0-conundrum-2.8.2\"><span class=\"toc-item-num\">2.8.2&nbsp;&nbsp;</span>The 0,0 conundrum</a></span></li></ul></li><li><span><a href=\"#Create-a-two-layer-map\" data-toc-modified-id=\"Create-a-two-layer-map-2.9\"><span class=\"toc-item-num\">2.9&nbsp;&nbsp;</span>Create a two layer map</a></span></li><li><span><a href=\"#Subplots-for-multi-layered-maps\" data-toc-modified-id=\"Subplots-for-multi-layered-maps-2.10\"><span class=\"toc-item-num\">2.10&nbsp;&nbsp;</span>Subplots for multi-layered maps</a></span></li><li><span><a href=\"#The-spatial-join\" data-toc-modified-id=\"The-spatial-join-2.11\"><span class=\"toc-item-num\">2.11&nbsp;&nbsp;</span>The spatial join</a></span></li><li><span><a href=\"#Join-the-value-counts-back-to-the-gdf\" data-toc-modified-id=\"Join-the-value-counts-back-to-the-gdf-2.12\"><span class=\"toc-item-num\">2.12&nbsp;&nbsp;</span>Join the value counts back to the gdf</a></span></li><li><span><a href=\"#Normalizing:-Arrests-per-1000-people\" data-toc-modified-id=\"Normalizing:-Arrests-per-1000-people-2.13\"><span class=\"toc-item-num\">2.13&nbsp;&nbsp;</span>Normalizing: Arrests per 1000 people</a></span></li><li><span><a href=\"#Choropleth-map-of-arrests\" data-toc-modified-id=\"Choropleth-map-of-arrests-2.14\"><span class=\"toc-item-num\">2.14&nbsp;&nbsp;</span>Choropleth map of arrests</a></span></li></ul></li><li><span><a href=\"#Global-Spatial-Autocorrelation\" data-toc-modified-id=\"Global-Spatial-Autocorrelation-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Global Spatial Autocorrelation</a></span><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><span><a href=\"#Spatial-Weights\" data-toc-modified-id=\"Spatial-Weights-3.0.1\"><span class=\"toc-item-num\">3.0.1&nbsp;&nbsp;</span>Spatial Weights</a></span></li><li><span><a href=\"#Spatial-lag\" data-toc-modified-id=\"Spatial-lag-3.0.2\"><span class=\"toc-item-num\">3.0.2&nbsp;&nbsp;</span>Spatial lag</a></span></li><li><span><a href=\"#The-donut-and-the-diamond\" data-toc-modified-id=\"The-donut-and-the-diamond-3.0.3\"><span class=\"toc-item-num\">3.0.3&nbsp;&nbsp;</span>The donut and the diamond</a></span></li></ul></li><li><span><a href=\"#Spatial-lag-map\" data-toc-modified-id=\"Spatial-lag-map-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Spatial lag map</a></span></li><li><span><a href=\"#Side-by-side-maps\" data-toc-modified-id=\"Side-by-side-maps-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Side-by-side maps</a></span></li><li><span><a href=\"#Interactive-spatial-lag-satellite-map\" data-toc-modified-id=\"Interactive-spatial-lag-satellite-map-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Interactive spatial lag satellite map</a></span></li><li><span><a href=\"#Moran's-Plot\" data-toc-modified-id=\"Moran's-Plot-3.4\"><span class=\"toc-item-num\">3.4&nbsp;&nbsp;</span>Moran's Plot</a></span></li></ul></li><li><span><a href=\"#Local-Spatial-Autocorrelation\" data-toc-modified-id=\"Local-Spatial-Autocorrelation-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Local Spatial Autocorrelation</a></span><ul class=\"toc-item\"><li><span><a href=\"#Moral-Local-Scatterplot\" data-toc-modified-id=\"Moral-Local-Scatterplot-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Moral Local Scatterplot</a></span></li><li><span><a href=\"#Spatial-Autocorrelation-Map\" data-toc-modified-id=\"Spatial-Autocorrelation-Map-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>Spatial Autocorrelation Map</a></span></li></ul></li><li><span><a href=\"#Resources\" data-toc-modified-id=\"Resources-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Resources</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "\n",
    "<h1>Take notice!</h1>\n",
    "<ul>\n",
    "    <li>This class will be recorded</li>\n",
    "</ul>\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-df0628a5-dfb4-4e7c-ac45-e80b52c54132",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Spatial Autocorrelation\n",
    "\n",
    "<img src=\"images/sa.png\" width=800>\n",
    "\n",
    "Visual interpretations are meaningful ways to determine spatial trends in our data. However, underlying factors—such as inconsistent geographies, scale, data gaps, overlapping data—have the potential to produce incorrect assumptions, as valuable information may be conveniently hidden from the visual output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-df0628a5-dfb4-4e7c-ac45-e80b52c54132",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<a href=\"https://en.wikipedia.org/wiki/Spatial_analysis#/media/File:Snow-cholera-map.jpg\" target=\"_blank\"><img src=\"https://upload.wikimedia.org/wikipedia/commons/c/c7/Snow-cholera-map.jpg\" width=400></a>\n",
    "\n",
    "One way to address this issue is to amend your visual output with geo-statistical validation. In this lab, we will look at one such approach: Spatial Autocorrelation. Spatial autocorrelation addresses the so-called \"First Law of Geography\":\n",
    "\n",
    "> “Everything is related to everything else. But near things are more related than distant things”. Waldo Tobler’s (1969) First Law of Geography\n",
    "\n",
    "In other words, things that happen somewhere are likely to also happen at nearby locations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-df0628a5-dfb4-4e7c-ac45-e80b52c54132",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Methodology\n",
    "In this lab, we take on the controversial topic of policing in Los Angeles, specifically looking at arrest records from the LAPD. Do arrest locations have a statistical significant tendency to cluster in certain communities? To answer this question, we not only look at the location of recorded arrests in the city, but compare these locations with other arrests. In short, we are seeking to see where spatial correlations occur based on the data. Our approach is:\n",
    "\n",
    "1. import census block group boundaries for Los Angeles\n",
    "1. import arrest data from the LA Open Data Portal\n",
    "1. spatially join the two datasets\n",
    "1. normalize the data to create arrests per 1000\n",
    "1. conduct [global spatial autocorrelation](https://geographicdata.science/book/notebooks/06_spatial_autocorrelation.html) using Moran's I\n",
    "1. conduct [local spatial autocorrelation](https://geographicdata.science/book/notebooks/07_local_autocorrelation.html) using Local Indicators of Spatial Association (LISAs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Libraries to use\n",
    "\n",
    "<img src=\"https://pandas.pydata.org/docs/_static/pandas.svg\" width=300>\n",
    "\n",
    "<img src=\"https://matplotlib.org/_static/logo2_compressed.svg\" width=300>\n",
    "\n",
    "<img src=\"https://geopandas.org/_static/geopandas_logo_web.svg\" width=300>\n",
    "\n",
    "- [Data vis with Pandas and Matplotlib](https://pandas.pydata.org/docs/user_guide/visualization.html)\n",
    "- [geopandas introduction](https://geopandas.org/getting_started/introduction.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Libraries to use\n",
    "\n",
    "<img src=\"images/esda.png\">\n",
    "\n",
    "<img src=\"images/splot.png\">\n",
    "\n",
    "- [ESDA](https://pysal.org/esda/)\n",
    "- [SPLOT](https://github.com/pysal/splot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00001-85a84e8b-13d6-4356-8656-408ce2a2630c",
    "deepnote_cell_type": "code",
    "execution_millis": 1181,
    "execution_start": 1605804400268,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "bc3fef15"
   },
   "outputs": [],
   "source": [
    "# to read and wrangle data\n",
    "import pandas as pd\n",
    "\n",
    "# to import data from LA Data portal\n",
    "from sodapy import Socrata\n",
    "\n",
    "# to create spatial data\n",
    "import geopandas as gpd\n",
    "\n",
    "# for basemaps\n",
    "import contextily as ctx\n",
    "\n",
    "# For spatial statistics\n",
    "import esda\n",
    "from esda.moran import Moran, Moran_Local\n",
    "\n",
    "import splot\n",
    "from splot.esda import moran_scatterplot, plot_moran, lisa_cluster,plot_moran_simulation\n",
    "\n",
    "import libpysal as lps\n",
    "\n",
    "# Graphics\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.express as px"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Data preparation\n",
    "\n",
    "<img src=\"images/method.png\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00002-8c57606b-b800-4154-99ce-86d0e45b3499",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Block Groups\n",
    "\n",
    "Our first task is to bring in a geography that will allow us to summarize the location of arrests. The smaller the geography, the better our spatial correlation results will be. Short of creating our own grid, the census block groups provides an easily accessible boundary layer at a human scale. Additionally, working with census geographies will allow for future analyses that may include census data.\n",
    "\n",
    "* Date source: \n",
    "   * [Census Reporter: ACS 2018 5 year: Table B01003: Total Population in Los Angeles: Census Block Groups](https://censusreporter.org/data/table/?table=B01003&geo_ids=16000US0644000,150|16000US0644000&primary_geo_id=16000US0644000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00003-27a9c38f-0d55-4f36-9268-f26b88e2ae29",
    "deepnote_cell_type": "code",
    "execution_millis": 3186,
    "execution_start": 1605804403202,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "b8c3febb"
   },
   "outputs": [],
   "source": [
    "# population block group data from census reporter\n",
    "gdf = gpd.read_file('data/acs2018_5yr_B01003_15000US060372711003.geojson')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Data cleanup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# show first 5 rows\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00004-e26db948-08d3-46f9-b0ce-1da8563c02d9",
    "deepnote_cell_type": "code",
    "execution_millis": 7,
    "execution_start": 1605804407207,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "494d51a4"
   },
   "outputs": [],
   "source": [
    "# show columns and data types\n",
    "gdf.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00005-65875226-9284-44aa-b62e-3a1e2212d84b",
    "deepnote_cell_type": "code",
    "execution_millis": 13,
    "execution_start": 1605804409430,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "7c13ef85"
   },
   "outputs": [],
   "source": [
    "# trim the data to the bare minimum columns\n",
    "gdf = gdf[['geoid','B01003001','geometry']]\n",
    "\n",
    "# rename the columns\n",
    "gdf.columns = ['FIPS','TotalPop','geometry']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00005-65875226-9284-44aa-b62e-3a1e2212d84b",
    "deepnote_cell_type": "code",
    "execution_millis": 13,
    "execution_start": 1605804409430,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "7c13ef85"
   },
   "outputs": [],
   "source": [
    "# last 5 rows\n",
    "gdf.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# delete last row which is for the entire city of LA\n",
    "gdf=gdf.drop(2515)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# fix FIPS code\n",
    "gdf['FIPS'] = gdf['FIPS'].str.replace('15000US','')\n",
    "gdf.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One more data cleanup: get rid of census blocks groups with less than 100 total population.\n",
    "\n",
    "- `sort_values()` [documentation](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.sort_values.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# sort by total pop\n",
    "gdf.sort_values(by='TotalPop').head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Subsetting the data\n",
    "\n",
    "- [Selecting a subset of a dataframe](https://pandas.pydata.org/docs/getting_started/intro_tutorials/03_subset_data.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# delete less than 100 population geographies\n",
    "gdf = gdf[gdf['TotalPop']>100]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Map the census block groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00006-c7db2dde-5f71-47f2-bfc2-0a0d94f1a03d",
    "deepnote_cell_type": "code",
    "execution_millis": 930,
    "execution_start": 1605804413958,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "e699b335"
   },
   "outputs": [],
   "source": [
    "# get the layers into a web mercator projection\n",
    "# reproject to web mercator\n",
    "gdf = gdf.to_crs(epsg=3857)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Using Matplotlib subplots\n",
    "\n",
    "- `plt.subplots()` [documentation](https://matplotlib.org/stable/gallery/subplots_axes_and_figures/subplots_demo.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00007-ffa16cca-a0b5-4e66-9a4c-73d90565821d",
    "deepnote_cell_type": "code",
    "execution_millis": 2725,
    "execution_start": 1605804417356,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "93ebb39"
   },
   "outputs": [],
   "source": [
    "# plot it!\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "gdf.plot(ax=ax,\n",
    "         color='black', \n",
    "         edgecolor='white',\n",
    "         lw=0.5,\n",
    "         alpha=0.4)\n",
    "\n",
    "# no axis\n",
    "ax.axis('off')\n",
    "\n",
    "# add a basemap\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00008-5ca674be-a4a5-4476-b7bd-3c86722b0246",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Get Arrest Data from LA Open Data Portal\n",
    "Next, we acquire the data using the socrata API. Use the socrata documentation to grab the code syntax for our arrests data.\n",
    "- [LA Data Portal](https://data.lacity.org)\n",
    "- [LA Data Portal arrest data](https://data.lacity.org/Public-Safety/Arrest-Data-from-2020-to-Present/amvf-fr72)\n",
    "- [Socrata endpoint documentation for this dataset](https://dev.socrata.com/foundry/data.lacity.org/amvf-fr72)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00009-684c0663-6a05-4f66-82c3-63b259473b81",
    "deepnote_cell_type": "code",
    "execution_millis": 4626,
    "execution_start": 1605804430393,
    "output_cleared": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "a077b4d5"
   },
   "outputs": [],
   "source": [
    "# connect to the data portal\n",
    "client = Socrata(\"data.lacity.org\", None)\n",
    "\n",
    "results = client.get(\"amvf-fr72\", \n",
    "                     limit=50000,\n",
    "                     where = \"arst_date between '2020-07-01T00:00:00' and '2021-01-31T00:00:00'\",\n",
    "                     order='arst_date desc')\n",
    "\n",
    "# Convert to pandas DataFrame\n",
    "arrests = pd.DataFrame.from_records(results)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00010-edd31267-0772-400b-9220-b4c2af28281d",
    "deepnote_cell_type": "code",
    "execution_millis": 2,
    "execution_start": 1605804438759,
    "output_cleared": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "d8a3f800"
   },
   "outputs": [],
   "source": [
    "arrests.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "arrests.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00011-87e77a3f-36e6-48db-bc94-ae8b2747e411",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convert data to a geodataframe\n",
    "\n",
    "Geopandas allows us to convert different types of data into a spatial format.\n",
    "- https://geopandas.org/gallery/create_geopandas_from_pandas.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00012-11bd538c-28c8-4c61-9994-e819c3feb7e8",
    "deepnote_cell_type": "code",
    "execution_millis": 409,
    "execution_start": 1605804442578,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "dcfb7ede"
   },
   "outputs": [],
   "source": [
    "# convert pandas dataframe to geodataframe\n",
    "arrests = gpd.GeoDataFrame(arrests, \n",
    "                         crs='EPSG:4326',\n",
    "                         geometry=gpd.points_from_xy(arrests.lon, arrests.lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00013-9c89a5ae-ae0f-4844-89e8-3b9105d30b6e",
    "deepnote_cell_type": "code",
    "execution_millis": 1327,
    "execution_start": 1605804596814,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "8eb5f93e"
   },
   "outputs": [],
   "source": [
    "# get the layers into a web mercator projection\n",
    "# reproject to web mercator\n",
    "arrests = arrests.to_crs(epsg=3857)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00014-5b46baf9-e870-4232-a425-28f0d65fb571",
    "deepnote_cell_type": "code",
    "execution_millis": 60,
    "execution_start": 1605804598151,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "ced73d7"
   },
   "outputs": [],
   "source": [
    "# convert lat/lon to floats\n",
    "arrests.lon = arrests.lon.astype('float')\n",
    "arrests.lat = arrests.lat.astype('float')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00015-37008c4f-0307-455b-b316-da6e65f084c5",
    "deepnote_cell_type": "code",
    "execution_millis": 3485,
    "execution_start": 1605804600104,
    "output_cleared": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "ec9920c6"
   },
   "outputs": [],
   "source": [
    "# map it!\n",
    "fig,ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "arrests.plot(ax=ax,\n",
    "             color='red',\n",
    "             markersize=1)\n",
    "\n",
    "# no axis\n",
    "ax.axis('off')\n",
    "\n",
    "# add a basemap\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00016-7294ac05-aa91-433c-aefd-bb0ee45e7ef5",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The 0,0 conundrum\n",
    "What is that red dot off the coast of Africa? Yes, that is the infamous spatial black hole, the [0,0] coordinate. There can be many reasons why those red dots get lost and find themselves there. No data may default to 0's, null values may be converted to 0's, or a human may have inadvertently entered 0's for unknown locations. Either which way, since these records do not have valid locations, they need to be deleted from our dataframe in order to proceed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00017-7d512d29-36e2-4a73-80ec-1a49b1845f1d",
    "deepnote_cell_type": "code",
    "execution_millis": 11,
    "execution_start": 1605804631750,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "dc9df796"
   },
   "outputs": [],
   "source": [
    "# subset the zero coordinate records\n",
    "arrests[arrests.lon==0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# drop the unmapped rows\n",
    "arrests = arrests[arrests.lon!=0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00015-37008c4f-0307-455b-b316-da6e65f084c5",
    "deepnote_cell_type": "code",
    "execution_millis": 3485,
    "execution_start": 1605804600104,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "ec9920c6"
   },
   "outputs": [],
   "source": [
    "# map it!\n",
    "fig,ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "arrests.plot(ax=ax,\n",
    "             color='red',\n",
    "             markersize=1)\n",
    "\n",
    "# no axis\n",
    "ax.axis('off')\n",
    "\n",
    "# add a basemap\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00020-c58096aa-ccb6-44c8-a710-06b464722ea2",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Create a two layer map\n",
    "\n",
    "- https://geopandas.org/mapping.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00021-54ac341c-54f2-4549-b9fd-97ba8c8f9628",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Since we want to zoom to the extent of the arrests layer (and not the block groups), get the bounding coordinates for our axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00022-530381e9-7db1-4821-bcee-a66b1b5983f4",
    "deepnote_cell_type": "code",
    "execution_millis": 489,
    "execution_start": 1605804643081,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "2933d808"
   },
   "outputs": [],
   "source": [
    "# get the bounding box coordinates for the arrest data\n",
    "minx, miny, maxx, maxy = arrests.geometry.total_bounds\n",
    "print(minx)\n",
    "print(maxx)\n",
    "print(miny)\n",
    "print(maxy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Subplots for multi-layered maps\n",
    "\n",
    "For our multi-layered maps, we are taking it one step further from our previous lab using matplotlib's `subplots`. `subplots` allows the creation of multiple plots on a gridded canvas. For our map, we only need a single subplot, but we are layering multiple datasets *on top of one another* on that subplot. To specify which subplot to put the layer on, you use the `ax` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00023-de1f44e5-dc3e-4a90-9423-884b716cdfa7",
    "deepnote_cell_type": "code",
    "execution_millis": 4318,
    "execution_start": 1605804911562,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "75d2d69a"
   },
   "outputs": [],
   "source": [
    "# set up the plot canvas with plt.subplots in one column, one row\n",
    "fig, ax = plt.subplots(1,1,figsize=(15, 15))\n",
    "\n",
    "# block groups\n",
    "gdf.plot(ax=ax, # this puts it in the ax plot\n",
    "        color='gray', \n",
    "        edgecolor='white',\n",
    "        alpha=0.5)\n",
    "\n",
    "# arrests\n",
    "arrests.plot(ax=ax, # this also puts it in the same ax plot\n",
    "            color='red',\n",
    "            markersize=1,\n",
    "            alpha=0.2)\n",
    "\n",
    "# use the bounding box coordinates to set the x and y limits\n",
    "ax.set_xlim(minx - 1000, maxx + 1000) # added/substracted value is to give some margin around total bounds\n",
    "ax.set_ylim(miny - 1000, maxy + 1000)\n",
    "\n",
    "# no axis\n",
    "ax.axis('off')\n",
    "\n",
    "# add a basemap\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00024-350b1c6f-5d68-4f37-841d-9b832e497dfe",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The spatial join\n",
    "\n",
    "* https://geopandas.org/mergingdata.html?highlight=spatial%20join\n",
    "\n",
    "In a Spatial Join, two geometry objects are merged based on their spatial relationship to one another."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "While the official documentation may seem confusing, consider the following as a rule of thumb. When you do a spatial join with `gpd.sjoin()`, you feed it three arguments: a left dataframe, a right dataframe, and a how statement.\n",
    "\n",
    "- **Left dataframe**: identify the layer you want *to* attach infomation that will come from the other layer\n",
    "- **Right dataframe**: identify the layer that you want to get information *from* to attach to the other layer\n",
    "\n",
    "Once you identify your left and right dataframes, use `how=\"left\"` to spatially join the two layers (think: \"I'm sending data from the right to the left\"). Note that this will result in a dataframe with the same records and datatype as the LEFT layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00025-83617000-2144-45cb-9145-3e070dc1da7e",
    "deepnote_cell_type": "code",
    "execution_millis": 1339,
    "execution_start": 1605806117634,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "8c53b4a0"
   },
   "outputs": [],
   "source": [
    "# Do the spatial join\n",
    "join = gpd.sjoin(arrests, gdf, how='left')\n",
    "join.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00026-a88db1eb-74f2-4c08-9db8-e964326dc862",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This creates a dataframe that has every arrest record with the corresponding FIPS code.\n",
    "\n",
    "Next, we create another dataframe that counts crime by their corresponding block group.\n",
    "\n",
    "- `value_counts()` [documentation](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.value_counts.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00027-ca714a55-7276-4cd5-8cfd-64ed87c51b47",
    "deepnote_cell_type": "code",
    "execution_millis": 0,
    "execution_start": 1605806144598,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "bbec6db9"
   },
   "outputs": [],
   "source": [
    "arrests_by_gdf = join.FIPS.value_counts().rename_axis('FIPS').reset_index(name='arrests_count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00028-0edb1f2f-bf0d-464e-91a5-8e00178c2fc7",
    "deepnote_cell_type": "code",
    "execution_millis": 1,
    "execution_start": 1605806249414,
    "output_cleared": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "1a26c77d"
   },
   "outputs": [],
   "source": [
    "arrests_by_gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- [Pandas bar plot documentation](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.plot.bar.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00029-c703c03b-ae27-4935-85ae-6f32896bcd05",
    "deepnote_cell_type": "code",
    "execution_millis": 256,
    "execution_start": 1605806267816,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "a09e3e30"
   },
   "outputs": [],
   "source": [
    "# make a bar chart of top 20 geographies\n",
    "arrests_by_gdf[:20].plot.bar(figsize=(20,4),\n",
    "                             x='FIPS',\n",
    "                             y='arrests_count')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00030-23e15688-b87f-4d89-8e04-78e98f59d8ce",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Join the value counts back to the gdf\n",
    "\n",
    "How many people know their census block number? The bar chart is nice, but it is not informative. Without spatial awareness, the data chart does little to convey knowledge. What we want is a choropleth map to accompany it. To do so, we merge the counts back to the block group gdf.\n",
    "\n",
    "- [merge documentation](https://geopandas.org/docs/user_guide/mergingdata.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00031-ef0c1b84-f54e-464b-9183-d7b6952c93ad",
    "deepnote_cell_type": "code",
    "execution_millis": 0,
    "execution_start": 1605806281114,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "6918564f"
   },
   "outputs": [],
   "source": [
    "# join the summary table back to the gdf\n",
    "gdf=gdf.merge(arrests_by_gdf,on='FIPS')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00032-de630962-35a2-45c8-a028-ad929bebbff7",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now the block group gdf has a new column for arrest counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00033-4405a8bc-92cd-44dc-9f79-e1756818d28f",
    "deepnote_cell_type": "code",
    "execution_millis": 10,
    "execution_start": 1605806285885,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "436bb116"
   },
   "outputs": [],
   "source": [
    "# our neighborhood table now has a count column\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Normalizing: Arrests per 1000 people\n",
    "Rather than proceeding with an absolute count of arrests, let's normalize it by number of people who live in the census block group.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gdf['arrests_per_1000'] = gdf['arrests_count']/gdf['TotalPop']*1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gdf.sort_values(by=\"arrests_per_1000\").tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And if you are curious, we can map a slice of the data. Here, we sort the values by descending arrest rate, and only show a slice of the data, the top 20 geographies using the handy `[:20]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00034-d60cf7ed-32db-4b0e-ac55-6dc3119c40c9",
    "deepnote_cell_type": "code",
    "execution_millis": 455,
    "execution_start": 1605806413842,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "c70d7036"
   },
   "outputs": [],
   "source": [
    "# map the top 20 geographies\n",
    "fig,ax = plt.subplots(figsize=(12,10))\n",
    "gdf.sort_values(by='arrests_per_1000',ascending=False)[:20].plot(ax=ax,\n",
    "                                                                 color='red',\n",
    "                                                                 edgecolor='white',\n",
    "                                                                 alpha=0.5,legend=True)\n",
    "\n",
    "\n",
    "# title\n",
    "ax.set_title('Top 20 locations of LAPD arrests per 1000 people (July 2020-January 2021)')\n",
    "\n",
    "# no axis\n",
    "ax.axis('off')\n",
    "\n",
    "# add a basemap\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00035-452962c1-1d0a-4b6e-9f25-bf2bac2eb75e",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Choropleth map of arrests\n",
    "\n",
    "Finally, we are ready to generate a choropleth map of arrests. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- [geopandas choropleths](https://geopandas.org/docs/user_guide/mapping.html#choropleth-maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00040-751a6cac-f335-4b67-b470-3e240904a0cf",
    "deepnote_cell_type": "code",
    "execution_millis": 1290,
    "execution_start": 1605806440098,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    },
    "source_hash": "cf9722a3"
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(15,15))\n",
    "\n",
    "gdf.plot(ax=ax,\n",
    "        column='arrests_per_1000', # this makes it a choropleth\n",
    "        legend=True,\n",
    "        alpha=0.8,\n",
    "        cmap='RdYlGn_r', # a diverging color scheme\n",
    "        scheme='quantiles') # how to break the data into bins\n",
    "\n",
    "ax.axis('off')\n",
    "ax.set_title('2020 July to January 2021 arrests per 1000 people',fontsize=22)\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00041-1c2e4521-6cae-4ca5-9731-1b9c505d6077",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "The map above is a good way to begin exploring spatial patterns in our data. What does this map tell you? Is it informative? Do you notice any significant clusters? What if you change the map? Notice the `scheme` argument is set to `naturalbreaks`. Experiment with other map classfications such as `equalinterval`, `quantiles`. How does each classification change the map? \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00042-b24787d4-0ab3-44d5-ac9e-91f78dd286e3",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Global Spatial Autocorrelation\n",
    "\n",
    "<img src=\"images/sa.png\" width=800>\n",
    "\n",
    "We have imported two datasets. Cleaned them up, spatialized them, and connected them spatially. We successfully mapped them to show the location of arrests per 1000 people by census block groups. The resulting map intuitively and visually tells us that there does appear to be spatial clusters of where arrests are more prevalent, but to what degree of certainty can we say so? Actually, very little, without statitistically backing up our determinations. Could this exact pattern be a matter of chance? Or is the pattern so distinct that there is no way it could have happened randomly?\n",
    "\n",
    "In order to answer this question, we conduct spatial autocorrelation, a process that determines to what degree an existing pattern is or is not random."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00042-b24787d4-0ab3-44d5-ac9e-91f78dd286e3",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Global Moran's I statistic is a way to *quantify* the degree to which similar geographies are clustered. To do so, we compare each geography based on a given value (in this case arrest counts) with that of its neighbors. The first step of this process is to define a \"spatial weight.\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00043-8ff210ab-9ca8-4ba3-a467-3b94e460e0d1",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spatial Weights\n",
    "\n",
    "<img src=\"https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-statistics-toolbox/GUID-A772FE33-6EE7-4B45-A654-A3C73235E9A2-web.gif\" width=400>\n",
    "\n",
    "Image source: [ESRI](https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-statistics-toolbox/generate-spatial-weights-matrix.htm)\n",
    "\n",
    "Spatial weights are how we determine the area’s neighborhood. There are different statistical methods that are used for determining spatial weights. Some use a contiguity model, assigning neighbors based on boundaries that touch each other. Others are based on distance, finding the closest neighbors based on the centroid of each geography. So which method should we use? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00043-8ff210ab-9ca8-4ba3-a467-3b94e460e0d1",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For this lab, we will use the KNN weight, where `k` is the number of \"nearest neighbors\" to count in the calculations. Let's proceed with `k=8` for our KNN spatial weights. \n",
    "\n",
    "We also **row standardize** the data: A technique for adjusting the weights in a spatial weights matrix. When weights are row standardized, each weight is divided by its row sum. The row sum is the sum of weights for a features neighbors.\n",
    "\n",
    "- https://geographicdata.science/book/notebooks/04_spatial_weights.html#distance-based-weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00044-c20d9e6b-953d-44f8-8f15-51052303c063",
    "deepnote_cell_type": "code",
    "execution_millis": 111,
    "execution_start": 1605806657728,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "3ee158c7"
   },
   "outputs": [],
   "source": [
    "# calculate spatial weight\n",
    "wq =  lps.weights.KNN.from_dataframe(gdf,k=8)\n",
    "\n",
    "# Row-standardization\n",
    "wq.transform = 'r'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00045-109ded14-a0b2-4c8c-a0a0-8454c76c230d",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spatial lag\n",
    "\n",
    "Now that we have our spatial weights assigned, we use it to calculate the spatial lag. While the mathematical operations are beyond the scope of this lab, you are welcome to check it out [here](https://geographicdata.science/book/notebooks/06_spatial_autocorrelation.html#spatial-lag). Simply put, the spatial lag is a calculated assignment to each geography in your data, which takes into account the data values from others in their  \"neighborhood\" as defined by the spatial weight. This operation can be done with a single line of code which is part of the pysal module, but the underlying calculations are not that difficult to understand: \n",
    "\n",
    "It takes the average of all the neighbors as defined by the spatial weight to come up with a single associated value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00046-020e64b9-d823-45ea-bdfd-6950215f410d",
    "deepnote_cell_type": "code",
    "execution_millis": 4,
    "execution_start": 1605806659021,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "dd2f79d2"
   },
   "outputs": [],
   "source": [
    "# create a new column for the spatial lag\n",
    "gdf['arrests_per_1000_lag'] = lps.weights.lag_spatial(wq, gdf['arrests_per_1000'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00047-80344eb5-7c68-4523-80fb-88e94e5ce690",
    "deepnote_cell_type": "code",
    "execution_millis": 0,
    "execution_start": 1605806663078,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "7aafb24b"
   },
   "outputs": [],
   "source": [
    "# sample gives us 10 random rows\n",
    "gdf.sample(10)[['TotalPop','arrests_count','arrests_per_1000','arrests_per_1000_lag']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "Take a moment to look at the values in the dataframe. What do they tell you?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The donut and the diamond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a column that calculates the difference betwen arrests and lag\n",
    "gdf['arrest_lag_diff'] = gdf['arrests_per_1000'] - gdf['arrests_per_1000_lag']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# output to get the head and tail\n",
    "gdf.sort_values(by='arrest_lag_diff')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00048-88965640-7413-4991-a51d-6abf7e886e78",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In order to better understand the significance of the spatial lag values, consider the following two geographies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# the FIPS with highest negative difference\n",
    "gdf_donut = gdf.sort_values(by='arrest_lag_diff').head(1)\n",
    "gdf_donut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# the FIPS with highest positive difference\n",
    "gdf_diamond = gdf.sort_values(by='arrest_lag_diff').tail(1)\n",
    "gdf_diamond"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To better illustrate our assumptions, let's display these locations using satellite imagery. To do so, we use Plotly Express and its Mapbox Satellite basemap (which provides high quality and up to date imagery). Note that the use of Mapbox satellite image tiles require a user token. \n",
    "\n",
    "- Grab your token from the [Mapbox account page](https://account.mapbox.com/) (after you have created an account)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# set the mapbox access token\n",
    "token = 'your_token'\n",
    "px.set_mapbox_access_token(token)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# subset donut, project to WGS84, and get its centroid\n",
    "gdf_donut = gdf_donut.to_crs('epsg:4326')\n",
    "\n",
    "# what's the centroid?\n",
    "minx, miny, maxx, maxy = gdf_donut.geometry.total_bounds\n",
    "center_lat_donut = (maxy-miny)/2+miny\n",
    "center_lon_donut = (maxx-minx)/2+minx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# subset diamond, project to WGS84, and get its centroid\n",
    "gdf_diamond = gdf_diamond.to_crs('epsg:4326')\n",
    "\n",
    "# what's the centroid?\n",
    "minx, miny, maxx, maxy = gdf_diamond.geometry.total_bounds\n",
    "center_lat_diamond = (maxy-miny)/2+miny\n",
    "center_lon_diamond = (maxx-minx)/2+minx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- [plotly choropleth maps](https://plotly.com/python/mapbox-county-choropleth/#)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "px.choropleth_mapbox(gdf_donut, \n",
    "                     geojson=gdf_donut.geometry, \n",
    "                     locations=gdf_donut.index, \n",
    "                     mapbox_style=\"satellite-streets\",\n",
    "                     zoom=14, \n",
    "                     center = {\"lat\": center_lat_donut, \"lon\": center_lon_donut},\n",
    "                     hover_data=['arrests_count','arrests_per_1000','arrests_per_1000_lag'],\n",
    "                     opacity=0.4,\n",
    "                     title='The Donut')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "px.choropleth_mapbox(gdf_diamond, \n",
    "                     geojson=gdf_diamond.geometry, \n",
    "                     locations=gdf_diamond.index, \n",
    "                     mapbox_style=\"satellite-streets\",\n",
    "                     zoom=12, \n",
    "                     center = {\"lat\": center_lat_diamond, \"lon\": center_lon_diamond},\n",
    "                     hover_data=['arrests_count','arrests_per_1000','arrests_per_1000_lag'],\n",
    "                     opacity=0.4,\n",
    "                     title='The Diamond')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What's the story here?\n",
    "\n",
    "Our donut is located in Venice Beach, in the famous enclave with homes adjacent to the grand canal known as part of the Venice Canal Historic District. It has a low arrests rate at just over 9 per 1000. It's spatial lag? 76! \n",
    "\n",
    "Our diamond, coincidentally, is the also famous Venice Beach Boardwalk. It is located in close proximity to our donut, with an over 400 arrests per 1000. It's spatial lag? Ninety, which is still high!\n",
    "\n",
    "What these numbers tell us is that the Venice Canal is like a donut. A prestine historical housing community which is surrounded by neighbors with high arrest rates. On the flip side, Venice Boardwalk is a diamond: a super high valued geography surrounded by lower neighbors. It indicates that the park itself is a magnet location for arrests, but its immediate neighbors have comparatively lower number of arrests."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Spatial lag map\n",
    "But we digress. Let's map the entire dataframe by the newly created spatial lag column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00051-c11fcc2b-df8f-420c-9365-3812f835d6ce",
    "deepnote_cell_type": "code",
    "execution_millis": 1107,
    "execution_start": 1605807561410,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "51b90aba"
   },
   "outputs": [],
   "source": [
    "# use subplots that make it easier to create multiple layered maps\n",
    "fig, ax = plt.subplots(1,1,figsize=(15, 15))\n",
    "\n",
    "# spatial lag choropleth\n",
    "gdf.plot(ax=ax,\n",
    "         figsize=(15,15),\n",
    "         column='arrests_per_1000_lag',\n",
    "         legend=True,\n",
    "         alpha=0.8,\n",
    "         cmap='RdYlGn_r',\n",
    "         scheme='quantiles')\n",
    "\n",
    "ax.axis('off')\n",
    "ax.set_title('July 2020-January 2021 Arrests per 1000 people',fontsize=22)\n",
    "\n",
    "ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00052-7673144c-5be9-47cd-9e0c-46d5b9fe324b",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Side-by-side maps\n",
    "We can now compare these two map outputs side by side. Notice that the syntax is a bit different from past labs where we have only worked with one figure at a time. This output produces 1 row, and 2 columns of figures in `subplots`.\n",
    "- [subplots documentation](https://matplotlib.org/3.3.0/gallery/subplots_axes_and_figures/subplots_demo.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00053-3d2f883f-e36b-42a8-b0e4-8d46ac893d67",
    "deepnote_cell_type": "code",
    "execution_millis": 1427,
    "execution_start": 1605807604042,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "26dd2021"
   },
   "outputs": [],
   "source": [
    "# create the 1x2 subplots\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 8))\n",
    "\n",
    "# two subplots produces ax[0] (left) and ax[1] (right)\n",
    "\n",
    "# regular count map on the left\n",
    "gdf.plot(ax=ax[0], # this assigns the map to the left subplot\n",
    "         column='arrests_per_1000', \n",
    "         cmap='RdYlGn_r', \n",
    "         scheme='quantiles',\n",
    "         k=5, \n",
    "         edgecolor='white', \n",
    "         linewidth=0, \n",
    "         alpha=0.75, \n",
    "           )\n",
    "\n",
    "\n",
    "ax[0].axis(\"off\")\n",
    "ax[0].set_title(\"Arrests per 1000\")\n",
    "\n",
    "# spatial lag map on the right\n",
    "gdf.plot(ax=ax[1], # this assigns the map to the right subplot\n",
    "         column='arrests_per_1000_lag', \n",
    "         cmap='RdYlGn_r', \n",
    "         scheme='quantiles',\n",
    "         k=5, \n",
    "         edgecolor='white', \n",
    "         linewidth=0, \n",
    "         alpha=0.75\n",
    "           )\n",
    "\n",
    "ax[1].axis(\"off\")\n",
    "ax[1].set_title(\"Arrests per 1000 - Spatial Lag\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Interactive spatial lag satellite map\n",
    "Building the equivalent map as an interactive javascript map is a bit more challenging. While there are several options to choose from, this lab will use plotly express's choropleth_mapbox feature.\n",
    "- https://plotly.com/python/mapbox-county-choropleth/#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# interactive version needs to be in WGS84\n",
    "gdf_web = gdf.to_crs('EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# what's the centroid?\n",
    "minx, miny, maxx, maxy = gdf_web.geometry.total_bounds\n",
    "center_lat_gdf_web = (maxy-miny)/2+miny\n",
    "center_lon_gdf_web = (maxx-minx)/2+minx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Unlike the matplotlib map, plotly's mapbox map only gives us a continuous scale option (there is no magical `scheme` option). To produce a similar quantile map, we need to calculate the values manually.\n",
    "\n",
    "As we want to produce a choropleth map based on our spatial lag column, let's get some simple stats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# some stats\n",
    "gdf_web.arrests_per_1000_lag.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# grab the median\n",
    "median = gdf_web.arrests_per_1000_lag.median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = px.choropleth_mapbox(gdf_web, \n",
    "                     geojson=gdf_web.geometry, # the geometry column\n",
    "                     locations=gdf_web.index, # the index\n",
    "                     mapbox_style=\"satellite-streets\",\n",
    "                     zoom=9, \n",
    "                     color='arrests_per_1000_lag',\n",
    "                     color_continuous_scale='RdYlGn_r',\n",
    "                     color_continuous_midpoint =median, # put the median as the midpoint\n",
    "                     range_color =(0,median*2),\n",
    "                     hover_data=['arrests_count','arrests_per_1000','arrests_per_1000_lag'],\n",
    "                     center = {\"lat\": center_lat_gdf_web, \"lon\": center_lon_gdf_web},\n",
    "                     opacity=0.8,\n",
    "                     width=1000,\n",
    "                     height=800,\n",
    "                     labels={\n",
    "                             'arrests_per_1000_lag':'Arrests per 1000 (Spatial Lag)',\n",
    "                             'arrests_per_1000':'Arrests per 1000',\n",
    "                     })\n",
    "fig.update_traces(marker_line_width=0.1, marker_line_color='white')\n",
    "fig.update_layout(margin={\"r\":0,\"t\":0,\"l\":0,\"b\":0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00054-cde72154-74e9-492d-9384-f438d860221d",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Moran's Plot\n",
    "\n",
    "We now have a spatial lag map: a map that displays geographies weighted against the values of its neighbors. The clusters are much clearer and cleaner than the original arrest count map. Downtown, Venice, South LA, Van Nuys... But we still have not *quantified* the degree of the spatial correlations. To begin this process, we test for global autocorrelation for a continuous attribute (arrest counts)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "y = gdf.arrests_per_1000\n",
    "moran = Moran(y, wq)\n",
    "moran.I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The moran's I value is nothing more than the calculated slope of the scatterplot of our \"arrests per 1000\" and \"arrests per 1000 spatial lag\" columns. It does indicate whether or not you have a positive or negative autocorrelation. Values will range from positive one, to negative one. \n",
    "\n",
    "- **Positive** spatial autocorrelation: high values are close to high values, and/or low values are close to low values\n",
    "- **Negative** spatial autocorrelation (less common): similar values are far from each other; high values are next to low values, low values are next to high values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can output a scatterplot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = moran_scatterplot(moran, aspect_equal=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So what is the significance of our Moran value of 0.3? In other words, **how likely is our observed pattern on the map generated by an entirely random process?** To find out, we compare our value with a simulation of 999 permutations that randomly shuffles the arrest data throughout the given geographies. The output is a sampling distribution of Moran’s I values under the (null) hypothesis that attribute values are randomly distributed across the study area. We then compare our observed Moran’s I value to this \"Reference Distribution.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00055-489abd0d-4dfc-4a43-ab07-631157aaea0d",
    "deepnote_cell_type": "code",
    "execution_millis": 47,
    "execution_start": 1605807625694,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "slide"
    },
    "source_hash": "707e03ac"
   },
   "outputs": [],
   "source": [
    "plot_moran_simulation(moran,aspect_equal=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00056-8d12869d-48fa-468e-9cb5-1babfcbf3bf9",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can compute the P-value:\n",
    "\n",
    "<img src=\"https://www.simplypsychology.org/p-value.png\" width=600>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00057-90a65e59-c521-42fa-bc23-4141b9369225",
    "deepnote_cell_type": "code",
    "output_cleared": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "moran.p_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00058-368e0922-3062-4b2a-9155-27c2332200fd",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The value is calculated as an empirical P-value that represents the proportion of realisations in the simulation under spatial randomness that are more extreme than the observed value. A small enough p-value associated with the Moran’s I of a map allows to reject the hypothesis that the map is random. In other words, we can conclude that the map displays more spatial pattern than we would expect if the values had been randomly allocated to a locations.\n",
    "\n",
    "That is a very low value, particularly considering it is actually the minimum value we could have obtained given the simulation behind it used 999 permutations (default in PySAL) and, by standard terms, it would be deemed statistically significant. We can ellaborate a bit further on the intuition behind the value of p_sim. If we generated a large number of maps with the same values but randomly allocated over space, and calculated the Moran’s I statistic for each of those maps, only 0.1% of them would display a larger (absolute) value than the one we obtain from the observed data, and the other 99.9% of the random maps would receive a smaller (absolute) value of Moran’s I. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Local Spatial Autocorrelation\n",
    "So far, we have only determined that there is a positive spatial autocorrelation between the price of properties in neighborhoods and their locations. But we have not detected where clusters are. Local Indicators of Spatial Association (LISA) is used to do that. LISA classifies areas into four groups: high values near to high values (HH), Low values with nearby low values (LL), Low values with high values in its neighborhood, and vice-versa.\n",
    "\n",
    "- HH: high arrest rate geographies near other high arrest rate neighbors\n",
    "- LL: low arrest rate geographies near other low arrest rate neighbors\n",
    "- LH (donuts): low arrest rate geographies surrounded by high arrest neighbors\n",
    "- HL (diamonds): high arrest rate geographies surrounded by low arrest neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Moral Local Scatterplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00060-4a2f2b79-53f9-44b4-ae86-31cc8fdce2bd",
    "deepnote_cell_type": "code",
    "execution_millis": 3398,
    "execution_start": 1605807811791,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "bf84c31e"
   },
   "outputs": [],
   "source": [
    "# calculate local moran values\n",
    "lisa = esda.moran.Moran_Local(y, wq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00063-d9fc9212-bf39-4525-90db-ff4cbcfd7783",
    "deepnote_cell_type": "code",
    "execution_millis": 214,
    "execution_start": 1605807831354,
    "output_cleared": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "2b643"
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig,ax = plt.subplots(figsize=(10,15))\n",
    "\n",
    "moran_scatterplot(lisa, ax=ax, p=0.05)\n",
    "ax.set_xlabel(\"Arrests\")\n",
    "ax.set_ylabel('Spatial Lag of Arrests')\n",
    "\n",
    "# add some labels\n",
    "plt.text(1.95, 0.5, \"HH\", fontsize=25)\n",
    "plt.text(1.95, -1, \"HL\", fontsize=25)\n",
    "plt.text(-2, 1, \"LH\", fontsize=25)\n",
    "plt.text(-1, -1, \"LL\", fontsize=25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the scatterplot above, the colored dots represents the rows (census block groups) that have a P-value less that 0.05 in each quadrant. In other words, these are the statisticaly significantly, spatially autocorrelated geographies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Spatial Autocorrelation Map\n",
    "Finally, you can visually these statistically significant clusters using the `lisa_cluster` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00064-5038673d-0806-4e6c-a49e-9938454eedd4",
    "deepnote_cell_type": "code",
    "execution_millis": 819,
    "execution_start": 1605807832811,
    "output_cleared": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "4f4bbffa"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(14,12))\n",
    "lisa_cluster(lisa, gdf, p=0.05, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And create a map comparing different p-values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "00053-3d2f883f-e36b-42a8-b0e4-8d46ac893d67",
    "deepnote_cell_type": "code",
    "execution_millis": 1427,
    "execution_start": 1605807604042,
    "output_cleared": false,
    "slideshow": {
     "slide_type": "fragment"
    },
    "source_hash": "26dd2021"
   },
   "outputs": [],
   "source": [
    "# create the 1x2 subplots\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 8))\n",
    "\n",
    "# regular count map on the left\n",
    "lisa_cluster(lisa, gdf, p=0.05, ax=ax[0])\n",
    "\n",
    "ax[0].axis(\"off\")\n",
    "ax[0].set_title(\"P-value: 0.05\")\n",
    "\n",
    "# spatial lag map on the right\n",
    "lisa_cluster(lisa, gdf, p=0.01, ax=ax[1])\n",
    "ax[1].axis(\"off\")\n",
    "ax[1].set_title(\"P-value: 0.01\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00073-51e53eb0-0154-43e8-835d-ce181f9c7186",
    "deepnote_cell_type": "markdown",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Resources\n",
    "\n",
    "- https://geographicdata.science/book/notebooks/06_spatial_autocorrelation.html\n",
    "- https://pysal.org/esda/notebooks/spatialautocorrelation.html\n",
    "- https://towardsdatascience.com/what-is-exploratory-spatial-data-analysis-esda-335da79026ee"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "deepnote_execution_queue": [],
  "deepnote_notebook_id": "e3ee145c-1ca6-4b9e-9637-f285c2264264",
  "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.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "349.047px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}