如何基于R语言terra包实现spatRaster的加权核密度估计？

问题背景 #

Terra包内置的density()函数可计算SpatRaster的核密度估计，但不支持加权核密度估计。比如在估算特定土地利用类型的年均温KDE时，该类型在单元格中的占比不同，必须通过权重体现这种差异，此时内置函数无法满足需求。

现有一个自定义函数dwr()实现了加权KDE，但仅通过采样优化大栅格场景的性能，效率仍有提升空间，原函数及示例代码如下：

原自定义函数与基础示例 #

library(terra)
# Terra内置density函数基础用法
x <- rast(matrix(rnorm(100,10,1), nr = 10, nc = 10))
density(x, plot = F)

# 自定义加权栅格密度函数
dwr <- function(x, w, samp = 1, ...){
  compareGeom(x, w, lyrs = TRUE) # 检查数据与权重栅格的一致性
  w <- subst(w, NA, 0) # 将缺失权重转为0
  if(as.numeric(global(w, min)) < 0) return("Warning: Negative weights") # 检查权重是否为负
  n <- global(x, "notNA")[1,] # 统计有效数据量
  sz <- if(samp <= 1) round(samp * n) else samp # 确定采样数量
  if(sz > n) return("Warning: Sample size is too large") # 检查采样量是否超标
  ss <- spatSample(c(x,w), sz, na.rm = TRUE) # 抽取样本
  ss$wt <- ss[,2]/sum(ss[,2]) # 将权重转为占比
  density(ss[,1], weights = ss$wt, ...) # 计算加权核密度
}

# 函数调用示例
x <- rast(matrix(rnorm(100,10,1), nr = 10, nc = 10))
w <- rast(matrix(rpois(100,5), nr = 10, nc = 10))
dwr(x,w)

高效改进方案 #

针对大栅格场景，可从以下方向优化性能与精度：

1. 跳过采样，直接批量提取有效数据 #

放弃采样步骤，直接提取所有有效栅格的数值与权重，减少采样带来的随机误差和额外开销：

dwr_opt1 <- function(x, w, ...){
  # 几何一致性检查
  if(!compareGeom(x, w, lyrs = TRUE, stopOnError = FALSE)){
    stop("x和w栅格的几何属性不匹配")
  }
  # 处理缺失值与负权重
  w <- subst(w, NA, 0)
  if(global(w, min, na.rm = TRUE)[[1]] < 0){
    stop("权重不能为负值")
  }
  # 批量提取有效数据
  vals <- cbind(values(x), values(w))
  vals <- vals[complete.cases(vals[,1]), ] # 仅保留x非NA的行
  if(nrow(vals) == 0){
    stop("无有效数据可计算")
  }
  # 归一化权重
  vals[,2] <- vals[,2]/sum(vals[,2])
  # 计算加权KDE
  density(vals[,1], weights = vals[,2], ...)
}

优势：保留所有有效数据的权重信息，结果更准确；批量提取数据的IO效率远高于采样，适合中大型栅格。

2. 分块处理超大栅格 #

当栅格数据量超出内存上限时，采用分块读取+加权聚合的方式，避免内存溢出：

dwr_opt2 <- function(x, w, blocksize = 1000, ...){
  # 几何一致性检查
  if(!compareGeom(x, w, lyrs = TRUE, stopOnError = FALSE)){
    stop("x和w栅格的几何属性不匹配")
  }
  w <- subst(w, NA, 0)
  if(global(w, min, na.rm = TRUE)[[1]] < 0){
    stop("权重不能为负值")
  }
  
  # 初始化权重总和与数据容器
  total_w <- 0
  data_list <- list()
  
  # 分块读取栅格
  blocks <- makeBlocks(x, n = blocksize)
  for(i in 1:nrow(blocks)){
    b <- blocks[i,]
    x_block <- readValues(x, row = b$row, nrows = b$nrows, col = b$col, ncols = b$ncols)
    w_block <- readValues(w, row = b$row, nrows = b$nrows, col = b$col, ncols = b$ncols)
    # 过滤x为NA的记录
    valid_idx <- !is.na(x_block)
    x_valid <- x_block[valid_idx]
    w_valid <- w_block[valid_idx]
    if(length(x_valid) == 0) next
    # 累加权重总和，暂存数据
    total_w <- total_w + sum(w_valid)
    data_list[[i]] <- data.frame(x = x_valid, w = w_valid)
  }
  
  if(total_w == 0){
    stop("有效权重总和为0，无法计算")
  }
  
  # 合并数据并归一化权重
  all_data <- do.call(rbind, data_list)
  all_data$w <- all_data$w / total_w
  # 计算加权KDE
  density(all_data$x, weights = all_data$w, ...)
}

优势：通过分块控制内存占用，可处理TB级超大栅格；分块读取的效率远高于全量加载，适合内存有限的环境。

3. 预计算全局权重总和，减少重复计算 #

原函数在采样后才归一化权重，可提前计算全局权重总和，避免重复求和：

dwr_opt3 <- function(x, w, samp = 1, ...){
  # 几何一致性检查
  if(!compareGeom(x, w, lyrs = TRUE, stopOnError = FALSE)){
    stop("x和w栅格的几何属性不匹配")
  }
  w <- subst(w, NA, 0)
  min_w <- global(w, min, na.rm = TRUE)[[1]]
  if(min_w < 0){
    stop("权重不能为负值")
  }
  total_w <- global(w, sum, na.rm = TRUE)[[1]]
  if(total_w == 0){
    stop("有效权重总和为0，无法计算")
  }
  
  n <- global(x, "notNA")[1,]
  sz <- if(samp <= 1) round(samp * n) else samp
  if(sz > n){
    stop("采样量超过有效数据量")
  }
  
  # 采样时直接使用预计算的全局权重总和归一化
  ss <- spatSample(c(x,w), sz, na.rm = TRUE)
  ss$wt <- ss[,2]/total_w
  density(ss[,1], weights = ss$wt, ...)
}

优势：减少重复计算量，同时保留采样功能，适合不需要全量精度的快速估算场景。

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