R中两个光栅堆栈之间每个位置的最佳拟合线

编程问答 2022-05-17

问题描述

我知道如何通过以下方式计算两个栅格堆栈之间每个位置的 r^2 相关性：

library(raster)
s1<-rasterstack1
s2<-rasterstack2
a <- length(s1)
b <- length(s2)
s <- stack(s1,s2)
fun=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[1:a] ~ x[(a+1):(a+b)]);summary(m)$r.squared }}
r.squared <- calc(s,fun)

但是，我试图找到最适合每个位置的线。因此，假设输出栅格/矩阵类似于：

[y=x^.8,y=X^.3]
[y=x^.5,y=x^.9]

其中 y 是给定 x 变量的估计因变量（r^2 值基于什么）。然后理论上我希望能够给出一个输入 x 并有一个估计的栅格或估计 y 值的矩阵...... 我假设输出将是这样的：

output<- function(x) A <- matrix(c(x**.8,x**.3,x**.5,x**.9),2,2)
end<-output(x)

我挣扎的部分是计算这些方程并将它们存储为矩阵形式。我得到的最远的是：

predict <-function(rasterstack){
  # Values for (layers,ncell,ncol,nrow,method,crs,extent) come straight from the input raster stack
  # e.g. nlayers(rasterstack),ncell(rasterstack)... etc.
  print(paste("Start:",Sys.time()))
  print("Loading parameters")
  layers=nlayers(rasterstack);ncell=ncell(rasterstack);
  ncol=ncol(rasterstack);nrow=nrow(rasterstack);crs=crs(rasterstack);
  extent=extent(rasterstack);pb = txtProgressBar(min = 0,max = ncell,initial = 0)
  print("Done loading parameters")
  mtrx <- as.matrix(rasterstack,ncol=layers)
  empt <- matrix(nrow=ncell,ncol=2)
  print("Initiating loop operation")
  numb=1
  
  for (i in 1:ncell){
    setTxtProgressBar(pb,i)
    if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){ 
      empt[i,1] <- NA 
    } else 
      if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
        empt[i,1] <- NA 
      } else {
        
        empt[i,1]<-as.expression(coef(lm( mtrx[i,((layers/2)+1):layers] ~ mtrx[i,1:(layers/2)]))[2])*x+(coef(lm( mtrx[i,1:(layers/2)]))[1])
        
        
        
  }
  print("Creating empty raster")
  eqmatrix <- raster(nrows=nrow,ncols=ncol,crs=crs)
  extent(eqmatrix) <- extent
  print("Populating correlation raster")
  values(eqmatrix) <- empt[,1]
  print(paste("End on",Sys.time()))
  eqmatrix
  
}
}

predict_matrix<-predict(rasterstack=s)

但这会产生错误： "as.expression(coef(lm(mtrx[i,((layers/2) + 1):layers] ~mtrx[i,: 二元运算符的非数字参数

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

rrr-rasterrasterstatistics