0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.