[R-sig-Geo] Mann-kendall

Nuno Sá nunocesardesa at gmail.com
Wed Oct 8 23:23:54 CEST 2014


Hello!

I've been trying for two days to implement the Mann-Kendall test. It is
available in IDRISI and it would be cool that it would be available here.

I believe I made part of it correctly but I am stuck in two steps.

1. Calculating the ties/NA coefficients. (Algorithm problem..)

2. Estimating the p-value on the last step (Lack of statistical knowledge
problem maybe)

So here it goes what I have so far. Bear with me that I have tried to keep
my variables names with a certain "instant" meaning.

---------------------- Problem here, is identifying the coefficients for
the ties/NA -------------------

dummy.stack <- raster(stack.ndvi,layer=1) #required to initiate stack
within loop
mtx <- matrix(0,ncol=nlayers(stack.ndvi),nrow=(nlayers(stack.ndvi)-1))
#this allows the construction of a matrix with tie values signalled to 1

for (i in 1:(nlayers(stack.ndvi)-1)){

  print("----- i ----")
  print(i)
  for (j in (i+1):nlayers(stack.ndvi)){

    print(j)

    dummy.stack <-
stack(dummy.stack,sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)))
    tie.break <-
sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)) == 0
    cc <- maxValue(tie.break)
    if( cc == 1){
      mtx[i,j-1]=1
    }
  }
}

S.stack <- dropLayer(dummy.stack,1) #removing dummy

S.stati <- raster(S.stack,layer=1) #for the somatory - initiating in the
second term (can be substituted by in-built overlay function)

for (i in 2:nlayers(S.stack)){

  S.stati <- S.stati + raster(S.stack,layer=i)
}

-----------------------------------------------------------------------

#Ignoring the existence of ties/na (which I've read is acceptable for large
enough datasets)
#then:

  var.mk.sqr <- sqrt(n*(n-1)*(2*n+5)/18)

  S.c1 <- ((S.stati > 0)*S.stati - 1 )/ var.mk.sqr
  S.c2 <- (S.stati ==0)*S.stati* 0
  S.c3 <- ((S.stati < 0)*S.stati + 1 )/ var.mk.sqr

  Z.combi <- overlay(S.c1,S.c2,S.c3,fun="sum")

#"The P value (probability value p) of the MK statistic S of sample data can
be estimated using the normal cumulative distribution function:" (Yue,2002
and Neeti, 2011)
# So, how to implement the CDF using raster data? Am I forced to do series
expansions? would this not be to slow with raster? Or is there a quicker
way to iterate pixel by pixel?

Sorry for the long email and thank you in advance for any
comments/corrections/recommendations!

Best regards to all,
Nuno Sá


-- 

Nuno César de Sá
+351 91 961 90 37

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list