在R中计算两组对应点绕陆地的通行距离

核心思路 #

针对gridDistance仅支持单原点的限制，我们可以通过单起点路径计算+终点距离提取的方式实现点对距离计算，并借助批量迭代工具处理500+组数据。同时使用terra包替代维护中的raster包，提升计算效率。

步骤1：准备环境与数据 #

# 安装并加载必要包
install.packages(c("terra", "sf", "rnaturalearth", "fasterize"))
library(terra)
library(sf)
library(rnaturalearth)
library(fasterize)

# 获取北美陆地矢量数据并转换为terra格式
na_land <- ne_countries(scale = "large", continent = "North America", returnclass = "sf")
na_land_vect <- vect(na_land)

# 创建栅格（调整分辨率平衡精度与速度，示例用1000x1000）
r <- rast(ext(na_land_vect), nrows = 1000, ncols = 1000, crs = crs(na_land_vect))
# 栅格化陆地：陆地设为NA（障碍），海洋设为1（可通行）
land_raster <- fasterize(na_land, r)
ocean_raster <- ifel(is.na(land_raster), 1, NA)

# 加载并整理点对数据
tagxy <- data.frame(Y_tagged = c(45.141, 45.311, 44.954),
                    X_tagged = c(-65.501, -65.221, -66.014))
recapxy <- data.frame(Y_returned = c(43.187, 43.295, 43.279),
                      X_returned = c(-65.754, -66.051, -66.348))

# 合并为点对数据集，转换为terra点对象
point_pairs <- cbind(tagxy, recapxy)
colnames(point_pairs) <- c("tag_y", "tag_x", "recap_y", "recap_x")

# 转换点为terra矢量，匹配栅格坐标系
point_pairs$tag_pt <- vect(st_as_sf(point_pairs, coords = c("tag_x", "tag_y"), crs = crs(ocean_raster)))
point_pairs$recap_pt <- vect(st_as_sf(point_pairs, coords = c("recap_x", "recap_y"), crs = crs(ocean_raster)))

步骤2：定义单对距离计算函数 #

该函数会：

1. 将数据转换为UTM投影（获得米级实际距离）
2. 计算起点到全栅格的绕陆最短路径
3. 提取终点对应的距离值

calc_land_avoid_distance <- function(tag_pt, recap_pt, ocean_rast) {
  # 转换为UTM投影（北美东部选zone19N，根据实际区域调整）
  utm_crs <- "+proj=utm +zone=19 +datum=WGS84 +units=m +no_defs"
  ocean_utm <- project(ocean_rast, utm_crs)
  tag_utm <- project(tag_pt, utm_crs)
  recap_utm <- project(recap_pt, utm_crs)
  
  # 计算绕陆路径距离（避开NA的陆地）
  dist_rast <- distance(ocean_utm, tag_utm)
  # 提取终点的距离值
  extract(dist_rast, recap_utm)[[1]]
}

步骤3：批量计算所有点对 #

使用mapply实现批量迭代，也可根据需求改用并行计算（parallel::mcmapply）提升速度

# 批量计算所有点对的绕陆距离（单位：米）
point_pairs$distance_m <- mapply(calc_land_avoid_distance,
                                 point_pairs$tag_pt,
                                 point_pairs$recap_pt,
                                 MoreArgs = list(ocean_rast = ocean_raster))

# 查看结果
print(point_pairs[, c("tag_x", "tag_y", "recap_x", "recap_y", "distance_m")])

关键优化建议 #

• 分辨率调整：若1000x1000精度不足，可逐步提升至2000x2000；500+点对不建议用10000x10000（计算时间过长）
• 并行处理：对于大量点对，可使用parallel包开启多线程：

  library(parallel)
  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, c("calc_land_avoid_distance", "ocean_raster", "utm_crs"))
  point_pairs$distance_m <- parMapply(cl, calc_land_avoid_distance,
                                      point_pairs$tag_pt,
                                      point_pairs$recap_pt,
                                      MoreArgs = list(ocean_rast = ocean_raster))
  stopCluster(cl)
  

• CRS匹配：确保所有数据（陆地、栅格、点）的坐标系完全一致，避免投影错误

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