Блог

Я знаю, как рассчитать корреляцию 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,((layers/2) 1):layers] ~ 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, ((слои/2) 1):слои] ~mtrx[i, :
нечисловой аргумент двоичного оператора