scripts/REFERENCE/2b_intrinsic_analysis_reference.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<img src=\"../../images/BikeDNA_logo.svg\" width=\"250\"  alt=\"BikeDNA logo\" style=\"display:block; margin-left: auto; margin-right: auto;\">\n",
    "<a href=\"https://github.com/anerv/BikeDNA\">Github</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2b. Intrinsic Analysis of Reference Data\n",
    "\n",
    "This notebook analyses the quality of a user-provided reference bicycle infrastructure data set for a given area. The quality assessment is *intrinsic*, i.e. based only on the one input data set, and making no use of information external to the data set. For an extrinsic quality assessment that compares the reference data set to corresponding OSM data, see the notebooks 3a and 3b.\n",
    "\n",
    "The analysis assesses the *fitness for purpose* ([Barron et al., 2014](https://onlinelibrary.wiley.com/doi/10.1111/tgis.12073)) of the reference data for a given area. Outcomes of the analysis can be relevant for bicycle planning and research - especially for projects that include a network analysis of bicycle infrastructure, in which case the topology of the geometries is of particular importance.\n",
    "\n",
    "Since the assessment does not make use of an external reference data set as the ground truth, no universal claims of data quality can be made. The idea is rather to enable those working with bicycle networks to assess whether their data are good enough for their particular use case. The analysis assists in finding potential data quality issues but leaves the final interpretation of the results to the user.\n",
    "\n",
    "The notebook makes use of quality metrics from a range of previous projects investigating OSM/VGI data quality, such as [Ferster et al. (2020)](https://www.tandfonline.com/doi/full/10.1080/15568318.2018.1519746), [Hochmair et al. (2015)](https://onlinelibrary.wiley.com/doi/abs/10.1111/tgis.12081), [Barron et al. (2014)](https://onlinelibrary.wiley.com/doi/10.1111/tgis.12073), and [Neis et al. (2012](https://www.mdpi.com/1999-5903/4/1/1))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Prerequisites &amp; Input/Output</b>\n",
    "<p>\n",
    "Run notebook 1b in advance.  \n",
    "    \n",
    "Output files of this notebook (data, maps, plots) are saved to the <span style=\"font-family:courier;\">../results/REFERENCE/[study_area]/</span> subfolders.\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Familiarity required</b>\n",
    "<p>\n",
    "For a correct interpretation of some of the metrics for spatial data quality, some familiarity with the area is necessary.\n",
    "</p>\n",
    "</div>\n",
    "\n",
    "* [Data completeness](#Data-completeness)\n",
    "    * [Network density](#Network-density)\n",
    "* [Network topology](#Network-topology)\n",
    "    * [Simplification outcome](#Simplification-outcome)\n",
    "    * [Dangling nodes](#Dangling-nodes)\n",
    "    * [Under/overshoots](#Under/overshoots)\n",
    "* [Network components](#Network-components)\n",
    "    * [Disconnected components](#Disconnected-components)\n",
    "    * [Components per grid cell](#Components-per-grid-cell)\n",
    "    * [Component length distribution](#Component-length-distribution)\n",
    "    * [Largest connected component](#Largest-connected-component)\n",
    "    * [Missing links](#Missing-links)\n",
    "    * [Component connectivity](#Component-connectivity)\n",
    "* [Summary](#Summary)\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Load libraries, settings and data\n",
    "\n",
    "import json\n",
    "import pickle\n",
    "import warnings\n",
    "from collections import Counter\n",
    "import math\n",
    "\n",
    "import contextily as cx\n",
    "import folium\n",
    "import geopandas as gpd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import osmnx as ox\n",
    "import pandas as pd\n",
    "import yaml\n",
    "from matplotlib import cm, colors\n",
    "\n",
    "from bikedna import evaluation_functions as eval_func\n",
    "from bikedna import plotting_functions as plot_func\n",
    "\n",
    "%run ../settings/yaml_variables.py\n",
    "%run ../settings/plotting.py\n",
    "%run ../settings/tiledict.py\n",
    "%run ../settings/load_refdata.py\n",
    "%run ../settings/df_styler.py\n",
    "%run ../settings/paths.py\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "grid = ref_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data completeness\n",
    "\n",
    "### Network density\n",
    "\n",
    "In this setting, network density refers to the length of edges or number of nodes per km2. This is the usual definition of network density in spatial (road) networks, which is distinct from the *structural* network density known more generally in network science. Without comparing to a reference data set, network density does not in itself indicate spatial data quality. For anyone familiar with the study area, network density can however indicate whether parts of the area appear to be under- or over-mapped.\n",
    "\n",
    "**Method**\n",
    "\n",
    "The density here is not based on the geometric length of edges, but instead on the computed length of the infrastructure. For example, a 100-meter-long bidirectional path contributes with 200 meters of bicycle infrastructure. This method is used to take into account different ways of mapping bicycle infrastructure, which otherwise can introduce large deviations in network density. With `compute_network_density`, the number of elements (nodes, dangling nodes, and total infrastructure length) per unit area is calculated. The density is computed twice: first for the study area for both the entire network ('global density'), then for each of the grid cells ('local density'). Both global and local densities are computed for the entire network and for protected and unprotected infrastructure.\n",
    "\n",
    "**Interpretation**\n",
    "\n",
    "Since the analysis conducted here is intrinsic, i.e. it makes no use of external information, it cannot be known whether a low-density value is due to incomplete mapping, or due to actual lack of infrastructure in the area. However, a comparison of the grid cell density values can provide some insights, for example:\n",
    "* lower-than-average infrastructure density indicates a locally sparser network\n",
    "* higher-than-average node density indicates that there are relatively many intersections in a grid cell\n",
    "* higher-than-average dangling node density indicates that there are relatively many dead ends in a grid cell\n",
    "\n",
    "#### Global network density\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entire study area\n",
    "edge_density, node_density, dangling_node_density = eval_func.compute_network_density(\n",
    "    (ref_edges_simplified, ref_nodes_simplified),\n",
    "    grid.unary_union.area,\n",
    "    return_dangling_nodes=True,\n",
    ")\n",
    "\n",
    "density_results = {}\n",
    "density_results[\"edge_density_m_sqkm\"] = edge_density\n",
    "density_results[\"node_density_count_sqkm\"] = node_density\n",
    "density_results[\"dangling_node_density_count_sqkm\"] = dangling_node_density\n",
    "\n",
    "ref_protected = ref_edges_simplified.loc[ref_edges_simplified.protected == \"protected\"]\n",
    "ref_unprotected = ref_edges_simplified.loc[\n",
    "    ref_edges_simplified.protected == \"unprotected\"\n",
    "]\n",
    "ref_mixed = ref_edges_simplified.loc[ref_edges_simplified.protected == \"mixed\"]\n",
    "\n",
    "ref_data = [ref_protected, ref_unprotected, ref_mixed]\n",
    "labels = [\"protected_density\", \"unprotected_density\", \"mixed_density\"]\n",
    "\n",
    "for data, label in zip(ref_data, labels):\n",
    "    if len(data) > 0:\n",
    "        ref_edge_density_type, _ = eval_func.compute_network_density(\n",
    "            (data, ref_nodes_simplified),\n",
    "            grid.unary_union.area,\n",
    "            return_dangling_nodes=False,\n",
    "        )\n",
    "        density_results[label + \"_m_sqkm\"] = ref_edge_density_type\n",
    "    else:\n",
    "        density_results[label + \"_m_sqkm\"] = 0\n",
    "\n",
    "protected_edge_density = density_results[\"protected_density_m_sqkm\"]\n",
    "unprotected_edge_density = density_results[\"unprotected_density_m_sqkm\"]\n",
    "mixed_protection_edge_density = density_results[\"mixed_density_m_sqkm\"]\n",
    "\n",
    "print(f\"For the entire study area, there are:\")\n",
    "print(f\"- {edge_density:.2f} meters of bicycle infrastructure per km2.\")\n",
    "print(f\"- {node_density:.2f} nodes in the bicycle network per km2.\")\n",
    "print(\n",
    "    f\"- {dangling_node_density:.2f} dangling nodes in the bicycle network per km2.\"\n",
    ")\n",
    "print(\n",
    "    f\"- {protected_edge_density:.2f} meters of protected bicycle infrastructure per km2.\"\n",
    ")\n",
    "print(\n",
    "    f\"- {unprotected_edge_density:.2f} meters of unprotected bicycle infrastructure per km2.\"\n",
    ")\n",
    "print(\n",
    "    f\"- {mixed_protection_edge_density:.2f} meters of mixed protection bicycle infrastructure per km2.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save stats to csv\n",
    "pd.DataFrame(\n",
    "    {\n",
    "        \"metric\": [\n",
    "            \"meters of bicycle infrastructure per square km\",\n",
    "            \"nodes in the bicycle network per square km\",\n",
    "            \"dangling nodes in the bicycle network per square km\",\n",
    "            \"meters of protected bicycle infrastructure per square km\",\n",
    "            \"meters of unprotected bicycle infrastructure per square km\",\n",
    "            \"meters of mixed protection bicycle infrastructure per square km\",\n",
    "        ],\n",
    "        \"value\": [\n",
    "            np.round(edge_density, 2),\n",
    "            np.round(node_density, 2),\n",
    "            np.round(dangling_node_density, 2),\n",
    "            np.round(protected_edge_density, 2),\n",
    "            np.round(unprotected_edge_density, 2),\n",
    "            np.round(mixed_protection_edge_density, 2),\n",
    "        ],\n",
    "    }\n",
    ").to_csv(ref_results_data_fp + \"stats_area.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Local network density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Per grid cell\n",
    "results_dict = {}\n",
    "data = (ref_edges_simp_joined, ref_nodes_simp_joined.set_index(\"osmid\"))\n",
    "\n",
    "[\n",
    "    eval_func.run_grid_analysis(\n",
    "        grid_id,\n",
    "        data,\n",
    "        results_dict,\n",
    "        eval_func.compute_network_density,\n",
    "        grid[\"geometry\"].loc[grid.grid_id == grid_id].area.values[0],\n",
    "        return_dangling_nodes=True,\n",
    "    )\n",
    "    for grid_id in grid_ids\n",
    "]\n",
    "\n",
    "results_df = pd.DataFrame.from_dict(results_dict, orient=\"index\")\n",
    "results_df.reset_index(inplace=True)\n",
    "results_df.rename(\n",
    "    columns={\n",
    "        \"index\": \"grid_id\",\n",
    "        0: \"ref_edge_density\",\n",
    "        1: \"ref_node_density\",\n",
    "        2: \"ref_dangling_node_density\",\n",
    "    },\n",
    "    inplace=True,\n",
    ")\n",
    "\n",
    "grid = eval_func.merge_results(grid, results_df, \"left\")\n",
    "\n",
    "ref_protected = ref_edges_simp_joined.loc[\n",
    "    ref_edges_simp_joined.protected == \"protected\"\n",
    "]\n",
    "ref_unprotected = ref_edges_simp_joined.loc[\n",
    "    ref_edges_simp_joined.protected == \"unprotected\"\n",
    "]\n",
    "ref_mixed = ref_edges_simp_joined.loc[ref_edges_simp_joined.protected == \"mixed\"]\n",
    "\n",
    "ref_data = [ref_protected, ref_unprotected, ref_mixed]\n",
    "\n",
    "ref_labels = [\"ref_protected_density\", \"ref_unprotected_density\", \"ref_mixed_density\"]\n",
    "\n",
    "for data, label in zip(ref_data, ref_labels):\n",
    "    if len(data) > 0:\n",
    "        results_dict = {}\n",
    "        data = (ref_edges_simp_joined.loc[data.index], ref_nodes_simp_joined)\n",
    "        [\n",
    "            eval_func.run_grid_analysis(\n",
    "                grid_id,\n",
    "                data,\n",
    "                results_dict,\n",
    "                eval_func.compute_network_density,\n",
    "                grid[\"geometry\"].loc[grid.grid_id == grid_id].area.values[0],\n",
    "            )\n",
    "            for grid_id in grid_ids\n",
    "        ]\n",
    "\n",
    "        results_df = pd.DataFrame.from_dict(results_dict, orient=\"index\")\n",
    "        results_df.reset_index(inplace=True)\n",
    "        results_df.rename(columns={\"index\": \"grid_id\", 0: label}, inplace=True)\n",
    "        results_df.drop(1, axis=1, inplace=True)\n",
    "\n",
    "        grid = eval_func.merge_results(grid, results_df, \"left\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Network density grid plots\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "plot_cols = [\"ref_edge_density\", \"ref_node_density\", \"ref_dangling_node_density\"]\n",
    "plot_titles = [\n",
    "    area_name + f\": {reference_name} edge density\",\n",
    "    area_name + f\": {reference_name} node density\",\n",
    "    area_name + f\": {reference_name} dangling node density\",\n",
    "]\n",
    "filepaths = [\n",
    "    ref_results_static_maps_fp + \"density_edge_reference\",\n",
    "    ref_results_static_maps_fp + \"density_node_reference\",\n",
    "    ref_results_static_maps_fp + \"density_danglingnode_reference\",\n",
    "]\n",
    "cmaps = [pdict[\"pos\"], pdict[\"pos\"], pdict[\"pos\"]]\n",
    "no_data_cols = [\"count_ref_edges\", \"count_ref_nodes\", \"count_ref_nodes\"]\n",
    "\n",
    "plot_func.plot_grid_results(\n",
    "    grid=grid,\n",
    "    plot_cols=plot_cols,\n",
    "    plot_titles=plot_titles,\n",
    "    filepaths=filepaths,\n",
    "    cmaps=cmaps,\n",
    "    alpha=pdict[\"alpha_grid\"],\n",
    "    cx_tile=cx_tile_2,\n",
    "    no_data_cols=no_data_cols,\n",
    "    attr=reference_name\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Densities of protected and unprotected infrastructure\n",
    "\n",
    "In BikeDNA, *protected infrastructure* refers to all bicycle infrastructure which is either separated from car traffic by for example an elevated curb, bollards, or other physical barriers, or for cycle tracks that are not adjacent to a street.\n",
    "\n",
    "*Unprotected infrastructure* are all other types of lanes that are dedicated for bicyclists, but which only are separated by car traffic by e.g., a painted line on the street."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Network density grid plots\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "plot_cols = [\"ref_protected_density\", \"ref_unprotected_density\", \"ref_mixed_density\"]\n",
    "\n",
    "plot_cols = [p for p in plot_cols if p in grid.columns]\n",
    "\n",
    "plot_titles = [\n",
    "    area_name + f\": {reference_name} protected infrastructure density (m/km2)\",\n",
    "    area_name + f\": {reference_name} unprotected infrastructure density (m/km2)\",\n",
    "    area_name + f\": {reference_name} mixed protection infrastructure density (m/km2)\"\n",
    "]\n",
    "\n",
    "# plot_titles = plot_titles[0:len(plot_cols)]\n",
    "\n",
    "filepaths = [\n",
    "    ref_results_static_maps_fp + \"density_protected_reference\",\n",
    "    ref_results_static_maps_fp + \"density_unprotected_reference\",\n",
    "    ref_results_static_maps_fp + \"density_mixed_reference\",\n",
    "]\n",
    "\n",
    "cmaps = [pdict[\"pos\"]] * len(plot_cols)\n",
    "no_data_cols = [\"ref_protected_density\", \"ref_unprotected_density\", \"ref_mixed_density\"]\n",
    "norm_min = [0] * len(plot_cols)\n",
    "norm_max = [max(grid[plot_cols].fillna(value=0).max())] * len(plot_cols)\n",
    "\n",
    "plot_func.plot_grid_results(\n",
    "    grid=grid,\n",
    "    plot_cols=plot_cols,\n",
    "    plot_titles=plot_titles,\n",
    "    filepaths=filepaths,\n",
    "    cmaps=cmaps,\n",
    "    alpha=pdict[\"alpha_grid\"],\n",
    "    cx_tile=cx_tile_2,\n",
    "    no_data_cols=no_data_cols,\n",
    "    use_norm=True,\n",
    "    norm_min=norm_min,\n",
    "    norm_max=norm_max,\n",
    "    attr=reference_name\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Network topology\n",
    "\n",
    "This section explores the geometric and topological features of the data. These are, for example, network density, disconnected components, and dangling (degree one) nodes. It also includes exploring whether there are nodes that are very close to each other but do not share an edge - a potential sign of edge undershoots - or if there are intersecting edges without a node at the intersection, which might indicate a digitizing error that will distort routing on the network.\n",
    "\n",
    "Due to the fragmented nature of most bicycle networks, many metrics, such as missing links or network gaps, can simply reflect the true extent of the infrastructure ([Natera Orozco et al., 2020](https://onlinelibrary.wiley.com/doi/pdfdirect/10.1111/gean.12324)). This is different for road networks, where e.g., disconnected components could more readily be interpreted as a data quality issue. Therefore, the analysis only takes very small network gaps into account as potential data quality issues.\n",
    "\n",
    "### Simplification outcome\n",
    "\n",
    "To compare the structure and true ratio between nodes and edges in the network, a simplified network representation which only includes nodes at endpoints and intersections was created in notebook `1b` by removing all interstitial nodes.\n",
    "\n",
    "Comparing the degree distribution for the networks before and after simplification is a quick sanity check for the simplification routine. Typically, the vast majority of nodes in the non-simplified network will be of degree two; in the simplified network, however, most nodes will have degrees other than two. Degree two nodes are retained in only two cases: if they represent a connection point between two different types of infrastructure; or if they are needed in order to avoid self-loops (edges whose start and end points are identical) or multiple edges between the same pair of nodes. \n",
    "\n",
    "<p align=\"center\">\n",
    "<img src='../../images/network_simplification_illustration.png' width=300/>\n",
    "\n",
    "*Non-simplified network (left) and simplified network (right)*.\n",
    "\n",
    "</p>\n",
    "\n",
    "**Method**\n",
    "\n",
    "The node degree distributions before and after simplification are plotted below. \n",
    "\n",
    "**Interpretation**\n",
    "\n",
    "Typically, the node degree distribution will go from high (before simplification) to low (after simplification) counts of degree two nodes, while it will not change for all other degrees (1, or 3 and higher). Further, the total number of nodes will see a strong decline. If the simplified graph still maintains a relatively high number of degree two nodes, or if the number of nodes with other degrees changes after the simplification, this might point to issues either with the graph conversion or with the simplification process.\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decrease in network elements after simplification\n",
    "\n",
    "edge_percent_diff = (len(ref_edges) - len(ref_edges_simplified)) / len(ref_edges) * 100\n",
    "node_percent_diff = (len(ref_nodes) - len(ref_nodes_simplified)) / len(ref_nodes) * 100\n",
    "\n",
    "simplification_results = {\n",
    "    \"edge_percent_diff\": edge_percent_diff,\n",
    "    \"node_percent_diff\": node_percent_diff,\n",
    "}\n",
    "\n",
    "print(\n",
    "    f\"Simplifying the network decreased the number of edges with {edge_percent_diff:.1f}% and the number of nodes with {node_percent_diff:.1f}%.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Node degree distribution\n",
    "\n",
    "set_renderer(renderer_plot)\n",
    "fig, ax = plt.subplots(1, 2, figsize=pdict[\"fsbar_short\"], sharey=True)\n",
    "\n",
    "degree_sequence_before = sorted((d for n, d in ref_graph.degree()), reverse=True)\n",
    "degree_sequence_after = sorted(\n",
    "    (d for n, d in ref_graph_simplified.degree()), reverse=True\n",
    ")\n",
    "\n",
    "\n",
    "ax[0].bar(*np.unique(degree_sequence_before, return_counts=True), tick_label = np.unique(degree_sequence_before), color=pdict[\"ref_base\"])\n",
    "ax[0].set_title(\"Before simplification\")\n",
    "ax[0].set_xlabel(\"Degree\")\n",
    "ax[0].set_ylabel(\"Nodes\")\n",
    "\n",
    "ax[1].bar(*np.unique(degree_sequence_after, return_counts=True), tick_label = np.unique(degree_sequence_after), color=pdict[\"ref_base\"])\n",
    "ax[1].set_title(\"After simplification\")\n",
    "ax[1].set_xlabel(\"Degree\")\n",
    "\n",
    "plt.suptitle(f\"{area_name}: {reference_name} degree distributions\")\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_plots_fp + \"degree_dist_reference\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dangling nodes\n",
    "\n",
    "Dangling nodes are nodes of degree one, i.e. they have only one single edge attached to them. Most networks will naturally contain a number of dangling nodes. Dangling nodes can occur at actual dead-ends (representing a cul-de-sac) or at the endpoints of certain features, e.g. when a bicycle path ends in the middle of a street. However, dangling nodes can also occur as a data quality issue in case of over/undershoots (see next section). The number of dangling nodes in a network does to some extent also depend on the digitization method, as shown in the illustration below. \n",
    "\n",
    "Therefore, the presence of dangling nodes is in itself not a sign of low data quality. However, a high number of dangling nodes in an area that is not known for containing many dead-ends can indicate digitization errors and problems with edge over/undershoots.\n",
    "\n",
    "<!-- <table><tr><td><img src='../../images/dangling_nodes_illustration.png' width=300 /></td><td><img src='../../images/no_dangling_nodes_illustration.png' width=295 /></td></tr></table> -->\n",
    "\n",
    "<p align=\"center\">\n",
    "\n",
    "<img src='../../images/dangling_nodes_illustration_new.png' width=350/>\n",
    "\n",
    "*Left: Dangling nodes occur where road features end. Right: However, when separate features are joined at the end, there will be no dangling nodes.*\n",
    "\n",
    "</p>\n",
    "\n",
    "**Method**\n",
    "\n",
    "Below, a list of all dangling nodes is obtained with the help of `get_dangling_nodes`. Then, the network with all its nodes is plotted. The dangling nodes are shown in color, all other nodes are shown in black.\n",
    "\n",
    "**Interpretation**\n",
    "\n",
    "We recommend a visual analysis in order to interpret the spatial distribution of dangling nodes, with particular attention to areas of high dangling node density. It is important to understand where dangling nodes come from: are they actual dead-ends or digitization errors (e.g., over/undershoots)? A higher number of digitization errors points to lower data quality.\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute number of dangling nodes\n",
    "dangling_nodes = eval_func.get_dangling_nodes(\n",
    "    ref_edges_simplified, ref_nodes_simplified\n",
    ")\n",
    "\n",
    "# Export results\n",
    "dangling_nodes.to_file(ref_results_data_fp + \"dangling_nodes.gpkg\", index=False)\n",
    "\n",
    "# Compute local count and pct of dangling nodes\n",
    "dn_ref_joined = gpd.overlay(\n",
    "    dangling_nodes, grid[[\"geometry\", \"grid_id\"]], how=\"intersection\"\n",
    ")\n",
    "df = eval_func.count_features_in_grid(dn_ref_joined, \"ref_dangling_nodes\")\n",
    "grid = eval_func.merge_results(grid, df, \"left\")\n",
    "\n",
    "grid[\"ref_dangling_nodes_pct\"] = np.round(\n",
    "    100 * grid.count_ref_dangling_nodes / grid.count_ref_simplified_nodes, 2\n",
    ")\n",
    "\n",
    "# set to zero where there are simplified nodes but no dangling nodes\n",
    "grid[\"ref_dangling_nodes_pct\"].loc[\n",
    "    grid.count_ref_simplified_nodes.notnull() & grid.ref_dangling_nodes_pct.isnull()\n",
    "] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot dangling nodes\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "fig, ax = plt.subplots(1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"3.5%\", pad=\"1%\")\n",
    "    \n",
    "grid.plot(\n",
    "    cax=cax,\n",
    "    column=\"ref_dangling_nodes_pct\",\n",
    "    ax=ax,\n",
    "    alpha=pdict[\"alpha_grid\"],\n",
    "    cmap=pdict[\"pos\"],\n",
    "    legend=True,\n",
    ")\n",
    "\n",
    "# add no data patches\n",
    "grid[grid[\"count_ref_simplified_nodes\"].isnull()].plot(\n",
    "    cax=cax,\n",
    "    ax=ax,\n",
    "    facecolor=pdict[\"nodata_face\"],\n",
    "    edgecolor=pdict[\"nodata_edge\"],\n",
    "    linewidth= pdict[\"line_nodata\"],\n",
    "    hatch=pdict[\"nodata_hatch\"],\n",
    "    alpha= pdict[\"alpha_nodata\"],\n",
    ")\n",
    "\n",
    "ax.legend(handles=[nodata_patch], loc=\"upper right\")\n",
    "ax.set_title(f\"{area_name}: {reference_name} percent of dangling nodes\")\n",
    "ax.set_axis_off()\n",
    "cx.add_basemap(ax=ax, crs=study_crs, source=cx_tile_2)\n",
    "cx.add_attribution(ax=ax, text=f\"(C) {reference_name}\")\n",
    "txt = ax.texts[-1]\n",
    "txt.set_position([1,0.00])\n",
    "txt.set_ha('right')\n",
    "txt.set_va('bottom')\n",
    "\n",
    "plot_func.save_fig(fig,ref_results_static_maps_fp + \"pct_dangling_nodes_reference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "interactive"
    ]
   },
   "outputs": [],
   "source": [
    "# Interactive plot of dangling nodes\n",
    "\n",
    "edges_simplified_folium = plot_func.make_edgefeaturegroup(\n",
    "    gdf=ref_edges_simplified,\n",
    "    mycolor=pdict[\"base\"],\n",
    "    myweight=pdict[\"line_base\"],\n",
    "    nametag=\"Edges\",\n",
    "    show_edges=True,\n",
    ")\n",
    "\n",
    "nodes_simplified_folium = plot_func.make_nodefeaturegroup(\n",
    "    gdf=ref_nodes_simplified,\n",
    "    mysize=pdict[\"mark_base\"],\n",
    "    mycolor=pdict[\"base\"],\n",
    "    nametag=\"All nodes\",\n",
    "    show_nodes=True,\n",
    ")\n",
    "\n",
    "dangling_nodes_folium = plot_func.make_nodefeaturegroup(\n",
    "    gdf=dangling_nodes,\n",
    "    mysize=pdict[\"mark_emp\"],\n",
    "    mycolor=pdict[\"ref_base\"],\n",
    "    nametag=\"Dangling nodes\",\n",
    "    show_nodes=True,\n",
    ")\n",
    "\n",
    "m = plot_func.make_foliumplot(\n",
    "    feature_groups=[\n",
    "        edges_simplified_folium,\n",
    "        nodes_simplified_folium,\n",
    "        dangling_nodes_folium,\n",
    "    ],\n",
    "    layers_dict=folium_layers,\n",
    "    center_gdf=ref_nodes_simplified,\n",
    "    center_crs=ref_nodes_simplified.crs,\n",
    "    attr=reference_name\n",
    ")\n",
    "\n",
    "bounds = plot_func.compute_folium_bounds(ref_nodes_simplified)\n",
    "m.fit_bounds(bounds)\n",
    "\n",
    "\n",
    "m.save(ref_results_inter_maps_fp + \"folium_danglingmap_reference.html\")\n",
    "\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Interactive map saved at \" + ref_results_inter_maps_fp.lstrip(\"../\") + \"folium_danglingmap_reference.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under/overshoots \n",
    "\n",
    "When two nodes in a simplified network are placed within a distance of a few meters, but do not share a common edge, it is often due to an edge over/undershoot or another digitizing error. An undershoot occurs when two features are supposed to meet, but instead are just in close proximity to each other. An overshoot occurs when two features meet and one of them extends beyond the other. See the image below for an illustration. For a more detailed explanation of over/undershoots, see the [GIS Lounge website](https://www.gislounge.com/digitizing-errors-in-gis/).\n",
    "\n",
    "<!-- <table>\n",
    "<tr><td><img src='../../images/undershoot_illustration.png' width=300 /></td><td><img src='../../images/overshoot_illustration.png' width=300  /></td></tr>\n",
    "</table>\n",
    "\n",
    "*Left: Undershoots happen when two line features are not properly joined, for example at intersection. Right: Overshoots refer to situations where a line feature extends too far beyond at intersecting line, rather than ending at the intersection.*  -->\n",
    "\n",
    "<p align=\"center\">\n",
    "\n",
    "<img src='../../images/over_undershoots2.png' width=350/>\n",
    "\n",
    "*Left: Undershoots happen when two line features are not properly joined, for example at an intersection. Right: Overshoots refer to situations where a line feature extends too far beyond at intersecting line, rather than ending at the intersection.* \n",
    "\n",
    "</p>\n",
    "\n",
    "**Method**\n",
    "\n",
    "*Undershoots:* First, the `length_tolerance` (in meters) is defined in the cell below. Then, with `find_undershoots`, all pairs of dangling nodes that have a maximum of `length_tolerance` distance between them, are identified as undershoots, and the results are plotted.\n",
    "\n",
    "*Overshoots:* First, the `length_tolerance` (in meters) is defined in the cell below. Then, with `find_overshoots`, all network edges that have a dangling node attached to them and that have a maximum length of `length_tolerance` are identifed as overshoots, and the results are plotted.\n",
    "\n",
    "The method for over/undershoot detection is inspired by [Neis et al. (2012)](https://www.mdpi.com/1999-5903/4/1/1).\n",
    "\n",
    "**Interpretation**\n",
    "\n",
    "Under/overshoots are not necessarily always a data quality issue - they might be instead an accurate representation of the network conditions or of the digitization strategy. For example, a cycle path might end abruptly soon after a turn, which results in an overshoot. Protected cycle paths are sometimes digitized as interrupted at intersections which results in intersection undershoots.\n",
    "\n",
    "The interpretation of the impact of over/undershoots on data quality is context dependent. For certain applications, such as routing, overshoots do not present a particular challenge; they can, however, pose an issue for other applications such as network analysis, given that they skew the network structure.  Undershoots, on the contrary, are a serious problem for routing applications, especially if only bicycle infrastructure is considered. They also pose a problem for network analysis, for example for any path-based metric, such as most centrality measures like betweenness centrality.\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>User configurations</b>\n",
    "<br>\n",
    "<br>\n",
    "In the analysis of over and undershoots, the user can modify the length tolerance for both over and undershoots. <br>\n",
    "For example, a length tolerance of 3 meters for overshoots means that only edge snippets with a length of 3 meters or less are considered overshoots.<br>\n",
    "A tolerance of 5 meters for undershoots means that only gaps of 5 meters or less are considered undershoots.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# USER INPUT: LENGTH TOLERANCE FOR OVER- AND UNDERSHOOTS\n",
    "length_tolerance_over = 3\n",
    "length_tolerance_under = 3\n",
    "\n",
    "for s in [length_tolerance_over, length_tolerance_under]:\n",
    "    assert isinstance(s, int) or isinstance(s, float), print(\n",
    "        \"Settings must be integer or float values!\"\n",
    "    )\n",
    "\n",
    "print(f\"Running overshoot analysis with a tolerance threshold of {length_tolerance_over} m.\")\n",
    "print(f\"Running undershoot analysis with a tolerance threshold of {length_tolerance_under} m.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Overshoots\n",
    "overshoots = eval_func.find_overshoots(\n",
    "    dangling_nodes,\n",
    "    ref_edges_simplified,\n",
    "    length_tolerance_over,\n",
    "    return_overshoot_edges=True,\n",
    ")\n",
    "\n",
    "print(\n",
    "    f\"{len(overshoots)} potential overshoots were identified with a length tolerance of {length_tolerance_over} m.\"\n",
    ")\n",
    "\n",
    "### Undershoots\n",
    "undershoot_dict, undershoot_nodes = eval_func.find_undershoots(\n",
    "    dangling_nodes,\n",
    "    ref_edges_simplified,\n",
    "    length_tolerance_under,\n",
    "    \"edge_id\",\n",
    "    return_undershoot_nodes=True,\n",
    ")\n",
    "\n",
    "print(\n",
    "    f\"{len(undershoot_nodes)} potential undershoots were identified with a length tolerance of {length_tolerance_under} m.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to csv\n",
    "overshoots[[\"edge_id\", \"length\"]].to_csv(\n",
    "    ref_results_data_fp + f\"overshoot_edges_{length_tolerance_over}.csv\", header = [\"edge_id\", \"length (m)\"], index = False\n",
    ")\n",
    "\n",
    "\n",
    "pd.DataFrame(undershoot_nodes[\"nodeID\"].to_list(), columns=[\"node_id\"]).to_csv(\n",
    "    ref_results_data_fp + f\"undershoot_nodes_{length_tolerance_under}.csv\", index=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "interactive"
    ]
   },
   "outputs": [],
   "source": [
    "# Interactive plot over/undershoots\n",
    "\n",
    "simplified_edges_folium = plot_func.make_edgefeaturegroup(\n",
    "    gdf=ref_edges_simplified,\n",
    "    mycolor=pdict[\"base\"],\n",
    "    myweight=pdict[\"line_base\"],\n",
    "    nametag=\"Edges\",\n",
    "    show_edges=True,\n",
    ")\n",
    "\n",
    "fg = [simplified_edges_folium]\n",
    "\n",
    "if len(overshoots) > 0 or len(undershoot_nodes) > 0:\n",
    "\n",
    "    if len(overshoots) > 0:\n",
    "\n",
    "        overshoots_folium = plot_func.make_edgefeaturegroup(\n",
    "            gdf=overshoots,\n",
    "            mycolor=pdict[\"ref_emp2\"],\n",
    "            myweight=pdict[\"line_emp2\"],\n",
    "            nametag=\"Overshoots\",\n",
    "            show_edges=True,\n",
    "        )\n",
    "\n",
    "        fg.append(overshoots_folium)\n",
    "\n",
    "    if len(undershoot_nodes) > 0:\n",
    "\n",
    "        undershoot_nodes_folium = plot_func.make_nodefeaturegroup(\n",
    "            gdf=undershoot_nodes,\n",
    "            mysize=pdict[\"mark_emp\"],\n",
    "            mycolor=pdict[\"ref_contrast\"],\n",
    "            nametag=\"Undershoot nodes\",\n",
    "            show_nodes=True,\n",
    "        )\n",
    "\n",
    "        fg.append(undershoot_nodes_folium)\n",
    "\n",
    "    m = plot_func.make_foliumplot(\n",
    "        feature_groups=fg,\n",
    "        layers_dict=folium_layers,\n",
    "        center_gdf=ref_nodes_simplified,\n",
    "        center_crs=ref_nodes_simplified.crs,\n",
    "    )\n",
    "\n",
    "    bounds = plot_func.compute_folium_bounds(ref_nodes_simplified)\n",
    "    m.fit_bounds(bounds)\n",
    "\n",
    "    m.save(\n",
    "        ref_results_inter_maps_fp\n",
    "        + f\"underovershoots_{length_tolerance_under}_{length_tolerance_over}_reference.html\"\n",
    "    )\n",
    "\n",
    "    display(m)\n",
    "\n",
    "if len(undershoot_nodes) == 0:\n",
    "    print(\"There are no undershoots to plot.\")\n",
    "if len(overshoots) == 0:\n",
    "    print(\"There are no overshoots to plot.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if len(overshoots) > 0 or len(undershoot_nodes) > 0:\n",
    "    print(\"Interactive map saved at \" + ref_results_inter_maps_fp.lstrip(\"../\") + f\"overundershoots_{length_tolerance_under}_{length_tolerance_over}_reference.html\")\n",
    "else:\n",
    "    print(\"There are no under/overshoots to plot.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Network components\n",
    "\n",
    "Disconnected components do not share any elements (nodes/edges). In other words, there is no network path that could lead from one disconnected component to the other. As mentioned above, most real-world networks of bicycle infrastructure do consist of many disconnected components ([Natera Orozco et al., 2020](https://onlinelibrary.wiley.com/doi/pdfdirect/10.1111/gean.12324)). However, when two disconnected components are very close to each other, it might be a sign of a missing edge or another digitizing error.\n",
    "\n",
    "**Method**\n",
    "\n",
    "First, with the help of `return_components`, a list of all (disconnected) components of the network is obtained. The total number of components is printed and all components are plotted in different colors for visual analysis. Next, the component size distribution (with components ordered by the network length they contain) is plotted, followed by a plot of the largest connected component. \n",
    "\n",
    "**Interpretation**\n",
    "\n",
    "As with many of the previous analysis steps, knowledge of the area is crucial for a correct interpretation of component analysis. Given that the data represents the actual infrastructure accurately, bigger components indicate coherent network parts, while smaller components indicate scattered infrastructure (e.g., one single bicycle path along a street that does not connect to any other bicycle infrastructure). A high number of disconnected components in near vicinity of each other indicates digitization errors or missing data. \n",
    "\n",
    "### Disconnected components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_components = eval_func.return_components(ref_graph_simplified)\n",
    "print(\n",
    "    f\"The network in the study area has {len(ref_components)} disconnected components.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot disconnected components\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "\n",
    "# set seed for colors\n",
    "np.random.seed(42)\n",
    "\n",
    "# generate enough random colors to plot all components\n",
    "randcols = np.random.rand(len(ref_components), 3)\n",
    "randcols[0, :] = col_to_rgb(pdict['ref_base'])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "ax.set_title(f\"{area_name}: {reference_name} disconnected components\")\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "for j, c in enumerate(ref_components):\n",
    "    if len(c.edges) > 0:\n",
    "        edges = ox.graph_to_gdfs(c, nodes=False)\n",
    "        edges.plot(ax=ax, color=randcols[j])\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=study_crs, source=cx_tile_2)\n",
    "cx.add_attribution(ax=ax, text=f\"(C) {reference_name}\")\n",
    "txt = ax.texts[-1]\n",
    "txt.set_position([1,0.00])\n",
    "txt.set_ha('right')\n",
    "txt.set_va('bottom')\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + \"all_components_reference\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Components per grid cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign component ids to grid\n",
    "grid = eval_func.assign_component_id_to_grid(\n",
    "    ref_edges_simplified,\n",
    "    ref_edges_simp_joined,\n",
    "    ref_components,\n",
    "    grid,\n",
    "    prefix=\"ref\",\n",
    "    edge_id_col=\"edge_id\",\n",
    ")\n",
    "\n",
    "fill_na_dict = {\"component_ids_ref\": \"\"}\n",
    "grid.fillna(value=fill_na_dict, inplace=True)\n",
    "\n",
    "grid[\"component_count_ref\"] = grid.component_ids_ref.apply(lambda x: len(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot number of components per grid cell\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "ncolors = grid[\"component_count_ref\"].max()\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"3.5%\", pad=\"1%\")\n",
    "\n",
    "mycm = cm.get_cmap(pdict[\"seq\"], ncolors)\n",
    "grid[grid.component_count_ref>0].plot(\n",
    "    cax=cax,\n",
    "    ax=ax,\n",
    "    column=\"component_count_ref\",\n",
    "    legend=True,\n",
    "    legend_kwds={'ticks': list(range(1, ncolors+1))},\n",
    "    cmap=mycm,\n",
    "    alpha=pdict[\"alpha_grid\"],\n",
    ")\n",
    "\n",
    "# add no data patches\n",
    "grid[grid[\"count_ref_edges\"].isnull()].plot(\n",
    "    cax=cax,\n",
    "    ax=ax,\n",
    "    facecolor=pdict[\"nodata_face\"],\n",
    "    edgecolor=pdict[\"nodata_edge\"],\n",
    "    linewidth= pdict[\"line_nodata\"],\n",
    "    hatch=pdict[\"nodata_hatch\"],\n",
    "    alpha=pdict[\"alpha_nodata\"],\n",
    ")\n",
    "\n",
    "ax.legend(handles=[nodata_patch], loc=\"upper right\")\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=study_crs, source=cx_tile_2)\n",
    "cx.add_attribution(ax=ax, text=f\"(C) {reference_name}\")\n",
    "txt = ax.texts[-1]\n",
    "txt.set_position([1,0.00])\n",
    "txt.set_ha('right')\n",
    "txt.set_va('bottom')\n",
    "ax.set_title(area_name + f\": {reference_name} number of components in grid cells\") \n",
    "ax.set_axis_off()\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + f\"number_of_components_in_grid_cells_reference\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Component length distribution\n",
    "\n",
    "The distribution of all network component lengths can be visualized in a so-called *Zipf plot*, which orders the lengths of each component by rank, showing the largest component's length on the left, then the second largest component's length, etc., until the smallest component's length on the right. When a Zipf plot follows a straight line in [log-log scale](https://en.wikipedia.org/wiki/Logarithmic_scale), it means that there is a much higher chance to find small disconnected components than expected from traditional distributions [(Clauset et al., 2009)]( https://epubs.siam.org/doi/abs/10.1137/070710111). This can mean that there has been no consolidation of the network, only piece-wise or random additions [(Szell et al., 2022)](https://www.nature.com/articles/s41598-022-10783-y), or that the data itself suffers from many gaps and topology errors resulting in small disconnected components.\n",
    "\n",
    "However, it can also happen that the largest connected component (the leftmost marker in the plot at rank $10^0$) is a clear outlier, while the rest of the plot follows a different shape. This can mean that at the infrastructure level, most of the infrastructure has been connected to one large component, and that the data reflects this - i.e. the data is not suffering from gaps and missing links to a large extent.\n",
    "\n",
    "Bicycle networks might also be somewhere inbetween, with several large components as outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Zipf plot of component lengths\n",
    "\n",
    "set_renderer(renderer_plot)\n",
    "\n",
    "components_length = {}\n",
    "\n",
    "for i, c in enumerate(ref_components):\n",
    "    c_length = 0\n",
    "    for (u, v, l) in c.edges(data=\"length\"):\n",
    "        c_length += l\n",
    "    components_length[i] = c_length\n",
    "\n",
    "components_df = pd.DataFrame.from_dict(components_length, orient=\"index\")\n",
    "components_df.rename(columns={0: \"component_length\"}, inplace=True)\n",
    "\n",
    "fig = plt.figure(figsize=pdict[\"fsbar_small\"])\n",
    "axes = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "axes.set_axisbelow(True)\n",
    "axes.grid(True,which=\"major\",ls=\"dotted\")\n",
    "yvals = sorted(list(components_df[\"component_length\"] / 1000), reverse = True)\n",
    "axes.scatter(\n",
    "    x=[i+1 for i in range(len(components_df))],\n",
    "    y=yvals,\n",
    "    s=18,\n",
    "    color=pdict[\"ref_base\"],\n",
    ")\n",
    "axes.set_ylim(ymin=10**math.floor(math.log10(min(yvals))), ymax=10**math.ceil(math.log10(max(yvals))))\n",
    "axes.set_xscale(\"log\")\n",
    "axes.set_yscale(\"log\")\n",
    "\n",
    "axes.set_ylabel(\"Component length [km]\")\n",
    "axes.set_xlabel(\"Component rank (largest to smallest)\")\n",
    "axes.set_title(area_name+\": \" + f\"{reference_name} component length distribution\")\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_plots_fp + \"component_length_distribution_reference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "largest_cc = max(ref_components, key=len)\n",
    "\n",
    "largest_cc_length = 0\n",
    "\n",
    "for (u, v, l) in largest_cc.edges(data=\"length\"):\n",
    "\n",
    "    largest_cc_length += l\n",
    "\n",
    "largest_cc_pct = largest_cc_length / components_df[\"component_length\"].sum() * 100\n",
    "\n",
    "print(\n",
    "    f\"The largest connected component contains {largest_cc_pct:.2f}% of the network length.\"\n",
    ")\n",
    "\n",
    "lcc_edges = ox.graph_to_gdfs(\n",
    "    G=largest_cc, nodes=False, edges=True, node_geometry=False, fill_edge_geometry=False\n",
    ")\n",
    "\n",
    "# Export to GPKG\n",
    "lcc_edges[[\"edge_id\", \"geometry\"]].to_file(\n",
    "    ref_results_data_fp + \"largest_connected_component.gpkg\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Largest connected component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot of largest connected component\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "fig, ax = plt.subplots(1, 1, figsize=pdict[\"fsmap\"])\n",
    "ref_edges_simplified.plot(ax=ax, color = pdict[\"base\"], linewidth = 1.5, label = \"All smaller components\")\n",
    "lcc_edges.plot(ax=ax, color=pdict[\"ref_base\"], linewidth = 2, label = \"Largest connected component\")\n",
    "grid.plot(ax=ax,alpha=0)\n",
    "ax.set_axis_off()\n",
    "ax.set_title(area_name + f\": {reference_name} largest connected component\")\n",
    "ax.legend()\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=study_crs, source=cx_tile_2)\n",
    "cx.add_attribution(ax=ax, text=f\"(C) {reference_name}\")\n",
    "txt = ax.texts[-1]\n",
    "txt.set_position([1,0.00])\n",
    "txt.set_ha('right')\n",
    "txt.set_va('bottom')\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + f\"largest_conn_comp_reference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Save plot without basemap for potential report titlepage\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "fig, ax = plt.subplots(1, 1, figsize=pdict[\"fsmap\"])\n",
    "ref_edges_simplified.plot(ax=ax, color = pdict[\"base\"], linewidth = 1.5, label = \"Disconnected components\")\n",
    "lcc_edges.plot(ax=ax, color=pdict[\"ref_base\"], linewidth = 2, label = \"Largest connected component\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + \"titleimage\",plot_res=\"high\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing links\n",
    "\n",
    "In the plot of potential missing links between components, all edges that are within the specified distance of an edge on another component are plotted. The gaps between disconnected edges are highlighted with a marker.\n",
    "The map thus highlights edges which, despite being in close proximity of each other, are disconnected and where it thus would not be possible to bike on cycling infrastructure between the edges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>User configurations</b>\n",
    "<br>\n",
    "<br>\n",
    "In the analysis of potential missing links between components, the user must define the threshold for when the distance between two components is considered to be low enough that a digitization error is suspected.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# DEFINE MAX BUFFER DISTANCE BETWEEN COMPONENTS CONSIDERED A GAP/MISSING LINK\n",
    "component_min_distance = 10\n",
    "\n",
    "assert isinstance(component_min_distance, int) or isinstance(\n",
    "    component_min_distance, float\n",
    "), print(\"Setting must be integer or float value!\")\n",
    "\n",
    "print(f\"Running analysis with component distance threshold of {component_min_distance} meters.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "component_gaps = eval_func.find_adjacent_components(\n",
    "    components=ref_components,\n",
    "    buffer_dist=component_min_distance,\n",
    "    crs=study_crs,\n",
    "    edge_id=\"edge_id\",\n",
    ")\n",
    "component_gaps_gdf = gpd.GeoDataFrame.from_dict(\n",
    "    component_gaps, orient=\"index\", geometry=\"geometry\", crs=study_crs\n",
    ")\n",
    "\n",
    "edge_ids = set(\n",
    "    component_gaps_gdf[\"edge_id\" + \"_left\"].to_list()\n",
    "    + component_gaps_gdf[\"edge_id\" + \"_right\"].to_list()\n",
    ")\n",
    "\n",
    "edge_ids = [int(i) for i in edge_ids]\n",
    "edges_with_gaps = ref_edges_simplified.loc[ref_edges_simplified.edge_id.isin(edge_ids)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to csv\n",
    "pd.DataFrame(edge_ids, columns=[\"edge_id\"]).to_csv(\n",
    "    ref_results_data_fp + f\"component_gaps_edges_{component_min_distance}.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "component_gaps_gdf.to_file(\n",
    "    ref_results_data_fp + f\"component_gaps_centroids_{component_min_distance}.gpkg\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "interactive"
    ]
   },
   "outputs": [],
   "source": [
    "# Interactive plot of adjacent components\n",
    "\n",
    "if len(component_gaps) > 0:\n",
    "\n",
    "    simplified_edges_folium = plot_func.make_edgefeaturegroup(\n",
    "        gdf=ref_edges_simplified,\n",
    "        mycolor=pdict[\"ref_base\"],\n",
    "        myweight=pdict[\"line_base\"],\n",
    "        nametag=\"All edges\",\n",
    "        show_edges=True,\n",
    "    )\n",
    "\n",
    "    component_issues_edges_folium = plot_func.make_edgefeaturegroup(\n",
    "        gdf=edges_with_gaps,\n",
    "        mycolor=pdict[\"ref_emp\"],\n",
    "        myweight=pdict[\"line_emp\"],\n",
    "        nametag=\"Adjacent disconnected edges\",\n",
    "        show_edges=True,\n",
    "    )\n",
    "\n",
    "    component_issues_gaps_folium = plot_func.make_markerfeaturegroup(\n",
    "        gdf=component_gaps_gdf, nametag=\"Component gaps\", show_markers=True\n",
    "    )\n",
    "\n",
    "    m = plot_func.make_foliumplot(\n",
    "        feature_groups=[\n",
    "            simplified_edges_folium,\n",
    "            component_issues_edges_folium,\n",
    "            component_issues_gaps_folium,\n",
    "        ],\n",
    "        layers_dict=folium_layers,\n",
    "        center_gdf=ref_nodes_simplified,\n",
    "        center_crs=ref_nodes_simplified.crs,\n",
    "    )\n",
    "\n",
    "    bounds = plot_func.compute_folium_bounds(ref_nodes_simplified)\n",
    "    m.fit_bounds(bounds)\n",
    "    m.save(ref_results_inter_maps_fp + f\"component_gaps_{component_min_distance}_reference.html\")\n",
    "\n",
    "    display(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if len(component_gaps) > 0:\n",
    "    print(\"Interactive map saved at \" + ref_results_inter_maps_fp.lstrip(\"../\") + f\"component_gaps_{component_min_distance}_reference.html\")\n",
    "else:\n",
    "    print(\"There are no component gaps to plot.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Component connectivity\n",
    "\n",
    "Here we visualize differences between how many cells can be reached from each cell. This is a crude measure for network connectivity but has the benefit of being computationally cheap and thus able to quickly highlight stark differences in network connectivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_components_cell_count = eval_func.count_component_cell_reach(\n",
    "    components_df, grid, \"component_ids_ref\"\n",
    ")\n",
    "grid[\"cells_reached_ref\"] = grid[\"component_ids_ref\"].apply(\n",
    "    lambda x: eval_func.count_cells_reached(x, ref_components_cell_count)\n",
    "    if x != \"\"\n",
    "    else 0\n",
    ")\n",
    "\n",
    "grid[\"cells_reached_ref_pct\"] = grid.apply(\n",
    "    lambda x: np.round((x.cells_reached_ref / len(grid)) * 100, 2), axis=1\n",
    ")\n",
    "\n",
    "grid.loc[grid[\"cells_reached_ref_pct\"] == 0, \"cells_reached_ref_pct\"] = np.NAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot percent of cells reachable\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "# norm for color bars\n",
    "cbnorm_reach = colors.Normalize(vmin=0, vmax=100)\n",
    "cbnorm_reach_diff = colors.Normalize(vmin=-100, vmax=100)\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"3.5%\", pad=\"1%\")\n",
    "\n",
    "grid.plot(\n",
    "    cax=cax,\n",
    "    ax=ax,\n",
    "    column=\"cells_reached_ref_pct\",\n",
    "    legend=True,\n",
    "    cmap=pdict[\"seq\"],\n",
    "    norm=cbnorm_reach,\n",
    "    alpha=pdict[\"alpha_grid\"],\n",
    ")\n",
    "\n",
    "ref_edges_simplified.plot(ax=ax, color=pdict[\"ref_emp\"], linewidth=pdict[\"line_base\"])\n",
    "\n",
    "# add no data patches\n",
    "grid[grid[\"count_ref_edges\"].isnull()].plot(\n",
    "    cax=cax,\n",
    "    ax=ax,\n",
    "    facecolor=pdict[\"nodata_face\"],\n",
    "    edgecolor=pdict[\"nodata_edge\"],\n",
    "    linewidth= pdict[\"line_nodata\"],\n",
    "    hatch=pdict[\"nodata_hatch\"],\n",
    "    alpha=pdict[\"alpha_nodata\"],\n",
    ")\n",
    "\n",
    "ax.legend(handles=[nodata_patch], loc=\"upper right\")\n",
    "\n",
    "cx.add_basemap(ax=ax, crs=study_crs, source=cx_tile_2)\n",
    "cx.add_attribution(ax=ax, text=f\"(C) {reference_name}\")\n",
    "txt = ax.texts[-1]\n",
    "txt.set_position([1,0.00])\n",
    "txt.set_ha('right')\n",
    "txt.set_va('bottom')\n",
    "ax.set_title(area_name+f\": {reference_name} percent of cells reachable\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "plot_func.save_fig(fig,ref_results_static_maps_fp + \"percent_cells_reachable_grid_reference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "components_results = {}\n",
    "components_results[\"component_count\"] = len(ref_components)\n",
    "components_results[\"largest_cc_pct_size\"] = largest_cc_pct\n",
    "components_results[\"largest_cc_length\"] = largest_cc_length\n",
    "components_results[\"count_component_gaps\"] = len(component_gaps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summarize_results = {**density_results, **components_results}\n",
    "\n",
    "summarize_results[\"count_dangling_nodes\"] = len(dangling_nodes)\n",
    "summarize_results[\"count_overshoots\"] = len(overshoots)\n",
    "summarize_results[\"count_undershoots\"] = len(undershoot_nodes)\n",
    "\n",
    "# Add total node count and total infrastructure length\n",
    "summarize_results[\"total_nodes\"] = len(ref_nodes_simplified)\n",
    "summarize_results[\"total_length\"] = ref_edges_simplified.infrastructure_length.sum() / 1000\n",
    "\n",
    "summarize_results_df = pd.DataFrame.from_dict(summarize_results, orient=\"index\")\n",
    "summarize_results_df.rename({0: \" \"}, axis=1, inplace=True)\n",
    "\n",
    "# Convert length to km\n",
    "summarize_results_df.loc[\"largest_cc_length\"] = (\n",
    "    summarize_results_df.loc[\"largest_cc_length\"] / 1000\n",
    ")\n",
    "\n",
    "summarize_results_df = summarize_results_df.reindex([\n",
    "    'total_length',\n",
    "    'protected_density_m_sqkm',\n",
    "    'unprotected_density_m_sqkm',\n",
    "    'mixed_density_m_sqkm',\n",
    "    'edge_density_m_sqkm',\n",
    "    'total_nodes',\n",
    "    'count_dangling_nodes',\n",
    "    'node_density_count_sqkm',\n",
    "    'dangling_node_density_count_sqkm',\n",
    "    'count_overshoots',\n",
    "    'count_undershoots',\n",
    "    'component_count',\n",
    "    'largest_cc_length',\n",
    "    'largest_cc_pct_size', \n",
    "    'count_component_gaps'\n",
    "     ])\n",
    "\n",
    "rename_metrics = {\n",
    "    \"total_length\": \"Total infrastructure length (km)\",\n",
    "    \"total_nodes\": \"Nodes\",\n",
    "    \"edge_density_m_sqkm\": \"Bicycle infrastructure density (m/km2)\",\n",
    "    \"node_density_count_sqkm\": \"Nodes per km2\",\n",
    "    \"dangling_node_density_count_sqkm\": \"Dangling nodes per km2\",\n",
    "    \"protected_density_m_sqkm\": \"Protected bicycle infrastructure density (m/km2)\",\n",
    "    \"unprotected_density_m_sqkm\": \"Unprotected bicycle infrastructure density (m/km2)\",\n",
    "    \"mixed_density_m_sqkm\": \"Mixed protection bicycle infrastructure density (m/km2)\",\n",
    "    \"component_count\": \"Components\",\n",
    "    \"largest_cc_pct_size\": \"Largest component's share of network length\",\n",
    "    \"largest_cc_length\": \"Length of largest component (km)\",\n",
    "    \"count_component_gaps\": \"Component gaps\",\n",
    "    \"count_dangling_nodes\": \"Dangling nodes\",\n",
    "    \"count_overshoots\": \"Overshoots\",\n",
    "    \"count_undershoots\": \"Undershoots\",\n",
    "}\n",
    "\n",
    "summarize_results_df.rename(rename_metrics, inplace=True)\n",
    "summarize_results_df.style.pipe(format_ref_style)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "## Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_results = {}\n",
    "\n",
    "all_results[\"network_density\"] = density_results\n",
    "all_results[\"count_overshoots\"] = len(overshoots)\n",
    "all_results[\"count_undershoots\"] = len(undershoot_nodes)\n",
    "all_results[\"dangling_node_count\"] = len(dangling_nodes)\n",
    "all_results[\"simplification_outcome\"] = simplification_results\n",
    "all_results[\"component_analysis\"] = components_results\n",
    "\n",
    "with open(ref_intrinsic_fp, \"w\") as outfile:\n",
    "    json.dump(all_results, outfile)\n",
    "\n",
    "# Save summary dataframe\n",
    "summarize_results_df.to_csv(\n",
    "    ref_results_data_fp + \"intrinsic_summary_results.csv\", index=True\n",
    ")\n",
    "\n",
    "# Save grid with results\n",
    "with open(ref_intrinsic_grid_fp, \"wb\") as f:\n",
    "    pickle.dump(grid, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "from time import strftime\n",
    "print(\"Time of analysis: \" + strftime(\"%a, %d %b %Y %H:%M:%S\"))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "bikedna",
   "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.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "25d89d047ff73a4f08e1f58df0505313a978ba3f16552d4fb4e98a48d36b765b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}