Linear regression between two raster images in R

Question

I need a linear regression for calculating an empirical parameter. L1 is a raster image, format .tif. L2 is a raster image as well, calculated beforehand. Both images have the same number of columns and rows.

The formula is: L1 = a + b * L2 which translates in R as:

lm(L1 ~ L2)

In a second formula I later need a nd b.

I am now facing the problem, that both raster contain NA values and I not sure how to build the function for the linear regression. I am not that familiar with R so I am stuck at this problem which might be rather simple. I think I have to use calc, but not sure how.

Edit: So far I have this code:

s = stack(L1,L2)
fun = function(x) {if (is.na(x[1])) { NA } else {lm(x[1] ~ x[2])$coefficients[2]}}

However, it takes a really long time calculating and does not come to an result

Answer 1

You would use calc if you wanted to do local regression, that is a separate regression for each grid cell (pixel). But that makes no sense in this case as you only have two rasters; and thus only one data point per grid cell.

In your case, you appear to want a global regression. You can that like this:

s <- stack(L1, L2)
v <- data.frame(na.omit(values(s)))
# this step may not be necessary
names(v) <- c('L1', 'L2')
m <- lm(L2 ~ L1, data=v)
m

if s is too large for that, you can do something like

v <- sampleRegular(s, 100000)  
v <- data.frame(na.omit(v))

etc.

Now with some data (and showing how to get residuals)

library(raster)
f <- stack(system.file("external/rlogo.grd", package="raster")) 
s <- stack(f)
names(s)
v <- data.frame(na.omit(values(s)))
m <- lm(red ~ green, data=v)
m

p <- predict(s, m)
residuals <- s$red - p