2

我有两组栅格，它们都具有相同的 x、y、z 范围。我做了两个堆栈：stacka 和 stackb。我想计算沿时间线的两个堆栈之间的每个网格单元中的 Pearson 相关系数 ( PCC )。我做了一个更简单的例子（请原谅我创建栅格的愚蠢方式）

a1<-c(1,1,1,1,1,1,1,1,NA)
a2<-c(2,2,2,2,1,2,2,NA,2)
a3<-c(3,3,3,3,3,2,NA,3,3)
b1<-c(2,2,2,2,2,2,2,2,2)
b2<-c(3,3,3,3,3,3,3,3,3)
b3<-c(4,4,4,4,4,4,4,4,4)
matrixa1<-matrix(a1,3,3)
matrixa2<-matrix(a2,3,3)
matrixa3<-matrix(a3,3,3)
matrixb1<-matrix(b1,3,3)
matrixb2<-matrix(b2,3,3)
matrixb3<-matrix(b3,3,3)
rastera1<-raster(matrixa1)
rastera2<-raster(matrixa2)
rastera3<-raster(matrixa3)
rasterb1<-raster(matrixb1)
rasterb2<-raster(matrixb2)
rasterb3<-raster(matrixb3)
stacka<-stack(rastera1,rastera2,rastera3)
stackb<-stack(rasterb1,rasterb2,rasterb3)

a_bar<-calc(stacka,mean,na.rm=TRUE)
b_bar<-calc(stackb,mean,na.rm=TRUE)
numerator<-setValues(rastera1,0)
denominator1<-numerator
denominator2<-numerator
for(i in 1:noflayers){
  numerator<-numerator+(stacka[[i]]-a_bar)*(stackb[[i]]-b_bar)
  denominator1<-denominator1+(stacka[[i]]-a_bar)^2
  denominator2<-denominator2+(stackb[[i]]-b_bar)^2
}
pearsoncoeff<-numerator/sqrt(denominator1*denominator2)

最后，我有一个栅格，每个网格单元都填充了 PCC。问题是，数据 a 是断断续续的，某些网格在某些层中是 NA。所以最终产品有一些空白。我的算法在遇到 NA 时会吐出“NA”。我需要像na.rm=TRUE计算中的一些选项，因此输出将计算任何月份的值。

我能想到的方法是is.na(stacka[[nlayers]][nrows,ncols]==FALSE在stackb中使用并找到对应的对，但这是基于单元格的，这需要大量的计算机时间。

4

2 回答 2

4

我编辑了 Paulo 推荐的在计算中处理 NA 的方法，它似乎在一系列测试中运行得很快，包括上面的数据集：

stack.correlation <- function(stack1, stack2, cor.method){
  # output template
  cor.map <- raster(stack1)
  # combine stacks
  T12 <- stack(stack1,stack2)
  rnlayers=nlayers(T12)
  # the function takes a vector, partitions it in half, then correlates
  # the two sections, returning the correlation coefficient. 
  stack.sequence.cor <- function(myvec,na.rm=T){
    myvecT1<-myvec[1:(length(myvec)/2)]
    myvecT2<-myvec[(length(myvec)/2+1):length(myvec)]
    return(cor(myvecT1,myvecT2, method =  cor.method, use="complete.obs"))
  }
  # apply the function above to each cell and write the correlation
  # coefficient to the output template. 
  cor.map <- stackApply(T12, indices = rep(1, rnlayers), 
                        fun = stack.sequence.cor, na.rm = FALSE)

  return(cor.map)
}
cor_r=stack.correlation(stacka, stackb, "pearson")
于 2013-05-23T02:08:15.290 回答
3

一种更简单的方法：

library(raster)
a1 <- raster(matrix(c(1,1,1,1,1,1,1,1,NA),3,3))
a2 <- raster(matrix(c(2,2,2,2,1,2,2,NA,2), 3, 3))
a3 <- raster(matrix(c(3,3,3,3,3,2,NA,3,3), 3, 3))
b1 <- raster(matrix(c(2,2,2,2,2,2,2,2,2), 3, 3))
b2 <- raster(matrix(c(3,3,3,3,3,3,3,3,3), 3, 3))
b3 <- raster(matrix(c(4,4,4,4,4,4,4,4,4), 3, 3))
sa <- stack(a1, a2, a3)
sb <- stack(b1, b2, b3)


funcal <- function(xy) {
    xy <- na.omit(matrix(xy, ncol=2))
    if (ncol(xy) < 2) {
        NA
    } else {
        cor(xy[, 1], xy[, 2])
    }
}

s <- stack(sa, sb)
calc(s, funcal)
于 2013-05-25T19:05:26.413 回答