Doing a pixel-wise regression between two raster time series can be useful for several reasons, for example:

find the relation between vegetation and rainfall for each pixel, e.g. a low correlation could be a sign of degradation

derive regression coefficients to model the depending variable using the independend variable (e.g. model NDVI with rainfall data)

check the consistency between different sensors, e.g. MODIS and SPOT VGT

adjust the data recorded by different sensors to each other

It’s easy and fast in R. Let’s say we have two raster stacks or bricks which contain 120 rasters each, one with NDVI and one with rainfall (named “NDVI” and “rainfall”). We put everything in one single stack and derive rasters showing the slope, intercept and R² of the two time series as pixel values:

library(raster)
s <- stack(NDVI, rainfall)
fun=function(x) { if (is.na(x[1])){ NA } else { lm(x[1:120] ~ x[121:240])$coefficients[2] }}
slope <- calc(s, fun)

fun=function(x) { if (is.na(x[1])){ NA } else { lm(x[1:120] ~ x[121:240])$coefficients[1] }}
intercept <- calc(s, fun)

fun=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[1:120] ~ x[121:240]);summary(m)$r.squared }}
r.squared <- calc(s, fun)

It only works if your two stacks have the same dimensions and NA values, if not, you can try something like that:

z = mask(NDVI, rainfall)
s <- stack(z, rainfall)

and do the analysis as described above.

This example shows R² between annually summed NDVI (SPOT VGT) and TRMM rainfall from 1998-2010 for the Linguere district in Senegal, and we see that major differences exist, which need to further examined.

9 thoughts on “Pixel-wise regression between two raster time series (e.g. NDVI and rainfall)”

Hi Martin,
Thanks for posting such a useful material. I was wondering how to get the significance values for this. I’m trying the following:
fun = function(x) { if (is.na(x[1])){ NA } else { lm(x[1:360] ~ x[361:720])$coefficients[8] }}
significance <- calc(s, fun)
plot(significance)
I'm getting a raster with NA's, however the other coefficients seem to be good. Do you know how to get a map with statistical significance?
Many thanks for this!
All the best,
Betsa

my problem still remain even using this modification

Hi Miles,
I am still struggling with this, I mean, this script have worked for me and I could plot the map. But I am not understanding what exactly it shows me; is this a map of p-value?
I will be very thankfull if you could help me.
Best regards,
Isis Dias

Hi Betsa,
I am also trying to get a map with statistical significance. Did you have any positive result applying a new script in your R code?
Best regards,
Isis.

Hi, Martin
This post is great. But still I am encountering a problem while executing the code. It was saying ‘Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases’. I think this code handles NA value issues. Can you tell me what may be the problem?

Hi Martin,
I am Bhavani , Research scholar.
The post is very informative and useful. I have executed the same procedure for my research area of interest NDVI with rainfall. Recently I am trying to do multiple regression at pixel level NDVI~ rainfall+ temperature. The model execute well. But how to pick the slopes values for two parameters. As $coefficients[2] works only for one parameter. Please kindly request you to sort the problem.
Thanking you
with regards
Bhavani

Hi Martin,
What a useful post!! I just discovered your page and it’s amazing. Just wondering if I can use a gam procedure instead lm approach, since few times normal assumptions are totally full filled.
Here’s my code but its seems there’s something wrong
fun=function(x) { if (is.na(x[1])){ NA } else { m <- gam(x[1:100] ~ s(x[101:200])); summary(m)$dev.expl}}
I will be very thankful if you could help me.
Best regards,
Arian CD

Hi Martin,

Thanks for posting such a useful material. I was wondering how to get the significance values for this. I’m trying the following:

fun = function(x) { if (is.na(x[1])){ NA } else { lm(x[1:360] ~ x[361:720])$coefficients[8] }}

significance <- calc(s, fun)

plot(significance)

I'm getting a raster with NA's, however the other coefficients seem to be good. Do you know how to get a map with statistical significance?

Many thanks for this!

All the best,

Betsa

Betsa,

I was struggling with this too. Try using:

fun=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[ ] ~ x[ ]);summary(m)$coefficients[8]}}

Note the inclusion of "summary", which seems to do the job!

my problem still remain even using this modification

Hi Miles,

I am still struggling with this, I mean, this script have worked for me and I could plot the map. But I am not understanding what exactly it shows me; is this a map of p-value?

I will be very thankfull if you could help me.

Best regards,

Isis Dias

Hi Betsa,

I am also trying to get a map with statistical significance. Did you have any positive result applying a new script in your R code?

Best regards,

Isis.

If you mean p-value which you got from r.square, you may try;

r.0 <- corLocal(r.stack.0L[[1:12]], r.stack.0L[[13:24]], test=TRUE )

when;

r.stack.0L[[1:12]] = NDVI raster stack

r.stack.0L[[13:24]] = Precipitation raster stack

Finally, you will get the map of R square and p-value

I hope it could help. 🙂

Hi, Martin

This post is great. But still I am encountering a problem while executing the code. It was saying ‘Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases’. I think this code handles NA value issues. Can you tell me what may be the problem?

Hi Martin,

I am Bhavani , Research scholar.

The post is very informative and useful. I have executed the same procedure for my research area of interest NDVI with rainfall. Recently I am trying to do multiple regression at pixel level NDVI~ rainfall+ temperature. The model execute well. But how to pick the slopes values for two parameters. As $coefficients[2] works only for one parameter. Please kindly request you to sort the problem.

Thanking you

with regards

Bhavani

Hi Martin,

What a useful post!! I just discovered your page and it’s amazing. Just wondering if I can use a gam procedure instead lm approach, since few times normal assumptions are totally full filled.

Here’s my code but its seems there’s something wrong

fun=function(x) { if (is.na(x[1])){ NA } else { m <- gam(x[1:100] ~ s(x[101:200])); summary(m)$dev.expl}}

I will be very thankful if you could help me.

Best regards,

Arian CD