регрессия фокального растра с RasterStack

Question

Я пытаюсь запустить многомерную регрессию с разными слоями в RasterStack с использованием focal {raster} или localFun {raster}. С помощью similar post и rasterreference manual мой код отлично работает с одним RasterLayers в качестве входного сигнала (см. Пример воспроизводимого, хотя, вероятно, «неуклюжий», пример ниже). Тем не менее, я хотел бы сделать это, используя разные уровни в RasterStack, как описано в {SECTION2} приведенного ниже кода. Я бы очень признателен за любые советы.

Спасибо

КОД:

library(raster) 

#%%%%%%%%%%%%%%%%%%%%% 
## SECTION1 
#%%%%%%%%%%%%%%%%%%%%% 

# create test data 
set.seed(0) 
resp = expl = raster(nrow=10, ncol=10) 
# response variable 
resp = setValues(resp,runif(100,min=15,max=45)) 
# explanatory variable 
expl = setValues(expl,runif(100,min=2,max=6)) 
expl = expl * resp 
resp[1:5] = NA; expl[1:5] = NA # add some NA values 
par(mfrow=c(1,2)) 
plot(resp); plot(expl) 
#.............................................................. 

# check global lm() results 
data1.df = na.omit(as.data.frame(stack(list(resp=resp,expl=expl)))) 
head(data1.df) 
data1.lm = lm(resp ~ expl, data=data1.df) 
(data1.lmSum = summary(data1.lm)) 
data1.lmSum$coefficients[1];data1.lmSum$coefficients[2];data1.lmSum$coefficients[8] 
data1.lmSum$r.squared 
data1.lmSum$sigma 
# pf(data1.lmSum$fstatistic[1], data1.lmSum$fstatistic[2], data1.lmSum$fstatistic[3],lower.tail = FALSE)  
#.............................................................. 

# lm function for focal {raster} with RasterLayers 
# output coefficients, r-squared, residual standard error and p-value(F stat) 

# Calculate focal ("moving window") weight 
fw = focalWeight(resp, 2, "Gauss") 

# focal regression: 
lm.focal = function(x, y, ...) { 
    if(all(is.na(x) & is.na(y))) {NA} 
    else { 
    m = lm(y~x) 
    summary(m)$r.squared #r-squared 
    # summary(m)$coefficients #intercept and slope together 
    #---> Error in setValues(x, value) : cannot use a matrix with these dimensions 
    # summary(m)$coefficients[1] #intercept 
    # summary(m)$coefficients[2] #slope 
    # summary(m)$coefficients[8] #p-value 
    # summary(m)$sigma #residual standard error 
    } 
} 
#---> How to output all at once? 

lm.focal.out1 = localFun(resp, expl, w=fw, fun=lm.focal, na.rm=TRUE) 
plot(lm.focal.out1) 

#%%%%%%%%%%%%%%%%%%%%% 
## SECTION2 
#%%%%%%%%%%%%%%%%%%%%% 

# create test data 
set.seed(1) 
resp = expl1 = expl2 = expl3 = expl4 = raster(nrow=10, ncol=10) 
# x1 response variable 
resp = setValues(resp,runif(100,min=15,max=45)) 
# x3 explanatory variables 
expl1 = setValues(expl,runif(100,min=2,max=6)) 
expl1 = expl1 * resp 
expl2 = expl1 * resp/runif(100,min=1,max=4) 
expl3 = ((expl1 * resp)/1.5)/10 
expl4 = ((expl1 * resp)/runif(100,min=0.5,max=2))/100 
# add some NA values 
resp[1:5] = NA; expl1[1:5] = NA; expl2[1:5] = NA; expl3[1:5] = NA; expl4[1:5] = NA 

#stack RasterLayers 
stack1 = stack(list(resp=resp,expl1=expl1,expl2=expl2,expl3=expl3,expl4=expl4)) 
# par(mfrow=c(1,1)) 
plot(stack1) 
#.............................................................. 

# check global lm() results 
stack1.df = na.omit(as.data.frame(stack1)) 
head(stack1.df) 
stack1.lm = lm(resp ~ expl1+expl2+expl3+expl4, data=stack1.df) 
(stack1.lmSum = summary(stack1.lm)) 
stack1.lmSum$coefficients[1] 
stack1.lmSum$coefficients[2];stack1.lmSum$coefficients[3];stack1.lmSum$coefficients[4];stack1.lmSum$coefficients[5] 
stack1.lmSum$r.squared 
stack1.lmSum$sigma 
pf(stack1.lmSum$fstatistic[1], stack1.lmSum$fstatistic[2], stack1.lmSum$fstatistic[3],lower.tail = FALSE) 
#.............................................................. 

# lm function for focal {raster} with RasterStack 
# output coefficients, r-squared, residual standard error and p-value(F stat) 

# Calculate focal ("moving window") weight 
fw.s = focalWeight(stack1, 2, "Gauss") 

# focal regression with raster stack: 
lm.focal.stack = function(x, ...) { 
    if(all(is.na(x))) {NA} 
    else { 
    m = lm(x[1]~x[2]+x[3]+x[4]+x[5]) 
    summary(m)$r.squared #r-squared 
    # summary(m)$coefficients #intercept and slope together 
    #---> Error in setValues(x, value) : cannot use a matrix with these dimensions 
    # summary(m)$coefficients[1] #intercept 
    # summary(m)$coefficients[2] #slope 
    # pf(summary(m)$fstatistic[1], summary(m)$fstatistic[2], summary(m)$fstatistic[3],lower.tail = FALSE) #p-value 
    # summary(m)$sigma #residual standard error 
    } 
} 
#---> How to output all at once? 

lm.focal.stack.out1 <- focal(stack1, w=fw.s, fun=lm.focal.stack, na.rm=TRUE) 
#---> unable to find an inherited method for function ‘focal’ for signature ‘"RasterStack"’ 
#plot(lm.focal.stack.out1) 

#----------------------------------------------------------- 
> sessionInfo() 
R version 3.3.1 (2016-06-21) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 
Running under: Windows >= 8 x64 (build 9200) 

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

other attached packages: 
[1] raster_2.5-8 sp_1.2-3  

loaded via a namespace (and not attached): 
[1] rgdal_1.1-10 tools_3.3.1  Rcpp_0.12.5  grid_3.3.1  lattice_0.20-33 

Answer 1

Не уверен, что если вы все еще это нужно ответить, но у меня был тот же вопрос, и сделал функцию, называемую localFunStack, чтобы сделать работу вектора выхода из локальная функция как объект растрового объекта, с небольшим взломом, чтобы получить правильные имена слоев:

# localFun modified to write out a layer stack 
localFunStack <- function(x, y, ngb=5, fun, ...) { 

    compareRaster(x,y) 
    rasterList <- list() 
    nc1 <- 1:(ngb*ngb) 
    nc2 <- ((ngb*ngb)+1):(2*(ngb*ngb)) 

    if (canProcessInMemory(x, n=2*ngb)) { 
    vx <- getValuesFocal(x, 1, nrow(x), ngb=ngb) 
    vy <- getValuesFocal(y, 1, nrow(y), ngb=ngb) 
    v <- apply(cbind(vx, vy), 1, function(x, ...) fun(x[nc1], x[nc2], ...)) 

    for (j in 1:nrow(v)) { 
     if (length(rasterList) < j) { 
     rasterList[[j]] <- raster(x) 
     } 
     values(rasterList[[j]]) <- v[j,] 
    } 
    } 

    else { 
    tr <- blockSize(out) 
    pb <- pbCreate(tr$n, label='localFun', ...) 

    for (i in 1:tr$n) { 
     vx <- getValuesFocal(x, tr$row[i], tr$nrows[i], ngb=ngb) 
     vy <- getValuesFocal(y, tr$row[i], tr$nrows[i], ngb=ngb) 
     v <- apply(cbind(vx, vy), 1, function(x, ...) fun(x[nc1], x[nc2], ...)) 
     for (j in 1:nrow(v)) { 
     if (length(rasterList) < j) { 
      rasterList[[j]] <- raster(x) 
     } 
     rasterList[[j]] <- writeValues(rasterList[[j]], v[j,], tr$row[i]) 
     } 
    } 
    } 
    return(stack(rasterList)) 
} 

# local regression function 
lm.focal <- function(x, y, ...) { 
    if(all(is.na(x) & is.na(y)) || all(is.na(x)) || all(is.na(y))) {rep(NA, 8)} 
    else { 
    m <- lm(y~x) 
    coef <- summary(m)$coef 
    if (nrow(coef) == 1) { # Add NAs for cases where the response is constant 
     coef <- rbind(coef, rep(NA, 4)) 
     rownames(coef) <- rownames(nm) 
    } 
    coef <- as.vector(coef) 
    names(coef) <- c(outer(rownames(nm), colnames(nm) ,FUN=paste ,sep=" ")) 
    coef 
    # summary(m)$r.squared #r-squared 
    # summary(m)$sigma #residual standard error 
    } 
} 

lm.focal.out = localFunStack(expl, resp, ngb=5, fun=lm.focal, na.rm=TRUE) 

m <- lm(resp ~ expl) 
nm <- summary(m)$coefficients  
names(lm.focal.out) <- c(outer(rownames(nm), colnames(nm), FUN=paste, sep=" ")) 
plot(lm.focal.out)