R 中的成对栅格比较：for 循环的替代方案？

Question

如何有效地比较分布栅格对（raster 层仅包含 0 和 1）？我需要衡量 ~6500 个单独的全球栅格之间的相似性。来自 SDMTools 的 Istat 应该可以完成这项工作。

这是我的代码：

library(raster)
library(SDMTools)

创建可重现的示例数据：值为 0 和 1 的栅格

# first raster
r1 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=18000000, ymn=-9000000, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=0)
r2 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=0, ymn=0, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=2)
r12 <- mosaic(r1, r2, fun=mean)

# second raster
r3 <- raster(nrow=1800, ncol=3600, xmn=-18000000, xmx=18000000, ymn=-9000000, ymx=9000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=0)
r4 <- raster(nrow=1800, ncol=3600, xmn=-12000000, xmx=15000000, ymn=2000000, ymx=3000000, 
             crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs", 
             resolution=10000, vals=2)
r34 <- mosaic(r3, r4, fun=mean)

列出栅格

files_list <- list(r12, r34)

创建空矩阵以填充来自循环的数据

ras_comp <- matrix(NA, nrow=length(files_list), ncol=length(files_list))
ras_comp
# label rows and columns of matrix
rownames(ras_comp) <- c("r12", "r34")
colnames(ras_comp) <- c("r12", "r34")
ras_comp

循环比较所有可能的 matrices/rasters

for (i in 1:length(files_list)) {
  # load raster i
  ras_i <- as.matrix(files_list[[i]])

  for (j in 1:length(files_list)) {
    # load raster j
    ras_j <- as.matrix(files_list[[j]])

    # compare both rasters
    ras_Istat <- Istat(ras_i, ras_j, old=F)

    # write value into matrix
    ras_comp[i,j] <- ras_Istat
  }
}

检查最终矩阵

ras_comp
> ras_comp
          r12       r34
r12 1.0000000 0.1814437
r34 0.1814437 1.0000000

使用 as.matrix 将栅格转换为矩阵显着减少了计算时间，最终得到的 table 是我需要的，但是对数千个栅格执行此操作需要永远完成。如何优化代码以便更有效地比较栅格？

Answer 1

Istat 在进行简单计算之前进行了一系列测试和缩放。如果您知道这些测试通过了，您可以一次性进行缩放并处理缩放后的值。它确实：

if (length(which(dim(x) == dim(y))) != 2) 
    stop("matrix / raster objects must be of the same extent")
if (min(c(x, y), na.rm = T) < 0) 
    stop("all values must be positive")

然后检查两个栅格的位置 "finite"，其中包括 NA 个值：

pos = which(is.finite(x) & is.finite(y))

然后计算栅格的缩放值：

px = x[pos]/sum(x[pos])
py = y[pos]/sum(y[pos])
H = sqrt(sum((sqrt(px) - sqrt(py))^2))

如果old=FALSE像你那样returns:

    return(1 - (H^2)/2)

> Istat(r12,r34)
[1] 0.1814437

如果我删除测试并编写一个适用于缩放值的函数，我可以将其归结为：

fIstat = function(px,py){
    1 - (sum((sqrt(px) - sqrt(py))^2))/2
}

通过缩放栅格进行测试运行:

r12px = r12[]/sum(r12[])
r34px = r34[]/sum(r34[])
fIstat(r12px, r34px)
# [1] 0.1814437

相同的值。很好，但是它更快吗？

> microbenchmark(fIstat(r12px, r34px), Istat(r12,r34))
Unit: milliseconds
                 expr        min         lq       mean     median         uq
 fIstat(r12px, r34px)   49.95867   78.28649   78.10863   79.45235   80.85234
      Istat(r12, r34) 1084.84825 1181.31116 1217.64122 1212.93180 1263.50811
       max neval
  106.6803   100
 1349.0239   100

是的，很大一部分。

所以... 如果您的数据没有缺失值或无穷大，请创建这些缩放栅格值的 files_list，调用我的 fIstat，只在上面的三角形上循环，你应该把它加快 10 倍。

R 中的成对栅格比较：for 循环的替代方案？

Pairwise raster comparison in R: alternative to for-loop?

for-loop

r

r-raster

pairwise-distance

sdmtools