高效将交叉路线矩阵转换为简化空间网络（图）

Question

我有一个详细的路线矩阵，想把它高效地转化成一个简单的空间网络。这里的“简单”是指我不太关心起点和终点附近的复杂交通情况和可能的交叉路口。我希望在网络中添加一些主要的交叉口作为节点，这些交叉口不在起点和终点之间。下面我给了一个简单的例子。我的真实数据有12,500条路线，连接500个起点和终点，大小大约是2GB。

library(fastverse)
#> -- Attaching packages --------------------------------------- fastverse 0.3.2 --
#> v data.table 1.15.0     v kit        0.0.13
#> v magrittr   2.0.3      v collapse   2.0.12
fastverse_extend(osrm, sf, sfnetworks, install = TRUE)
#> -- Attaching extension packages ----------------------------- fastverse 0.3.2 --
#> v osrm       4.1.1      v sfnetworks 0.6.3 
#> v sf         1.0.16

largest_20_german_cities <- data.frame(
  city = c("Berlin", "Stuttgart", "Munich", "Hamburg", "Cologne", "Frankfurt",
    "Duesseldorf", "Leipzig", "Dortmund", "Essen", "Bremen", "Dresden",
    "Hannover", "Nuremberg", "Duisburg", "Bochum", "Wuppertal", "Bielefeld", "Bonn", "Muenster"),
  lon = c(13.405, 9.18, 11.575, 10, 6.9528, 8.6822, 6.7833, 12.375, 7.4653, 7.0131,
          8.8072, 13.74, 9.7167, 11.0775, 6.7625, 7.2158, 7.1833, 8.5347, 7.1, 7.6256),
  lat = c(52.52, 48.7775, 48.1375, 53.55, 50.9364, 50.1106, 51.2333, 51.34, 51.5139,
          51.4508, 53.0758, 51.05, 52.3667, 49.4539, 51.4347, 51.4819, 51.2667, 52.0211, 50.7333, 51.9625))

# Unique routes
m <- matrix(1, 20, 20)
diag(m) <- NA
m[upper.tri(m)] <- NA
routes_ind <- which(!is.na(m), arr.ind = TRUE)
rm(m)

# Routes DF
routes <- data.table(from_city = largest_20_german_cities$city[routes_ind[, 1]], 
                     to_city = largest_20_german_cities$city[routes_ind[, 2]], 
                     duration = NA_real_, 
                     distance = NA_real_, 
                     geometry = list())
# Fetch Routes
i = 1L
for (r in mrtl(routes_ind)) {
  route <- osrmRoute(ss(largest_20_german_cities, r[1], c("lon", "lat")),
                     ss(largest_20_german_cities, r[2], c("lon", "lat")), overview = "full")
  set(routes, i, 3:5, fselect(route, duration, distance, geometry))
  i <- i + 1L
}
routes %<>% st_as_sf(crs = st_crs(route))

routes_net = as_sfnetwork(routes, directed = FALSE)
print(routes_net)
#> # A sfnetwork with 20 nodes and 190 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # An undirected simple graph with 1 component with spatially explicit edges
#> #
#> # A tibble: 20 × 1
#>              geometry
#>           <POINT [°]>
#> 1  (9.179999 48.7775)
#> 2      (13.405 52.52)
#> 3 (11.57486 48.13675)
#> 4 (10.00001 53.54996)
#> 5   (6.95285 50.9364)
#> 6   (8.68202 50.1109)
#> # ℹ 14 more rows
#> #
#> # A tibble: 190 × 7
#>    from    to from_city to_city duration distance                       geometry
#>   <int> <int> <chr>     <chr>      <dbl>    <dbl>               <LINESTRING [°]>
#> 1     1     2 Stuttgart Berlin      390.     633. (9.179999 48.7775, 9.18005 48…
#> 2     2     3 Munich    Berlin      356.     586. (11.57486 48.13675, 11.57486 …
#> 3     2     4 Hamburg   Berlin      176.     288. (10.00001 53.54996, 10.0002 5…
#> # ℹ 187 more rows
plot(routes_net)

^{创建于2024-03-28，使用了 reprex v2.0.2}

关于可能的解决方案，我对任何软件都持开放态度（比如R、Python、QGIS等）。我知道在R中有一个叫 tidygraph 的工具，可以让我做一些类似的事情

library(tidygraph)
routes_net_subdiv = convert(routes_net, to_spatial_subdivision)

但是即使是这个简单的例子，它运行起来似乎也要花很长时间。我还看到有人建议使用 GRASS的v.clean工具 来拆分几何形状，但我还没有尝试过，而且有点不想安装GRASS。

我觉得也许最好的性能解决方案是转换为S2格式，然后使用 s2_intersection() 分别比较所有的线段，最后把这些信息转化成图形。不过我希望能找到更优雅和高效的解决方案。

Answer 1

我对任何软件（R、Python、QGIS等）都很开放。

使用Python，如果我理解正确，你可以在获取可用的directions之后，使用这个基本的方法：

cities = gpd.GeoDataFrame(
    pd.DataFrame(data), crs="EPSG:4326",
    geometry=gpd.points_from_xy(data["lon"], data["lat"])
)

client = openrouteservice.Client(key="") # put here your key

def get_route(coords, as_ls=True):
    res = client.directions(coords, profile="driving-car", format="geojson")
    geom = res["features"][0]["geometry"]
    return shape(geom) if as_ls else geom

combs_lonlat = combinations(cities[["lon", "lat"]].agg(tuple, axis=1), 2)

# maybe we should use a standard loop and add time.sleep(n) ?
lines = [get_route(c) for c in combs_lonlat] # will eventually show warnings

routes = gpd.GeoDataFrame(geometry=lines, crs=4326)

注意：你需要在openrouteservice上注册，以便生成/获取你的API密钥。

现在，为了创建图，你可以使用momepy：

splines = set()
for ll, rl in combinations(routes["geometry"], 2):
    try:
        for line in split(ll, rl).geoms:
            splines.add(line)
    except ValueError:
        pass

G = gdf_to_nx(
    gpd.GeoDataFrame(geometry=list(splines)),
    multigraph=False, directed=False,
    approach="primal",
)

len(G.nodes) # 250
len(G.edges) # 1296

绘图（可选）：

fig, (axl, axr) = plt.subplots(1, 2)

for ax, inp in zip((axl, axr), (cities, G)):
    routes.plot(color="k", ax=ax)
    cities.plot(color="r", marker="x", zorder=3, ax=ax)
    _ = ax.set_title(f"{len(inp)} nodes", fontweight="bold")
    _ = ax.axis("off")

for x, y, name in zip(
    cities["geometry"].x,
    cities["geometry"].y,
    cities["city"]
):
    for ax in (axl, axr):
        _ = ax.annotate(
            name, xy=(x, y), xytext=(3, 3), textcoords="offset points"
        )

(gpd.GeoSeries([Point(n) for n in G.nodes], crs=4326)
     .plot(color="b", marker="o", ax=axr))

使用的导入/输入：

import pandas as pd
import geopandas as gpd
import openrouteservice
import matplotlib.pyplot as plt
from itertools import combinations
from shapely.geometry import shape, Point
from shapely.ops import split
from momepy import gdf_to_nx

data = {
    "city": [
        "Berlin", "Stuttgart", "Munich", "Hamburg", "Cologne",
        "Frankfurt", "Duesseldorf", "Leipzig", "Dortmund", "Essen",
        "Bremen", "Dresden", "Hannover", "Nuremberg", "Duisburg",
        "Bochum", "Wuppertal", "Bielefeld", "Bonn", "Muenster"
    ],
    "lon": [
        13.405, 9.18, 11.575, 10, 6.9528, 8.6822, 6.7833, 12.375,
        7.4653, 7.0131, 8.8072, 13.74, 9.7167, 11.0775, 6.7625, 7.2158,
        7.1833, 8.5347, 7.1, 7.6256
    ],
    "lat": [
        52.52, 48.7775, 48.1375, 53.55, 50.9364, 50.1106, 51.2333, 51.34,
        51.5139, 51.4508, 53.0758, 51.05, 52.3667, 49.4539, 51.4347, 51.4819,
        51.2667, 52.0211, 50.7333, 51.9625
    ]
}