点数据快速提取栅格值：如何加速raster::extract()函数运行

原并行失效原因 #

raster包的beginCluster仅对projectRaster等少数内置函数做了并行适配，extract处理点坐标时未接入该并行框架，因此只会占用单核心运行，符合你观察到的现象。

下面给出3种可落地的提速方案，优先级从上到下： #

方案1：手动计算行列号直接索引（性能最高，适配你的场景） #

你的XY坐标与DTM的CRS完全一致，无需重投影校验，可直接通过栅格参数计算对应行列号索引取值，跳过所有extract的冗余逻辑，性能比调用官方函数高10倍以上。

library(data.table)
library(raster)

# 提取DTM固定参数
xmin <- DTM@extent@xmin
ymax <- DTM@extent@ymax
res_x <- DTM@file@res[1]
res_y <- DTM@file@res[2]
ncol_dtm <- ncol(DTM)

# 读取DTM所有值到内存（如果内存不足可转用terra包分块读取）
dtm_vals <- values(DTM)

# data.table批量计算行列号，无循环耗时极短
XY[, `:=`(
  col = floor((x - xmin)/res_x) + 1,
  row = floor((ymax - y)/res_y) + 1
)]

# 过滤超出栅格范围的异常点，避免索引报错
XY <- XY[col >= 1 & col <= ncol(DTM) & row >= 1 & row <= nrow(DTM)]

# 计算栅格向量索引位置，直接取值
XY[, idx := (row - 1)*ncol_dtm + col]
XY[, altitude := dtm_vals[idx]]

• 注意：如果DTM总大小超过服务器内存，可将XY按空间范围拆分为多块，逐块读取对应DTM区域值计算即可。

方案2：换用terra包（改造成本最低） #

terra是raster包作者推出的下一代栅格处理包，底层用C++重写，点提取效率是raster的3-10倍，原生优化无需手动开并行。

library(terra)
library(data.table)

# 读取DTM为terra对象
dtm_terra <- rast("dtm_16x16_utm33.tif")

# 直接提取，兼容data.table格式的XY坐标
altitude <- extract(dtm_terra, XY)

# 合并结果到原XY表
XY <- cbind(XY, altitude)

方案3：手动分块并行（适配必须使用raster包的场景） #

如果受环境限制只能用raster包，可手动拆分XY为多块，用parallel包实现并行提取，避免单核心瓶颈。

library(raster)
library(data.table)
library(parallel)

# 配置核心数，预留2核给系统进程
use_core <- detectCores() - 2

# 拆分XY为对应核心数的分块
xy_block <- split(XY, cut(1:nrow(XY), use_core, labels = F))

# 初始化集群
cl <- makeCluster(use_core)
clusterExport(cl, "DTM")

# 并行执行提取
alt_list <- parLapply(cl, xy_block, function(b) raster::extract(DTM, b))
stopCluster(cl)

# 合并结果
XY$altitude <- unlist(alt_list)

• 注意：如果DTM过大导出到集群节点内存占用过高，可将extract逻辑改为每个节点独立读取DTM文件，避免内存溢出。

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