R语言doParallel并行计算通视性时遇[readValues]读取错误求助

问题背景 #

使用R语言doParallel+foreach并行化，基于DEM的rasterprofile函数计算地形观测点间的通视性。曾成功处理1000×1000规模数据，处理时间压缩至原1/10；但处理邻区341×1200规模数据时出现异常。

异常现象 #

• 计算规模超过10×1200时，并行计数器显示完成，但无法输出结果，报错：task 1 failed - "[readValues] cannot read values"
• 缩减origin数据至10行（10×1200规模）可完成计算，但随机出现closing unused connection警告
• 使用投影坐标系EPSG:32635，DEM栅格覆盖所有点位

数据示例 #

origin数据结构：

> origin
           X       Y   h
 [1,] 453727.1 4706737 1.5
 [2,] 452999.1 4707739 2.0
 ...

核心代码 #

并行计算逻辑 #

cl <- makeCluster(10) # 16核机器，未过载
registerDoParallel(cl)

ie_table <- foreach (i = 1:nrow(origin), # 当前测试用40行
                     .combine = 'rbind',
                     .packages = c("data.table", "raster")) %do% {  
  print(i)
  foreach(e = 1:nrow(target),  # 1200行
          .combine = 'rbind',.packages = c("data.table", "raster")) %dopar% {
  print(e)
  result <-  cansee(Y35, origin[i,1:2], target[e,1:2], h1 = origin[i,3], h2 = target[e,3])
  result_row <- c(i,e, result)
  }
}

stopCluster(cl)

通视性判断函数 #

cansee <- function(r, xy1, xy2, h1=0, h2=0){
### xy1能否在DEM r上看到xy2？
### r为DEM（Y35），与x,y,z单位一致
### xy1和xy2为长度为2的x,y坐标向量
### h1和h2为额外高度偏移（即特征高度）
### （例如观测点的观察者高度等）
    xyz = rasterprofile(r, xy1, xy2)
    np = nrow(xyz)-1
    h1 = xyz$z[1] + h1
    h2 = xyz$z[np] + h2
    hpath = h1 + (0:np)*(h2-h1)/np
    return(!any(hpath < xyz$z))
} 

解决思路 #

1. 修复栅格对象的并行传递 #

并行集群启动时，DEM栅格Y35可能未被正确传递到子进程，导致读取失败。显式导出栅格到所有子进程：

cl <- makeCluster(10)
registerDoParallel(cl)
clusterExport(cl, "Y35") # 将栅格对象传递给每个子进程

也可改用raster包的clusterR函数，它针对栅格并行操作优化了数据传递逻辑，避免手动处理的连接问题。

2. 解决连接资源泄漏 #

closing unused connection警告提示栅格文件句柄未正确释放，修改cansee函数显式关闭连接（仅针对磁盘存储的栅格有效）：

cansee <- function(r, xy1, xy2, h1=0, h2=0){
    xyz = rasterprofile(r, xy1, xy2)
    np = nrow(xyz)-1
    h1 = xyz$z[1] + h1
    h2 = xyz$z[np] + h2
    hpath = h1 + (0:np)*(h2-h1)/np
    close(r) # 释放栅格文件连接
    return(!any(hpath < xyz$z))
}

若Y35是内存栅格，可跳过此操作，转而将磁盘栅格读入内存：Y35 <- readAll(Y35)，减少IO竞争。

3. 优化并行结构 #

当前嵌套foreach（外层%do%+内层%dopar%）易引发资源竞争，改为单层并行处理所有任务组合：

cl <- makeCluster(10)
registerDoParallel(cl)
clusterExport(cl, c("Y35", "origin", "target", "cansee"))

# 生成所有(i,e)任务组合
task_grid <- expand.grid(i=1:nrow(origin), e=1:nrow(target))

ie_table <- foreach(k=1:nrow(task_grid),
                    .combine='rbind',
                    .packages=c("data.table", "raster")) %dopar% {
  i <- task_grid$i[k]
  e <- task_grid$e[k]
  result <- cansee(Y35, origin[i,1:2], target[e,1:2], h1=origin[i,3], h2=target[e,3])
  c(i,e, result)
}
stopCluster(cl)

4. 栅格数据预处理 #

• 用summary(Y35)检查DEM是否存在NA或异常值，若有则用na.omit()或mask()清理
• 确认所有origin/target坐标均在栅格范围内：用extract(Y35, origin[,1:2])验证是否返回有效值

5. 定位具体错误点 #

在并行任务中添加错误捕获，定位出错的观测点对：

result_row <- tryCatch({
  result <- cansee(Y35, origin[i,1:2], target[e,1:2], h1=origin[i,3], h2=target[e,3])
  c(i,e, result)
}, error=function(err){
  c(i,e, NA, err$message)
})

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