Can't Calculate pixel-wise regression in R on raster stack with fun

user_8990027 0
Viewed 1038

I am working with rasters and I've a RasterStack with 7n layers. I would like to calculate pixel-wise regression, using formula beneath. I was trying to do it with raster::calc, but my function failed with message :

'Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases.'

But all rasters are OK, and contain numbers (not only NAs), I can plot it, and I can calculate general linear regression with formula

 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))

But when I stack layers with

allData = stack(MK_trend.EVI.sig_Only,stack.dtr,stack.wet,stack.pre,stack.vap,stack.tmp,stack.pet)

and try calc function

    # 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)
    }
}

I see the error message.
I am pretty new in R and progranning, so, maybe, there is some silly mistake? I would appreciate any hint in order to make the processing work.

calc raster regression r function
1 Answers

You can use calc for pixel-wise (local) regression, but your formula seems to suggest you want something else (a global model).

If the regression were pixel wise, you would have an equal number of x and y values for each cell, and you can use calc. See ?calc for examples.

Instead you have 1 y (independent) and 6 x (dependent) variable values for each cell. This suggests you want a global model. For that, you can do something like this:

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)

This requires loading all values into memory. If you can not do that, you could take a large regular sample. See sampleRegular

And to show why your approach does not work:

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

Even though I do not get the error message you report (but all R^2 values are zero).

user_6352450
Related