Pixel-wise time series trend anaylsis with NDVI (GIMMS) and R

The GIMMS dataset is currently offline and the new GIMMS3g will soon be released, but it does not really matter which type of data is used for this purpose. It can also be SPOT VGT, MODIS or anything else as long as the temporal resolution is high and the time frame is long enough to detect significant trends. The purpose is to do a pixelwise trend analyis and extract only significant trends over a certain period of time fora selected region. Everything is done using open and free R software. Input data are continuous NDVI images, in this case it’s GIMMS with bi-monthly images, so 24 per year for the period 1982-2006.

so let’s load all Geotiffs lying in the GIMMS directory (should be  600 in this case):

sg = stack(list.files(pattern='*.tif'))
gimms = brick(sg)

first, let’s create annual sums, as autocorrelation might be present with monthly values. To keep NDVI as the unit, the results are devided by 24:

fun <- function(x) { 
 gimms.ts = ts(x, start=c(1982,1), end=c(2006,24), frequency=24)
 x <- aggregate(gimms.ts) 
gimms.sum <- calc(gimms, fun)

then the slope is calculated to get the direction and magnitude of trends, multiplied by the number of years to get the change in NDVI units:

time <- 1:nlayers(gimms.sum) 
fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[2] }}
gimms.slope=calc(gimms.sum, fun)

now we need to see which trends are significant. Thus we first extract the p-value:

fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[8] }}
p <- calc(gimms.sum, fun=fun)
plot(p, main="p-Value")

then mask all values >0.05 to get a confidence level of 95%:

m = c(0, 0.05, 1, 0.05, 1, 0)
rclmat = matrix(m, ncol=3, byrow=TRUE)
p.mask = reclassify(p, rclmat)
fun=function(x) { x[x<1] <- NA; return(x)}
p.mask.NA = calc(p.mask, fun)

and finaly mask all insignificant values in the trend map, so we only get NDVI change significant at the 95% level:

trend.sig = mask(gimms.slope, p.mask.NA)
plot(trend.sig, main="significant NDVI change")

The result could look like that:

significant NDVI change (1982-2006) using integrated GIMMS

significant NDVI change (1982-2006) using integrated GIMMS


25 thoughts on “Pixel-wise time series trend anaylsis with NDVI (GIMMS) and R

      • Hi Martin,
        i appreciated your good work, it helps me a lot. But i do have some more question, since i am new for R. So could you respond me through my email:s6genega@uni-bonn.de
        i have been working on NDVI for assessing desertification in Ethiopia for 27 yrs (1982-2008) so after i got slope (Trend) each NDVI pixel value of the study area , i would like to test how much significance it is spatially.

  1. Hi Martin,
    First of all, thanks for sharing your work on web. I am also working on time series analysis of NDVI data. I have gone through your code, I understand everything except ‘$coefficients[2] and $coefficients[8]‘. Could you please explain me what are these two things? My email ID is ekta.geo@gmail.com.
    Thanks for your time!

    • if you do a linear regression with the “lm” function, it stores several coefficients, like the slope, the intercept, the significance, etc, in our case they are stored in “m” and can be shown by “summary(m)”. If you want to extract a specific coefficient, like the slope, you can do that with “summary(m)$coefficient[2]”, where [2] means slope. [8] is the significance (p), and you can also extract the other coefficients with other numbers….

      hope this helps,

  2. hi Martin,
    first of all thanks for all those important codes that you have used in your study. I am also a new user in R. i am working on classification of crop rotation system over Punjab and Haryana in India using MODIS time series EVI image. I have stacked them, and now i want to run SAVGOL filter for smoothing the temporal signal for each pixel in the entire scene. i have used the code
    sg<- sgolayfilt(EVItimeseries, 3, 5)
    but it returns error.
    I also want to know that how i can execute the FFT on the same stacked image after smoothing.
    is it possible to reply in my mail (prasenjit.ac@gmail.com) also? thanks in advance for your help

    • it’s possible to use the SavGol filter as a funtion, try this, assuming you have a raster brick or stack called MODIS for one year. Here I interpolate NA values with mean values:


      fun <- function(x) {
      z=substituteNA(v, type="mean")
      MODIS.ts2 = ts(z, start=c(2011,1), end=c(2012,23), frequency=23)
      x=sgolayfilt(MODIS.ts2, p=1, n=3, ts=30)

      MODIS.filtered <- calc(MODIS, fun)

      • Thank you Martin for your reply.
        Things become weird after i run the code. Let me tell you how I stack those time series images. I have used ERDAS 9.1 image processing software to stack the MODIS 16 day NDVI time series images from May 9 2011 to Apr 22 2012. then I subset the image according to my desired area. After doing this, I have import the .img file into R by using “stack(myfile.img)” command. This is followed by the steps that you had given me. It fails to return the result I want. Instead of smoothing it breaks the periodic nature of the NDVI signals. Can you please help me to solve the issue.
        with regards,

      • i think you need to adjust the time frame in this line:
        start=c(2011,1), end=c(2012,23), frequency=23
        where frequency is the number of images per year, start is the beginning image and end the last one. So assuming you have 16 day composites, in your case I think
        start=c(2011,9), end=c(2012,8), frequency=23
        should do it

      • Hello Martin,
        Still facing the same issue. can I send you a small patch of the the image file that I am using?
        it would be better for me to understand then. my email…. prasenjit.ac@gmail.com
        with regards,

  3. Thank you for your article,

    I am currently working with several NDVI products from different satellite and trying to asses the consistency of those data. Do you think it it worth to apply savitzky golay filter to GIMMS time series data as they have no NA value?

    thanks before.

      • Hi Martin!
        Many thanks for sharing these codes which are absolutely helpful.

        Just wondering why you said GIMMS images do not need smoothing? What about if you process them in TIMESAT? They will have a smoothing process there, aren’t they?

        Do you think is better use NDVI3g data as downloaded? why?

        All the best!

      • Hi Betsa,
        the GIMMS images are not that noisy, but if you process them in TIMESAT it will work for sure and you can extract some metrics, but the smoothing effect will be rather low. Don’t use the NDVI3g yet, wait for the new version 1 which will be available soon…

  4. Thanks for putting this together Martin, this is very helpful to me.

    I am new to R and am reaching a memory error: “Reached total allocation of 7168Mb: see help(memory.size)”. It gets the error when running: “gimms.sum <- calc(gimms, fun)".
    I was wondering if you have any tricks to get around this.

    I am working with the NOAA STAR NDVI data. It is also based on AVHRR but at 4km2 resolution and from mid 1981-present. I ran your script, after averaging to yearly myself for 34 years (1982-2015).


  5. Hi Martin,

    Thanks for sharing this code. I’m trying to use this code for an “irrigulay Space time-series”. Meaning that I have some missing data. I have NDVI data at annual level and have some gaps between the years.

    How should I modify the code?

    Thank you

  6. Hi Martin, I’m following ur example code and getting error message: “invalid time series
    parameters specified” Pls suggest a solution.

    (gimmsstack is: 48 gimms images in an ordered stack)

    > gimmsTS = ts(gimmsstack, start=c(2008,1), end=c(2009,2), frequency=24)
    Error in ts(gimmsstack, start = c(2008, 1), end = c(2009, 2), frequency = 24) :
    invalid time series parameters specified

    (600 images)
    > gimmsTS = ts(gimmsstack, start=c(1982,1), end=c(2006,24), frequency=24)
    Error in ts(gimmsstack, start = c(1982, 1), end = c(2006, 24), frequency = 24) :
    invalid time series parameters specified

  7. Very interesting Martin. Thank you.
    I haven’t understood this line of code:
    fun=function(x) { x[x<1] <- NA; return(x)}
    p.mask.NA = calc(p.mask, fun)
    You want to mask all values <0.05, so why " x[x<1]"?

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