如何高效提取每个点位置的最近非NA栅格值？

针对你遇到的两个性能瓶颈，这里提供几个层级的优化方案：

1. 替换全栅格距离计算，用邻域逐步搜索 #

原函数计算全栅格的距离矩阵，这在栅格较大时完全没必要。我们可以从点所在像元开始，逐步扩大8邻域范围，直到找到最近的非NA像元，大幅减少计算量：

# 优化后的找最近非NA值函数
get_nearest_non_na_fast <- function(point, raster_layer) {
  # 将sf点转为terra的SpatVector
  p_vect <- vect(point)
  # 获取点对应的栅格像元索引
  cell_idx <- cellFromXY(raster_layer, p_vect)
  
  # 如果当前像元非NA，直接返回值
  if (!is.na(raster_layer[cell_idx])) {
    return(raster_layer[cell_idx])
  }
  
  # 从1开始逐步扩大搜索邻域的步长
  step <- 1
  while(TRUE) {
    # 获取当前步长下的所有邻域像元
    adj_cells <- adjacent(raster_layer, cell_idx, directions = 8, pairs = FALSE, include = FALSE, distance = step)
    # 筛选邻域中的非NA像元
    valid_cells <- adj_cells[!is.na(raster_layer[adj_cells])]
    
    if (length(valid_cells) > 0) {
      # 计算有效像元到目标点的距离，取最近的一个
      cell_coords <- xyFromCell(raster_layer, valid_cells)
      p_coords <- xyFromCell(raster_layer, cell_idx)
      distances <- sqrt((cell_coords[,1] - p_coords[1])^2 + (cell_coords[,2] - p_coords[2])^2)
      nearest_cell <- valid_cells[which.min(distances)]
      return(raster_layer[nearest_cell])
    }
    # 没找到就扩大步长继续搜
    step <- step + 1
  }
}

2. 批量处理栅格栈，取消逐图层循环 #

利用terra::extract的批量提取能力，一次性获取所有图层的点值，再统一处理NA，避免循环带来的开销：

# 一次性提取所有栅格图层的点值
extracted_all <- terra::extract(r, pnts, method = "simple")
# 拆分点ID和变量值列
pnts_id <- extracted_all[,1]
values_all <- extracted_all[,-1]

# 批量处理所有变量的NA值
corrected_all <- apply(values_all, 2, function(col) {
  var_name <- names(col)
  # 找到当前变量中NA的位置
  na_idx <- which(is.na(col))
  if (length(na_idx) == 0) {
    return(col)
  }
  # 对每个NA点调用优化后的函数
  col[na_idx] <- sapply(na_idx, function(j) {
    get_nearest_non_na_fast(pnts[j, ], r[[var_name]])
  })
  return(col)
})

# 将修正后的值合并到点数据集中
pnts <- cbind(pnts, corrected_all)
# 查看结果
print(pnts)

3. 超大规模数据进阶优化：预存非NA像元+KD树搜索 #

如果你的栅格非NA区域固定，且点数量极多，可以预存所有非NA像元的坐标和值，用KD树做快速最近邻搜索（不用每次搜索都遍历邻域）：

# 预计算单个栅格图层的非NA像元坐标和值
prep_raster_data <- function(raster_layer) {
  valid_cells <- which(!is.na(raster_layer[]))
  valid_coords <- xyFromCell(raster_layer, valid_cells)
  valid_values <- raster_layer[valid_cells]
  list(coords = valid_coords, values = valid_values)
}

# 为栅格栈的每个图层预存数据
raster_prepped <- lapply(1:nlyr(r), function(i) prep_raster_data(r[[i]]))
names(raster_prepped) <- names(r)

# 用预存数据快速找最近非NA值
get_nearest_kdtree <- function(point, prep_data) {
  p_coords <- st_coordinates(point)[1,]
  # 计算所有预存坐标到点的距离，取最近的
  distances <- sqrt((prep_data$coords[,1] - p_coords[1])^2 + (prep_data$coords[,2] - p_coords[2])^2)
  nearest_idx <- which.min(distances)
  prep_data$values[nearest_idx]
}

# 批量处理NA值
corrected_all_kdtree <- apply(values_all, 2, function(col) {
  var_name <- names(col)
  na_idx <- which(is.na(col))
  if (length(na_idx) == 0) return(col)
  col[na_idx] <- sapply(na_idx, function(j) {
    get_nearest_kdtree(pnts[j, ], raster_prepped[[var_name]])
  })
  return(col)
})

# 合并结果
pnts <- cbind(pnts, corrected_all_kdtree)

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