Linear regression between every 3×3 pixels between two rasters using R

Question

I have two rasters of precipitation and elevation, PPT and ELEV. They have equal number of rows and columns.I want to create a new raster of precipitation lapse rate (GRADIENT) with same number of rows and columns as PPT and ELEV. I want to compute this raster by performing linear regression between every 3×3 pixels of the two rasters (PPT and ELEV) such that each pixel of the GRADIENT will hold the lapse rate value obtained from linear regression of the corresponding 3×3 pixels in PPT and ELEV that surrounds that pixel.

Answer 1

I assume that by "lapse rate" you mean the regression slope. You are essentially describing a focal regression and since there is no multivariate version of raster::focal, there is no canned way to do what you are after. This could be accomplished through a series of focal matrix algebra functions. However, you can also leverage the read/write functions that the raster package has via gdal.

First, we add some packages and create some data (all-be-it, somewhat nonsensical). The final object we are after is the r raster stack.

library(raster)
library(gstat)
library(sp)

data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y

v1 <- variogram(log(copper) ~ 1, meuse)
x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))
copper <- krige(copper ~ 1, meuse, meuse.grid, x1, nmax = 30)
gridded(copper) <- TRUE
copper@data = as.data.frame(copper@data[,-2])

v2 <- variogram(log(elev) ~ 1, meuse)
x2 <- fit.variogram(v2, vgm(.1, "Sph", 1000, .6))
elev <- krige(elev ~ 1, meuse, meuse.grid, x2, nmax = 30)
gridded(elev) <- TRUE
elev@data <- as.data.frame(elev@data[,-2])
elev@data[,1] <- elev@data[,1]

r <- stack( raster(copper), raster(elev) )

Now that we have some data, we can set a focal window size and implement a loop that will iterate through the rows of the rasters. The basic idea is that we take the focal values from each raster in the stack, regress the resulting vectors and assign the regression slope back to an empty raster. The workhorse function here is getValuesFocal which returns a list object containing data.frames for each raster. The rows of the data.frames are the focal values for each cell in a given row for the specified NxN window. I provided the basic code here but, this could certainly be optimized by extracting more data, than just a row at a time, and using mapply.

w <- c(3,3)
lapse <- r[[1]]
  lapse[] <- NA
  for( rl in 1:nrow(r) ) { 
    v <- getValuesFocal(r[[1:2]], row=rl, nrows=1, ngb = w, array = FALSE)
      fit <- rep(NA,nrow(v[[1]]))
        for(i in 1:nrow(v[[1]]) ) {
          xy <- na.omit( data.frame(x=v[[1]][i,],y=v[[2]][i,]) )    
          if( nrow(xy) > 4 ) {
            fit[i] <- coefficients(lm(as.numeric(xy$y) ~ as.numeric(xy$x)))[2]
              if(is.null(fit)) fit = 1
          } else {
            fit[i] <- NA 
         }
     }
    lapse[rl,] <- fit
  }

lapse
plot(lapse)