6

I'm interested in doing a pixel-wize trend analysis of a series of NDVI images from Landsat in R.

The problem is that my images are not exactly happening at regular intervals. But decent imagery of my study area is sporadic and not regular at all (July-21-1992, May-20-1996, ...). Therefore, I can't directly use any of the timeseries packages (ts, greenbrown) to do this. The default assumption of these packages/functions is that your raster stack has a regular time interval (daily, monthly, ...). Here is an example:

Find and stack all the layers:

Library(greenbrown)
all <- list.files("C:/Users/R_test", full.names = TRUE, pattern = "*.tif")
st <- stack(all)

Perform the pixel-wize analysis: 

trends <- TrendRaster(st, start = c(1992, 1), freq=1, method=c("ATT"), breaks = 1)

Therefore, my question is how can I implement these gaps into the timeseries before doing any analysis?

I have tried the zoo package and managed to create a vector that has NA for the dates that I dont have imagery, but I'm not sure how to insert the zoo object back into the trend analysis function (here the greenbrown package).

Regards

Improve this question

1 Answer 1

Reset to default
6

I guess I found a way to my own problem:

Basically, I parsed the time variable from the filename and added it back to regression model. Below is the peice of code I used. 

library(raster)


all <- list.files("/home/R_test/", full.names = TRUE, pattern = "*.tif")
fn <- list.files("/home/R_test/", full.names = FALSE, pattern = "*.tif")

#Stack rasters
r <- stack(all)

## Parsing out the date. My file names are
## like this: 1993154.tif which is the year and Julian day
Date <- strptime(fn, "%Y%j")
Date <- as.Date(Date)


## Here is to insert the time component back to the
## linear model. This one is to calculate the slope
fun_slope <- function(y) { 
    if(all(is.na(y))) {
      NA
    } else {
    m = lm(y ~ Date); summary(m)$coefficients[2] 
   }
  }
## and this one is to calculate the p-value
fun_pvalue <- function(y) { 
    if(all(is.na(y))) {
    NA
    } else {
     m = lm(y ~ Date); summary(m)$coefficients[8] 
  }
}

slope <- calc(r, fun_slope)
pvalue <- calc(r,fun_pvalue)

The output is slope and p-value of the fitted line to each pixel enter image description here enter image description here

Improve this answer
1
  • 3
    As long as it is a constant vector length, calc() will recycle a resulting vector of >1 into a raster stack. As such, you do not need separate functions return slope and p-value. You can simply use 'summary(m)$coefficients[c(2,8)]' This will return a raster stack where the first raster is the slope and the second the p-value. This avoids running calc() twice which, on large rasters, is quite a bit of overhead.
    – Jeffrey Evans
    Commented Feb 1, 2018 at 16:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.