modules/06-spatial-data/lecture.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatial Data\n",
    "\n",
    "Overview of today's topics:\n",
    "\n",
    "  - Working with shapefiles, GeoPackages, CSV files, and rasters\n",
    "  - Projection\n",
    "  - Geometric operations\n",
    "  - Spatial joins\n",
    "  - Web mapping\n",
    "  - Spatial indexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ast\n",
    "\n",
    "import contextily as cx\n",
    "import folium\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import rasterio\n",
    "import rasterio.features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Loading a shapefile or GeoPackage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tell geopandas to read shapefile with its read_file() function, passing in shapefile folder\n",
    "# this produces a GeoDataFrame\n",
    "gdf_tracts = gpd.read_file(\"../../data/tl_2020_06_tract/\")\n",
    "gdf_tracts.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just like regular pandas, see the first 5 rows of the GeoDataFrame\n",
    "# this is a shapefile of polygon geometries, that is, tract boundaries\n",
    "gdf_tracts.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rudimentary mapping is as easy as calling the GeoDataFrame's plot method\n",
    "ax = gdf_tracts.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what is the CRS?\n",
    "# this derives from the shapefile's .prj file\n",
    "# always make sure the shapefile you load has prj info so you get a CRS attribute!\n",
    "gdf_tracts.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading a GeoPackage works the same way\n",
    "gdf_stations = gpd.read_file(\"../../data/rail_stations.gpkg\")\n",
    "gdf_stations.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_stations.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Loading a CSV file\n",
    "\n",
    "Often, you won't have a shapefile or GeoPackage (which is explicitly spatial), but rather a CSV file which is implicitly spatial (contains lat-lng columns). If you're loading a CSV file (or other non-explicitly spatial file type) of lat-lng data:\n",
    "\n",
    "  1. first load the CSV file as a DataFrame the usual way with pandas\n",
    "  2. then create a new geopandas GeoDataFrame from your DataFrame\n",
    "  3. manually create a geometry column\n",
    "  4. set the CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# load business location data as a regular pandas dataframe\n",
    "df = pd.read_csv(\"../../data/Listing_of_Active_Businesses.csv\")\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clean up the data (same code from the data cleaning lecture)\n",
    "df.columns = df.columns.str.lower().str.replace(\" \", \"_\").str.strip(\"_#\")\n",
    "df = df.set_index(\"location_account\").sort_index()\n",
    "df[\"location_start_date\"] = pd.to_datetime(df[\"location_start_date\"])\n",
    "slicer = pd.IndexSlice[:, \"business_name\":\"mailing_city\"]\n",
    "df.loc[slicer] = df.loc[slicer].apply(lambda col: col.str.title(), axis=\"rows\")\n",
    "mask = pd.notnull(df[\"location\"])\n",
    "latlng = df.loc[mask, \"location\"].map(ast.literal_eval)\n",
    "df.loc[mask, [\"lat\", \"lng\"]] = pd.DataFrame(\n",
    "    latlng.to_list(), index=latlng.index, columns=[\"lat\", \"lng\"]\n",
    ")\n",
    "df = df.drop(columns=[\"location\"]).dropna(subset=[\"lat\", \"lng\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# examine first five rows\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a geometry array to contain shapely geometry for geopandas to use\n",
    "# notice the shapely points are lng, lat so that they are equivalent to x, y\n",
    "geometry = gpd.points_from_xy(x=df[\"lng\"], y=df[\"lat\"])\n",
    "geometry[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a geopandas geodataframe from the pandas dataframe\n",
    "# also notice that we set the CRS explicitly\n",
    "gdf_business = gpd.GeoDataFrame(df, geometry=geometry, crs=\"epsg:4326\")\n",
    "gdf_business.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Always define the CRS** if you are manually creating a GeoDataFrame! Earlier, when we loaded the shapefile, geopandas loaded the CRS from the shapefile itself. But our CSV file is not explicitly spatial and it contains no CRS data, so we have to tell it what it is. In our case, the CRS is EPSG:4326, which is WGS84 lat-lng data, such as for GPS. Your data source should always tell you what CRS their coordinates are in. If they don't, ask! Don't just guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_business.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what's the CRS\n",
    "gdf_business.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Loading a raster\n",
    "\n",
    "So far we've worked with vector data. We can also work with raster data. Raster datasets are grids of pixels, where each pixel has a value (or multiple values if multiple bands of data), while vector datasets contain geometry objects with attributes, where each geometry is represented by mathematical coordinates for points, lines, polygons, etc. \"Raster is faster but vector is corrector.\"\n",
    "\n",
    "Common raster data include:\n",
    "  - tree cover\n",
    "  - urbanization footprints\n",
    "  - land use\n",
    "  - elevation\n",
    "\n",
    "In this example we load the SRTM 30m elevation raster, downloaded from https://dwtkns.com/srtm30m/, and cropped (to make a small dataset that can fit in laptop memory) via [raster-crop-bbox.ipynb](raster-crop-bbox.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load raster file and view its band count, pixel width/height, null value, and geographic bounds\n",
    "raster = rasterio.open(\"../../data/la-elevation.tif\")\n",
    "print(raster.count, raster.width, raster.height)\n",
    "print(raster.nodata)\n",
    "print(raster.bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the raster data\n",
    "df = pd.DataFrame(raster.read(1))\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram of elevations (meters above sea level) around downtown LA\n",
    "ax = df[df != raster.nodata].stack().hist(bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get shapes representing groups of adjacent pixels with same values\n",
    "# affine transformation maps pixel row/col -> spatial x/y\n",
    "shapes = rasterio.features.shapes(source=raster.read(1), transform=raster.transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert raster to GeoJSON-like vector features and create a gdf from them\n",
    "# pro-tip: use generator comprehension for memory efficiency\n",
    "features = ({\"geometry\": polygon, \"properties\": {\"elevation\": value}} for polygon, value in shapes)\n",
    "gdf_srtm = gpd.GeoDataFrame.from_features(features, crs=raster.crs)\n",
    "\n",
    "# drop any null rows\n",
    "gdf_srtm = gdf_srtm[gdf_srtm[\"elevation\"] != raster.nodata]\n",
    "gdf_srtm.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the gdf\n",
    "gdf_srtm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check its crs\n",
    "gdf_srtm.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the elevation pixels and identify pershing square\n",
    "fig, ax = plt.subplots(facecolor=\"#111111\")\n",
    "ax = gdf_srtm.plot(ax=ax, column=\"elevation\", cmap=\"inferno\")\n",
    "_ = ax.axis(\"off\")\n",
    "_ = ax.scatter(y=34.048097, x=-118.253233, c=\"w\", marker=\"x\", s=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# change the colors and also show the location of city hall on the map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Projection\n",
    "\n",
    "Your datasets need to be in the same CRS if you want to work with them together. If they're not, then project one or more of them so they're in the same CRS.\n",
    "\n",
    "Take note of the important difference here between *setting* a CRS (i.e., identifying a dataset's current CRS) and *projecting* to a CRS (i.e., mathematically transforming your coordinates from their current CRS to a different one). Projection lets you transform, for example, from lat-lng coordinates on the surface of the round Earth to a flat two-dimensional plane for mapping and analysis in intuitive units like meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if all our datasets have the same CRS\n",
    "gdf_tracts.crs == gdf_stations.crs == gdf_business.crs == gdf_srtm.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# project them all to UTM zone 11N (see http://epsg.io/32611)\n",
    "utm_crs = \"epsg:32611\"\n",
    "gdf_tracts = gdf_tracts.to_crs(utm_crs)\n",
    "gdf_stations = gdf_stations.to_crs(utm_crs)\n",
    "gdf_business = gdf_business.to_crs(utm_crs)\n",
    "gdf_srtm = gdf_srtm.to_crs(utm_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if all our datasets have the same CRS\n",
    "gdf_tracts.crs == gdf_stations.crs == gdf_business.crs == gdf_srtm.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Be careful**: heed the difference between the `gdf.crs` attribute and the `gdf.to_crs()` method. The former is the geodataframe's current CRS, whereas the latter projects the geodataframe to a new CRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# pick a different CRS and re-project the data to it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Geometric operations\n",
    "\n",
    "GIS and spatial analysis use common \"computational geometry\" operations like intersects, within, and dissolve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# takes a few seconds...\n",
    "# dissolve lets you aggregate (merge geometries together) by shared attribute values\n",
    "# this is the spatial equivalent of pandas's groupby function\n",
    "gdf_counties = gdf_tracts.dissolve(by=\"COUNTYFP\", aggfunc=\"sum\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now that we've dissolved tracts -> counties and summed their attributes,\n",
    "# plot the counties by land area\n",
    "fig, ax = plt.subplots(facecolor=\"#111111\")\n",
    "ax = gdf_counties.plot(ax=ax, column=\"ALAND\", cmap=\"Blues_r\")\n",
    "_ = ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just like in regular pandas, we can filter and subset the GeoDataFrame\n",
    "# retain only tracts in LA county (FIPS code 037)\n",
    "mask = gdf_tracts[\"COUNTYFP\"] == \"037\"\n",
    "gdf_tracts_la = gdf_tracts[mask]\n",
    "ax = gdf_tracts_la.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discard the channel islands' tracts to retain only the mainland\n",
    "# how? sort by centroids' y-coord and discard the two southern-most\n",
    "labels = gdf_tracts_la.centroid.y.sort_values().iloc[2:].index\n",
    "gdf_tracts_la = gdf_tracts_la.loc[labels]\n",
    "ax = gdf_tracts_la.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unary union merges all geometries in gdf into one\n",
    "la_geom = gdf_tracts_la.unary_union\n",
    "la_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convex hull generates the minimal convex polygon around feature(s)\n",
    "la_geom.convex_hull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# envelope generates the minimal rectangular polygon around feature(s)\n",
    "la_geom.envelope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a bounding box around our elevation data\n",
    "elev_bounds = gdf_srtm.unary_union.envelope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many spatial operations, such as intersects/within, scale in time complexity as a function of 1) the number of objects, and 2) the number of vertices in the reference polygon. Using a simplified reference polygon, such as a bounding box, can drastically speed up your operation at the cost of imprecision. In this case, our raster is already approximately square, and we don't need to do precise matching, so let's use the bounding box for intersects/within to filter our tracts and businesses to those that lie within the area covered by our elevation data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all the tracts that intersect those bounds\n",
    "# intersects tells you if each geometry in one dataset intersects with some other (single) geometry\n",
    "mask = gdf_tracts_la.intersects(elev_bounds)\n",
    "gdf_tracts_dtla = gdf_tracts_la[mask]\n",
    "gdf_tracts_dtla.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all business points within those bounds\n",
    "# within tells you if each geometry in one dataset is within some other (single) geometry\n",
    "mask = gdf_business.within(elev_bounds)\n",
    "gdf_business_dtla = gdf_business[mask]\n",
    "gdf_business_dtla.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# euclidean buffers let you analyze the area around features (use projected CRS!)\n",
    "# buffer the rail stations by a half km (5-10 minute walk)\n",
    "gdf_stations[\"geometry\"] = gdf_stations.buffer(500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 8), facecolor=\"#111111\")\n",
    "ax = gdf_tracts_dtla.plot(ax=ax, color=\"k\")\n",
    "ax = gdf_stations.plot(ax=ax, color=\"w\", alpha=0.3)\n",
    "ax = gdf_business_dtla.plot(\n",
    "    ax=ax, color=\"#ffff66\", marker=\".\", linewidth=0, markersize=20, alpha=0.05\n",
    ")\n",
    "_ = ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# you can do set operations like union, intersection, and difference\n",
    "# get all the portions of tracts >0.5km from a rail station\n",
    "gdf_diff = gpd.overlay(gdf_tracts_dtla, gdf_stations, how=\"difference\")\n",
    "ax = gdf_diff.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Spatial join\n",
    "\n",
    "Joins two geodataframes based on some shared spatial location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remember (again): always double-check CRS before any spatial operations\n",
    "gdf_tracts_dtla.crs == gdf_stations.crs == gdf_business.crs == gdf_srtm.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6a. How hilly is it around the stations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# join stations to elevation data\n",
    "gdf = gpd.sjoin(gdf_srtm, gdf_stations, how=\"inner\", predicate=\"intersects\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# counts vary because these aren't elevation pixels, but regions of same value\n",
    "gdf_elev_desc = gdf.groupby(\"name\")[\"elevation\"].describe().astype(int)\n",
    "gdf_elev_desc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_stations_elev = gdf_stations.merge(gdf_elev_desc, left_on=\"name\", right_index=True)\n",
    "gdf_stations_elev.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# which stations have the greatest elevation variation around them?\n",
    "ax = gdf_stations_elev.plot(column=\"std\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# which station buffer covers the largest elevation range?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6b. Which stations have the most businesses in their catchment areas?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# join stations to businesses data\n",
    "gdf = gpd.sjoin(gdf_business, gdf_stations, how=\"inner\", predicate=\"intersects\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# counts vary because these aren't elevation pixels, but regions of same value\n",
    "gdf_business_desc = gdf.groupby(\"name\").size().sort_values(ascending=False)\n",
    "gdf_business_desc.name = \"count\"\n",
    "gdf_business_desc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beware artificial peripheries! Some station buffers extend beyond the spatially-cropped business locations. How would you fix this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# change earlier parts of the notebook to make sure\n",
    "# our station buffers capture all the businesses within them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_stations_business = gdf_stations.merge(gdf_business_desc, left_on=\"name\", right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# which stations have the most businesses around them?\n",
    "ax = gdf_stations_business.plot(column=\"count\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works ok as a quick and dirty way to visually inspect our results. But it only works because we're analyzing/visualizing counts across study sites (i.e., station buffers) that *are all the same size* as each other. If the study site sizes varied (such as tracts or counties), counts might be correlated with area! Then you're just visualizing which study sites are the largest. In such cases, make sure you normalize. For example, use densities instead of counts.\n",
    "\n",
    "### 6c. Which tracts have the most businesses?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# join tracts to business data\n",
    "gdf = gpd.sjoin(gdf_business, gdf_tracts_dtla, how=\"inner\", predicate=\"intersects\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# count businesses per tract\n",
    "counts = gdf.groupby(\"GEOID\").size()\n",
    "counts.name = \"count\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge in the counts then calculate density (businesses per km^2)\n",
    "gdf_tracts_dtla_business = gdf_tracts_dtla.merge(counts, left_on=\"GEOID\", right_index=True)\n",
    "gdf_tracts_dtla_business[\"density\"] = (\n",
    "    gdf_tracts_dtla_business[\"count\"] / gdf_tracts_dtla_business[\"ALAND\"] * 1e6\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot tracts as choropleth plus station buffers\n",
    "fig, ax = plt.subplots(figsize=(8, 8), facecolor=\"#111111\")\n",
    "ax = gdf_tracts_dtla_business.plot(ax=ax, column=\"density\", cmap=\"viridis\")\n",
    "ax = gdf_stations.plot(ax=ax, alpha=0.2, linewidth=3, edgecolor=\"w\", color=\"none\")\n",
    "_ = ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this time, let's add a basemap for context\n",
    "fig, ax = plt.subplots(figsize=(8, 8), facecolor=\"#111111\")\n",
    "ax = gdf_tracts_dtla_business.plot(\n",
    "    ax=ax, column=\"density\", cmap=\"viridis\", alpha=0.7, linewidth=0.3, edgecolor=\"k\"\n",
    ")\n",
    "ax = gdf_stations.plot(ax=ax, alpha=0.3, linewidth=3, edgecolor=\"w\", color=\"none\")\n",
    "_ = ax.axis(\"off\")\n",
    "\n",
    "# add the basemap with contextily, choosing a tile provider\n",
    "# or try cx.providers.Stamen.TonerBackground, etc\n",
    "cx.add_basemap(ax, crs=gdf_stations.crs.to_string(), source=cx.providers.CartoDB.DarkMatter)\n",
    "\n",
    "ax.figure.savefig(\"map.png\", dpi=600, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# change the tile provider, the tract colors, the alphas, etc to find a plot your like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about an interactive web map instead?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optionally bin the data into quintiles\n",
    "bins = list(gdf_tracts_dtla_business[\"density\"].quantile([0, 0.2, 0.4, 0.6, 0.8, 1]))\n",
    "\n",
    "# create leaflet choropleth web map\n",
    "m = folium.Map(location=(34.047223, -118.253555), zoom_start=15, tiles=\"cartodbdark_matter\")\n",
    "c = folium.Choropleth(\n",
    "    geo_data=gdf_tracts_dtla_business,\n",
    "    data=gdf_tracts_dtla_business,\n",
    "    # bins=bins,\n",
    "    columns=[\"GEOID\", \"density\"],\n",
    "    key_on=\"feature.properties.GEOID\",\n",
    "    highlight=True,\n",
    "    fill_color=\"YlOrRd_r\",\n",
    "    legend_name=\"Businesses per square km\",\n",
    ").add_to(m)\n",
    "\n",
    "# add mouseover tooltip to the tracts\n",
    "c.geojson.add_child(folium.features.GeoJsonTooltip([\"GEOID\", \"density\"]))\n",
    "\n",
    "# save web map to disk\n",
    "m.save(\"webmap.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now it's your turn\n",
    "# try binning the data in different ways. how would you do it?\n",
    "# try changing the colors, basemap, and what variable you're visualizing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Spatial Indexing\n",
    "\n",
    "When you need to find which page a topic appears on in a book, do you search through every word, page by page, until you find it? When you need to find which polygon a point lies in, do you search through every polygon, one at a time, until you find it? Sometimes. But you can avoid that slow brute-force search if you use an index.\n",
    "\n",
    "A spatial index such as an [r-tree](https://en.wikipedia.org/wiki/R-tree) can drastically speed up spatial operations like intersects and joins. In computer science, a tree data structure represents parent and children objects like the branches of a tree. For example, a k-d tree lets you partition space for fast nearest-neighbor search. But an r-tree is particularly useful for finding what geometries intersect with some other geometry, such as point-in-polygon queries.\n",
    "\n",
    "An r-tree represents individual objects and their bounding boxes (\"r\" is for \"rectangle\") as the lowest level of the spatial index. It then aggregates nearby objects and represents them with their aggregate bounding box in the next higher level of the index. At yet higher levels, the r-tree aggregates bounding boxes and represents them by their bounding box, iteratively, until everything is nested into one top-level bounding box.\n",
    "\n",
    "To search, the r-tree takes a query box and, starting at the top level, sees which (if any) bounding boxes intersect it. It then expands each intersecting bounding box and sees which of the child bounding boxes inside it intersect the query box. This proceeds recursively until all intersecting boxes are searched down to the lowest level, and returns the matching objects from the lowest level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# geopandas uses r-tree spatial indexes\n",
    "# if a spatial index doesn't already exist,\n",
    "# it will be created the first time the sindex attribute is accessed\n",
    "sindex = gdf_business.sindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# count all the businesses within the station buffers\n",
    "polygon = gdf_stations.unary_union\n",
    "gdf_business.within(polygon).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can break this out into a two-step process. First find approximate matches with spatial index, then precise matches from those approximate ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(gdf_business)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# get the positions of possible matches, to use with iloc\n",
    "positions = sindex.intersection(polygon.bounds)\n",
    "possible_matches = gdf_business.iloc[positions]\n",
    "len(possible_matches)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was fast! And we're nearly there. We intersected the spatial index with the bounds of our polygon. This returns a set of possible matches. That is, there are no false negatives but there may be some false positives if an r-tree rectangle within the bounds contains some points outside the tracts' true borders.\n",
    "\n",
    "Unfortunately, the heavy lifting remains in filtering down the possible matches within the bounds to figure out which are within the polygon itself. To identify the precise matches (those points exactly within our polygon), we intersect the possible matches with the polygon itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "mask = possible_matches.intersects(polygon)\n",
    "precise_matches = possible_matches[mask]\n",
    "len(precise_matches)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the r-tree lets us filter out ~90% of the points (from the rest of the county) nearly instantly, but then the final precise point-in-polygon search (i.e., of the remaining points within the station buffers' bounding box, which are within the station buffers themselves?) consumes nearly all the runtime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what false positives appeared among the possible matches?\n",
    "labels = possible_matches.index.difference(precise_matches.index)\n",
    "false_positives = possible_matches.loc[labels]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the precise matches vs the false positives\n",
    "ax = gdf_stations.plot(color=\"gray\")\n",
    "ax = false_positives.plot(ax=ax, c=\"r\", markersize=0.1)\n",
    "ax = precise_matches.plot(ax=ax, c=\"b\", markersize=0.1)\n",
    "_ = ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (ppde642)",
   "language": "python",
   "name": "ppde642"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}