[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