calc_mean.R

# Calculate mean 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)) #  "Sunday"=1,

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


#### zApply method
z <-
  raster::zApply(mystack,
                 by = lubridate::wday,
                 fun = mean,
                 na.rm = T)


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




# Print the results
stackapply_mean
ver_mean
z