r3.to.rast: add multiply and add coeff options

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

raster3d/r3.to.rast/main.c

@@ -27,6 +27,8 @@
 typedef struct {
     struct Option *input, *output;
     struct Option *type;
+    struct Option *coeff_a;
+    struct Option *coeff_b;
     struct Flag *mask;
     struct Flag *res; /*If set, use the same resolution as the input map */
 } paramType;
@@ -37,7 +39,8 @@ paramType param; /*Parameters */
 void fatal_error(void *map, int *fd, int depths, char *errorMsg); /*Simple Error message */
 void set_params(); /*Fill the paramType structure */
 void g3d_to_raster(void *map, RASTER3D_Region region, int *fd,
-                   int output_type); /*Write the raster */
+                   int output_type, int use_coeffs, double coeff_a,
+                   double coeff_b); /*Write the raster */
 int open_output_map(const char *name, int res_type); /*opens the outputmap */
 void close_output_map(int fd); /*close the map */
 
@@ -106,6 +109,20 @@ void set_params()
     param.type = G_define_standard_option(G_OPT_R_TYPE);
     param.type->required = NO;
 
+    param.coeff_a = G_define_option();
+    param.coeff_a->key = "multiply";
+    param.coeff_a->type = TYPE_DOUBLE;
+    param.coeff_a->required = NO;
+    param.coeff_a->label = _("Value to multiply the raster values with");
+    param.coeff_a->description = _("Coefficient a in the equation y = ax + b");
+
+    param.coeff_b = G_define_option();
+    param.coeff_b->key = "add";
+    param.coeff_b->type = TYPE_DOUBLE;
+    param.coeff_b->required = NO;
+    param.coeff_b->label = _("Value to add to the raster values");
+    param.coeff_b->description = _("Coefficient b in the equation y = ax + b");
+
     param.mask = G_define_flag();
     param.mask->key = 'm';
     param.mask->description = _("Use 3D raster mask (if exists) with input map");
@@ -121,8 +138,11 @@ void set_params()
 /* Write the slices to seperate raster maps ******************************** */
 
 /* ************************************************************************* */
+/* coefficients are used only when needed, otherwise the original values
+ * is preserved as well as possible */
 void g3d_to_raster(void *map, RASTER3D_Region region, int *fd,
-                   int output_type)
+                   int output_type, int use_coeffs, double coeff_a,
+                   double coeff_b)
 {
     CELL c1 = 0;
     FCELL f1 = 0;
@@ -164,16 +184,24 @@ void g3d_to_raster(void *map, RASTER3D_Region region, int *fd,
             for (x = 0; x < cols; x++) {
                 if (typeIntern == FCELL_TYPE) {
                     Rast3d_get_value(map, x, y, z, &f1, typeIntern);
-                    if (Rast3d_is_null_value_num(&f1, FCELL_TYPE))
+                    if (Rast3d_is_null_value_num(&f1, FCELL_TYPE)) {
                         Rast_set_null_value(ptr, 1, output_type);
-                    else
+                    }
+                    else {
+                        if (use_coeffs)
+                            f1 = coeff_a * f1 + coeff_b;
                         Rast_set_f_value(ptr, f1, output_type);
+                    }
                 } else {
                     Rast3d_get_value(map, x, y, z, &d1, typeIntern);
-                    if (Rast3d_is_null_value_num(&d1, DCELL_TYPE))
+                    if (Rast3d_is_null_value_num(&d1, DCELL_TYPE)) {
                         Rast_set_null_value(ptr, 1, output_type);
-                    else
+                    }
+                    else {
+                        if (use_coeffs)
+                            d1 = coeff_a * d1 + coeff_b;
                         Rast_set_d_value(ptr, d1, output_type);
+                    }
                 }
                 ptr = G_incr_void_ptr(ptr, cell_size);
             }
@@ -219,6 +247,9 @@ int main(int argc, char *argv[])
     int *fd = NULL, output_type, cols, rows;
     char *RasterFileName;
     int overwrite = 0;
+    int use_coeffs = 0;  /* bool */
+    double coeff_a = 1;
+    double coeff_b = 0;
 
     /* Initialize GRASS */
     G_gisinit(argv[0]);
@@ -243,6 +274,14 @@ int main(int argc, char *argv[])
         Rast3d_fatal_error(_("3D raster map <%s> not found"),
                        param.input->answer);
 
+    /* coefficients to modify the map */
+    if (param.coeff_a->answer || param.coeff_b->answer)
+        use_coeffs = 1;
+    if (param.coeff_a->answer)
+        coeff_a = atof(param.coeff_a->answer);
+    if (param.coeff_b->answer)
+        coeff_b = atof(param.coeff_b->answer);
+
     /*Set the defaults */
     Rast3d_init_defaults();
 
@@ -313,7 +352,6 @@ int main(int argc, char *argv[])
         output_type = Rast3d_file_type_map(map);
     }
 
-
     /*prepare the filehandler */
     fd = (int *) G_malloc(region.depths * sizeof (int));
 
@@ -350,7 +388,7 @@ int main(int argc, char *argv[])
     }
 
     /*Create the Rastermaps */
-    g3d_to_raster(map, region, fd, output_type);
+    g3d_to_raster(map, region, fd, output_type, use_coeffs, coeff_a, coeff_b);
 
 
     /*Loop over all output maps! close */

raster3d/r3.to.rast/testsuite/test_a_b_coeff.py

@@ -0,0 +1,160 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE:        test_small_data
+# AUTHOR:        Vaclav Petras
+# PURPOSE:       Fast test using a small example
+# COPYRIGHT:     (C) 2016 by Vaclav Petras and 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.
+#
+#############################################################################
+
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+from grass.script import list_strings
+
+# generated by
+# g.region n=12 s=9 e=21 w=18 t=8 b=4 res=1 res3=1 -p3
+# r3.mapcalc "x = rand(0,10)" seed=100 && r3.out.ascii x prec=0
+INPUT = """\
+version: grass7
+order: nsbt
+north: 12
+south: 9
+east: 21
+west: 18
+top: 8
+bottom: 4
+rows: 3
+cols: 3
+levels: 4
+6 5 1
+0 7 5
+1 7 1
+8 2 1
+3 4 2
+8 5 6
+1 2 8
+1 5 5
+1 1 3
+1 8 3
+6 5 1
+5 1 7
+"""
+
+# created from the above and template from
+# r.mapcalc "x = rand(0,10)" seed=100 && r.out.ascii x prec=0
+OUTPUTS = [
+"""\
+north: 12
+south: 9
+east: 21
+west: 18
+rows: 3
+cols: 3
+12.5 10.5 2.5
+0.5 14.5 10.5
+2.5 14.5 2.5
+""",
+"""\
+north: 12
+south: 9
+east: 21
+west: 18
+rows: 3
+cols: 3
+16.5 4.5 2.5
+6.5 8.5 4.5
+16.5 10.5 12.5
+""",
+"""\
+north: 12
+south: 9
+east: 21
+west: 18
+rows: 3
+cols: 3
+2.5 4.5 16.5
+2.5 10.5 10.5
+2.5 2.5 6.5
+""",
+"""\
+north: 12
+south: 9
+east: 21
+west: 18
+rows: 3
+cols: 3
+2.5 16.5 6.5
+12.5 10.5 2.5
+10.5 2.5 14.5
+""",
+]
+
+
+class TestR3ToRast(TestCase):
+    # TODO: replace by unified handing of maps
+    # mixing class and object attributes
+    to_remove_3d = []
+    to_remove_2d = []
+    rast3d = 'r3_to_rast_test_a_b_coeff'
+    rast2d = 'r3_to_rast_test_a_b_coeff'
+    rast2d_ref = 'r3_to_rast_test_a_b_coeff_ref'
+    rast2d_refs = []
+
+    def setUp(self):
+        self.use_temp_region()
+        self.runModule('r3.in.ascii', input='-', stdin_=INPUT,
+                       output=self.rast3d)
+        self.to_remove_3d.append(self.rast3d)
+        self.runModule('g.region', raster_3d=self.rast3d)
+
+        for i, data in enumerate(OUTPUTS):
+            rast = "%s_%d" % (self.rast2d_ref, i)
+            self.runModule('r.in.ascii', input='-', stdin_=data,
+                           output=rast)
+            self.to_remove_2d.append(rast)
+            self.rast2d_refs.append(rast)
+
+    def tearDown(self):
+        if self.to_remove_3d:
+            self.runModule('g.remove', flags='f', type='raster_3d',
+                           name=','.join(self.to_remove_3d), verbose=True)
+        if self.to_remove_2d:
+            self.runModule('g.remove', flags='f', type='raster',
+                           name=','.join(self.to_remove_2d), verbose=True)
+        self.del_temp_region()
+
+    def test_a_b_coeff(self):
+        self.assertModule('r3.to.rast', input=self.rast3d,
+                          output=self.rast2d, multiply=2, add=0.5)
+        rasts = list_strings('raster', mapset=".",
+                             pattern="%s_*" % self.rast2d,
+                             exclude="%s_*" % self.rast2d_ref)
+        self.assertEquals(len(rasts), 4,
+                          msg="Wrong number of 2D rasters present"
+                              " in the mapset")
+        ref_info = dict(cells=9)
+        ref_univar = dict(cells=9, null_cells=0)
+        for rast in rasts:
+            self.assertRasterExists(rast)
+            # the following doesn't make much sense because we just listed them
+            self.to_remove_2d.append(rast)
+            self.assertRasterFitsInfo(raster=rast, reference=ref_info,
+                                      precision=0)
+            self.assertRasterFitsUnivar(raster=rast, reference=ref_univar,
+                                        precision=0)
+
+        # check the actual values
+        for rast_ref, rast in zip(self.rast2d_refs, rasts):
+            self.assertRastersNoDifference(actual=rast,
+                                           reference=rast_ref,
+                                           precision=0.1)
+
+
+if __name__ == '__main__':
+    test()