批量使用R geosphere::distVincentyEllipsoid计算图加权边长度报错求助

我需要创建一个包含顶点和边的图，边权重为大圆距离。真实数据有14316个顶点，直接创建会出现Error: cannot allocate vector of size 21858.7 Gb的内存错误，因此尝试用批量方式处理。

示例顶点数据 #

Lon <- c(179.5,  31.5,  32.5,  33.5,  34.5,  30.5,  31.5,  32.5,  33.5,  34.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5)
Lat <- c(-16.402778,  -1.402778,  -1.402778,  -1.402778,  -1.402778,  -2.402778,  -2.402778,  -2.402778,  -2.402778,  -2.402778, -3.402778,  -3.402778,  -3.402778,  -3.402778,  -3.402778,  -3.402778)
Elevation <- c(99, 1194, 1134, 1134, 1358, 1425, 1307, 1135, 1134, 1452, 1206, 1199, 1194, 1357, 1478, 1333)
vertices <- data.frame(Lon, Lat, Elevation)

运行后顶点数据预览：

Lon        Lat Elevation
1  179.5 -16.402778        99
2   31.5  -1.402778      1194
3   32.5  -1.402778      1134
4   33.5  -1.402778      1134
5   34.5  -1.402778      1358
6   30.5  -2.402778      1425
7   31.5  -2.402778      1307
8   32.5  -2.402778      1135
9   33.5  -2.402778      1134
10  34.5  -2.402778      1452
11  31.5  -3.402778      1206
12  32.5  -3.402778      1199
13  33.5  -3.402778      1194
14  34.5  -3.402778      1357
15  35.5  -3.402778      1478
16  36.5  -3.402778      1333

使用igraph和geosphere包，定义calculate_distances函数通过distVincentyEllipsoid计算大圆距离，批量处理时出现错误：Error in out[, i] <- r : number of items to replace is not a multiple of replacement length。

原批量处理代码 #

library(igraph)
library(geosphere)

# Create an empty graph
g <- graph.empty()

# Define the batch size (adjust as needed)
batch_size <- 4

# Calculate the total number of batches
num_batches <- ceiling(nrow(vertices) / batch_size)

# Function to calculate great circle distances
calculate_distances <- function(vertices1, vertices2) {
  # Extract the Lon and Lat columns for each set of vertices
  coords1 <- vertices1[, c("Lon", "Lat")]
  coords2 <- vertices2[, c("Lon", "Lat")]
  
  # Convert coordinates to a matrix with numeric values
  coords1 <- as.matrix(coords1)
  coords2 <- as.matrix(coords2)
  
  # Calculate distances using distVincentyEllipsoid
  distances <- distVincentyEllipsoid(coords1, coords2,
                                     a = 6378137, b = 6356752.3142, f = 1/298.257223563)
  
  return(distances)
}

# Iterate through batches
for (i in 1:num_batches) {
  # Determine the start and end indices for the current batch
  start_index <- (i - 1) * batch_size + 1
  end_index <- min(i * batch_size, nrow(vertices))
  
  # Subset the vertices for the current batch
  batch_vertices <- vertices[start_index:end_index, c("Lon", "Lat")]
  
  # Add vertices to the existing graph and get their indices
  batch_vertex_indices <- add.vertices(g, nv = nrow(batch_vertices))
  
  # Calculate great circle distances for the current batch
  batch_distances <- calculate_distances(batch_vertices, vertices)
  
  # Generate combinations of vertex pairs within the current batch
  edge_combinations <- combn(batch_vertex_indices, 2)
  
  # Ensure edge combinations and weights match correctly
  num_edges <- ncol(edge_combinations)
  num_weights <- length(batch_distances)
  
  # Adjust the number of weights to match the number of edges
  adjusted_weights <- rep(batch_distances, each = num_edges)
  
  # Trim excess weights if necessary
  adjusted_weights <- adjusted_weights[1:num_edges]
  
  # Add edges to the existing graph with edge distances
  g <- add_edges(g, edges = edge_combinations, edge.attr = adjusted_weights)
}

错误原因分析 #

1. 距离计算结果维度不匹配：distVincentyEllipsoid返回的是nrow(coords1)*nrow(coords2)的矩阵，原代码错误地将其当作一维向量处理，后续权重调整逻辑完全混乱。
2. 边创建逻辑错误：原代码仅创建批量内部的顶点对，但未明确是否是需求；同时add.vertices返回的新顶点ID管理不当，容易导致边对索引错误。
3. 内存逻辑冗余：calculate_distances(batch_vertices, vertices)计算了批量到所有顶点的距离，但后续未正确使用这些值，造成资源浪费。

修正后的代码 #

场景1：创建全连接图（任意两点间都有边） #

library(igraph)
library(geosphere)

# 先添加全部顶点，避免多次add.vertices的索引混乱
g <- graph_from_data_frame(directed = FALSE, vertices = vertices)

batch_size <- 4
num_vertices <- nrow(vertices)
num_batches <- ceiling(num_vertices / batch_size)

# 计算大圆距离的函数，返回矩阵
calculate_distances <- function(coords1, coords2) {
  distVincentyEllipsoid(as.matrix(coords1), as.matrix(coords2),
                        a = 6378137, b = 6356752.3142, f = 1/298.257223563)
}

# 批量处理：新批量顶点与所有已存在顶点创建边
for (i in 2:num_batches) {
  start_idx <- (i-1)*batch_size + 1
  end_idx <- min(i*batch_size, num_vertices)
  current_batch_idx <- start_idx:end_idx
  existing_idx <- 1:(start_idx-1)
  
  # 获取坐标
  current_coords <- vertices[current_batch_idx, c("Lon", "Lat")]
  existing_coords <- vertices[existing_idx, c("Lon", "Lat")]
  
  # 计算距离矩阵
  dist_matrix <- calculate_distances(current_coords, existing_coords)
  
  # 生成边对并提取对应权重
  edge_pairs <- expand.grid(from = existing_idx, to = current_batch_idx)
  edge_weights <- as.vector(t(dist_matrix))
  
  # 添加边到图中
  g <- add_edges(g, edges = as.vector(t(edge_pairs)), attr = list(weight = edge_weights))
}

# 处理每个批量内部的边
for (i in 1:num_batches) {
  start_idx <- (i-1)*batch_size + 1
  end_idx <- min(i*batch_size, num_vertices)
  batch_idx <- start_idx:end_idx
  
  if (length(batch_idx) >=2) {
    edge_pairs <- combn(batch_idx, 2)
    coords_batch <- vertices[batch_idx, c("Lon", "Lat")]
    dist_matrix <- calculate_distances(coords_batch, coords_batch)
    # 提取下三角距离（对应combn的无序对）
    edge_weights <- dist_matrix[lower.tri(dist_matrix)]
    
    g <- add_edges(g, edges = as.vector(edge_pairs), attr = list(weight = edge_weights))
  }
}

场景2：仅创建批量内部的边 #

如果需求是每个批量内部的顶点互相连接，使用以下代码：

library(igraph)
library(geosphere)

g <- graph.empty()

batch_size <- 4
num_vertices <- nrow(vertices)
num_batches <- ceiling(num_vertices / batch_size)

calculate_distances <- function(coords1, coords2) {
  distVincentyEllipsoid(as.matrix(coords1), as.matrix(coords2),
                        a = 6378137, b = 6356752.3142, f = 1/298.257223563)
}

for (i in 1:num_batches) {
  start_idx <- (i-1)*batch_size + 1
  end_idx <- min(i*batch_size, num_vertices)
  batch_vertices_sub <- vertices[start_idx:end_idx, ]
  
  # 添加顶点并获取新ID
  new_v_ids <- add.vertices(g, nv = nrow(batch_vertices_sub), attr = as.list(batch_vertices_sub))
  
  # 生成批量内部边对
  if (length(new_v_ids) >=2) {
    edge_pairs <- combn(new_v_ids, 2)
    coords_batch <- batch_vertices_sub[, c("Lon", "Lat")]
    dist_matrix <- calculate_distances(coords_batch, coords_batch)
    edge_weights <- dist_matrix[lower.tri(dist_matrix)]
    
    g <- add_edges(g, edges = as.vector(edge_pairs), attr = list(weight = edge_weights))
  }
}

关键修正点 #

• 正确处理distVincentyEllipsoid返回的矩阵，提取对应边对的距离值，避免盲目重复向量。
• 明确边的创建逻辑，按需选择全连接或批量内部连接。
• 优化顶点索引管理，确保边对ID与图中顶点ID对应。
• 通过批量处理避免一次性生成全量距离矩阵，减少内存占用。

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