如何在R中优化超11000个多边形质心的地理距离计算？

核心思路：利用规则网格特性跳过冗余计算 #

你的网格是固定111km×111km的全球规则网格，完全不需要依赖地理空间包的通用函数，直接基于坐标的简化计算是最高效的路径。

1. 直接推导质心坐标，跳过st_centroid #

规则网格的质心坐标可以直接从网格的经纬度范围或行列号算出，无需调用空间转换函数：

• 若网格数据包含lon_min/lon_max、lat_min/lat_max字段，质心坐标直接取(mean(lon_min, lon_max), mean(lat_min, lat_max))
• 若有网格行列号，可直接用行列号换算成近似经纬度（比如每列对应1°经度，每行对应1°纬度，正好匹配111km的网格尺度），彻底省去空间对象的创建和转换开销

2. 用简化球面距离公式替代st_distance #

由于精度要求不高，无需st_distance的严格椭球计算，**哈弗辛公式（Haversine）**的简化版本完全满足需求，且计算速度远快于通用地理函数：

• 哈弗辛公式基于球面近似，对于111km级别的网格，误差可忽略
• 向量化实现的公式能批量处理百万级距离对

3. 向量化计算全距离矩阵，避免循环 #

如果必须计算所有质心对的距离，用R的向量化操作代替逐对循环，效率提升显著：

# 假设已生成质心经纬度数据框
centroids_df <- data.frame(
  lon = c(0.5, 1.5, 2.5, -0.5, -1.5),
  lat = c(0.5, 0.5, 0.5, -0.5, -0.5)
)

# 简化哈弗辛距离函数（返回公里数）
haversine_dist <- function(lon1, lat1, lon2, lat2) {
  # 转换为弧度
  lon1 <- lon1 * pi / 180
  lat1 <- lat1 * pi / 180
  lon2 <- lon2 * pi / 180
  lat2 <- lat2 * pi / 180
  
  dlon <- lon2 - lon1
  dlat <- lat2 - lat1
  
  a <- sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  
  # 地球半径取6371公里
  6371 * c
}

# 向量化生成全距离矩阵
n <- nrow(centroids_df)
distance_matrix <- outer(1:n, 1:n, function(i,j) {
  haversine_dist(centroids_df$lon[i], centroids_df$lat[i],
                 centroids_df$lon[j], centroids_df$lat[j])
})

4. 按需计算：只保留必要的距离对 #

如果不需要完整的距离矩阵（比如仅需相邻网格或特定范围内的距离），可进一步优化：

• 基于网格行列号筛选：只计算行列差在±N范围内的对
• 用空间索引快速筛选：通过sf::st_is_within_distance设置距离阈值，只提取符合条件的质心对，避免全矩阵计算的内存和时间浪费

5. 底层加速：用Rcpp或并行计算 #

如果向量化仍无法满足速度需求，可：

• 用Rcpp实现哈弗辛公式的底层C++版本，速度提升数倍
• 用data.table结合parallel包实现并行化计算，利用多核CPU资源

内容的提问来源于stack exchange，提问作者Rubén Bernardo-Madrid