Calculating shortest path between two points along a rail network

Question

I would like to calculate the shortest path from mines to ports, along the rail network. Since the mines and ports are not directly on the network, I first worked out the nearest rail point for each mine and port.

However, I am now struggling to convert the rail network into a NetworkX graph and from there, calculating the shortest distance between each mine and port. When I ran my code, it is unable to work out the nearest path since all of the distances are zero. Please see below for the code I have written thus far. The rail network can be downloaded from here. I would kindly appreciate any assistance that can be provided. Thank you.

# Ports
ports = 
[{'port_name': 'Geraldton', 'geometry': POINT (114.59786 -28.77688)},
{'port_name': 'Bunbury', 'geometry': POINT(115.673447 -33.318797)},
{'port_name': 'Albany', 'geometry': POINT(117.895025 -35.032831)}, 
{'port_name': 'Esperance', 'geometry': POINT(121.897114 -33.871834)}]  

# Mine sites
mines = 
[{'mine_name': 'Gold', 'geometry': POINT (117.94568 -34.93467), 
{'mine_name': 'Silver', 'geometry': POINT (115.16923 -29.65613)},
{'mine_name': 'Bronze', 'geometry': POINT (115.11039 -29.51287)}, 
{'mine_name': 'Platinum', 'geometry': POINT (115.11130 -29.42621)}]

# Convert to GeoDataFrame
crs = 'EPSG:4326'
ports = gpd.GeoDataFrame(ports, crs=crs, geometry = 'geometry')
mine_sites = gpd.GeoDataFrame(ports, crs=crs, geometry = 'geometry')

# Column for nearest rail point
ports['nearest_rail_point'] = None
mine_sites['nearest_rail_point'] = None

# Nearest rail point
for index, row in ports.iterrows():
    port_point = row['geometry']
    nearest_point = nearest_points(port_point, rail_network.unary_union)[1] 
    ports.at[index, 'nearest_rail_point'] = nearest_point

for index, row in mine_sites.iterrows():
    mine_site_point = row['geometry']
    nearest_point = nearest_points(mine_site_point, rail_network.unary_union)[1] 
    mine_sites.at[index, 'nearest_rail_point'] = nearest_point

G = momepy.gdf_to_nx(rail_network)

shortest_distances = {}

# Calculate shortest path
for mine_index, mine_row in mines.iterrows():
    mine_name = mine_row['mine_name']
    mine_point = mine_row['geometry']
    
    shortest_distances[mine_name] = {}
    
    for port_index, port_row in ports.iterrows():
        port_name = port_row['port_name']
        port_point = port_row['geometry']

        try:
            shortest_distance = nx.shortest_path_length(G, 
                                                        source=(mine_point.x, mine_point.y), 
                                                        target=(port_point.x, port_point.y), 
                                                        weight='length')
            shortest_distances[mine_name][port_name] = shortest_distance
        except nx.NetworkXNoPath:
            shortest_distances[mine_name][port_name] = float('inf')  # No path found

Answer 1

It's unclear if you need the shortest path between each mine and its corresponding port (at the same index) or among all ports. Either way, one option would be to use a primal approach (with momempy) to create the graph and compute the shortest_path_length (or shortest_path, depending on your expected output). However, before that, you need to snap the mines to the nearest rail boundary :

spots = (pd.concat([
    ports.rename(columns={"port_name": "spot_name"}),
    mine_sites.rename(columns={"mine_name": "spot_name"})])
             .to_crs("EPSG:3857"))

import momepy
G = (momepy.gdf_to_nx(
    rail_network[["name", "geometry"]].explode(),
    approach="primal", multigraph=False, length="dis"))

def nearest_bnd(p, ls, shp=False):
    sp, *_, ep = ls.coords
    nb = min(map(Point, [sp, ep]), key=p.distance)
    return nb if shp else (nb.x, nb.y)
        
coos = (spots.rename_geometry("geom_sp")
        .sjoin_nearest(rail_network.assign(geom_rn=rail_network["geometry"]))
        .set_geometry("geom_rn").assign(coo=lambda x: [nearest_bnd(p, ls)
        for p, ls in x[["geom_sp", "geom_rn"]].to_numpy()])
        .set_index("spot_name")["coo"])

all_spl = {}
for mn in mine_sites["mine_name"]:
    for pn in ports["port_name"]:
        try:
            uv = map(coos.get, [mn, pn])
            spl = nx.shortest_path_length(G, *uv, weight="dis")
            spl = round(spl / 1000, 2)
        except nx.NetworkXNoPath:
            spl = float("inf")
        all_spl.setdefault(mn, []).append({pn: spl})

---

{mn: min(d, key=lambda x: list(x.values())[0]) for mn, d in all_spl.items()}

# {
#     "Gold": {"Albany": 0.0},
#     "Silver": {"Geraldton": 117.85},
#     "Bronze": {"Geraldton": 117.85},
#     "Platinum": {"Geraldton": 117.85},
# }