r.colors: add offset and scale (#954)

adding offset and scale enables r.colors to add e.g. the celsius rules to a map with units = Kelvin * 50

raster/r.colors/edit_colors.c

@@ -59,13 +59,16 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
     struct GModule *module;
     struct maps_info input_maps;
     int stat = -1;
+    int do_scale, rule_is_percent;
+    double scale, offset;
 
     struct {
         struct Flag *r, *w, *l, *d, *g, *a, *n, *e;
     } flag; 
 
     struct {
-        struct Option *maps, *colr, *rast, *volume, *rules, *file;
+        struct Option *maps, *colr, *rast, *volume, *rules, *file,
+		      *scale, *offset;
     } opt;
     
     G_gisinit(argv[0]);
@@ -127,6 +130,24 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
     opt.rules->description = _("\"-\" to read rules from stdin");
     opt.rules->guisection = _("Define");
 
+    opt.offset = G_define_option();
+    opt.offset->key = "offset";
+    opt.offset->type = TYPE_DOUBLE;
+    opt.offset->required = NO;
+    opt.offset->answer = "0";
+    opt.offset->label = _("Offset for color rule values");
+    opt.offset->description = _("New value = (old value + offset) * scale");
+    opt.offset->guisection = _("Define");
+
+    opt.scale = G_define_option();
+    opt.scale->key = "scale";
+    opt.scale->type = TYPE_DOUBLE;
+    opt.scale->required = NO;
+    opt.scale->answer = "1";
+    opt.scale->label = _("Scale for color rule values");
+    opt.scale->description = _("New value = (old value + offset) * scale");
+    opt.scale->guisection = _("Define");
+
     flag.r = G_define_flag();
     flag.r->key = 'r';
     flag.r->description = _("Remove existing color table");
@@ -199,6 +220,9 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
     rules = opt.rules->answer;
     file = opt.file->answer;
 
+    offset = atof(opt.offset->answer);
+    scale = atof(opt.scale->answer);
+
     if (opt.rast->answer)
         cmap = opt.rast->answer;
     if (opt.volume->answer)
@@ -360,7 +384,7 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
 	}
 
 	if (i > 0) {
-	    if(has_fcell_type && has_cell_type) {
+	    if (has_fcell_type && has_cell_type) {
 		G_fatal_error("Input maps must have the same cell type. "
 				"Mixing of integer and floating point maps is not supported.");
 	    }
@@ -381,9 +405,12 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
 	}
     }
 
+    rule_is_percent = 0;
+    do_scale = 0;
     if (is_from_stdin) {
-        if (!read_color_rules(stdin, &colors, min, max, has_fcell_type))
+        if (!read_color_rules(stdin, &colors, min, max, has_fcell_type, &rule_is_percent))
             exit(EXIT_FAILURE);
+	do_scale = 1;
     }
     else if (style) {
         /* 
@@ -409,13 +436,27 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
             Rast_make_histogram_log_colors(&colors, &statf, (CELL) min,
                                            (CELL) max);
         }
-	else if (G_find_color_rule(style))
+	else if (G_find_color_rule(style)) {
+	    char path[GPATH_MAX];
+
             Rast_make_fp_colors(&colors, style, min, max);
+
+	    /* check if this style is a percentage style */
+            /* don't bother with native dirsep as not needed for backwards compatibility */
+            G_snprintf(path, GPATH_MAX, "%s/etc/colors/%s", G_gisbase(), style);
+	    rule_is_percent = check_percent_rule(path);
+	    do_scale = 1;
+	}
         else
             G_fatal_error(_("Unknown color request '%s'"), style);
     }
     else if (rules) {
-        if (!Rast_load_fp_colors(&colors, rules, min, max)) {
+	do_scale = 1;
+	/* check if these rules are percentage rules */
+        if (Rast_load_fp_colors(&colors, rules, min, max)) {
+	    rule_is_percent = check_percent_rule(rules);
+	}
+	else {
             /* for backwards compatibility try as std name; remove for GRASS 7 */
             char path[GPATH_MAX];
 
@@ -424,6 +465,7 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
 
             if (!Rast_load_fp_colors(&colors, path, min, max))
                 G_fatal_error(_("Unable to load rules file <%s>"), rules);
+	    rule_is_percent = check_percent_rule(path);
         }
     }
     else {
@@ -449,6 +491,11 @@ int edit_colors(int argc, char **argv, int type, const char *maptype,
     if (has_fcell_type)
         Rast_mark_colors_as_fp(&colors);
 
+    if (do_scale && !rule_is_percent && (offset != 0 || scale != 1)) {
+	rescale_colors(&colors_tmp, &colors, offset, scale);
+	colors = colors_tmp;
+    }
+
     if (flag.n->answer)
         Rast_invert_colors(&colors);
 

raster/r.colors/local_proto.h

@@ -47,6 +47,8 @@ void get_fp_stats(struct maps_info *, struct FP_stats *statf,
 int edit_colors(int, char **, int, const char *, const char*);
 
 /* rules.c */
-int read_color_rules(FILE *, struct Colors *, DCELL, DCELL, int);
+int read_color_rules(FILE *, struct Colors *, DCELL, DCELL, int, int *);
+int check_percent_rule(const char *);
+void rescale_colors(struct Colors *, struct Colors *, double, double);
 
 #endif /* __LOCAL_PROTO_H__ */

raster/r.colors/r.colors.html

@@ -38,6 +38,13 @@ solution.
 
 <p>The <b>-e</b> and <b>-g</b> flags are not mutually exclusive.
 
+<p><b>offset</b> and <b>scale</b> modify color rules to match the
+values of a raster map using the formula <em>new_value = (old_value +
+offset) x scale</em>. For example, if the units of a raster map are
+Kelvin x 50, the color rules <em>celsius</em> can be applied with
+<em>offset=273.15</em> and <em>scale=50</em>.
+
+
 <p>If the user specifies the <b>-w</b> flag, the current color table file for
 the input map will not be overwritten. This means that the color table is
 created only if the <i>map</i> does not already have a color table. If this

raster/r.colors/rules.c

@@ -30,8 +30,10 @@ static int read_rule(void *, DCELL, DCELL, DCELL *, int *, int *, int *,
 static void badrule(int, const char *, int);
 static int show_colors(FILE *);
 
+static int rule_is_percent = 0;
+
 int read_color_rules(FILE * fp, struct Colors *colors, DCELL min, DCELL max,
-		     int is_fp)
+		     int is_fp, int *is_percent)
 {
     DCELL rulemin, rulemax;
 
@@ -64,9 +66,78 @@ int read_color_rules(FILE * fp, struct Colors *colors, DCELL min, DCELL max,
 	G_warning(_("Your color rules do not cover the whole range of data!\n (rules %f to %f but data %f to %f)"),
 		  rulemin, rulemax, min, max);
 
+    *is_percent = rule_is_percent;
+
     return 1;
 }
 
+int check_percent_rule(const char *path)
+{
+    FILE *fp;
+    char buf[1024];
+
+    fp = fopen(path, "r");
+    if (!fp)
+	G_fatal_error(_("Unable to open color rule"));
+
+    rule_is_percent = 0;
+    while (G_getl2(buf, sizeof(buf), fp)) {
+	char value[80], color[80];
+	double x;
+	char c;
+
+	G_strip(buf);
+
+	if (*buf == '\0')
+	    continue;
+	if (*buf == '#')
+	    continue;
+
+	if (sscanf(buf, "%s %[^\n]", value, color) != 2)
+	    continue;
+
+	if (sscanf(value, "%lf%c", &x, &c) == 2 && c == '%') {
+	    rule_is_percent = 1;
+	    break;
+	}
+    }
+    fclose(fp);
+
+    return rule_is_percent;
+}
+
+void rescale_colors(struct Colors *colors_tmp, struct Colors *colors,
+                    double offset, double scale)
+{
+    int i, rcount;
+    unsigned char r1, g1, b1, r2, g2, b2;
+    int red, grn, blu;
+    DCELL dmin, dmax;
+
+    Rast_init_colors(colors_tmp);
+
+    Rast_get_default_color(&red, &grn, &blu, colors);
+    Rast_set_default_color(red, grn, blu, colors_tmp);
+
+    Rast_get_null_value_color(&red, &grn, &blu, colors);
+    Rast_set_null_value_color(red, grn, blu, colors_tmp);
+
+    rcount = Rast_colors_count(colors);
+    for (i = 0; i < rcount; i++) {
+
+	Rast_get_fp_color_rule(&dmin, &r1, &g1, &b1,
+			       &dmax, &r2, &g2, &b2,
+			       colors, i);
+
+	dmin = (dmin + offset) * scale;
+	dmax = (dmax + offset) * scale;
+
+	Rast_add_d_color_rule(&dmin, r1, g1, b1,
+			      &dmax, r2, g2, b2, colors_tmp);
+
+    }
+}
+
 static int read_rule(void *closure, DCELL min, DCELL max,
 		     DCELL * val, int *r, int *g, int *b,
 		     int *norm, int *nval, int *dflt)
@@ -79,6 +150,9 @@ static int read_rule(void *closure, DCELL min, DCELL max,
     for (;;) {
 	char buf[1024];
 	int ret, i;
+	char value[80], color[80];
+	double x;
+	char c;
 
 	if (tty)
 	    fprintf(stderr, "> ");
@@ -113,6 +187,12 @@ static int read_rule(void *closure, DCELL min, DCELL max,
 	    continue;
 	}
 
+	if (sscanf(buf, "%s %[^\n]", value, color) == 2) {
+	    if (sscanf(value, "%lf%c", &x, &c) == 2 && c == '%') {
+		rule_is_percent = 1;
+	    }
+	}
+
 	ret =
 	    Rast_parse_color_rule(min, max, buf, val, r, g, b, norm, nval, dflt);
 	if (ret == 0)

raster/r.colors/stats.c

@@ -21,7 +21,8 @@
 #include <stdlib.h>
 #include "local_proto.h"
 
-int get_stats(struct maps_info *input_maps, struct Cell_stats *statf) {
+int get_stats(struct maps_info *input_maps, struct Cell_stats *statf)
+{
     CELL *cell;
     int row, nrows, ncols;
     int fd;
@@ -29,25 +30,25 @@ int get_stats(struct maps_info *input_maps, struct Cell_stats *statf) {
 
     Rast_init_cell_stats(statf);
 
-    for(i = 0; i < input_maps->num; i++) {
-		fd = Rast_open_old(input_maps->names[i], input_maps->mapsets[i]);
+    for (i = 0; i < input_maps->num; i++) {
+	fd = Rast_open_old(input_maps->names[i], input_maps->mapsets[i]);
 
-	    cell = Rast_allocate_c_buf();
-	    nrows = Rast_window_rows();
-	    ncols = Rast_window_cols();
+	cell = Rast_allocate_c_buf();
+	nrows = Rast_window_rows();
+	ncols = Rast_window_cols();
 
-		G_verbose_message(_("(%i/%i) Reading raster map <%s>..."),
-				i + 1, input_maps->num,
-				G_fully_qualified_name(input_maps->names[i], input_maps->mapsets[i]));
+	G_verbose_message(_("(%i/%i) Reading raster map <%s>..."),
+			i + 1, input_maps->num,
+			G_fully_qualified_name(input_maps->names[i], input_maps->mapsets[i]));
 
-		for (row = 0; row < nrows; row++) {
-			G_percent(row, nrows, 2);
-			Rast_get_c_row(fd, cell, row);
-			Rast_update_cell_stats(cell, ncols, statf);
-		}
-		G_percent(row, nrows, 2);
-		Rast_close(fd);
-	    G_free(cell);
+	for (row = 0; row < nrows; row++) {
+	    G_percent(row, nrows, 2);
+	    Rast_get_c_row(fd, cell, row);
+	    Rast_update_cell_stats(cell, ncols, statf);
+	}
+	G_percent(nrows, nrows, 2);
+	Rast_close(fd);
+	G_free(cell);
     }
 
     return 1;
@@ -55,7 +56,8 @@ int get_stats(struct maps_info *input_maps, struct Cell_stats *statf) {
 
 void get_fp_stats(struct maps_info *input_maps,
         struct FP_stats *statf,
-        DCELL min, DCELL max, int geometric, int geom_abs, int type) {
+        DCELL min, DCELL max, int geometric, int geom_abs, int type)
+{
     DCELL *dcell = NULL;
     int row, col, depth, nrows, ncols, ndepths = 1;
     int fd;
@@ -64,108 +66,113 @@ void get_fp_stats(struct maps_info *input_maps,
     char *mapset;
     RASTER3D_Map *map3d = NULL;
 
-	statf->geometric = geometric;
-	statf->geom_abs = geom_abs;
-	statf->flip = 0;
+    statf->geometric = geometric;
+    statf->geom_abs = geom_abs;
+    statf->flip = 0;
 
-	if (statf->geometric) {
-		if (min * max < 0)
-			G_fatal_error(_("Unable to use logarithmic scaling if range includes zero"));
+    if (statf->geometric) {
+	if (min * max < 0)
+	    G_fatal_error(_("Unable to use logarithmic scaling if range includes zero"));
 
-		if (min < 0) {
-			statf->flip = 1;
-			min = -min;
-			max = -max;
-		}
-		min = log(min);
-		max = log(max);
+	if (min < 0) {
+	    statf->flip = 1;
+	    min = -min;
+	    max = -max;
 	}
+	min = log(min);
+	max = log(max);
+    }
+
+    if (statf->geom_abs) {
+	double a = log(fabs(min) + 1);
+	double b = log(fabs(max) + 1);
+	int has_zero = min * max < 0;
+
+	min = a < b ? a : b;
+	max = a > b ? a : b;
+	if (has_zero)
+	    min = 0;
+    }
+
+    statf->count = 1000;
+    statf->min = min;
+    statf->max = max;
+    statf->stats = G_calloc(statf->count + 1, sizeof (unsigned long));
+    statf->total = 0;
+
+    /* Loop over all input maps */
+    for(i = 0; i < input_maps->num; i++) {
+	name = input_maps->names[i];
+	mapset = input_maps->mapsets[i];
 
-	if (statf->geom_abs) {
-		double a = log(fabs(min) + 1);
-		double b = log(fabs(max) + 1);
-		int has_zero = min * max < 0;
-		min = a < b ? a : b;
-		max = a > b ? a : b;
-		if (has_zero)
-			min = 0;
+	fd = -1;
+	if (type == RASTER_TYPE) {
+	    fd = Rast_open_old(name, mapset);
+	    dcell = Rast_allocate_d_buf();
+	    nrows = Rast_window_rows();
+	    ncols = Rast_window_cols();
 	}
+	else {
+	    /* Initiate the default settings */
+	    Rast3d_init_defaults();
 
-	statf->count = 1000;
-	statf->min = min;
-	statf->max = max;
-	statf->stats = G_calloc(statf->count + 1, sizeof (unsigned long));
-	statf->total = 0;
-
-	/* Loop over all input maps */
-	for(i = 0; i < input_maps->num; i++) {
-		name = input_maps->names[i];
-		mapset = input_maps->mapsets[i];
-
-		if (type == RASTER_TYPE) {
-			fd = Rast_open_old(name, mapset);
-			dcell = Rast_allocate_d_buf();
-			nrows = Rast_window_rows();
-			ncols = Rast_window_cols();
-		} else {
-			/* Initiate the default settings */
-			Rast3d_init_defaults();
-
-			map3d = Rast3d_open_cell_old(name, mapset, RASTER3D_DEFAULT_WINDOW,
-					RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
-
-			if (map3d == NULL)
-				Rast3d_fatal_error(_("Error opening 3d raster map"));
-
-			nrows = map3d->window.rows;
-			ncols = map3d->window.cols;
-			ndepths = map3d->window.depths;
-		}
+	    map3d = Rast3d_open_cell_old(name, mapset, RASTER3D_DEFAULT_WINDOW,
+			    RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
 
-		G_verbose_message(_("(%i/%i) Reading map <%s>..."), i, input_maps->num,
-				G_fully_qualified_name(name, mapset));
+	    if (map3d == NULL)
+		Rast3d_fatal_error(_("Error opening 3d raster map"));
 
-		for (depth = 0; depth < ndepths; depth++) {
-			for (row = 0; row < nrows; row++) {
-				G_percent(row, nrows, 2);
+	    nrows = map3d->window.rows;
+	    ncols = map3d->window.cols;
+	    ndepths = map3d->window.depths;
+	}
 
-				if (type == RASTER_TYPE)
-					Rast_get_d_row(fd, dcell, row);
+	G_verbose_message(_("(%i/%i) Reading map <%s>..."), i, input_maps->num,
+			G_fully_qualified_name(name, mapset));
 
-				for (col = 0; col < ncols; col++) {
-					DCELL x;
-					int j;
+	for (depth = 0; depth < ndepths; depth++) {
+	    for (row = 0; row < nrows; row++) {
+		G_percent(row, nrows, 2);
 
-					if (type == RASTER_TYPE)
-						x = dcell[col];
-					else
-						x = Rast3d_get_double(map3d, col, row, depth);
+		if (type == RASTER_TYPE)
+		    Rast_get_d_row(fd, dcell, row);
 
-					if (Rast_is_d_null_value(&x))
-						continue;
+		for (col = 0; col < ncols; col++) {
+		    DCELL x;
+		    int j;
 
-					if (statf->flip)
-						x = -x;
-					if (statf->geometric)
-						x = log(x);
-					if (statf->geom_abs)
-						x = log(fabs(x) + 1);
+		    if (type == RASTER_TYPE)
+			x = dcell[col];
+		    else
+			x = Rast3d_get_double(map3d, col, row, depth);
 
-					j = (int) floor(statf->count * (x - statf->min) / (statf->max - statf->min));
-					statf->stats[j]++;
-					statf->total++;
-				}
-			}
-		}
+		    if (Rast_is_d_null_value(&x))
+			continue;
 
-		G_percent(row, nrows, 2);
+		    if (statf->flip)
+			x = -x;
+		    if (statf->geometric)
+			x = log(x);
+		    if (statf->geom_abs)
+			x = log(fabs(x) + 1);
 
-		if (type == RASTER_TYPE) {
-			Rast_close(fd);
-		if(dcell)
-				G_free(dcell);
-		} else {
-			Rast3d_close(map3d);
+		    j = (int) floor(statf->count * (x - statf->min) / (statf->max - statf->min));
+		    statf->stats[j]++;
+		    statf->total++;
 		}
+	    }
+	}
+
+	G_percent(nrows, nrows, 2);
+
+	if (type == RASTER_TYPE) {
+	    Rast_close(fd);
+	if(dcell)
+	    G_free(dcell);
+	}
+	else {
+	    Rast3d_close(map3d);
+	    map3d = NULL;
 	}
+    }
 }