0

What I want:

I have a list of rasterstacks which are timeseries (300+ layers) for different bands. The timeseries are irregular, and therefore I want to create annual composites based on median ndvi. Thus from the available images for one year, the pixel with the (close to) median ndvi is chosen. For the other bands I want to create annual composites based on this ndvi composite. I try to make a mask in which each pixel has the value of the index of the image used for the median ndvi composite. I will apply this on the other bands, so I have 'the same' annual composite for each band per year.

Where I am at

Unfortunately, I am stuck on making the mask. I created some dummy data and somehow it returns two indexrasters (one with values 1-3, the other 2-4), while I expected one (with values 1-4). Also my function cannot handle NA values and adding na.rm to the calc function does not solve this. I am wondering whats going wrong and how to solve this.

#dummy data:
r <- raster(nrow=5, ncol=5)
set.seed(20181114)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))
names(s) <- c("X2001", "X2002a", "X2002b", "X2002c", "X2002d")
#s$X2002a[2] <- NA 

AnnualMask <- function(ts, year){
  year <- as.character(year)
  ts_year <- subset(ts, (grep(year, names(ts))))
  indexraster <- calc(ts_year, function(x){
    medianval <- median(x, na.rm = TRUE)
    which(abs(x - medianval) == min(abs(x - medianval)))
           })
  return(indexraster)
}

mask2002 <- AnnualMask(s, 2002)
plot(mask2002)

What do I need to do?

1

In a median, if the number of cases in consideration is even, then the result is the mean of two middle results. As such, your which() will always return two values, hence you get two layers, one with the index of one closest layer to the median, and one with the other.

It works as you expected if you make the code a bit simpler using which.min(), which only returns the first minimum value:

indexraster <- calc(ts_year, function(x){
    medianval <- median(x, na.rm = TRUE)
    distances = abs(x - medianval)
    return(which.min(distances))
       })

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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