CODEWITHSUNDEEP
rdatemoving-average
foo
foo

Reputation: 43

Moving mean within certain window while preserving Date, R

I have a large rasterstack of 50 layers and I want to calculate the mean and max withing a moving window of 10. What I want for example is to calculate the mean for layers 1:10 and for 11:20 and so on (so in total 5 mean rasters). At the same time, I want to keep the dates of my layers based on the function that I use (e.g. the date of the rasters with the max value or the the mean date).

So far I have tried the following but is very slow. Can anyone hep me to do this more efficiently?

Considering s_all as my raster stack: 

for(i in 1:5){
  t[[i]]<-calc(
            s_all[[((i-1)*10 + 1):((i-1)*10 + 10)]], 
            fun = mean, na.rm = T)
  t@z$Date[[i]]<-mean.Date(as.Date(
                              c(s_all@z$Date[[((i-1)*10 + 1)]],
                              s_all@z$Date[[((i-1)*10 + 10)]])))
}

EDIT Sample data 

r <- raster(ncol=10, nrow=10, vals=1:100)
s_all <- stack(replicate(50, r))
d<-sample(seq(as.Date('1999/01/01'), as.Date('2000/01/01'), by="day"), 50)
s_all<-setZ(s_all,d[],"Date")

Upvotes: 1

Views: 278

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47091

You say you want a "moving average", but that is not what you describe (a moving average would return the same number of layers, with locally smoothed values). What you describe is an aggregation of layers with a step of 10. You can get that with stackApply or aggregate.

I changed your example data to have some variation in the values, and I do not randomize the dates.

library(raster)
r <- raster(ncol=10, nrow=10)
s <- stack(lapply(1:50, function(i) setValues(r, i)))
d <- seq(as.Date('1999/01/01'), as.Date('2000/01/01'), by="day")[1:50]
s <- setZ(s, d, "Date")

For stackApply you need an index:

idx <- rep(1:5, each=10)
ss <- stackApply(s, idx, mean)

Given that the number of layers aggregated is constant, you can also use aggregate

sss <- aggregate(s, c(1,1,10), mean)

Now, set a new date, here I use the maximum.

newd <- tapply(getZ(s), idx, max)
newd <- as.Date(newd, origin="1970-01-01")
ss <- setZ(ss, newd, "Date")

ss
#class      : RasterBrick 
#dimensions : 10, 10, 100, 5  (nrow, ncol, ncell, nlayers)
#resolution : 36, 18  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : memory
#names      : index_1, index_2, index_3, index_4, index_5 
#min values :     5.5,    15.5,    25.5,    35.5,    45.5 
#max values :     5.5,    15.5,    25.5,    35.5,    45.5 
#Date       : 1999-01-10, 1999-01-20, 1999-01-30, 1999-02-09, 1999-02-19

Upvotes: 1

Related Questions