try sjoin shapes to road segments

Author: tiffanychu90

gtfs_schedule/06_shapes_on_roads.ipynb

@@ -0,0 +1,186 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "2e515aeb-4ac8-4ecb-a09c-23fb86c7bcab",
+   "metadata": {},
+   "source": [
+    "# Spatial join shapes to roads\n",
+    "\n",
+    "* Parallel routes\n",
+    "\n",
+    "Instead of using points (stops, stop arrivals), let's join a bunch of shapes to road segments and see how that turns out."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "876a573f-a163-4218-b6a0-a8d0f23fc677",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dask.dataframe as dd\n",
+    "import dask_geopandas as dg\n",
+    "import geopandas as gpd\n",
+    "import intake\n",
+    "import pandas as pd\n",
+    "\n",
+    "from shared_utils import rt_dates, rt_utils\n",
+    "from segment_speed_utils import helpers\n",
+    "from segment_speed_utils.project_vars import (SCHED_GCS, \n",
+    "                                              SHARED_GCS, \n",
+    "                                              PROJECT_CRS\n",
+    "                                             )\n",
+    "                                             \n",
+    "catalog = intake.open_catalog(\n",
+    "    \"../_shared_utils/shared_utils/shared_data_catalog.yml\")\n",
+    "\n",
+    "analysis_date = rt_dates.DATES[\"sep2023\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "42c8a487-03bb-48ae-8833-032962b7f62c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "primary_secondary = dg.read_parquet(\n",
+    "    f\"{SHARED_GCS}segmented_roads_2020_primarysecondary.parquet\"\n",
+    ")\n",
+    "\n",
+    "local = dg.read_parquet(\n",
+    "    f\"{SHARED_GCS}segmented_roads_2020_local.parquet\",\n",
+    "    columns = [i for i in primary_secondary.columns]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e8fb012d-c94b-4d47-9a2c-7af88045d6c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "road_segments = dd.multi.concat(\n",
+    "    [primary_secondary, local], \n",
+    "    axis=0\n",
+    ").reset_index(drop=True).repartition(npartitions=10)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d3ef08ab-ba37-47d6-a767-9515aa657336",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "shapes = helpers.import_scheduled_shapes(\n",
+    "    analysis_date,\n",
+    "    columns = [\"shape_array_key\", \"n_trips\", \"geometry\"],\n",
+    "    get_pandas = True,\n",
+    "    crs = PROJECT_CRS\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c93729ab-fb86-41f3-86f1-ca021e506b14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def spatial_join_by_partition(road_segments, shapes):\n",
+    "    road_segments_buff = road_segments.assign(\n",
+    "        geometry = road_segments.geometry.buffer(35)\n",
+    "    )\n",
+    "    \n",
+    "    s1 = gpd.sjoin(\n",
+    "        road_segments_buff,\n",
+    "        shapes,\n",
+    "        how = \"inner\",\n",
+    "        predicate = \"intersects\"\n",
+    "    ).drop(columns = [\"index_right\", \"geometry\"])\n",
+    "        \n",
+    "    return s1\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "63711923-9035-4597-b162-19b2dbece88d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "road_cols = [\n",
+    "    \"linearid\", \"mtfcc\", \n",
+    "    \"fullname\", \"segment_sequence\", \n",
+    "    \"primary_direction\"\n",
+    "]\n",
+    "\n",
+    "road_dtypes = road_segments[road_cols].dtypes.to_dict()\n",
+    "shape_dtypes = shapes[[\"shape_array_key\", \"n_trips\"]].dtypes.to_dict()\n",
+    "\n",
+    "sjoin_results = road_segments.map_partitions(\n",
+    "    spatial_join_by_partition,\n",
+    "    shapes,\n",
+    "    meta = {\n",
+    "        **road_dtypes,\n",
+    "        **shape_dtypes\n",
+    "    },\n",
+    "    align_dataframes = False,\n",
+    ").persist()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cbd685c6-22ba-4d29-bfaa-c46a6e7c5b7e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "trips_per_segment = (sjoin_results.groupby(road_cols, \n",
+    "                                           observed=True, group_keys=False)\n",
+    "                     .agg({\n",
+    "                         \"shape_array_key\": \"count\",\n",
+    "                         \"n_trips\": \"sum\"\n",
+    "                     }).reset_index()\n",
+    "                    ).compute()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd3e3963-78f3-489e-8bc2-8d25185d3c1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "trips_per_segment.to_parquet(\n",
+    "    f\"{SCHED_GCS}trips_road_segments_{analysis_date}.parquet\"\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}