Merge pull request #276 from arup-group/snap-facilities

facility link snapping

Author: Theodore-Chatziioannou

CHANGELOG.md

@@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 
 ### Added
+* Facility link snapping (#276).
 
 ### Changed
 

src/pam/cli.py

@@ -9,6 +9,7 @@
 from pam import read, write
 from pam.operations.combine import pop_combine
 from pam.operations.cropping import simplify_population
+from pam.operations.snap import run_facility_link_snapping
 from pam.report.benchmarks import benchmarks as bms
 from pam.report.stringify import stringify_plans
 from pam.report.summary import pretty_print_summary, print_summary
@@ -667,3 +668,29 @@ def plan_filter(plan):
 
     logger.info("Population wipe complete")
     logger.info(f"Output saved at {path_population_output}")
+
+
+@cli.command()
+@click.argument("path_population_in", type=click.Path(exists=True))
+@click.argument("path_population_out", type=click.Path(exists=False, writable=True))
+@click.argument("path_network_geometry", type=click.Path(exists=True))
+@click.option(
+    "--link_id_field",
+    "-f",
+    type=str,
+    default="id",
+    help="The link ID field to use in the network shapefile. Defaults to 'id'.",
+)
+def snap_facilities(
+    path_population_in: str,
+    path_population_out: str,
+    path_network_geometry: int,
+    link_id_field: Optional[str],
+):
+    """Snap facilities to a network geometry."""
+    run_facility_link_snapping(
+        path_population_in=path_population_in,
+        path_population_out=path_population_out,
+        path_network_geometry=path_network_geometry,
+        link_id_field=link_id_field,
+    )

src/pam/operations/snap.py

@@ -0,0 +1,59 @@
+""" Methods for snapping elements to the network or facilities. """
+
+from pathlib import Path
+
+import geopandas as gp
+import numpy as np
+from scipy.spatial import cKDTree
+
+from pam.core import Population
+from pam.read import read_matsim
+from pam.write import write_matsim
+
+
+def snap_facilities_to_network(
+    population: Population, network: gp.GeoDataFrame, link_id_field: str = "id"
+) -> None:
+    """Snaps activity facilities to a network geometry (in-place).
+
+    Args:
+        population (Population): A PAM population.
+        network (gp.GeoDataFrame): A network geometry shapefile.
+        link_id_field (str, optional): The link ID field to use in the network shapefile. Defaults to "id".
+    """
+    if network.geometry.geom_type[0] == "Point":
+        coordinates = np.array(list(zip(network.geometry.x, network.geometry.y)))
+    else:
+        coordinates = np.array(list(zip(network.geometry.centroid.x, network.geometry.centroid.y)))
+
+    tree = cKDTree(coordinates)
+    link_ids = network[link_id_field].values
+
+    for _, _, person in population.people():
+        for act in person.activities:
+            point = act.location.loc
+            distance, index = tree.query([(point.x, point.y)])
+            act.location.link = link_ids[index[0]]
+
+
+def run_facility_link_snapping(
+    path_population_in: str,
+    path_population_out: str,
+    path_network_geometry: str,
+    link_id_field: str = "id",
+) -> None:
+    """Reads a population, snaps activity facilities to a network geometry, and saves the results.
+
+    Args:
+        path_population_in (str): Path to a PAM population.
+        path_population_out (str): The path to save the output population.
+        path_network_geometry (str): Path to the network geometry file.
+        link_id_field (str, optional): The link ID field to use in the network shapefile. Defaults to "id".
+    """
+    population = read_matsim(path_population_in)
+    if ".parquet" in Path(path_network_geometry).suffixes:
+        network = gp.read_parquet(path_network_geometry)
+    else:
+        network = gp.read_file(path_network_geometry)
+    snap_facilities_to_network(population=population, network=network, link_id_field=link_id_field)
+    write_matsim(population=population, plans_path=path_population_out)

tests/test_29_snap.py

@@ -0,0 +1,41 @@
+import os
+
+import geopandas as gp
+import pytest
+from pam.operations.snap import run_facility_link_snapping, snap_facilities_to_network
+from pam.read import read_matsim
+
+
+def test_add_snapping_adds_link_attribute(population_heh):
+    network = gp.read_file(pytest.test_data_dir / "test_link_geometry.geojson")
+    for _, _, person in population_heh.people():
+        for act in person.activities:
+            assert act.location.link is None
+
+    snap_facilities_to_network(population=population_heh, network=network)
+    for _, _, person in population_heh.people():
+        for act in person.activities:
+            assert act.location.link is not None
+
+            # check that the link is indeed the nearest one
+            link_distance = (
+                network.set_index("id")
+                .loc[act.location.link, "geometry"]
+                .distance(act.location.loc)
+            )
+            min_distance = network.distance(act.location.loc).min()
+            assert link_distance == min_distance
+
+
+def test_links_resnapped(tmpdir):
+    path_out = os.path.join(tmpdir, "pop_snapped.xml")
+    run_facility_link_snapping(
+        path_population_in=pytest.test_data_dir / "1.plans.xml",
+        path_population_out=path_out,
+        path_network_geometry=pytest.test_data_dir / "test_link_geometry.geojson",
+    )
+    assert os.path.exists(path_out)
+    pop_snapped = read_matsim(path_out)
+    for _, _, person in pop_snapped.people():
+        for act in person.activities:
+            assert "link-" in act.location.link

tests/test_data/test_link_geometry.geojson

@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2157" } },
+"features": [
+{ "type": "Feature", "properties": { "id": "link-1" }, "geometry": { "type": "LineString", "coordinates": [ [ 10000.0, 5000.0 ], [ 10000.0, 10000.0 ] ] } },
+{ "type": "Feature", "properties": { "id": "link-2" }, "geometry": { "type": "LineString", "coordinates": [ [ 10000.0, 10000.0 ], [ 5000.0, 5000.0 ] ] } }
+]
+}