0

I want to perform a moving window regression on every pixel of two raster stacks representing Band3 and Band4 of Landsat data. The result should be two additional stacks, one representing the Intercept and the other one representing the slope of the regression. So layer 1 of stack "B3" and stack "B4" result in layer 1 of stack "intercept" and stack "slope". Layer 2 of stack B3 and stack B4 result in layer 2,.... and so on.

I already came along the gwrfunction, but want to stay in the raster package. I somehow know that focal must be included in order to set my moving window (which should be 3x3 pixels) and somehow a linear model like: lm(as.matrix(b3)~as.matrix(b4)) although I don't think that this gets me the values pixelwise...

Instead of a rasterstack a layer by layer approach is also possible. (So it must not necessarily be a rasterstack of Band3.

Has anyone a glue how to program this in R?

phil652
  • 1,484
  • 1
  • 23
  • 48
user2978751
  • 57
  • 1
  • 10

1 Answers1

2

Here is one approach, adapted from ?raster::localFun

set.seed(0)
b <- stack(system.file("external/rlogo.grd", package="raster"))
x <- flip(b[[2]], 'y') + runif(ncell(b))
y <- b[[1]] + runif(ncell(b))

# local regression:
rfun <- function(x, y, ...) {
    d <- na.omit(data.frame(x, y))
    if (nrow(d) < 3) return(NA)
    m <- lm(y~x, data=d)
    # return slope
    coefficients(m)[2]
}

ff <- localFun(x, y, fun=rfun)
plot(ff)

Unfortunately you will have to run this twice to get both the slope and intercept (coefficients(m)[1]).

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • 1
    When setting `ngb` as 3 rows and 3 columns it users the intersection (i.e. 9 cells). In your case you can use `ngb=3` or `ngb=c(3,3)`. – Robert Hijmans May 22 '15 at 16:22
  • Thanks Robert! - I would have one more question: as my actual raster stacks of B3 and B4 include NA values (which are varying in the layers but are on the same position between the rasterstacks - layer 1 of B3 has same NA positions as layer 1 in B4). I thought this shouldn't be any problem via `m <- lm(x~y, na.action=na.exclude)` but it gives me following error: `Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases Called from: eval(substitute(browser(skipCalls = pos), list(pos = 9 - frame)), envir = sys.frame(frame))` – user2978751 May 23 '15 at 10:29
  • A bit more than half of the raster is usually NA. I suppose its not the best condition for performing a moving window regression - I guess it should be possible. – user2978751 May 23 '15 at 10:34
  • Do you know if it is possible to apply a "circular" moving window. (Like with `focal()` - specifying a matrix of {0,1} or using something like `focalWeight()` - both unfortunately didn't work as an input for argument `ngb`). Thanks a lot! – user2978751 Jun 01 '15 at 11:52