I want to run a mann-kendall trend analysis with pre-whitening on a stacked raster and then save the resulting tau, sen-slope and p_value. I am currently doing this by looping through every pixel but it is incredibly slow so I am wondering if there is a better approach to this? I am using the modifiedmk package for the trend assessment. I would be ok with any package that applies pre-whitening, and returns a sen slope, a tau value and a p-value though.

This is my current code:


#create a fake raster
r = raster(nrow = 50, ncol = 50)
r[] = rnorm(n = ncell(r))

#make raster stack
stacked = stack(x = c(r, r*2, r))
stacked2 = stack(x = c(r, r*2, r))
stacked3 = stack(x = c(r, r*2, r))

#get a depth of 9
final_stacked = stack(stacked, stacked2, stacked3)

#some empty data to repopulate
coef = raster(nrow = 50, ncol = 50)
p = raster(nrow = 50, ncol = 50)
t = raster(nrow = 50, ncol = 50)

#loop through every pixel stack and runb a mann kendall test
for (r in 1:nrow(coef)){
  for (c in 1:ncol(coef)){

    #one pixel stack
    pixel = data.frame(final_stacked[r, c][1,])
    names(pixel) = 'Values'

    #run mann kendall
    result = pwmk(pixel$Values)

    #save the tau, p-value and sen slope
    tau = result[[7]]
    p_value = result[[4]]
    sens = result[[3]]

    #save the value to new rasters
    coef[r,c] = sens
    p[r,c] = p_value
    t[r,c] = tau


I also tried this which should be faster than looping if it worked:

tsfun <- function(x) { 
  } else { 
    r <- data.frame(pwmk(x))
    names(r) = c("Z", "Sen", "Old_Sen", "p", "s", "var(s)", "Tau")
    a <- r$Sen

raster.mkt <- calc(final_stacked, fun=tsfun, na.rm = TRUE) 

but it returns:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function. Perhaps add '...' or 'na.rm' to the function arguments?

Always check that your function for calc works before using it:

> tsfun(runif(3))
Error in names(r) = c("Z", "Sen", "Old_Sen", "p", "s", "var(s)", "Tau") : 
  'names' attribute [7] must be the same length as the vector [1]

calc tries to trap errors and gives a less informative error message. Running it separately can give an insight.

This error is because you've got a data frame with one column and seven rows and trying to give the column names. You can more easily extract the element you want from the pwmk function:

 pwmk(runif(3))[["Sen's Slope"]]

So rewrite tsfun like this:

> tsfun = function(x,na.rm){
      return(pwmk(x)[["Sen's Slope"]])

and then this works:

> raster.mkt <- calc(final_stacked, fun=tsfun, na.rm = TRUE) 

Problem here might be if you want to get more than just one value out of the pwmk function - then you either have to call this several times or try something else. Here's something else:

First make some test data

r1 = raster(matrix(runif(24),4,6))
r2 = raster(matrix(runif(24),4,6))
r3 = raster(matrix(runif(24),4,6))
z = stack(r1,r2,r3)

Then convert the stack into an array, and apply the pwmk function, converting to vector because pwmk complains about getting a matrix:

pw = apply(as.array(z),c(1,2),function(x){pwmk(as.vector(x))})

That is now an array with the values from pwmk in. Convert to a stack requires a bit of permutation but after trial and error this seems to do it:

s = stack(brick(aperm(pw,c(2,3,1))))

To test, the pwmk value from cell (2,3) is:

> pwmk(c(z[2,3,]))
         Z-Value      Sen's Slope old. Sen's Slope          P-value 
       0.0000000       -0.2621047       -0.2515844        1.0000000 
               S           Var(S)              Tau 
      -1.0000000        1.0000000       -1.0000000 

Which should be element (2,3,) from the stack:

> s[2,3,]
     Z.Value Sen.s.Slope old..Sen.s.Slope P.value  S Var.S. Tau
[1,]       0  -0.2621047       -0.2515844       1 -1      1  -1

So now you have a stack with 7 layers, one from each of the quantities output by pwmk:

> s
class       : RasterStack 
dimensions  : 4, 6, 24, 7  (nrow, ncol, ncell, nlayers)
resolution  : 0.1666667, 0.25  (x, y)
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
names       :    Z.Value, Sen.s.Slope, old..Sen.s.Slope,    P.value,          S,     Var.S.,        Tau 
min values  :  0.0000000,  -0.5586920,       -0.4114719,  1.0000000, -1.0000000,  1.0000000, -1.0000000 
max values  :  0.0000000,   0.5530541,        0.2235587,  1.0000000,  1.0000000,  1.0000000,  1.0000000 

Then doing image(s[[2]]) should be a map of "Sen's Slope".

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.