From 7d4416274a7b11909227528ee9310f3e2b65d744 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Tue, 28 Apr 2026 11:29:32 +0200
Subject: [PATCH 1/8] Copy-paste changes from folcoexistingnw branch

Does not run yet, needs fixes
---
 examples/mwe.py            |   3 +-
 growbikenet/functions.py   | 108 +++++++++++++++++++++++++-
 growbikenet/growbikenet.py | 152 ++++++++++++++++++++++---------------
 3 files changed, 200 insertions(+), 63 deletions(-)

diff --git a/examples/mwe.py b/examples/mwe.py
index 3c4ecf6..a9f6c39 100644
--- a/examples/mwe.py
+++ b/examples/mwe.py
@@ -3,9 +3,10 @@
 import growbikenet as gbn
 
 a_edges = gbn.growbikenet(
-    city_name="Bath",
+    city_name="Turin",
     proj_crs="3857",
     ranking="betweenness_centrality",
+    existing_network_spacing=500,
     export_data=True,
     export_plots=False,
     export_video=False,
diff --git a/growbikenet/functions.py b/growbikenet/functions.py
index 7fe3598..55c2f30 100644
--- a/growbikenet/functions.py
+++ b/growbikenet/functions.py
@@ -2,6 +2,7 @@
 import pandas as pd
 import geopandas as gpd
 import networkx as nx
+import osmnx as ox
 from scipy.spatial import Delaunay
 from shapely.prepared import prep
 from shapely.geometry import Point
@@ -27,9 +28,77 @@ def intersects_properly(geom1, geom2):
     return geom1.intersects(geom2) and not geom1.touches(geom2)
 
 
+def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None):
+    """Download and prepare a street network from OSM via OSMnx
+    Downloads a network with a given network_type and custom_filter using ox.graph_from_place.
+    Then, stores the undirected OSM data in gdfs and projects using proj_crs.
+    Parameters
+    ----------
+    city_name : 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).
+    network_type : {“all”, “all_public”, “bike”, “drive”, “drive_service”, “walk”} 
+        What type of street network to retrieve if custom_filter is None.
+    custom_filter : (str | list[str] | None)
+        A custom ways filter to be used instead of the network_type presets
+    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
+    """
+    # Fetch street network data from osmnx
+    g = ox.graph_from_place(
+    city_name, network_type=network_type, custom_filter=custom_filter, retain_all=True
+    )
+    g_undir = g.to_undirected().copy() # convert to undirected (dropping OSMnx keys!)
+
+    # Export osmnx data to gdfs
+    nodes, edges = nx_to_nodes_edges(g_undir, proj_crs)
+    return nodes, edges, g_undir
+
+def nx_to_nodes_edges(G, proj_crs='3857'):
+    """Get nodes and projected edges from networkX graph
+    Parameters
+    ----------
+    G : networkx.classes.multigraph.MultiGraph
+        networkX graph, undirected
+    proj_crs : str, default '3857'
+        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+    Returns
+    -------
+    nodes : geopandas.geodataframe.GeoDataFrame
+        Extracted OSM nodes, projected, osmid is index
+    edges : geopandas.geodataframe.GeoDataFrame
+        Extracted OSM edges, projected
+    """
+    nodes, edges = ox.graph_to_gdfs(
+    G,
+    nodes=True,
+    edges=True,
+    node_geometry=True,
+    fill_edge_geometry=True
+    )
+
+    # Replace after dropping edges with key = 1
+    edges = edges.loc[:,:,0].copy()
+    # This also means we are dropping the "key" level from edge index (u,v,key becomes: u,v)
+
+    # Project geometries of nodes, edges, seed points
+    edges = edges.to_crs(proj_crs)
+    nodes = nodes.to_crs(proj_crs)
+
+    # Add osm ID as column to node gdf
+    nodes["osmid"] = nodes.index
+    return nodes, edges
+    
 def get_correct_edgetuples(edge_gdf, nodelist):
     """
-    helper function that maps a node list (output of nx.shortest_paths)
+    Helper function that maps a node list (output of nx.shortest_paths)
     to the correct set of edge tuples that can be used for INDEXING THE EDGE GDF
 
     Parameters
@@ -54,6 +123,43 @@ def get_correct_edgetuples(edge_gdf, nodelist):
     return edgelist_final
 
 
+def get_existing_network_seed_points(nodes_exnw, existing_network_spacing):
+    """Get seed points on an existing bicycle network
+    
+    Start with the first (arbitrary) node from nodes_exnw. Then, for each node: Delete all other nodes closer than existing_network_spacing, proceed with the closest of the remaining nodes. Finish once all nodes are found or deleted.
+    
+    Parameters
+    ----------
+    nodes_exnw: geopandas.geodataframe.GeoDataFrame
+        Nodes of the existing bicycle network, in a projected coordinate reference system.
+    existing_network_spacing: int
+        Distance between seed points, in meters.
+    Returns
+    -------
+    seed_points_exnw: geopandas.geodataframe.GeoDataFrame
+        Seed points, already part of the network, in the same projected coordinate reference system as edges
+    """
+    # Start with the first (arbitrary) node from nodes_exnw
+    node_current = nodes_exnw.iloc[[0]]
+
+    seed_points_exnw = gpd.GeoDataFrame()
+    while len(node_current)>0 and len(nodes_exnw)>1:
+        # Find all too close nodes to the current nodes
+        nodes_too_close = nodes_exnw.loc[(nodes_exnw.geometry.distance(Point(node_current.iloc[0].geometry)) <= existing_network_spacing)]
+        nodes_too_close = nodes_too_close.iloc[:, :-1] # osmid is there twice now (once in the end), so it needs to be dropped
+
+        # Delete the nodes that are too close to nodes_exnw
+        nodes_exnw = nodes_exnw.overlay(nodes_too_close, how='difference')
+
+        # Add current node to seed_points_exnw
+        seed_points_exnw = pd.concat([seed_points_exnw, node_current], ignore_index=True)
+
+        # Find the node in nodes_exnw that is closest to the existing seed points
+        node_current = seed_points_exnw.sjoin_nearest(nodes_exnw, how="inner") 
+        node_current = nodes_exnw[nodes_exnw.osmid == node_current["osmid_right"].values[0]]
+
+    return seed_points_exnw
+    
 def get_grid_seed_points(edges, seed_point_spacing, principal_bearing):
     """Get grid seed points for street network, rotated by principal bearing
 
diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index f1ec7d1..1ab4255 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -11,6 +11,7 @@ def growbikenet(
     seed_point_type="grid",
     seed_point_grid_spacing=1707,
     seed_point_delta=500,
+    existing_network_spacing=None,
     export_data=True,
     export_data_slug=None,
     export_plots=False,
@@ -18,33 +19,33 @@ def growbikenet(
 ):
     """Creates a list of edges ordered by a specified ranking method.
 
-    The edges form a subnetwork of a city's street network, interpreted as a growing bicycle network following [1]_.
-    Note that the original paper [1]_ uses minimum weight triangulation, but Delaunay triangulation is much faster due to the Delaunay scipy function and gives in most cases identical results.
-    Triangulation is calculated for the abstract network, but metrics (betweenness, closeness) are calculated for the routed network accounting for lengths.
+    The edges form a subnetwork of a city's street network, interpreted as a growing bicycle network following [1]_. By default, growth is from scratch, but the existing bicycle network can also be used as a starting point[2]_. Note that the original paper [1]_ uses minimum weight triangulation, but Delaunay triangulation is much faster due to the Delaunay scipy function and gives in most cases identical results. Triangulation is calculated for the abstract network, but metrics (betweenness, closeness) are calculated for the routed network accounting for lengths.
 
     Parameters
     ----------
     city_name : str
-        name of the city that the analysis should be performed on
+        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)
+        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator)
     ranking : str, default 'betweenness_centrality'
         Method used to rank edges. Must be 'betweenness_centrality' (default), 'closeness_centrality', or 'all'. If 'all', will also add a random ranking.
     seed_point_type : str, optional, default 'grid'
-        if set to 'grid', creates a square grid
-        if set to 'rail', uses rail stations
+        If set to 'grid', creates a square grid
+        If set to 'rail', uses rail stations
     seed_point_grid_spacing : int, optional, default 1707
-        if seed_point_type is set to 'grid', this is the spacing between seed points, in meters
+        If seed_point_type is set to 'grid', this is the spacing between seed points, in meters
     seed_point_delta : int, optional, default 500
-        maximum distance between generated seed points and osm nodes for snapping
+        Maximum distance between generated seed points and osm nodes for snapping
+    existing_network_spacing : int, optional, default None
+        Spacing between seed points, in meters, only on the existing bicycle network. If not set to a positive integer, the existing network is ignored.
     export_data : bool, optional, default True
-        if set to True, data will be saved to a file. The filename is [slug]-[ranking]-[seed_point_type].gpkg, where slug is a string id made out of city_name
+        If set to True, data will be saved to a file. The filename is [slug]-[ranking]-[seed_point_type].gpkg, where slug is a string id made out of city_name
     export_data_slug : string, optional, default None
-        if not set to None, it will be slugified and used as the slug in the filename of the data export
+        If not set to None, it will be slugified and used as the slug in the filename of the data export
     export_plots : bool, optional, default False
-        if set to True, plots will be saved to a file
+        If set to True, plots will be saved to a file
     export_video : bool, optional, default False
-        if set to True, video will be saved to a file (only possible if export_plots is set to True)
+        If set to True, video will be saved to a file (only possible if export_plots is set to True)
 
     Returns
     -------
@@ -54,9 +55,11 @@ def growbikenet(
     References
     ----------
     .. [1] M. Szell, S. Mimar, T. Perlman, G. Ghoshal, R. Sinatra, "Growing urban bicycle networks", Scientific Reports 12, 6765 (2022)
+    .. [2] P. Folco, L. Gauvin, M. Tizzoni, M. Szell, "Data-driven micromobility network planning for demand and safety", Environment and planning B: Urban analytics and city science 50(8), 2087-2102 (2023)
 
     """
-    # check if user input is valid
+    
+    # Check if user input is valid
     if type(city_name) != str:
         raise TypeError("city_name must be a string")
     if type(proj_crs) != str:
@@ -71,8 +74,16 @@ def growbikenet(
         raise ValueError("seed_point_type must be 'grid' or 'rail'")
     if seed_point_type == "grid" and type(seed_point_grid_spacing) != int:
         raise TypeError("seed_point_grid_spacing must be an integer")
+    if seed_point_type == 'grid' and type(seed_point_grid_spacing) == int and seed_point_grid_spacing <= 0:  
+        raise ValueError("seed_point_grid_spacing must be a positive integer")
     if type(seed_point_delta) != int:
         raise TypeError("seed_point_delta must be an integer")
+    if type(seed_point_delta) == int and seed_point_delta <= 0:
+        raise ValueError("seed_point_delta must be a positive integer")
+    if type(existing_network_spacing) != int and existing_network_spacing is not None:
+        raise TypeError("existing_network_spacing must be None or a positive integer")
+    if type(existing_network_spacing) == int and existing_network_spacing <= 0:
+        raise ValueError("existing_network_spacing must be None or a positive integer")
     if type(export_data) != bool:
         raise TypeError("export_data must be a boolean")
     if export_data_slug is not None and type(export_data_slug) != str:
@@ -90,34 +101,29 @@ def growbikenet(
 
     np.random.seed(42)  # Set random number generator seed for reproducibility
 
-    ### downloading and preprocessing data from OSM
+    ### Download and preprocess data from OSM
     print("Downloading OSM data..")
 
-    # fetch street network data from osmnx
-    g = ox.graph_from_place(city_name, network_type="all")
-    g_undir = g.to_undirected().copy()  # convert to undirected (dropping OSMnx keys!)
-
-    # export osmnx data to gdfs
-    nodes, edges = ox.graph_to_gdfs(
-        g_undir, nodes=True, edges=True, node_geometry=True, fill_edge_geometry=True
-    )
-
-    # # save "original" graph data (in orig_crs)
-    # nodes.to_file("nodes.gpkg", driver='GPKG')
-    # edges.to_file("edges.gpkg", driver='GPKG')
-
-    # replace after dropping edges with key = 1
-    edges = edges.loc[:, :, 0].copy()
-    # this also means we are dropping the "key" level from edge index (u,v,key becomes: u,v)
-
-    # project geometries of nodes, edges, seed points
-    edges = edges.to_crs(proj_crs)
-    nodes = nodes.to_crs(proj_crs)
-
-    # add osm ID as column to node gdf
-    nodes["osmid"] = nodes.index
-
-    ### creating seed points
+    # Fetch street network data from osmnx
+    nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public')
+
+    if existing_network_spacing:
+        cf = ['["cycleway"~"track"]',
+              '["highway"~"cycleway"]',
+              '["highway"~"path"]["bicycle"~"designated"]',
+              '["cycleway:right"~"track"]',
+              '["cycleway:left"~"track"]',
+              '["cyclestreet"]',
+              '["highway"~"living_street"]'
+             ]
+        for custom_tag in ["cycleway", "bicycle", "cycleway:right", "cycleway:left", "cyclestreet"]:
+            if custom_tag not in ox.settings.useful_tags_way:
+                ox.settings.useful_tags_way.extend(custom_tag)
+        nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf)
+        g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
+        _, edges = nx_to_nodes_edges(g_undir, proj_crs)
+
+    ### Create seed points
     print("Creating " + seed_point_type + " seed points..")
 
     if seed_point_type == "grid":
@@ -140,25 +146,50 @@ def growbikenet(
     seed_points_snapped = snap_seed_points(seed_points, nodes)
     seed_points_snapped = filter_seed_points(seed_points_snapped, seed_point_delta)
 
+    if existing_network_spacing:
+        # If the existing bicycle network is used, create extra seed points on it. They are by construction already snapped.
+        seed_points_exnw = get_existing_network_seed_points(nodes_exnw, existing_network_spacing)
+        seed_points_exnw.to_crs(edges.crs, inplace=True)
+
+        # Afterwards, drop all previously determined seed points (grid or rail) that are now too close to these extra points.
+        buffer_seed_points_exnw = gpd.GeoDataFrame(seed_points_exnw.buffer(existing_network_spacing))
+        buffer_seed_points_exnw = buffer_seed_points_exnw.rename(columns={0:'geometry'}).set_geometry('geometry') # https://gis.stackexchange.com/questions/266098/how-to-convert-a-geoseries-to-a-geodataframe-with-geopandas
+        buffer_seed_points_exnw.to_crs(edges.crs, inplace=True)
+
+        # Delete the seed points that are too close to seed_points_exnw via its buffer
+        seed_points_snapped = seed_points_snapped.overlay(buffer_seed_points_exnw, how='difference')
+
+        # Merge original snapped points with new existing network points (=already snapped)
+        seed_points_snapped = seed_points_snapped.overlay(seed_points_exnw, how='union')
+        # seed_points_snapped.to_file(seed_points_snapped, driver="GPKG")
+
+        # Bring back to original form (pandas df, columns, osmid index)
+        # This is a bit of a mess but it works. Simplify it in the future.
+        seed_points_snapped = pd.DataFrame(seed_points_snapped)
+        seed_points_snapped.loc[seed_points_snapped['osmid_1'].isnull(), 'osmid_1'] = seed_points_snapped['osmid_2'] # _1 comes from one side, _2 from the other. One has NaNs, the other too. https://stackoverflow.com/a/60132614
+        seed_points_snapped.drop(["y","x","street_count", "highway", "osmid_2"], axis=1, inplace=True)
+        seed_points_snapped.rename(columns={"osmid_1": "osmid"}, inplace=True)
+        seed_points_snapped.set_index("osmid", drop=False, inplace=True)
+        
     # Abort if only 0 or 1 seed points
     if len(seed_points_snapped) < 2:
         raise RuntimeError("Found less than 2 seed points")
 
-    ### running greedy triangulation
+    ### Triangulate
     # Triangulation is calculated for the abstract network, but metrics (betweenness, closeness) are calculated for the routed network accounting for lengths.
     print("Greedy triangulation..")
 
-    # create df with delaunay edges
+    # Create df with delaunay edges
     df = create_delaunay_edges(seed_points_snapped)
 
-    # map each abstract edge to a merged geometry of corresponding osmnx edges (routed on g_undir)
+    # Map each abstract edge to a merged geometry of corresponding osmnx edges (routed on g_undir)
     df = add_path_to_df(df, edges, g_undir)
 
-    # get "routed" geometry (LineString) for each abstract edge (row)
+    # Get "routed" geometry (LineString) for each abstract edge (row)
     print("Routing..")
     gdf = create_gdf_with_geoms(df, edges)
 
-    # add distances between source and target from geometry
+    # Add distances between source and target from geometry
     gdf["dist"] = gdf["geometry"].length
 
     edge_list = gdf["pair"]
@@ -166,40 +197,39 @@ def growbikenet(
     dist_dict = dict(zip(edge_list, dist_list))
     geom_dict = dict(zip(edge_list, gdf["geometry"].tolist()))
 
-    # make graph object from edge list
+    # Make graph object from edge list
     A = nx.Graph()
     A.add_nodes_from(seed_points_snapped.index)
     A.add_edges_from(edge_list)
     nx.set_edge_attributes(A, dist_dict, "distance")
     nx.set_edge_attributes(A, geom_dict, "geometry")
 
-    ### compute edge attributes
+    ### Compute edge attributes
     print("Computing edge attributes..")
-    # metric_dict =
     if ranking == "betweenness_centrality" or ranking == "all":
-        # add betweenness attributes to edges
+        # Add betweenness attributes to edges
         bc_values = nx.edge_betweenness_centrality(
             A, weight="distance", normalized=True
         )
         nx.set_edge_attributes(A, bc_values, name="betweenness_centrality")
     if ranking == "closeness_centrality" or ranking == "all":
-        # add closeness attributes to nodes and edges
+        # Add closeness attributes to nodes and edges
         cc_values_nodes = nx.closeness_centrality(A, distance="distance")
         nx.set_node_attributes(A, cc_values_nodes, name="closeness_centrality")
         cc_values = node_to_edge_attributes(cc_values_nodes, A.edges)
         nx.set_edge_attributes(A, cc_values, name="closeness_centrality")
 
-    # export attributes to gdfs:
+    ### Export attributes to gdfs:
 
-    # create dataframe and add method as edge attribute
+    # Create dataframe and add method as edge attribute
     a_edges = df_from_graph(A, ranking)
 
-    # rank edges by specified method
+    # Rank edges by specified method
     a_edges = rank_df(a_edges, ranking)
 
     a_edges = gpd.GeoDataFrame(a_edges, crs=edges.crs, geometry="geometry")
 
-    # generate export data filename
+    # Generate export data filename
     if export_data or export_plots or export_video:
         if export_data_slug is None:
             city_string = city_name
@@ -209,32 +239,32 @@ def growbikenet(
             slugify(city_string) + "-" + ranking + "-" + seed_point_type + ".gpkg"
         )
 
-    # save to file
+    # Save to file
     if export_data:
         ### save data
         print("Saving data..")
         a_edges.to_file(export_data_filename, driver="GPKG")
 
     if export_plots or export_video:
-        ### Visualization
+        ### Visualize
         print("Creating visualizations..")
 
-        # create directories
+        # Create directories
         os.makedirs("./results/", exist_ok=True)
         os.makedirs("./results/plots/", exist_ok=True)
         os.makedirs("./results/plots/video/", exist_ok=True)
 
-        # read in file to plot
+        # Read in file to plot
         routed_edges_gdf = gpd.read_file(export_data_filename)
 
-        # viz/plot settings (move to config file later)
+        # Viz/plot settings (move to config file later)
 
-        # define color palette (from Michael's project: https://github.com/mszell/bikenwgrowth/blob/main/parameters/parameters.py)
+        # Define color palette (from Michael's project: https://github.com/mszell/bikenwgrowth/blob/main/parameters/parameters.py)
         streetcolor = "#999999"
         edgecolor = "#0EB6D2"
         seedcolor = "#ff7338"
 
-        # define linewidths
+        # Define linewidths
 
         lws = {"street": 0.75, "bike": 2}
 

From 144420f706ffe33c50dee8c83e29e9c3c1802985 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Tue, 28 Apr 2026 20:52:34 +0200
Subject: [PATCH 2/8] Fix seed_points_snapped, start work on LCC

---
 examples/mwe.py            |  2 +-
 growbikenet/functions.py   | 15 ++++++++++-----
 growbikenet/growbikenet.py | 28 ++++++++++++++++++++--------
 3 files changed, 31 insertions(+), 14 deletions(-)

diff --git a/examples/mwe.py b/examples/mwe.py
index a9f6c39..759acc0 100644
--- a/examples/mwe.py
+++ b/examples/mwe.py
@@ -6,7 +6,7 @@
     city_name="Turin",
     proj_crs="3857",
     ranking="betweenness_centrality",
-    existing_network_spacing=500,
+    existing_network_spacing=200,
     export_data=True,
     export_plots=False,
     export_video=False,
diff --git a/growbikenet/functions.py b/growbikenet/functions.py
index 55c2f30..2ed1b4f 100644
--- a/growbikenet/functions.py
+++ b/growbikenet/functions.py
@@ -10,7 +10,7 @@
 
 def intersects_properly(geom1, geom2):
     """
-    helper function to check whether newly to be added edge intersects with already added edges
+    Helper function to check whether newly to be added edge intersects with already added edges
     for 2 shapely geometries, check whether they "properly intersect" (i.e. intersect but not touch, i.e. don't share endpoints)
 
     Parameters
@@ -28,7 +28,7 @@ def intersects_properly(geom1, geom2):
     return geom1.intersects(geom2) and not geom1.touches(geom2)
 
 
-def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None):
+def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None, retain_all=True):
     """Download and prepare a street network from OSM via OSMnx
     Downloads a network with a given network_type and custom_filter using ox.graph_from_place.
     Then, stores the undirected OSM data in gdfs and projects using proj_crs.
@@ -42,6 +42,8 @@ def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None)
         What type of street network to retrieve if custom_filter is None.
     custom_filter : (str | list[str] | None)
         A custom ways filter to be used instead of the network_type presets
+    retain_all : bool, default True
+        If True, return the entire graph even if it is not connected, useful for disconnected bicycle networks. If False, retain only the largest weakly connected component, useful for road networks.
     Returns
     -------
     nodes : geopandas.geodataframe.GeoDataFrame
@@ -53,7 +55,7 @@ def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None)
     """
     # Fetch street network data from osmnx
     g = ox.graph_from_place(
-    city_name, network_type=network_type, custom_filter=custom_filter, retain_all=True
+    city_name, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
     )
     g_undir = g.to_undirected().copy() # convert to undirected (dropping OSMnx keys!)
 
@@ -63,12 +65,14 @@ def prepare_network(city_name, proj_crs, network_type='all', custom_filter=None)
 
 def nx_to_nodes_edges(G, proj_crs='3857'):
     """Get nodes and projected edges from networkX graph
+    
     Parameters
     ----------
     G : networkx.classes.multigraph.MultiGraph
         networkX graph, undirected
     proj_crs : str, default '3857'
         Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+        
     Returns
     -------
     nodes : geopandas.geodataframe.GeoDataFrame
@@ -143,7 +147,7 @@ def get_existing_network_seed_points(nodes_exnw, existing_network_spacing):
     node_current = nodes_exnw.iloc[[0]]
 
     seed_points_exnw = gpd.GeoDataFrame()
-    while len(node_current)>0 and len(nodes_exnw)>1:
+    while len(node_current)>0 and len(nodes_exnw)>0:
         # Find all too close nodes to the current nodes
         nodes_too_close = nodes_exnw.loc[(nodes_exnw.geometry.distance(Point(node_current.iloc[0].geometry)) <= existing_network_spacing)]
         nodes_too_close = nodes_too_close.iloc[:, :-1] # osmid is there twice now (once in the end), so it needs to be dropped
@@ -156,7 +160,8 @@ def get_existing_network_seed_points(nodes_exnw, existing_network_spacing):
 
         # Find the node in nodes_exnw that is closest to the existing seed points
         node_current = seed_points_exnw.sjoin_nearest(nodes_exnw, how="inner") 
-        node_current = nodes_exnw[nodes_exnw.osmid == node_current["osmid_right"].values[0]]
+        if len(node_current)>0: # Current nodes could already be depleted here. Then loop will stop.
+            node_current = nodes_exnw[nodes_exnw.osmid == node_current["osmid_right"].values[0]]
 
     return seed_points_exnw
     
diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index 1ab4255..47fa0ec 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -105,7 +105,7 @@ def growbikenet(
     print("Downloading OSM data..")
 
     # Fetch street network data from osmnx
-    nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public')
+    nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public', retain_all=False)
 
     if existing_network_spacing:
         cf = ['["cycleway"~"track"]',
@@ -119,8 +119,16 @@ def growbikenet(
         for custom_tag in ["cycleway", "bicycle", "cycleway:right", "cycleway:left", "cyclestreet"]:
             if custom_tag not in ox.settings.useful_tags_way:
                 ox.settings.useful_tags_way.extend(custom_tag)
-        nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf)
+        nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True)
         g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
+        # Take largest connected component lcc
+        lcc = max(nx.connected_components(g_undir), key=len)
+        g_undir = g_undir.subgraph(lcc).copy() 
+        # TO DO: Restrict nodes and edges of existing net to the lcc
+        # print(g_undir.nodes())
+        # print(nodes_exnw)
+        # print(edges_exnw)
+        # sys.exit()
         _, edges = nx_to_nodes_edges(g_undir, proj_crs)
 
     ### Create seed points
@@ -145,9 +153,11 @@ def growbikenet(
     # Snap seed points to OSM nodes
     seed_points_snapped = snap_seed_points(seed_points, nodes)
     seed_points_snapped = filter_seed_points(seed_points_snapped, seed_point_delta)
-
+    
     if existing_network_spacing:
         # If the existing bicycle network is used, create extra seed points on it. They are by construction already snapped.
+        # print(nodes_exnw)
+        # sys.exit()
         seed_points_exnw = get_existing_network_seed_points(nodes_exnw, existing_network_spacing)
         seed_points_exnw.to_crs(edges.crs, inplace=True)
 
@@ -161,19 +171,21 @@ def growbikenet(
 
         # Merge original snapped points with new existing network points (=already snapped)
         seed_points_snapped = seed_points_snapped.overlay(seed_points_exnw, how='union')
-        # seed_points_snapped.to_file(seed_points_snapped, driver="GPKG")
 
-        # Bring back to original form (pandas df, columns, osmid index)
+        # Bring back to original form (geometry and osmid columns, osmid index)
         # This is a bit of a mess but it works. Simplify it in the future.
-        seed_points_snapped = pd.DataFrame(seed_points_snapped)
         seed_points_snapped.loc[seed_points_snapped['osmid_1'].isnull(), 'osmid_1'] = seed_points_snapped['osmid_2'] # _1 comes from one side, _2 from the other. One has NaNs, the other too. https://stackoverflow.com/a/60132614
-        seed_points_snapped.drop(["y","x","street_count", "highway", "osmid_2"], axis=1, inplace=True)
+        seed_points_snapped.drop(["y","x","street_count", "highway", "railway", "osmid_2"], axis=1, inplace=True)
         seed_points_snapped.rename(columns={"osmid_1": "osmid"}, inplace=True)
         seed_points_snapped.set_index("osmid", drop=False, inplace=True)
         
+        # seed_points_snapped.drop(["osmid"], axis=1, inplace=True)
+        # print(seed_points_snapped)
+        # seed_points_snapped.to_file("test.gpkg", driver="GPKG")
+        
     # Abort if only 0 or 1 seed points
     if len(seed_points_snapped) < 2:
-        raise RuntimeError("Found less than 2 seed points")
+        raise RuntimeError("Found less than 2 seed points, but more are needed.")
 
     ### Triangulate
     # Triangulation is calculated for the abstract network, but metrics (betweenness, closeness) are calculated for the routed network accounting for lengths.

From 2bb372bd93add6f7f276ccb922e70f2c732e02a0 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 08:32:31 +0200
Subject: [PATCH 3/8] Fix ruff style issues with new type checks

---
 growbikenet/growbikenet.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index 47fa0ec..4b152fa 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -80,9 +80,9 @@ def growbikenet(
         raise TypeError("seed_point_delta must be an integer")
     if type(seed_point_delta) == int and seed_point_delta <= 0:
         raise ValueError("seed_point_delta must be a positive integer")
-    if type(existing_network_spacing) != int and existing_network_spacing is not None:
+    if type(existing_network_spacing) is not int and existing_network_spacing is not None:
         raise TypeError("existing_network_spacing must be None or a positive integer")
-    if type(existing_network_spacing) == int and existing_network_spacing <= 0:
+    if type(existing_network_spacing) is int and existing_network_spacing <= 0:
         raise ValueError("existing_network_spacing must be None or a positive integer")
     if type(export_data) != bool:
         raise TypeError("export_data must be a boolean")

From d47ef63f0c6b61d38e0c7cccdd9595970fb0e130 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 09:32:41 +0200
Subject: [PATCH 4/8] Restrict nodes and edges of existing bikenw to lcc

---
 growbikenet/growbikenet.py | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index 4b152fa..ecbfa72 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -105,6 +105,7 @@ def growbikenet(
     print("Downloading OSM data..")
 
     # Fetch street network data from osmnx
+    # Due to retain_all=False, this fetches the largest connected component
     nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public', retain_all=False)
 
     if existing_network_spacing:
@@ -119,16 +120,22 @@ def growbikenet(
         for custom_tag in ["cycleway", "bicycle", "cycleway:right", "cycleway:left", "cyclestreet"]:
             if custom_tag not in ox.settings.useful_tags_way:
                 ox.settings.useful_tags_way.extend(custom_tag)
+        # Fetch protected bike network data from osmnx
+        # Due to retain_all=True, this fetches all the connected components
         nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True)
         g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
-        # Take largest connected component lcc
+
+        # Now we could have some leftover bike infra that is disconnected from the street network and thus not routable.
+        # We delete those parts next:
+        # Take largest connected component lcc of the merged network
         lcc = max(nx.connected_components(g_undir), key=len)
         g_undir = g_undir.subgraph(lcc).copy() 
-        # TO DO: Restrict nodes and edges of existing net to the lcc
-        # print(g_undir.nodes())
-        # print(nodes_exnw)
-        # print(edges_exnw)
-        # sys.exit()
+        # Restrict nodes and edges of the existing bike net to this lcc
+        valid_node_osmids = g_undir.nodes()
+        nodes_exnw = nodes_exnw[nodes_exnw['osmid'].isin(valid_node_osmids)]
+        # edges_exnw has a MultiIndex ('u','v'), so we must use get_level_values, see https://stackoverflow.com/a/18835121
+        edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('u').isin(valid_node_osmids)]
+        edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('v').isin(valid_node_osmids)]
         _, edges = nx_to_nodes_edges(g_undir, proj_crs)
 
     ### Create seed points

From c086e0a53523aa3a30a76cb62a1eb0648b9e76ac Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 11:19:59 +0200
Subject: [PATCH 5/8] Refactor update_with_existing_bike_network

---
 growbikenet/functions.py   | 58 +++++++++++++++++++++++++++
 growbikenet/growbikenet.py | 82 +++++++++++++-------------------------
 2 files changed, 85 insertions(+), 55 deletions(-)

diff --git a/growbikenet/functions.py b/growbikenet/functions.py
index 2ed1b4f..91019b3 100644
--- a/growbikenet/functions.py
+++ b/growbikenet/functions.py
@@ -165,6 +165,64 @@ def get_existing_network_seed_points(nodes_exnw, existing_network_spacing):
 
     return seed_points_exnw
     
+def update_with_existing_bike_network(city_name, proj_crs, g_undir):
+    """Update street network with existing bike network
+
+    Downloads a network of protected bike infrastructure from OSM (retaining all connected components) and merges it to a given street network graph g_undir.
+    
+    Parameters
+    ----------
+    city_name : str
+        Name of the city that the analysis should be performed on.
+    proj_crs : str
+        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+    g_undir : networkx.classes.multigraph.MultiGraph
+        Street network networkX graph, undirected
+
+    Returns
+    -------
+    nodes : geopandas.geodataframe.GeoDataFrame
+        Updated OSM nodes of the street network, projected
+    edges : geopandas.geodataframe.GeoDataFrame
+        Updated OSM edges of the street network, projected
+    g_undir : networkx.classes.multigraph.MultiGraph
+        Updated street networkX graph, undirected
+    nodes_exnw : geopandas.geodataframe.GeoDataFrame
+        OSM nodes of the corresponding bike network, projected
+    edges_exnw : geopandas.geodataframe.GeoDataFrame
+        OSM edges of the corresponding bike network, projected
+    """
+    cf = ['["cycleway"~"track"]',
+            '["highway"~"cycleway"]',
+          '["highway"~"path"]["bicycle"~"designated"]',
+          '["cycleway:right"~"track"]',
+          '["cycleway:left"~"track"]',
+          '["cyclestreet"]',
+          '["highway"~"living_street"]'
+        ]
+    for custom_tag in ["cycleway", "bicycle", "cycleway:right", "cycleway:left", "cyclestreet"]:
+        if custom_tag not in ox.settings.useful_tags_way:
+            ox.settings.useful_tags_way.extend(custom_tag)
+    # Fetch protected bike network data from osmnx
+    # Due to retain_all=True, this fetches all the connected components
+    nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True)
+    g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
+
+    # Now we could have some leftover bike infra that is disconnected from the street network and thus not routable.
+    # We delete those parts next:
+    # Take largest connected component lcc of the merged network
+    lcc = max(nx.connected_components(g_undir), key=len)
+    g_undir = g_undir.subgraph(lcc).copy() 
+    # Restrict nodes and edges of the existing bike net to this lcc
+    valid_node_osmids = g_undir.nodes()
+    nodes_exnw = nodes_exnw[nodes_exnw['osmid'].isin(valid_node_osmids)]
+    # edges_exnw has a MultiIndex ('u','v'), so we must use get_level_values, see https://stackoverflow.com/a/18835121
+    edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('u').isin(valid_node_osmids)]
+    edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('v').isin(valid_node_osmids)]
+    nodes, edges = nx_to_nodes_edges(g_undir, proj_crs)
+
+    return nodes, edges, g_undir, nodes_exnw, edges_exnw
+
 def get_grid_seed_points(edges, seed_point_spacing, principal_bearing):
     """Get grid seed points for street network, rotated by principal bearing
 
diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index ecbfa72..e5f623b 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -2,7 +2,7 @@
 from slugify import slugify
 from growbikenet.functions import *
 from growbikenet.visualizations import *
-
+import warnings
 
 def growbikenet(
     city_name,
@@ -84,6 +84,8 @@ def growbikenet(
         raise TypeError("existing_network_spacing must be None or a positive integer")
     if type(existing_network_spacing) is int and existing_network_spacing <= 0:
         raise ValueError("existing_network_spacing must be None or a positive integer")
+    if type(existing_network_spacing) is int and existing_network_spacing >= seed_point_grid_spacing:
+        warnings.warn("existing_network_spacing is recommended to be smaller than seed_point_grid_spacing to ensure that the existing bicycle network is built first.")
     if type(export_data) != bool:
         raise TypeError("export_data must be a boolean")
     if export_data_slug is not None and type(export_data_slug) != str:
@@ -109,34 +111,7 @@ def growbikenet(
     nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public', retain_all=False)
 
     if existing_network_spacing:
-        cf = ['["cycleway"~"track"]',
-              '["highway"~"cycleway"]',
-              '["highway"~"path"]["bicycle"~"designated"]',
-              '["cycleway:right"~"track"]',
-              '["cycleway:left"~"track"]',
-              '["cyclestreet"]',
-              '["highway"~"living_street"]'
-             ]
-        for custom_tag in ["cycleway", "bicycle", "cycleway:right", "cycleway:left", "cyclestreet"]:
-            if custom_tag not in ox.settings.useful_tags_way:
-                ox.settings.useful_tags_way.extend(custom_tag)
-        # Fetch protected bike network data from osmnx
-        # Due to retain_all=True, this fetches all the connected components
-        nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True)
-        g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
-
-        # Now we could have some leftover bike infra that is disconnected from the street network and thus not routable.
-        # We delete those parts next:
-        # Take largest connected component lcc of the merged network
-        lcc = max(nx.connected_components(g_undir), key=len)
-        g_undir = g_undir.subgraph(lcc).copy() 
-        # Restrict nodes and edges of the existing bike net to this lcc
-        valid_node_osmids = g_undir.nodes()
-        nodes_exnw = nodes_exnw[nodes_exnw['osmid'].isin(valid_node_osmids)]
-        # edges_exnw has a MultiIndex ('u','v'), so we must use get_level_values, see https://stackoverflow.com/a/18835121
-        edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('u').isin(valid_node_osmids)]
-        edges_exnw = edges_exnw.iloc[edges_exnw.index.get_level_values('v').isin(valid_node_osmids)]
-        _, edges = nx_to_nodes_edges(g_undir, proj_crs)
+        nodes, edges, g_undir, nodes_exnw, edges_exnw = update_with_existing_bike_network(city_name, proj_crs, g_undir)
 
     ### Create seed points
     print("Creating " + seed_point_type + " seed points..")
@@ -163,8 +138,6 @@ def growbikenet(
     
     if existing_network_spacing:
         # If the existing bicycle network is used, create extra seed points on it. They are by construction already snapped.
-        # print(nodes_exnw)
-        # sys.exit()
         seed_points_exnw = get_existing_network_seed_points(nodes_exnw, existing_network_spacing)
         seed_points_exnw.to_crs(edges.crs, inplace=True)
 
@@ -186,17 +159,13 @@ def growbikenet(
         seed_points_snapped.rename(columns={"osmid_1": "osmid"}, inplace=True)
         seed_points_snapped.set_index("osmid", drop=False, inplace=True)
         
-        # seed_points_snapped.drop(["osmid"], axis=1, inplace=True)
-        # print(seed_points_snapped)
-        # seed_points_snapped.to_file("test.gpkg", driver="GPKG")
-        
     # Abort if only 0 or 1 seed points
     if len(seed_points_snapped) < 2:
         raise RuntimeError("Found less than 2 seed points, but more are needed.")
 
     ### Triangulate
     # Triangulation is calculated for the abstract network, but metrics (betweenness, closeness) are calculated for the routed network accounting for lengths.
-    print("Greedy triangulation..")
+    print("Triangulation..")
 
     # Create df with delaunay edges
     df = create_delaunay_edges(seed_points_snapped)
@@ -250,6 +219,7 @@ def growbikenet(
 
     # Generate export data filename
     if export_data or export_plots or export_video:
+        os.makedirs("./results/", exist_ok=True)
         if export_data_slug is None:
             city_string = city_name
         else:
@@ -262,41 +232,43 @@ def growbikenet(
     if export_data:
         ### save data
         print("Saving data..")
-        a_edges.to_file(export_data_filename, driver="GPKG")
+        a_edges.to_file("./results/"+export_data_filename, driver="GPKG")
 
     if export_plots or export_video:
         ### Visualize
         print("Creating visualizations..")
 
-        # Create directories
-        os.makedirs("./results/", exist_ok=True)
-        os.makedirs("./results/plots/", exist_ok=True)
-        os.makedirs("./results/plots/video/", exist_ok=True)
-
         # Read in file to plot
         routed_edges_gdf = gpd.read_file(export_data_filename)
 
         # Viz/plot settings (move to config file later)
-
         # Define color palette (from Michael's project: https://github.com/mszell/bikenwgrowth/blob/main/parameters/parameters.py)
         streetcolor = "#999999"
         edgecolor = "#0EB6D2"
         seedcolor = "#ff7338"
-
         # Define linewidths
-
         lws = {"street": 0.75, "bike": 2}
 
-        create_plots(
-            routed_edges_gdf,
-            seed_points_snapped,
-            streetcolor,
-            edgecolor,
-            seedcolor,
-            lws,
-        )
+        if ranking == "all": 
+            ranking_list = ["betweenness_centrality", "closeness_centrality", "random"]
+        else:
+            ranking_list = [ranking]
+        for ranking_this in ranking_list:
+            os.makedirs("./results/plots/ordering_"+ranking_this+"/", exist_ok=True)
+            create_plots(
+                routed_edges_gdf,
+                seed_points_snapped,
+                streetcolor,
+                edgecolor,
+                seedcolor,
+                lws,
+                ranking_this,
+            )
+
         if export_video:
-            print("Generating video..")
-            make_video(img_folder_name="./results/plots/", fps=1)
+            for ranking_this in ranking_list:
+                print("Generating video..")
+                os.makedirs("./results/plots/ordering_"+ranking_this+"/video/", exist_ok=True)
+                make_video(img_folder_name="./results/plots/ordering_"+ranking_this+"/", fps=5)
 
     return a_edges

From 66fbee6bdf2dae951e88d565b846a9ee824d315c Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 11:20:51 +0200
Subject: [PATCH 6/8] Fix create_plot and results folder structure

---
 .gitignore                    |  2 +-
 growbikenet/visualizations.py | 34 +++++++++++++++++++---------------
 2 files changed, 20 insertions(+), 16 deletions(-)

diff --git a/.gitignore b/.gitignore
index da0b8df..edd770c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,6 @@
 # Folders
 data/
-plots/
+results/
 cache/
 examples/*.gpkg
 /*.gpkg
diff --git a/growbikenet/visualizations.py b/growbikenet/visualizations.py
index 961ce1e..cfad76e 100644
--- a/growbikenet/visualizations.py
+++ b/growbikenet/visualizations.py
@@ -58,28 +58,32 @@ def make_video(
 
 
 def create_plots(
-    routed_edges_gdf, seed_points_snapped, streetcolor, edgecolor, seedcolor, lws
+    routed_edges_gdf, seed_points_snapped, streetcolor, edgecolor, seedcolor, lws, ranking
 ):
-    for ordering in sorted(routed_edges_gdf["ordering"].unique()):
-        fig, ax = plt.subplots(1, 1, figsize=(10, 10))
+    if ranking == "all": 
+        ranking_list = ["betweenness_centrality", "closeness_centrality", "random"]
+    else:
+        ranking_list = [ranking]
 
-        # first, plot street network as "base line"
-        routed_edges_gdf.plot(ax=ax, color=streetcolor, lw=lws["street"], zorder=0)
+    for ranking_this in ranking_list:
+        for ordering in sorted(routed_edges_gdf["ordering_"+ranking_this].unique()):
+            fig, ax = plt.subplots(1, 1, figsize=(10, 10))
 
-        # plot all edges up to current rank
+            # first, plot street network as "base line"
+            routed_edges_gdf.plot(ax=ax, color=streetcolor, lw=lws["street"], zorder=0)
 
-        routed_edges_gdf[routed_edges_gdf["ordering"] <= ordering].plot(
-            ax=ax, color=edgecolor, lw=lws["bike"], zorder=1
-        )
+            # plot all edges up to current rank
 
-        seed_points_snapped.plot(ax=ax, color=seedcolor, zorder=2)
+            routed_edges_gdf[routed_edges_gdf["ordering_"+ranking_this] <= ordering].plot(
+                ax=ax, color=edgecolor, lw=lws["bike"], zorder=1
+            )
 
-        ax.set_axis_off()
+            seed_points_snapped.plot(ax=ax, color=seedcolor, zorder=2)
 
-        plot_id = "{:03d}".format(ordering)  # format plot ID with leading zeros
+            ax.set_axis_off()
 
-        fig.savefig(f"./results/plots/{plot_id}.png", dpi=300)
+            plot_id = "{:03d}".format(ordering)  # format plot ID with leading zeros
 
-        plt.close()
+            fig.savefig(f"./results/plots/ordering_{ranking_this}/{plot_id}.png", dpi=150, bbox_inches='tight')
 
-    return None
+            plt.close()

From fd52bb4e8c448f8dcdf8428606ec4f7376707b3c Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 11:32:12 +0200
Subject: [PATCH 7/8] refactor update_seed_points_with_existing_bike_network

---
 examples/mwe.py            |  2 +-
 growbikenet/functions.py   | 45 ++++++++++++++++++++++++++++++++++++++
 growbikenet/growbikenet.py | 28 ++++--------------------
 3 files changed, 50 insertions(+), 25 deletions(-)

diff --git a/examples/mwe.py b/examples/mwe.py
index 759acc0..4b0e073 100644
--- a/examples/mwe.py
+++ b/examples/mwe.py
@@ -6,7 +6,7 @@
     city_name="Turin",
     proj_crs="3857",
     ranking="betweenness_centrality",
-    existing_network_spacing=200,
+    existing_network_spacing=600,
     export_data=True,
     export_plots=False,
     export_video=False,
diff --git a/growbikenet/functions.py b/growbikenet/functions.py
index 91019b3..430f4ba 100644
--- a/growbikenet/functions.py
+++ b/growbikenet/functions.py
@@ -223,6 +223,51 @@ def update_with_existing_bike_network(city_name, proj_crs, g_undir):
 
     return nodes, edges, g_undir, nodes_exnw, edges_exnw
 
+def update_seed_points_with_existing_bike_network(seed_points_snapped, nodes_exnw, existing_network_spacing, proj_crs):
+    """Update seed points with existing bike network
+
+    Updates given snapped seed points by incorporating seed points from an existing bike network.
+    
+    Parameters
+    ----------
+    seed_points_snapped : geopandas.geodataframe.GeoDataFrame
+        Snapped seed points on the street network, constructed with seed_point_grid_spacing
+    nodes_exnw : geopandas.geodataframe.GeoDataFrame
+        Nodes of the existing bike network
+    existing_network_spacing : int
+        Positive integer denoting spacing between seed points, in meters, only on the existing bicycle network.
+    proj_crs : str
+        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+
+    Returns
+    -------
+    seed_points_snapped : geopandas.geodataframe.GeoDataFrame
+        Snapped seed points incorporating both street grid and existing bike network
+    """
+
+    # If the existing bicycle network is used, create extra seed points on it. They are by construction already snapped.
+    seed_points_exnw = get_existing_network_seed_points(nodes_exnw, existing_network_spacing)
+    seed_points_exnw.to_crs(proj_crs, inplace=True)
+
+    # Afterwards, drop all previously determined seed points (grid or rail) that are now too close to these extra points.
+    buffer_seed_points_exnw = gpd.GeoDataFrame(seed_points_exnw.buffer(existing_network_spacing))
+    buffer_seed_points_exnw = buffer_seed_points_exnw.rename(columns={0:'geometry'}).set_geometry('geometry') # https://gis.stackexchange.com/questions/266098/how-to-convert-a-geoseries-to-a-geodataframe-with-geopandas
+    buffer_seed_points_exnw.to_crs(proj_crs, inplace=True)
+
+    # Delete the seed points that are too close to seed_points_exnw via its buffer
+    seed_points_snapped = seed_points_snapped.overlay(buffer_seed_points_exnw, how='difference')
+
+    # Merge original snapped points with new existing network points (=already snapped)
+    seed_points_snapped = seed_points_snapped.overlay(seed_points_exnw, how='union')
+
+    # Bring back to original form (geometry and osmid columns, osmid index)
+    # This is a bit of a mess but it works. Simplify it in the future.
+    seed_points_snapped.loc[seed_points_snapped['osmid_1'].isnull(), 'osmid_1'] = seed_points_snapped['osmid_2'] # _1 comes from one side, _2 from the other. One has NaNs, the other too. https://stackoverflow.com/a/60132614
+    seed_points_snapped.drop(["y","x","street_count", "highway", "railway", "osmid_2"], axis=1, inplace=True)
+    seed_points_snapped.rename(columns={"osmid_1": "osmid"}, inplace=True)
+    seed_points_snapped.set_index("osmid", drop=False, inplace=True)
+    return seed_points_snapped
+
 def get_grid_seed_points(edges, seed_point_spacing, principal_bearing):
     """Get grid seed points for street network, rotated by principal bearing
 
diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index e5f623b..f0d3df5 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -130,35 +130,15 @@ def growbikenet(
             city_name, {"railway": ["station", "halt"]}
         )
         seed_points = seed_points[seed_points["geometry"].type == "Point"]
-        seed_points.to_crs(edges.crs, inplace=True)
+        seed_points.to_crs(proj_crs, inplace=True)
 
     # Snap seed points to OSM nodes
     seed_points_snapped = snap_seed_points(seed_points, nodes)
     seed_points_snapped = filter_seed_points(seed_points_snapped, seed_point_delta)
     
     if existing_network_spacing:
-        # If the existing bicycle network is used, create extra seed points on it. They are by construction already snapped.
-        seed_points_exnw = get_existing_network_seed_points(nodes_exnw, existing_network_spacing)
-        seed_points_exnw.to_crs(edges.crs, inplace=True)
-
-        # Afterwards, drop all previously determined seed points (grid or rail) that are now too close to these extra points.
-        buffer_seed_points_exnw = gpd.GeoDataFrame(seed_points_exnw.buffer(existing_network_spacing))
-        buffer_seed_points_exnw = buffer_seed_points_exnw.rename(columns={0:'geometry'}).set_geometry('geometry') # https://gis.stackexchange.com/questions/266098/how-to-convert-a-geoseries-to-a-geodataframe-with-geopandas
-        buffer_seed_points_exnw.to_crs(edges.crs, inplace=True)
-
-        # Delete the seed points that are too close to seed_points_exnw via its buffer
-        seed_points_snapped = seed_points_snapped.overlay(buffer_seed_points_exnw, how='difference')
-
-        # Merge original snapped points with new existing network points (=already snapped)
-        seed_points_snapped = seed_points_snapped.overlay(seed_points_exnw, how='union')
-
-        # Bring back to original form (geometry and osmid columns, osmid index)
-        # This is a bit of a mess but it works. Simplify it in the future.
-        seed_points_snapped.loc[seed_points_snapped['osmid_1'].isnull(), 'osmid_1'] = seed_points_snapped['osmid_2'] # _1 comes from one side, _2 from the other. One has NaNs, the other too. https://stackoverflow.com/a/60132614
-        seed_points_snapped.drop(["y","x","street_count", "highway", "railway", "osmid_2"], axis=1, inplace=True)
-        seed_points_snapped.rename(columns={"osmid_1": "osmid"}, inplace=True)
-        seed_points_snapped.set_index("osmid", drop=False, inplace=True)
-        
+        seed_points_snapped = update_seed_points_with_existing_bike_network(seed_points_snapped, nodes_exnw, existing_network_spacing, proj_crs)
+          
     # Abort if only 0 or 1 seed points
     if len(seed_points_snapped) < 2:
         raise RuntimeError("Found less than 2 seed points, but more are needed.")
@@ -215,7 +195,7 @@ def growbikenet(
     # Rank edges by specified method
     a_edges = rank_df(a_edges, ranking)
 
-    a_edges = gpd.GeoDataFrame(a_edges, crs=edges.crs, geometry="geometry")
+    a_edges = gpd.GeoDataFrame(a_edges, crs=proj_crs, geometry="geometry")
 
     # Generate export data filename
     if export_data or export_plots or export_video:

From 9692fb41367379cbc117e3a22e9f8f33630c3f52 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Wed, 29 Apr 2026 11:44:05 +0200
Subject: [PATCH 8/8] Fix imports

---
 growbikenet/growbikenet.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index fc4ebdb..1b3b772 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -16,6 +16,9 @@
     node_to_edge_attributes,
     df_from_graph,
     rank_df,
+    prepare_network,
+    update_with_existing_bike_network,
+    update_seed_points_with_existing_bike_network,
 )
 from growbikenet.visualizations import make_video, create_plots