0

I have a netcdf raster object with multiple layers (time: 30 years * 12 month = 360 steps).

The raster is constituted of multiple grid-cells (say 1000). I want to linearly detrend each month separately and each grid cell separately to obtain a residual raster object with same dimensions as the input.

For a long format data frame with columns [x,y,value,month,year], I would do it this way:

library(tidyverse)
data %>%
group_by(x, y, month) %>%
      nest() %>%
      mutate(model_trend = map(data,  ~ lm(value ~ year, data = .))) %>%
      mutate(augment_trend  = map(model_trend,  ~ broom::augment(.))) %>%
      unnest(cols = c(augment_trend)) %>%
      ungroup() %>%
      dplyr::select(x, y, year, month, .resid)

I am sure there are much more efficient and similarly easy ways of doing this operation using spatial packages such as terra, stars and so forth. Could you please advice on how to do that ?

Thanks

Omniswitcher
  • 368
  • 3
  • 13
Raed Hamed
  • 327
  • 1
  • 10

1 Answers1

0

I am not quite sure what you are doing as you provide no data and you do not use standard R. But here is how you can detrend cells in a SpatRaster using linear regression.

library(terra)
s <- rast(system.file("ex/logo.tif", package="terra"))   

n <- 1:nlyr(s)
f <- function(y) {
   residuals( lm(y ~ n) )
}
a <- app(s, f)

And a much faster approach with linear algebra (may not work with your data)

X <- cbind(1, n)
invXtX <- solve(t(X) %*% X) %*% t(X)

qfun <- function(y) {
    p <- invXtX %*% y
    y - (p[1] + p[2] * n)
}

app(s, qfun) 

#class       : SpatRaster 
#dimensions  : 77, 101, 3  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#names       :       red,      green,      blue 
#min values  : -3.833333, -17.000000, -3.833333 
#max values  :  8.500000,   7.666667,  8.500000 

You can use the above for all months by using tapp instead of app

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63