From 33b21562c65b3202297c2c45c3c666d004ff71bb Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 08:37:53 +0200
Subject: [PATCH 1/7] Remove release-please

---
 .github/workflows/release-please.yml | 19 -------------------
 .release-please-manifest.json        |  3 ---
 bikenetlib/_version.py               |  2 +-
 release-please-config.json           | 12 ------------
 4 files changed, 1 insertion(+), 35 deletions(-)
 delete mode 100644 .github/workflows/release-please.yml
 delete mode 100644 .release-please-manifest.json
 delete mode 100644 release-please-config.json

diff --git a/.github/workflows/release-please.yml b/.github/workflows/release-please.yml
deleted file mode 100644
index 1a2c43a..0000000
--- a/.github/workflows/release-please.yml
+++ /dev/null
@@ -1,19 +0,0 @@
-name: release-please
-
-on:
-  push:
-    branches:
-      - main
-
-permissions:
-  contents: write
-  pull-requests: write
-
-jobs:
-  release-please:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: googleapis/release-please-action@v4
-        with:
-          token: ${{ secrets.GITHUB_TOKEN }}
-          config-file: release-please-config.json
diff --git a/.release-please-manifest.json b/.release-please-manifest.json
deleted file mode 100644
index f1c1e58..0000000
--- a/.release-please-manifest.json
+++ /dev/null
@@ -1,3 +0,0 @@
-{
-  ".": "0.5.0"
-}
diff --git a/bikenetlib/_version.py b/bikenetlib/_version.py
index f75962d..553f839 100644
--- a/bikenetlib/_version.py
+++ b/bikenetlib/_version.py
@@ -1,3 +1,3 @@
 """bikenetlib package version."""
 
-__version__ = "0.5.0"  # x-release-please-version
+__version__ = "0.5.0"
diff --git a/release-please-config.json b/release-please-config.json
deleted file mode 100644
index f4d36d0..0000000
--- a/release-please-config.json
+++ /dev/null
@@ -1,12 +0,0 @@
-{
-  "packages": {
-    ".": {
-      "release-type": "python",
-      "changelog-path": "docs/changelog.md",
-      "extra-files": [
-        "bikenetlib/_version.py",
-        "docs/conf.py"
-      ]
-    }
-  }
-}

From 5f78f913213aea11081eb6cbdc084e9490350c4e Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 08:56:37 +0200
Subject: [PATCH 2/7] Update project and workflows

---
 .github/workflows/publish.yml |  2 +-
 .github/workflows/test.yml    | 15 +++++++-------
 .pre-commit-config.yaml       |  2 +-
 bikenetlib/osm.py             |  0
 environment-dev.yml           | 36 ++++++++++++++------------------
 environment.yml               | 20 +++---------------
 pyproject.toml                | 39 +++++++++++++++++++----------------
 tests/test_osm.py             |  1 +
 8 files changed, 50 insertions(+), 65 deletions(-)
 create mode 100644 bikenetlib/osm.py
 create mode 100644 tests/test_osm.py

diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 62cd726..e496489 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -22,7 +22,7 @@ jobs:
       url: https://pypi.org/project/bikenetlib/
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v7
 
       - name: Build package
         run: |
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 2656162..6aaf719 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -21,7 +21,7 @@ jobs:
   caching:
     strategy:
       matrix:
-        python-version: ["3.12"]
+        python-version: ["3.12", "3.15"]
         os: [ubuntu-latest]
       fail-fast: false
     env:
@@ -34,13 +34,13 @@ jobs:
     name: Cache for ${{ matrix.python-version }} on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v7
 
       - name: Get week number
         run: echo "WEEK=$(date +'%V')" >> $GITHUB_ENV
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: mamba-org/setup-micromamba@v2
+        uses: mamba-org/setup-micromamba@v3
         with:
           create-args: python=${{ matrix.python-version }}
           environment-file: environment-dev.yml
@@ -51,12 +51,13 @@ jobs:
     needs: caching
     strategy:
       matrix:
-        python-version: ["3.12"]
+        python-version: ["3.12", "3.15"]
         os: [ubuntu-latest]
         submodule:
-          #- { name: "io", pytest_args: "tests/test_io.py" } # To do
+          # - { name: "io", pytest_args: "tests/test_io.py" } # To do
+          # - { name: "osm", pytest_args: "tests/test_osm.py" } # To do
           - { name: "utils", pytest_args: "tests/test_utils.py" }
-          #- { name: "viz", pytest_args: "tests/test_viz.py" } # To do
+          # - { name: "viz", pytest_args: "tests/test_viz.py" } # To do
       fail-fast: false
 
     name: ${{ matrix.submodule.name }} (${{ matrix.python-version }} on ${{ matrix.os }})
@@ -73,7 +74,7 @@ jobs:
         run: echo "WEEK=$(date +'%V')" >> $GITHUB_ENV
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: mamba-org/setup-micromamba@v2
+        uses: mamba-org/setup-micromamba@v3
         with:
           create-args: python=${{ matrix.python-version }}
           environment-file: environment-dev.yml
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a041caa..8427de9 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -23,7 +23,7 @@ repos:
         types_or: [markdown, yaml]
   - repo: https://github.com/astral-sh/ruff-pre-commit
     # Ruff version.
-    rev: v0.15.12
+    rev: v0.16.7
     hooks:
       # Run the linter.
       - id: ruff
diff --git a/bikenetlib/osm.py b/bikenetlib/osm.py
new file mode 100644
index 0000000..e69de29
diff --git a/environment-dev.yml b/environment-dev.yml
index 80c452a..4e7617b 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -2,32 +2,26 @@ name: bnlenvdev
 channels:
   - conda-forge
 dependencies:
-  - python=3.12
-  - geopandas>=0.14
-  - jupyter
-  - jupyterlab
-  - ipywidgets
-  - ipykernel
-  - matplotlib
-  - osmnx>=1.9.4
-  - tqdm
+  - python>=3.12
+  - matplotlib=3.11.0 # The exact version is important for visualization tests
+  - osmnx>=2
+  - scipy
   - geojson
-  - pip
-  - momepy
+
+  # additional to the non-dev env
+  - ipykernel
+  - ipywidgets
   - pytest
-  - contextily
+  - pytest-mpl
+  - pytest-cov
   - sphinx
   - numpydoc
   - myst-nb
-  - sphinx-book-theme
-  - nbsphinx
-  - myst-parser
-  - ipython
+  - pre-commit
+  - jupyterlab
+  - pixi-kernel
   - pip:
       # not available via conda-forge:
-      - opencv-python
-      - pre-commit
-      - sphinxcontrib-bibtex
       - sphinx-copybutton
-      - sphinx-gallery
-      - python-slugify
+      - sphinxext-opengraph
+      - furo
diff --git a/environment.yml b/environment.yml
index 1dc39db..b14423a 100644
--- a/environment.yml
+++ b/environment.yml
@@ -2,22 +2,8 @@ name: bnlenv
 channels:
   - conda-forge
 dependencies:
-  - python=3.12
-  - geopandas>=0.14
-  - jupyter
-  - jupyterlab
-  - ipywidgets
-  - ipykernel
+  - python>=3.12
   - matplotlib
-  - osmnx>=1.9.4
-  - tqdm
+  - osmnx>=2
+  - scipy
   - geojson
-  - pip
-  - momepy
-  - pytest
-  - contextily
-  - pip:
-      # not available via conda-forge:
-      - opencv-python
-      - pre-commit
-      - python-slugify
diff --git a/pyproject.toml b/pyproject.toml
index ce8d591..0a8fc1b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,14 +4,15 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "bikenetlib"
-# dynamic = ["version"]
 version = "0.5.0"
 authors = [
-  { name = "Manuel Knepper" },
-  { name = "Anastassia Vybornova" },
-  { name = "Michael Szell" },
+    { name = "Michael Szell", email = "michael@szell.net" },
+    { name = "Manuel Knepper" },
+    { name = "Anastassia Vybornova" },
+]
+maintainers = [
+    { name = "Michael Szell", email = "michael@szell.net" },
 ]
-maintainers = [{ name = "Michael Szell", email = "michael@szell.net" }]
 license = "AGPL-3.0-or-later"
 readme = "README.md"
 description = "BikeNetKit Python core utility package"
@@ -20,40 +21,42 @@ classifiers = [
     "Development Status :: 3 - Alpha",
     "Intended Audience :: Developers",
     "Intended Audience :: Science/Research",
+    "Intended Audience :: Other Audience",
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
     "Programming Language :: Python :: 3 :: Only",
-    "Programming Language :: Python :: 3.10",
-    "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
+    "Programming Language :: Python :: 3.15",
     "Natural Language :: English",
+    "Topic :: Scientific/Engineering :: GIS",
 ]
-requires-python = ">=3.10"
+requires-python = ">=3.12"
 dependencies = [
-    "geopandas>=0.14",
+    "geopandas>=1",
     "matplotlib",
-    "osmnx>=1.9.4",
+    "osmnx>=2",
+    "scipy",
     "geojson",
-    "tqdm",
-    "pip",
-    "momepy",
-    "opencv-python",
-    "python-slugify",
 ]
 
 [project.urls]
+Homepage = "https://bikenetkit.org/"
+Documentation = "https://docs.bikenetkit.org/BikeNetLib"
 Repository = "https://github.com/BikeNetKit/BikeNetLib"
 Issues = "https://github.com/BikeNetKit/BikeNetLib/issues"
+Changelog = "https://docs.bikenetkit.org/BikeNetKit/changelog"
 
 [project.optional-dependencies]
 test = ["pytest"]
-doc = ["sphinx", "momepy"]
+doc = ["sphinx"]
 
 [tool.ruff.lint.extend-per-file-ignores]
 "__init__.py" = ["F401"]  # "Imported but unused: happens with packages
 
-# [tool.setuptools.dynamic]
-# version = {attr = "bikenetlib.__version__"}
+[tool.ruff.lint.pydocstyle]
+convention = "numpy"
 
 [tool.setuptools.packages.find]
 exclude = [
diff --git a/tests/test_osm.py b/tests/test_osm.py
new file mode 100644
index 0000000..17efb5f
--- /dev/null
+++ b/tests/test_osm.py
@@ -0,0 +1 @@
+# To do
\ No newline at end of file

From 9d10add394ed9266668c2032b3354b9b4cfcc38a Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 08:58:08 +0200
Subject: [PATCH 3/7] Add Python 3.11 support

---
 .github/workflows/test.yml | 4 ++--
 pyproject.toml             | 3 ++-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 6aaf719..db810ea 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -21,7 +21,7 @@ jobs:
   caching:
     strategy:
       matrix:
-        python-version: ["3.12", "3.15"]
+        python-version: ["3.11", "3.15"]
         os: [ubuntu-latest]
       fail-fast: false
     env:
@@ -51,7 +51,7 @@ jobs:
     needs: caching
     strategy:
       matrix:
-        python-version: ["3.12", "3.15"]
+        python-version: ["3.11", "3.15"]
         os: [ubuntu-latest]
         submodule:
           # - { name: "io", pytest_args: "tests/test_io.py" } # To do
diff --git a/pyproject.toml b/pyproject.toml
index 0a8fc1b..2374003 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -25,6 +25,7 @@ classifiers = [
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
     "Programming Language :: Python :: 3 :: Only",
+    "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
     "Programming Language :: Python :: 3.14",
@@ -32,7 +33,7 @@ classifiers = [
     "Natural Language :: English",
     "Topic :: Scientific/Engineering :: GIS",
 ]
-requires-python = ">=3.12"
+requires-python = ">=3.11"
 dependencies = [
     "geopandas>=1",
     "matplotlib",

From 6b74f023bc4017b2ac14ee90151a3bc7d7c3af56 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 09:06:03 +0200
Subject: [PATCH 4/7] Update actions versions

---
 .github/workflows/test.yml | 2 +-
 .pre-commit-config.yaml    | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index db810ea..d0811ff 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -68,7 +68,7 @@ jobs:
         shell: bash -el {0}
 
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v7
 
       - name: Get week number
         run: echo "WEEK=$(date +'%V')" >> $GITHUB_ENV
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8427de9..c05ea4a 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -16,14 +16,14 @@ repos:
       - id: mixed-line-ending
       - id: check-illegal-windows-names
 
-  - repo: https://github.com/pre-commit/mirrors-prettier
-    rev: "v3.1.0"
+  - repo: https://github.com/pre-commit/pre-commit
+    rev: "v4.6.2"
     hooks:
       - id: prettier
         types_or: [markdown, yaml]
   - repo: https://github.com/astral-sh/ruff-pre-commit
     # Ruff version.
-    rev: v0.16.7
+    rev: v0.17.0
     hooks:
       # Run the linter.
       - id: ruff

From 79a6c7ba36ce5da603226664570915c7cc06492d Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 10:11:00 +0200
Subject: [PATCH 5/7] Add osm module and polished assign_pbi_to_edges()

---
 bikenetlib/osm.py   | 78 +++++++++++++++++++++++++++++++++++++++++++++
 bikenetlib/utils.py | 29 +++++++++++++++++
 2 files changed, 107 insertions(+)

diff --git a/bikenetlib/osm.py b/bikenetlib/osm.py
index e69de29..71c89f9 100644
--- a/bikenetlib/osm.py
+++ b/bikenetlib/osm.py
@@ -0,0 +1,78 @@
+"""Global constants for `bikenetlib` related to OpenStreetMap(OSM).
+
+GROWABLE_NETWORK_CUSTOM_FILTER : list[str] or None
+    Custom filter for all infrastructure elements that are considered as 
+    growable by `growbikenet`. By default, `growbikenet` uses a custom filter 
+    to retrieve the combined drive and pbi (protected bicycle infrastructure) 
+    network. To only consider the drive network, set 
+    `GROWABLE_NETWORK_CUSTOM_FILTER` to None and `GROWABLE_NETWORK_TYPE` to 
+    'drive'. However, doing so can lead to issues: 
+    https://github.com/BikeNetKit/GrowBikeNet/issues/255. 
+GROWABLE_NETWORK_TYPE : {'drive', 'all', 'all_public', 'bike', 'drive_service', 'walk'}, default 'drive' 
+        What type of street network to retrieve for the growable network if `GROWABLE_NETWORK_CUSTOM_FILTER` is None.
+PBI_CUSTOM_FILTER : list[str]
+    Custom filter for protected bicycle infrastructure (pbi).
+PBI_DICT : dict
+    Dictionary for protected bicycle infrastructure (pbi), for single criteria.
+PBI_DICT_SUB : dict
+    Dictionary for protected bicycle infrastructure (pbi), with one 
+    sub-criterion. For example, highway~path AND bicycle~designated. Always 
+    adds a second criterion access!~private.
+ROUTING_PENALTY : dict, default {0: 1.5, 1: 1}
+    Factor to multiply length of non-pbi/pbi for routing, to avoid routing
+    through parallel streets when slightly longer pbi is available. By 
+    default, non-pbi counts as 50% longer than pbi.
+"""
+
+# Populate ox.settings.useful_tags_way to make application of custom filter possible
+import osmnx as ox
+for custom_tag in ["highway", "cycleway", "bicycle", "cycleway:right", "cycleway:left", "cycleway:both", "cyclestreet", "access", "area", "service", "motor_vehicle", "motorcar"]: # This list should contain all tags used in any custom filters
+    if custom_tag not in ox.settings.useful_tags_way:
+        ox.settings.useful_tags_way.extend(custom_tag)
+
+
+GROWABLE_NETWORK_CUSTOM_FILTER = [ # adapted from https://github.com/gboeing/osmnx/blob/2fc39cb2792ff869881b99f724a2d97f0e958667/osmnx/_overpass.py#L77
+    # Car infra - Instead of excluding what we don't want, we only include what we want. Note: we are OK with alleys and driveways
+    '["highway"~"motorway|trunk|primary|secondary|tertiary|residential|motorway_link|trunk_link|primary_link|secondary_link|tertiary_link|service|unclassified"]["highway"!~"abandoned|construction|no|planned|platform|proposed|raceway|razed|rest_area|services"]["service"!~"emergency_access|parking|parking_aisle|private"]["area"!~"yes"]["access"!~"private"]["motor_vehicle"!~"no"]["motorcar"!~"no"]',
+    # Bike infra, copied from PBI_CUSTOM_FILTER
+    '["highway"~"cycleway|living_street"]',
+    '["highway"~"path|pedestrian|footway"]["bicycle"~"designated|yes|permissive"]["access"!~"private"]',
+    '["cyclestreet"~"yes"]',
+    '["bicycle_road"~"yes"]',
+    '["cycleway"~"track"]',
+    '["cycleway:right"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+    '["cycleway:left"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+    '["cycleway:both"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+]
+
+GROWABLE_NETWORK_TYPE = 'drive'
+
+PBI_CUSTOM_FILTER = [
+    '["highway"~"cycleway|living_street"]',
+    '["highway"~"path|pedestrian|footway"]["bicycle"~"designated|yes|permissive"]["access"!~"private"]',
+    '["cyclestreet"~"yes"]',
+    '["bicycle_road"~"yes"]',
+    '["cycleway"~"track"]',
+    '["cycleway:right"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+    '["cycleway:left"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+    '["cycleway:both"~"track|opposite_track"]', # opposite_track is deprecated, but could still exist
+]
+
+PBI_DICT = {
+    'highway': ['cycleway','living_street'],
+    'cyclestreet': ['yes'],
+    'bicycle_road': ['yes'],
+    'cycleway': ['track'],
+    'cycleway:right': ['track','opposite_track'],
+    'cycleway:left': ['track','opposite_track'],
+    'cycleway:both': ['track','opposite_track'],
+}
+PBI_DICT_SUB = {
+    'highway': [
+        ['path','pedestrian','footway'],
+        {'bicycle': ['designated','yes','permissive']},
+    ]
+}
+
+ROUTING_PENALTY = {0: 1.5, 1: 1}
+
diff --git a/bikenetlib/utils.py b/bikenetlib/utils.py
index 90cddf5..f62ed26 100644
--- a/bikenetlib/utils.py
+++ b/bikenetlib/utils.py
@@ -7,6 +7,35 @@
 from shapely.prepared import prep
 from shapely.geometry import Point, MultiLineString
 from tqdm import tqdm
+from . import osm
+
+
+def assign_pbi_to_edges(g):
+    """Assign a boolean pbi attribute to all edges in a graph depending if they 
+    are considered protected bike infrastructure (pbi).
+
+    Parameters
+    ----------
+    g : networkx.MultiDiGraph
+        Simplified graph representing the street network.
+
+    Returns
+    -------
+    g : networkx.MultiDiGraph
+        Simplified graph representing the street network, with added binary edge attribute "pbi".
+    """
+    for edge in g.edges(keys=True):
+        assigned_pbi = 0
+        for tag in osm.PBI_DICT:
+            if g.edges[edge].get(tag) in osm.PBI_DICT[tag]:
+                assigned_pbi = 1
+        if not assigned_pbi: # Check if the edge fulfils a subcriterion
+            for tag in osm.PBI_DICT_SUB:
+                subtag = osm.PBI_DICT_SUB[tag][1].keys()[0]
+                if g.edges[edge].get(tag) in osm.PBI_DICT_SUB[tag][0] and g.edges[edge].get(subtag) in osm.PBI_DICT_SUB[tag][1][subtag] and g.edges[edge].get("access") != 'private':
+                    assigned_pbi = 1
+        g.edges[edge]["pbi"] = assigned_pbi
+    return g
 
 
 def intersects_properly(geom1, geom2):

From c6f5e14aad7ea34bba0e7c673fb353cdbed814cb Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 12:00:53 +0200
Subject: [PATCH 6/7] Cut and polish functions

---
 bikenetlib/io.py    |  58 +---
 bikenetlib/utils.py | 666 +++-----------------------------------------
 bikenetlib/viz.py   |   1 +
 environment-dev.yml |   1 -
 environment.yml     |   1 -
 pyproject.toml      |   1 -
 tests/test_utils.py |  60 ----
 7 files changed, 47 insertions(+), 741 deletions(-)

diff --git a/bikenetlib/io.py b/bikenetlib/io.py
index b03003a..992d1fb 100644
--- a/bikenetlib/io.py
+++ b/bikenetlib/io.py
@@ -1,57 +1 @@
-import geopandas as gpd
-import osmnx as ox
-
-def prepare_network(city_name, proj_crs, network_type='all_public', custom_filter=None, retain_all=True, city_boundary_file=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. Overruled (for data fetching) if city_boundary_file is set.
-    proj_crs : str
-        Coordinate reference system that is used to project osm data.
-    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
-    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.
-    city_boundary_file : (str | None), default None
-        If not set to None, the study area will be selected from the (Multi)Polygon provided in the city_boundary_file shape file. For example, "copenhagen.shp".
-
-    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
-    if city_boundary_file is None:
-        g = ox.graph_from_place(
-        city_name, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
-        )
-    else:
-        shp = gpd.read_file(city_boundary_file)
-        city_boundary_polygon = shp.iloc[0].geometry
-        g = ox.graph_from_polygon(
-        city_boundary_polygon, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
-        )
-
-    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 download_pois(): # To develop
-    pass
-
-def save_to_file(): # To develop
-    pass
\ No newline at end of file
+"""Import/Output functions of `bikenetlib`."""
\ No newline at end of file
diff --git a/bikenetlib/utils.py b/bikenetlib/utils.py
index f62ed26..f5b9495 100644
--- a/bikenetlib/utils.py
+++ b/bikenetlib/utils.py
@@ -1,12 +1,10 @@
+"""Utility functions of `bikenetlib`."""
+
 import numpy as np
 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, MultiLineString
-from tqdm import tqdm
 from . import osm
 
 
@@ -37,235 +35,37 @@ def assign_pbi_to_edges(g):
         g.edges[edge]["pbi"] = assigned_pbi
     return g
 
-
-def intersects_properly(geom1, geom2):
-    """
-    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
-    ----------
-    geom1 : shapely geometry
-        A shapely geometry, for example shapely.geometry.Point() or shapely.geometry.LineString()
-    geom2 : shapely geometry
-        A shapely geometry, for example shapely.geometry.Point() or shapely.geometry.LineString()
-
-    Returns
-    -------
-    boolean
-        Returns true if the two provided geometries intersect but do not touch
-    """
-    return geom1.intersects(geom2) and not geom1.touches(geom2)
-
-
-def nx_to_nodes_edges(G, proj_crs):
-    """Get nodes and projected edges from networkX graph
     
-    Parameters
-    ----------
-    G : networkx.classes.multigraph.MultiGraph
-        networkX graph, undirected
-    proj_crs : str
-        Coordinate reference system that is used to project osm data.
-        
-    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)
+def nodelist_to_edgetuples(edge_gdf, nodelist):
+    """Map a list of ndoes (from `nx.shortest_paths()`) to a list of edge 
+    tuples that can be used for indexing an edge geodataframe.
 
-    # 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)
-    to the correct set of edge tuples that can be used for INDEXING THE EDGE GDF
+    Useful to apply as a lambda function to turn path nodes into path edges.
 
     Parameters
     ----------
-    edge_gdf: geopandas.geodataframe.GeoDataFrame
-        The street network, in a projected coordinate reference system
-    nodelist: list
-        A list of nodes that make up source and targets of edges
+    edge_gdf : geopandas.geodataframe.GeoDataFrame
+        The edges of a spatial network, in a projected CRS.
+    nodelist : list
+        A list of nodes that make up source and targets of edges.
 
     Returns
     -------
-    edgelist_final: list
-        List of edge tuples that can be used for INDEXING THE EDGE GDF
+    edgelist_final : list
+        List of edge tuples that can be used for indexing the edge 
+        geodataframe.
     """
     edgelist_prelim = zip(nodelist, nodelist[1:])
     edgelist_final = []
+    temp_gdf = edge_gdf.sort_index() # To circumvent PerformanceWarning, see https://stackoverflow.com/questions/54307300/what-causes-indexing-past-lexsort-depth-warning-in-pandas
     for edge_prelim in edgelist_prelim:
-        if edge_prelim in edge_gdf.index:
+        if edge_prelim in temp_gdf.index:
             edgelist_final.append(edge_prelim)
         else:
             edgelist_final.append(tuple([edge_prelim[1], edge_prelim[0]]))
     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)>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
-
-        # 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") 
-        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
-    
-def update_with_existing_bike_network(city_name, proj_crs, g_undir, city_boundary_file=None):
-    """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. Overruled (for data fetching) if city_boundary_file is set.
-    proj_crs : str
-        Coordinate reference system that is used to project osm data.
-    g_undir : networkx.classes.multigraph.MultiGraph
-        Street network networkX graph, undirected
-    city_boundary_file : (str | None), default None
-        If not set to None, the study area will be selected from the (Multi)Polygon provided in the city_boundary_file shape file. For example, "copenhagen.shp".
-
-    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, city_boundary_file=city_boundary_file)
-    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 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.
-
-    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 = seed_points_snapped[['osmid_1','geometry']]
-    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_principal_bearing(G):
     """Determine the most common (principal) bearing, for the best grid orientation.
 
@@ -296,35 +96,35 @@ def get_principal_bearing(G):
         except:  # noqa (To do: make specific and remove noqa)
             pass  # Bearings cannot be calculated in rare edge cases.
     b = pd.Series(city_bearings)
-    bearings = pd.concat([b, b.map(reverse_bearing)]).reset_index(drop="True")
+    bearings = pd.concat([b, b.map(_reverse_bearing)]).reset_index(drop="True")
     bins = np.arange(bearingbins + 1) * 360 / bearingbins
-    count = count_and_merge(bearingbins, bearings)
+    count = _count_and_merge(bearingbins, bearings)
     principal_bearing = bins[np.where(count == max(count))][0]
 
     return principal_bearing
 
 
-def reverse_bearing(x):
+def _reverse_bearing(x):
     """Reverse bearing.
 
     Adapted from: https://github.com/gboeing/osmnx-examples/blob/v0.11/notebooks/17-street-network-orientations.ipynb
 
     Parameters
     ----------
-    x: float
-        The bearing to reverse
+    x : float
+        The bearing to reverse.
 
     Returns
     -------
-    x_rev: float
-        The reversed bearing
+    x_rev : float
+        The reversed bearing.
     """
     x_rev = x + 180 if x < 180 else x - 180
     return x_rev
 
 
-def count_and_merge(n, bearings):
-    """Double, then merge bins to avoid edge effects
+def _count_and_merge(n, bearings):
+    """Double, then merge bins to avoid edge effects.
 
     Make twice as many bins as desired, then merge them in pairs.
     Prevents bin-edge effects around common values like 0° and 90°.
@@ -332,15 +132,15 @@ def count_and_merge(n, bearings):
 
     Parameters
     ----------
-    n: int
-        Number of bins
-    bearings: pandas.Series
-        Series of bearings
+    n : int
+        Number of bins.
+    bearings : pandas.Series
+        Series of bearings.
 
     Returns
     -------
-    bearings_merged: numpy.ndarray, dtype=int
-        The frequencies of the new merged bearings
+    bearings_merged : numpy.ndarray, dtype=int
+        The frequencies of the new merged bearings.
     """
     n *= 2
     bins = np.arange(n + 1) * 360 / n
@@ -352,261 +152,6 @@ def count_and_merge(n, bearings):
     return bearings_merged
 
 
-def snap_seed_points(seed_points, nodes):
-    """
-    snap generated seed_points to actual osm nodes
-    Parameters
-    ----------
-    seed_points: geopandas.geodataframe.GeoDataFrame
-        Seed points that were created within city area, to be snapped to actual osm nodes
-    nodes: geopandas.geodataframe.GeoDataFrame
-        actual osm nodes, downloaded from osmnx
-
-    Returns
-    -------
-    seed_points_snapped: geopandas.geodataframe.GeoDataFrame
-        seed_points with additional information about geometries of osm nodes that seed nodes were snapped to
-
-    """
-    # Ensure same CRS
-    if seed_points.crs != nodes.crs:
-        seed_points = seed_points.to_crs(nodes.crs)
-
-    # Find nearest nodes (returns indices)
-    idx_seed, idx_nodes = nodes.sindex.nearest(seed_points.geometry, return_all=False)
-
-    # Assign osmid safely
-    seed_points = seed_points.copy()
-    seed_points["osmid"] = nodes.iloc[idx_nodes]["osmid"].values
-
-    # Keep original geometry
-    seed_points = seed_points.rename(columns={"geometry": "geometry_generated"})
-
-    # Attach node geometry + attributes
-    nodes_subset = nodes[["osmid", "geometry"]].rename(
-        columns={"geometry": "geometry_osm"}
-    )
-
-    seed_points = seed_points.reset_index(drop=True)
-    nodes_subset = nodes_subset.reset_index(drop=True)
-
-    seed_points_snapped = seed_points.merge(nodes_subset, on="osmid", how="left")
-
-    seed_points_snapped.set_geometry("geometry_osm")
-
-    return seed_points_snapped
-
-
-def filter_seed_points(seed_points_snapped, seed_point_delta):
-    """
-    remove seed_points that are further than delta away from an actual osm node
-    Parameters
-    ----------
-    seed_points_snapped: geopandas.geodataframe.GeoDataFrame
-        seed_points with additional information about geometries of osm nodes that seed nodes were snapped to
-    seed_point_delta: int
-        maximum distance a seed_point may be removed from an actual osm node
-
-    Returns
-    -------
-    seed_points_snapped: geopandas.geodataframe.GeoDataFrame
-        seed_points within delta away from an actual osm node, only columns are osmid and the associated osm geometry
-    """
-    gdf = seed_points_snapped.copy()
-
-    # Compute distance
-    gdf["snap_dist"] = gdf.geometry_generated.distance(gdf.geometry_osm)
-
-    # Filter by threshold
-    gdf = gdf[gdf["snap_dist"] <= seed_point_delta].copy()
-
-    # Drop duplicates: one row per osmid
-    gdf = gdf.sort_values("snap_dist").drop_duplicates("osmid")
-
-    # Keep only node geometry
-    gdf = gdf[["osmid", "geometry_osm"]].rename(columns={"geometry_osm": "geometry"})
-
-    # Set geometry + index
-    gdf = gdf.set_geometry("geometry")
-    gdf = gdf.set_index("osmid", drop=False)
-
-    seed_points_snapped = gdf.copy()
-
-    return seed_points_snapped
-
-
-def create_delaunay_edges(nodes_gdf):
-    """Create df with edges that are part of Delaunay triangulation
-
-    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 and metrics (betweenness, closeness) are calculated for the abstract network for which egde lengths are taken from the routed network.
-
-    Parameters
-    ----------
-    nodes_gdf: geopandas.geodataframe.GeoDataFrame
-        seed points with osmid and corresponding point geometry
-
-    Returns
-    -------
-    df : pandas.DataFrame
-        DataFrame with Edge pairs and singled out source and target nodes
-
-    References
-    ----------
-    .. [1] M. Szell, S. Mimar, T. Perlman, G. Ghoshal, R. Sinatra, "Growing urban bicycle networks", Scientific Reports 12, 6765 (2022)
-    """
-    # Ensure projected CRS
-    if nodes_gdf.crs.is_geographic:
-        raise ValueError("CRS must be projected for triangulation.")
-
-    # Extract coordinates
-    coords = np.array([(geom.x, geom.y) for geom in nodes_gdf.geometry])
-    osmids = nodes_gdf["osmid"].values
-
-    # Compute triangulation
-    tri = Delaunay(coords)
-
-    edges_set = set()
-
-    # Each triangle has 3 edges
-    for simplex in tri.simplices:
-        i, j, k = simplex
-
-        edges_set.add(tuple(sorted((i, j))))
-        edges_set.add(tuple(sorted((j, k))))
-        edges_set.add(tuple(sorted((i, k))))
-
-    pairs = []
-    sources = []
-    targets = []
-
-    for i, j in edges_set:
-        pairs.append((osmids[i], osmids[j]))
-        sources.append(osmids[i])
-        targets.append(osmids[j])
-
-    df = pd.DataFrame(
-        {
-            "pair": pairs,
-            "source": sources,
-            "target": targets,
-        }
-    )
-
-    return df
-
-
-def remove_edge_overlaps(edges_in):
-    """In the grown network, remove edge overlaps stepwise
-
-    Parameters
-    ----------
-    edges_in: geopandas.geodataframe.GeoDataFrame
-        The grown bike network, in a projected coordinate reference system
-
-    Returns
-    -------
-    edges_out: geopandas.geodataframe.GeoDataFrame
-        The grown bike network without edge overlaps, in a projected coordinate reference system
-    """
-    # import time
-    # start = time.time()
-
-    edges_out = edges_in.copy()
-    grown_net = MultiLineString()
-    for row in tqdm(
-        edges_in.itertuples(),
-        desc="{:<26}".format("Removing edge overlaps"),
-        leave=True,
-        unit="edge",
-        total=len(list(edges_in.itertuples())),
-        bar_format='{l_bar}{bar:20}{r_bar}',
-        ):
-        grown_net_new = grown_net | row.geometry # Calculate union
-        if grown_net_new.length > grown_net.length: # Something was added
-            grown_net_diff = row.geometry - grown_net # Calculate difference
-            # print(row.geometry, grown_net, grown_net_diff)
-            edges_out.loc[row.Index, ['geometry']] = grown_net_diff # Add difference
-            grown_net = grown_net_new # Only update if something was added
-        else: # There was nothing added, so we delete the row
-            edges_out.drop(index=row.Index, inplace=True)
-    edges_out.drop_duplicates(inplace=True) # How can duplicates happen??
-    edges_out.reset_index(drop=True, inplace=True)
-
-    # end = time.time()
-    # print(end - start)
-    return edges_out
-
-def df_from_graph(A, method):
-    """Create a dataframe from an input graph
-
-    Parameters
-    ----------
-    A: networkx.graph
-        Graph created from triangulation edge list
-    method: str
-        Method used to rank edges. Must be 'betweenness_centrality' (default), 'closeness_centrality', or 'random'.
-
-    Returns
-    -------
-    df: pandas.DataFrame
-        Dataframe with source and target information for each edge, as well as edge attributes as columns
-    """
-
-    if method != "random":
-        attrs = {
-            edge: {
-                method: data.get(method),
-                "geometry": data.get("geometry"),
-            }
-            for edge, data in A.edges.items()
-        }
-        df = pd.DataFrame.from_dict(attrs, orient="index", columns=[method, "geometry"])
-    else:
-        attrs = {
-            edge: {
-                "geometry": data.get("geometry"),
-            }
-            for edge, data in A.edges.items()
-        }
-        df = pd.DataFrame.from_dict(attrs, orient="index", columns=["geometry"])
-    df["node_tuple"] = df.index
-    df["source"] = [t[0] for t in df.node_tuple]
-    df["target"] = [t[1] for t in df.node_tuple]
-    df.drop(columns=["node_tuple"], inplace=True)
-    return df
-
-
-def rank_df(df, method):
-    """Rank dataframe by specified method
-
-    Parameters
-    ----------
-    df: pandas.DataFrame
-        Dataframe with source and target information for each edge, as well as edge attributes as columns
-    method: str
-        Method used to rank edges. Must be 'betweenness_centrality' (default), 'closeness_centrality', or 'random'.
-
-    Results
-    -------
-    df: pandas.DataFrame
-        Dataframe sorted by specified ranking method.
-    """
-    if method == "random": # ranking is random
-        df["ordering_random"] = np.random.permutation(np.arange(df.shape[0]))
-        df = df.sort_values(by="ordering_random", ascending=False)
-        df.reset_index(drop=True, inplace=True)
-        df["ordering_" + method] = (
-            df.index
-        )  
-    else: # ranking is the order of appearance in the method's ranking
-        df = df.sort_values(by=method, ascending=False)
-        df.reset_index(drop=True, inplace=True)
-        df["ordering_" + method] = (
-            df.index
-        )  
-    return df
-
-
 def node_to_edge_attributes(values_nodes, edges):
     """Map node to edge attributes.
 
@@ -616,7 +161,7 @@ def node_to_edge_attributes(values_nodes, edges):
     ----------
     values_nodes : dict
         Keys: node ids, Values: Node attributes (for example a scalar)
-    edges : networkx.classes.reportviews.EdgeView
+    edges : networkx.classes.reportviews.EdgeView 
         A view of edge attributes of a networkx graph. Could also be a list of tuples of node ids.
 
     Returns
@@ -630,155 +175,34 @@ def node_to_edge_attributes(values_nodes, edges):
     return values_edges
 
 
-def add_path_to_df(df, edges, g):
-    """
+def route_nodepairs(nodepairs, edges, g_undir):
+    """Route pairs of nodes on the graph `g_undir` to a merged geometry of corresponding OSM nodes and edges.
+
     Parameters
     ----------
-    df: pandas.DataFrame
-        Dataframe with information about edges
-    edges: geopandas.geodataframe.GeoDataFrame
-        The street network, in a projected coordinate reference system
-    g: networkx.graph undirected
-        graph to use for routing
+    df : pandas.DataFrame
+        Dataframe with node pair integer IDs, columns `source` and `target`.
+    edges : geopandas.geodataframe.GeoDataFrame
+        The street network, in a projected CRS.
+    g_undir : networkx.graph undirected
+        Graph to use for routing. The weight variable is "weight".
 
     Returns
     -------
-    df: pandas.DataFrame
-        Dataframe with added path nodes and path edges
+    df : pandas.DataFrame
+        Dataframe with added path nodes and path edges.
     """
     paths = []
     for _, row in df.iterrows():
         paths.append(
             nx.shortest_path(
-                G=g,  # !! use undirected graph here
+                G=g_undir,
                 source=int(row.source),
                 target=int(row.target),
-                weight="length",
+                weight="weight",
             )
         )
     df["path_nodes"] = paths
-    df["path_edges"] = df.path_nodes.apply(lambda x: get_correct_edgetuples(edges, x))
+    df["path_edges"] = df.path_nodes.apply(lambda x: nodelist_to_edgetuples(edges, x))
     return df
 
-
-def create_gdf_with_geoms(df, edges):
-    """
-    Parameters
-    ----------
-    df: pandas.DataFrame
-        Dataframe with path nodes and path edges
-    edges: geopandas.GeoDataFrame
-        The street network, in a projected coordinate reference system
-
-    Returns
-    -------
-    gdf: geopandas.GeoDataFrame
-        projected GeoDataFrame with path nodes and path edges and merged geometries
-    """
-    # get geometry by merging all geoms from edge gdf
-    df["geometry"] = df.path_edges.apply(lambda x: edges.loc[x].geometry.union_all())
-    # convert edges into a gdf
-    gdf = gpd.GeoDataFrame(df, crs=edges.crs, geometry="geometry")
-    # merge multilinestring into linestring where possible (should be possible everywhere)
-    gdf["geometry"] = gdf.line_merge()
-    return gdf
-
-
-
-def map_edges_to_bike_infrastructure(g):
-    """
-    map if edges in graph have bike infrastructure as specified in config_osm.yml
-
-    Parameters
-    ----------
-    g :networkx.MultiDiGraph
-        simplified graph representing the street network
-
-    Returns
-    -------
-    g : networkx.MultiDiGraph
-        simplified graph representing the street network, with added binary edge attribute "pbi"
-    """
-    # first step: map all highway attributes
-    protected_bike_infra = yaml.load(
-        open("./config/config_osm.yml"),
-        Loader=yaml.FullLoader)["protected_bike_infra"]
-
-    # add binary edge attribute "pbi" (protected bike infra: True/False)
-    for edge in g.edges(keys=True):
-        g.edges[edge]["pbi"] = protected_bike_infra[
-            g.edges[edge]["highway"]
-        ]
-    return g
-
-def bike_infra_mapping_gdf(g, edges_gdf):
-    """
-    add binary edge attribute pbi to edges_gdf
-
-    Parameters
-    ----------
-    g : networkx.MultiDiGraph
-        simplified graph representing the street network, with added binary edge attribute "pbi"
-    edges_gdf: geopandas.GeoDataFrame
-        edges representing the street network
-
-    Returns
-    -------
-    edges_gdf: geopandas.GeoDataFrame
-        edges representing the street network with added binary attribute "pbi"
-    """
-    # Build dict of edge attribute
-    attr_dict = {
-        (u, v, k): data['pbi']
-        for u, v, k, data in g.edges(keys=True, data=True)
-    }
-
-    # Map to GeoDataFrame
-    edges_gdf['pbi'] = edges_gdf.index.map(attr_dict)
-    return edges_gdf
-
-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 300 nodes
-    random.seed(1312)
-    random_nodes = random.choices(list(G.nodes), k=300)
-    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
\ No newline at end of file
diff --git a/bikenetlib/viz.py b/bikenetlib/viz.py
index e69de29..da87717 100644
--- a/bikenetlib/viz.py
+++ b/bikenetlib/viz.py
@@ -0,0 +1 @@
+"""Visualization functions of `bikenetlib`."""
\ No newline at end of file
diff --git a/environment-dev.yml b/environment-dev.yml
index 4e7617b..4dd3cba 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -5,7 +5,6 @@ dependencies:
   - python>=3.12
   - matplotlib=3.11.0 # The exact version is important for visualization tests
   - osmnx>=2
-  - scipy
   - geojson
 
   # additional to the non-dev env
diff --git a/environment.yml b/environment.yml
index b14423a..739b777 100644
--- a/environment.yml
+++ b/environment.yml
@@ -5,5 +5,4 @@ dependencies:
   - python>=3.12
   - matplotlib
   - osmnx>=2
-  - scipy
   - geojson
diff --git a/pyproject.toml b/pyproject.toml
index 2374003..f2791f3 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -38,7 +38,6 @@ dependencies = [
     "geopandas>=1",
     "matplotlib",
     "osmnx>=2",
-    "scipy",
     "geojson",
 ]
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index de6f05f..b6b12e3 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -4,71 +4,11 @@
 import geopandas as gpd
 from pandas.testing import assert_frame_equal
 from bikenetlib.utils import ( # List them in the order appearing in utils.py
-    intersects_properly,
     get_principal_bearing,
-    filter_seed_points,
 )
 from shapely.geometry import Point, LineString, MultiLineString
 
 
-# intersects_properly
-
-@pytest.fixture
-def geom_1():
-    linestring = LineString([(0, 0), (1, 1), (2, 2)])
-    return linestring
-
-
-@pytest.fixture
-def geom_2():
-    linestring = LineString([(3, 3), (4, 4), (5, 5)])
-    return linestring
-
-
-def test_intersects_properly(geom_1, geom_2):
-    assert intersects_properly(geom_1, geom_2) is False
-
-
-# filter_seed_points
-
-@pytest.fixture
-def seed_point_delta():
-    return 500
-
-
-@pytest.fixture
-def snapped_seed_points():
-    d = {
-        "osmid": ["1", "2", "3"],
-        "geometry_generated": [Point(1000, 1000), Point(2000, 2000), Point(3000, 3000)],
-    }
-    gdf = gpd.GeoDataFrame(d, geometry="geometry_generated", crs="EPSG:3857")
-    gdf["geometry_osm"] = gpd.GeoSeries(
-        [Point(1001, 1001), Point(10000, 10000), Point(3001, 3001)], crs="EPSG:3857"
-    )
-    return gdf
-
-
-@pytest.fixture
-def filtered_seed_points():
-    d = {"osmid": ["1", "3"], "geometry": [Point(1001, 1001), Point(3001, 3001)]}
-    gdf = gpd.GeoDataFrame(d, geometry="geometry", crs="EPSG:3857")
-    gdf = gdf.set_index("osmid")
-    gdf["osmid"] = gdf.index
-    gdf = gdf.iloc[:, [1, 0]]
-    return gdf
-
-
-def test_filter_seed_points(
-    snapped_seed_points, filtered_seed_points, seed_point_delta
-):
-    assert_frame_equal(
-        filter_seed_points(snapped_seed_points, seed_point_delta),
-        filtered_seed_points,
-        check_dtype=False,
-    )
-
-
 # get_principal_bearing
 
 @pytest.fixture

From 1ed94370ed068b4bbe3045c5966a5477bc1bff10 Mon Sep 17 00:00:00 2001
From: Michael Szell <michael@szell.net>
Date: Sat, 10 Oct 2026 13:49:08 +0200
Subject: [PATCH 7/7] Fix new functions

---
 bikenetlib/osm.py   |  6 ------
 bikenetlib/utils.py | 20 ++++++++++----------
 2 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/bikenetlib/osm.py b/bikenetlib/osm.py
index 71c89f9..ba78d98 100644
--- a/bikenetlib/osm.py
+++ b/bikenetlib/osm.py
@@ -18,10 +18,6 @@
     Dictionary for protected bicycle infrastructure (pbi), with one 
     sub-criterion. For example, highway~path AND bicycle~designated. Always 
     adds a second criterion access!~private.
-ROUTING_PENALTY : dict, default {0: 1.5, 1: 1}
-    Factor to multiply length of non-pbi/pbi for routing, to avoid routing
-    through parallel streets when slightly longer pbi is available. By 
-    default, non-pbi counts as 50% longer than pbi.
 """
 
 # Populate ox.settings.useful_tags_way to make application of custom filter possible
@@ -74,5 +70,3 @@
     ]
 }
 
-ROUTING_PENALTY = {0: 1.5, 1: 1}
-
diff --git a/bikenetlib/utils.py b/bikenetlib/utils.py
index f5b9495..7c651a6 100644
--- a/bikenetlib/utils.py
+++ b/bikenetlib/utils.py
@@ -27,9 +27,9 @@ def assign_pbi_to_edges(g):
         for tag in osm.PBI_DICT:
             if g.edges[edge].get(tag) in osm.PBI_DICT[tag]:
                 assigned_pbi = 1
-        if not assigned_pbi: # Check if the edge fulfils a subcriterion
+        if not assigned_pbi: # Check if the edge fulfills a subcriterion
             for tag in osm.PBI_DICT_SUB:
-                subtag = osm.PBI_DICT_SUB[tag][1].keys()[0]
+                subtag = list(osm.PBI_DICT_SUB[tag][1].keys())[0]
                 if g.edges[edge].get(tag) in osm.PBI_DICT_SUB[tag][0] and g.edges[edge].get(subtag) in osm.PBI_DICT_SUB[tag][1][subtag] and g.edges[edge].get("access") != 'private':
                     assigned_pbi = 1
         g.edges[edge]["pbi"] = assigned_pbi
@@ -37,7 +37,7 @@ def assign_pbi_to_edges(g):
 
     
 def nodelist_to_edgetuples(edge_gdf, nodelist):
-    """Map a list of ndoes (from `nx.shortest_paths()`) to a list of edge 
+    """Map a list of nodes (from `nx.shortest_paths()`) to a list of edge 
     tuples that can be used for indexing an edge geodataframe.
 
     Useful to apply as a lambda function to turn path nodes into path edges.
@@ -180,7 +180,7 @@ def route_nodepairs(nodepairs, edges, g_undir):
 
     Parameters
     ----------
-    df : pandas.DataFrame
+    nodepairs : pandas.DataFrame
         Dataframe with node pair integer IDs, columns `source` and `target`.
     edges : geopandas.geodataframe.GeoDataFrame
         The street network, in a projected CRS.
@@ -189,11 +189,11 @@ def route_nodepairs(nodepairs, edges, g_undir):
 
     Returns
     -------
-    df : pandas.DataFrame
-        Dataframe with added path nodes and path edges.
+    nodepairs : pandas.DataFrame
+        Dataframe of node pairs with added path nodes and path edges.
     """
     paths = []
-    for _, row in df.iterrows():
+    for _, row in nodepairs.iterrows():
         paths.append(
             nx.shortest_path(
                 G=g_undir,
@@ -202,7 +202,7 @@ def route_nodepairs(nodepairs, edges, g_undir):
                 weight="weight",
             )
         )
-    df["path_nodes"] = paths
-    df["path_edges"] = df.path_nodes.apply(lambda x: nodelist_to_edgetuples(edges, x))
-    return df
+    nodepairs["path_nodes"] = paths
+    nodepairs["path_edges"] = nodepairs.path_nodes.apply(lambda x: nodelist_to_edgetuples(edges, x))
+    return nodepairs