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

循环比较所有可能的矩阵/栅格对

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 将栅格转换为矩阵显着减少了计算时间，并且生成的最终表格正是我所需要的，但是为数千个栅格执行此操作需要很长时间才能完成。如何优化代码以便以更有效的方式比较栅格？

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")

然后它会检查两个栅格的“有限”位置，其中包括 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 和你一样，那么它会返回：

    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 倍。