高效地将相交的路线矩阵转化为简化的空间网络（图）

Question

我有一个真实详细路线的矩阵，我想将其有效地转变为一个简单的空间网络。简单意味着我不关心当地交通的复杂性以及起点和终点附近可能的路线交叉点。我确实想将起点和终点之外的主要交叉点添加为网络的节点。下面我举一个简单的例子。我的真实数据有 500 个起点和终点之间的 12,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 和 IIUC，您可以在获得可用的 directions

 后使用这种 

primal 方法 :

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 = list(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(
    pd.DataFrame(
        list(combinations(cities["city"], 2)), columns=["city_l", "city_r"]
    ).join(pd.DataFrame(combs_lonlat, columns=["coords_l", "coords_r"])),
    geometry=lines, crs=4326,
)

注意：您需要注册 openrouteservice 才能生成/获取您的 api 密钥。

现在，要制作图，您可以使用

momepy

:

G = gdf_to_nx(
    routes.to_crs(4839), # optional conversion ? 
    multigraph=False, directed=False,
    approach="primal",
)

len(G.nodes) # 20
len(G.edges) # 190

情节（可选）：

ax = routes.plot(color="k")
cities.plot(color="r", marker="x", ax=ax)

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

_ = ax.axis("off")

使用过的进口/输入：

import geopandas as gpd
import openrouteservice
from itertools import combinations
from shapely.geometry import shape
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
    ]
}