[GRASS-SVN] r63407 - grass-addons/grass7/imagery/i.theilsen

[svn_grass at osgeo.org]

Author: ychemin
Date: 2014-12-07 01:19:11 -0800 (Sun, 07 Dec 2014)
New Revision: 63407

Added:
   grass-addons/grass7/imagery/i.theilsen/mk.c
Modified:
   grass-addons/grass7/imagery/i.theilsen/main.c
Log:
Added Mann-Kendall test

Modified: grass-addons/grass7/imagery/i.theilsen/main.c
===================================================================
--- grass-addons/grass7/imagery/i.theilsen/main.c	2014-12-07 05:45:06 UTC (rev 63406)
+++ grass-addons/grass7/imagery/i.theilsen/main.c	2014-12-07 09:19:11 UTC (rev 63407)
@@ -27,6 +27,9 @@
 int *cellfd;
 int open_files(void);
 /*-------------------------------------*/
+/*Mann-Kendall test*/
+double mk_test(double *signal, int t);
+/*-------------------------------------*/
 
 int main(int argc, char *argv[])
 {
@@ -47,7 +50,8 @@
     DCELL mk_max=-10000.0;/*Mann-Kendall total max for colour palette */
     DCELL mk_min=100000.0;/*Mann-Kendall total min for colour palette */
     DCELL temp=0.0;/*swapping temp value*/
-    
+    DCELL pvalue=0.0;/*Mann-Kendall trend test p-value*/    
+
     int outfd0, outfd1;
     DCELL *outrast0, *outrast1;
 
@@ -160,10 +164,12 @@
             if (sorted[n0n1/2]>ts_max)
                 ts_max=sorted[n0n1/2];
             /* Mann-Kendall Trend Test */
-            /* Not yet implemented     */
-            mk_min=ts_min;
-            mk_max=ts_max;
-            outrast1[col] = sorted[n0n1/2];
+            pvalue=mk_test(signal,nfiles);
+            if(pvalue<mk_min)
+                mk_min=pvalue;
+            if(pvalue>mk_max)
+                mk_max=pvalue;
+            outrast1[col] = pvalue;
             /*-------------------------*/
         }
         Rast_put_d_row(outfd0, outrast0);

Added: grass-addons/grass7/imagery/i.theilsen/mk.c
===================================================================
--- grass-addons/grass7/imagery/i.theilsen/mk.c	                        (rev 0)
+++ grass-addons/grass7/imagery/i.theilsen/mk.c	2014-12-07 09:19:11 UTC (rev 63407)
@@ -0,0 +1,37 @@
+#include <stdio.h>
+#include <math.h>
+#define ISNEG(X) (!((X) > 0) && ((X) != 0))
+
+/* sign() returns -1 (negative) or +1 (positive) or 0 (zero)*/
+int sign(double signedvalue){
+    if (signedvalue == 0.0) return 0;
+    else if (ISNEG(signedvalue)==1) return -1;
+    else return 1;
+}
+
+/*normalized cumulative distribution function*/
+double normcdf(double x, double mu, double sigma){
+    double y;
+    double l[10]={-1.26551223,1.00002368,0.37409196,0.09678418,-0.18628806,0.27886807,-1.13520398,1.48851587,-0.82215223,0.17087277};
+    double val = -(x-mu)/(sigma*sqrt(2.0));
+    double t = 1.0/(1.0+fabs(val)/2.0);
+    double r = t*exp(-fabs(val)*fabs(val)+l[0]+t*(l[1]+t*(l[2]+t*(l[3]+t*(l[4]+t*(l[5]+t*(l[6]+t*(l[7]+t*(l[8]+t*l[9])))))))));
+    if (x >= 0.) y=r/2.0;
+    else y=(2.0-r)/2.0;
+    if (y>1.0) y = 1.0;
+    return(y);
+}
+
+/*Mann-Kendall test input signal and its length (t)*/
+double mk_test(double *signal, int t){
+    double value=0.0, z=0.0;
+    int i, j;
+    for (i=1;i<t-1;i++)
+        for (j=i+1;j<t;j++) 
+            value += sign(signal[j]-signal[i]); 
+    double variance = (t*(t-1)*(2*t+5))/18.0;
+    double stddev = sqrt(variance);
+    if (value > 0) z=((value-1)/stddev)*(value);
+    else z=(value+1)/stddev;
+    return(2*(1-normcdf(fabs(z),0,1)));
+}