analysis/analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "import psutil\n",
    "import sys\n",
    "\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "from railrailrail.network.stage import Stage\n",
    "from railrailrail.railgraph import RailGraph\n",
    "\n",
    "from analysis.analysis_helpers import (\n",
    "    get_stage_journeys,\n",
    "    make_stage_journeys_dataframe,\n",
    "    contains_point,\n",
    "    get_station_agg_stats,\n",
    "    get_planning_area_to_region_df,\n",
    ")\n",
    "\n",
    "# Do not truncate DataFrames\n",
    "pd.set_option(\"display.max_rows\", None)\n",
    "pd.set_option(\"display.max_columns\", None)\n",
    "pd.set_option(\"display.width\", None)\n",
    "pd.set_option(\"display.max_colwidth\", None)\n",
    "\n",
    "minimum_available_ram_GB = 3\n",
    "if psutil.virtual_memory().available < minimum_available_ram_GB * 1e9:\n",
    "    print(\n",
    "        \"Warning: Less than %dGB RAM available. This program may crash with an Out-Of-Memory error.\"\n",
    "        % minimum_available_ram_GB,\n",
    "        file=sys.stderr,\n",
    "    )\n",
    "\n",
    "\n",
    "planning_areas = gpd.read_file(\"MasterPlan2019PlanningAreaBoundaryNoSea.geojson\")\n",
    "planning_areas[\"Planning Area\"] = planning_areas[\"Description\"].apply(\n",
    "    lambda description: description.split(\"<td>\", 1)[-1].split(\"</td>\", 1)[0].strip()\n",
    ")\n",
    "\n",
    "planning_area_to_region_df = get_planning_area_to_region_df()\n",
    "\n",
    "\n",
    "stages = list(Stage.stages)\n",
    "journeys_by_stage: dict[str, dict[str, pd.DataFrame | RailGraph]] = dict()\n",
    "station_agg_stats_by_stage: dict[str, pd.DataFrame] = dict()\n",
    "\n",
    "\n",
    "with multiprocessing.get_context(\"fork\").Pool(\n",
    "    None\n",
    ") as pool:  # https://bugs.python.org/issue33725 https://stackoverflow.com/a/72605750\n",
    "    for stage, journeys, rail_graph in pool.imap_unordered(get_stage_journeys, stages):\n",
    "        journeys_by_stage[stage] = {\n",
    "            \"df\": make_stage_journeys_dataframe(stage, journeys),\n",
    "            \"rail_graph\": rail_graph,\n",
    "        }\n",
    "        station_agg_stats_by_stage[stage] = get_station_agg_stats(\n",
    "            df=journeys_by_stage[stage][\"df\"],\n",
    "            rail_graph=journeys_by_stage[stage][\"rail_graph\"],\n",
    "        )\n",
    "        # Stations within interchanges have different geographic coordinates.\n",
    "        # This causes them to have different average circuity. For example, see BP1/NS4/JS1 Choa Chu Kang.\n",
    "\n",
    "        station_agg_stats_by_stage[stage][\"planning_area\"] = station_agg_stats_by_stage[\n",
    "            stage\n",
    "        ].apply(\n",
    "            lambda row: next(\n",
    "                (\n",
    "                    feature[\"Planning Area\"]\n",
    "                    for _, feature in planning_areas.iterrows()\n",
    "                    if contains_point(\n",
    "                        feature[\"geometry\"],\n",
    "                        row[\"latitude\"],\n",
    "                        row[\"longitude\"],\n",
    "                    )\n",
    "                ),\n",
    "                None,\n",
    "            ),\n",
    "            axis=1,\n",
    "        )\n",
    "\n",
    "        station_agg_stats_by_stage[stage] = station_agg_stats_by_stage[stage].merge(\n",
    "            planning_area_to_region_df, on=\"planning_area\", how=\"left\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Journeys between only stations that were around in 2003.\n",
    "old_stations = list(journeys_by_stage[\"nel\"][\"df\"][\"start\"].unique())\n",
    "old_stations_df = pd.DataFrame({\"station\": old_stations})\n",
    "\n",
    "recent_stages = set()\n",
    "start = False\n",
    "for stage in journeys_by_stage.keys():\n",
    "    if stage == \"nel\":\n",
    "        start = True\n",
    "    if not start:\n",
    "        continue\n",
    "    recent_stages.add(stage)\n",
    "\n",
    "journeys_by_stage_from_nel: dict[str, dict[str, pd.DataFrame | RailGraph]] = dict()\n",
    "station_agg_stats_by_stage_from_nel: dict[str, pd.DataFrame] = dict()\n",
    "\n",
    "\n",
    "def get_stage_journeys_cumulative(stage: str):\n",
    "    if stage not in recent_stages:\n",
    "        return None\n",
    "    return (\n",
    "        journeys_by_stage[stage][\"df\"]\n",
    "        .merge(old_stations_df, left_on=\"start\", right_on=\"station\")\n",
    "        .merge(old_stations_df, left_on=\"end\", right_on=\"station\")\n",
    "    )\n",
    "\n",
    "\n",
    "for stage in stages:\n",
    "    df = get_stage_journeys_cumulative(stage)\n",
    "    if df is None:\n",
    "        continue\n",
    "    journeys_by_stage_from_nel[stage] = dict()\n",
    "    journeys_by_stage_from_nel[stage][\"rail_graph\"] = journeys_by_stage[stage][\n",
    "        \"rail_graph\"\n",
    "    ]\n",
    "    journeys_by_stage_from_nel[stage][\"df\"] = df\n",
    "\n",
    "    station_agg_stats_by_stage_from_nel[stage] = get_station_agg_stats(\n",
    "        df=journeys_by_stage_from_nel[stage][\"df\"],\n",
    "        rail_graph=journeys_by_stage_from_nel[stage][\"rail_graph\"],\n",
    "    )\n",
    "    # Stations within interchanges have different geographic coordinates.\n",
    "    # This causes them to have different average circuity. For example, see BP1/NS4 Choa Chu Kang.\n",
    "\n",
    "    station_agg_stats_by_stage_from_nel[stage][\"planning_area\"] = (\n",
    "        station_agg_stats_by_stage_from_nel[stage].apply(\n",
    "            lambda row: next(\n",
    "                (\n",
    "                    feature[\"Planning Area\"]\n",
    "                    for _, feature in planning_areas.iterrows()\n",
    "                    if contains_point(\n",
    "                        feature[\"geometry\"],\n",
    "                        row[\"latitude\"],\n",
    "                        row[\"longitude\"],\n",
    "                    )\n",
    "                ),\n",
    "                None,\n",
    "            ),\n",
    "            axis=1,\n",
    "        )\n",
    "    )\n",
    "\n",
    "    station_agg_stats_by_stage_from_nel[stage] = station_agg_stats_by_stage_from_nel[\n",
    "        stage\n",
    "    ].merge(planning_area_to_region_df, on=\"planning_area\", how=\"left\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_planning_area_agg_stats_scatter(\n",
    "    station_agg_stats_by_stage: dict[str, pd.DataFrame], stage: str, show=False\n",
    "):\n",
    "    stage_description, stage_timestamp = Stage.stages_info[stage]\n",
    "    planning_area_agg_stats = (\n",
    "        station_agg_stats_by_stage[stage]\n",
    "        .copy()\n",
    "        .groupby(by=[\"planning_area\", \"region\"])\n",
    "        .mean(numeric_only=True)\n",
    "        .reset_index()\n",
    "        .sort_values(by=[\"path_distance_mean\"])\n",
    "    )\n",
    "    planning_area_agg_stats[\"size\"] = 10  # Chart dot size\n",
    "    # print(\n",
    "    #     \"Number of active planning areas:\",\n",
    "    #     planning_area_agg_stats[\"planning_area\"].nunique(),\n",
    "    # )\n",
    "\n",
    "    fig = px.scatter(\n",
    "        planning_area_agg_stats,\n",
    "        title=(\n",
    "            \"Mean Circuity against Mean Fastest Path Length (Averaged By Planning Area)\"\n",
    "            \"<br />\"\n",
    "            f\"{stage_description} @ {stage_timestamp.strftime('%d %B %Y')}\"\n",
    "        ),\n",
    "        x=\"path_distance_mean\",\n",
    "        y=\"circuity_mean\",\n",
    "        range_x=(\n",
    "            int(planning_area_agg_stats[\"path_distance_mean\"].min() * 0.95),\n",
    "            int(planning_area_agg_stats[\"path_distance_mean\"].max() * 1.05),\n",
    "        ),\n",
    "        range_y=(\n",
    "            planning_area_agg_stats[\"circuity_mean\"].min() - 0.1,\n",
    "            planning_area_agg_stats[\"circuity_mean\"].max() + 0.1,\n",
    "        ),\n",
    "        labels={\n",
    "            \"path_distance_mean\": \"Mean Fastest Path Length\",\n",
    "            \"circuity_mean\": \"Mean Circuity\",\n",
    "            \"planning_area\": \"Planning Area\",\n",
    "            \"region\": \"Region\",\n",
    "        },\n",
    "        size=\"size\",\n",
    "        color=\"region\",\n",
    "        color_discrete_map={\n",
    "            \"CENTRAL\": \"#ffa15a\",\n",
    "            \"WEST\": \"#636efa\",\n",
    "            \"EAST\": \"#00cc96\",\n",
    "            \"NORTH-EAST\": \"fuchsia\",\n",
    "            \"NORTH\": \"#ef553b\",\n",
    "        },\n",
    "        text=\"planning_area\",\n",
    "        hover_name=\"planning_area\",\n",
    "        hover_data={\"size\": False},\n",
    "        # trendline=\"lowess\",  # Locally WEighted Scatterplot Smoothing (LOWESS)\n",
    "        # trendline_scope=\"overall\",\n",
    "        log_x=False,\n",
    "        log_y=False,\n",
    "        size_max=10,\n",
    "        template=\"plotly_dark\",\n",
    "    )\n",
    "\n",
    "    fig.add_hline(\n",
    "        y=planning_area_agg_stats[\"circuity_mean\"].mean(),\n",
    "        line_width=2,\n",
    "        line_dash=\"dash\",\n",
    "        line_color=\"red\",\n",
    "    )\n",
    "    fig.add_vline(\n",
    "        x=planning_area_agg_stats[\"path_distance_mean\"].mean(),\n",
    "        line_width=2,\n",
    "        line_dash=\"dash\",\n",
    "        line_color=\"red\",\n",
    "    )\n",
    "\n",
    "    fig.add_annotation(\n",
    "        x=planning_area_agg_stats[\"path_distance_mean\"].max(),\n",
    "        y=planning_area_agg_stats[\"circuity_mean\"].mean(),\n",
    "        text=f\"Mean: {planning_area_agg_stats['circuity_mean'].mean():.2f}\",\n",
    "        showarrow=False,\n",
    "        yshift=10,\n",
    "    )  # Median Mean Circuity across all Planning Areas\n",
    "\n",
    "    fig.add_annotation(\n",
    "        x=planning_area_agg_stats[\"path_distance_mean\"].mean(),\n",
    "        y=planning_area_agg_stats[\"circuity_mean\"].max(),\n",
    "        text=f\"Mean: {planning_area_agg_stats['path_distance_mean'].mean() / 1000:.1f}km\",\n",
    "        showarrow=False,\n",
    "        xshift=50,\n",
    "    )  # Median Mean Fastest Path Length across all Planning Areas\n",
    "\n",
    "    fig.update_layout(\n",
    "        title={\"x\": 0.5},\n",
    "        dragmode=\"pan\",\n",
    "        xaxis=dict(\n",
    "            tickformat=\"\",\n",
    "            ticksuffix=\"m\",\n",
    "            tickfont=dict(color=\"white\", family=\"Arial Black\"),\n",
    "        ),\n",
    "        yaxis=dict(\n",
    "            tickfont=dict(color=\"white\", family=\"Arial Black\"),\n",
    "        ),\n",
    "    )\n",
    "    fig.update_traces(\n",
    "        textposition=\"top center\",\n",
    "        textfont=dict(size=13, color=\"white\", family=\"Arial\"),\n",
    "    )\n",
    "\n",
    "    if show:\n",
    "        fig.show()\n",
    "    return planning_area_agg_stats, fig\n",
    "\n",
    "\n",
    "for stage in [\"teck_lee\", \"cg_tel_c\"]:\n",
    "    planning_area_agg_stats, fig = plot_planning_area_agg_stats_scatter(\n",
    "        station_agg_stats_by_stage, stage, show=True\n",
    "    )\n",
    "    print(fig.to_html(full_html=False, include_plotlyjs=\"cdn\"))\n",
    "    print()\n",
    "\n",
    "for stage in [\"teck_lee\", \"cg_tel_c\"]:\n",
    "    planning_area_agg_stats, fig = plot_planning_area_agg_stats_scatter(\n",
    "        station_agg_stats_by_stage_from_nel, stage, show=True\n",
    "    )\n",
    "    print(fig.to_html(full_html=False, include_plotlyjs=\"cdn\"))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = journeys_by_stage[\"phase_1_2\"][\"df\"]\n",
    "print(\n",
    "    \"%.1f%% of journeys as of Phase 1b are between adjacent stations.\"\n",
    "    % (len(df.query(\"total_cost != 0 and circuity == 1.0\")) / len(df) * 100,)\n",
    ")\n",
    "\n",
    "df = journeys_by_stage[\"teck_lee\"][\"df\"]\n",
    "print(\n",
    "    \"%.1f%% of journeys as of Teck Lee Station launch are between adjacent stations.\"\n",
    "    % (len(df.query(\"total_cost != 0 and circuity == 1.0\")) / len(df) * 100,)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_agg_stats_map_multi(\n",
    "    station_agg_stats_by_stage: dict[str, pd.DataFrame],\n",
    "    stages: list[str],\n",
    "    color: str,\n",
    "    hover_data: list[str],\n",
    "    show: bool = False,\n",
    "):\n",
    "    station_agg_stats = pd.DataFrame()\n",
    "    for stage in stages:\n",
    "        df = station_agg_stats_by_stage[stage].copy()\n",
    "        df[\"full_name\"] = df.apply(\n",
    "            lambda row: row[\"start\"] + \" \" + row[\"station_name\"], axis=1\n",
    "        )\n",
    "        df[\"stage_description\"] = Stage.stages_info[stage][0]\n",
    "        df[\"stage_timestamp\"] = Stage.stages_info[stage][1].strftime(\"%d %B %Y\")\n",
    "        station_agg_stats = pd.concat([station_agg_stats, df], axis=0)\n",
    "\n",
    "    max_circuity_mean = station_agg_stats[\"circuity_mean\"].max()\n",
    "    max_circuity_median = station_agg_stats[\"circuity_median\"].max()\n",
    "    max_circuity = max(max_circuity_mean, max_circuity_median)\n",
    "\n",
    "    max_path_distance_mean = station_agg_stats[\"path_distance_mean\"].max()\n",
    "    max_path_distance_median = station_agg_stats[\"path_distance_median\"].max()\n",
    "    max_path_distance = max(max_path_distance_mean, max_path_distance_median)\n",
    "\n",
    "    station_agg_stats[\"marker_size\"] = 1\n",
    "    fig = px.scatter_mapbox(\n",
    "        station_agg_stats,\n",
    "        title=f\"{'Mean Circuity' if 'circuity' in color else 'Mean Fastest Path Length'}\",\n",
    "        lat=\"latitude\",\n",
    "        lon=\"longitude\",\n",
    "        animation_frame=\"stage_timestamp\",\n",
    "        hover_name=\"full_name\",\n",
    "        hover_data=hover_data,\n",
    "        color=color,\n",
    "        color_continuous_scale=px.colors.sequential.Jet,\n",
    "        range_color=(\n",
    "            (1, max_circuity) if \"circuity\" in color else (8000, max_path_distance)\n",
    "        ),\n",
    "        mapbox_style=\"carto-darkmatter\",\n",
    "        template=\"plotly_dark\",\n",
    "        zoom=9,\n",
    "        labels={\n",
    "            \"path_distance_mean\": \"Mean Fastest Path Length\",\n",
    "            \"circuity_mean\": \"Mean Circuity\",\n",
    "            \"planning_area\": \"Planning Area\",\n",
    "            \"region\": \"Region\",\n",
    "        },\n",
    "        size=\"marker_size\",\n",
    "        size_max=9,\n",
    "        opacity=1,\n",
    "        center={\"lat\": 1.330270, \"lon\": 103.851959},\n",
    "    )\n",
    "    fig[\"layout\"].pop(\"updatemenus\")\n",
    "    fig.update_layout(\n",
    "        title={\"x\": 0.5, \"y\": 0.95},\n",
    "        margin=dict(t=50, b=0),\n",
    "        sliders=[\n",
    "            {\n",
    "                \"currentvalue\": {\n",
    "                    \"font\": {\"size\": 20},\n",
    "                    \"prefix\": None,\n",
    "                    \"visible\": True,\n",
    "                    \"xanchor\": \"center\",\n",
    "                },\n",
    "                \"transition\": {\"duration\": 0, \"easing\": \"linear\"},\n",
    "                \"pad\": {\"t\": 0, \"l\": 0},\n",
    "            }\n",
    "        ],\n",
    "    )\n",
    "    if \"circuity\" in color:\n",
    "        fig.update_coloraxes(colorbar_title=\"Mean Circuity\")\n",
    "    else:\n",
    "        fig.update_coloraxes(\n",
    "            colorbar_ticksuffix=\"m\", colorbar_title=\"Mean Fastest Path Length\"\n",
    "        )\n",
    "    fig.update_coloraxes(showscale=True)\n",
    "\n",
    "    if show:\n",
    "        fig.show()\n",
    "    return station_agg_stats, fig\n",
    "\n",
    "\n",
    "for agg_stats_by_stage in [\n",
    "    station_agg_stats_by_stage,\n",
    "    station_agg_stats_by_stage_from_nel,\n",
    "]:\n",
    "    for color in [\"path_distance_mean\", \"circuity_mean\"]:\n",
    "        station_agg_stats, fig = plot_agg_stats_map_multi(\n",
    "            agg_stats_by_stage,\n",
    "            [\n",
    "                \"nel\",\n",
    "                \"pglrt_east_loop_and_sklrt_west_loop\",\n",
    "                \"dtl_2\",\n",
    "                \"teck_lee\",\n",
    "                \"cg_tel_c\",\n",
    "            ],\n",
    "            color=color,\n",
    "            hover_data={\n",
    "                \"circuity_mean\": True,\n",
    "                \"path_distance_mean\": True,\n",
    "                \"region\": True,\n",
    "                \"planning_area\": True,\n",
    "                \"marker_size\": False,\n",
    "                \"stage_timestamp\": False,\n",
    "            },\n",
    "            show=True,\n",
    "        )\n",
    "        print(fig.to_html(full_html=False, include_plotlyjs=\"cdn\", auto_play=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_agg_stats_map(\n",
    "    stage: str, color: str, hover_data: list[str], show: bool = False\n",
    "):\n",
    "    station_agg_stats = station_agg_stats_by_stage[stage].copy()\n",
    "    station_agg_stats[\"full_name\"] = station_agg_stats.apply(\n",
    "        lambda row: row[\"start\"] + \" \" + row[\"station_name\"], axis=1\n",
    "    )\n",
    "    stage_description = Stage.stages_info[stage][0]\n",
    "    stage_timestamp = Stage.stages_info[stage][1]\n",
    "\n",
    "    max_circuity_mean = station_agg_stats[\"circuity_mean\"].max()\n",
    "    max_circuity_median = station_agg_stats[\"circuity_median\"].max()\n",
    "    max_circuity = max(max_circuity_mean, max_circuity_median)\n",
    "\n",
    "    max_path_distance_mean = station_agg_stats[\"path_distance_mean\"].max()\n",
    "    max_path_distance_median = station_agg_stats[\"path_distance_median\"].max()\n",
    "    max_path_distance = max(max_path_distance_mean, max_path_distance_median)\n",
    "\n",
    "    station_agg_stats[\"marker_size\"] = 1\n",
    "    fig = px.scatter_mapbox(\n",
    "        station_agg_stats,\n",
    "        title=(\n",
    "            f\"{'Mean Circuity' if 'circuity' in color else 'Mean Fastest Path Length'}\"\n",
    "            \"<br />\"\n",
    "            f\"{stage_description} @ {stage_timestamp.strftime('%d %B %Y')}\"\n",
    "        ),\n",
    "        lat=\"latitude\",\n",
    "        lon=\"longitude\",\n",
    "        hover_name=\"full_name\",\n",
    "        hover_data=hover_data,\n",
    "        color=color,\n",
    "        color_continuous_scale=px.colors.sequential.Jet,\n",
    "        range_color=(\n",
    "            (1, max_circuity) if \"circuity\" in color else (8000, max_path_distance)\n",
    "        ),\n",
    "        mapbox_style=\"carto-darkmatter\",\n",
    "        template=\"plotly_dark\",\n",
    "        zoom=9,\n",
    "        labels={\n",
    "            \"path_distance_mean\": \"Mean Fastest Path Length\",\n",
    "            \"circuity_mean\": \"Mean Circuity\",\n",
    "            \"planning_area\": \"Planning Area\",\n",
    "            \"region\": \"Region\",\n",
    "        },\n",
    "        size=\"marker_size\",\n",
    "        size_max=9,\n",
    "        opacity=1,\n",
    "        center={\"lat\": 1.330270, \"lon\": 103.851959},\n",
    "    )\n",
    "    fig[\"layout\"].pop(\"updatemenus\")\n",
    "    fig.update_layout(\n",
    "        title={\"x\": 0.5},\n",
    "    )\n",
    "    if \"circuity\" in color:\n",
    "        fig.update_coloraxes(colorbar_title=\"Mean Circuity\")\n",
    "    else:\n",
    "        fig.update_coloraxes(\n",
    "            colorbar_ticksuffix=\"m\", colorbar_title=\"Mean Fastest Path Length\"\n",
    "        )\n",
    "    fig.update_coloraxes(showscale=True)\n",
    "\n",
    "    if show:\n",
    "        fig.show()\n",
    "    return station_agg_stats, fig\n",
    "\n",
    "\n",
    "station_agg_stats, fig = plot_agg_stats_map(\n",
    "    \"teck_lee\",\n",
    "    color=\"circuity_mean\",\n",
    "    hover_data={\n",
    "        \"circuity_mean\": True,\n",
    "        \"path_distance_mean\": True,\n",
    "        \"region\": True,\n",
    "        \"planning_area\": True,\n",
    "        \"marker_size\": False,\n",
    "    },\n",
    "    show=True,\n",
    ")\n",
    "print(fig.to_html(full_html=False, include_plotlyjs=\"cdn\", auto_play=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_agg_stats_grid(color: str, hover_data: list[str], show: bool = False):\n",
    "    station_agg_stats = pd.DataFrame()\n",
    "    for stage in stages:\n",
    "        df = station_agg_stats_by_stage[stage].copy()\n",
    "        df[\"full_name\"] = df.apply(\n",
    "            lambda row: row[\"start\"] + \" \" + row[\"station_name\"], axis=1\n",
    "        )\n",
    "        df[\"stage_description\"] = Stage.stages_info[stage][0]\n",
    "        df[\"stage_timestamp\"] = Stage.stages_info[stage][1].strftime(\"%b %Y\")\n",
    "        station_agg_stats = pd.concat([station_agg_stats, df], axis=0)\n",
    "\n",
    "    max_circuity_mean = station_agg_stats[\"circuity_mean\"].max()\n",
    "    max_circuity_median = station_agg_stats[\"circuity_median\"].max()\n",
    "    max_circuity = max(max_circuity_mean, max_circuity_median)\n",
    "\n",
    "    max_path_distance_mean = station_agg_stats[\"path_distance_mean\"].max()\n",
    "    max_path_distance_median = station_agg_stats[\"path_distance_median\"].max()\n",
    "    max_path_distance = max(max_path_distance_mean, max_path_distance_median)\n",
    "\n",
    "    station_agg_stats[\"marker_size\"] = 1\n",
    "    fig = px.scatter(\n",
    "        station_agg_stats,\n",
    "        title=f\"{'Mean Circuity' if 'circuity' in color else 'Mean Fastest Path Length'}\",\n",
    "        y=\"latitude\",\n",
    "        x=\"longitude\",\n",
    "        range_y=[1.25, 1.46],\n",
    "        range_x=[103.625, 104.02],\n",
    "        animation_frame=\"stage_timestamp\",\n",
    "        hover_name=\"full_name\",\n",
    "        hover_data=hover_data,\n",
    "        color=color,\n",
    "        color_continuous_scale=px.colors.sequential.Jet,\n",
    "        range_color=(\n",
    "            (1, max_circuity) if \"circuity\" in color else (8000, max_path_distance)\n",
    "        ),\n",
    "        template=\"plotly_dark\",\n",
    "        size=\"marker_size\",\n",
    "        size_max=8,\n",
    "        opacity=1,\n",
    "    )\n",
    "    fig[\"layout\"].pop(\"updatemenus\")\n",
    "    fig.update_layout(\n",
    "        title={\"x\": 0.5},\n",
    "        xaxis=dict(scaleanchor=\"y\", scaleratio=1),\n",
    "        dragmode=\"pan\",\n",
    "        sliders=[\n",
    "            {\n",
    "                \"currentvalue\": {\n",
    "                    \"font\": {\"size\": 20},\n",
    "                    \"prefix\": None,\n",
    "                    \"visible\": True,\n",
    "                    \"xanchor\": \"center\",\n",
    "                },\n",
    "                \"transition\": {\"duration\": 0, \"easing\": \"linear\"},\n",
    "                \"pad\": {\"t\": 50, \"l\": 5},\n",
    "            }\n",
    "        ],\n",
    "    )\n",
    "    if \"circuity\" in color:\n",
    "        fig.update_coloraxes(colorbar_title=\"Mean Circuity\")\n",
    "    else:\n",
    "        fig.update_coloraxes(\n",
    "            colorbar_ticksuffix=\"m\", colorbar_title=\"Mean Fastest Path Length\"\n",
    "        )\n",
    "    fig.update_coloraxes(showscale=True)\n",
    "\n",
    "    # Change initial frame\n",
    "    frame_num = next(\n",
    "        i for i, frame in enumerate(fig.frames) if frame[\"name\"] == \"Aug 2024\"\n",
    "    )\n",
    "    fig.layout[\"sliders\"][0][\"active\"] = frame_num\n",
    "    fig = go.Figure(\n",
    "        data=fig[\"frames\"][frame_num][\"data\"], frames=fig[\"frames\"], layout=fig.layout\n",
    "    )\n",
    "\n",
    "    if show:\n",
    "        fig.show(config={\"scrollZoom\": True})\n",
    "    return station_agg_stats, fig\n",
    "\n",
    "\n",
    "for color in [\"path_distance_mean\", \"circuity_mean\"]:\n",
    "    station_agg_stats, fig = plot_agg_stats_grid(\n",
    "        color=color,\n",
    "        hover_data=[\"circuity_mean\", \"path_distance_mean\", \"region\", \"planning_area\"],\n",
    "        show=True,\n",
    "    )\n",
    "    print(fig.to_html(full_html=False, include_plotlyjs=\"cdn\", auto_play=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for stage in journeys_by_stage:\n",
    "    stage_df = journeys_by_stage[stage][\"df\"].query(\"total_cost != 0\")\n",
    "    stage_mean_total_cost = int(stage_df[\"total_cost\"].mean())\n",
    "    stage_max_total_cost = int(stage_df[\"total_cost\"].max())\n",
    "    stage_total_cost_std_dev = int(stage_df[\"total_cost\"].std())\n",
    "    print(stage, stage_mean_total_cost, stage_max_total_cost, stage_total_cost_std_dev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "median_price_per_sq_ft_per_lease_year = pd.DataFrame.from_dict(\n",
    "    {\n",
    "        \"ANG MO KIO\": 9.79,\n",
    "        \"BEDOK\": 9.32,\n",
    "        \"BISHAN\": 11.27,\n",
    "        \"BUKIT BATOK\": 7.33,\n",
    "        \"BUKIT MERAH\": 10.8,\n",
    "        \"BUKIT PANJANG\": 7.08,\n",
    "        \"BUKIT TIMAH\": 12.86,\n",
    "        \"CENTRAL AREA\": 13.2,  # Not a planning area.\n",
    "        \"CHOA CHU KANG\": 6.44,\n",
    "        \"CLEMENTI\": 9.95,\n",
    "        \"GEYLANG\": 10.53,\n",
    "        \"HOUGANG\": 8,\n",
    "        \"JURONG EAST\": 8.08,\n",
    "        \"JURONG WEST\": 6.8,\n",
    "        \"KALLANG\": 10.69,\n",
    "        \"MARINE PARADE\": 12.5,\n",
    "        \"PASIR RIS\": 7.63,\n",
    "        \"PUNGGOL\": 7.11,\n",
    "        \"QUEENSTOWN\": 11.22,\n",
    "        \"SEMBAWANG\": 6.54,\n",
    "        \"SENGKANG\": 6.97,\n",
    "        \"SERANGOON\": 9.74,\n",
    "        \"TAMPINES\": 8.34,\n",
    "        \"TOA PAYOH\": 10.77,\n",
    "        \"WOODLANDS\": 6.76,\n",
    "        \"YISHUN\": 7.69,\n",
    "    },\n",
    "    orient=\"index\",\n",
    ").reset_index()\n",
    "median_price_per_sq_ft_per_lease_year.columns = [\n",
    "    \"planning_area\",\n",
    "    \"median_price_per_sq_ft_per_lease_year\",\n",
    "]\n",
    "\n",
    "planning_area_agg_stats = (\n",
    "    station_agg_stats_by_stage[\"teck_lee\"]\n",
    "    .copy()\n",
    "    .groupby(by=[\"planning_area\", \"region\"])\n",
    "    .mean(numeric_only=True)\n",
    "    .reset_index()\n",
    "    .sort_values(by=[\"path_distance_mean\"])\n",
    "    .merge(median_price_per_sq_ft_per_lease_year, on=\"planning_area\", how=\"inner\")\n",
    ")\n",
    "\n",
    "hover_data = {\n",
    "    \"haversine_distance_median\": True,\n",
    "    \"path_distance_median\": True,\n",
    "    \"region\": True,\n",
    "    \"planning_area\": True,\n",
    "}\n",
    "\n",
    "fig = px.scatter_3d(\n",
    "    planning_area_agg_stats,\n",
    "    x=\"haversine_distance_median\",\n",
    "    y=\"path_distance_median\",\n",
    "    z=\"median_price_per_sq_ft_per_lease_year\",\n",
    "    color=\"region\",\n",
    "    range_x=[\n",
    "        planning_area_agg_stats[\"haversine_distance_median\"].max(),\n",
    "        planning_area_agg_stats[\"haversine_distance_median\"].min(),\n",
    "    ],\n",
    "    text=\"planning_area\",\n",
    "    hover_name=\"planning_area\",\n",
    "    hover_data=hover_data,\n",
    "    height=850,\n",
    ")\n",
    "fig.update_layout(\n",
    "    scene=dict(\n",
    "        xaxis_title=\"Median Haversine (km)\",\n",
    "        yaxis_title=\"Median Path (km)\",\n",
    "        zaxis_title=\"$/sqft/yr\",\n",
    "    )\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}