New temporal accumulation modules.

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@58651 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/Makefile

@@ -88,6 +88,7 @@ SUBDIRS = \
 	r.rescale.eq \
 	r.ros \
 	r.series \
+	r.series.accumulate \
 	r.series.interp \
 	r.shaded.relief \
 	r.slope.aspect \

raster/r.series.accumulate/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.series.accumulate
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+EXTRA_LIBS = $(OMPLIB)
+EXTRA_CFLAGS = $(OMPCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.series.accumulate/main.c

@@ -0,0 +1,458 @@
+/****************************************************************************
+ * 
+ * MODULE:       r.series.accumulate
+ * AUTHOR(S):    Markus Metz
+ *               Soeren Gebbert
+ *               based on r.series
+ * PURPOSE:      Calculates (accumulated) raster value means, growing degree days
+ *               (GDDs) or Winkler indices from several input maps.
+ * COPYRIGHT:    (C) 2012 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#define METHOD_GDD 1
+#define METHOD_MEAN 2
+#define METHOD_WINKLER 3
+
+struct map_info
+{
+    const char *name;
+    int fd;
+    DCELL *buf;
+};
+
+struct map_info_out
+{
+    const char *name;
+    int fd;
+    void *buf;
+};
+
+int main(int argc, char *argv[])
+{
+    struct GModule *module;
+    struct
+    {
+        struct Option *input, *basemap, *file, *output,
+        *range, *scale, *shift, *lower,
+        *upper, *limits, *method;
+    } parm;
+    struct
+    {
+        struct Flag *nulls, *lazy, *float_output;
+    } flag;
+    int i;
+    int num_inputs, max_inputs;
+    int method;
+    struct map_info *inputs = NULL;
+    struct map_info_out *out = NULL;
+    struct map_info *basemap = NULL;
+    struct map_info *map_lower = NULL;
+    struct map_info *map_upper = NULL;
+    struct History history;
+    int nrows, ncols;
+    int row, col;
+    DCELL lo, hi, tscale, tshift, lower = 10.0, upper = 30.0;
+    DCELL dcell_null;
+    RASTER_MAP_TYPE out_type;
+    int out_size;
+    
+    G_gisinit(argv[0]);
+    
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("series"));
+    G_add_keyword(_("accumulation"));
+    module->description = _("Calculates (accumulated) raster value means, growing degree days (GDDs) or Winkler indices from several input maps.");
+    
+    parm.basemap = G_define_standard_option(G_OPT_R_INPUT);
+    parm.basemap->key = "basemap";
+    parm.basemap->description = _("Existing map to be added to output");
+    parm.basemap->required = NO;
+    
+    parm.input = G_define_standard_option(G_OPT_R_INPUTS);
+    parm.input->required = NO;
+    
+    parm.file = G_define_standard_option(G_OPT_F_INPUT);
+    parm.file->key = "file";
+    parm.file->description = _("Input file with raster map names, one per line");
+    parm.file->required = NO;
+    
+    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+    parm.output->multiple = NO;
+    
+    parm.scale = G_define_option();
+    parm.scale->key = "scale";
+    parm.scale->type = TYPE_DOUBLE;
+    parm.scale->answer = "1.0";
+    parm.scale->required = NO;
+    parm.scale->description = _("Scale factor for input");
+    
+    parm.shift = G_define_option();
+    parm.shift->key = "shift";
+    parm.shift->type = TYPE_DOUBLE;
+    parm.shift->answer = "0.0";
+    parm.shift->required = NO;
+    parm.shift->description = _("Shift factor for input");
+    
+    parm.lower = G_define_standard_option(G_OPT_R_INPUT);
+    parm.lower->key = "lower";
+    parm.lower->required = NO;
+    parm.lower->description = _("The raster map specifying the lower accumulation limit");
+    
+    parm.upper = G_define_standard_option(G_OPT_R_INPUT);
+    parm.upper->key = "upper";
+    parm.upper->required = NO;
+    parm.upper->description = _("The raster map specifying the upper accumulation limit");
+    
+    parm.range = G_define_option();
+    parm.range->key = "range";
+    parm.range->type = TYPE_DOUBLE;
+    parm.range->key_desc = "min,max";
+    parm.range->description = _("Ignore values outside this range");
+    
+    parm.limits = G_define_option();
+    parm.limits->key = "limits";
+    parm.limits->type = TYPE_DOUBLE;
+    parm.limits->key_desc = "lower,upper";
+    parm.limits->answer = "10,30";
+    parm.limits->description = _("Use these limits in case lower and/or upper input maps are not defined");
+    
+    parm.method = G_define_option();
+    parm.method->key = "method";
+    parm.method->type = TYPE_STRING;
+    parm.method->options = "mean,gdd,winkler";
+    parm.method->answer = "mean";
+    parm.method->description = _("This method will be applied to compute the accumulative values from the input maps");
+    
+    flag.nulls = G_define_flag();
+    flag.nulls->key = 'n';
+    flag.nulls->description = _("Propagate NULLs");
+    
+    flag.lazy = G_define_flag();
+    flag.lazy->key = 'z';
+    flag.lazy->description = _("Don't keep files open");
+    
+    flag.float_output = G_define_flag();
+    flag.float_output->key = 'f';
+    flag.float_output->description = _("Create a FCELL map (floating point single precision) as output");
+    
+    if (G_parser(argc, argv))
+        exit(EXIT_FAILURE);
+    
+    lo = -1.0 / 0.0;    /* -inf */
+    hi =  1.0 / 0.0;    /* inf */
+    
+    if (G_strncasecmp(parm.method->answer, "gdd", 3) == 0)
+        method = METHOD_GDD;
+    
+    if (G_strncasecmp(parm.method->answer, "mean", 4) == 0)
+        method = METHOD_MEAN;
+    
+    if (G_strncasecmp(parm.method->answer, "winkler", 7) == 0)
+        method = METHOD_WINKLER;
+    
+    if (parm.range->answer) {
+        lo = atof(parm.range->answers[0]);
+        hi = atof(parm.range->answers[1]);
+    }
+    
+    if (parm.limits->answer) {
+        lower = atof(parm.limits->answers[0]);
+        upper = atof(parm.limits->answers[1]);
+    }
+    
+    if (parm.scale->answer)
+        tscale = atof(parm.scale->answer);
+    else
+        tscale = 1.;
+    
+    if (parm.shift->answer)
+        tshift = atof(parm.shift->answer);
+    else
+        tshift = 0.;
+    
+    if (parm.input->answer && parm.file->answer)
+        G_fatal_error(_("input= and file= are mutually exclusive"));
+    
+    if (!parm.input->answer && !parm.file->answer)
+        G_fatal_error(_("Please specify input= or file="));
+    
+    max_inputs = 0;
+    
+    /* process the input maps from the file */
+    if (parm.file->answer) {
+        FILE *in;
+        
+        in = fopen(parm.file->answer, "r");
+        if (!in)
+            G_fatal_error(_("Unable to open input file <%s>"),
+                          parm.file->answer);
+            
+            num_inputs = 0;
+        
+        for (;;) {
+            char buf[GNAME_MAX];
+            char *name;
+            struct map_info *p;
+            
+            if (!G_getl2(buf, sizeof(buf), in))
+                break;
+            
+            name = G_chop(buf);
+            
+            /* Ignore empty lines */
+            if (!*name)
+                continue;
+            
+            if (num_inputs >= max_inputs) {
+                max_inputs += 100;
+                inputs =
+                G_realloc(inputs, max_inputs * sizeof(struct map_info));
+            }
+            p = &inputs[num_inputs++];
+            
+            p->name = G_store(name);
+            G_verbose_message(_("Reading raster map <%s>..."), p->name);
+            p->buf = Rast_allocate_d_buf();
+            if (!flag.lazy->answer)
+                p->fd = Rast_open_old(p->name, "");
+        }
+        
+        if (num_inputs < 1)
+            G_fatal_error(_("No raster map name found in input file"));
+        
+        fclose(in);
+    }
+    else {
+        for (i = 0; parm.input->answers[i]; i++)
+            ;
+        num_inputs = i;
+        
+        if (num_inputs < 1)
+            G_fatal_error(_("Raster map not found"));
+        
+        inputs = G_malloc(num_inputs * sizeof(struct map_info));
+        
+        for (i = 0; i < num_inputs; i++) {
+            struct map_info *p = &inputs[i];
+            
+            p->name = parm.input->answers[i];
+            G_verbose_message(_("Reading raster map <%s>..."), p->name);
+            p->buf = Rast_allocate_d_buf();
+            if (!flag.lazy->answer)
+                p->fd = Rast_open_old(p->name, "");
+        }
+        max_inputs = num_inputs;
+    }
+    
+    if (parm.basemap->answer) {
+        basemap = G_malloc(1 * sizeof(struct map_info));
+        basemap->name = parm.basemap->answer;
+        G_verbose_message(_("Reading raster map <%s>..."), basemap->name);
+        basemap->buf = Rast_allocate_d_buf();
+        basemap->fd = Rast_open_old(basemap->name, "");
+    }
+    
+    if (parm.lower->answer) {
+        map_lower = G_malloc(1 * sizeof(struct map_info));
+        map_lower->name = parm.lower->answer;
+        G_verbose_message(_("Reading raster map <%s>..."), map_lower->name);
+        map_lower->buf = Rast_allocate_d_buf();
+        map_lower->fd = Rast_open_old(map_lower->name, "");
+    }
+    
+    if (parm.upper->answer) {
+        map_upper = G_malloc(1 * sizeof(struct map_info));
+        map_upper->name = parm.upper->answer;
+        G_verbose_message(_("Reading raster map <%s>..."), map_upper->name);
+        map_upper->buf = Rast_allocate_d_buf();
+        map_upper->fd = Rast_open_old(map_upper->name, "");
+    }
+    
+    out = G_calloc(1, sizeof(struct map_info_out));
+    out->name = parm.output->answer;
+    out_type = flag.float_output->answer ? FCELL_TYPE : DCELL_TYPE;
+    out->buf = Rast_allocate_buf(out_type);
+    out_size = Rast_cell_size(out_type);
+    out->fd = Rast_open_new(out->name, out_type);
+    
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+    
+    Rast_set_d_null_value(&dcell_null, 1);
+    
+    /* process the data */
+    G_verbose_message(_("Percent complete..."));
+    
+    for (row = 0; row < nrows; row++) {
+        G_percent(row, nrows, 2);
+        
+        if (basemap)
+            Rast_get_d_row(basemap->fd, basemap->buf, row);
+        if (map_lower)
+            Rast_get_d_row(map_lower->fd, map_lower->buf, row);
+        if (map_upper)
+            Rast_get_d_row(map_upper->fd, map_upper->buf, row);
+        
+        if (flag.lazy->answer) {
+            /* Open the files only on run time */
+            for (i = 0; i < num_inputs; i++) {
+                inputs[i].fd = Rast_open_old(inputs[i].name, "");
+                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
+                Rast_close(inputs[i].fd);
+            }
+        }
+        else {
+            for (i = 0; i < num_inputs; i++)
+                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
+        }
+        
+        #pragma omp for schedule (static) private (col)
+        for (col = 0; col < ncols; col++) {
+            int null = 0, non_null = 0;
+            DCELL min, max, avg, value;
+            
+            if (map_lower)
+                lower = map_lower->buf[col];
+            if (map_upper)
+                upper = map_upper->buf[col];
+            
+            if (upper <= lower)
+                G_fatal_error(_("'%s'=%f must be > '%s'=%f"), parm.upper->key, upper,
+                              parm.lower->key, lower);
+                
+                min = dcell_null;
+            max = dcell_null;
+            avg = 0;
+            
+            for (i = 0; i < num_inputs; i++) {
+                DCELL v = inputs[i].buf[col];
+                
+                if (Rast_is_d_null_value(&v))
+                    null = 1;
+                else {
+                    v = v * tscale + tshift;
+                    if (parm.range->answer && (v < lo || v > hi)) {
+                        null = 1;
+                    }
+                    else  {
+                        if (method == METHOD_MEAN || method == METHOD_WINKLER)
+                            avg += v;
+                        else if (method == METHOD_GDD) {
+                            if (min > v || Rast_is_d_null_value(&min))
+                                min = v;
+                            if (max < v || Rast_is_d_null_value(&max))
+                                max = v;
+                        }
+                        non_null++;
+                    }
+                }
+            }
+            
+            if (!non_null || (null && flag.nulls->answer)) {
+                if (basemap)
+                    value = basemap->buf[col];
+                else
+                    value = dcell_null;
+            }
+            else {
+                switch(method) {
+                    case METHOD_WINKLER:
+                        
+                        avg /= non_null;
+                        
+                        if (avg < lower)
+                            avg = lower;
+                        else if (avg > upper)
+                            avg = upper;
+                        
+                        value = avg - lower;
+                        
+                        if (value < 0.)
+                            value = 0.;
+                        
+                        break;
+                    case METHOD_GDD:
+                        
+                        if (min < lower)
+                            min = lower;
+                        else if (min > upper)
+                            min = upper;
+                        if (max < lower)
+                            max = lower;
+                        else if (max > upper)
+                            max = upper;
+                        
+                        value = (min + max) / 2. - lower;
+                        
+                        if (value < 0.)
+                            value = 0.;
+                        
+                        break;
+                    default:
+                        
+                        avg /= non_null;
+                        
+                        if (avg < lower)
+                            avg = 0.0;
+                        else if (avg > upper)
+                            avg = upper;
+                        
+                        value = avg;
+                }
+                
+                if (basemap)
+                    value += basemap->buf[col];
+            }
+            Rast_set_d_value((char *)out->buf + col * out_size, value, out_type);
+        }
+        Rast_put_row(out->fd, out->buf, out_type);
+    }
+    
+    G_percent(row, nrows, 2);
+    
+    /* close output map */
+    Rast_close(out->fd);
+    
+    Rast_short_history(out->name, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(out->name, &history);
+    
+    /* Close input maps */
+    if (basemap)
+        Rast_close(basemap->fd);
+    if (map_lower)
+        Rast_close(map_lower->fd);
+    if (map_upper)
+        Rast_close(map_upper->fd);
+    
+    if (!flag.lazy->answer) {
+        for (i = 0; i < num_inputs; i++)
+            Rast_close(inputs[i].fd);
+    }
+    
+    if (method == METHOD_GDD) {
+        struct Colors colors;
+        
+        Rast_init_colors(&colors);
+        Rast_make_colors(&colors, "gdd", 0, 6000);
+        Rast_write_colors(out->name, G_mapset(), &colors);
+    }
+    
+    exit(EXIT_SUCCESS);
+}

raster/r.series.accumulate/r.series.accumulate.html

@@ -0,0 +1,109 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.series.accumulate</em> calculates (accumulated) raster value means, growing degree days (GDDs) or
+Winkler indices's from several input maps.
+<p>
+The formula for calculating the raster value means
+<div class="code"><pre>
+    value = average
+</pre></div>
+The formula for calculating GDDs is
+<div class="code"><pre>
+    value = (max + min) / 2 - lower
+</pre></div>
+The formula for calculating the Winkler index is
+<div class="code"><pre>
+    value = average - lower
+</pre></div>
+with <em>min</em> being the minimum value, <em>max</em>
+the maximum value and <em>average</em> the average value.
+The <em>min</em>, <em>max</em> and <em>average</em> values
+are automatically calculated from the input maps.
+<p>
+The <em>shift</em> and <em>scale</em> values are applied directly to the 
+input values. The <em>lower</em> and <em>upper</em> maps, as well as the <em>range</em>
+options are applied to constrain the accumulation. In case the <em>lower</em> and <em>upper</em>
+maps are not provided the <em>limits</em> option with default values will be applied.
+<p>
+Any <em>min</em>, <em>max</em> and <em>average</em> values
+above the <em>upper</em> values are set to <em>upper</em>, and any
+<em>min</em>, <em>max</em> and <em>average</em> values below the 
+<em>lower</em> values are set to <em>lower</em>, or in case of <em>mean</em>
+computation to 0 (zero). Negative results are set to 0 (zero). 
+<p>
+If an existing map is provided with the <em>basemap</em> option, the
+values of this map are added to the output, thus accumulating means, GDDs or winkler indices's.
+
+<h2>NOTES</h2>
+
+The <em>scale</em> and <em>shift</em> parameters are used to transform 
+input values with
+<div class="code"><pre>
+    new = old * scale + shift
+</pre></div>
+<p>
+With the <em>-n</em> flag, any cell for which any of the 
+corresponding input cells are NULL is automatically set to NULL 
+(NULL propagation) and the accumulated value is not calculated.
+<p>
+Without the <em>-n</em> flag, all non-NULL cells are used for calculation.
+<p>
+If the <em>range=</em> option is given, any values which fall outside 
+that range will be treated as if they were NULL. Note that the range is 
+applied to the scaled and shifted input data. The <em>range</em> 
+parameter can be set to <em>low,high</em> thresholds: 
+values outside of this range are treated as NULL (i.e., they will be 
+ignored by most aggregates, or will cause the result to be NULL if -n 
+is given). The <em>low,high</em> thresholds are floating point, so use 
+<em>-inf</em> or <em>inf</em> for a single threshold (e.g., 
+<em>range=0,inf</em> to ignore negative values, or 
+<em>range=-inf,-200.4</em> to ignore values above -200.4).
+<p>
+The number of input raster maps to be processed is given by the limit of the 
+operating system. For example, both the hard and soft limits are 
+typically 1024. The soft limit can be changed with e.g. 
+<tt>ulimit -n 1500</tt> (UNIX-based operating systems) but not higher 
+than the hard limit. If it is too low, you can as superuser add an 
+entry in
+
+<div class="code"><pre>
+/etc/security/limits.conf
+# &lt;domain&gt;      &lt;type&gt;  &lt;item&gt;         &lt;value&gt;
+your_username  hard    nofile          1500
+</pre></div>
+
+This would raise the hard limit to 1500 file. Be warned that more
+files open need more RAM.
+
+<p>
+Use the <em>file</em> option to analyze large amount of raster maps 
+without hitting open files limit and the size limit of command line 
+arguments. The computation is slower than the <em>input</em> option 
+method. For every sinlge row in the output map(s) all input maps are 
+opened and closed. The amount of RAM will rise linear with the 
+number of specified input maps. The input and file options are 
+mutually exclusive. Input is a text file with a new line separated 
+list of raster map names and optional weights. As separator between 
+the map name and the weight the character | must be used.
+
+<h2>EXAMPLES</h2>
+
+Example with MODIS Land Surface Temperature, transforming values from 
+Kelvin * 50 to degrees Celsius:
+<div class="code"><pre>
+r.series.accumulate in=MOD11A1.Day,MOD11A1.Night,MYD11A1.Day,MYD11A1.Night out=MCD11A1.GDD \
+      scale=0.02 shift=-273.15 limits=10,30
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.series.html">r.series</a>,
+<a href="g.mlist.html">g.mlist</a>,
+<a href="g.region.html">g.region</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Markus Metz and Soeren Gebbert (based on r.series)
+
+<p><i>Last changed: $Date: 2012-07-03 23:47:59 +0200 (Di, 03. Jul 2012) $</i>

raster/r.series.accumulate/test_suite/test.r.series.accumulate.sh

@@ -0,0 +1,49 @@
+#!/bin/sh -e
+# Test r.series.accumulate
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster maps with r.mapcalc
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+export GRASS_OVERWRITE=1
+
+r.mapcalc expr="basemap = if(row() == 3, null(), 10)"
+r.mapcalc expr="lower = 5"
+r.mapcalc expr="upper = 10"
+
+r.mapcalc expr="map_a = rand(0, 15)"
+r.mapcalc expr="map_b = rand(1, 14)"
+r.mapcalc expr="map_c = rand(2, 13)"
+
+r.series.accumulate basemap=basemap input=map_a lower=lower \
+                    output=test_accu_1 upper=upper method=gdd -f --verbose 
+
+r.series.accumulate basemap=basemap input=map_a lower=lower \
+                    output=test_accu_2 upper=upper method=winkler -f --verbose 
+
+r.series.accumulate basemap=basemap input=map_a lower=lower \
+                    output=test_accu_3 upper=upper method=mean --verbose 
+
+r.series.accumulate basemap=basemap input=map_a,map_b,map_c limits=5,10 \
+                    output=test_accu_4 method=gdd -f --verbose 
+
+r.series.accumulate basemap=basemap input=map_a,map_b,map_c lower=lower limits=5,10 \
+                    output=test_accu_5 method=winkler -f --verbose 
+
+r.series.accumulate basemap=basemap input=map_a,map_b,map_c lower=lower limits=5,10 \
+                    output=test_accu_6 range=6,9 method=mean --verbose 
+
+# Test for correct results
+for map in `g.mlist type=rast pattern=test_accu_*` ; do
+    r.out.ascii input=${map} output=${map}.txt dp=2
+done
+
+for i in `ls test_accu_*.txt` ; do
+    diff $i "`basename $i .txt`.ref" >> out.diff
+done
+
+CHAR_NUM=`cat out.diff | wc -c`
+
+# Return as exit status 0 in case no diffs are found
+exit $CHAR_NUM

raster/r.series.accumulate/test_suite/test_accu_1.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 10 15 12 10 12 15 12 15 10 15 10 
+13 15 11 14 10 12 14 12 15 15 14 15 
+* * * * * * * * * * * * 
+15 12 15 15 10 13 15 14 13 10 15 13 
+14 15 11 13 15 10 15 10 10 15 11 15 
+15 15 10 15 13 13 14 10 10 10 10 10 
+14 14 10 10 15 10 10 11 10 10 12 14 
+11 10 13 10 15 11 15 15 14 15 11 10 

raster/r.series.accumulate/test_suite/test_accu_2.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 10 15 12 10 12 15 12 15 10 15 10 
+13 15 11 14 10 12 14 12 15 15 14 15 
+* * * * * * * * * * * * 
+15 12 15 15 10 13 15 14 13 10 15 13 
+14 15 11 13 15 10 15 10 10 15 11 15 
+15 15 10 15 13 13 14 10 10 10 10 10 
+14 14 10 10 15 10 10 11 10 10 12 14 
+11 10 13 10 15 11 15 15 14 15 11 10 

raster/r.series.accumulate/test_suite/test_accu_3.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 10 20 17 10 17 20 17 20 10 20 10 
+18 20 16 19 10 17 19 17 20 20 19 20 
+* * * * * * * * * * * * 
+20 17 20 20 10 18 20 19 18 10 20 18 
+19 20 16 18 20 10 20 15 10 20 16 20 
+20 20 10 20 18 18 19 10 15 10 15 10 
+19 19 10 10 20 10 10 16 10 15 17 19 
+16 10 18 15 20 16 20 20 19 20 16 10 

raster/r.series.accumulate/test_suite/test_accu_4.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 12.5 12.5 11 10 11 14.5 11 12.5 12.5 12.5 12.5 
+12.5 12.5 13 14 11 13.5 12.5 11 12.5 14 14.5 14 
+* * * * * * * * * * * * 
+12.5 13 13.5 12.5 12.5 14 12.5 12 13 11.5 13 11.5 
+12.5 12.5 10.5 12.5 14.5 10 15 12.5 11 14 11 14.5 
+12.5 13 12 12.5 13 14 12.5 12.5 12.5 12 11.5 12.5 
+12 12.5 12.5 10 14 10 11 10.5 12.5 12.5 12.5 12 
+12.5 12 11.5 12.5 12.5 12.5 12.5 13 12.5 12.5 13 10.5 

raster/r.series.accumulate/test_suite/test_accu_5.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 12.67 12.33 10.33 10 10.33 15 10.33 14.67 13 11.67 13 
+13 11.67 13 14.33 10 13 12 10 12.67 15 15 14.67 
+* * * * * * * * * * * * 
+15 12.67 15 12 11.67 15 12.33 12.33 13.33 10.67 15 11 
+12.33 12.33 10 13 15 10 15 13.67 10 15 10 15 
+11 13.67 10.67 12.33 13.33 14 13.67 10.33 12.67 10 11.67 10 
+10 12.33 12 10 15 10 10.67 10.33 11.67 14 12 11 
+13 10 10 13 12.67 12.67 10.33 13.67 12.67 12.33 13 10 

raster/r.series.accumulate/test_suite/test_accu_6.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 120
+west: 0
+rows: 8
+cols: 12
+10 10 10 17 10 16.5 19 16.5 10 10 17 10 
+18 10 17 18.5 17 17 19 17 10 18.5 19 18 
+* * * * * * * * * * * * 
+10 16.5 17 19 17 18 10 18.5 17 18 16 18 
+19 16 16 18 19 10 10 19 16.5 18 16.5 19 
+16 17.5 19 17 17 18 19 10 16 19 17.5 10 
+19 17.33 18 10 18 10 17 16 10 10 17 19 
+18 19 18 16 16 17.67 10 16.5 19 17 16.5 16 

temporal/Makefile

@@ -12,6 +12,8 @@ SUBDIRS = \
 	t.sample \
 	t.register \
 	t.unregister \
+	t.rast.accumulate \
+	t.rast.accdetect \
 	t.rast.aggregate \
 	t.rast.aggregate.ds \
 	t.rast.colors \

temporal/t.rast.accdetect/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.accdetect
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)

temporal/t.rast.accdetect/t.rast.accdetect.html

@@ -0,0 +1,63 @@
+<h2>DESCRIPTION</h2>
+
+<b>t.rast.accdetect</b> is designed to detect accumulation pattern in temporally accumulated
+space time raster datasets created by <a href="t.rast.accumulate.html">t.rast.accumulate</a>.
+
+
+This module expects a space time raster dataset as input that is the result 
+of a <a href="t.rast.accumulate.html">t.rast.accumulate</a> run.
+<p>
+The <b>start</b> time and the <b>end</b> time of the pattern detection process must be set, 
+eg. <b>start="2000-03-01" end="2011-01-01"</b>. The <b>start</b> and <b>end</b> time do not need to be the same
+as for the accumulation run that produced the input space time raster dataset. 
+In addition a <b>cycle</b>, eg. "8 months", can be specified, that defines
+after which time interval the accumulation pattern detection process restarts. The <b>offset</b> option specifies the
+time between two cycles that should be skipped, eg. "4 months". Please make sure that the <b>cycle</b> and
+<b>offset</b> options are same as in the accumulation process that produces the input space time raster dataset, otherwise the accumulation pattern detection will produce wrong results.
+<p>
+The <b>minimum</b> and <b>maximum</b> values of the pattern detection 
+process can be set, either by using space time raster datasets or 
+by using fixed values for all raster cells and time steps. 
+<p>
+Using space time raster datasets allow to specify minimum and maximum values for each raster cell and each time step. For example,
+we want to detect the germination (minimum value) and harvesting (maximum value) dates for different crops in Germany using the growing-degree-day (GDD) method for several years. Different crops may grow in different raster cells and change with time because of crop rotation. 
+Hence we need to specify different GDD germination/harvesting (minimum/maximum) values for different raster cells and different years.
+<p>
+The raster maps that specifies the
+minimum and maximum values of the actual granule will be detected using the following temporal relations: 
+equals, during, overlaps, overlapped and contains. 
+First all maps with equal time stamps to the current granule of the input STRDS will be detected, 
+the first minimum map and the first maximum map that were found are used as range definitions. 
+If no equal maps are found then maps with a temporal during relation are detected, 
+then maps that temporally overlap the actual granules, until
+maps are detected that have a temporal contain relation. 
+If no maps are found or minimum/maximum STRDS are not set, then the <b>range</b> option is used, eg. <b>range=480,730</b>.
+<p>
+The <b>base</b> name of of the generated maps must always be set.
+<p>
+This module produces two output space time raster datasets. The <b>occurrence</b> output STRDS stores the time in days
+from the begin of a cycle for each raster cell and time step that has a value within the minimum and maximum definition.
+These values can be used to compute the duration of the recognized accumulation pattern. The <b>indicator</b> output STRDS
+uses three values, that can be set using the <b>staend</b> option, to mark raster cells with integer values that indicate the start, the intermediate state and the end of a accumulation pattern. As default specifies the value 1 the start, the value 2 the intermediate state and the value 3 the end of the accumulation pattern in a cycle.
+
+<h2>EXAMPLE</h2>
+
+Please have a look at the <a href="t.rast.accumulate.html">t.rast.accumulate</a> example.
+
+<h2>SEE ALSO</h2>
+
+<b>
+<a href="t.rast.accumulate.html">t.rast.accumulate</a><br>
+<a href="t.rast.aggregate.html">t.rast.aggregate</a><br>
+<a href="t.rast.mapcalc.html">t.rast.mapcalc</a><br>
+<a href="t.info.html">t.info</a><br>
+<a href="r.series.accumulate.html">r.series.accumulate</a><br>
+<a href="g.region.html">g.region</a><br>
+</b>
+
+
+<h2>AUTHOR</h2>
+
+S&ouml;ren Gebbert
+
+<p><i>Last changed: $Date: 2012-08-29 14:55:21 +0200 (Mi, 29. Aug 2012) $</i>

temporal/t.rast.accdetect/t.rast.accdetect.py

@@ -0,0 +1,551 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:       t.rast.accdetect
+# AUTHOR(S):    Soeren Gebbert
+#
+# PURPOSE:      Detect accumulation pattern in temporally accumulated space time raster datasets created by t.rast.accumulate.
+# COPYRIGHT:    (C) 2013 - 2014 by the GRASS Development Team
+#
+#               This program is free software under the GNU General Public
+#               License (version 2). Read the file COPYING that comes with GRASS
+#               for details.
+#
+#############################################################################
+
+#%module
+#% description: Detect accumulation pattern in temporally accumulated space time raster datasets created by t.rast.accumulate.
+#% keywords: temporal
+#% keywords: accumulation
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#% key: minimum
+#% description: Input space time raster dataset that specifies the minimum values to detect the accumulation pattern
+#% required: no
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#% key: maximum
+#% description: Input space time raster dataset that specifies the maximum values to detect the accumulation pattern
+#% required: no
+#%end
+
+#%option G_OPT_STRDS_OUTPUT
+#% key: occurrence
+#% description: The output space time raster dataset that stores the occurrence of the the accumulation pattern using the provided data range
+#% required: yes
+#%end
+
+#%option G_OPT_STRDS_OUTPUT
+#% key: indicator
+#% description: The output space time raster dataset that stores the indication of the start, intermediate and end of the specified data range
+#% required: no
+#%end
+
+#%option
+#% key: start
+#% type: string
+#% description: The temporal starting point to begin the accumulation, eg '2001-01-01'
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: stop
+#% type: string
+#% description: The temporal date to stop the accumulation, eg '2009-01-01'
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: cycle
+#% type: string
+#% description: The temporal cycle to restart the accumulation, eg '12 months'
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: offset
+#% type: string
+#% description: The temporal offset to the begin of the next cycle, eg '6 months'
+#% required: no
+#% multiple: no
+#%end
+
+#%option G_OPT_R_BASE
+#%end
+
+#%option
+#% key: range
+#% type: double
+#% key_desc: min,max
+#% description: The minimum and maximum value of the occurrence of accumulated values, these values will be used if the min/max space time raster datasets are not specified
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: staend
+#% type: integer
+#% key_desc: start,intermediate,end
+#% description: The user defined values that indicate start, intermediate and end status in the indicator output space time raster dataset
+#% answer: 1,2,3
+#% required: no
+#% multiple: no
+#%end
+
+#%flag
+#% key: n
+#% description: Register empty maps in the output space time raster datasets, otherwise they will be deleted
+#%end
+
+#%flag
+#% key: r
+#% description: Reverse time direction in cyclic accumulation
+#%end
+
+import grass.script as grass
+import grass.temporal as tgis
+
+
+############################################################################
+
+range_relations = ["EQUALS", "DURING", "OVERLAPS", "OVERLAPPING", "CONTAINS"]
+
+def main():
+    # Get the options
+    input = options["input"]
+    start = options["start"]
+    stop = options["stop"]
+    base = options["base"]
+    cycle = options["cycle"]
+    offset = options["offset"]
+    minimum = options["minimum"]
+    maximum = options["maximum"]
+    occurrence = options["occurrence"]
+    range = options["range"]
+    indicator = options["indicator"]
+    staend = options["staend"]
+    register_null = flags["n"]
+    reverse = flags["r"]
+
+    grass.set_raise_on_error(True)
+
+    # Make sure the temporal database exists
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    mapset = tgis.get_current_mapset()
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        id = input + "@" + mapset
+
+    input_strds = tgis.SpaceTimeRasterDataset(id)
+
+    if input_strds.is_in_db() == False:
+        dbif.close()
+        grass.fatal(_("Space time %s dataset <%s> not found") % (
+            input_strds.get_output_map_instance(None).get_type(), id))
+
+    input_strds.select(dbif)
+    dummy = input_strds.get_new_map_instance(None)
+
+    # The occurrence space time raster dataset
+    if occurrence:
+        if not minimum or not maximum:
+            if not range:
+                dbif.close()
+                grass.fatal(_("You need to set the range to compute the occurrence"
+                              " space time raster dataset"))
+
+        if occurrence.find("@") >= 0:
+            occurrence_id = occurrence
+        else:
+            occurrence_id = occurrence + "@" + mapset
+
+        occurrence_strds = tgis.SpaceTimeRasterDataset(occurrence_id)
+        if occurrence_strds.is_in_db(dbif):
+            if not grass.overwrite():
+                dbif.close()
+                grass.fatal(_("Space time raster dataset <%s> is already in the "
+                              "database, use overwrite flag to overwrite") % occurrence_id)
+
+    # The indicator space time raster dataset
+    if indicator:
+        if not occurrence:
+            dbif.close()
+            grass.fatal(_("You need to set the occurrence to compute the indicator"
+                          " space time raster dataset"))
+        if not staend:
+            dbif.close()
+            grass.fatal(_("You need to set the staend options to compute the indicator"
+                          " space time raster dataset"))
+        if indicator.find("@") >= 0:
+            indicator = indicator
+        else:
+            indicator_id = indicator + "@" + mapset
+
+        indicator_strds = tgis.SpaceTimeRasterDataset(indicator_id)
+        if indicator_strds.is_in_db(dbif):
+            if not grass.overwrite():
+                dbif.close()
+                grass.fatal(_("Space time raster dataset <%s> is already in the "
+                              "database, use overwrite flag to overwrite") % indicator_id)
+        staend = staend.split(",")
+        indicator_start = int(staend[0])
+        indicator_mid = int(staend[1])
+        indicator_end = int(staend[2])
+
+    # The minimum threshold space time raster dataset
+    minimum_strds = None
+    if minimum:
+        if minimum.find("@") >= 0:
+            minimum_id = minimum
+        else:
+            minimum_id = minimum + "@" + mapset
+
+        minimum_strds = tgis.SpaceTimeRasterDataset(minimum_id)
+        if minimum_strds.is_in_db() == False:
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> not found") % (minimum_strds.get_id()))
+
+        if minimum_strds.get_temporal_type() != input_strds.get_temporal_type():
+            dbif.close()
+            grass.fatal(_("Temporal type of input strds and minimum strds must be equal"))
+
+        minimum_strds.select(dbif)
+
+    # The maximum threshold space time raster dataset
+    maximum_strds = None
+    if maximum:
+        if maximum.find("@") >= 0:
+            maximum_id = maximum
+        else:
+            maximum_id = maximum + "@" + mapset
+
+        maximum_strds = tgis.SpaceTimeRasterDataset(maximum_id)
+        if maximum_strds.is_in_db() == False:
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> not found") % (maximum_strds.get_id()))
+
+        if maximum_strds.get_temporal_type() != input_strds.get_temporal_type():
+            dbif.close()
+            grass.fatal(_("Temporal type of input strds and maximum strds must be equal"))
+
+        maximum_strds.select(dbif)
+
+    input_strds_start, input_strds_end = input_strds.get_temporal_extent_as_tuple()
+
+    if input_strds.is_time_absolute():
+        start = tgis.string_to_datetime(start)
+        if stop:
+            stop = tgis.string_to_datetime(stop)
+        else:
+            stop = input_strds_end
+    else:
+        start = int(start)
+        if stop:
+            stop = int(stop)
+        else:
+            stop = input_strds_end
+
+    if input_strds.is_time_absolute():
+        end = tgis.increment_datetime_by_string(start, cycle)
+    else:
+        end = start + cycle
+
+    count = 1
+    indi_count = 1
+    occurrence_maps = {}
+    indicator_maps = {}
+
+    while input_strds_end > start and stop > start:
+
+        # Make sure that the cyclic computation will stop at the correct time
+        if stop and end > stop:
+            end = stop
+
+        where = "start_time >= \'%s\' AND start_time < \'%s\'"%(str(start),
+                                                                str(end))
+        input_maps = input_strds.get_registered_maps_as_objects(where=where,
+                                                                dbif=dbif)
+
+        print len(input_maps)
+
+        input_topo = tgis.SpatioTemporalTopologyBuilder()
+        input_topo.build(input_maps, input_maps)
+
+        if len(input_maps) == 0:
+            continue
+
+        grass.message(_("Processing cycle %s - %s"%(str(start), str(end))))
+
+        count = compute_occurrence(occurrence_maps, input_strds, input_maps,
+                                   start, base, count, mapset, where, reverse,
+                                   range, minimum_strds, maximum_strds, dbif)
+
+        # Indicator computation is based on the occurrence so we need to start it after
+        # the occurrence cycle
+        if indicator:
+            num_maps = len(input_maps)
+            for i in xrange(num_maps):
+                if reverse:
+                    map = input_maps[num_maps - i - 1]
+                else:
+                    map = input_maps[i]
+
+                indicator_map_name = "%s_indicator_%i" % (base, indi_count)
+                indicator_map_id = dummy.build_id(indicator_map_name, mapset)
+                indicator_map = input_strds.get_new_map_instance(indicator_map_id)
+
+                # Check if new map is in the temporal database
+                if indicator_map.is_in_db(dbif):
+                    if grass.overwrite():
+                        # Remove the existing temporal database entry
+                        indicator_map.delete(dbif)
+                        indicator_map = input_strds.get_new_map_instance(indicator_map_id)
+                    else:
+                        grass.fatal(_("Map <%s> is already registered in the temporal"
+                                     " database, use overwrite flag to overwrite.") %
+                                    (indicator_map.get_map_id()))
+
+                curr_map = occurrence_maps[map.get_id()].get_name()
+
+                # Reverse time
+                if reverse:
+                    if i ==  0:
+                        prev_map = curr_map
+                        subexpr1 = "null()"
+                        subexpr3 = "%i"%(indicator_start)
+                    elif i > 0 and i < num_maps - 1:
+                        prev_map = occurrence_maps[map.next().get_id()].get_name()
+                        next_map = occurrence_maps[map.prev().get_id()].get_name()
+                        # In case the previous map is null() set null() or the start indicator
+                        subexpr1 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_start)
+                        # In case the previous map was not null() if the current map is null() set null()
+                        # if the current map is not null() and the next map is not null() set
+                        # intermediate indicator, if the next map is null set the end indicator
+                        subexpr2 = "if(isnull(%s), %i, %i)"%(next_map, indicator_end, indicator_mid)
+                        subexpr3 = "if(isnull(%s), null(), %s)"%(curr_map, subexpr2)
+                        expression = "%s = if(isnull(%s), %s, %s)"%(indicator_map_name,
+                                                                    prev_map, subexpr1,
+                                                                    subexpr3)
+                    else:
+                        prev_map = occurrence_maps[map.next().get_id()].get_name()
+                        subexpr1 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_start)
+                        subexpr3 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_mid)
+                else:
+                    if i == 0:
+                        prev_map = curr_map
+                        subexpr1 = "null()"
+                        subexpr3 = "%i"%(indicator_start)
+                    elif i > 0 and i < num_maps - 1:
+                        prev_map = occurrence_maps[map.prev().get_id()].get_name()
+                        next_map = occurrence_maps[map.next().get_id()].get_name()
+                        # In case the previous map is null() set null() or the start indicator
+                        subexpr1 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_start)
+                        # In case the previous map was not null() if the current map is null() set null()
+                        # if the current map is not null() and the next map is not null() set
+                        # intermediate indicator, if the next map is null set the end indicator
+                        subexpr2 = "if(isnull(%s), %i, %i)"%(next_map, indicator_end, indicator_mid)
+                        subexpr3 = "if(isnull(%s), null(), %s)"%(curr_map, subexpr2)
+                        expression = "%s = if(isnull(%s), %s, %s)"%(indicator_map_name,
+                                                                    prev_map, subexpr1,
+                                                                    subexpr3)
+                    else:
+                        prev_map = occurrence_maps[map.prev().get_id()].get_name()
+                        subexpr1 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_start)
+                        subexpr3 = "if(isnull(%s), null(), %i)"%(curr_map, indicator_mid)
+
+                expression = "%s = if(isnull(%s), %s, %s)"%(indicator_map_name,
+                                                            prev_map, subexpr1,
+                                                            subexpr3)
+                print expression
+                grass.mapcalc(expression, overwrite=True)
+
+                map_start, map_end = map.get_temporal_extent_as_tuple()
+
+                if map.is_time_absolute():
+                    indicator_map.set_absolute_time(map_start, map_end)
+                else:
+                    indicator_map.set_relative_time(map_start, map_end,
+                                                 map.get_relative_time_unit())
+
+                indicator_maps[map.get_id()] = indicator_map
+                indi_count += 1
+
+        # Increment the cycle
+        start = end
+        if input_strds.is_time_absolute():
+            start = end
+            if offset:
+                start = tgis.increment_datetime_by_string(end, offset)
+            end = tgis.increment_datetime_by_string(start, cycle)
+        else:
+            if offset:
+                start = end + offset
+            end = start + cycle
+
+    empty_maps = []
+
+    create_strds_register_maps(input_strds, occurrence_strds, occurrence_maps,
+                               register_null, empty_maps, dbif)
+
+    if indicator:
+        create_strds_register_maps(input_strds, indicator_strds, indicator_maps,
+                                   register_null, empty_maps, dbif)
+
+    dbif.close()
+
+    # Remove empty maps
+    if len(empty_maps) > 0:
+        for map in empty_maps:
+            grass.run_command("g.remove", rast=map.get_name(), quiet=True)
+
+############################################################################
+
+def create_strds_register_maps(in_strds, out_strds, out_maps, register_null,
+                    empty_maps, dbif):
+
+    out_id = out_strds.get_id()
+
+    if out_strds.is_in_db(dbif):
+        if grass.overwrite():
+            out_strds.delete(dbif)
+            out_strds = in_strds.get_new_instance(out_id)
+
+    temporal_type, semantic_type, title, description = in_strds.get_initial_values()
+    out_strds.set_initial_values(temporal_type, semantic_type, title,
+                                    description)
+    out_strds.insert(dbif)
+
+    # Register the maps in the database
+    count = 0
+    for map in out_maps.values():
+        count += 1
+        if count%10 == 0:
+            grass.percent(count, len(out_maps), 1)
+        # Read the raster map data
+        map.load()
+        # In case of a empty map continue, do not register empty maps
+        if not register_null:
+            if map.metadata.get_min() is None and \
+                map.metadata.get_max() is None:
+                empty_maps.append(map)
+                continue
+
+        # Insert map in temporal database
+        map.insert(dbif)
+        out_strds.register_map(map, dbif)
+
+    out_strds.update_from_registered_maps(dbif)
+    grass.percent(1, 1, 1)
+
+############################################################################
+
+def compute_occurrence(occurrence_maps, input_strds, input_maps, start, base,
+               count, mapset, where, reverse, range, minimum_strds,
+               maximum_strds, dbif):
+
+    if minimum_strds:
+        input_maps_minimum = input_strds.get_registered_maps_as_objects(where=where,
+                                                                      dbif=dbif)
+        minimum_maps = minimum_strds.get_registered_maps_as_objects(dbif=dbif)
+        minimum_topo = tgis.SpatioTemporalTopologyBuilder()
+        minimum_topo.build(input_maps_minimum, minimum_maps)
+
+    if maximum_strds:
+        input_maps_maximum = input_strds.get_registered_maps_as_objects(where=where,
+                                                                      dbif=dbif)
+        maximum_maps = maximum_strds.get_registered_maps_as_objects(dbif=dbif)
+        maximum_topo = tgis.SpatioTemporalTopologyBuilder()
+        maximum_topo.build(input_maps_maximum, maximum_maps)
+
+    # Aggregate
+    num_maps = len(input_maps)
+    for i in xrange(num_maps):
+        if reverse:
+            map = input_maps[num_maps - i - 1]
+        else:
+            map = input_maps[i]
+
+        # Compute the days since start
+        input_start, input_end = map.get_temporal_extent_as_tuple()
+
+        td = input_start - start
+        if map.is_time_absolute():
+            days = tgis.time_delta_to_relative_time(td)
+        else:
+            days = td
+
+        occurrence_map_name = "%s_%i" % (base, count)
+        occurrence_map_id = map.build_id(occurrence_map_name, mapset)
+        occurrence_map = input_strds.get_new_map_instance(occurrence_map_id)
+
+        # Check if new map is in the temporal database
+        if occurrence_map.is_in_db(dbif):
+            if grass.overwrite():
+                # Remove the existing temporal database entry
+                occurrence_map.delete(dbif)
+                occurrence_map = input_strds.get_new_map_instance(occurrence_map_id)
+            else:
+                grass.fatal(_("Map <%s> is already registered in the temporal"
+                             " database, use overwrite flag to overwrite.") %
+                            (occurrence_map.get_map_id()))
+
+        range_vals = range.split(",")
+        min = range_vals[0]
+        max = range_vals[1]
+
+        if minimum_strds:
+            relations = input_maps_minimum[i].get_temporal_relations()
+            for relation in range_relations:
+                if relation in relations:
+                    min = str(relations[relation][0].get_id())
+                    break
+
+        if maximum_strds:
+            relations = input_maps_maximum[i].get_temporal_relations()
+            for relation in range_relations:
+                if relation in relations:
+                    max = str(relations[relation][0].get_id())
+                    break
+
+        expression = "%s = if(%s > %s && %s < %s, %s, null())"%(occurrence_map_name,
+                                                                map.get_name(),
+                                                                min, map.get_name(),
+                                                                max, days)
+        print expression
+        grass.mapcalc(expression, overwrite=True)
+
+        map_start, map_end = map.get_temporal_extent_as_tuple()
+
+        if map.is_time_absolute():
+            occurrence_map.set_absolute_time(map_start, map_end)
+        else:
+            occurrence_map.set_relative_time(map_start, map_end,
+                                         map.get_relative_time_unit())
+
+        # Store the new maps
+        occurrence_maps[map.get_id()] = occurrence_map
+
+        count += 1
+
+    return count
+
+############################################################################
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()

temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.reverse.sh

@@ -0,0 +1,43 @@
+#!/bin/sh -e
+# Space time raster dataset neighborhood operations
+# We need to set a specific region in the
+# @preprocess step of this test.
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 -p
+
+export GRASS_OVERWRITE=1
+
+# Generate data
+r.mapcalc expr="temp_6 = 5"
+r.mapcalc expr="temp_5 = 10"
+r.mapcalc expr="temp_4 = 15"
+r.mapcalc expr="temp_3 = 20"
+r.mapcalc expr="temp_2 = 25"
+r.mapcalc expr="temp_1 = 30"
+
+t.create type=strds temporaltype=absolute output=temp_abs1 title="A test" descr="A test"
+t.register -i type=rast input=temp_abs1 maps=temp_1,temp_2,temp_3,temp_4,temp_5,temp_6 \
+    start="2001-01-01" increment="2 months"
+
+# The first @test
+
+t.rast.accumulate -r input=temp_abs1 output=temp_accumulation base=temp_acc \
+    limits=10,25 start="2001-01-01" gran="2 months" cycle="12 months"
+
+t.rast.accdetect -r input=temp_accumulation occurrence=temp_occ base=temp_occ \
+    range=20,80 start="2001-01-01" cycle="12 months" staend=1,2,3 indi=temp_indi
+
+t.rast.list temp_accumulation col=name,start_time,min,max > test_2_temp_accumulation.txt
+t.rast.list temp_occ col=name,start_time,min,max  > test_2_temp_occ.txt
+t.rast.list temp_indi col=name,start_time,min,max > test_2_temp_indi.txt
+
+t.remove -rf type=strds input=temp_abs1,temp_accumulation,temp_indi
+
+for i in `ls test_2_*.txt` ; do
+    diff $i "`basename $i .txt`.ref" >> out.diff
+done
+
+CHAR_NUM=`cat out.diff | wc -c`
+
+# Return as exit status 0 in case no diffs are found
+exit $CHAR_NUM

temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.sh

@@ -0,0 +1,63 @@
+#!/bin/sh -e
+# Space time raster dataset neighborhood operations
+# We need to set a specific region in the
+# @preprocess step of this test.
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 -p
+
+export GRASS_OVERWRITE=1
+
+# Generate data
+r.mapcalc expr="temp_1 = 5"
+r.mapcalc expr="temp_2 = 10"
+r.mapcalc expr="temp_3 = 15"
+r.mapcalc expr="temp_4 = 20"
+r.mapcalc expr="temp_5 = 25"
+r.mapcalc expr="temp_6 = 30"
+
+t.create type=strds temporaltype=absolute output=temp_abs1 title="A test" descr="A test"
+t.register -i type=rast input=temp_abs1 maps=temp_1,temp_2,temp_3,temp_4,temp_5,temp_6 \
+    start="2001-01-01" increment="2 months"
+
+# The first @test
+
+t.rast.accumulate input=temp_abs1 output=temp_accumulation base=temp_acc \
+    limits=10,25 start="2001-01-01" gran="2 months" cycle="12 months"
+
+t.rast.accdetect input=temp_accumulation occurrence=temp_occ base=temp_occ \
+    range=20,80 start="2001-01-01" cycle="12 months"
+
+# Check the registered maps metadata
+t.rast.list temp_accumulation col=name,start_time,min,max > test_1_temp_accumulation.txt
+t.rast.list temp_occ col=name,start_time,min,max          > test_1_temp_occ_a.txt
+
+# Leets test the minimum and maximum STRDS implementation
+
+r.mapcalc expr="minimum = 18"
+r.mapcalc expr="maximum = 78"
+
+# We use different
+t.create type=strds temporaltype=absolute output=minimum title="minimum limit" descr="minimum limit"
+t.register -i type=rast input=minimum maps=minimum start="2001-01-01" increment="8 months"
+
+t.create type=strds temporaltype=absolute output=maximum title="maximum limit" descr="maximum limit"
+t.register -i type=rast input=maximum maps=maximum start="2001-01-01" increment="10 months"
+
+t.rast.accdetect input=temp_accumulation occurrence=temp_occ base=temp_occ \
+    range=20,80 start="2001-01-01" cycle="12 months" min=minimum \
+    max=maximum staend=1,2,3 indi=temp_indi
+
+# Check the registered maps metadata
+t.rast.list temp_occ col=name,start_time,min,max          > test_1_temp_occ_b.txt
+t.rast.list temp_indi col=name,start_time,min,max         > test_1_temp_indi.txt
+
+t.remove -rf type=strds input=temp_abs1,temp_accumulation,temp_indi,minimum,maximum
+
+for i in `ls test_1_*.txt` ; do
+    diff $i "`basename $i .txt`.ref" >> out.diff
+done
+
+CHAR_NUM=`cat out.diff | wc -c`
+
+# Return as exit status 0 in case no diffs are found
+exit $CHAR_NUM

temporal/t.rast.accdetect/test_suite/test_1_temp_accumulation.ref

@@ -0,0 +1,6 @@
+temp_acc_1	2001-01-01 00:00:00	0.0	0.0
+temp_acc_2	2001-03-01 00:00:00	10.0	10.0
+temp_acc_3	2001-05-01 00:00:00	25.0	25.0
+temp_acc_4	2001-07-01 00:00:00	45.0	45.0
+temp_acc_5	2001-09-01 00:00:00	70.0	70.0
+temp_acc_6	2001-11-01 00:00:00	95.0	95.0

temporal/t.rast.accdetect/test_suite/test_1_temp_indi.ref

@@ -0,0 +1,3 @@
+temp_occ_indicator_3	2001-05-01 00:00:00	1.0	1.0
+temp_occ_indicator_4	2001-07-01 00:00:00	2.0	2.0
+temp_occ_indicator_5	2001-09-01 00:00:00	3.0	3.0

temporal/t.rast.accdetect/test_suite/test_1_temp_occ_a.ref

@@ -0,0 +1,3 @@
+temp_occ_3	2001-05-01 00:00:00	120.0	120.0
+temp_occ_4	2001-07-01 00:00:00	181.0	181.0
+temp_occ_5	2001-09-01 00:00:00	243.0	243.0

temporal/t.rast.accdetect/test_suite/test_1_temp_occ_b.ref

@@ -0,0 +1,3 @@
+temp_occ_3	2001-05-01 00:00:00	120.0	120.0
+temp_occ_4	2001-07-01 00:00:00	181.0	181.0
+temp_occ_5	2001-09-01 00:00:00	243.0	243.0

temporal/t.rast.accdetect/test_suite/test_2_temp_accumulation.ref

@@ -0,0 +1,6 @@
+temp_acc_6	2001-01-01 00:00:00	95.0	95.0
+temp_acc_5	2001-03-01 00:00:00	70.0	70.0
+temp_acc_4	2001-05-01 00:00:00	45.0	45.0
+temp_acc_3	2001-07-01 00:00:00	25.0	25.0
+temp_acc_2	2001-09-01 00:00:00	10.0	10.0
+temp_acc_1	2001-11-01 00:00:00	0.0	0.0

temporal/t.rast.accdetect/test_suite/test_2_temp_indi.ref

@@ -0,0 +1,3 @@
+temp_occ_indicator_5	2001-03-01 00:00:00	3.0	3.0
+temp_occ_indicator_4	2001-05-01 00:00:00	2.0	2.0
+temp_occ_indicator_3	2001-07-01 00:00:00	1.0	1.0

temporal/t.rast.accdetect/test_suite/test_2_temp_occ.ref

@@ -0,0 +1,3 @@
+temp_occ_5	2001-03-01 00:00:00	59.0	59.0
+temp_occ_4	2001-05-01 00:00:00	120.0	120.0
+temp_occ_3	2001-07-01 00:00:00	181.0	181.0

temporal/t.rast.accumulate/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.accumulate
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)

temporal/t.rast.accumulate/t.rast.accumulate.html

@@ -0,0 +1,271 @@
+<h2>DESCRIPTION</h2>
+
+<b>t.rast.accumulate</b> is designed to perform temporal accumulations of space time raster datasets.
+
+This module expects a space time raster dataset as input that will be sampled by a given <b>granularity</b>.
+All maps that have the start time during the actual granule will be accumulated with the predecessor 
+granule accumulation result using the raster module <a href="r.series.accumulate.html">r.series.accumulate</a>. 
+The default granularity is 1 day, but any temporal granularity can be set.
+<p>
+The <b>start</b> time and the <b>end</b> time of the accumulation process must be set, 
+eg. <b>start="2000-03-01" end="2011-01-01"</b>. 
+In addition a <b>cycle</b>, eg. <b>cycle="8 months"</b>, can be specified, that defines
+after which time interval the accumulation process restarts. The <b>offset</b> option specifies the
+time between two cycles that should be skipped, eg. <b>offset="4 months"</b>.
+<p>
+The <b>lower</b> and <b>upper</b> <b>limits</b> of the accumulation 
+process can be set, either by using space time raster datasets or 
+by using fixed values for all raster cells and time steps.  The raster maps that specifies the
+lower and upper limits of the actual granule will be detected using the following temporal relations: 
+equals, during, overlaps, overlapped and contains. 
+First all maps with equal time stamps to the current granule will be detected, 
+the first lower map and the first upper map that were found are used as limit definitions. 
+If no equal maps are found then maps with a temporal during relation are detected, 
+then maps that temporally overlap the actual granules, until
+maps are detected that have a temporal contain relation. 
+If no maps are found or lower/upper STRDS are not defined, then the <b>limits</b> option is used, eg. <b>limits=10,30</b>.
+<p>
+The options <b>shift</b>, <b>scale</b> and <b>method</b> are passed to 
+<a href="r.series.accumulate.html">r.series.accumulate</a>. 
+Please refer to the manual page of <a href="r.series.accumulate.html">r.series.accumulate</a>
+for detailed option description.
+<p>
+The <b>output</b> is a new space time raster dataset with the provided start time, end time and granularity
+containing the accumulated raster maps. The <b>base</b> name of of the generated maps must always be set.
+The <b>output</b> space time raster dataset can then be analyzed 
+using <a href="t.rast.accdetect.html">t.rast.accdetect</a> to detect specific accumulation patterns.
+
+<h2>EXAMPLE</h2>
+
+This is an example how to accumulate the daily mean temperature of Europe from 1990 to 2000 using
+the growing-degree-day method to detect grass hopper reproduction cycles that are critical 
+to agriculture. 
+
+<div class="code"><pre>
+# Get the temperature data
+wget www-pool.math.tu-berlin.de/~soeren/grass/temperature_mean_1990_2000_daily_celsius.tar.gz
+
+# Create a temporary location directory
+mkdir -p /tmp/grassdata/LL
+
+# Start GRASS and create a new location with PERMANENT mapset
+grass70  -c EPSG:4326 /tmp/grassdata/LL/PERMANENT
+
+# Import the temperature data
+t.rast.import in=temperature_mean_1990_2000_daily_celsius.tar.gz \
+      out=temperature_mean_1990_2000_daily_celsius extr=/tmp
+
+# We need to set the region correctly
+g.region -p rast=`t.rast.list input=temperature_mean_1990_2000_daily_celsius column=name | tail -1`
+# We can zoom to the raster map
+g.region -p zoom=`t.rast.list input=temperature_mean_1990_2000_daily_celsius column=name | tail -1`
+
+#############################################################################
+#### ACCUMULATION USING GDD METHOD ##########################################
+#############################################################################
+# The computation of grashopper pest control cycles is based on:
+#
+#   Using Growing Degree Days For Insect Management
+#   Nancy E. Adams
+#   Extension Educator, Agricultural Resources
+#
+# available here: http://extension.unh.edu/agric/gddays/docs/growch.pdf
+
+# Now we compute the GDD from 1990 - 2000 for each year (12 month cycle) with
+# a granularity of one day. Base temperature is 10°C, upper limit is 30°C.
+# Hence the accumulation starts at 10°C and does not accumulate values above 30°C. 
+t.rast.accumulate input="temperature_mean_1990_2000_daily_celsius" \
+      output="temperature_mean_1990_2000_daily_celsius_accumulated_10_30" \
+      limits="10,30" start="1990-01-01" stop="2000-01-01" cycle="12 months" \
+      base="temp_acc_daily_10_30" method="gdd"
+
+#############################################################################
+#### ACCUMULATION PATTERN DETECTION #########################################
+#############################################################################
+# Now we detect the three grasshopper pest control cycles
+
+# First cycle at 325°C - 427°C GDD
+t.rast.accdetect input=temperature_mean_1990_2000_daily_celsius_accumulated_10_30@PERMANENT \
+      occ=leafhopper_occurrence_c1_1990_2000 start="1990-01-01" stop="2000-01-01" \
+      cycle="12 months" range=325,427 base=lh_c1 indicator=leafhopper_indicator_c1_1990_2000
+
+# Second cycle at 685°C - 813°C GDD
+t.rast.accdetect input=temperature_mean_1990_2000_daily_celsius_accumulated_10_30@PERMANENT \
+      occ=leafhopper_occurrence_c2_1990_2000 start="1990-01-01" stop="2000-01-01" \
+      cycle="12 months" range=685,813 base=lh_c2 indicator=leafhopper_indicator_c2_1990_2000
+
+# Third cycle at 1047°C - 1179°C GDD
+t.rast.accdetect input=temperature_mean_1990_2000_daily_celsius_accumulated_10_30@PERMANENT \
+      occ=leafhopper_occurrence_c3_1990_2000 start="1990-01-01" stop="2000-01-01" \
+      cycle="12 months" range=1047,1179 base=lh_c3 indicator=leafhopper_indicator_c3_1990_2000
+
+
+#############################################################################
+#### YEARLY SPATIAL OCCURRENCE COMPUTATION OF ALL CYCLES ####################
+#############################################################################
+
+# Extract the areas that have full cycles
+t.rast.aggregate input=leafhopper_indicator_c1_1990_2000 gran="1 year" \
+      output=leafhopper_cycle_1_1990_2000_yearly method=maximum base=li_c1
+
+t.rast.mapcalc input=leafhopper_cycle_1_1990_2000_yearly base=lh_clean_c1 \
+               output=leafhopper_cycle_1_1990_2000_yearly_clean \
+               expr="if(leafhopper_cycle_1_1990_2000_yearly == 3, 1, null())"
+
+t.rast.aggregate input=leafhopper_indicator_c2_1990_2000 gran="1 year" \
+      output=leafhopper_cycle_2_1990_2000_yearly method=maximum base=li_c2
+
+t.rast.mapcalc input=leafhopper_cycle_2_1990_2000_yearly base=lh_clean_c2 \
+               output=leafhopper_cycle_2_1990_2000_yearly_clean \
+               expr="if(leafhopper_cycle_2_1990_2000_yearly == 3, 2, null())"
+
+t.rast.aggregate input=leafhopper_indicator_c3_1990_2000 gran="1 year" \
+      output=leafhopper_cycle_3_1990_2000_yearly method=maximum base=li_c3
+
+t.rast.mapcalc input=leafhopper_cycle_3_1990_2000_yearly base=lh_clean_c3 \
+               output=leafhopper_cycle_3_1990_2000_yearly_clean \
+               expr="if(leafhopper_cycle_3_1990_2000_yearly == 3, 3, null())"
+
+
+t.rast.mapcalc input=leafhopper_cycle_1_1990_2000_yearly_clean,leafhopper_cycle_2_1990_2000_yearly_clean,leafhopper_cycle_3_1990_2000_yearly_clean \
+               base=lh_cleann_all_cycles \
+               output=leafhopper_all_cycles_1990_2000_yearly_clean \
+               expr="if(isnull(leafhopper_cycle_3_1990_2000_yearly_clean), if(isnull(leafhopper_cycle_2_1990_2000_yearly_clean), if(isnull(leafhopper_cycle_1_1990_2000_yearly_clean), null() ,1),2),3)"
+
+cat > color.table << EOF
+3 yellow
+2 blue
+1 red
+EOF
+
+t.rast.colors input=leafhopper_cycle_1_1990_2000_yearly_clean rules=color.table
+t.rast.colors input=leafhopper_cycle_2_1990_2000_yearly_clean rules=color.table
+t.rast.colors input=leafhopper_cycle_3_1990_2000_yearly_clean rules=color.table
+t.rast.colors input=leafhopper_all_cycles_1990_2000_yearly_clean rules=color.table
+
+#############################################################################
+################ DURATION COMPUTATION #######################################
+#############################################################################
+
+# Extract the duration in days of the first cycle
+t.rast.aggregate input=leafhopper_occurrence_c1_1990_2000 gran="1 year" \
+      output=leafhopper_min_day_c1_1990_2000 method=minimum base=occ_min_day_c1
+t.rast.aggregate input=leafhopper_occurrence_c1_1990_2000 gran="1 year" \
+      output=leafhopper_max_day_c1_1990_2000 method=maximum base=occ_max_day_c1
+t.rast.mapcalc input=leafhopper_min_day_c1_1990_2000,leafhopper_max_day_c1_1990_2000 \
+               base=occ_duration_c1 \
+               output=leafhopper_duration_c1_1990_2000 \
+               expr="leafhopper_max_day_c1_1990_2000 - leafhopper_min_day_c1_1990_2000"
+
+
+# Extract the duration in days of the second cycle
+t.rast.aggregate input=leafhopper_occurrence_c2_1990_2000 gran="1 year" \
+      output=leafhopper_min_day_c2_1990_2000 method=minimum base=occ_min_day_c2
+t.rast.aggregate input=leafhopper_occurrence_c2_1990_2000 gran="1 year" \
+      output=leafhopper_max_day_c2_1990_2000 method=maximum base=occ_max_day_c2
+t.rast.mapcalc input=leafhopper_min_day_c2_1990_2000,leafhopper_max_day_c2_1990_2000 \
+               base=occ_duration_c2 \
+               output=leafhopper_duration_c2_1990_2000 \
+               expr="leafhopper_max_day_c2_1990_2000 - leafhopper_min_day_c2_1990_2000"
+
+
+# Extract the duration in days of the third cycle
+t.rast.aggregate input=leafhopper_occurrence_c3_1990_2000 gran="1 year" \
+      output=leafhopper_min_day_c3_1990_2000 method=minimum base=occ_min_day_c3
+t.rast.aggregate input=leafhopper_occurrence_c3_1990_2000 gran="1 year" \
+      output=leafhopper_max_day_c3_1990_2000 method=maximum base=occ_max_day_c3
+t.rast.mapcalc input=leafhopper_min_day_c3_1990_2000,leafhopper_max_day_c3_1990_2000 \
+               base=occ_duration_c3 \
+               output=leafhopper_duration_c3_1990_2000 \
+               expr="leafhopper_max_day_c3_1990_2000 - leafhopper_min_day_c3_1990_2000"
+
+t.rast.colors input=leafhopper_duration_c1_1990_2000 color=rainbow
+t.rast.colors input=leafhopper_duration_c2_1990_2000 color=rainbow
+t.rast.colors input=leafhopper_duration_c3_1990_2000 color=rainbow
+
+#############################################################################
+################ MONTHLY CYCLES OCCURRENCE ##################################
+#############################################################################
+
+# Extract the monthly indicator that show the start and end of a cycle
+
+# First cycle
+t.rast.aggregate input=leafhopper_indicator_c1_1990_2000 gran="1 month" \
+      output=leafhopper_indi_min_month_c1_1990_2000 method=minimum base=occ_indi_min_month_c1
+t.rast.aggregate input=leafhopper_indicator_c1_1990_2000 gran="1 month" \
+      output=leafhopper_indi_max_month_c1_1990_2000 method=maximum base=occ_indi_max_month_c1
+t.rast.mapcalc input=leafhopper_indi_min_month_c1_1990_2000,leafhopper_indi_max_month_c1_1990_2000 \
+               base=indicator_monthly_c1 \
+               output=leafhopper_monthly_indicator_c1_1990_2000 \
+               expr="if(leafhopper_indi_min_month_c1_1990_2000 == 1, 1, if(leafhopper_indi_max_month_c1_1990_2000 == 3, 3, 2))"
+
+# Second cycle
+t.rast.aggregate input=leafhopper_indicator_c2_1990_2000 gran="1 month" \
+      output=leafhopper_indi_min_month_c2_1990_2000 method=minimum base=occ_indi_min_month_c2
+t.rast.aggregate input=leafhopper_indicator_c2_1990_2000 gran="1 month" \
+      output=leafhopper_indi_max_month_c2_1990_2000 method=maximum base=occ_indi_max_month_c2
+t.rast.mapcalc input=leafhopper_indi_min_month_c2_1990_2000,leafhopper_indi_max_month_c2_1990_2000 \
+               base=indicator_monthly_c2 \
+               output=leafhopper_monthly_indicator_c2_1990_2000 \
+               expr="if(leafhopper_indi_min_month_c2_1990_2000 == 1, 1, if(leafhopper_indi_max_month_c2_1990_2000 == 3, 3, 2))"
+
+# Third cycle
+t.rast.aggregate input=leafhopper_indicator_c3_1990_2000 gran="1 month" \
+      output=leafhopper_indi_min_month_c3_1990_2000 method=minimum base=occ_indi_min_month_c3
+t.rast.aggregate input=leafhopper_indicator_c3_1990_2000 gran="1 month" \
+      output=leafhopper_indi_max_month_c3_1990_2000 method=maximum base=occ_indi_max_month_c3
+t.rast.mapcalc input=leafhopper_indi_min_month_c3_1990_2000,leafhopper_indi_max_month_c3_1990_2000 \
+               base=indicator_monthly_c3 \
+               output=leafhopper_monthly_indicator_c3_1990_2000 \
+               expr="if(leafhopper_indi_min_month_c3_1990_2000 == 1, 1, if(leafhopper_indi_max_month_c3_1990_2000 == 3, 3, 2))"
+
+cat > color.table << EOF
+3 red
+2 yellow
+1 green
+EOF
+
+t.rast.colors input=leafhopper_monthly_indicator_c1_1990_2000 rules=color.table
+t.rast.colors input=leafhopper_monthly_indicator_c2_1990_2000 rules=color.table
+t.rast.colors input=leafhopper_monthly_indicator_c3_1990_2000 rules=color.table
+
+#############################################################################
+################ VISUALIZATION ##############################################
+#############################################################################
+# Now we use g.gui.animation to visualize the yearly occurrence, the duration and the monthly occurrence
+
+# Yearly occurrence of all reproduction cycles 
+g.gui.animation strds=leafhopper_all_cycles_1990_2000_yearly_clean
+
+# Yearly duration of reproduction cycle 1
+g.gui.animation strds=leafhopper_duration_c1_1990_2000
+# Yearly duration of reproduction cycle 2
+g.gui.animation strds=leafhopper_duration_c2_1990_2000
+# Yearly duration of reproduction cycle 3
+g.gui.animation strds=leafhopper_duration_c3_1990_2000
+
+# Monthly occurrence of reproduction cycle 1
+g.gui.animation strds=leafhopper_monthly_indicator_c1_1990_2000
+# Monthly occurrence of reproduction cycle 2
+g.gui.animation strds=leafhopper_monthly_indicator_c2_1990_2000
+# Monthly occurrence of reproduction cycle 3
+g.gui.animation strds=leafhopper_monthly_indicator_c3_1990_2000
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<b>
+<a href="t.rast.accdetect.html">t.rast.accdetect</a><br>
+<a href="t.rast.aggregate.html">t.rast.aggregate</a><br>
+<a href="t.rast.mapcalc.html">t.rast.mapcalc</a><br>
+<a href="t.info.html">t.info</a><br>
+<a href="g.region.html">g.region</a><br>
+<a href="r.series.accumulate.html">r.series.accumulate</a><br>
+</b>
+
+
+<h2>AUTHOR</h2>
+
+S&ouml;ren Gebbert
+
+<p><i>Last changed: $Date: 2012-08-29 14:55:21 +0200 (Mi, 29. Aug 2012) $</i>

temporal/t.rast.accumulate/t.rast.accumulate.py

@@ -0,0 +1,500 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:       t.rast.accumulate
+# AUTHOR(S):    Soeren Gebbert
+#
+# PURPOSE:      Compute cyclic accumulations of a space time raster dataset
+# COPYRIGHT:    (C) 2013 - 2014 by the GRASS Development Team
+#
+#               This program is free software under the GNU General Public
+#               License (version 2). Read the file COPYING that comes with GRASS
+#               for details.
+#
+#############################################################################
+
+#%module
+#% description: Compute cyclic accumulations of a space time raster dataset
+#% keywords: temporal
+#% keywords: accumulation
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_STRDS_OUTPUT
+#%end
+
+
+#%option G_OPT_STRDS_INPUT
+#% key: lower
+#% description: Input space time raster dataset that defines the lower threshold, values lower this threshold are excluded from accumulation
+#% required: no
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#% key: upper
+#% description: Input space time raster dataset that defines the upper threshold, values upper this threshold are excluded from accumulation
+#% required: no
+#%end
+
+#%option
+#% key: start
+#% type: string
+#% description: The temporal starting point to begin the accumulation, eg '2001-01-01'
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: stop
+#% type: string
+#% description: The temporal date to stop the accumulation, eg '2009-01-01'
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: cycle
+#% type: string
+#% description: The temporal cycle to restart the accumulation, eg '12 months'
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: offset
+#% type: string
+#% description: The temporal offset to the begin of the next cycle, eg '6 months'
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: granularity
+#% type: string
+#% description: The granularity for accumulation '1 day'
+#% answer: 1 day
+#% required: no
+#% multiple: no
+#%end
+
+#%option G_OPT_R_BASE
+#%end
+
+#%option
+#% key: limits
+#% type: double
+#% key_desc: lower,upper
+#% description: Use these limits in case lower and/or upper input  space time raster datasets are not defined
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: shift
+#% type: double
+#% description: Scale factor for input space time raster dataset
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: scale
+#% type: double
+#% description: Shift factor for input space time raster dataset
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: method
+#% type: string
+#% description: This method will be applied to compute the accumulative values from the input maps
+#% options: mean,gdd,winkler
+#% answer: mean
+#% required: no
+#% multiple: no
+#%end
+
+#%flag
+#% key: n
+#% description: Register empty maps in the output space time raster dataset, otherwise they will be deleted
+#%end
+
+#%flag
+#% key: r
+#% description: Reverse time direction in cyclic accumulation
+#%end
+
+import grass.script as grass
+import grass.temporal as tgis
+from grass.pygrass.modules import Module
+from copy import copy
+
+############################################################################
+
+def main():
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    start = options["start"]
+    stop = options["stop"]
+    base = options["base"]
+    cycle = options["cycle"]
+    lower = options["lower"]
+    upper = options["upper"]
+    offset = options["offset"]
+    limits = options["limits"]
+    shift = options["shift"]
+    scale = options["scale"]
+    method = options["method"]
+    granularity = options["granularity"]
+    register_null = flags["n"]
+    reverse = flags["r"]
+
+    # Make sure the temporal database exists
+    tgis.init()
+
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    mapset = tgis.get_current_mapset()
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        id = input + "@" + mapset
+
+    input_strds = tgis.SpaceTimeRasterDataset(id)
+
+    if input_strds.is_in_db() == False:
+        dbif.close()
+        grass.fatal(_("Space time raster dataset <%s> not found") % (id))
+
+    input_strds.select(dbif)
+
+    if output.find("@") >= 0:
+        out_id = output
+    else:
+        out_id = output + "@" + mapset
+
+    # The output space time raster dataset
+    output_strds = tgis.SpaceTimeRasterDataset(out_id)
+    if output_strds.is_in_db(dbif):
+        if not grass.overwrite():
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> is already in the "
+                          "database, use overwrite flag to overwrite") % out_id)
+
+    if tgis.check_granularity_string(granularity,
+                                     input_strds.get_temporal_type()) == False:
+            dbif.close()
+            grass.fatal(_("Invalid granularity"))
+
+    if tgis.check_granularity_string(cycle,
+                                     input_strds.get_temporal_type()) == False:
+            dbif.close()
+            grass.fatal(_("Invalid cycle"))
+
+    if offset:
+        if tgis.check_granularity_string(offset,
+                                         input_strds.get_temporal_type()) == False:
+                dbif.close()
+                grass.fatal(_("Invalid offset"))
+
+    # The lower threshold space time raster dataset
+    if lower:
+        if not range:
+            dbif.close()
+            grass.fatal(_("You need to set the range to compute the occurrence"
+                          " space time raster dataset"))
+
+        if lower.find("@") >= 0:
+            lower_id = lower
+        else:
+            lower_id = lower + "@" + mapset
+
+        lower_strds = tgis.SpaceTimeRasterDataset(lower_id)
+        if lower_strds.is_in_db() == False:
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> not found") % (lower_strds.get_id()))
+
+        if lower_strds.get_temporal_type() != input_strds.get_temporal_type():
+            dbif.close()
+            grass.fatal(_("Temporal type of input strds and lower strds must be equal"))
+
+        lower_strds.select(dbif)
+
+    # The upper threshold space time raster dataset
+    if upper:
+        if not lower:
+            dbif.close()
+            grass.fatal(_("The upper option works only in conjunction with the lower option"))
+
+        if upper.find("@") >= 0:
+            upper = upper
+        else:
+            upper_id = upper + "@" + mapset
+
+        upper_strds = tgis.SpaceTimeRasterDataset(upper_id)
+        if upper_strds.is_in_db() == False:
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> not found") % (upper_strds.get_id()))
+
+        if upper_strds.get_temporal_type() != input_strds.get_temporal_type():
+            dbif.close()
+            grass.fatal(_("Temporal type of input strds and upper strds must be equal"))
+
+        upper_strds.select(dbif)
+
+    input_strds_start, input_strds_end = input_strds.get_temporal_extent_as_tuple()
+
+    if input_strds.is_time_absolute():
+        start = tgis.string_to_datetime(start)
+        if stop:
+            stop = tgis.string_to_datetime(stop)
+        else:
+            stop = input_strds_end
+        start = tgis.adjust_datetime_to_granularity(start, granularity)
+    else:
+        start = int(start)
+        if stop:
+            stop = int(stop)
+        else:
+            stop = input_strds_end
+
+    if input_strds.is_time_absolute():
+        end = tgis.increment_datetime_by_string(start, cycle)
+    else:
+        end = start + cycle
+
+    limit_relations = ["EQUALS", "DURING", "OVERLAPS", "OVERLAPPING", "CONTAINS"]
+
+    count = 1
+    output_maps = []
+
+
+    while input_strds_end > start and stop > start:
+
+        # Make sure that the cyclic computation will stop at the correct time
+        if stop and end > stop:
+            end = stop
+
+        where = "start_time >= \'%s\' AND start_time < \'%s\'"%(str(start),
+                                                                str(end))
+        input_maps = input_strds.get_registered_maps_as_objects(where=where,
+                                                                dbif=dbif)
+
+        grass.message(_("Processing cycle %s - %s"%(str(start), str(end))))
+
+        if len(input_maps) == 0:
+            continue
+
+        # Lets create a dummy list of maps with granularity conform intervals
+        gran_list = []
+        gran_list_low = []
+        gran_list_up = []
+        gran_start = start
+        while gran_start < end:
+            map = input_strds.get_new_map_instance("%i@%i"%(count, count))
+            if input_strds.is_time_absolute():
+                gran_end = tgis.increment_datetime_by_string(gran_start,
+                                                             granularity)
+                map.set_absolute_time(gran_start, gran_end)
+                gran_start = tgis.increment_datetime_by_string(gran_start,
+                                                               granularity)
+            else:
+                gran_end = gran_start + granularity
+                map.set_relative_time(gran_start, gran_end,
+                                      input_strds.get_relative_time_unit())
+                gran_start = gran_start + granularity
+            gran_list.append(copy(map))
+            gran_list_low.append(copy(map))
+            gran_list_up.append(copy(map))
+        # Lists to compute the topology with upper and lower datasets
+
+        # Create the topology between the granularity conform list and all maps
+        # of the current cycle
+        gran_topo = tgis.SpatioTemporalTopologyBuilder()
+        gran_topo.build(gran_list, input_maps)
+
+        if lower:
+            lower_maps = lower_strds.get_registered_maps_as_objects(dbif=dbif)
+            gran_lower_topo = tgis.SpatioTemporalTopologyBuilder()
+            gran_lower_topo.build(gran_list_low, lower_maps)
+
+        if upper:
+            upper_maps = upper_strds.get_registered_maps_as_objects(dbif=dbif)
+            gran_upper_topo = tgis.SpatioTemporalTopologyBuilder()
+            gran_upper_topo.build(gran_list_up, upper_maps)
+
+        old_map_name = None
+
+        # Aggregate
+        num_maps = len(gran_list)
+
+        for i in xrange(num_maps):
+            if reverse:
+                map = gran_list[num_maps - i - 1]
+            else:
+                map = gran_list[i]
+            # Select input maps based on temporal topology relations
+            input_maps = []
+            if map.get_equal():
+                input_maps += map.get_equal()
+            elif map.get_contains():
+                input_maps += map.get_contains()
+            elif map.get_overlaps():
+                input_maps += map.get_overlaps()
+            elif map.get_overlapped():
+                input_maps += map.get_overlapped()
+            elif map.get_during():
+                input_maps += map.get_during()
+
+            # Check input maps
+            if len(input_maps) == 0:
+                continue
+
+            # New output map
+            output_map_name = "%s_%i" % (base, count)
+            output_map_id = map.build_id(output_map_name, mapset)
+            output_map = input_strds.get_new_map_instance(output_map_id)
+
+            # Check if new map is in the temporal database
+            if output_map.is_in_db(dbif):
+                if grass.overwrite():
+                    # Remove the existing temporal database entry
+                    output_map.delete(dbif)
+                    output_map = input_strds.get_new_map_instance(output_map_id)
+                else:
+                    grass.fatal(_("Map <%s> is already registered in the temporal"
+                                 " database, use overwrite flag to overwrite.") %
+                                (output_map.get_map_id()))
+
+            map_start, map_end = map.get_temporal_extent_as_tuple()
+
+            if map.is_time_absolute():
+                output_map.set_absolute_time(map_start, map_end)
+            else:
+                output_map.set_relative_time(map_start, map_end,
+                                             map.get_relative_time_unit())
+
+            limits_vals = limits.split(",")
+            limits_lower = float(limits_vals[0])
+            limits_upper = float(limits_vals[1])
+
+            lower_map_name = None
+            if lower:
+                relations = gran_list_low[i].get_temporal_relations()
+                for relation in limit_relations:
+                    if relation in relations:
+                        lower_map_name = str(relations[relation][0].get_id())
+                        break
+
+            upper_map_name = None
+            if upper:
+                relations = gran_list_up[i].get_temporal_relations()
+                for relation in limit_relations:
+                    if relation in relations:
+                        upper_map_name = str(relations[relation][0].get_id())
+                        break
+
+            input_map_names = []
+            for input_map in input_maps:
+                input_map_names.append(input_map.get_id())
+
+            # Set up the module
+            accmod = Module("r.series.accumulate", input=input_map_names,
+                            output=output_map_name, run_=False)
+
+            if old_map_name:
+                accmod.inputs["basemap"].value = old_map_name
+            if lower_map_name:
+                accmod.inputs["lower"].value = lower_map_name
+            if upper_map_name:
+                accmod.inputs["upper"].value = upper_map_name
+
+            accmod.inputs["limits"].value = (limits_lower, limits_upper)
+
+            if shift:
+                accmod.inputs["shift"].value = float(shift)
+
+            if scale:
+                accmod.inputs["scale"].value = float(scale)
+
+            if method:
+                accmod.inputs["method"].value = method
+
+            print accmod
+            accmod.run()
+
+            if accmod.popen.returncode != 0:
+                dbif.close()
+                grass.fatal(_("Error running r.series.accumulate"))
+
+            output_maps.append(output_map)
+            old_map_name = output_map_name
+            count += 1
+
+        # Increment the cycle
+        start = end
+        if input_strds.is_time_absolute():
+            start = end
+            if offset:
+                start = tgis.increment_datetime_by_string(end, offset)
+
+            end = tgis.increment_datetime_by_string(start, cycle)
+        else:
+            if offset:
+                start = end + offset
+            end = start + cycle
+
+    # Insert the maps into the output space time dataset
+    if output_strds.is_in_db(dbif):
+        if grass.overwrite():
+            output_strds.delete(dbif)
+            output_strds = input_strds.get_new_instance(out_id)
+
+    temporal_type, semantic_type, title, description = input_strds.get_initial_values()
+    output_strds.set_initial_values(temporal_type, semantic_type, title,
+                                    description)
+    output_strds.insert(dbif)
+
+    empty_maps = []
+    # Register the maps in the database
+    count = 0
+    for output_map in output_maps:
+        count += 1
+        if count%10 == 0:
+            grass.percent(count, len(output_maps), 1)
+        # Read the raster map data
+        output_map.load()
+        # In case of a empty map continue, do not register empty maps
+
+        if not register_null:
+            if output_map.metadata.get_min() is None and \
+                output_map.metadata.get_max() is None:
+                empty_maps.append(output_map)
+                continue
+
+        # Insert map in temporal database
+        output_map.insert(dbif)
+        output_strds.register_map(output_map, dbif)
+
+    # Update the spatio-temporal extent and the metadata table entries
+    output_strds.update_from_registered_maps(dbif)
+    grass.percent(1, 1, 1)
+
+    dbif.close()
+
+    # Remove empty maps
+    if len(empty_maps) > 0:
+        for map in empty_maps:
+            grass.run_command("g.remove", rast=map.get_name(), quiet=True)
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()

temporal/t.rast.accumulate/test_suite/acc_1.ref

@@ -0,0 +1,39 @@
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ |                                                                            |
+ +-------------------- Basic information -------------------------------------+
+ | Name: ...................... precip_abs2
+ | Creator: ................... soeren
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2001-01-01 00:00:00
+ | End time:................... 2002-01-01 00:00:00
+ | Granularity:................ 4 months
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 0.0
+ | Maximum value min:.......... 21.5
+ | Maximum value max:.......... 60.0
+ |
+ | Title:
+ | A test
+ | Description:
+ | A test
+ | Command history:
+ | t.rast.accumulate input="precip_abs1"
+ |     output="precip_abs2" base="prec_acc" limits="10,30" start="2001-01-01"
+ |     gran="4 months" cycle="12 months"
+ | 
+ +----------------------------------------------------------------------------+

temporal/t.rast.accumulate/test_suite/acc_2.ref

@@ -0,0 +1,39 @@
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ |                                                                            |
+ +-------------------- Basic information -------------------------------------+
+ | Name: ...................... precip_abs2
+ | Creator: ................... soeren
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2001-01-01 00:00:00
+ | End time:................... 2001-09-01 00:00:00
+ | Granularity:................ 4 months
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 0.0
+ | Maximum value min:.......... 21.5
+ | Maximum value max:.......... 39.5
+ |
+ | Title:
+ | A test
+ | Description:
+ | A test
+ | Command history:
+ | t.rast.accumulate input="precip_abs1"
+ |     output="precip_abs2" base="prec_acc" limits="10,30" start="2001-01-01"
+ |     stop="2002-01-01" gran="4 months" cycle="7 months"
+ | 
+ +----------------------------------------------------------------------------+

temporal/t.rast.accumulate/test_suite/acc_3.ref

@@ -0,0 +1,39 @@
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ |                                                                            |
+ +-------------------- Basic information -------------------------------------+
+ | Name: ...................... precip_abs2
+ | Creator: ................... soeren
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2001-01-01 00:00:00
+ | End time:................... 2002-01-01 00:00:00
+ | Granularity:................ 2 months
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 0.0
+ | Maximum value min:.......... 24.0
+ | Maximum value max:.......... 120.0
+ |
+ | Title:
+ | A test
+ | Description:
+ | A test
+ | Command history:
+ | t.rast.accumulate input="precip_abs1"
+ |     output="precip_abs2" base="prec_acc" limits="10,30" start="2001-01-01"
+ |     gran="2 months" cycle="12 months"
+ | 
+ +----------------------------------------------------------------------------+

temporal/t.rast.accumulate/test_suite/acc_4.ref

@@ -0,0 +1,39 @@
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ |                                                                            |
+ +-------------------- Basic information -------------------------------------+
+ | Name: ...................... precip_abs2
+ | Creator: ................... soeren
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2001-01-01 00:00:00
+ | End time:................... 2002-01-01 00:00:00
+ | Granularity:................ 1 month
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 0.0
+ | Maximum value min:.......... 24.0
+ | Maximum value max:.......... 240.0
+ |
+ | Title:
+ | A test
+ | Description:
+ | A test
+ | Command history:
+ | t.rast.accumulate input="precip_abs1"
+ |     output="precip_abs2" base="prec_acc" limits="10,30" start="2001-01-01"
+ |     stop="2002-01-01" gran="1 months" cycle="12 months"
+ | 
+ +----------------------------------------------------------------------------+

temporal/t.rast.accumulate/test_suite/acc_5.ref

@@ -0,0 +1,40 @@
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ |                                                                            |
+ +-------------------- Basic information -------------------------------------+
+ | Name: ...................... precip_abs2
+ | Creator: ................... soeren
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2001-01-01 00:00:00
+ | End time:................... 2002-01-01 00:00:00
+ | Granularity:................ 1 month
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 0.0
+ | Maximum value min:.......... 24.0
+ | Maximum value max:.......... 240.0
+ |
+ | Title:
+ | A test
+ | Description:
+ | A test
+ | Command history:
+ | t.rast.accumulate input="precip_abs1"
+ |     output="precip_abs2" base="prec_acc" limits="8,33" lower="lower"
+ |     upper="upper" start="2001-01-01" stop="2002-01-01" gran="1 months"
+ |     cycle="12 months"
+ | 
+ +----------------------------------------------------------------------------+

temporal/t.rast.accumulate/test_suite/test.t.rast.accumulate.sh

@@ -0,0 +1,78 @@
+#!/bin/sh -e
+# Space time raster dataset neighborhood operations
+# We need to set a specific region in the
+# @preprocess step of this test. 
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 -p
+
+export GRASS_OVERWRITE=1
+
+# Generate data
+r.mapcalc expr="prec_1 = rand(0, 25)"
+r.mapcalc expr="prec_2 = rand(0, 24)"
+r.mapcalc expr="prec_3 = rand(0, 23)"
+r.mapcalc expr="prec_4 = rand(0, 25)"
+r.mapcalc expr="prec_5 = rand(0, 23)"
+r.mapcalc expr="prec_6 = rand(0, 26)"
+
+t.create type=strds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
+t.register -i type=rast input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 \
+    start="2001-01-01" increment="2 months"
+
+# The first @test
+
+t.rast.accumulate input=precip_abs1 output=precip_abs2 base=prec_acc \
+    limits=10,30 start="2001-01-01" gran="4 months" \
+    cycle="12 months"
+
+# We need to filter the identifier that may be different on different systems
+# like creation time and mapset 
+t.info input=precip_abs2 | grep -v \# | grep -v Creation | grep -v register | grep -v Id | grep -v Mapset > acc_1.txt
+
+t.rast.accumulate input=precip_abs1 output=precip_abs2 base=prec_acc \
+    limits=10,30 start="2001-01-01" stop="2002-01-01" gran="4 months" \
+    cycle="7 months"
+
+t.info input=precip_abs2 | grep -v \# | grep -v Creation | grep -v register | grep -v Id | grep -v Mapset > acc_2.txt
+
+t.rast.accumulate input=precip_abs1 output=precip_abs2 base=prec_acc \
+    limits=10,30 start="2001-01-01" gran="2 months" \
+    cycle="12 months"
+
+t.info input=precip_abs2 | grep -v \# | grep -v Creation | grep -v register | grep -v Id | grep -v Mapset > acc_3.txt
+
+t.rast.accumulate input=precip_abs1 output=precip_abs2 base=prec_acc \
+    limits=10,30 start="2001-01-01" stop="2002-01-01" gran="1 months" \
+    cycle="12 months"
+
+t.info input=precip_abs2 | grep -v \# | grep -v Creation | grep -v register | grep -v Id | grep -v Mapset > acc_4.txt
+
+# Second test
+
+r.mapcalc expr="lower = 10"
+r.mapcalc expr="upper = 35"
+
+t.create type=strds temporaltype=absolute output=lower title="lower limit" descr="lower limit"
+t.register -i type=rast input=lower maps=lower start="2001-01-01" increment="8 months"
+
+t.create type=strds temporaltype=absolute output=upper title="upper limit" descr="upper limit"
+t.register -i type=rast input=upper maps=upper start="2001-01-01" increment="10 months"
+
+t.rast.accumulate input=precip_abs1 output=precip_abs2 base=prec_acc \
+    limits=8,33 lower=lower upper=upper start="2001-01-01" stop="2002-01-01" gran="1 months" \
+    cycle="12 months"
+
+t.info input=precip_abs2 | grep -v \# | grep -v Creation | grep -v register | grep -v Id | grep -v Mapset > acc_5.txt
+
+t.remove -rf type=strds input=precip_abs1,precip_abs2,lower,upper
+
+for i in `ls acc_*.txt` ; do
+    diff $i "`basename $i .txt`.ref" >> out.diff
+done
+
+CHAR_NUM=`cat out.diff | wc -c`
+
+# Return as exit status 0 in case no diffs are found
+exit $CHAR_NUM
+
+