如何有效比较分布栅格对(raster
仅包含0和1的图层)?我需要衡量约6500个单个全局栅格之间的相似性。Istat
从SDMTools
应该做的工作。
这是我的代码:
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
可以显着减少计算时间,最终的最终表正是我所需要的,但是要成千上万个栅格执行此操作将永远需要完成。为了更有效地比较栅格,如何优化代码?
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
这些缩放的栅格值中的一个,调用my fIstat
,仅在上三角处循环,您应将其加快10倍。
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句