diff --git a/FixBikeMVP.ipynb b/FixBikeMVP.ipynb index f6e51f6..c927ef9 100644 --- a/FixBikeMVP.ipynb +++ b/FixBikeMVP.ipynb @@ -15,17 +15,14 @@ "id": "1b6a69660258f567", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:44:35.560984Z", - "start_time": "2026-03-05T14:44:34.251482Z" + "end_time": "2026-03-06T13:24:09.401596Z", + "start_time": "2026-03-06T13:24:08.014186Z" } }, "source": [ "# import packages\n", "import pandas as pd\n", "import osmnx as ox\n", - "import networkx as nx\n", - "import random\n", - "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import os\n", "\n", @@ -48,8 +45,8 @@ "id": "9173dd015d3280f3", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:44:43.976210Z", - "start_time": "2026-03-05T14:44:43.954651Z" + "end_time": "2026-03-06T13:24:10.609718Z", + "start_time": "2026-03-06T13:24:10.596152Z" } }, "source": [ @@ -79,8 +76,8 @@ "id": "a9a51035deb92425", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:45:07.394672Z", - "start_time": "2026-03-05T14:44:45.668432Z" + "end_time": "2026-03-06T13:24:34.025385Z", + "start_time": "2026-03-06T13:24:12.224519Z" } }, "source": [ @@ -101,8 +98,8 @@ "id": "a8b4eed5", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:45:15.988833Z", - "start_time": "2026-03-05T14:45:13.924067Z" + "end_time": "2026-03-06T13:24:46.446529Z", + "start_time": "2026-03-06T13:24:44.423659Z" } }, "source": [ @@ -137,8 +134,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:45:22.101582Z", - "start_time": "2026-03-05T14:45:21.982420Z" + "end_time": "2026-03-06T13:24:47.841310Z", + "start_time": "2026-03-06T13:24:47.728797Z" } }, "cell_type": "code", @@ -164,8 +161,8 @@ { "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:47:06.726898Z", - "start_time": "2026-03-05T14:45:23.501207Z" + "end_time": "2026-03-06T13:26:35.471171Z", + "start_time": "2026-03-06T13:24:49.019530Z" } }, "cell_type": "code", @@ -179,8 +176,8 @@ "id": "e2670fe4", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:47:09.810662Z", - "start_time": "2026-03-05T14:47:09.732303Z" + "end_time": "2026-03-06T13:31:14.903844Z", + "start_time": "2026-03-06T13:31:14.725758Z" } }, "source": [ @@ -222,8 +219,8 @@ "id": "ac84d5bc", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:47:11.808025Z", - "start_time": "2026-03-05T14:47:11.457935Z" + "end_time": "2026-03-06T13:31:16.805584Z", + "start_time": "2026-03-06T13:31:16.458067Z" } }, "source": [ @@ -248,8 +245,8 @@ "id": "c6bc1fa5", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:47:13.884562Z", - "start_time": "2026-03-05T14:47:13.788483Z" + "end_time": "2026-03-06T13:31:18.314233Z", + "start_time": "2026-03-06T13:31:18.226835Z" } }, "source": "G = weigh_edges(G, penalty)", @@ -271,8 +268,8 @@ "id": "602e50f5", "metadata": { "ExecuteTime": { - "end_time": "2026-03-05T14:47:22.674995Z", - "start_time": "2026-03-05T14:47:22.515390Z" + "end_time": "2026-03-06T13:31:20.049312Z", + "start_time": "2026-03-06T13:31:19.891469Z" } }, "source": [ @@ -288,7 +285,7 @@ ] } ], - "execution_count": 11 + "execution_count": 10 }, { "cell_type": "markdown", @@ -297,39 +294,32 @@ "source": [ "### Create list of all potential gaps\n", "\n", - "Here defined as: contact node pair combinations **within `maxgap` euclidean distance from each other** " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a02eeca1", - "metadata": {}, - "outputs": [], - "source": [ - "potential_gaps = []\n", - "\n", - "for node in contact_nodes:\n", - " node_buffer = nodes_gdf.loc[node, \"geometry\"].buffer(maxgap)\n", - " q = nodes_gdf.sindex.query(node_buffer, predicate=\"intersects\")\n", - " neighbours = list(nodes_gdf.iloc[q].index)\n", - " # convention: sort by ascending OSMID...\n", - " node_pairs = [tuple(sorted(z)) for z in zip([node]*len(neighbours), neighbours)]\n", - " potential_gaps += node_pairs\n", - "\n", - "# ... so that we can easily deduplicate\n", - "potential_gaps = list(set(potential_gaps))" + "Here defined as: contact node pair combinations **within `maxgap` euclidean distance from each other**" ] }, { "cell_type": "code", - "execution_count": null, "id": "a01da088", - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2026-03-06T13:31:24.244710Z", + "start_time": "2026-03-06T13:31:21.777852Z" + } + }, "source": [ + "potential_gaps = find_potential_gaps(contact_nodes, nodes_gdf, maxgap)\n", "print(\"potential gaps found:\", len(potential_gaps))" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "potential gaps found: 142482\n" + ] + } + ], + "execution_count": 11 }, { "cell_type": "markdown", @@ -345,35 +335,16 @@ }, { "cell_type": "code", - "execution_count": null, "id": "43c578c3", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2026-03-06T13:31:31.246886Z", + "start_time": "2026-03-06T13:31:26.417817Z" + } + }, + "source": "found_gaps, found_gaps_nsp = find_actual_gaps(G, potential_gaps)", "outputs": [], - "source": [ - "pbi_dict = nx.get_edge_attributes(G, \"pbi\")\n", - "\n", - "found_gaps = []\n", - "found_gaps_nsp = [] # naive shortest paths (by length, in node list format)\n", - "\n", - "for i, gap in enumerate(potential_gaps):\n", - " u, v = gap\n", - " nodelist = nx.shortest_path(\n", - " G=G,\n", - " source=u,\n", - " target=v, \n", - " weight=\"length\"\n", - " )\n", - " pbis = set([pbi_dict[tuple(sorted(z))] for z in zip(nodelist, nodelist[1:])])\n", - " \n", - " # confirm that it is an actual gap if it consists only of pbi==0 infra:\n", - " if pbis==set([0]): \n", - " found_gaps.append(gap)\n", - " found_gaps_nsp.append(nodelist)\n", - " \n", - " # # (manual timer)\n", - " # if i % 100000 == 0:\n", - " # print(i)\n" - ] + "execution_count": 12 }, { "cell_type": "markdown", @@ -385,17 +356,22 @@ }, { "cell_type": "code", - "execution_count": null, "id": "b5ca9fdf440935d8", - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2026-03-06T13:31:33.093955Z", + "start_time": "2026-03-06T13:31:33.006586Z" + } + }, "source": [ "edgelist = []\n", "for nodelist in found_gaps_nsp:\n", " edgelist += [tuple(sorted(z)) for z in zip(nodelist, nodelist[1:])]\n", "# deduplicate\n", "edgelist = list(set(edgelist))" - ] + ], + "outputs": [], + "execution_count": 13 }, { "cell_type": "markdown", @@ -407,38 +383,16 @@ }, { "cell_type": "code", - "execution_count": null, "id": "b1f0c683", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2026-03-06T13:32:08.385205Z", + "start_time": "2026-03-06T13:31:34.994480Z" + } + }, + "source": "ebc = compute_local_betweenness_centrality(G, nodes_gdf, radius)", "outputs": [], - "source": [ - "# set current ebc value of all G edges to 0\n", - "for edge in G.edges:\n", - " G.edges[edge][\"ebc\"] = 0\n", - "\n", - "# create dict that will be updated at each step\n", - "ebc = nx.get_edge_attributes(G, \"ebc\")\n", - "\n", - "# for each node, compute \"local\" ebc (buffered with radius!)\n", - "# for comp feas, now only subset of randomly drawn 100 nodes\n", - "random.seed(1312)\n", - "random_nodes = random.choices(list(G.nodes), k = 100)\n", - "for node in random_nodes:\n", - " node_buffer = nodes_gdf.loc[node, \"geometry\"].buffer(radius)\n", - " q = nodes_gdf.sindex.query(node_buffer, predicate=\"intersects\")\n", - " neighbours = list(nodes_gdf.iloc[q].index)\n", - " local_ebc = nx.edge_betweenness_centrality_subset(\n", - " G=G,\n", - " sources=[node],\n", - " targets=neighbours,\n", - " normalized=False, # important! otherwise the addition makes no sense\n", - " weight=\"weight\"# using penalty for non-pbi\n", - " )\n", - "\n", - " # update ebc dictionary\n", - " for k, v in local_ebc.items():\n", - " ebc[k] += v # updating ebc!!" - ] + "execution_count": 14 }, { "cell_type": "markdown", @@ -450,19 +404,16 @@ }, { "cell_type": "code", - "execution_count": null, "id": "4c920f49e57db5c", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2026-03-06T13:32:13.235834Z", + "start_time": "2026-03-06T13:32:12.972818Z" + } + }, + "source": "Bs = rank_gaps_by_b(found_gaps_nsp, G, ebc)", "outputs": [], - "source": [ - "Bs = []\n", - "for nodelist in found_gaps_nsp:\n", - " edgelist = [tuple(sorted(z)) for z in zip(nodelist, nodelist[1:])]\n", - " lengths = np.array([G.edges[edge][\"length\"] for edge in edgelist])\n", - " ebcs = np.array([ebc[edge] for edge in edgelist])\n", - " B = sum(lengths * ebcs)/sum(lengths)\n", - " Bs.append(B)" - ] + "execution_count": 15 }, { "cell_type": "markdown", diff --git a/fixbikenet/functions.py b/fixbikenet/functions.py index fb6336a..f7b47eb 100644 --- a/fixbikenet/functions.py +++ b/fixbikenet/functions.py @@ -1,4 +1,7 @@ import yaml +import networkx as nx +import random +import numpy as np def map_edges_to_bike_infrastructure(g): """ @@ -106,4 +109,149 @@ def find_contact_nodes(G): pbis = set([G.edges[edge]["pbi"] for edge in G.edges(node)]) if len(pbis) == 2: contact_nodes.append(node) - return contact_nodes \ No newline at end of file + return contact_nodes + +def find_potential_gaps(contact_nodes, nodes_gdf, maxgap): + """ + finds potential gaps in protected bicycle network, corresponding to two contact nodes that are within maxgap euclidean distance of each other + + Parameters + ---------- + contact_nodes: list + list of all nodes that fulfill criteria to be a contact node + nodes_gdf: geopandas.GeoDataFrame + all nodes in street network + maxgap: int + user defined maximal euclidean distance between two contact nodes + + Returns + ------- + potential_gaps: list + all unique potential gaps in protected bicycle network + """ + potential_gaps = [] + + for node in contact_nodes: + node_buffer = nodes_gdf.loc[node, "geometry"].buffer(maxgap) + q = nodes_gdf.sindex.query(node_buffer, predicate="intersects") + neighbours = list(nodes_gdf.iloc[q].index) + neighbours.remove(node) + # convention: sort by ascending OSMID... + node_pairs = [tuple(sorted(z)) for z in zip([node] * len(neighbours), neighbours)] + potential_gaps += node_pairs + + # ... so that we can easily deduplicate + potential_gaps = list(set(potential_gaps)) + return potential_gaps + +def find_actual_gaps(G, potential_gaps): + """ + determines which potential gaps are actual gaps by finding paths between all contact nodes and only keeping the gaps that have no protected bike infrastructure + + Parameters + ---------- + G: networkx.Graph + undirected simple graph representing the street network with weighted edges + potential_gaps: list + all unique potential gaps in protected bicycle network + + Returns + ------- + found_gaps: list + list of all gaps in protected bicycle network + found_gaps_nsp: list + list of paths in network for all gaps in protected bicycle network + """ + pbi_dict = nx.get_edge_attributes(G, "pbi") + + found_gaps = [] + found_gaps_nsp = [] # naive shortest paths (by length, in node list format) + + for i, gap in enumerate(potential_gaps): + u, v = gap + nodelist = nx.shortest_path( + G=G, + source=u, + target=v, + weight="length" + ) + pbis = set([pbi_dict[tuple(sorted(z))] for z in zip(nodelist, nodelist[1:])]) + + # confirm that it is an actual gap if it consists only of pbi==0 infra: + if pbis == set([0]): + found_gaps.append(gap) + found_gaps_nsp.append(nodelist) + return found_gaps, found_gaps_nsp + +def compute_local_betweenness_centrality(G, nodes_gdf, radius): + """ + computes weighted betweenness centrality for paths within radius + + Parameters + ---------- + G: networkx.Graph + undirected simple graph representing the street network with weighted edges + nodes_gdf: geopandas.GeoDataFrame + all nodes in street network + radius: int + maximum length of path for betweennessn centrality calculation, set by user + + Returns + ------- + ebc: dict + local betweenness centrality values for all edges in network + """ + # set current ebc value of all G edges to 0 + for edge in G.edges: + G.edges[edge]["ebc"] = 0 + + # create dict that will be updated at each step + ebc = nx.get_edge_attributes(G, "ebc") + + # for each node, compute "local" ebc (buffered with radius!) + # for comp feas, now only subset of randomly drawn 100 nodes + random.seed(1312) + random_nodes = random.choices(list(G.nodes), k=100) + for node in random_nodes: + node_buffer = nodes_gdf.loc[node, "geometry"].buffer(radius) + q = nodes_gdf.sindex.query(node_buffer, predicate="intersects") + neighbours = list(nodes_gdf.iloc[q].index) + local_ebc = nx.edge_betweenness_centrality_subset( + G=G, + sources=[node], + targets=neighbours, + normalized=False, # important! otherwise the addition makes no sense + weight="weight" # using penalty for non-pbi + ) + + # update ebc dictionary + for k, v in local_ebc.items(): + ebc[k] += v # updating ebc!! + return ebc + +def rank_gaps_by_b(found_gaps_nsp, G, ebc): + """ + calculates b for all gaps + + Parameters + ---------- + found_gaps_nsp: list + list of paths in network for all gaps in protected bicycle network + G: networkx.Graph + undirected simple graph representing the street network with weighted edges + ebc: dict + local betweenness centrality values for all edges in network + + Returns + ------- + Bs: list + list of values of b for all gaps in protected bicycle network + """ + Bs = [] + for nodelist in found_gaps_nsp: + edgelist = [tuple(sorted(z)) for z in zip(nodelist, nodelist[1:])] + lengths = np.array([G.edges[edge]["length"] for edge in edgelist]) + ebcs = np.array([ebc[edge] for edge in edgelist]) + B = sum(lengths * ebcs) / sum(lengths) + Bs.append(B) + return Bs \ No newline at end of file diff --git a/tests/test_functions.py b/tests/test_functions.py index 802587a..a34e9fd 100644 --- a/tests/test_functions.py +++ b/tests/test_functions.py @@ -1,6 +1,7 @@ import pytest -import networkx as nx from networkx.utils.misc import graphs_equal +import geopandas as gpd +from shapely.geometry import Point from fixbikenet.functions import * @@ -83,4 +84,86 @@ def create_validation_contact_nodes(): return contact_nodes def test_find_contact_nodes(create_weighted_graph, create_validation_contact_nodes): - assert find_contact_nodes(create_weighted_graph) == create_validation_contact_nodes \ No newline at end of file + assert find_contact_nodes(create_weighted_graph) == create_validation_contact_nodes + +@pytest.fixture +def create_test_contact_nodes(): + contact_nodes = [1,3,4] + return contact_nodes + +@pytest.fixture +def create_maxgap(): + maxgap = 50 + return maxgap + +@pytest.fixture +def create_test_nodes_gdf(): + nodes ={'osmid':[1,2,3,4],'geometry':[Point(1,1), Point(1000,1000), Point(3,3), Point(2000,2000)]} + nodes_gdf = gpd.GeoDataFrame(nodes) + nodes_gdf.set_index('osmid',inplace=True) + return nodes_gdf + +@pytest.fixture +def create_validation_gaps(): + gaps = [(1,3)] + return gaps + +def test_find_potential_gaps(create_test_contact_nodes, create_maxgap, create_test_nodes_gdf, create_validation_gaps): + assert find_potential_gaps(create_test_contact_nodes,create_test_nodes_gdf, create_maxgap) == create_validation_gaps + +@pytest.fixture +def create_graph_for_routing(): + G = nx.Graph() + G.add_edges_from([(1, 2), (1, 3), (2, 3), (3,4)]) + attributes = {(1, 2): {'length': 10, 'pbi': True, 'weight': 10}, + (1, 3): {'length': 4, 'pbi': False, 'weight': 20}, + (2, 3): {'length': 5, 'pbi': True, 'weight': 5}, + (3, 4): {'length': 6, 'pbi': False, 'weight': 30}} + nx.set_edge_attributes(G, attributes) + return G + +@pytest.fixture +def create_potential_gaps(): + potential_gaps = [(1,4), (2,3), (3,4)] + return potential_gaps + +@pytest.fixture +def create_validation_found_gap_paths(): + found_gap_paths = ([(1,4), (3,4)], [[1,3,4],[3,4]]) + return found_gap_paths + +def test_find_actual_gaps(create_graph_for_routing, create_potential_gaps, create_validation_found_gap_paths): + assert find_actual_gaps(create_graph_for_routing, create_potential_gaps) == create_validation_found_gap_paths + +@pytest.fixture +def create_betweenness_nodes(): + nodes = {'osmid': [1, 2, 3, 4], 'geometry': [Point(1, 1), Point(5000, 5000), Point(3, 3), Point(1000, 1000)]} + nodes_gdf = gpd.GeoDataFrame(nodes) + nodes_gdf.set_index('osmid', inplace=True) + return nodes_gdf + +@pytest.fixture +def create_radius(): + radius = 2000 + return radius + +@pytest.fixture +def create_validation_ebc(): + ebc = {(1,2):54.5, (1,3): 0, (2,3): 54.5, (3,4): 56.5} + return ebc + +def test_compute_local_betweenness_centrality(create_graph_for_routing, create_betweenness_nodes, create_radius, create_validation_ebc): + assert compute_local_betweenness_centrality(create_graph_for_routing, create_betweenness_nodes, create_radius) == create_validation_ebc + +@pytest.fixture +def create_found_paths(): + found_paths = [(1,3,4),(3,4)] + return found_paths + +@pytest.fixture +def create_validation_bs(): + validation_bs = [33.9,56.5] + return validation_bs + +def test_rank_gaps_by_b(create_found_paths, create_graph_for_routing, create_validation_ebc, create_validation_bs): + assert rank_gaps_by_b(create_found_paths, create_graph_for_routing, create_validation_ebc) == create_validation_bs \ No newline at end of file