Calculate median per day of week

calc_median.R

# Calculate median per day of week

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

r <- raster(nrows = 10, ncols = 10, res = 1, 
            xmn = 0, 
            xmx = 10, 
            ymn = 0, 
            ymx = 10)


n=2000

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