scripts/REFERENCE/2a_initialize_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": [
    "# 2a. Initialize reference data\n",
    "\n",
    "This notebook:\n",
    "- Loads the polygon defining the study area and then creates a grid overlay for the study area.\n",
    "- Loads the reference data.\n",
    "- Processes the reference data to create the network structure and attributes needed in the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Prerequisites &amp; Input/Output</b>\n",
    "<p>\n",
    "<span style=\"font-family:courier;\">config.yml</span> must be set in advance.\n",
    "    \n",
    "Output files of this notebook are saved to the <span style=\"font-family:courier;\">../data/REFERENCE/[study_area]/processed/</span> subfolder and used in notebooks 2b, 3a and 3b.\n",
    "\n",
    "When the network is created from the reference data, each edge is assigned a unique edge ID. To relate the edge ID back to the ID in the original data, the dataframes with network edges, <span style=\"font-family:courier;\">ref_edges</span> and <span style=\"font-family:courier;\">ref_edges_simplified</span>, can be used.\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sections**\n",
    "* [Load data for study area and create analysis grid](#Load-data-for-study-area-and-create-analysis-grid)\n",
    "* [Load and preprocess reference data](#Load-and-preprocess-reference-data)\n",
    "\n",
    "<br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data for study area and create analysis grid\n",
    "\n",
    "This step:\n",
    "\n",
    "- Loads settings for the analysis from the configuration file `config.yml`.\n",
    "- Reads data for the study area.\n",
    "- Creates a grid overlay of the study area, with grid cell size defined in `config.yml`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Troubleshooting</b>\n",
    "<p>\n",
    "If in the cell below the functions located in the <i>bikedna</i> folder do not load properly: Check that <span style=\"font-family:courier;\">pip install -e .</span> was run successfully (see instructions in <span style=\"font-family:courier;\">README.md</span>).\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Load libraries, settings and data\n",
    "\n",
    "import pickle\n",
    "import warnings\n",
    "\n",
    "import contextily as cx\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "import numpy as np\n",
    "\n",
    "from bikedna import evaluation_functions as eval_func\n",
    "from bikedna import graph_functions as graph_func\n",
    "from bikedna import simplification_functions as simp_func\n",
    "from bikedna import plotting_functions as plot_func\n",
    "\n",
    "%run ../settings/plotting.py\n",
    "%run ../settings/yaml_variables.py\n",
    "%run ../settings/paths.py\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data for study area\n",
    "\n",
    "The study area is defined by the user-provided polygon. It will be used for the computation of **global** results, i.e. quality metrics based on all data in the study area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read polygon for study area\n",
    "\n",
    "study_area_poly = gpd.read_file(study_area_poly_fp)\n",
    "\n",
    "if study_area_poly.crs == None:\n",
    "    print(\"Please assign a crs to the study area polygon!\")\n",
    "\n",
    "if study_area_poly.crs != study_crs:\n",
    "    study_area_poly = study_area_poly.to_crs(study_crs)\n",
    "\n",
    "assert study_area_poly.crs == study_crs\n",
    "\n",
    "area = study_area_poly.area.values[0]\n",
    "print(f\"The size of the study area is {area / 1000000:.2f} km2.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Plot study area\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "print(\"This study area will be used for the analysis:\")\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=pdict[\"fsmap\"])\n",
    "    \n",
    "study_area_poly.plot(ax=ax, edgecolor=pdict[\"base\"], facecolor='none',alpha=pdict[\"alpha_back\"],hatch='/')\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "cx.add_basemap(ax, crs=study_area_poly.crs, source=cx_tile_1)\n",
    "\n",
    "ax.set_title(f\"{area_name} study area ({area / 1000000:.2f} sqkm)\")\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + \"area_size_reference\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create analysis grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>User configurations</b>\n",
    "<p>\n",
    "The grid is used to compute <b>local</b> results for many of the quality metrics, i.e. on grid cell level. The cell size of the grid cells must be specified in the configuration file. Smaller grid cells will increase the computation time but give a more granular result. The grid cell size should not be smaller than the scale at which you would expect actual differences in data quality.\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create and plot grid\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "grid = eval_func.create_grid_geometry(study_area_poly, grid_cell_size)\n",
    "\n",
    "grid[\"grid_id\"] = grid.index\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "grid_cell_area = (grid_cell_size / 1000) ** 2  # in km2\n",
    "print(\n",
    "    f\"The grid contains {len(grid)} square cells with a side length of {grid_cell_size} m and an area of {grid_cell_area:.2f} km2.\"\n",
    ")\n",
    "print(\"This grid will be used for local (grid cell level) analysis:\")\n",
    "\n",
    "grid.plot(ax=ax, facecolor=\"none\", edgecolor=pdict[\"base\"], alpha=pdict[\"alpha_back\"],linewidth=0.5)\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "cx.add_basemap(ax, crs=study_area_poly.crs, source=cx_tile_1)\n",
    "\n",
    "ax.set_title(\n",
    "    f\"{area_name} study area ({len(grid)} grid cells, side length {grid_cell_size}m)\"\n",
    ")\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + \"area_grid_reference\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and preprocess reference data\n",
    "\n",
    "This step:\n",
    "\n",
    "- Creates a network from the reference data.\n",
    "- Projects it to the chosen CRS.\n",
    "- Clips the data to the polygon defining the study area.\n",
    "- Measures the infrastructure length of the edges based on the geometry type and whether they allow for bidirectional travel or not.\n",
    "- Simplifies the network.\n",
    "- Creates copies of all edge and node data sets indexed by their intersecting grid cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Network data model**\n",
    "\n",
    "In BikeDNA, all input data are converted to a network structure consisting of *nodes* and *edges*. Edges represents the actual infrastructure, such as bike lanes and paths, while nodes represents the start and end points for the edges, as well as all intersections. For further details, read more about the [network data model](https://networkx.org/documentation/stable/reference/introduction.html).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details><summary><b>Network simplification</b></summary>\n",
    "\n",
    "When converting a set of linestrings to network (graph) format, not all vertices (nodes) are of equal meaning. For the geometry of edges, all nodes are used to define the shape of the edge. For the topology of the network, however, only those nodes that are endpoints or intersection points with other edges are needed, while all other (so-called 'interstitial') nodes do not add any information. To compare the structure and true ratio between nodes and edges in a network, a simplified network representation which only includes nodes at endpoints and intersections, or where the value of important attributes changes, is required. An additional advantage of simplifying the network is the resulting substantial reduction of the number of nodes and edges, which makes computational routines much faster.\n",
    "\n",
    "Note that the simplification function only simplifies and merges edges that have identical values for important attributes such as geometry type, bicycle direction, and protection level. \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",
    "As part of the simplification routine, in cases where there are several edges between the same pair of nodes ('parallel edges' or 'multiedges'), only one of the edges is retained.\n",
    "\n",
    "To read more about the modified OSMnx simplification [(Boeing, 2017)](https://geoffboeing.com/publications/osmnx-complex-street-networks/) used here, we refer to this [GitHub repository](https://github.com/anerv/NERDS_osmnx) which contains both the simplification functions, explanation of the logic and a demonstration.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read reference data\n",
    "org_ref_data = gpd.read_file(reference_fp)\n",
    "\n",
    "# Create data with LineStrings only defined by end and start coordinate\n",
    "ref_data = graph_func.unzip_linestrings(org_ref_data, reference_id_col)\n",
    "\n",
    "# Create column with id called 'osmid'\n",
    "ref_data[\"osmid\"] = ref_data[\"new_edge_id\"]\n",
    "ref_data.drop([\"new_edge_id\"], axis=1, inplace=True)\n",
    "\n",
    "# Reproject\n",
    "if ref_data.crs == None:\n",
    "    print(\"Please assign a crs to the study area polygon!\")\n",
    "\n",
    "if ref_data.crs != study_crs:\n",
    "    ref_data = ref_data.to_crs(study_crs)\n",
    "\n",
    "assert ref_data.crs == study_crs\n",
    "\n",
    "# Clip reference data to study area poly\n",
    "ref_data = ref_data.clip(study_area_poly)\n",
    "\n",
    "ref_data = ref_data[ref_data.geometry.length!=0.0]\n",
    "\n",
    "# Convert to osmnx graph object\n",
    "graph_ref = graph_func.create_osmnx_graph(ref_data)\n",
    "\n",
    "ref_nodes, ref_edges = ox.graph_to_gdfs(graph_ref)\n",
    "\n",
    "# Overview of data\n",
    "graph_area = eval_func.get_graph_area(ref_nodes, study_area_poly, study_crs)\n",
    "print(\n",
    "    f\"The {reference_name} data covers an area of {graph_area / 1000000:.2f} km2.\"\n",
    ")\n",
    "print(\"\\n\")\n",
    "\n",
    "# Classify edges as protected or unprotected\n",
    "ref_edges = eval_func.define_protected_unprotected(\n",
    "    ref_edges, ref_bicycle_infrastructure_type\n",
    ")\n",
    "for key, value in ref_edges.protected.value_counts().items():\n",
    "    perc = np.round(100*value/len(ref_edges),2)\n",
    "    print(f\"Edges where the protection level is '{key}': {value} out of {len(ref_edges)} ({perc}%)\")\n",
    "print(\"\\n\")\n",
    "\n",
    "# Summarize attribute values - if there are columns with values for each row instead of a global setting\n",
    "if bicycle_bidirectional in ref_edges.columns:\n",
    "    for key, value in ref_edges[bicycle_bidirectional].value_counts().items():\n",
    "        perc = np.round(100*value/len(ref_edges),2)\n",
    "        print(f\"Edges where 'bicycle_bidirectional' is {key}: {value} out of {len(ref_edges)} ({perc}%)\")\n",
    "else: \n",
    "    print(\"Using global settings for cycling direction.\")\n",
    "print(\"\\n\")\n",
    "\n",
    "if reference_geometries in ref_edges.columns:\n",
    "    for key, value in ref_edges[reference_geometries].value_counts().items():\n",
    "        perc = np.round(100*value/len(ref_edges),2)\n",
    "        print(f\"Edges where the geometry type is '{key}': {value} out of {len(ref_edges)}({perc}%)\")\n",
    "else:\n",
    "    print(\"Using global settings for geometry type.\")\n",
    "print(\"\\n\")\n",
    "\n",
    "# Add attributes to graph\n",
    "protected_dict = ref_edges[\"protected\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref, protected_dict, \"protected\")\n",
    "\n",
    "# Add 'multiple edge' attribute to each edge of the graph\n",
    "for e in graph_ref.edges:\n",
    "    graph_ref.edges[e][\"multiedge\"] = 0\n",
    "\n",
    "# Count multiple edges and store info at key 0\n",
    "# (simplification routine keeps only key 0)\n",
    "for u, v, k in graph_ref.edges:\n",
    "    if k > 0:\n",
    "        graph_ref.edges[u, v, 0][\"multiedge\"] += 1\n",
    "\n",
    "# Remove geometry attribute (required by simplification function)\n",
    "for n1, n2, d in graph_ref.edges(data=True):\n",
    "    d.pop(\"geometry\", None)\n",
    "\n",
    "# Simplify - the graph is simplified as much as possible, but without combining edges with conflicting values of bicycle infrastructure types\n",
    "simplify_cols = [reference_geometries, bicycle_bidirectional, \"protected\"]\n",
    "simplify_cols = [s for s in simplify_cols if s in ref_edges.columns]\n",
    "graph_ref_simplified = simp_func.simplify_graph(\n",
    "    graph_ref, attributes=simplify_cols, remove_rings=False\n",
    ")\n",
    "\n",
    "graph_ref = ox.get_undirected(graph_ref)\n",
    "ref_nodes, ref_edges = ox.graph_to_gdfs(graph_ref)\n",
    "graph_ref_simplified = ox.get_undirected(graph_ref_simplified)\n",
    "\n",
    "ref_nodes_simplified, ref_edges_simplified = ox.graph_to_gdfs(graph_ref_simplified)\n",
    "\n",
    "ref_edges[\"length\"] = ref_edges.geometry.length\n",
    "ref_edges_simplified[\"length\"] = ref_edges_simplified.geometry.length\n",
    "\n",
    "# Measure the length of bicycle infrastructure (!=from length of edges)\n",
    "if (\n",
    "    reference_geometries not in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == str\n",
    "):\n",
    "\n",
    "    ref_edges_simplified[\"infrastructure_length\"] = ref_edges_simplified.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=x[reference_geometries],\n",
    "            bidirectional=x[bicycle_bidirectional],\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "elif (\n",
    "    reference_geometries in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == bool\n",
    "):\n",
    "\n",
    "    ref_edges_simplified[\"infrastructure_length\"] = ref_edges_simplified.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=reference_geometries,\n",
    "            bidirectional=bicycle_bidirectional,\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "\n",
    "elif (\n",
    "    reference_geometries in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == str\n",
    "):\n",
    "\n",
    "    ref_edges_simplified[\"infrastructure_length\"] = ref_edges_simplified.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=reference_geometries,\n",
    "            bidirectional=x[bicycle_bidirectional],\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "\n",
    "elif (\n",
    "    reference_geometries not in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == bool\n",
    "):\n",
    "\n",
    "    ref_edges_simplified[\"infrastructure_length\"] = ref_edges_simplified.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=x[reference_geometries],\n",
    "            bidirectional=bicycle_bidirectional,\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "else:\n",
    "    print(\n",
    "        f\"Invalid data types for the variables reference_geometetries and/or bicycle_bidirectional!\"\n",
    "    )\n",
    "\n",
    "print(\n",
    "    f\"The length of the {reference_name} network is {ref_edges_simplified.infrastructure_length.sum()/1000 :.2f} km.\"\n",
    ")\n",
    "\n",
    "\n",
    "# Creating coulumns with unique feature id (required by some functions)\n",
    "ref_edges[\"edge_id\"] = ref_edges.reset_index().index\n",
    "ref_edges_simplified[\"edge_id\"] = ref_edges_simplified.reset_index().index\n",
    "\n",
    "ref_id_dict = ref_edges[\"edge_id\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref, ref_id_dict, \"edge_id\")\n",
    "\n",
    "ref_id_dict_s = ref_edges_simplified[\"edge_id\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref_simplified, ref_id_dict_s, \"edge_id\")\n",
    "\n",
    "# Saving data to graphs\n",
    "ref_l_dict = ref_edges[\"length\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref, ref_l_dict, \"length\")\n",
    "\n",
    "ref_l_dict_s = ref_edges_simplified[\"length\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref_simplified, ref_l_dict_s, \"length\")\n",
    "\n",
    "ref_il_dict = ref_edges_simplified[\"infrastructure_length\"].to_dict()\n",
    "nx.set_edge_attributes(graph_ref_simplified, ref_il_dict, \"infrastructure_length\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot global statistics (by network length)\n",
    "\n",
    "set_renderer(renderer_plot)\n",
    "ymax = ref_edges.length.sum() * 1.2 /1000 # conversion m to km\n",
    "\n",
    "# attribute: bidirectional\n",
    "if bicycle_bidirectional in ref_edges.columns:\n",
    "    \n",
    "    x = []\n",
    "    y = []\n",
    "\n",
    "    for value in np.unique(ref_edges[bicycle_bidirectional]):\n",
    "        x.append(str(value))\n",
    "        y.append(ref_edges[ref_edges[bicycle_bidirectional]==value].length.sum()/1000) # conversion m to km\n",
    "\n",
    "    plot_func.make_bar_plot(\n",
    "        data = y,\n",
    "        bar_labels = x,\n",
    "        y_label = \"Network edge lengths [km]\",\n",
    "        x_positions = [i for i in range(len(x))],\n",
    "        title = f\"{area_name}\\n {reference_name} edges by bidirectionality\",\n",
    "        bar_colors = len(x)*[pdict[\"ref_base\"]],\n",
    "        filepath = ref_results_plots_fp + \"edge_global_bidirectional\",\n",
    "        ylim = ymax,\n",
    "        figsize=pdict[\"fsbar_small\"]\n",
    "    );\n",
    "\n",
    "# attribute: protection\n",
    "x = []\n",
    "y = []\n",
    "\n",
    "for value in np.unique(ref_edges.protected):\n",
    "    x.append(value)\n",
    "    y.append(ref_edges[ref_edges[\"protected\"]==value].length.sum()/1000) # conversion m to km\n",
    "    \n",
    "plot_func.make_bar_plot(\n",
    "    data = y,\n",
    "    bar_labels = x,\n",
    "    y_label = \"Network edge lengths [km]\",\n",
    "    x_positions = [i for i in range(len(x))],\n",
    "    title = f\"{area_name}\\n {reference_name} edges by infrastructure type\",\n",
    "    bar_colors = len(x)*[pdict[\"ref_base\"]],\n",
    "    filepath = ref_results_plots_fp + \"edge_global_infrastructuretype\",\n",
    "    ylim = ymax,\n",
    "    figsize=pdict[\"fsbar_small\"]\n",
    ");\n",
    "\n",
    "# attribute: bicycle geometries (mapping type)\n",
    "if reference_geometries in ref_edges.columns:\n",
    "\n",
    "    x = []\n",
    "    y = []\n",
    "\n",
    "    for value in np.unique(ref_edges[reference_geometries]):\n",
    "        x.append(value)\n",
    "        y.append(ref_edges[ref_edges[reference_geometries]==value].length.sum()/1000) # conversion m to km\n",
    "\n",
    "    plot_func.make_bar_plot(\n",
    "        data = y,\n",
    "        bar_labels = x,\n",
    "        y_label = \"Network edge lengths [km]\",\n",
    "        x_positions = [i for i in range(len(x))],\n",
    "        title = f\"{area_name}\\n {reference_name} edges by mapping type\",\n",
    "        bar_colors = len(x)*[pdict[\"ref_base\"]],\n",
    "        filepath = ref_results_plots_fp + \"edge_global_mappingtype\",\n",
    "        ylim = ymax,\n",
    "        figsize=pdict[\"fsbar_small\"]\n",
    "    );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot network\n",
    "\n",
    "set_renderer(renderer_map)\n",
    "fig, ax = plt.subplots(1, figsize=pdict[\"fsmap\"])\n",
    "\n",
    "ref_edges_simplified.plot(ax=ax, color=pdict[\"ref_base\"], linewidth=pdict[\"line_base\"])\n",
    "ref_nodes_simplified.plot(ax=ax, color=pdict[\"ref_emp\"], markersize=pdict[\"mark_emp\"])\n",
    "\n",
    "study_area_poly.plot(ax=ax, edgecolor=pdict[\"base\"], facecolor=\"None\", linewidth=1)\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "cx.add_basemap(ax, crs=study_area_poly.crs, source=cx_tile_2)\n",
    "\n",
    "ax.set_title(f\"{area_name}, {reference_name} network\")\n",
    "\n",
    "plot_func.save_fig(fig, ref_results_static_maps_fp + \"area_network_reference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Joining grid cell id to network elements\n",
    "\n",
    "ref_edges_simp_joined = gpd.overlay(\n",
    "    ref_edges_simplified.reset_index(), grid, how=\"intersection\", keep_geom_type=True\n",
    ")\n",
    "ref_edges_joined = gpd.overlay(\n",
    "    ref_edges.reset_index(), grid, how=\"intersection\", keep_geom_type=True\n",
    ")\n",
    "\n",
    "ref_nodes_simp_joined = gpd.overlay(\n",
    "    ref_nodes_simplified.reset_index(), grid, how=\"intersection\", keep_geom_type=True\n",
    ")\n",
    "ref_nodes_joined = gpd.overlay(\n",
    "    ref_nodes.reset_index(), grid, how=\"intersection\", keep_geom_type=True\n",
    ")\n",
    "\n",
    "# Count features in each grid cell\n",
    "all_data_ref = [\n",
    "    ref_edges_joined,\n",
    "    ref_nodes_joined,\n",
    "    ref_edges_simp_joined,\n",
    "    ref_nodes_simp_joined,\n",
    "]\n",
    "labels_ref = [\"ref_edges\", \"ref_nodes\", \"ref_simplified_edges\", \"ref_simplified_nodes\"]\n",
    "\n",
    "for data, label in zip(all_data_ref, labels_ref):\n",
    "\n",
    "    df = eval_func.count_features_in_grid(data, label)\n",
    "\n",
    "    grid = eval_func.merge_results(grid, df, \"left\")\n",
    "\n",
    "\n",
    "# Recomputing infrastructure length for data joined to grid\n",
    "ref_edges_simp_joined[\"length\"] = ref_edges_simp_joined[\"geometry\"].length\n",
    "\n",
    "if (\n",
    "    reference_geometries not in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == str\n",
    "):\n",
    "\n",
    "    ref_edges_simp_joined[\"infrastructure_length\"] = ref_edges_simp_joined.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=x[reference_geometries],\n",
    "            bidirectional=x[bicycle_bidirectional],\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "elif (\n",
    "    reference_geometries in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == bool\n",
    "):\n",
    "\n",
    "    ref_edges_simp_joined[\"infrastructure_length\"] = ref_edges_simp_joined.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=reference_geometries,\n",
    "            bidirectional=bicycle_bidirectional,\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "elif (\n",
    "    reference_geometries in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == str\n",
    "):\n",
    "\n",
    "    ref_edges_simp_joined[\"infrastructure_length\"] = ref_edges_simp_joined.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=reference_geometries,\n",
    "            bidirectional=x[bicycle_bidirectional],\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "\n",
    "elif (\n",
    "    reference_geometries not in [\"true_geometries\", \"centerline\"]\n",
    "    and type(bicycle_bidirectional) == bool\n",
    "):\n",
    "\n",
    "    ref_edges_simp_joined[\"infrastructure_length\"] = ref_edges_simp_joined.apply(\n",
    "        lambda x: eval_func.measure_infrastructure_length(\n",
    "            edge=x.geometry,\n",
    "            geometry_type=[reference_geometries],\n",
    "            bidirectional=bicycle_bidirectional,\n",
    "            bicycle_infrastructure=\"yes\",\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "else:\n",
    "    print(\n",
    "        f\"Invalid data types for the variables reference_geometries and/or bicycle_bidirectional!\"\n",
    "    )\n",
    "\n",
    "assert round(ref_edges_simplified.infrastructure_length.sum() / 1000, 1) == round(\n",
    "    ref_edges_simp_joined.infrastructure_length.sum() / 1000, 1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "noex"
    ]
   },
   "outputs": [],
   "source": [
    "# Save reference data\n",
    "exec(open(\"../settings/save_refdata.py\").read())"
   ]
  },
  {
   "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
}