[GRASS-SVN] r58651 - in grass/trunk: raster raster/r.series.accumulate raster/r.series.accumulate/test_suite temporal temporal/t.rast.accdetect temporal/t.rast.accdetect/test_suite temporal/t.rast.accumulate temporal/t.rast.accumulate/test_suite
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jan 8 18:10:13 PST 2014
Author: huhabla
Date: 2014-01-08 18:10:13 -0800 (Wed, 08 Jan 2014)
New Revision: 58651
Added:
grass/trunk/raster/r.series.accumulate/
grass/trunk/raster/r.series.accumulate/Makefile
grass/trunk/raster/r.series.accumulate/main.c
grass/trunk/raster/r.series.accumulate/r.series.accumulate.html
grass/trunk/raster/r.series.accumulate/test_suite/
grass/trunk/raster/r.series.accumulate/test_suite/test.r.series.accumulate.sh
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_1.ref
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_2.ref
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_3.ref
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_4.ref
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_5.ref
grass/trunk/raster/r.series.accumulate/test_suite/test_accu_6.ref
grass/trunk/temporal/t.rast.accdetect/
grass/trunk/temporal/t.rast.accdetect/Makefile
grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.html
grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.py
grass/trunk/temporal/t.rast.accdetect/test_suite/
grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.reverse.sh
grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.sh
grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_accumulation.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_indi.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_a.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_b.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_accumulation.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_indi.ref
grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_occ.ref
grass/trunk/temporal/t.rast.accumulate/
grass/trunk/temporal/t.rast.accumulate/Makefile
grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.html
grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.py
grass/trunk/temporal/t.rast.accumulate/test_suite/
grass/trunk/temporal/t.rast.accumulate/test_suite/acc_1.ref
grass/trunk/temporal/t.rast.accumulate/test_suite/acc_2.ref
grass/trunk/temporal/t.rast.accumulate/test_suite/acc_3.ref
grass/trunk/temporal/t.rast.accumulate/test_suite/acc_4.ref
grass/trunk/temporal/t.rast.accumulate/test_suite/acc_5.ref
grass/trunk/temporal/t.rast.accumulate/test_suite/test.t.rast.accumulate.sh
Modified:
grass/trunk/raster/Makefile
grass/trunk/temporal/Makefile
Log:
New temporal accumulation modules.
Modified: grass/trunk/raster/Makefile
===================================================================
--- grass/trunk/raster/Makefile 2014-01-09 02:03:39 UTC (rev 58650)
+++ grass/trunk/raster/Makefile 2014-01-09 02:10:13 UTC (rev 58651)
@@ -88,6 +88,7 @@
r.rescale.eq \
r.ros \
r.series \
+ r.series.accumulate \
r.series.interp \
r.shaded.relief \
r.slope.aspect \
Added: grass/trunk/raster/r.series.accumulate/Makefile
===================================================================
--- grass/trunk/raster/r.series.accumulate/Makefile (rev 0)
+++ grass/trunk/raster/r.series.accumulate/Makefile 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/main.c
===================================================================
--- grass/trunk/raster/r.series.accumulate/main.c (rev 0)
+++ grass/trunk/raster/r.series.accumulate/main.c 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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);
+}
\ No newline at end of file
Added: grass/trunk/raster/r.series.accumulate/r.series.accumulate.html
===================================================================
--- grass/trunk/raster/r.series.accumulate/r.series.accumulate.html (rev 0)
+++ grass/trunk/raster/r.series.accumulate/r.series.accumulate.html 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
+# <domain> <type> <item> <value>
+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>
Added: grass/trunk/raster/r.series.accumulate/test_suite/test.r.series.accumulate.sh
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test.r.series.accumulate.sh (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test.r.series.accumulate.sh 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
\ No newline at end of file
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_1.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_1.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_1.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_2.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_2.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_2.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_3.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_3.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_3.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_4.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_4.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_4.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_5.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_5.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_5.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/raster/r.series.accumulate/test_suite/test_accu_6.ref
===================================================================
--- grass/trunk/raster/r.series.accumulate/test_suite/test_accu_6.ref (rev 0)
+++ grass/trunk/raster/r.series.accumulate/test_suite/test_accu_6.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Modified: grass/trunk/temporal/Makefile
===================================================================
--- grass/trunk/temporal/Makefile 2014-01-09 02:03:39 UTC (rev 58650)
+++ grass/trunk/temporal/Makefile 2014-01-09 02:10:13 UTC (rev 58651)
@@ -12,6 +12,8 @@
t.sample \
t.register \
t.unregister \
+ t.rast.accumulate \
+ t.rast.accdetect \
t.rast.aggregate \
t.rast.aggregate.ds \
t.rast.colors \
Added: grass/trunk/temporal/t.rast.accdetect/Makefile
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/Makefile (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/Makefile 2014-01-09 02:10:13 UTC (rev 58651)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.accdetect
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)
Added: grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.html
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.html (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.html 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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ören Gebbert
+
+<p><i>Last changed: $Date: 2012-08-29 14:55:21 +0200 (Mi, 29. Aug 2012) $</i>
Added: grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.py
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.py (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/t.rast.accdetect.py 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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()
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.reverse.sh
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.reverse.sh (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.reverse.sh 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.sh
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.sh (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test.t.rast.accdetect.sh 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_accumulation.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_accumulation.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_accumulation.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_indi.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_indi.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_indi.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_a.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_a.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_a.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_b.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_b.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_1_temp_occ_b.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_accumulation.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_accumulation.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_accumulation.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_indi.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_indi.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_indi.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_occ.ref
===================================================================
--- grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_occ.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accdetect/test_suite/test_2_temp_occ.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
Added: grass/trunk/temporal/t.rast.accumulate/Makefile
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/Makefile (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/Makefile 2014-01-09 02:10:13 UTC (rev 58651)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.accumulate
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)
Added: grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.html
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.html (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.html 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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 at 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 at 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 at 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ören Gebbert
+
+<p><i>Last changed: $Date: 2012-08-29 14:55:21 +0200 (Mi, 29. Aug 2012) $</i>
Added: grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.py
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.py (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/t.rast.accumulate.py 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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()
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/acc_1.ref
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/acc_1.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/acc_1.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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"
+ |
+ +----------------------------------------------------------------------------+
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/acc_2.ref
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/acc_2.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/acc_2.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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"
+ |
+ +----------------------------------------------------------------------------+
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/acc_3.ref
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/acc_3.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/acc_3.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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"
+ |
+ +----------------------------------------------------------------------------+
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/acc_4.ref
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/acc_4.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/acc_4.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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"
+ |
+ +----------------------------------------------------------------------------+
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/acc_5.ref
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/acc_5.ref (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/acc_5.ref 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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"
+ |
+ +----------------------------------------------------------------------------+
Added: grass/trunk/temporal/t.rast.accumulate/test_suite/test.t.rast.accumulate.sh
===================================================================
--- grass/trunk/temporal/t.rast.accumulate/test_suite/test.t.rast.accumulate.sh (rev 0)
+++ grass/trunk/temporal/t.rast.accumulate/test_suite/test.t.rast.accumulate.sh 2014-01-09 02:10:13 UTC (rev 58651)
@@ -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
+
+
More information about the grass-commit
mailing list