{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4d5266c4-0da0-43bd-8e8f-01dd39a0488e",
   "metadata": {},
   "source": [
    "# Introduction to Geospatial Python\n",
    "\n",
    "---\n",
    "\n",
    "**A Tufts University Data Lab Tutorial**  \n",
    "Written by [Uku-Kaspar Uustalu](https://directory.tufts.edu/user/view/90E8E773F8EC92B23679584546E5E321/)\n",
    "\n",
    "Contact: <uku-kaspar.uustalu@tufts.edu>\n",
    "\n",
    "Last updated: `2023-02-27`\n",
    "\n",
    "***WORK IN PROGRESS***"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c147f789-d25a-4b66-b49b-a6448d327b16",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Reading Data with Geographic Coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45c85a64-fb05-43fd-835a-dea12d1b5bde",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cc08761-a252-49dc-9be2-0290156bf4fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops = pd.read_csv('data/mbta-transit-stations.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f47e6c9d-5ece-4bc7-a253-0ef41f5b3844",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28fd812d-d6a9-403b-84f7-4bc031424b30",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops = gpd.GeoDataFrame(data=stops,\n",
    "                         geometry=gpd.points_from_xy(x=stops.stop_lon,\n",
    "                                                     y=stops.stop_lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64d2871d-7d52-4683-9144-97dd29be86c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e5c4edb-49ca-45aa-bbc9-44d0611d35b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(stops)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85e5ac93-9b5a-4f97-9dbf-5cc21d293b10",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(stops.geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6981b17-e451-4df8-97dc-06fb218e9f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(stops.geometry[0])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5d2173a3-115f-4902-8a72-33e56f039b5e",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Coordinate Reference Systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b64bbb2d-7cff-4a1f-9e3d-10769c153fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "736ff1b0-7661-434b-ab80-ba28a97ee64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops = stops.set_crs(epsg=4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b73c518b-ace8-49f1-9f64-5dfdf555d95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56354346-aa9c-400c-967d-b937dbdfcafc",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.crs.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38db0691-28a8-46e7-ad33-493423675ded",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.to_crs(epsg=26986, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c6ae7f4-0ca1-4d50-adb6-5d46a5c84f48",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e782e33-86ff-4a1c-baaf-bd393ebb6e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a710b01-2991-404e-8a8c-ddc1a9a2923b",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.crs.name"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "499f0265-d817-4483-b703-d8ac7bd58ea2",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Creating Static Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b11b7777",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import contextily as cx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fab82c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.geometry.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "874740fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.geometry.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "959f0f4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(stops.geometry.x, stops.geometry.y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6110a27c-0c10-4397-ad81-546f7f5a6007",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7326c62-622b-45dc-a8c3-bd99c4390cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = stops.plot(figsize=(10, 10), color='black', markersize=50)\n",
    "cx.add_basemap(ax=ax, crs=stops.crs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06a0a951-f3b2-4461-9270-5c1b583b266d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = stops.plot(figsize=(10, 10), color='black', markersize=50)\n",
    "cx.add_basemap(ax=ax, crs=stops.crs, source=cx.providers.OpenStreetMap.Mapnik)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e30a80c-6e1e-4307-a45c-0150b2ea5a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "sns.kdeplot(x=stops.geometry.x,\n",
    "            y=stops.geometry.y,\n",
    "            fill=True,\n",
    "            alpha=0.75,\n",
    "            cmap='magma',\n",
    "            ax=ax)\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=stops.crs)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d06e5da-d281-4f6f-8df9-1c616f063066",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = stops[stops.accessible == 'yes'].plot(figsize=(10, 10),\n",
    "                                           color='blue',\n",
    "                                           label='accessible',\n",
    "                                           marker='o',\n",
    "                                           markersize=50)\n",
    "\n",
    "stops[stops.accessible == 'no'].plot(color='red',\n",
    "                                     label='not accessible',\n",
    "                                     marker='x',\n",
    "                                     markersize=100,\n",
    "                                     ax=ax)\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=stops.crs)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02b0fdea-7960-417b-bcdc-3077c19ca427",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Creating Interactive Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7521d23b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "import hvplot.pandas\n",
    "import plotly.express as px"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "452c2e3c-8fdc-4769-82fe-8d8a386ddc11",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.explore(column='accessible',\n",
    "              tooltip='stop_name',\n",
    "              tooltip_kwds=dict(labels=False),\n",
    "              marker_kwds=dict(radius=5, fill=True),\n",
    "              popup=True,\n",
    "              legend=True,\n",
    "              cmap='bwr_r',\n",
    "              tiles='CartoDB positron')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "277bd4da-f212-4d55-8c09-f68a532a95dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops_wgs84 = stops.to_crs(epsg=4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58614d97-b4db-4c0a-ae9e-d49d098987c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops_wgs84_ctr = stops_wgs84.unary_union.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80be72b2-15ca-4133-9637-8b2560915ade",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=(stops_wgs84_ctr.y, stops_wgs84_ctr.x),\n",
    "               zoom_start=12)\n",
    "\n",
    "def style_function(feature):\n",
    "    color = 'blue' if feature['properties']['accessible']=='yes' else 'red'\n",
    "    return {'color': color, 'fillColor': color}\n",
    "\n",
    "folium.GeoJson(data=stops_wgs84.to_json(),\n",
    "               style_function = style_function,\n",
    "               tooltip=folium.GeoJsonTooltip(fields=['stop_name'],\n",
    "                                             labels=False),\n",
    "               marker=folium.CircleMarker(radius=5,\n",
    "                                          fill=True,\n",
    "                                          fill_opacity=0.5)).add_to(m)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4695a9b-f6dd-4a2d-947a-543b5ca1d6aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops.hvplot(crs=stops.crs.to_epsg(),\n",
    "             tiles='StamenTerrainRetina',\n",
    "             hover_cols=['stop_name'],\n",
    "             c='accessible',\n",
    "             cmap='bwr',\n",
    "             alpha=0.7,\n",
    "             size=200,\n",
    "             width=1200,\n",
    "             height=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0531ea-60e2-4158-b4b8-f15a53ec7987",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter_mapbox(stops_wgs84,\n",
    "                        lat=stops_wgs84.geometry.y,\n",
    "                        lon=stops_wgs84.geometry.x,\n",
    "                        color='accessible',\n",
    "                        color_discrete_map={'yes': 'blue', 'no': 'red'},\n",
    "                        hover_name='stop_name',\n",
    "                        mapbox_style='carto-positron',\n",
    "                        opacity=0.7,\n",
    "                        zoom=10,\n",
    "                        width=1200,\n",
    "                        height=600)\n",
    "\n",
    "fig.update_traces(marker={'size': 15})\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d56ef1-7baf-4f61-a9ae-0a12ecd7b01c",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Exploring Shapely Objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57b51f36-1189-4302-a843-55f511dc4761",
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon, LineString\n",
    "from shapely.ops import transform\n",
    "import pyproj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dce10e8-f978-4ed1-ac5d-8e939a87027e",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops[stops.stop_name == 'Medford/Tufts']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d688de18-d5f2-4644-a541-706cda2c06e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops[stops.stop_name == 'Medford/Tufts'].geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e68cb9d2-7360-4be8-bfa7-d871f2130a11",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford = stops[stops.stop_name == 'Medford/Tufts'].geometry.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "666c6aa7-d1e8-4d90-bae8-5441523d9f04",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b662b4d-d855-4bad-a221-ab4eeb28459e",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(medford)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a62635",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd60db65",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06c331ed-9b23-4546-a929-3054e1d01615",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(medford.x, medford.y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d8f0377-d19f-4a83-a963-068a8547082c",
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = stops[stops.stop_name == 'Tufts Medical Center'].geometry.values[0]\n",
    "smfa = stops[stops.stop_name == 'Museum of Fine Arts'].geometry.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e257a54-4e08-4ab7-af68-f9c551be4b70",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = stops.plot(figsize=(10, 10), color='black', markersize=25)\n",
    "ax.scatter(medford.x, medford.y, color='blue', marker='*', s=200)\n",
    "ax.scatter(boston.x, boston.y, color='red', marker='X', s=100)\n",
    "ax.scatter(smfa.x, smfa.y, color='green', marker='s', s=50)\n",
    "cx.add_basemap(ax, crs=stops.crs, source=cx.providers.CartoDB.Positron)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6df85dd2-b551-4114-9e4b-308be6234f96",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle = Polygon([medford, boston, smfa])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50b2b0be-e751-4ea0-ae99-ecd2b2e1cab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91d6f7a9-3c8c-494e-a0e0-472b04fca110",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(triangle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7862ebc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.exterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd19123",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(triangle.exterior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4181544-5612-469c-b775-229fe684e0f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "019b85f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.exterior.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ca1cb68",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62009c0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.exterior.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e31e1295-9697-4391-93a2-54bd5859074b",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford.distance(boston)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3759c82",
   "metadata": {},
   "outputs": [],
   "source": [
    "medford.distance(smfa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63a9a6ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle.exterior.xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b86e43b-e4ff-4227-8386-93cd87ee368d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = triangle.exterior.xy\n",
    "plt.plot(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1241ae77",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = triangle.exterior.xy\n",
    "plt.fill(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3540b61-b13c-42a3-966f-342f19f7bf59",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = triangle.exterior.xy\n",
    "ax = stops.plot(figsize=(10, 10), color='black', markersize=50, zorder=3)\n",
    "plt.plot(x, y, color='blue', linewidth=3, zorder=2)\n",
    "plt.fill(x, y, color='blue', alpha=0.5, zorder=1)\n",
    "cx.add_basemap(ax, crs=stops.crs, source=cx.providers.CartoDB.Positron)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5b6ef21-fe4f-4b14-962f-b541b907df0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformer = pyproj.Transformer.from_crs(stops.crs, 'epsg:4326').transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8c8ba8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=transform(transformer, triangle).centroid.coords[0],\n",
    "               tiles='CartoDB positron',\n",
    "               zoom_start=12)\n",
    "folium.Polygon(locations=transform(transformer, triangle.exterior).coords,\n",
    "                color='blue', fill=True).add_to(m)\n",
    "folium.Marker(location=transform(transformer, medford).coords[0],\n",
    "              tooltip='Medford/Tufts').add_to(m)\n",
    "folium.Marker(location=transform(transformer, boston).coords[0],\n",
    "              tooltip='Tufts Medical Center').add_to(m)\n",
    "folium.Marker(location=transform(transformer, smfa).coords[0],\n",
    "              tooltip='Museum of Fine Arts').add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "077ddbf8-0ca1-4e08-847b-8857acf99e14",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Geocoding with GeoPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27580416-c87a-4a7b-9030-ea83e50ab78a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "from geopy.geocoders import Nominatim\n",
    "from geopy.extra.rate_limiter import RateLimiter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ed4fc74-4513-4d3d-b18c-722542fdbf9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx_stops = ['Medford/Tufts', 'Ball Square', 'Magoun Square', 'Gilman Square',\n",
    "             'East Somerville', 'Lechmere', 'Union Square']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a6dac58",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx = pd.read_csv('data/mbta-transit-stations.csv')\n",
    "glx = glx[glx.stop_name.isin(glx_stops)].reset_index(drop=True)\n",
    "glx.drop(columns=['stop_lat', 'stop_lon', 'accessible'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c05c9547",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75c8d3e2-9ce1-4a68-86a3-6528d15a3378",
   "metadata": {},
   "outputs": [],
   "source": [
    "token = 'uep239-python-spatial-{:04}'.format(random.randint(0,9999))\n",
    "print(token)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a717a802-271e-4f43-aa7c-c1aecf49ca4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "geolocator = Nominatim(user_agent=token)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06e5d749-e001-4e01-94f9-beb308695f28",
   "metadata": {},
   "outputs": [],
   "source": [
    "geocode = RateLimiter(geolocator.geocode, min_delay_seconds=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bae66ea-7303-489b-a6f0-08d4acb32d27",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx['location'] = glx.stop_address.apply(geocode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9689f568-8995-444e-ae50-e8385a035780",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aa6d7e2-df40-48bc-9659-1cdb33661e68",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx['stop_lat'] = glx.location.apply(lambda loc: loc.latitude if loc else None)\n",
    "glx['stop_lon'] = glx.location.apply(lambda loc: loc.longitude if loc else None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6d0e098-e4ce-4d2f-9c38-efef8b5ecef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx.drop(columns='location', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5ed3f13-7f23-4161-8c7d-910accdd8b24",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b9931fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx = gpd.GeoDataFrame(data=glx,\n",
    "                       geometry=gpd.points_from_xy(x=glx.stop_lon,\n",
    "                                                   y=glx.stop_lat,\n",
    "                                                   crs='epsg:4326'))\n",
    "glx.to_crs(epsg=26986, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7ffb559",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glx.plot(color='darkgreen', markersize=50)\n",
    "cx.add_basemap(ax=ax, crs=glx.crs, source=cx.providers.CartoDB.Positron)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a5b7d81",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx_line = LineString(glx.set_index('stop_name').loc[glx_stops, 'geometry'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35d0c4dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "glx_line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1648eb03",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = glx_line.xy\n",
    "ax = glx.plot(color='darkgreen', markersize=50, zorder=2)\n",
    "plt.plot(x, y, color='green', zorder=1)\n",
    "cx.add_basemap(ax=ax, crs=glx.crs, source=cx.providers.CartoDB.Positron)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53d9eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=transform(transformer, glx_line).centroid.coords[0],\n",
    "               tiles='CartoDB positron',\n",
    "               zoom_start=13)\n",
    "folium.PolyLine(locations=transform(transformer, glx_line).coords,\n",
    "                color='green').add_to(m)\n",
    "folium.GeoJson(data=glx.to_crs(epsg=4326).to_json(),\n",
    "               tooltip=folium.GeoJsonTooltip(fields=['stop_name'],\n",
    "                                             labels=False),\n",
    "               marker=folium.Marker(icon=folium.Icon(color='green',\n",
    "                                                     icon='train-tram',\n",
    "                                                     prefix='fa',))).add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfe8c993-7506-4999-ba1b-13858103c8c2",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "128eba7f-a94a-4575-9a9e-42ec9d2d7e09",
   "metadata": {},
   "outputs": [],
   "source": [
    "entries = (pd.read_csv('data/mbta-gated-entries-2022.csv')\n",
    "             .groupby('stop_id')\n",
    "             .gated_entries.sum()\n",
    "             .to_frame('total_entries')\n",
    "             .reset_index())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a15a8526-bbb3-444d-8b92-f498d2343ac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "entries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c51bf27-756a-4569-86e4-6199e55f0250",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops = stops.merge(entries, how='left', on='stop_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b609170b",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcbd205b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = stops.plot(figsize=(10, 10),\n",
    "                column='total_entries',\n",
    "                cmap='Greens',\n",
    "                alpha=0.8,\n",
    "                edgecolor='darkgreen',\n",
    "                legend=True,\n",
    "                markersize=100,\n",
    "                missing_kwds=dict(color='gray',\n",
    "                                  edgecolor=None,\n",
    "                                  alpha=0.5,\n",
    "                                  marker='x',\n",
    "                                  markersize=50))\n",
    "cx.add_basemap(ax, crs=stops.crs, source=cx.providers.CartoDB.Positron)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo",
   "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "471545a5e0d746e0c5e432ab2b855d187b4082cc93405d38fc83f4bf06d1f6fe"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}