Modified the cross section algorithm to use north, east and top to get

the volume data. Implemented a test script and added reference data.


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

raster3d/r3.cross.rast/main.c

@@ -1,45 +1,44 @@
 
 /****************************************************************************
-*
-* MODULE:       r3.cross.rast 
-*   	    	
-* AUTHOR(S):    Original author 
-*               Soeren Gebbert soerengebbert at gmx de
-* 		23 Feb 2006 Berlin
-* PURPOSE:      Creates a cross section 2D map from one G3D raster map based on a 2D elevation map  
-*
-* COPYRIGHT:    (C) 2005 by the GRASS Development Team
-*
-*               This program is free software under the GNU General Public
-*   	    	License (>=v2). Read the file COPYING that comes with GRASS
-*   	    	for details.
-*
-*****************************************************************************/
+ *
+ * MODULE:       r3.cross.rast 
+ *   	    	
+ * AUTHOR(S):    Original author 
+ *               Soeren Gebbert soerengebbert at gmx de
+ * 		23 Feb 2006 Berlin
+ * PURPOSE:      Creates a cross section 2D map from one G3D raster map based on a 2D elevation map  
+ *
+ * COPYRIGHT:    (C) 2005 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *   	    	License (>=v2). Read the file COPYING that comes with GRASS
+ *   	    	for details.
+ *
+ *****************************************************************************/
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/G3d.h>
 #include <grass/glocale.h>
 
-
 /*- param and global variables -----------------------------------------*/
-typedef struct
-{
+typedef struct {
     struct Option *input, *output, *elevation;
     struct Flag *mask;
 } paramType;
 
-paramType param;		/*param */
+paramType param; /*param */
 int globalElevMapType;
 
 
 /*- prototypes --------------------------------------------------------------*/
-void fatal_error(void *map, int elevfd, int outfd, char *errorMsg);	/*Simple Error message */
-void set_params();		/*Fill the paramType structure */
-void rast3d_cross_section(void *map, G3D_Region region, int elevfd, int outfd);	/*Write the raster */
-void close_output_map(int fd);	/*close the map */
+void fatal_error(void *map, int elevfd, int outfd, char *errorMsg); /*Simple Error message */
+void set_params(); /*Fill the paramType structure */
+void rast3d_cross_section(void *map, G3D_Region region, int elevfd, int outfd); /*Write the raster */
+void close_output_map(int fd); /*close the map */
 
 
 
@@ -52,16 +51,16 @@ void fatal_error(void *map, int elevfd, int outfd, char *errorMsg)
     /* Close files and exit */
 
     if (map != NULL) {
-	if (!G3d_closeCell(map))
-	    G3d_fatalError(_("Could not close G3D map"));
+        if (!G3d_closeCell(map))
+            G3d_fatalError(_("Could not close G3D map"));
     }
 
     /*unopen the output map */
     if (outfd != -1)
-	Rast_unopen(outfd);
+        Rast_unopen(outfd);
 
     if (elevfd != -1)
-	close_output_map(elevfd);
+        close_output_map(elevfd);
 
     G3d_fatalError(errorMsg);
     exit(EXIT_FAILURE);
@@ -95,7 +94,7 @@ void set_params()
     param.elevation->type = TYPE_STRING;
     param.elevation->required = YES;
     param.elevation->description =
-	_("2D elevation map used to create the cross section map");
+        _("2D elevation map used to create the cross section map");
     param.elevation->gisprompt = "old,cell,raster";
 
     param.output = G_define_option();
@@ -113,12 +112,11 @@ void set_params()
 
 
 /* ************************************************************************* */
-/* Calulates the resulting raster map ************************************** */
+/* Compute the cross section raster map ************************************ */
 /* ************************************************************************* */
-void rast3d_cross_section(void *map, G3D_Region region, int elevfd, int outfd)
+void rast3d_cross_section(void *map,G3D_Region region, int elevfd, int outfd)
 {
-    double d1 = 0, f1 = 0;
-    int x, y, z;
+    int col, row;
     int rows, cols, depths, typeIntern;
     FCELL *fcell = NULL;
     DCELL *dcell = NULL;
@@ -127,133 +125,100 @@ void rast3d_cross_section(void *map, G3D_Region region, int elevfd, int outfd)
     int intvalue;
     float fvalue;
     double dvalue;
-    int isnull;
     double elevation = 0;
-    double top, tbres, bottom;
-
-    /*shoter names ;) */
+    double north, east;
+    struct Cell_head window;
+ 
+    Rast_get_window(&window);
+    
     rows = region.rows;
     cols = region.cols;
     depths = region.depths;
-    top = region.top;
-    bottom = region.bottom;
-
-    /*Calculate the top-bottom resolution */
-    tbres = (top - bottom) / depths;
-
+    
     /*Typ of the G3D Tile */
     typeIntern = G3d_tileTypeMap(map);
 
     /*Allocate mem for the output maps row */
     if (typeIntern == FCELL_TYPE)
-	fcell = Rast_allocate_f_buf();
+        fcell = Rast_allocate_f_buf();
     else if (typeIntern == DCELL_TYPE)
-	dcell = Rast_allocate_d_buf();
+        dcell = Rast_allocate_d_buf();
 
     /*Mem for the input map row */
     elevrast = Rast_allocate_buf(globalElevMapType);
 
-    for (y = 0; y < rows; y++) {
-	G_percent(y, rows - 1, 10);
-
-	/*Read the input map row */
-	Rast_get_row(elevfd, elevrast, y, globalElevMapType);
-
-	for (x = 0, ptr = elevrast; x < cols; x++, ptr =
-	     G_incr_void_ptr(ptr, Rast_cell_size(globalElevMapType))) {
-
-	    /*we guess the elevation input map has no null values */
-	    isnull = 0;
-
-	    if (Rast_is_null_value(ptr, globalElevMapType)) {
-		isnull = 1;	/*input map has nulls */
-	    }
-
-	    /*Read the elevation value */
-	    if (globalElevMapType == CELL_TYPE) {
-		intvalue = *(CELL *) ptr;
-		elevation = intvalue;
-	    }
-	    else if (globalElevMapType == FCELL_TYPE) {
-		fvalue = *(FCELL *) ptr;
-		elevation = fvalue;
-	    }
-	    else if (globalElevMapType == DCELL_TYPE) {
-		dvalue = *(DCELL *) ptr;
-		elevation = dvalue;
-	    }
-
-	    /*Check if the elevation is located in the 3d raster map */
-	    if (isnull == 0) {
-		/*we guess, we have no hit */
-		isnull = 1;
-		for (z = 0; z < depths; z++) {	/*From the bottom to the top */
-		    if (elevation >= z * tbres + bottom && elevation <= (z + 1) * tbres + bottom) {	/*if at the border, choose the value from the top */
-			/*Read the value and put it in the output map row */
-            /* Because we read raster rows from north to south, but the coordinate system
-             of the g3d cube read from south to north we need to adjust the
-             Cube coordinates row = rows - y - 1.
-             */
-			if (typeIntern == FCELL_TYPE) {
-			    G3d_getValue(map, x, rows - y - 1, z, &f1, typeIntern);
-			    if (G3d_isNullValueNum(&f1, FCELL_TYPE))
-				Rast_set_null_value(&fcell[x], 1, FCELL_TYPE);
-			    else
-				fcell[x] = (FCELL) f1;
-			}
-			else {
-			    G3d_getValue(map, x, rows - y - 1, z, &d1, typeIntern);
-			    if (G3d_isNullValueNum(&d1, DCELL_TYPE))
-				Rast_set_null_value(&dcell[x], 1, DCELL_TYPE);
-			    else
-				dcell[x] = (DCELL) d1;
-
-			}
-			/*no NULL value should be set */
-			isnull = 0;
-		    }
-		}
-	    }
-
-	    /*Set the NULL values */
-	    if (isnull == 1) {
-		if (typeIntern == FCELL_TYPE)
-		    Rast_set_null_value(&fcell[x], 1, FCELL_TYPE);
-		else if (typeIntern == DCELL_TYPE)
-		    Rast_set_null_value(&dcell[x], 1, DCELL_TYPE);
-	    }
-	}
-
-	/*Write the data to the output map */
-	if (typeIntern == FCELL_TYPE)
-	    Rast_put_f_row(outfd, fcell);
-
-	if (typeIntern == DCELL_TYPE)
-	    Rast_put_d_row(outfd, dcell);
+    for (row = 0; row < rows; row++) {
+        G_percent(row, rows - 1, 10);
+
+        /*Read the input map row */
+        Rast_get_row(elevfd, elevrast, row, globalElevMapType);
+
+        for (col = 0, ptr = elevrast; col < cols; col++, ptr =
+            G_incr_void_ptr(ptr, Rast_cell_size(globalElevMapType))) {
+
+            if (Rast_is_null_value(ptr, globalElevMapType)) {
+                if (typeIntern == FCELL_TYPE)
+                    Rast_set_null_value(&fcell[col], 1, FCELL_TYPE);
+                else if (typeIntern == DCELL_TYPE)
+                    Rast_set_null_value(&dcell[col], 1, DCELL_TYPE);
+                continue;
+            }
+
+            /*Read the elevation value */
+            if (globalElevMapType == CELL_TYPE) {
+                intvalue = *(CELL *) ptr;
+                elevation = intvalue;
+            } else if (globalElevMapType == FCELL_TYPE) {
+                fvalue = *(FCELL *) ptr;
+                elevation = fvalue;
+            } else if (globalElevMapType == DCELL_TYPE) {
+                dvalue = *(DCELL *) ptr;
+                elevation = dvalue;
+            }
+
+            /* Compute the coordinates */
+            north = Rast_row_to_northing((double)row + 0.5, &window);
+            east = Rast_col_to_easting((double)col + 0.5, &window);
+
+            /* Get the voxel value */
+            if (typeIntern == FCELL_TYPE)
+                G3d_getRegionValue(map, north, east, elevation, &fcell[col], FCELL_TYPE);
+
+            if (typeIntern == DCELL_TYPE)
+                G3d_getRegionValue(map, north, east, elevation, &dcell[col], DCELL_TYPE);
+        }
+
+        /*Write the data to the output map */
+        if (typeIntern == FCELL_TYPE)
+            Rast_put_f_row(outfd, fcell);
+
+        if (typeIntern == DCELL_TYPE)
+            Rast_put_d_row(outfd, dcell);
     }
     G_debug(3, "\nDone\n");
 
     /*Free the mem */
     if (elevrast)
-	G_free(elevrast);
+        G_free(elevrast);
     if (dcell)
-	G_free(dcell);
+        G_free(dcell);
     if (fcell)
-	G_free(fcell);
+        G_free(fcell);
 }
 
 
 /* ************************************************************************* */
 /* Main function, open the G3D map and create the cross section map ******** */
+
 /* ************************************************************************* */
 int main(int argc, char *argv[])
 {
     G3D_Region region;
     struct Cell_head window2d;
     struct GModule *module;
-    void *map = NULL;		/*The 3D Rastermap */
+    void *map = NULL; /*The 3D Rastermap */
     int changemask = 0;
-    int elevfd = -1, outfd = -1;	/*file descriptors */
+    int elevfd = -1, outfd = -1; /*file descriptors */
     int output_type, cols, rows;
 
     /* Initialize GRASS */
@@ -264,20 +229,20 @@ int main(int argc, char *argv[])
     G_add_keyword(_("raster"));
     G_add_keyword(_("voxel"));
     module->description =
-	_("Creates cross section 2D raster map from 3d raster map based on 2D elevation map");
+        _("Creates cross section 2D raster map from 3d raster map based on 2D elevation map");
 
     /* Get parameters from user */
     set_params();
 
     /* Have GRASS get inputs */
     if (G_parser(argc, argv))
-	exit(EXIT_FAILURE);
+        exit(EXIT_FAILURE);
 
     G_debug(3, "Open 3D raster map %s", param.input->answer);
 
     if (NULL == G_find_grid3(param.input->answer, ""))
-	G3d_fatalError(_("3d raster map <%s> not found"),
-		       param.input->answer);
+        G3d_fatalError(_("3d raster map <%s> not found"),
+                       param.input->answer);
 
     /* Figure out the region from the map */
     G3d_initDefaults();
@@ -289,14 +254,14 @@ int main(int argc, char *argv[])
 
     /*If not equal, set the 2D windows correct */
     if (rows != region.rows || cols != region.cols) {
-	G_message
-	    (_("The 2d and 3d region settings are different. I will use the g3d settings to adjust the 2d region."));
-	G_get_set_window(&window2d);
-	window2d.ns_res = region.ns_res;
-	window2d.ew_res = region.ew_res;
-	window2d.rows = region.rows;
-	window2d.cols = region.cols;
-	Rast_set_window(&window2d);
+        G_message
+            (_("The 2d and 3d region settings are different. I will use the g3d settings to adjust the 2d region."));
+        G_get_set_window(&window2d);
+        window2d.ns_res = region.ns_res;
+        window2d.ew_res = region.ew_res;
+        window2d.rows = region.rows;
+        window2d.cols = region.cols;
+        Rast_set_window(&window2d);
     }
 
 
@@ -305,78 +270,77 @@ int main(int argc, char *argv[])
 
     /*******************/
     map = G3d_openCellOld(param.input->answer,
-			  G_find_grid3(param.input->answer, ""),
-			  &region, G3D_TILE_SAME_AS_FILE,
-			  G3D_USE_CACHE_DEFAULT);
+                          G_find_grid3(param.input->answer, ""),
+                          &region, G3D_TILE_SAME_AS_FILE,
+                          G3D_USE_CACHE_DEFAULT);
 
     if (map == NULL)
-	G3d_fatalError(_("Error opening 3d raster map <%s>"),
-		       param.input->answer);
+        G3d_fatalError(_("Error opening 3d raster map <%s>"),
+                       param.input->answer);
 
     /*Get the output type */
     output_type = G3d_fileTypeMap(map);
 
     if (output_type == FCELL_TYPE || output_type == DCELL_TYPE) {
 
-	/********************************/
-	/*Open the elevation raster map */
+        /********************************/
+        /*Open the elevation raster map */
 
-	/********************************/
+        /********************************/
 
-	elevfd = Rast_open_old(param.elevation->answer, "");
+        elevfd = Rast_open_old(param.elevation->answer, "");
 
-	globalElevMapType = Rast_get_map_type(elevfd);
+        globalElevMapType = Rast_get_map_type(elevfd);
 
-	/**********************/
-	/*Open the Outputmap */
+        /**********************/
+        /*Open the Outputmap */
 
-	/**********************/
+        /**********************/
 
-	if (G_find_raster2(param.output->answer, ""))
-	    G_message(_("Output map already exists. Will be overwritten!"));
+        if (G_find_raster2(param.output->answer, ""))
+            G_message(_("Output map already exists. Will be overwritten!"));
 
-	if (output_type == FCELL_TYPE)
-	    outfd = Rast_open_new(param.output->answer, FCELL_TYPE);
-	else if (output_type == DCELL_TYPE)
-	    outfd = Rast_open_new(param.output->answer, DCELL_TYPE);
+        if (output_type == FCELL_TYPE)
+            outfd = Rast_open_new(param.output->answer, FCELL_TYPE);
+        else if (output_type == DCELL_TYPE)
+            outfd = Rast_open_new(param.output->answer, DCELL_TYPE);
 
-	/*if requested set the Mask on */
-	if (param.mask->answer) {
-	    if (G3d_maskFileExists()) {
-		changemask = 0;
-		if (G3d_maskIsOff(map)) {
-		    G3d_maskOn(map);
-		    changemask = 1;
-		}
-	    }
-	}
+        /*if requested set the Mask on */
+        if (param.mask->answer) {
+            if (G3d_maskFileExists()) {
+                changemask = 0;
+                if (G3d_maskIsOff(map)) {
+                    G3d_maskOn(map);
+                    changemask = 1;
+                }
+            }
+        }
 
-	/************************/
-	/*Create the Rastermaps */
+        /************************/
+        /*Create the Rastermaps */
 
-	/************************/
-	rast3d_cross_section(map, region, elevfd, outfd);
+        /************************/
+        rast3d_cross_section(map, region, elevfd, outfd);
 
-	/*We set the Mask off, if it was off before */
-	if (param.mask->answer) {
-	    if (G3d_maskFileExists())
-		if (G3d_maskIsOn(map) && changemask)
-		    G3d_maskOff(map);
-	}
+        /*We set the Mask off, if it was off before */
+        if (param.mask->answer) {
+            if (G3d_maskFileExists())
+                if (G3d_maskIsOn(map) && changemask)
+                    G3d_maskOff(map);
+        }
 
-	Rast_close(outfd);
-	Rast_close(elevfd);
+        Rast_close(outfd);
+        Rast_close(elevfd);
 
-    }
-    else {
-	fatal_error(map, -1, -1,
-		    _("Wrong G3D Datatype! Cannot create raster map."));
+    } else {
+        fatal_error(map, -1, -1,
+                    _("Wrong G3D Datatype! Cannot create raster map."));
     }
 
     /* Close files and exit */
     if (!G3d_closeCell(map))
-	G3d_fatalError(_("Could not close G3D map <%s>"),
-		       param.input->answer);
+        G3d_fatalError(_("Could not close G3D map <%s>"),
+                       param.input->answer);
 
     return (EXIT_SUCCESS);
 }

raster3d/r3.cross.rast/test.r3.cross.rast.sh

@@ -0,0 +1,35 @@
+# This script tests the r3.cross.rast to compute
+# cross section raster maps based on a raster3d and elevation map
+
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster and voxel data with r.mapcalc and r3.mapcalc
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=100 b=0 t=50 res=10 res3=10 -p3
+# We create several evlevation maps to create slices of the voxel map
+# We start from bottom and raise to the top
+r.mapcalc --o expr="elev_0 = 0"
+r.mapcalc --o expr="elev_1 = 5"
+r.mapcalc --o expr="elev_2 = 15"
+r.mapcalc --o expr="elev_3 = 25"
+r.mapcalc --o expr="elev_4 = 35"
+r.mapcalc --o expr="elev_5 = 45"
+r.mapcalc --o expr="elev_NAN = 50"
+r.mapcalc --o expr="elev_cross = float(col()* 5)"
+# Now create a voxel map with value = col + row + depth. Beware the 
+# raster3d module count from south to north.
+r3.mapcalc --o expr="volume = col() + row() + depth()"
+# Add null value information
+r3.mapcalc --o expr="volume_null = if(row() == 1 || row() == 5, null(), volume)"
+
+# We @test the creation of slices and a cross section of the voxel map. Reference data
+# for @raster map validation is located in the r3.cross.rast source directory.
+# Slice 0 and 1 should be identical. The last slice should be NAN.
+r3.cross.rast --o input=volume_null elevation=elev_0 output=test_cross_section_slice_0
+r3.cross.rast --o input=volume_null elevation=elev_1 output=test_cross_section_slice_1
+r3.cross.rast --o input=volume_null elevation=elev_2 output=test_cross_section_slice_2
+r3.cross.rast --o input=volume_null elevation=elev_3 output=test_cross_section_slice_3
+r3.cross.rast --o input=volume_null elevation=elev_4 output=test_cross_section_slice_4
+r3.cross.rast --o input=volume_null elevation=elev_5 output=test_cross_section_slice_5
+r3.cross.rast --o input=volume_null elevation=elev_NAN output=test_cross_section_slice_NAN
+r3.cross.rast --o input=volume_null elevation=elev_cross output=test_cross_section_result

raster3d/r3.cross.rast/test_cross_section_result.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+10 12 13 15 16 18 19 21 22 * 
+9 11 12 14 15 17 18 20 21 * 
+8 10 11 13 14 16 17 19 20 * 
+* * * * * * * * * * 
+6 8 9 11 12 14 15 17 18 * 
+5 7 8 10 11 13 14 16 17 * 
+4 6 7 9 10 12 13 15 16 * 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_0.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+10 11 12 13 14 15 16 17 18 19 
+9 10 11 12 13 14 15 16 17 18 
+8 9 10 11 12 13 14 15 16 17 
+* * * * * * * * * * 
+6 7 8 9 10 11 12 13 14 15 
+5 6 7 8 9 10 11 12 13 14 
+4 5 6 7 8 9 10 11 12 13 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_1.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+10 11 12 13 14 15 16 17 18 19 
+9 10 11 12 13 14 15 16 17 18 
+8 9 10 11 12 13 14 15 16 17 
+* * * * * * * * * * 
+6 7 8 9 10 11 12 13 14 15 
+5 6 7 8 9 10 11 12 13 14 
+4 5 6 7 8 9 10 11 12 13 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_2.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+11 12 13 14 15 16 17 18 19 20 
+10 11 12 13 14 15 16 17 18 19 
+9 10 11 12 13 14 15 16 17 18 
+* * * * * * * * * * 
+7 8 9 10 11 12 13 14 15 16 
+6 7 8 9 10 11 12 13 14 15 
+5 6 7 8 9 10 11 12 13 14 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_3.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+12 13 14 15 16 17 18 19 20 21 
+11 12 13 14 15 16 17 18 19 20 
+10 11 12 13 14 15 16 17 18 19 
+* * * * * * * * * * 
+8 9 10 11 12 13 14 15 16 17 
+7 8 9 10 11 12 13 14 15 16 
+6 7 8 9 10 11 12 13 14 15 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_4.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+13 14 15 16 17 18 19 20 21 22 
+12 13 14 15 16 17 18 19 20 21 
+11 12 13 14 15 16 17 18 19 20 
+* * * * * * * * * * 
+9 10 11 12 13 14 15 16 17 18 
+8 9 10 11 12 13 14 15 16 17 
+7 8 9 10 11 12 13 14 15 16 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_5.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+14 15 16 17 18 19 20 21 22 23 
+13 14 15 16 17 18 19 20 21 22 
+12 13 14 15 16 17 18 19 20 21 
+* * * * * * * * * * 
+10 11 12 13 14 15 16 17 18 19 
+9 10 11 12 13 14 15 16 17 18 
+8 9 10 11 12 13 14 15 16 17 
+* * * * * * * * * * 

raster3d/r3.cross.rast/test_cross_section_slice_NAN.ref

@@ -0,0 +1,14 @@
+north: 80
+south: 0
+east: 100
+west: 0
+rows: 8
+cols: 10
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * * 
+* * * * * * * * * *