CODEWITHSUNDEEP
rraster
nate-m
nate-m

Reputation: 557

Remove linear trend from raster stack R

Trying remove the linear trend (detrend) from a monthly precipitation raster stack for the US from 1979-2015 (https://www.northwestknowledge.net/metdata/data/monthly/pr_gridMET.nc). These data are large enough that using those data as an example would be a bit unruly here so I am going to use the data from the raster package for sake of efficiency. The working model I have currently is to use `raster"::calc`` on a linear model and pull the residuals. My understanding is that those residuals are the detrended series, but I am not 100% sure that is correct. The code I am using is as follows:

library(raster)

fn <- raster(system.file("external/test.grd", package="raster"))
fn2 <- fn+1000
fn3 <- fn +500
fn4 <- fn +750
fn5 <- fn+100

fns <- stack(fn, fn2, fn3, fn4, fn5)

time <- 1:nlayers(fns)

# Get residuals to detrend the raw data
get_residuals <- function(x) {
  if (is.na(x[1])){ 
    rep(NA, length(x)) } 
  else {
    m <- lm(x~time)
    q <- residuals(m)
    return(q)
  }
}

detrended_fns <- calc(fns, get_residuals) # Create our residual (detrended) time series stack

I feel like I'm missing something here. Can anyone confirm that I'm on the right track here? If I'm not any suggestions on how to properly detrend these data would be helpful! thanks!

Upvotes: 0

Views: 1462

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47156

The residuals remove the slope and the intercept and you get anomalies. Perhaps you only want to remove the slope? In that case you could add the intercept to the residuals in get_residuals

q <- residuals(m) + coefficients(m)[1]

Or better: 

q <- residuals(m) + predict(m)[1]

So that you use year 1 (and not year 0) as the base, and it would also work if time is, say, 2000:2004

You could also take the last year, mid year, or average as base.

Upvotes: 1

Related Questions