2017-11-22 16 views
0

ラスタを使用していますが、7nレイヤのラスタスタックがあります。私は下の数式を使ってピクセルワイズ回帰を計算したいと思います。私はraster::calcでそれをやろうとしましたが、私の機能は、メッセージで失敗しました:できませんラスタスタックでRのピクセル単位の回帰を計算できません

「エラーlm.fitで(X、Y、オフセット=オフセット、singular.ok = singular.ok、 ...) :0(非NA)の場合。

しかし、すべてのラスタOKです、と数字(NASだけではなく)を含有し、私はそれをプロットすることができ、 と私は式

cr.sig=lm (raster::as.array(MK_trend.EVI.sig_Only) ~ raster::as.array(stack.pet)+raster::as.array(stack.tmp)+raster::as.array(stack.vap)+raster::as.array(stack.pre)+raster::as.array(stack.wet)+raster::as.array(stack.dtr)) 

で一般的な線形回帰を計算することができます。しかし、私は

で層を積み重ねるとき
allData = stack(MK_trend.EVI.sig_Only,stack.dtr,stack.wet,stack.pre,stack.vap,stack.tmp,stack.pet) 

# Regression Function, R2 
lmFun=function(x){ 
    x1=as.vector(x); 
    if (is.na(x1[1])){ 
     NA 
    } else { 
     m = lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7]) 
     return(summary(m)$r.squared) 
    } 
} 
カルク機能を試してみてください

エラーメッセージが表示されます。
私はRでかなり新しいです、そして、そうかもしれない、多分、いくつかの愚かな間違いがありますか? 処理を行うためには、何かヒントをいただきたいと思います。

+0

エラーを再現するサンプルデータを作成するコードを提供できますか? – RobertH

+0

@RobertH Robertさん、私のデータセットに問題があると思われるので、ZIPアーカイブ(1,7Gb)の形式でアップロードしてください。http://dropmefiles.com/PDZL0 –

答えて

0

ピクセルワイズ(ローカル)回帰にはcalcを使用できますが、何か他のもの(グローバルモデル)がほしいと思っているようです。

回帰がピクセル単位の場合は、各セルのx値とy値が同じで、calcを使用できます。例については、?calcを参照してください。

代わりに、各セルに1つのy(独立)変数と6 x(依存)変数値があります。これは、グローバルモデルが必要であることを示唆しています。そのためには、次のようなことができます:

library(raster) 
# example data 
r <- raster(nrow=10, ncol=10) 
set.seed(0) 
s <- stack(lapply(1:7, function(i) setValues(r, rnorm(ncell(r), i, 3)))) 
x <- values(s) 

# model 
m <- lm(layer.1 ~ ., data=data.frame(x)) 

# prediction 
p <- predict(s, m) 

これはすべての値をメモリにロードする必要があります。それができない場合は、大きな正規のサンプルを取ることができます。 sampleRegular

を参照してください。そして、あなたのアプローチが動作しない理由を示すこと:

testFun=function(x1){ 
    lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7]) 
} 

# first cell 
v <- s[1] 
v 
#  layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 layer.7 
#[1,] 4.788863 4.345578 -0.137153 3.626695 3.829971 4.120895 1.936597 

m <- testFun(v) 
m 
#Call: 
#lm(formula = x1[1] ~ x1[2] + x1[3] + x1[4] + x1[5] + x1[6] + x1[7]) 

#Coefficients: 
#(Intercept)  x1[2]  x1[3]  x1[4]  x1[5]  x1[6]  x1[7] 
#  4.789   NA   NA   NA   NA   NA   NA 


summary(m)$r.squared 
# 0 

私はあなたが報告したエラーメッセージ(ただし、すべてのR^2つの値がゼロである)を得ることはありませんが。

+0

親愛なるRobert、ありがとう解明、私は傍受(return(summary(m)$ coefficients [1])を出そうとしましたが、前と同じエラーがあります –

関連する問題