# R circular statistics on a raster time series

I have a `raster time-series` with values from 1 to 36 representing 10 days (`dekads`) timesteps and NA values. I need to calculate the `average` and the `standard deviation`, taking into account that dekads are `circular` (36 is last dekad of the year, 1 is the first ones).

Here below is the code for the average (it seems working) but I am not able to write a code for the standard deviation.

Can you check if the code for the average calculation produces the right outputs and help me in writing the ones for the standard deviation.

``````require(raster)
r <- raster(ncol=50, nrow=50)
s <- stack(lapply(1:length(r), function(x) setValues(r, values = sample(x = c(0:36),size = ncell(r), replace = T))))

# 0 is set to NA

s[s==0]<-NA

#  function
conv <- 2*pi/36  #SET HERE THE time steps

fun1 <- function(m,na.rm=T){
x1 = Arg(mean(exp(conv*(m-1)*1i),na.rm=na.rm))
x2 = x1/conv
x3 = (x2 + 36) %% 36
return(x3)
}

# MEAN calculation
s_avg <- calc(s,fun1)
s_avg <- round(s_avg,0)
s_avg <- s_avg + 1
s_avg [s_avg==37] <- 1
s_avg

plot(s_avg)
``````
• Have you tested your `fun1` function? For example, `fun1(c(1,2))` returns 0.5, and `fun1(c(37,38))` also returns 0.5. Is that expected and correct? If not, what is the correct answer given some sample inputs? – Spacedman Dec 29 '18 at 13:33
• Hi, you have got the point, 1 and 2 are the first two dekads of year n while 37 and 38 the first two dekads of year n+1 (a year has 36 dekads), so the result is correct. – Gianca Jan 2 at 14:28