基于栅格直接计算连续型数据Jaccard相似性指数的高效实现方法问询

嘿，我刚好处理过类似的连续型栅格Jaccard计算需求，你的随机采样平均方法确实在数据量大的时候会慢到让人崩溃，这里有几个直接基于栅格操作的高效方案，完全不用采样，也不需要二值化，应该能解决你的问题：

首先明确：连续型数据的Jaccard相似性通常采用逐像素取min求和除以逐像素取max求和的定义，也就是：
$$J = \frac{\sum \min(x_i, y_i)}{\sum \max(x_i, y_i)}$$
这个定义完全适配连续数据，不需要转二值图，而且可以用栅格包的向量化操作直接计算，速度远快于采样。

用terra包（现在是R中栅格处理的主流工具，比旧的raster包高效很多）的代码示例：

library(terra)

# 加载你的两个连续型栅格（确保投影、分辨率、范围一致）
r1 <- rast("your_first_raster.tif")
r2 <- rast("your_second_raster.tif")

# 如果栅格不对齐，先投影匹配
r2 <- project(r2, r1)

# 逐像素计算min和max
min_raster <- min(r1, r2)
max_raster <- max(r1, r2)

# 对所有非NA像素求和（自动忽略NA）
sum_min <- global(min_raster, sum, na.rm = TRUE)[[1]]
sum_max <- global(max_raster, sum, na.rm = TRUE)[[1]]

# 得到最终的连续型Jaccard指数
jaccard_index <- sum_min / sum_max

这个方法全程都是底层优化的栅格操作，不会循环每个像素，处理大栅格也毫无压力。

如果你需要计算多个栅格之间的两两Jaccard指数，可以把所有栅格的非NA像素提取为矩阵，然后自定义距离函数计算。注意如果栅格太大，直接提取全矩阵可能占内存，这时候可以结合分块处理：

library(terra)
library(vegan)

# 加载多个栅格为一个栅格栈
r_stack <- rast(c("raster1.tif", "raster2.tif", "raster3.tif"))

# 提取所有非NA像素为矩阵（每列对应一个栅格，每行对应一个像素）
pixel_matrix <- as.matrix(r_stack, na.rm = TRUE)

# 自定义连续型Jaccard距离函数（距离=1-相似性，适配vegan的dist格式）
continuous_jaccard_dist <- function(x) {
  n_grids <- ncol(x)
  dist_mat <- matrix(0, n_grids, n_grids)
  
  for (i in 1:(n_grids-1)) {
    for (j in (i+1):n_grids) {
      sum_min <- sum(pmin(x[,i], x[,j]), na.rm = TRUE)
      sum_max <- sum(pmax(x[,i], x[,j]), na.rm = TRUE)
      dist_mat[i,j] <- 1 - (sum_min / sum_max)
      dist_mat[j,i] <- dist_mat[i,j]
    }
  }
  as.dist(dist_mat)
}

# 计算所有栅格对的Jaccard距离矩阵
jaccard_dist_matrix <- continuous_jaccard_dist(pixel_matrix)

如果栅格超大，提取矩阵内存不够，就用上面的分块累加思路，先计算每个栅格对的sum_min和sum_max，再构建距离矩阵。

如果你的栅格是几十GB级别的，完全加载到内存不现实，可以用分块逐段计算sum_min和sum_max，最后累加得到结果：

library(terra)

r1 <- rast("extremely_large_raster1.tif")
r2 <- rast("extremely_large_raster2.tif")
r2 <- project(r2, r1)

# 设置分块大小（根据你的内存调整，比如1000x1000像素）
block_dim <- c(1000, 1000)

sum_min_total <- 0
sum_max_total <- 0

# 遍历栅格的所有块
for (block in writeStart(r1, filename = "temp.tif", overwrite = TRUE)) {
  # 读取当前块的两个栅格数据
  block_r1 <- readBlock(r1, block)
  block_r2 <- readBlock(r2, block)
  
  # 计算当前块的min和max求和
  block_sum_min <- sum(pmin(block_r1, block_r2, na.rm = TRUE))
  block_sum_max <- sum(pmax(block_r1, block_r2, na.rm = TRUE))
  
  # 累加到总结果
  sum_min_total <- sum_min_total + block_sum_min
  sum_max_total <- sum_max_total + block_sum_max
  
  # 占位写入（不需要实际输出，随便写原数据即可）
  writeBlock(r1, block, block_r1)
}
writeStop(r1)

# 计算最终Jaccard指数
jaccard_index <- sum_min_total / sum_max_total

还可以结合并行计算加速分块处理，用doParallel包注册并行集群后，用terra::app批量处理每个块，效率会更高。

这些方案都完全避开了随机采样的低效问题，直接利用所有非NA像素计算，而且全程不需要转二值图，应该能完美解决你的场景。

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