如何解决R中两个光栅堆栈之间每个位置的最佳拟合线
我知道如何通过以下方式计算两个栅格堆栈之间每个位置的 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 举报,一经查实,本站将立刻删除。