[GRASS-SVN] r71360 - in grass/trunk/raster/r.mapcalc: . testsuite
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 9 23:33:42 PDT 2017
Author: huhabla
Date: 2017-08-09 23:33:41 -0700 (Wed, 09 Aug 2017)
New Revision: 71360
Modified:
grass/trunk/raster/r.mapcalc/evaluate.c
grass/trunk/raster/r.mapcalc/globals.h
grass/trunk/raster/r.mapcalc/main.c
grass/trunk/raster/r.mapcalc/map.c
grass/trunk/raster/r.mapcalc/map3.c
grass/trunk/raster/r.mapcalc/mapcalc.h
grass/trunk/raster/r.mapcalc/r.mapcalc.html
grass/trunk/raster/r.mapcalc/testsuite/test_r3_mapcalc.py
grass/trunk/raster/r.mapcalc/testsuite/test_r_mapcalc.py
Log:
raster modules: Implemented computational region settings based on disjoint union and intersection computation from all input raster maps in an expression
Modified: grass/trunk/raster/r.mapcalc/evaluate.c
===================================================================
--- grass/trunk/raster/r.mapcalc/evaluate.c 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/evaluate.c 2017-08-10 06:33:41 UTC (rev 71360)
@@ -1,6 +1,7 @@
#include <stdlib.h>
#include <unistd.h>
+#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
@@ -15,13 +16,59 @@
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;
+ }
+}
+
+/****************************************************************************/
+
static void allocate_buf(expression * e)
{
e->buf = G_malloc(columns * Rast_cell_size(e->res_type));
@@ -47,8 +94,9 @@
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 @@
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(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 @@
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_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(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 @@
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 @@
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);
+ }
- initialize(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);
- if (e->type != expr_type_binding)
- continue;
+ setup_region();
- var = e->data.bind.var;
- val = e->data.bind.val;
+ /* 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;
- e->data.bind.fd = open_output_map(var, val->res_type);
+ initialize(e);
+
+ if (e->type != expr_type_binding)
+ continue;
+
+ 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 @@
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;
+ expression *e = l->exp;
+ const char *var;
+ expression *val;
+ int fd;
- if (e->type != expr_type_binding)
- continue;
+ if (e->type != expr_type_binding)
+ continue;
- var = e->data.bind.var;
- val = e->data.bind.val;
- fd = e->data.bind.fd;
+ var = e->data.bind.var;
+ val = e->data.bind.val;
+ fd = e->data.bind.fd;
- close_output_map(fd);
- e->data.bind.fd = -1;
+ 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);
- }
+ 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);
+ copy_history(var, val->data.map.idx);
+ }
+ else
+ create_history(var, val);
}
G_unset_error_routine();
@@ -369,20 +432,20 @@
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 @@
fprintf(fp, "\n");
}
+
/****************************************************************************/
Modified: grass/trunk/raster/r.mapcalc/globals.h
===================================================================
--- grass/trunk/raster/r.mapcalc/globals.h 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/globals.h 2017-08-10 06:33:41 UTC (rev 71360)
@@ -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_ */
Modified: grass/trunk/raster/r.mapcalc/main.c
===================================================================
--- grass/trunk/raster/r.mapcalc/main.c 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/main.c 2017-08-10 06:33:41 UTC (rev 71360)
@@ -28,6 +28,7 @@
long seed_value;
long seeded;
+int region_approach;
/****************************************************************************/
@@ -41,11 +42,11 @@
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 @@
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 @@
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 @@
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 @@
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;
Modified: grass/trunk/raster/r.mapcalc/map.c
===================================================================
--- grass/trunk/raster/r.mapcalc/map.c 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/map.c 2017-08-10 06:33:41 UTC (rev 71360)
@@ -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 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);
+ }
+
+}
+/****************************************************************************/
Modified: grass/trunk/raster/r.mapcalc/map3.c
===================================================================
--- grass/trunk/raster/r.mapcalc/map3.c 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/map3.c 2017-08-10 06:33:41 UTC (rev 71360)
@@ -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 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);
+ }
+
+}
+
+/****************************************************************************/
Modified: grass/trunk/raster/r.mapcalc/mapcalc.h
===================================================================
--- grass/trunk/raster/r.mapcalc/mapcalc.h 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/mapcalc.h 2017-08-10 06:33:41 UTC (rev 71360)
@@ -63,6 +63,9 @@
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_ */
Modified: grass/trunk/raster/r.mapcalc/r.mapcalc.html
===================================================================
--- grass/trunk/raster/r.mapcalc/r.mapcalc.html 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/r.mapcalc.html 2017-08-10 06:33:41 UTC (rev 71360)
@@ -59,10 +59,52 @@
<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>
The following operators are supported:
@@ -208,7 +250,7 @@
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 @@
<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, &&& and ||| are treated specially, as described below).
<li>The &&& and ||| operators observe the following axioms even when x is NULL:
@@ -389,25 +431,25 @@
<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
+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.
+<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
+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 @@
<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 @@
<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 @@
</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 @@
<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 @@
<h2>AUTHORS</h2>
-Michael Shapiro, U.S.Army Construction Engineering
+Michael Shapiro, U.S.Army Construction Engineering
Research Laboratory
<p>Glynn Clements
Modified: grass/trunk/raster/r.mapcalc/testsuite/test_r3_mapcalc.py
===================================================================
--- grass/trunk/raster/r.mapcalc/testsuite/test_r3_mapcalc.py 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/testsuite/test_r3_mapcalc.py 2017-08-10 06:33:41 UTC (rev 71360)
@@ -16,7 +16,7 @@
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 @@
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__':
Modified: grass/trunk/raster/r.mapcalc/testsuite/test_r_mapcalc.py
===================================================================
--- grass/trunk/raster/r.mapcalc/testsuite/test_r_mapcalc.py 2017-08-09 17:40:37 UTC (rev 71359)
+++ grass/trunk/raster/r.mapcalc/testsuite/test_r_mapcalc.py 2017-08-10 06:33:41 UTC (rev 71360)
@@ -9,16 +9,16 @@
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 @@
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 @@
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 @@
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()
More information about the grass-commit
mailing list