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

我有一个包含真实详细路径的数据集，需要高效转换为简化空间网络：

• 简化规则：无需关注起点、终点附近的本地交通细节及路径交叉，但要将起止点外的主要交叉点作为节点加入网络
• 真实数据规模：500个起止点间的12500条路径，大小约2GB

我已经用R语言通过osrm、sfnetworks等工具获取并构建了初始网络，但使用tidygraph的convert(routes_net, to_spatial_subdivision)时，哪怕是测试示例都运行极慢。我知道可以用GRASS的v.clean工具拆分几何图形，但暂时不想安装GRASS。考虑过转换为S2并通过s2_intersection()逐一比较线串来构建网络，但希望有更优雅、高效的解决方案，接受R、Python、QGIS等任何工具。

示例代码（R） #

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)

尝试过的代码（R） #

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

高效解决方案 #

1. R语言：sf空间索引+批量拆分 #

利用sf的空间索引加速相交查询，避免逐一线串比对的低效：

library(sf)
library(data.table)
library(sfnetworks)

# 读取路径数据（替换为你的数据路径/对象）
routes_sf <- st_read("your_routes_data.shp") # 或直接使用已有的routes对象

# 创建空间索引，大幅提升相交查询速度
st_geometry(routes_sf) <- st_geometry(routes_sf)
routes_sf <- st_sf(routes_sf)

# 提取所有原始起止点（需根据你的数据结构调整字段）
orig_dest_points <- st_union(c(
  st_as_sf(routes_sf, coords = c("from_lon", "from_lat")),
  st_as_sf(routes_sf, coords = c("to_lon", "to_lat"))
)$geometry)

# 计算所有线的交点，排除原始起止点
all_lines <- st_union(routes_sf$geometry)
all_intersections <- st_intersection(all_lines, all_lines)
split_points <- st_difference(all_intersections, orig_dest_points)

# 用交点批量拆分所有线
split_routes <- st_split(routes_sf$geometry, split_points)
split_routes_sf <- st_sf(geometry = st_collection_extract(split_routes, "LINESTRING"))

# 转换为简化后的空间网络
simplified_network <- as_sfnetwork(split_routes_sf, directed = FALSE)

2. Python语言：Geopandas+Shapely批量处理 #

Python的空间处理库在大数据量下性能更优，结合空间索引和批量操作：

import geopandas as gpd
from shapely.ops import split, unary_union
import networkx as nx

# 读取路径数据
routes_gdf = gpd.read_file("your_routes_data.shp")

# 创建空间索引
routes_gdf.sindex

# 提取所有原始起止点（根据数据结构调整字段）
origins = gpd.GeoSeries(gpd.points_from_xy(routes_gdf.from_lon, routes_gdf.from_lat))
destinations = gpd.GeoSeries(gpd.points_from_xy(routes_gdf.to_lon, routes_gdf.to_lat))
orig_dest_union = unary_union(origins.union(destinations))

# 计算所有线的交点并排除起止点
all_lines_union = unary_union(routes_gdf.geometry)
all_intersections = all_lines_union.intersection(all_lines_union)
split_points = all_intersections.difference(orig_dest_union)

# 拆分所有线
split_lines = []
for line in routes_gdf.geometry:
    split_parts = split(line, split_points)
    split_lines.extend(list(split_parts))

# 转换为GeoDataFrame并构建网络
split_gdf = gpd.GeoDataFrame(geometry=split_lines, crs=routes_gdf.crs)
# 构建NetworkX网络（如需保留属性可扩展）
simplified_graph = nx.from_edgelist([
    (tuple(line.coords[0]), tuple(line.coords[-1])) 
    for line in split_gdf.geometry
])

3. QGIS可视化操作（无代码） #

适合非编程用户，操作步骤清晰：

1. 导入路径线图层到QGIS
2. 打开处理工具箱，搜索并运行提取交点工具，输入路径图层，得到所有线的交点
3. 运行删除重复要素工具清理交点图层，再用选择要素工具排除属于原始起止点的交点
4. 运行线拆分工具，用清理后的交点图层拆分原始路径线
5. 最后运行构建网络工具，将拆分后的线转换为空间网络，可通过删除重复节点进一步简化

内容的提问来源于stack exchange，提问作者Sebastian