使用%dopar%并行化生成路径返回NULL，疑与terra有关？

问题现象 #

用foreach的%do%串行执行create_routes1时，能正常生成所有点对的最小成本路径网络；但改用%dopar%并行执行create_routes2时，返回结果为NULL。代码中未直接调用terra，但依赖leastcostpath包生成的成本表面。

核心原因 #

1. terra对象并行序列化失效：slope_cs是基于terra栅格生成的成本表面对象，默认的进程间对象传递机制无法正确解析terra栅格数据，导致子进程中create_lcp无法获取有效输入，返回空结果。
2. 依赖包未完全加载：create_lcp内部依赖terra和sf的函数，但foreach的.packages仅指定了leastcostpath，子进程缺少必要库导致执行失败；加上.errorhandling = "remove"，所有失败的迭代被移除后最终返回NULL。

解决方案 #

1. 在foreach的.packages中添加terra和sf，确保子进程加载全部依赖库；
2. 用terra::wrap()将成本表面序列化为可安全传递的二进制格式，子进程中用terra::unwrap()还原；
3. 集群创建时使用PSOCK类型（跨平台兼容），保证对象传递机制正常。

修正后的代码 #

library(devtools)
install_github("josephlewis/leastcostpath@dev")
library(leastcostpath)
library(foreach)
library(doParallel)
library(terra)
library(sf)

r <- terra::rast(system.file("extdata/SICILY_1000m.tif", package="leastcostpath"))

slope_cs <- create_slope_cs(x = r, cost_function = "tobler", neighbours = 4)

locs <- sf::st_sf(geometry = sf::st_sfc(
 sf::st_point(c(839769, 4199443)),
 sf::st_point(c(1038608, 4100024)),
 sf::st_point(c(907695, 4145478)),
 sf::st_point(c(1054446, 4232288)),
 sf::st_point(c(957989, 4208863)),
 crs = terra::crs(r)))

# 原串行版本（保留）
create_routes1 <- function (x, locations, cost_distance = FALSE) 
{
 network <- expand.grid(1:nrow(locations), 1:nrow(locations))
 network <- network[network[, 1] != network[, 2], ]
 lcp_network <- foreach::foreach(i = 1:nrow(network), .packages = c("leastcostpath", "terra", "sf"), 
                                 .errorhandling = "remove", .combine = "rbind") %do% 
   {
     lcp <- create_lcp(x = x, origin = locations[network[i, 1], ], destination = locations[network[i, 2],], cost_distance = cost_distance)
     lcp$origin_ID <- network[i, 1]
     lcp$destination_ID <- network[i, 2]
     return(lcp)
   }
 
 return(lcp_network)
}

lcps1 <- create_routes1(x = slope_cs, locations = locs, cost_distance = FALSE)

# 修正后的并行版本
create_routes2 <- function (x, locations, cost_distance = FALSE, ncores = 2)
{
 # 序列化成本表面对象
 wrapped_x <- terra::wrap(x)
 
 myCluster <- makeCluster(ncores, type = "PSOCK")
 registerDoParallel(myCluster)
 
 network <- expand.grid(1:nrow(locations), 1:nrow(locations))
 network <- network[network[, 1] != network[, 2], ]
 lcp_network <- foreach::foreach(i = 1:nrow(network), 
                                 .packages = c("leastcostpath", "terra", "sf"), 
                                 .errorhandling = "remove", 
                                 .combine = "rbind") %dopar% 
   {
     # 子进程中反序列化对象
     unwrapped_x <- terra::unwrap(wrapped_x)
     lcp <- create_lcp(x = unwrapped_x, origin = locations[network[i, 1], ], destination = locations[network[i, 2],], cost_distance = cost_distance)
     lcp$origin_ID <- network[i, 1]
     lcp$destination_ID <- network[i, 2]
     return(lcp)
   }
 
 stopCluster(myCluster)
 
 return(lcp_network)
}

lcps2 <- create_routes2(x = slope_cs, locations = locs, cost_distance = FALSE, ncores = 2)

验证说明 #

运行修正后的create_routes2后，lcps2会和lcps1结构一致，包含所有有效点对的最小成本路径及对应ID。若仍有结果缺失，可暂时移除.errorhandling = "remove"，查看具体错误信息排查个别点对的路径生成问题。

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