# Pixel-wise regression between two raster time series (e.g. NDVI and rainfall)

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)){ NA } else { lm(x[1:120] ~ x[121:240])\$coefficients }}
slope <- calc(s, fun)```
```fun=function(x) { if (is.na(x)){ NA } else { lm(x[1:120] ~ x[121:240])\$coefficients }}
intercept <- calc(s, fun)```
```fun=function(x) { if (is.na(x)){ 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:

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. R² annual sums SPOT VGT NDVI and TRMM rainfall, 1998-2010 for Linguere, Senegal

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

1. Betsa says:

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)){ NA } else { lm(x[1:360] ~ x[361:720])\$coefficients }}
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

1. Miles says:

Betsa,

I was struggling with this too. Try using:

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

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

1. Mohammad Aref says:

my problem still remain even using this modification

2. Isis Dias says:

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

2. Isis Dias says:

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.

3. Tawatchai NA-U-DOM says:

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. 🙂

2. Sujit says:

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?

3. Bhavani says:

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 works only for one parameter. Please kindly request you to sort the problem.
Thanking you
with regards
Bhavani

4. Arian Correa says:

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)){ 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