R: Annual composite raster based on the median: How to get the index of the original layer for each pixel?

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.

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 what I need to adjust to get one output layer with values from 1-4, and how to let the 'which'-function deal with NAs.

#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)

Solution

Thanks for providing good example data:

library(raster)
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[[2]][1:3] <- NA

Here are two functions that do the same. #1 is perhaps simpler but probably less efficient.

annualIndex1 <- function(ts, year) {
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))
    medianval <- calc(ts_year, median, na.rm = TRUE)
    dif <- abs(ts_year - medianval)
    which.min(dif)
}
x1 <- annualIndex1(s, 2002)

This one uses calc to combine multiple steps

annualIndex2 <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}  
x2 <- annualIndex2(s, 2002)

The results are the same (and there are no NAs)

all(values(x1) == values(x2))
#[1] TRUE

However, if you make a cell NA in all layers

s[5] <- NA

function 2 fails. It could be adjusted like this

annualIndex2b <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        if (all(is.na(x))) return( NA )
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}

x2b <- annualIndex2b(s, 2002)

You use the term "mask", but really you are creating an index, can you can use the appropriate cells with strackSelect like this:

ts_year <- subset(ts, (grep(year, names(ts))))
v <- stackSelect(ts_year, x2)