Search code examples
python-3.xroutesnetworkxgeopandasosmnx

Networkx Shortest Path Analysis on multiple source and target nodes


What I have:

  • A geodataframe of school points (source - a total of 18)
  • A geodataframe of hospital pts (target - a total of 27)
  • A projeceted Osmnx graph (nodes + edges)

What I want:

  • A geodataframe containing the shortest route geometries to each hospital from each school (a total of 486 [18*27] features in the table each with a route) ie
school id hospital id route
xxxxxxxxx xxxxxxxxxxx LineString(x,x,x)

After reading in schools/hospitals, pulling and projecting the osmnx street graph

I am able to define a function to get the neareset osm node for both source and target points

# import neceessary modules
import pandas as pd
import geopandas as gpd
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
from pyproj import CRS
from shapely.geometry import Polygon, Point, LineString
)

# read in file
hosp_fp = r'vtData/hospitals.shp'
school_fp = r'vtData/schools.shp'
​
# read files
hospitals = gpd.read_file(hosp_fp)
schools = gpd.read_file(school_fp)

#### skip the reading osmnx part (the error isn't here and, yes,everything in same crs) ######

# Create function to find nearest node
def get_nearest_node(points, graph):
    # a function that finds the nearest node
    # params: points (a gdf of points with an x and y column); graph (an osm network graph)
    points['nearest_osm'] = None
    for i in tqdm(points.index, desc='find nearest osm node from input points', position=0):
        points.loc[i, 'nearest_osm'] = ox.get_nearest_node(graph, [points.loc[i, 'y'], points.loc[i, 'x']], method='euclidean') # find the nearest node from hospital location
    return(points)

# use the function to find each destination point nearest node
## returns the original gdfs with corresponding osmid column
source = get_nearest_node(schools, graph)
target = get_nearest_node(hospitals, graph)

# extract osmid's from list
src_list = list(source['nearest_osm'])
trg_list = list(target['nearest_osm'])

### WHERE I AM STUCK ####
# a function meant to construct shortest path routes to each target from each source

def get_routes(graph, src_list, trg_list):
    # a function that constructs a shortest routes to each target from the source
    # params: graph_proj (a projected osmnx graph); src (source pts); trg (target pts)

    # an empty list to append route linestring geometries to
    routes = []

    # a loop to construct all shortest path geometries
    for src, trg in zip(src_list, trg_list):
        sp = nx.shortest_path(graph, source=src, target=trg, 
        weight='length')
        routes.append(sp)
        print(len(routes))

Instead of returning 486 routes (one for each source and target), I am only getting a list of of 18 points (basically, it is only calculating routes based on the corresponding indices of the source and target points, rather than calculating 27 shortest routes (total num hospitals) for each school

From here, I would append the list into a new geodataframe called routes, but I only have 18 of my 486 routes


Solution

  • You are looking for the cartesian product of your origins and destinations, rather than zipping them together. Example:

    import numpy as np
    import osmnx as ox
    from itertools import product
    ox.config(log_console=True)
    
    # get a graph and add edge travel times
    G = ox.graph_from_place('Piedmont, CA, USA', network_type='drive')
    G = ox.add_edge_travel_times(ox.add_edge_speeds(G))
    
    # randomly choose 10 origins and 10 destinations
    n = 10
    origs = np.random.choice(G.nodes, size=n)
    dests = np.random.choice(G.nodes, size=n)
    
    # calculate 100 (10 origins x 10 destinations) shortest paths
    paths = []
    for o, d in product(origs, dests):
        path = ox.shortest_path(G, o, d, weight='travel_time')
        paths.append(path)
        
    len(paths) #100