The purpose is to compare whether the indexes of names from the result of the stackApply function correspond to the indexes that were defined as parameter input. An alternative method (ver_median function) is used for comparison.

calc_median_wrong_results.R

# Calculate median per day of week. 

library(lubridate)
library(raster)
library(magrittr)

r <- raster(nrows = 300, ncols = 300, res = 500, 
            xmn = 0, 
            xmx = 150000, 
            ymn = 0, 
            ymx = 150000)

n=3000

mydates <- seq(as.Date("2012-01-19"), by = "day", length.out = n)
set.seed(5);
mystack <- raster::stack(lapply(1:n, function(i) {r<-raster::setValues(r,floor(runif(ncell(r), min=1, max=1000)));return(r)})) %>% 
  setNames(mydates) %>% 
  raster::setZ(mydates)



myindices <- lubridate::wday(raster::getZ(mystack), week_start = 1) #  "Monday"=1, 

### stackApply method 
### 
stackapply_median<- raster::stackApply(mystack, indices=myindices,fun=function(x,...){(median(x, na.rm = T))})




#### my verification
#### 
ver <- function(myday, s){
  sub_s <- raster::subset(s,which(lubridate::wday(raster::getZ(s),week_start = 1) == myday))
  m <- calc(sub_s, median,na.rm = T)
  return(m)
}
ver_median <- raster::stack(sapply(c(1:7), ver, s=mystack))

# Print the results
stackapply_median
ver_median