Feature/file import

examples/mwe.py

@@ -3,8 +3,8 @@
 import fixbikenet as fbn
 
 gaps = fbn.fixbikenet(
-    city_name="Frederiksberg municipality",
+    city_query="Frederiksberg municipality",
     export_file_format="geojson",
 )
 
-# data is saved in current working directory, as gaps.gpkg
+# data is saved in directory ./results

fixbikenet/fixbikenet.py

@@ -2,40 +2,53 @@
 import osmnx as ox
 import os
 import matplotlib.pyplot as plt
+from collections import defaultdict
 
 # import functions
 from fixbikenet.functions import *
 
 def fixbikenet(
-    city_name,
+    city_query,
     proj_crs = "3857",
     radius = 2500,
     maxgap = 1000,
-    penalty = {0: 5, 1: 1},
+    penalty = {0: 1.5, 1: 1},
     export_data = True,
+    city_id = None,
     export_file_format="geojson",
     export_plot=False,
+    import_files={},
 ):
     """
     Finds gaps in bicycle networks and returns the 100 that are the most important to fill.
     Parameters
     ----------
-    city_name : str
+    city_query : str
         name of the city that the analysis should be performed on
     proj_crs : str, default '3857'
         coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator)
     radius : int, default 2500
         cut-off length for computation of local betweenness centrality, in meters
     maxgap : int, default 1000
         maximum distance between node pairs to be considered as a potential gap
-    penalty : dict, default {0:5, 1: 1}
+    penalty : dict, default {0:1.5, 1: 1}
         weighing for shortest path calculations, where streets without protected bike infrastructure (pbi) get penalized
     export_data : bool, optional, default True
         If set to True, data will be saved to a file. The filename is [slug].gpkg, where slug is a string id made out of city_name
+    city_id : str | None, default None
+        If set, the slugified city_id is used in the filename of the data export. For example, a city_id "Athens" will slugify into "athens" in filenames. If set to None, the slugified city_query is used in the filename of the data export. It is useful to set a city_id for cities where the city_query is not the city name, for example to set for a city_query "Municipality of Athens" the city_id to "Athens".
     export_file_format : str, optional, default "geojson"
         File format for the data export, relevant if export_data set to True. Default "geojson", also possible "gpkg". If exporting as geojson, generates extra files for street network and city boundary. If exporting as gkpg, these are added all in one file as extra layers.
     export_plot : bool, optional, default False
         If set to True, plot will be saved to a file
+    import_files: dict, default {}
+        The following key:value entries can be set:
+            "street_network" : str | None, default None
+                If not set to None, the street network is loaded from this file. Must be a gpkg file in unprojected crs EPSG:4326 with layers nodes and edges, with the structure that an undirected osmnx street network g has after saved via ox.io.save_graph_geopackage(). For example:
+                >>> ox.settings.useful_tags_way = ["highway", "cycleway", "cycleway:right", "cycleway:left", "cycleway:both", "cyclestreet"]
+                >>> g = ox.graph_from_place("Barcelona", network_type='all', simplify=False)
+                >>> g = nx.MultiGraph(ox.convert.to_digraph(g))
+                >>> ox.io.save_graph_geopackage(g, "Barcelona_streets.gpkg").
     Returns
     -------
     gdf : geopandas.geodataframe.GeoDataFrame
@@ -46,8 +59,8 @@ def fixbikenet(
     [1] Vybornova, A., Cunha, T., Gühnemann, A. and Szell, M. (2023), Automated Detection of Missing Links in Bicycle Networks. Geogr Anal, 55: 239-267. https://doi.org/10.1111/gean.12324
     """
     # check if user input is valid
-    if type(city_name) != str:
-        raise TypeError("city_name must be a string")
+    if type(city_query) != str:
+        raise TypeError("city_query must be a string")
     if type(proj_crs) != str:
         raise TypeError("proj_crs must be a string")
     if type(radius) != int:
@@ -58,38 +71,33 @@ def fixbikenet(
         raise TypeError("export_data must be a boolean")
     if export_file_format != "geojson" and export_file_format != "gpkg":
         raise ValueError("export_file_format must be 'geojson' or 'gpkg'")
+    if type(import_files) is not dict:
+        raise TypeError("import_files must be a dictionary")
+        # Prepare special case import_files. Turn it into a defaultdict where missing keys are None.
+    import_files = defaultdict(lambda: None, import_files)
 
-    ### downloading and preprocessing data from OSM
-    print("Downloading OSM data..")
+    if import_files['street_network'] is not None:
+        print("Importing street network..")
+        g = import_network(import_files['street_network'])
 
-    ox.settings.useful_tags_way = ["highway", "cycleway", "cycleway:right", "cycleway:left", "cycleway:both", "cyclestreet"]
+    else:
+        ### downloading and preprocessing data from OSM
+        print("Downloading OSM data..")
+
+        ox.settings.useful_tags_way = ["highway", "cycleway", "cycleway:right", "cycleway:left", "cycleway:both", "cyclestreet"]
+
+        # fetch street network data from osmnx
+        g = ox.graph_from_place(
+            city_query, network_type='all', simplify=False
+        )
 
-    # fetch street network data from osmnx
-    g = ox.graph_from_place(
-        city_name, network_type='all', simplify=False
-    )
     g = ox.simplify_graph(
         g,
         edge_attrs_differ=['cycleway', 'highway', 'cycleway:right', 'cycleway:left', 'cycleway:both']
     )
 
-    # export osmnx data to gdfs
-    nodes_gdf, edges_gdf = ox.graph_to_gdfs(
-        g,
-        nodes=True,
-        edges=True,
-        node_geometry=True,
-        fill_edge_geometry=True
-    )
-
-    # project to proj_crs
-    nodes_gdf = nodes_gdf.to_crs(proj_crs)
-    edges_gdf = edges_gdf.to_crs(proj_crs)
-
-
     # check which edges have existing bike infrastructure as defined in config/config_osm.yml and assign boolean value to edges
     g = map_edges_to_bike_infrastructure(g)
-    edges_gdf = bike_infra_mapping_gdf(g, edges_gdf)
 
     print("Dropping parallel edges..")
     edges_to_drop = find_edges_to_drop(g)
@@ -156,34 +164,44 @@ def fixbikenet(
     # add actual geometries in network to each gap
     gdf = create_gdf_with_geoms(gap_df, edges_gdf)
 
+    gdf['ordering'] = gdf.index
+    gdf['length'] = gdf['geometry'].length
+
     edges_pbi_gdf = edges_gdf[edges_gdf["pbi"] == 1]
 
+    # Back to unprojected (potentially). No more calculations after here.
+    gdf.to_crs(epsg=4326, inplace=True)
+    edges_pbi_gdf.to_crs(epsg=4326, inplace=True)
+    edges_gdf.to_crs(epsg=4326, inplace=True)
+
     # Generate export data filename
     if export_data:
-        os.makedirs("./results/", exist_ok=True)
+        os.makedirs(settings.export_path, exist_ok=True)
+        if city_id is None:
+            city_string = city_query
+        else:
+            city_string = city_id
         export_data_filename = (
-                city_name + "." + export_file_format
+                slugify(city_string) + "-fixbikenet-gaps" + "." + export_file_format
         )
 
     if export_data:
         ### save data
         print("Saving data..")
         edges_pbi_gdf.drop(["osmid"], axis=1, inplace=True)
-        city_boundary = ox.geocoder.geocode_to_gdf(city_name)
-        city_boundary.to_crs(epsg=proj_crs, inplace=True)
-        # We have meter precision, so rounding to integers is fine. Better would be to
-        # change dtypes to int, but this does not seem possible without manual looping.
-        city_boundary.geometry = city_boundary.geometry.set_precision(grid_size=1)
-        edges_pbi_gdf.geometry = edges_pbi_gdf.geometry.set_precision(grid_size=1)
-        gdf.geometry = gdf.geometry.set_precision(grid_size=1)
+        edges_gdf.drop(["osmid"], axis=1, inplace=True)
+        city_boundary = ox.geocoder.geocode_to_gdf(city_query)
+        city_boundary.to_crs(epsg=4326, inplace=True)
         if export_file_format == "geojson":
-            gdf.to_file("./results/"+export_data_filename+".geojson", driver="GeoJSON")
-            edges_pbi_gdf.to_file("./results/"+city_name+"-existing_bike_network.geojson", driver="GeoJSON")
-            city_boundary.to_file("./results/"+city_name+"-city_boundary.geojson", driver="GeoJSON")
+            gdf.to_file(settings.export_path + export_data_filename, driver="GeoJSON", RFC7946="YES")
+            edges_pbi_gdf.to_file(settings.export_path + slugify(city_string) + "-fixbikenet" +  "-existing_bike_network.geojson", driver="GeoJSON", RFC7946="YES")
+            edges_gdf.to_file(settings.export_path + slugify(city_string) + "-fixbikenet" + "-existing_street_network.geojson", driver="GeoJSON", RFC7946="YES")
+            city_boundary.to_file(settings.export_path + slugify(city_string) + "-city_boundary.geojson", driver="GeoJSON", RFC7946="YES")
         elif export_file_format == "gpkg":
-            gdf.to_file("./results/"+export_data_filename, driver="GPKG", layer="Identified gaps")
-            edges_pbi_gdf.to_file("./results/"+export_data_filename, driver="GPKG", layer="Existing bike network", append=True)
-            city_boundary.to_file("./results/"+export_data_filename, driver="GPKG", layer="City boundary", append=True)
+            gdf.to_file(settings.export_path + export_data_filename, driver="GPKG", layer="Identified gaps")
+            edges_pbi_gdf.to_file(settings.export_path + export_data_filename, driver="GPKG", layer="Existing bike network", append=True)
+            edges_gdf.to_file(settings.export_path + export_data_filename, driver="GPKG", layer="Existing street network", append=True)
+            city_boundary.to_file(settings.export_path + export_data_filename, driver="GPKG", layer="City boundary", append=True)
 
         if export_plot:
             print("Saving plot..")

fixbikenet/functions.py

@@ -3,9 +3,54 @@
 import numpy as np
 import pandas as pd
 import geopandas as gpd
+import osmnx as ox
+import re
 import itertools
 from shapely.geometry import Point, LineString
 from . import config
+from . import settings
+
+def import_network(street_network, import_path=settings.import_path):
+    """Import and project a street network from gpkg file
+
+    For all edges between a pair of nodes u and v there must be one edge with key 0.
+
+    Parameters
+    ----------
+    street_network : str
+        The street network will be loaded from this file. Must be a gpkg file in unprojected crs EPSG:4326 with layers nodes and edges, with the structure that a osmnx street network g has after saving its undirected version via ox.io.save_graph_geopackage(). For example:
+        >>> g = ox.graph_from_place("Barcelona", network_type='all')
+        >>> g = nx.MultiGraph(ox.convert.to_digraph(g))
+        >>> ox.io.save_graph_geopackage(g, "Barcelona_streets.gpkg")
+    import_path : str, default settings.import_path
+        Path to import files.
+
+    Returns
+    -------
+    nodes : geopandas.geodataframe.GeoDataFrame
+        Extracted OSM nodes, projected
+    edges : geopandas.geodataframe.GeoDataFrame
+        Extracted OSM edges, projected
+    g_undir : networkx.classes.multigraph.MultiGraph
+        Extracted networkX graph, undirected
+    city_boundary_gdf : geopandas.geodataframe.GeoDataFrame
+        Convex hull of the street network
+    """
+
+    nodes = gpd.read_file(import_path+street_network, layer='nodes')
+    edges = gpd.read_file(import_path+street_network, layer='edges')
+
+    # Set indices as required by osmnx.convert.graph_from_gdfs
+    # See: https://osmnx.readthedocs.io/en/stable/user-reference.html#osmnx.utils_graph.graph_from_gdfs
+    nodes = nodes.set_index(['osmid'])
+    edges = edges.set_index(['u', 'v', 'key'])
+
+    g = ox.convert.graph_from_gdfs(nodes, edges)
+
+    #city_boundary_gdf = gpd.GeoDataFrame(gpd.GeoSeries(nodes.union_all().convex_hull), geometry=0, crs=nodes.crs) # We do this before the projection of nodes below
+    # To do: To be super-correct, the hull should be buffered by settings.seed_point_snap_distance (in degrees due to being unprojected)
+
+    return g
 
 def map_edges_to_bike_infrastructure(g):
     """
@@ -298,7 +343,8 @@ def rank_gaps_by_b(found_gaps_nsp, G, ebc):
     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])
+        #ebcs = np.array([ebc[edge] for edge in edgelist])
+        ebcs = np.array([ebc.get(edge, ebc.get(edge[::-1])) for edge in edgelist])
         B = sum(lengths * ebcs) / sum(lengths)
         Bs.append(B)
     return Bs
@@ -439,7 +485,8 @@ def compute_benefit_metric(comp, node_path, ebc):
     """
     edgelist = [tuple(sorted(z)) for z in zip(node_path, node_path[1:])]
     lengths = np.array([comp.edges[edge]["length"] for edge in edgelist])
-    ebcs = np.array([ebc[edge] for edge in edgelist])
+    #ebcs = np.array([ebc[edge] for edge in edgelist])
+    ebcs = np.array([ebc.get(edge, ebc.get(edge[::-1])) for edge in edgelist])
     B = sum(lengths * ebcs) / sum(lengths)
     return B
 
@@ -542,4 +589,31 @@ def gap_declustering(gaps_df, G, ebc):
             "benefit": selected_scores,
         }
     )
-    return result
+    return result
+
+def slugify(s):
+    """Slugify a string
+
+    Source: https://github.com/Chalarangelo/30-seconds-of-code/blob/master/content/snippets/python/s/slugify.md
+    Note: A clean global solution would be using unidecode, but we do not want extra dependencies for this. We assume European city names in latin alphabet, some special letters like Hungarian long ö already mapped.
+
+    Parameters
+    ----------
+    s : str
+        String to slufigy
+
+    Returns
+    -------
+    s : str
+        Slugified string
+    """
+    s = s.lower().strip()
+    s = re.sub(r'[\s-]+', '', s)  # Remove white spaces, -
+    s = re.sub(r'[^\w\s-]', '', s)
+    s = re.sub(r'^-+|-+$', '', s)
+    tab = str.maketrans(
+        "áéíóúàèìùòâêîôûäëïöüǎěǐǒǔãẽĩõũăåæçčıłñňøœřßșşšůŷÿźž",
+        "aeiouaeiouaeiouaeiouaeiouaeiouaaaccilnnoorssssuyyzz"
+    )
+    s = s.translate(tab)
+    return s

fixbikenet/settings.py

@@ -0,0 +1,5 @@
+"""Global settings for fixbikenet that can be configured by the user."""
+
+import_path = "./"
+export_path = "./results/"
+crs_projected = '3857'