Be cautious using GIMMS3g v0 for trend analysis!

GIMMS3g NDVI is widely used to assess vegetation trends from local to global scale. And even though it is the best long term dataset available right now (July 2015), people should be aware that a serious error affects the dataset in semi arid areas, like the Sahel, which has severe impacts on trend analysis. The dry season values around 2004 suddenly drop, which is not happening in reality. It is definitely related to the sensor change from NOAA16 to NOAA17, and also the last change to NOAA18 has an impact on the time series. The good thing is that the wet season values should be usable. The developers are aware about this problem and a new version is on it’s way. To illustrate the issue, find attached the GIMMS3g NDVI (v0) averaged over the Sahel belt, note the drop at the sensor change (first blue line) and the recovery at the next sensor change:


To further illustrate how serious this problem affects trend analysis, the same area is shown in the VOD dataset, which is found to be reliable:vod

a rough overlay highlights the different directions of annual trends:


The problem can be solved by excluding the dry season and using the small integral, calculated in TIMESAT:


see further literature:

  • Tian, F., Fensholt, R., Verbesselt, J., Grogan, K., Horion, S., Wang, Y., 2015. Evaluating temporal consistency of long-term global NDVI datasets for trend analysis. Remote Sensing of Environment 163, 326–340. doi:10.1016/j.rse.2015.03.031
  • Horion, S., Fensholt, R., Tagesson, T., Ehammer, A., 2014. Using earth observation-based dry season NDVI trends for assessment of changes in tree cover in the Sahel. International Journal of Remote Sensing 35, 2493–2515. doi:10.1080/01431161.2014.883104
  • Pinzon, J.E., Tucker, C.J., 2014. A Non-Stationary 1981–2012 AVHRR NDVI3g Time Series. Remote Sensing 6, 6929–6960. doi:10.3390/rs6086929
  • Liu, Y.Y., van Dijk, A.I.J.M., McCabe, M.F., Evans, J.P., de Jeu, R.A.M., 2013. Global vegetation biomass change (1988–2008) and attribution to environmental and human drivers. Global Ecology and Biogeography 22, 692–705. doi:10.1111/geb.12024
  • Jonsson, P., Eklundh, L., 2004. TIMESAT–a program for analyzing time-series of satellite sensor data* 1. Computers & Geosciences 30, 833–845.

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