2016-06-21 13 views
0

Ich versuche, einen pixelweisen granger kausalen Test auf zwei Rasterstapel mit je 60 Raster durchzuführen. Das folgende Beispiel hat nur 20 Rastern:Granger kausale Tests auf Raster-Stacks in R

library(raster) 
library(lmtest) 

r <- raster(ncol=10, nrow=10) 
r[]=1:ncell(r) 
S <- stack(r,r,r,r,r,r,r,r,r,r,r,r) 
R <- stack(r,r,r,r,r,r,r,r,r,r,r,r) 
FNO2<-stack(S,R) 

Die ursprüngliche Funktion th „lmtest“ Paket mit ist:

D<- grangertest(degg ~ dchick, order=4) 

Hier ist eine Modifikation Ich habe die ursprüngliche grangertest Funktion auf Raster Stapel laufen zu lassen?

funG <- function(x) { if (is.na(x[1])){ NA } else { grangertest(x[13:24] ~ x[1:12],order=1)}} 
granger<-calc(FNO2,funG) 

Wo FNO2 ist der Stapel von beiden Raster-Stacks. Ich erhalte den folgenden Fehler:

Error in `colnames<-`(`*tmp*`, value = c("x", "y", "x_1", "y_1")) : 
length of 'dimnames' [2] not equal to array extent 

Wie ändere ich diese Funktion für Raster bitte?

Antwort

1

Sie müssen prüfen, was grangertest kehrt:

library(lmtest) 
x <- grangertest(egg ~ chicken, order = 3, data = ChickEgg) 

class(x) 
#[1] "anova"  "data.frame" 

x 
#Granger causality test 
# 
#Model 1: egg ~ Lags(egg, 1:3) + Lags(chicken, 1:3) 
#Model 2: egg ~ Lags(egg, 1:3) 
# Res.Df Df  F Pr(>F) 
#1  44     
#2  47 -3 0.5916 0.6238 

Das ist nicht etwas, das wir in einer Rasterebene

str(x) 
#Classes ‘anova’ and 'data.frame':  2 obs. of 4 variables: 
# $ Res.Df: num 44 47 
# $ Df : num NA -3 
# $ F  : num NA 0.592 
# $ Pr(>F): num NA 0.624 
# - attr(*, "heading")= chr "Granger causality test\n" "Model 1: egg ~ Lags(egg, 1:3) + Lags(chicken, 1:3)\nModel 2: egg ~ Lags(egg, 1:3)" 
> 

Ich bin nicht sicher halten können, was Sie nach, aber wenn es das war p-Wert, vielleicht

x$'Pr(>F)'[2] 
#[1] 0.6237862 

Dann können wir die Funktion zu etwas li ändern ke:

funG <- function(x) { if (is.na(x[1])){ 
       NA 
      } else { 
       grangertest(x[13:24] ~ x[1:12],order=1)$'Pr(>F)'[2] 
      } 
    } 

Beispiel RasterStack:

library(raster) 
r <- raster(ncol=10, nrow=10) 
set.seed(9) 
FNO2 <- stack(sapply(1:24, function(i) setValues(r, runif(ncell(r))))) 

Testen Sie die Funktion:

d <- as.vector(FNO2[1]) 
funG(d) 
#[1] 0.03248222 

Jetzt:

granger<-calc(FNO2,funG) 

granger 
#class  : RasterLayer 
#dimensions : 10, 10, 100 (nrow, ncol, ncell) 
#resolution : 36, 18 (x, y) 
#extent  : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) 
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : in memory 
#names  : layer 
#values  : 0.007425836, 0.9895796 (min, max) 
+0

Dank @RobertH. das hat funktioniert. –