raster modules: Implemented computational region settings based on disjoint union and intersection computation from all input raster maps in an expression

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

raster/r.mapcalc/evaluate.c

@@ -1,6 +1,7 @@
 
 #include <stdlib.h>
 #include <unistd.h>
+#include <string.h>
 
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -15,10 +16,56 @@
 int current_depth, current_row;
 int depths, rows;
 
+/* Local variables for map management */
+static expression **map_list = NULL;
+static int num_maps = 0;
+static int max_maps = 0;
+
 /****************************************************************************/
 
+static void extract_maps(expression *e);
 static void initialize(expression *e);
 static void evaluate(expression *e);
+static int append_map(expression *e);
+
+/****************************************************************************/
+
+int append_map(expression *e)
+{
+    /* \brief Append a map to the global map list and reallocate if necessary
+     */
+    const char *mapset;
+
+    if (num_maps >= max_maps) {
+        max_maps += 10;
+        map_list = G_realloc(map_list, max_maps * sizeof(struct expression*));
+    }
+
+    map_list[num_maps] = e;
+
+    return num_maps++;
+}
+
+void extract_maps(expression * e)
+{
+    /* \brief Search for map names in the expression and add them to the global map list */
+    int i;
+
+    switch (e->type) {
+        case expr_type_map:
+            G_debug(1, "Found map %s", e->data.map.name);
+            append_map(e);
+            break;
+        case expr_type_function:
+            for (i = 1; i <= e->data.func.argc; i++) {
+                extract_maps(e->data.func.args[i]);
+            }
+            break;
+        case expr_type_binding:
+            extract_maps(e->data.bind.val);
+            break;
+    }
+}
 
 /****************************************************************************/
 
@@ -47,8 +94,9 @@ static void initialize_variable(expression * e)
 static void initialize_map(expression * e)
 {
     allocate_buf(e);
+
     e->data.map.idx = open_map(e->data.map.name, e->data.map.mod,
-			       e->data.map.row, e->data.map.col);
+                               e->data.map.row, e->data.map.col);
 }
 
 static void initialize_function(expression * e)
@@ -61,8 +109,8 @@ static void initialize_function(expression * e)
     e->data.func.argv[0] = e->buf;
 
     for (i = 1; i <= e->data.func.argc; i++) {
-	initialize(e->data.func.args[i]);
-	e->data.func.argv[i] = e->data.func.args[i]->buf;
+        initialize(e->data.func.args[i]);
+        e->data.func.argv[i] = e->data.func.args[i]->buf;
     }
 }
 
@@ -74,24 +122,25 @@ static void initialize_binding(expression * e)
 
 static void initialize(expression * e)
 {
+
     switch (e->type) {
-    case expr_type_constant:
-	initialize_constant(e);
-	break;
-    case expr_type_variable:
-	initialize_variable(e);
-	break;
-    case expr_type_map:
-	initialize_map(e);
-	break;
-    case expr_type_function:
-	initialize_function(e);
-	break;
-    case expr_type_binding:
-	initialize_binding(e);
-	break;
-    default:
-	G_fatal_error(_("Unknown type: %d"), e->type);
+        case expr_type_constant:
+            initialize_constant(e);
+            break;
+        case expr_type_variable:
+            initialize_variable(e);
+            break;
+        case expr_type_map:
+            initialize_map(e);
+            break;
+        case expr_type_function:
+            initialize_function(e);
+            break;
+        case expr_type_binding:
+            initialize_binding(e);
+            break;
+        default:
+            G_fatal_error(_("Unknown type: %d"), e->type);
     }
 }
 
@@ -123,21 +172,21 @@ static void evaluate_constant(expression * e)
 
     switch (e->res_type) {
     case CELL_TYPE:
-	for (i = 0; i < columns; i++)
-	    ibuf[i] = e->data.con.ival;
-	break;
+        for (i = 0; i < columns; i++)
+            ibuf[i] = e->data.con.ival;
+        break;
 
     case FCELL_TYPE:
-	for (i = 0; i < columns; i++)
-	    fbuf[i] = e->data.con.fval;
-	break;
+        for (i = 0; i < columns; i++)
+            fbuf[i] = e->data.con.fval;
+        break;
 
     case DCELL_TYPE:
-	for (i = 0; i < columns; i++)
-	    dbuf[i] = e->data.con.fval;
-	break;
+        for (i = 0; i < columns; i++)
+            dbuf[i] = e->data.con.fval;
+        break;
     default:
-	G_fatal_error(_("Invalid type: %d"), e->res_type);
+        G_fatal_error(_("Invalid type: %d"), e->res_type);
     }
 }
 
@@ -149,10 +198,10 @@ static void evaluate_variable(expression * e)
 static void evaluate_map(expression * e)
 {
     get_map_row(e->data.map.idx,
-		e->data.map.mod,
-		current_depth + e->data.map.depth,
-		current_row + e->data.map.row,
-		e->data.map.col, e->buf, e->res_type);
+                e->data.map.mod,
+                current_depth + e->data.map.depth,
+                current_row + e->data.map.row,
+                e->data.map.col, e->buf, e->res_type);
 }
 
 static void evaluate_function(expression * e)
@@ -215,23 +264,23 @@ static void evaluate_binding(expression * e)
 static void evaluate(expression * e)
 {
     switch (e->type) {
-    case expr_type_constant:
-	evaluate_constant(e);
-	break;
-    case expr_type_variable:
-	evaluate_variable(e);
-	break;
-    case expr_type_map:
-	evaluate_map(e);
-	break;
-    case expr_type_function:
-	evaluate_function(e);
-	break;
-    case expr_type_binding:
-	evaluate_binding(e);
-	break;
-    default:
-	G_fatal_error(_("Unknown type: %d"), e->type);
+        case expr_type_constant:
+            evaluate_constant(e);
+            break;
+        case expr_type_variable:
+            evaluate_variable(e);
+            break;
+        case expr_type_map:
+            evaluate_map(e);
+            break;
+        case expr_type_function:
+            evaluate_function(e);
+            break;
+        case expr_type_binding:
+            evaluate_binding(e);
+            break;
+        default:
+            G_fatal_error(_("Unknown type: %d"), e->type);
     }
 }
 
@@ -246,11 +295,11 @@ static void error_handler(void *p)
     expr_list *l;
 
     for (l = exprs; l; l = l->next) {
-	expression *e = l->exp;
-	int fd = e->data.bind.fd;
+        expression *e = l->exp;
+        int fd = e->data.bind.fd;
 
-	if (fd >= 0)
-	    unopen_output_map(fd);
+        if (fd >= 0)
+            unopen_output_map(fd);
     }
 }
 
@@ -260,42 +309,55 @@ void execute(expr_list * ee)
     expr_list *l;
     int count, n;
 
-    setup_region();
-
     exprs = ee;
     G_add_error_handler(error_handler, NULL);
 
     for (l = ee; l; l = l->next) {
-	expression *e = l->exp;
-	const char *var;
+        expression *e = l->exp;
+        const char *var;
 
-	if (e->type != expr_type_binding && e->type != expr_type_function)
-	    G_fatal_error("internal error: execute: invalid type: %d",
-			  e->type);
+        if (e->type != expr_type_binding && e->type != expr_type_function)
+            G_fatal_error("internal error: execute: invalid type: %d",
+                  e->type);
 
-	if (e->type != expr_type_binding)
-	    continue;
+        if (e->type != expr_type_binding)
+            continue;
 
-	var = e->data.bind.var;
+        var = e->data.bind.var;
 
-	if (!overwrite_flag && check_output_map(var))
-	    G_fatal_error(_("output map <%s> exists. To overwrite, use the --overwrite flag"), var);
+        if (!overwrite_flag && check_output_map(var))
+            G_fatal_error(_("output map <%s> exists. To overwrite, "
+                            "use the --overwrite flag"), var);
     }
 
+    /* Parse each expression and extract all raster maps */
     for (l = ee; l; l = l->next) {
-	expression *e = l->exp;
-	const char *var;
-	expression *val;
+        expression *e = l->exp;
+        extract_maps(e);
+    }
+
+    /* Set the region from the input maps*/
+    if (region_approach == 2)
+        prepare_region_from_maps_union(map_list, num_maps);
+    if (region_approach == 3)
+        prepare_region_from_maps_intersect(map_list, num_maps);
 
-	initialize(e);
+    setup_region();
+
+    /* Parse each expression and initialize the maps, buffers and variables */
+    for (l = ee; l; l = l->next) {
+        expression *e = l->exp;
+        const char *var;
+        expression *val;
 
-	if (e->type != expr_type_binding)
-	    continue;
+        initialize(e);
 
-	var = e->data.bind.var;
-	val = e->data.bind.val;
+        if (e->type != expr_type_binding)
+            continue;
 
-	e->data.bind.fd = open_output_map(var, val->res_type);
+        var = e->data.bind.var;
+        val = e->data.bind.val;
+        e->data.bind.fd = open_output_map(var, val->res_type);
     }
 
     setup_maps();
@@ -305,58 +367,59 @@ void execute(expr_list * ee)
 
     G_init_workers();
 
-    for (current_depth = 0; current_depth < depths; current_depth++)
-	for (current_row = 0; current_row < rows; current_row++) {
-	    if (verbose)
-		G_percent(n, count, 2);
+    for (current_depth = 0; current_depth < depths; current_depth++) {
+        for (current_row = 0; current_row < rows; current_row++) {
+            if (verbose)
+            G_percent(n, count, 2);
 
-	    for (l = ee; l; l = l->next) {
-		expression *e = l->exp;
-		int fd;
+            for (l = ee; l; l = l->next) {
+            expression *e = l->exp;
+            int fd;
 
-		evaluate(e);
+            evaluate(e);
 
-		if (e->type != expr_type_binding)
-		    continue;
+            if (e->type != expr_type_binding)
+                continue;
 
-		fd = e->data.bind.fd;
-		put_map_row(fd, e->buf, e->res_type);
-	    }
+            fd = e->data.bind.fd;
+            put_map_row(fd, e->buf, e->res_type);
+            }
 
-	    n++;
-	}
+            n++;
+        }
+    }
 
     G_finish_workers();
 
     if (verbose)
-	G_percent(n, count, 2);
+    G_percent(n, count, 2);
 
     for (l = ee; l; l = l->next) {
-	expression *e = l->exp;
-	const char *var;
-	expression *val;
-	int fd;
-
-	if (e->type != expr_type_binding)
-	    continue;
-
-	var = e->data.bind.var;
-	val = e->data.bind.val;
-	fd = e->data.bind.fd;
-
-	close_output_map(fd);
-	e->data.bind.fd = -1;
-
-	if (val->type == expr_type_map) {
-	    if (val->data.map.mod == 'M') {
-		copy_cats(var, val->data.map.idx);
-		copy_colors(var, val->data.map.idx);
-	    }
-
-	    copy_history(var, val->data.map.idx);
-	}
-	else
-	    create_history(var, val);
+        expression *e = l->exp;
+        const char *var;
+        expression *val;
+        int fd;
+
+        if (e->type != expr_type_binding)
+            continue;
+
+        var = e->data.bind.var;
+        val = e->data.bind.val;
+        fd = e->data.bind.fd;
+
+        close_output_map(fd);
+        e->data.bind.fd = -1;
+
+        if (val->type == expr_type_map) {
+            if (val->data.map.mod == 'M') {
+            copy_cats(var, val->data.map.idx);
+            copy_colors(var, val->data.map.idx);
+            }
+
+            copy_history(var, val->data.map.idx);
+        }
+        else
+        create_history(var, val);
     }
 
     G_unset_error_routine();
@@ -369,20 +432,20 @@ void describe_maps(FILE *fp, expr_list *ee)
     fprintf(fp, "output=");
 
     for (l = ee; l; l = l->next) {
-	expression *e = l->exp;
-	const char *var;
+        expression *e = l->exp;
+        const char *var;
 
-	if (e->type != expr_type_binding && e->type != expr_type_function)
-	    G_fatal_error("internal error: execute: invalid type: %d",
-			  e->type);
+        if (e->type != expr_type_binding && e->type != expr_type_function)
+            G_fatal_error("internal error: execute: invalid type: %d",
+                  e->type);
 
-	initialize(e);
+        initialize(e);
 
-	if (e->type != expr_type_binding)
-	    continue;
+        if (e->type != expr_type_binding)
+            continue;
 
-	var = e->data.bind.var;
-        fprintf(fp, "%s%s", l != ee ? "," : "", var);
+        var = e->data.bind.var;
+            fprintf(fp, "%s%s", l != ee ? "," : "", var);
     }
 
     fprintf(fp, "\n");
@@ -392,4 +455,5 @@ void describe_maps(FILE *fp, expr_list *ee)
     fprintf(fp, "\n");
 }
 
+
 /****************************************************************************/

raster/r.mapcalc/globals.h

@@ -5,8 +5,9 @@
 extern int overwrite_flag;
 extern long seed_value;
 extern long seeded;
+extern int region_approach;
 
 extern int current_depth, current_row;
-extern int depths, rows;
+extern int depths, rows, columns;
 
 #endif /* __GLOBALS_H_ */

raster/r.mapcalc/main.c

@@ -28,6 +28,7 @@ int overwrite_flag;
 
 long seed_value;
 long seeded;
+int region_approach;
 
 /****************************************************************************/
 
@@ -41,11 +42,11 @@ static expr_list *parse_file(const char *filename)
     FILE *fp;
 
     if (strcmp(filename, "-") == 0)
-	return parse_stream(stdin);
+        return parse_stream(stdin);
 
     fp = fopen(filename, "r");
     if (!fp)
-	G_fatal_error(_("Unable to open input file <%s>"), filename);
+        G_fatal_error(_("Unable to open input file <%s>"), filename);
 
     res = parse_stream(fp);
 
@@ -59,7 +60,7 @@ static expr_list *parse_file(const char *filename)
 int main(int argc, char **argv)
 {
     struct GModule *module;
-    struct Option *expr, *file, *seed;
+    struct Option *expr, *file, *seed, *region;
     struct Flag *random, *describe;
     int all_ok;
 
@@ -78,6 +79,19 @@ int main(int argc, char **argv)
     expr->description = _("Expression to evaluate");
     expr->guisection = _("Expression");
     
+    region = G_define_option();
+    region->key = "region";
+    region->type = TYPE_STRING;
+    region->required = NO;
+    region->answer = "current";
+    region->options = "current,intersect,union";
+    region->description = _("The computational region that should be used.\n"
+                            "               - current uses the current region of the mapset.\n"
+                            "               - intersect computes the intersection region between\n"
+                            "                 all input maps and uses the smallest resolution\n"
+                            "               - union computes the union extent of all map regions\n"
+                            "                 and uses the smallest resolution");
+
     file = G_define_standard_option(G_OPT_F_INPUT);
     file->key = "file";
     file->required = NO;
@@ -100,50 +114,62 @@ int main(int argc, char **argv)
 
     if (argc == 1)
     {
-	char **p = G_malloc(3 * sizeof(char *));
-	p[0] = argv[0];
-	p[1] = G_store("file=-");
-	p[2] = NULL;
-	argv = p;
-	argc = 2;
+        char **p = G_malloc(3 * sizeof(char *));
+        p[0] = argv[0];
+        p[1] = G_store("file=-");
+        p[2] = NULL;
+        argv = p;
+        argc = 2;
     }
 
     if (G_parser(argc, argv))
-	exit(EXIT_FAILURE);
+        exit(EXIT_FAILURE);
 
     overwrite_flag = module->overwrite;
 
     if (expr->answer && file->answer)
-	G_fatal_error(_("%s= and %s= are mutually exclusive"),
-			expr->key, file->key);
+        G_fatal_error(_("%s= and %s= are mutually exclusive"),
+                        expr->key, file->key);
 
     if (seed->answer && random->answer)
-	G_fatal_error(_("%s= and -%c are mutually exclusive"),
-		      seed->key, random->key);
+        G_fatal_error(_("%s= and -%c are mutually exclusive"),
+                        seed->key, random->key);
 
     if (expr->answer)
-	result = parse_string(expr->answer);
+        result = parse_string(expr->answer);
     else if (file->answer)
-	result = parse_file(file->answer);
+        result = parse_file(file->answer);
     else
-	result = parse_stream(stdin);
+        result = parse_stream(stdin);
 
     if (!result)
-	G_fatal_error(_("parse error"));
+        G_fatal_error(_("parse error"));
 
     if (seed->answer) {
-	seed_value = atol(seed->answer);
-	G_srand48(seed_value);
-	seeded = 1;
-	G_debug(3, "Read random seed from seed=: %ld", seed_value);
+        seed_value = atol(seed->answer);
+        G_srand48(seed_value);
+        seeded = 1;
+        G_debug(3, "Read random seed from seed=: %ld", seed_value);
     }
 
     if (random->answer) {
-	seed_value = G_srand48_auto();
-	seeded = 1;
-	G_debug(3, "Generated random seed (-s): %ld", seed_value);
+        seed_value = G_srand48_auto();
+        seeded = 1;
+        G_debug(3, "Generated random seed (-s): %ld", seed_value);
     }
 
+    /* Set the global variable of the region setup approach */ 
+    region_approach = 1;
+
+    if (G_strncasecmp(region->answer, "union", 5) == 0)
+        region_approach = 2;
+
+    if (G_strncasecmp(region->answer, "intersect", 9) == 0)
+        region_approach = 3;
+
+    G_debug(1, "Region answer %s region approach %i", region->answer,
+                                                      region_approach);
+    
     if (describe->answer) {
         describe_maps(stdout, result);
         return EXIT_SUCCESS;
@@ -156,8 +182,8 @@ int main(int argc, char **argv)
     all_ok = 1;
 
     if (floating_point_exception_occurred) {
-	G_warning(_("Floating point error(s) occurred in the calculation"));
-	all_ok = 0;
+        G_warning(_("Floating point error(s) occurred in the calculation"));
+        all_ok = 0;
     }
 
     return all_ok ? EXIT_SUCCESS : EXIT_FAILURE;

raster/r.mapcalc/map.c

@@ -22,6 +22,8 @@
 
 /****************************************************************************/
 
+static void prepare_region_from_maps(expression **, int, int);
+int columns;
 struct Cell_head current_region2;
 
 void setup_region(void)
@@ -777,3 +779,85 @@ void create_history(const char *dst, expression * e)
 }
 
 /****************************************************************************/
+
+void prepare_region_from_maps_union(expression **map_list, int num_maps) {
+    prepare_region_from_maps(map_list, 1, num_maps);
+}
+
+void prepare_region_from_maps_intersect(expression **map_list, int num_maps) {
+    prepare_region_from_maps(map_list, 2, num_maps);
+}
+
+void prepare_region_from_maps(expression **map_list, int type, int num_maps) {
+    /* \brief Setup the computational region from all raster maps
+     *
+     * This function computes the disjoint union or intersection extent
+     * that is covered by all raster maps and sets the smallest
+     * resolution that it can find.
+     */
+
+    int first = 0;
+    struct Cell_head window, temp_window;
+    int i;
+    const char *mapset;
+
+    for (i = 0; i < num_maps; i++) {
+        expression *e = map_list[i];
+
+        char rast_name[GNAME_MAX];
+
+        strcpy(rast_name, e->data.map.name);
+        mapset = G_find_raster2(rast_name, "");
+
+        if (!mapset)
+            G_fatal_error(_("Raster map <%s> not found"), rast_name);
+
+        G_debug(1, "Setting region from raster map: %s@%s", rast_name, mapset);
+
+        Rast_get_cellhd(rast_name, mapset, &temp_window);
+
+        if (!first) {
+            window = temp_window;
+            first = 1;
+        }
+        else {
+            if (type == 1) {
+                /* Union: find the largest extent */
+                window.north = (window.north > temp_window.north) ?
+                    window.north : temp_window.north;
+                window.south = (window.south < temp_window.south) ?
+                    window.south : temp_window.south;
+                window.east = (window.east > temp_window.east) ?
+                    window.east : temp_window.east;
+                window.west = (window.west < temp_window.west) ?
+                    window.west : temp_window.west;
+            }
+            else {
+                /* Intersect: Find the smallest extent */
+                window.north = (window.north < temp_window.north) ?
+                    window.north : temp_window.north;
+                window.south = (window.south > temp_window.south) ?
+                    window.south : temp_window.south;
+                window.east = (window.east < temp_window.east) ?
+                    window.east : temp_window.east;
+                window.west = (window.west > temp_window.west) ?
+                    window.west : temp_window.west;
+            }
+            /* Find the smallest resolution */
+            window.ns_res = (window.ns_res < temp_window.ns_res) ?
+                window.ns_res : temp_window.ns_res;
+            window.ew_res = (window.ew_res < temp_window.ew_res) ?
+                window.ew_res : temp_window.ew_res;
+        }
+    }
+
+    /* Set the region only if a map was found in the expression */
+    if (first == 1) {
+        G_adjust_Cell_head3(&window, 0, 0, 0);
+        G_debug(1, "Region was set to n %g s %g e %g w %g ns_res %g ew_res %g", window.north,
+                  window.south, window.east, window.west, window.ns_res, window.ew_res);
+        G_put_window(&window);
+    }
+
+}
+/****************************************************************************/

raster/r.mapcalc/map3.c

@@ -18,6 +18,7 @@
 
 /****************************************************************************/
 
+static void prepare_region_from_maps(expression **, int, int);
 RASTER3D_Region current_region3;
 
 void setup_region(void)
@@ -671,3 +672,90 @@ void create_history(const char *dst, expression * e)
 }
 
 /****************************************************************************/
+
+void prepare_region_from_maps_union(expression **map_list, int num_maps) {
+    prepare_region_from_maps(map_list, 1, num_maps);
+}
+
+void prepare_region_from_maps_intersect(expression **map_list, int num_maps) {
+    prepare_region_from_maps(map_list, 2, num_maps);
+}
+
+void prepare_region_from_maps(expression **map_list, int type, int num_maps) {
+   /* \brief Setup the computational region from all 3d raster maps
+    *
+    * This function computes the disjoint union or intersection extent
+    * that is covered by all 3d raster maps and sets the smallest
+    * resolution that it can find.
+    *
+    * TODO: Must be implemented
+    */
+
+    G_fatal_error("Union or intersection of regions is not implemented");
+
+    int first = 0;
+    struct Cell_head window, temp_window;
+    int i;
+    const char *mapset;
+
+    for (i = 0; i < num_maps; i++) {
+        expression *e = map_list[i];
+
+        char rast_name[GNAME_MAX];
+
+        strcpy(rast_name, e->data.map.name);
+        mapset = G_find_raster2(rast_name, "");
+
+        if (!mapset)
+            G_fatal_error(_("Raster map <%s> not found"), rast_name);
+
+        G_debug(1, "Setting region from raster map: %s@%s", rast_name, mapset);
+
+        Rast_get_cellhd(rast_name, mapset, &temp_window);
+
+        if (!first) {
+            window = temp_window;
+            first = 1;
+        }
+        else {
+            if (type == 1) {
+                /* Union: find the largest extent */
+                window.north = (window.north > temp_window.north) ?
+                    window.north : temp_window.north;
+                window.south = (window.south < temp_window.south) ?
+                    window.south : temp_window.south;
+                window.east = (window.east > temp_window.east) ?
+                    window.east : temp_window.east;
+                window.west = (window.west < temp_window.west) ?
+                    window.west : temp_window.west;
+            }
+            else {
+                /* Intersect: Find the smallest extent */
+                window.north = (window.north < temp_window.north) ?
+                    window.north : temp_window.north;
+                window.south = (window.south > temp_window.south) ?
+                    window.south : temp_window.south;
+                window.east = (window.east < temp_window.east) ?
+                    window.east : temp_window.east;
+                window.west = (window.west > temp_window.west) ?
+                    window.west : temp_window.west;
+            }
+            /* Find the smallest resolution */
+            window.ns_res = (window.ns_res < temp_window.ns_res) ?
+                window.ns_res : temp_window.ns_res;
+            window.ew_res = (window.ew_res < temp_window.ew_res) ?
+                window.ew_res : temp_window.ew_res;
+        }
+    }
+
+    /* Set the region only if a map was found in the expression */
+    if (first == 1) {
+        G_adjust_Cell_head3(&window, 0, 0, 0);
+        G_debug(1, "Region was set to n %g s %g e %g w %g ns_res %g ew_res %g", window.north,
+                  window.south, window.east, window.west, window.ns_res, window.ew_res);
+        G_put_window(&window);
+    }
+
+}
+
+/****************************************************************************/

raster/r.mapcalc/mapcalc.h

@@ -63,6 +63,9 @@ extern void copy_colors(const char *dst, int idx);
 extern void copy_history(const char *dst, int idx);
 extern void create_history(const char *dst, expression * e);
 
+extern void prepare_region_from_maps_union(expression **, int);
+extern void prepare_region_from_maps_intersect(expression **, int);
+
 /****************************************************************************/
 
 #endif /* _MAPCALC_H_ */

raster/r.mapcalc/r.mapcalc.html

@@ -59,9 +59,51 @@ EOF
 <em>result</em> map title (which appears in the category file for <em>result</em>)
 and in the history file for <em>result</em>.
 <p>Some characters have special meaning to the command shell. If the user
-is entering input to <em>r.mapcalc</em> on the command line, expressions
-should be enclosed within single quotes.  See NOTES, below.
-
+    is entering input to <em>r.mapcalc</em> on the command line, expressions
+    should be enclosed within single quotes. See NOTES, below.
+</p>
+<h3>Computational regions in r.mapcalc</h3>
+<p>
+    By default <em>r.mapcalc</em> uses the current region
+    as computational region that
+    was set with <a href="g.region.html">g.region</a> for processing.
+    Sometimes it is necessary to use a region that is derived from the
+    raster maps in the expression to set the computational region.
+    This is of high importance for modules that use r.mapcalc internally
+    to process time series of satellite images that all have different
+    spatial extents. A module that requires this feature
+    is <a href="t.rast.algebra.html">t.rast.algebra</a>.
+    The <em>region</em> option of <em>r.mapcalc</em>
+    was implemented to address this requirement.
+    It allows to compute and use a region based on
+    all raster maps in an expression. Three modes are supported:
+</p>
+<ul>
+    <li>
+        Setting the <em>region</em> parameter to <em>current</em>
+        will result in the use of the current region
+        as computational region. This is the default.
+        The current region can be set with <a href="g.region.html">g.region</a>.
+    </li>
+    <li>
+        The parameter <em>union</em> will force r.mapcalc
+        to compute the disjoint union of all regions from raster
+        maps specified in the expression. This
+        computed region will then be used as computational region at runtime.
+        The region of the mapset will not be modified.
+        The smallest spatial resolution
+        of all raster maps will be used for processing.
+    </li>
+    <li>
+        The parameter <em>intersect</em> will force r.mapcalc
+        to compute the intersection of all regions from raster
+        maps specified in the expression. This
+        computed region will then be used as computational region at runtime.
+        The region of the mapset will not be modified.
+        The smallest spatial resolution
+        of all raster maps will be used for processing.
+    </li>
+</ul>
 
 <h3>Operators and order of precedence</h3>
 
@@ -208,7 +250,7 @@ Therefore, the category label must start with a valid number.
 If the category label is integer, it will be represented by
 a floating point number. I the category label does not start
 with a number or is missing, it will be represented by NULL
-(no data) in the resulting raster map. 
+(no data) in the resulting raster map.
 
 
 <h3>Grey scale equivalents and color separates</h3>
@@ -373,8 +415,8 @@ produce a floating-point result, as will using float():
 <h3>NULL support</h3>
 
 <ul>
-<li>Division by zero should result in NULL. 
-<li>Modulus by zero should result in NULL. 
+<li>Division by zero should result in NULL.
+<li>Modulus by zero should result in NULL.
 <li>NULL-values in any arithmetic or logical operation should result
 in NULL. (however, &amp;&amp;&amp; and ||| are treated specially, as described below).
 <li>The &amp;&amp;&amp; and ||| operators observe the following axioms even when x is NULL:
@@ -389,25 +431,25 @@ if(), eval() and isnull() are treated specially, as described below).
 <li>The eval() function always returns its last argument
 <li>The situation for if() is:
 <div class="code"><pre>
-if(x) 
-	NULL if x is NULL; 0 if x is zero; 1 otherwise 
-if(x,a) 
-	NULL if x is NULL; a if x is non-zero; 0 otherwise 
-if(x,a,b) 
-	NULL if x is NULL; a if x is non-zero; b otherwise 
-if(x,n,z,p) 
-	NULL if x is NULL; n if x is negative; 
-z if x is zero; p if x is positive 
-</pre></div>
-<li>The (new) function isnull(x) returns: 1 if x is NULL; 
-0 otherwise. The (new) function null() 
-(which has no arguments) returns an integer NULL. 
-<li>Non-NULL, but invalid, arguments to functions should result in NULL. 
-<div class="code"><pre>
-Examples: 
-log(-2) 
-sqrt(-2) 
-pow(a,b) where a is negative and b is not an integer 
+if(x)
+	NULL if x is NULL; 0 if x is zero; 1 otherwise
+if(x,a)
+	NULL if x is NULL; a if x is non-zero; 0 otherwise
+if(x,a,b)
+	NULL if x is NULL; a if x is non-zero; b otherwise
+if(x,n,z,p)
+	NULL if x is NULL; n if x is negative;
+z if x is zero; p if x is positive
+</pre></div>
+<li>The (new) function isnull(x) returns: 1 if x is NULL;
+0 otherwise. The (new) function null()
+(which has no arguments) returns an integer NULL.
+<li>Non-NULL, but invalid, arguments to functions should result in NULL.
+<div class="code"><pre>
+Examples:
+log(-2)
+sqrt(-2)
+pow(a,b) where a is negative and b is not an integer
 </pre></div>
 </ul>
 <p>NULL support: Please note that any math performed with NULL cells always
@@ -510,7 +552,7 @@ r.mapcalc command. E.g. rather than:
 
 <h3>Backwards compatibility</h3>
 
-For the backwards compatibility with GRASS 6, 
+For the backwards compatibility with GRASS 6,
 <!-- check wording: -->
 if no options are given, it manufactures <tt>file=-</tt> (which reads from
 stdin), so you can continue to use e.g.:
@@ -618,7 +660,7 @@ g.rename raster=newmap,oldmap
 
 <h3>Random number generator initialization</h3>
 <p>The pseudo-random number generator used by the rand() function can
-be initialised to a specific value using the <b>seed</b> option. 
+be initialised to a specific value using the <b>seed</b> option.
 This can be used to replicate a previous calculation.
 <p>Alternatively, it can be initialised from the system time and the
 PID using the <b>-r</b> flag. This should result in a different seed
@@ -688,7 +730,7 @@ random_fcell = float(rand(-100.0,100.0))
 </pre></div>
 
 <p>
-The graph() function allows users to specify a x-y conversion using 
+The graph() function allows users to specify a x-y conversion using
 pairs of x,y coordinates.
 In some situations a transformation from one value to another is not
 easily established mathematically, but can be represented by a 2-D
@@ -770,7 +812,7 @@ NULL cells. It is left to the user to utilize the <tt>isnull()</tt> function.
 <b><a href="http://grass.osgeo.org/uploads/grass/history_docs/mapcalc-algebra.pdf">r.mapcalc: An Algebra for GIS and Image
 Processing</a></b>, by Michael Shapiro and Jim Westervelt, U.S. Army
 Construction Engineering Research Laboratory (March/1991).
-<p> 
+<p>
 <b><a href="http://grass.osgeo.org/uploads/grass/history_docs/mapcalc.pdf">Performing Map Calculations on GRASS Data:
 r.mapcalc Program Tutorial</a></b>, by Marji Larson, Michael Shapiro and Scott
 Tweddale, U.S. Army Construction Engineering Research Laboratory (December
@@ -781,7 +823,7 @@ by Anil K. Jain (Prentice Hall, NJ, 1989; p 67).
 
 
 <h2>AUTHORS</h2>
-Michael Shapiro, U.S.Army Construction Engineering 
+Michael Shapiro, U.S.Army Construction Engineering
 Research Laboratory
 <p>Glynn Clements
 

raster/r.mapcalc/testsuite/test_r3_mapcalc.py

@@ -16,7 +16,7 @@ class TestBasicOperations(TestCase):
         cls.runModule('g.region',
             n=85, s=5, e=85, w=5,
             b=0, t=2000,
-            res=1e-07)
+            res=1, res3=1)
 
     @classmethod
     def tearDownClass(cls):
@@ -51,14 +51,14 @@ class TestBasicOperations(TestCase):
             expression='diff_e_e = 3 * x() * y() * z() - 3 * x() * y() * z()')
         self.to_remove.append('diff_e_e')
         self.assertRaster3dMinMax('diff_e_e', refmin=0, refmax=0)
-    
+
     def test_nrows_ncols_ndepths_sum(self):
         """Test if sum of nrows, ncols and ndepths matches one
         expected from current region settigs"""
         self.assertModule('r3.mapcalc',
             expression='nrows_ncols_ndepths_sum = nrows() + ncols() + ndepths()')
         self.to_remove.append('nrows_ncols_ndepths_sum')
-        self.assertRasterMinMax('nrows_ncols_ndepths_sum', refmin=2160, refmax=2160)
+        self.assertRaster3dMinMax('nrows_ncols_ndepths_sum', refmin=2160, refmax=2160)
 
 
 if __name__ == '__main__':

raster/r.mapcalc/testsuite/test_r_mapcalc.py

@@ -9,16 +9,16 @@ east: 25
 west: 15
 rows: 10
 cols: 10
-121 12 183 55 37 96 138 117 182 40 
-157 70 115 1 149 125 42 193 108 24 
-83 66 82 84 186 182 179 122 67 113 
-151 93 144 173 128 196 61 125 64 193 
-180 175 14 41 44 27 165 27 90 60 
-97 57 12 104 98 13 87 24 83 107 
-174 133 146 114 115 60 78 154 49 130 
-55 138 144 25 32 58 47 137 139 32 
-143 193 155 190 131 124 87 81 160 154 
-56 45 48 66 9 182 69 12 154 19 
+121 12 183 55 37 96 138 117 182 40
+157 70 115 1 149 125 42 193 108 24
+83 66 82 84 186 182 179 122 67 113
+151 93 144 173 128 196 61 125 64 193
+180 175 14 41 44 27 165 27 90 60
+97 57 12 104 98 13 87 24 83 107
+174 133 146 114 115 60 78 154 49 130
+55 138 144 25 32 58 47 137 139 32
+143 193 155 190 131 124 87 81 160 154
+56 45 48 66 9 182 69 12 154 19
 """
 
 dcell_seed_600 = """\
@@ -93,7 +93,7 @@ class TestRandFunction(TestCase):
 
     def test_seed_required(self):
         """Test that seed is required when rand() is used
-        
+
         This test can, and probably should, generate an error message.
         """
         self.assertModuleFail('r.mapcalc', expression='rand_x = rand(1, 200)')
@@ -218,7 +218,7 @@ class TestBasicOperations(TestCase):
             expression='diff_e_e = 3 * x() * y() - 3 * x() * y()')
         self.to_remove.append('diff_e_e')
         self.assertRasterMinMax('diff_e_e', refmin=0, refmax=0)
-    
+
     def test_nrows_ncols_sum(self):
         """Test if sum of nrows and ncols matches one
         expected from current region settigs"""
@@ -228,5 +228,54 @@ class TestBasicOperations(TestCase):
         self.assertRasterMinMax('nrows_ncols_sum', refmin=20, refmax=20)
 
 
+class TestRegionOperations(TestCase):
+
+    # TODO: replace by unified handing of maps
+    to_remove = []
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule('g.region', n=30, s=15, e=30, w=15, res=5)
+        cls.runModule('r.mapcalc', expression="test_region_1 = 1", seed=1)
+        cls.runModule('g.region', n=25, s=10, e=25, w=10, res=5)
+        cls.runModule('r.mapcalc', expression="test_region_2 = 2", seed=1)
+        cls.runModule('g.region', n=20, s=5, e=20, w=5, res=1)
+        cls.runModule('r.mapcalc', expression="test_region_3 = 3", seed=1)
+
+        cls.to_remove.append("test_region_1")
+        cls.to_remove.append("test_region_2")
+        cls.to_remove.append("test_region_3")
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        if cls.to_remove:
+            cls.runModule('g.remove', flags='f', type='raster',
+                name=','.join(cls.to_remove), verbose=True)
+
+    def test_union(self):
+        """Test the union region option"""
+        self.assertModule('r.mapcalc', region="union", seed=1,
+                          expression='test_region_4 = test_region_1 + test_region_2 + test_region_3')
+        self.to_remove.append('test_region_4')
+
+        self.assertModuleKeyValue('r.info', map='test_region_4', flags='gr',
+                                  reference=dict(min=6, max=6, cells=625, north=30,
+                                  south=5, west=5, east=30, nsres=1, ewres=1),
+                                  precision=0.01, sep='=')
+
+    def test_intersect(self):
+        """Test the intersect region option"""
+        self.assertModule('r.mapcalc', region="intersect", seed=1,
+                          expression='test_region_5 = test_region_1 + test_region_2 + test_region_3')
+        self.to_remove.append('test_region_5')
+
+        self.assertModuleKeyValue('r.info', map='test_region_5', flags='gr',
+                                  reference=dict(min=6, max=6, cells=25, north=20,
+                                  south=15, west=15, east=20, nsres=1, ewres=1),
+                                  precision=0.01, sep='=')
+
+
 if __name__ == '__main__':
     test()