Using a custom regression model with terra

Question

I would like to fit a custom regression model (mblm package) to each cell in a time series raster (each layer is a year). My raster manipulation is all happening with terra, so I'd prefer to continue in that package world if possible.

Example dataset:

d<-rast(lapply(1:10, FUN = function(x){
  rast(matrix(rnorm(100, mean=x), ncol=10))
}))
names(d)<-1:10

I've looked through the regress function, but I think I might be misunderstanding something about how it works. For example, if I create a hypothetical 10-layer year raster:

d1<-rast(lapply(names(d)[1:10], function(x) {rast(matrix(x, nrow=10, ncol=10))
}))
names(d1)<-1:10

and use the built-in lm within regress:

plot(regress(d, d1))

The slope raster (labelled 'x' comes out as all zero).

Eonema · Accepted Answer · 2024-05-17 06:02:42Z

You use names(d)[1:10] to fill the layers of d1. However, these values are characters*, not numbers:

 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"

One way to resolve this is by changing x to as.numeric(x):

d1<-rast(lapply(names(d)[1:10], function(x) {rast(matrix(as.numeric(x), nrow=10, ncol=10))
}))
names(d1)<-1:10

This gives you model parameters like:

names       : (Intercept),         x 
min values  :   -1.159427, 0.6790784 
max values  :    1.454451, 1.2057739

Note though, that since your time values are the same for all cells, you can avoid creating d1 altogether, and just use a single vector instead:

regress(d, as.numeric(names(d)))

And get the same thing.

* If you ran one of these cells through lm, it would treat the independent variable as categorical, and give you nine different slopes. I don't know why regress doesn't do the same thing. Maybe it isn't built for categorical variables, or maybe it's something about how the values are stored in a SpatRaster.

Robert Hijmans · Accepted Answer · 2024-05-21 05:56:36Z

0

Here are two ways in which you can use terra::regress.

A single SpatRaster and an index (that is, "x" in y~x is fixed)

library(terra)
set.seed(1)
r <- rast(ncol=10, nrow=10, nlyr=10, vals=runif(1000))

regress(r, 1:nlyr(r))
#class       : SpatRaster 
#dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#names       : (Intercept),           x 
#min values  :  0.02779796, -0.07366238 
#max values  :  0.90875803,  0.06692299

And with another SpatRaster (that is, (that is, "x" in y~x is changes by cell)

p <- rast(ncol=10, nrow=10, nlyr=10, vals=runif(1000))

regress(r, p)
#class       : SpatRaster 
#dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#names       : (Intercept),          x 
#min values  :  0.08374815, -0.7327174 
#max values  :  0.87454291,  0.6017231

You can achieve the same with lapp and lm (or other regression method) like this

x <- lapp(sds(p, r), \(x, y) {
    sapply(1:nrow(x), \(i) {
        coefficients(lm(b~a, data=data.frame(a=x[i,], b=y[i,])))
        }) |> t()
    })

x
#class       : SpatRaster 
#dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#names       : (Intercept),          a 
#min values  :  0.08374815, -0.7327174 
#max values  :  0.87454291,  0.6017231

The code is a bit gnarly because lm is not a vectorized function.

It is advisable to check some cells. E.g.,

x[1]
#  (Intercept)          a
#1   0.6239804 0.06438924

coefficients(lm(y~x, data=data.frame(x=unlist(p[1]), y=unlist(r[1]))))
#(Intercept)           x 
# 0.62398039  0.06438924

edited May 21, 2024 at 5:56

answered May 17, 2024 at 12:24

Robert Hijmans

49.5k4 gold badges62 silver badges77 bronze badges

2 Comments

FrostyBottom Over a year ago

Thank you for those suggestions. That works perfectly. The second part of the question, which in hindsight I did not state clearly enough: is it possible to use a completely different regression model in place of lm? Looking at the source code for regress, lm fits via ols look pretty baked in. In my case, if I wanted to apply a Theil Sen regression model from the mblm package, is there a way to apply that to each cell instead?

Robert Hijmans Over a year ago

I have now expanded my answer to capture that.

Collectives™ on Stack Overflow

Using a custom regression model with terra

2 Answers 2

Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related