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

R² annual sums SPOT VGT NDVI and TRMM rainfall, 1998-2010 for Linguere, Senegal

R² annual sums SPOT VGT NDVI and TRMM rainfall, 1998-2010 for Linguere, Senegal

Advertisements

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

  1. 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!

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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s