add graph from travel network (#755)

* add graph from travel network

* typo graphsummary

* add network tests

* test_method

* skip net test on windows

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* format

* picky ass linter

* distance-->cost

* ugh. lint.

* docstring; test length

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* validate geometry; remove unused options

* validate before

* imports

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* document node_ids param

* add example in docstring

* Update libpysal/graph/base.py

Co-authored-by: Martin Fleischmann <[email protected]>

* Update libpysal/graph/base.py

Co-authored-by: Martin Fleischmann <[email protected]>

* Update libpysal/graph/base.py

Co-authored-by: Martin Fleischmann <[email protected]>

* Update libpysal/graph/base.py

Co-authored-by: Martin Fleischmann <[email protected]>

* continue formatting

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* Update libpysal/graph/_network.py

Co-authored-by: James Gaboardi <[email protected]>

* curlies if were being picky; im gonna fight this linter--this one was your fault james!!! :P

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Martin Fleischmann <[email protected]>
Co-authored-by: James Gaboardi <[email protected]>

ci/310-latest.yaml

@@ -31,6 +31,7 @@ dependencies:
   - sqlalchemy
   - xarray
   - zstd
+  - pandana
   - pip
   - pip:
       - pulp

ci/310-oldest.yaml

@@ -29,6 +29,7 @@ dependencies:
   - sqlalchemy=2.0
   - zstd
   - xarray=2022.3
+  - pandana
   - pip
   - pip:
       - pulp

ci/311-latest.yaml

@@ -31,6 +31,7 @@ dependencies:
   - sqlalchemy
   - xarray
   - zstd
+  - pandana
   - pip
   - pip:
       - pulp

ci/312-dev.yaml

@@ -28,6 +28,7 @@ dependencies:
   - pyproj
   - sqlalchemy
   - zstd
+  - pandana
   - pip
   - pip:
       # dev versions of packages

ci/312-latest.yaml

@@ -30,6 +30,7 @@ dependencies:
   - scikit-learn
   - sqlalchemy
   - zstd
+  - pandana
   - xarray
   # for docs build action (this env only)
   - myst-parser

libpysal/graph/_network.py

@@ -0,0 +1,119 @@
+import numpy as np
+
+from ._utils import _induce_cliques, _validate_geometry_input
+
+
+def _build_coplanarity_node_lookup(geoms):
+    """
+    Identify coplanar points and create a look-up table for the coplanar geometries.
+    Same function as in ``graph._utils``, but need to keep the index to use as graph ids
+    """
+    # geoms = geoms.reset_index(drop=True)
+    coplanar = []
+    nearest = []
+    r = geoms.groupby(geoms).groups
+    for g in r.values():
+        if len(g) == 2:
+            coplanar.append(g[0])
+            nearest.append(g[1])
+        elif len(g) > 2:
+            for n in g[1:]:
+                coplanar.append(n)
+                nearest.append(g[0])
+    return np.asarray(coplanar), np.asarray(nearest)
+
+
+def pdna_to_adj(origins, network, node_ids, threshold):
+    """Create an adjacency list of shortest network-based travel between
+       origins and destinations in a pandana.Network.
+
+    Parameters
+    ----------
+    origins : geopandas.GeoDataFrame
+        Geodataframe of origin geometries to begin routing.
+    network : pandana.Network
+        pandana.Network instance that stores the local travel network
+    node_ids:
+        array of node_ids in the pandana.Network aligned with the input
+        observations in ``origins``. This is created via a call like
+        ``pandana.Network.get_node_ids(df.geometry.x, df.geometry.y)``
+    threshold : int
+        maximum travel distance (inclusive)
+
+    Returns
+    -------
+    pandas.DataFrame
+        adjacency list with columns 'origin', 'destination', and 'cost'
+    """
+
+    # map node ids in the network to index in the gdf
+    mapper = dict(zip(node_ids, origins.index.values, strict=False))
+
+    namer = {"source": "origin", network.impedance_names[0]: "cost"}
+
+    adj = network.nodes_in_range(node_ids, threshold)
+    adj = adj.rename(columns=namer)
+    # swap osm ids for gdf index
+    adj = adj.set_index("destination").rename(index=mapper).reset_index()
+    adj = adj.set_index("origin").rename(index=mapper).reset_index()
+    adj = adj[adj.destination.isin(origins.index.values)]
+
+    return adj
+
+
+def build_travel_graph(
+    df,
+    network,
+    threshold,
+):
+    """Compute the shortest path between gdf centroids via a pandana.Network
+    and return an adjacency list with weight=cost. Note unlike distance_band,
+    :math:`G_{ij}` and :math:`G_{ji}` are often different because travel networks
+    may be directed.
+
+    Parameters
+    ----------
+    df : geopandas.GeoDataFrame
+        geodataframe of observations. CRS should be the same as the locations
+        of ``node_x`` and ``node_y`` in the pandana.Network (usually 4326 if network
+        comes from OSM, but sometimes projected to improve snapping quality).
+    network : pandana.Network
+        Network that encodes travel costs. See <https://udst.github.io/pandana/>
+    threshold : int
+        maximum travel cost to consider neighbors
+
+    Returns
+    -------
+    pandas.Series
+        adjacency formatted as multiindexed (focal, neighbor) series
+    """
+    df = df.copy()
+    _validate_geometry_input(df.geometry, ids=None, valid_geometry_types="Point")
+    df["node_ids"] = network.get_node_ids(df.geometry.x, df.geometry.y)
+
+    # depending on density of the graph nodes / observations, it is common to have
+    # multiple observations snapped to the same network node, so use the clique
+    # expansion logic to handle these cases
+
+    # get indices of multi-observations at unique nodes
+    coplanar, nearest = _build_coplanarity_node_lookup(df["node_ids"])
+    # create adjacency on unique nodes
+    adj = pdna_to_adj(df, network, df["node_ids"].values, threshold)
+    # add clique members back to adjacency
+    adj_cliques = _induce_cliques(
+        adj.rename(
+            columns={"origin": "focal", "destination": "neighbor", "cost": "weight"}
+        ),
+        coplanar=coplanar,
+        nearest=nearest,
+    )
+    # reorder, drop induced dupes, and return
+    adj_cliques = (
+        adj_cliques.groupby(["focal", "neighbor"])
+        .first()
+        .reindex(df.index, level=0)
+        .reindex(df.index, level=1)
+        .reset_index()
+    )
+
+    return adj_cliques

libpysal/graph/base.py

@@ -19,6 +19,7 @@
 from ._indices import _build_from_h3
 from ._kernel import _distance_band, _kernel
 from ._matching import _spatial_matching
+from ._network import build_travel_graph as _build_travel_graph
 from ._plotting import _explore_graph, _plot
 from ._raster import _generate_da, _raster_contiguity
 from ._set_ops import SetOpsMixin
@@ -1454,6 +1455,89 @@ def build_h3(cls, ids, order=1, weight="distance"):
         else:
             raise ValueError("weight must be one of 'distance', 'binary', or 'inverse'")
 
+    @classmethod
+    def build_travel_cost(cls, df, network, threshold, kernel=None):
+        """Generate a Graph based on shortest travel costs from a pandana.Network
+
+        Parameters
+        ----------
+        df : geopandas.GeoDataFrame
+            geodataframe representing observations which are snapped to the nearest
+            node in the pandana.Network. If passing polygon geometries, the spatial
+            support will be reduced to Points (via centroid) before snapping.
+        network : pandana.Network
+            pandana Network object describing travel costs between nodes in the study
+            area
+        threshold : int
+            threshold representing maximum cost distances. This is measured in the same
+            units as the pandana.Network (not influenced by the df.crs in any way). For
+            travel modes with relatively constant speeds like walking or biking, this is
+            usually distance (e.g. meters if the Network is constructed from OSM). For a
+            a multimodal or auto network with variable travel speeds, this is usually
+            some measure of travel time
+        kernel : str or callable, optional
+            kernel transformation applied to the weights. See
+            libpysal.graph.Graph.build_kernel for more information on kernel
+            transformation options. Default is None, in which case the Graph weight
+            is pure distance between focal and neighbor
+
+        Returns
+        -------
+        Graph
+
+        Examples
+        ---------
+        >>> import geodatasets
+        >>> import geopandas as gpd
+        >>> import osmnx as ox
+        >>> import pandana as pdna
+
+        Read an example geodataframe:
+
+        >>> df = gpd.read_file(geodatasets.get_path("geoda Cincinnati")).to_crs(4326)
+
+        Download a walk network using osmnx
+
+        >>> osm_graph = ox.graph_from_polygon(df.union_all(), network_type="walk")
+        >>> nodes, edges = ox.utils_graph.graph_to_gdfs(osm_graph)
+        >>> edges = edges.reset_index()
+
+        Generate a routable pandana network from the OSM nodes and edges
+
+        >>> network = pdna.Network(
+        >>>     edge_from=edges["u"],
+        >>>     edge_to=edges["v"],
+        >>>     edge_weights=edges[["length"]],
+        >>>     node_x=nodes["x"],
+        >>>     node_y=nodes["y"],)
+
+        Use the pandana network to compute shortest paths between gdf centroids and
+        generate a Graph
+
+        >>> G = Graph.build_travel_cost(df.set_geometry(df.centroid), network, 500)
+        >>> G.adjacency.head()
+        focal  neighbor
+        0       62          385.609009
+                65          309.471985
+                115         346.858002
+                116           0.000000
+                117         333.639008
+        Name: weight, dtype: float64
+        """
+        adj = _build_travel_graph(df, network, threshold)
+        g = cls.from_adjacency(adj)
+        if kernel is not None:
+            arrays = _kernel(
+                g.sparse,
+                metric="precomputed",
+                kernel=kernel,
+                bandwidth=threshold,
+                resolve_isolates=False,
+                ids=df.index.values,
+            )
+            return cls.from_arrays(*arrays)
+        return g
+
     @cached_property
     def neighbors(self):
         """Get neighbors dictionary

libpysal/graph/tests/test_builders.py

@@ -1,3 +1,5 @@
+import sys
+
 import geodatasets
 import geopandas as gpd
 import numpy as np
@@ -495,3 +497,65 @@ def test_matching_strids(self):
         assert pd.api.types.is_string_dtype(g._adjacency.index.dtypes["focal"])
         assert pd.api.types.is_string_dtype(g._adjacency.index.dtypes["neighbor"])
         assert pd.api.types.is_numeric_dtype(g._adjacency.dtype)
+
+
+@pytest.mark.network
+@pytest.mark.skipif(
+    sys.platform.startswith("win"), reason="pandana has dtype issues on windows"
+)
+class TestTravelNetwork:
+    def setup_method(self):
+        pandana = pytest.importorskip("pandana")
+        import pooch
+
+        self.net_path = pooch.retrieve(
+            "https://spatial-ucr.s3.amazonaws.com/osm/metro_networks_8k/17140.h5",
+            known_hash=None,
+        )
+        df = gpd.read_file(geodatasets.get_path("geoda cincinnati")).to_crs(4326)
+        self.df = df.set_geometry(df.centroid)
+        self.network = pandana.Network.from_hdf5(self.net_path)
+
+    def test_build_travel_network(self):
+        g = graph.Graph.build_travel_cost(self.df, self.network, 500)
+        assert_array_almost_equal(
+            g.adjacency.head(10).to_numpy(),
+            np.array(
+                [
+                    418.28601074,
+                    228.23899841,
+                    196.0269928,
+                    0.0,
+                    341.73498535,
+                    478.47799683,
+                    298.91699219,
+                    445.60501099,
+                    174.64199829,
+                    0.0,
+                ]
+            ),
+        )
+        assert g.n == self.df.shape[0]
+
+    def test_build_travel_network_kernel(self):
+        g = graph.Graph.build_travel_cost(
+            self.df, self.network, 500, kernel="triangular"
+        )
+        assert_array_almost_equal(
+            g.adjacency.head(10).to_numpy(),
+            np.array(
+                [
+                    0.16342798,
+                    0.543522,
+                    0.60794601,
+                    1.0,
+                    0.31653003,
+                    0.04304401,
+                    0.40216602,
+                    0.10878998,
+                    0.650716,
+                    1.0,
+                ]
+            ),
+        )
+        assert g.n == self.df.shape[0]