1

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.

4

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)
  • Would this script take into account those pixels that lie at the border of the two rasters PPT and ELEV and are not surrounded exactly by 3×3 pixels or in other words, 8 pixels with them in the center? The ones at the edge would require linear regression of 2×3 pixels or 3×2 pixels or sometimes even 2×2 pixels in the two rasters? Thanks a lot Jeffrey Evans for your quick and considerate reply. I have worked with MATLAB but I am a beginner in R programming. Therefore, I am trying to understand and apply it to my analysis but I had to clear this issue with pixels at the border. – – Alejandro Apr 11 '18 at 18:20
  • This would be subject to edge effect and is an inherent limitation of focal functions. You could, in theory, pad edge values but this would invalidate the regression. It is common practice to apply a focal function and then pull the extent back equal to 1/2 the focal window size thus, excluding edge pixels. This is why it is good pratice to account for this when defining your analysis extent. – Jeffrey Evans Apr 11 '18 at 18:27
  • I wonder if I may take the two rasters PPT and ELEV occupying a larger area than my study area and apply the function to calculate the regression slope and obtain the raster of regression slope (GRADIENT) and later crop this raster GRADIENT using shapefile of my study area to avoid this inherent limitation of focal function. Forgive me if you meant the same in your comment and I could not get it.Thanks again for your prompt reply and your kind consideration. – Alejandro Apr 11 '18 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.