# [GRASS-SVN] r48510 - grass/trunk/raster/r.surf.area

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Sep 27 18:12:47 EDT 2011

```Author: martinl
Date: 2011-09-27 15:12:46 -0700 (Tue, 27 Sep 2011)
New Revision: 48510

grass/trunk/raster/r.surf.area/area.c
grass/trunk/raster/r.surf.area/local_proto.h
Modified:
grass/trunk/raster/r.surf.area/main.c
grass/trunk/raster/r.surf.area/r.surf.area.html
Log:
r.surf.area: major clean up for G7

===================================================================
--- grass/trunk/raster/r.surf.area/area.c	                        (rev 0)
+++ grass/trunk/raster/r.surf.area/area.c	2011-09-27 22:12:46 UTC (rev 48510)
@@ -0,0 +1,124 @@
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+
+#include "local_proto.h"
+
+void add_row_area(DCELL * top, DCELL * bottom, double sz, struct Cell_head *w,
+		  double *low, double *high)
+{
+    double guess1, guess2, mag, tedge1[3], tedge2[3], crossp[3];
+    int col;
+
+    for (col = 0; col < w->cols - 1; col++) {
+
+	/*
+	   For each cell**, we triangulate the four corners in
+	   two different ways, 1) UppperLeft to LowerRight diagonal
+	   and 2) LowerLeft to UpperRight diagonal.  Then we add the
+	   smaller of the two areas to "low" and the greater of
+	   the two areas to "high".
+
+	   ** here, the "cell" is actually the quadrangle formed by
+	   the center point of four cells, since these are the
+	   known elevation points.
+	 */
+
+	/* If NAN go to next or we get NAN for everything */
+	if (Rast_is_d_null_value(&(bottom[col + 1])) ||
+	    Rast_is_d_null_value(&(top[col])) ||
+	    Rast_is_d_null_value(&(top[col + 1])) ||
+	    Rast_is_d_null_value(&(bottom[col]))
+	    )
+	    continue;
+
+	/* guess1 --- ul to lr diag */
+	{
+	    tedge1[X] = w->ew_res;
+	    tedge1[Y] = -w->ns_res;
+	    tedge1[Z] = sz * (bottom[col + 1] - top[col]);
+
+	    /* upper */
+	    tedge2[X] = 0.0;
+	    tedge2[Y] = w->ns_res;
+	    tedge2[Z] = sz * (top[col + 1] - bottom[col + 1]);
+
+	    v3cross(tedge1, tedge2, crossp);
+	    v3mag(crossp, &mag);
+	    guess1 = .5 * mag;
+
+	    /* lower */
+	    tedge2[X] = -w->ew_res;
+	    tedge2[Y] = 0.0;
+	    tedge2[Z] = sz * (bottom[col] - bottom[col + 1]);
+
+	    v3cross(tedge1, tedge2, crossp);
+	    v3mag(crossp, &mag);
+	    guess1 += .5 * mag;
+	}
+
+	/* guess2 --- ll to ur diag */
+	{
+	    tedge1[X] = w->ew_res;
+	    tedge1[Y] = w->ns_res;
+	    tedge1[Z] = sz * (top[col + 1] - bottom[col]);
+
+	    /* upper */
+	    tedge2[X] = -w->ew_res;
+	    tedge2[Y] = 0.0;
+	    tedge2[Z] = sz * (top[col + 1] - top[col + 1]);
+
+	    v3cross(tedge1, tedge2, crossp);
+	    v3mag(crossp, &mag);
+	    guess2 = .5 * mag;
+
+	    /* lower */
+	    tedge2[X] = 0.0;
+	    tedge2[Y] = -w->ns_res;
+	    tedge2[Z] = sz * (bottom[col + 1] - top[col + 1]);
+
+	    v3cross(tedge1, tedge2, crossp);
+	    v3mag(crossp, &mag);
+	    guess2 += .5 * mag;
+	}
+	*low += (guess1 < guess2) ? guess1 : guess2;
+	*high += (guess1 < guess2) ? guess2 : guess1;
+
+    }				/* ea col */
+
+}
+
+/* calculate the running area of null data cells */
+{
+    int col;
+
+    for (col = 0; col < region->cols; col++) {
+	if (Rast_is_d_null_value(&(rast[col]))) {
+	    *area += region->ew_res * region->ns_res;
+	}
+    }
+}
+
+/* return the cross product v3 = v1 cross v2 */
+void v3cross(double v1[3], double v2[3], double v3[3])
+{
+    v3[X] = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
+    v3[Y] = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
+    v3[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
+}
+
+/* magnitude of vector */
+void v3mag(double v1[3], double *mag)
+{
+    *mag = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
+}
+
+double conv_value(double value, int units)
+{
+    if (units == U_UNDEFINED)
+	return value;
+
+    return value * G_units_to_meters_factor_sq(units);
+}

Property changes on: grass/trunk/raster/r.surf.area/area.c
___________________________________________________________________
+ text/x-csrc
+ native

===================================================================
--- grass/trunk/raster/r.surf.area/local_proto.h	                        (rev 0)
+++ grass/trunk/raster/r.surf.area/local_proto.h	2011-09-27 22:12:46 UTC (rev 48510)
@@ -0,0 +1,10 @@
+#define X 0
+#define Y 1
+#define Z 2
+
+		  double *, double *);
+void v3cross(double[], double[], double[]);
+void v3mag(double[], double *);
+double conv_value(double, int);

Property changes on: grass/trunk/raster/r.surf.area/local_proto.h
___________________________________________________________________
+ text/x-chdr
+ native

Modified: grass/trunk/raster/r.surf.area/main.c
===================================================================
--- grass/trunk/raster/r.surf.area/main.c	2011-09-27 21:30:04 UTC (rev 48509)
+++ grass/trunk/raster/r.surf.area/main.c	2011-09-27 22:12:46 UTC (rev 48510)
@@ -52,32 +52,22 @@

#include <stdlib.h>
#include <stdio.h>
-#include <math.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>

+#include "local_proto.h"

-#define X 0
-#define Y 1
-#define Z 2
-
-
-		  double *, double *);
-void v3cross(double v1[3], double v2[3], double v3[3]);
-void v3mag(double v1[3], double *mag);
-
-
int main(int argc, char *argv[])
{
struct GModule *module;
-    struct Option *surf, *vscale;
+    struct {
+	struct Option *surf, *vscale, *units;
+    } opt;
DCELL *cell_buf[2];
int row;
-    int cellfile = -1;
+    int cellfile, units;
double minarea, maxarea, sz, nullarea;

G_gisinit(argv[0]);
@@ -86,40 +76,44 @@
-    module->description = _("Surface area estimation for rasters.");
+    module->description = _("Prints estimation of surface area for raster map.");

-    surf = G_define_option();
-    surf->key = "input";
-    surf->type = TYPE_STRING;
-    surf->required = YES;
-    surf->multiple = NO;
-    surf->gisprompt = "old,cell,Raster";
-    surf->description = _("Raster file for surface");
-
-    vscale = G_define_option();
-    vscale->key = "vscale";
-    vscale->type = TYPE_DOUBLE;
-    vscale->required = NO;
-    vscale->multiple = NO;
-    vscale->description = _("Vertical scale");
-
+    opt.surf = G_define_standard_option(G_OPT_R_MAP);
+
+    opt.vscale = G_define_option();
+    opt.vscale->key = "vscale";
+    opt.vscale->type = TYPE_DOUBLE;
+    opt.vscale->required = NO;
+    opt.vscale->multiple = NO;
+    opt.vscale->description = _("Vertical scale");
+
+    opt.units = G_define_standard_option(G_OPT_M_UNITS);
+    opt.units->label = _("Output units");
+    opt.units->description = _("Default: square map units");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);

-    else
-	sz = 1.0;
-
+	G_verbose_message(_("Output in '%s'"), G_get_units_name(units, TRUE, TRUE));
+    }
+    else {
+	units = U_UNDEFINED;
+	G_verbose_message(_("Output in 'square map units'"));
+    }
+
G_get_set_window(&w);

/* open raster map for reading */

cell_buf[0] = (DCELL *) G_malloc(w.cols * Rast_cell_size(DCELL_TYPE));
cell_buf[1] = (DCELL *) G_malloc(w.cols * Rast_cell_size(DCELL_TYPE));

-    fprintf(stdout, "\n");
{
DCELL *top, *bottom;

@@ -146,144 +140,29 @@
G_free(cell_buf[1]);
Rast_close(cellfile);

-    {				/* report */
+    {	/* report */
double reg_area, flat_area, estavg;

flat_area = (w.cols - 1) * (w.rows - 1) * w.ns_res * w.ew_res;
reg_area = w.cols * w.rows * w.ns_res * w.ew_res;
estavg = (minarea + maxarea) / 2.0;

-	fprintf(stdout, "Null value area ignored in calculation %e\n",
-		nullarea);
-	fprintf(stdout, "Plan area used in calculation: %e\n", flat_area);
-	fprintf(stdout,
-		"Surface Area Calculation(low, high, avg):\n\t%e %e %e\n",
-		minarea, maxarea, estavg);
+	fprintf(stdout, "%s %f\n",
+		_("Null value area ignored in calculation:"),
+		conv_value(nullarea, units));
+	fprintf(stdout, "%s %f\n", _("Plan area used in calculation:"),
+		conv_value(flat_area, units));
+	fprintf(stdout, "%s\n\t%f %f %f\n",
+		_("Surface area calculation(low, high, avg):"),
+		conv_value(minarea, units),
+		conv_value(maxarea, units),
+		conv_value(estavg, units));

-	fprintf(stdout, "Current Region plan area: %e\n", reg_area);
-	fprintf(stdout, "Estimated Region Surface Area: %e\n",
-		reg_area * estavg / flat_area);
-	fprintf(stdout, "\nDone.\n");
+	fprintf(stdout, "%s %f\n", _("Current region plan area:"),
+		conv_value(reg_area, units));
+	fprintf(stdout, "%s %f\n", _("Estimated region Surface Area:"),
+		conv_value(reg_area * estavg / flat_area, units));
}

-    return (0);			/* Zero means success */
+    exit(EXIT_SUCCESS);
}
-
-/************************************************************************/
-
-void
-	     double *low, double *high)
-{
-    double guess1, guess2, mag, tedge1[3], tedge2[3], crossp[3];
-    int col;
-
-    for (col = 0; col < w->cols - 1; col++) {
-
-	/*
-	   For each cell**, we triangulate the four corners in
-	   two different ways, 1) UppperLeft to LowerRight diagonal
-	   and 2) LowerLeft to UpperRight diagonal.  Then we add the
-	   smaller of the two areas to "low" and the greater of
-	   the two areas to "high".
-
-	   ** here, the "cell" is actually the quadrangle formed by
-	   the center point of four cells, since these are the
-	   known elevation points.
-	 */
-
-	/* If NAN go to next or we get NAN for everything */
-	if (Rast_is_d_null_value(&(bottom[col + 1])) ||
-	    Rast_is_d_null_value(&(top[col])) ||
-	    Rast_is_d_null_value(&(top[col + 1])) ||
-	    Rast_is_d_null_value(&(bottom[col]))
-	    )
-	    continue;
-
-	/* guess1 --- ul to lr diag */
-	{
-	    tedge1[X] = w->ew_res;
-	    tedge1[Y] = -w->ns_res;
-	    tedge1[Z] = sz * (bottom[col + 1] - top[col]);
-
-	    /* upper */
-	    tedge2[X] = 0.0;
-	    tedge2[Y] = w->ns_res;
-	    tedge2[Z] = sz * (top[col + 1] - bottom[col + 1]);
-
-	    v3cross(tedge1, tedge2, crossp);
-	    v3mag(crossp, &mag);
-	    guess1 = .5 * mag;
-
-	    /* lower */
-	    tedge2[X] = -w->ew_res;
-	    tedge2[Y] = 0.0;
-	    tedge2[Z] = sz * (bottom[col] - bottom[col + 1]);
-
-	    v3cross(tedge1, tedge2, crossp);
-	    v3mag(crossp, &mag);
-	    guess1 += .5 * mag;
-	}
-
-	/* guess2 --- ll to ur diag */
-	{
-	    tedge1[X] = w->ew_res;
-	    tedge1[Y] = w->ns_res;
-	    tedge1[Z] = sz * (top[col + 1] - bottom[col]);
-
-	    /* upper */
-	    tedge2[X] = -w->ew_res;
-	    tedge2[Y] = 0.0;
-	    tedge2[Z] = sz * (top[col + 1] - top[col + 1]);
-
-	    v3cross(tedge1, tedge2, crossp);
-	    v3mag(crossp, &mag);
-	    guess2 = .5 * mag;
-
-	    /* lower */
-	    tedge2[X] = 0.0;
-	    tedge2[Y] = -w->ns_res;
-	    tedge2[Z] = sz * (bottom[col + 1] - top[col + 1]);
-
-	    v3cross(tedge1, tedge2, crossp);
-	    v3mag(crossp, &mag);
-	    guess2 += .5 * mag;
-	}
-	*low += (guess1 < guess2) ? guess1 : guess2;
-	*high += (guess1 < guess2) ? guess2 : guess1;
-
-    }				/* ea col */
-
-}
-
-/************************************************************************/
-/* calculate the running area of null data cells */
-{
-    int col;
-
-    for (col = 0; col < region->cols; col++) {
-	if (Rast_is_d_null_value(&(rast[col]))) {
-	    *area += region->ew_res * region->ns_res;
-	}
-    }
-}
-
-/************************************************************************/
-/* return the cross product v3 = v1 cross v2 */
-void v3cross(double v1[3], double v2[3], double v3[3])
-{
-    v3[X] = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
-    v3[Y] = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
-    v3[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
-}
-
-/************************************************************************/
-/* magnitude of vector */
-void v3mag(double v1[3], double *mag)
-{
-    *mag = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
-}
-
-/************************************************************************/
-/* vim: set softtabstop=4 shiftwidth=4 expandtab: */

Modified: grass/trunk/raster/r.surf.area/r.surf.area.html
===================================================================
--- grass/trunk/raster/r.surf.area/r.surf.area.html	2011-09-27 21:30:04 UTC (rev 48509)
+++ grass/trunk/raster/r.surf.area/r.surf.area.html	2011-09-27 22:12:46 UTC (rev 48510)
@@ -1,66 +1,83 @@
<h2>DESCRIPTION</h2>

-<em>r.surf.area</em> Calculates area of regular 3D triangulated points
-(centers of cells) in current region by adding areas of triangles.
-Therefore, area of a flat surface will be reported as
-(rows + cols -1)*(area of cell) less than area of
-flat region due to a half row and half column missing around the
-perimeter.
+<em>r.surf.area</em> calculates area of regular 3D triangulated points
+(centers of cells) in current region by adding areas of triangles.
+Therefore, area of a flat surface will be reported as (<tt>rows + cols
+-1) * (area of cell)</tt> less than area of flat region due to a half
+row and half column missing around the perimeter.

-<h2>NOTE</h2>
+<h2>NOTES</h2>

-This calculation is heavily dependent on
-data resolution (think of it as a fractal shoreline problem, the more
-resolution the more detail, the more area, etc).  This program uses the
-<em>CURRENT GRASS REGION</em>, not the resolution of the map.  This is
-especially important for surfaces with <tt>NULL</tt> values and highly
-irregular edges.  The program does not [currently] attempt to correct for
-the error introduced by this <em>edge effect</em>.
+This calculation is heavily dependent on data resolution (think of it
+as a fractal shoreline problem, the more resolution the more detail,
+the more area, etc). This module uses the <b>current region
+settings</b>, not the resolution of the raster map. This is especially
+important for surfaces with <tt>NULL</tt> values and highly irregular
+edges. The module does not [currently] attempt to correct for the
+error introduced by this <em>edge effect</em>.

<p>
-
This version actually calculates area twice for each triangle pair,
-keeping a running minimum and maximum area depending on the
-direction of the diagonal used.
+keeping a running minimum and maximum area depending on the direction
+of the diagonal used.

<p>
-
Reported totals are:
<ol>
-<li>"Plan" area of <tt>NULL</tt> values within the current GRASS region
-<li>"Plan" area within calculation region (rows-1 * cols-1 * cellarea)
-<li>Average of the minimum and maximum calculated 3d triangle area
-within this region
-<li>"Plan" area within current GRASS region (rows * cols * cellarea)
-<li>Scaling of calculated area to current GRASS region (see
-<a href="#note">NOTE</a>)
+<li>"Plan" area of <tt>NULL</tt> values within the current GRASS
+region</li>
+<li>"Plan" area within calculation region (<tt>rows-1 * cols-1 *
+cellarea</tt>)</li>
+<li>Average of the minimum and maximum calculated 3d triangle area
+within this region</li>
+<li>"Plan" area within current computational region (<tt>rows * cols *
+cellarea</tt>)</li>
+<li>Scaling of calculated area to current region</li>
</ol>

-<h2>NOTES</h2>
+<p>
+<em>r.surf.area</em> works best when the surface being evaluated
+extends to the edges of the current region and the cell resolution is
+small. Surfaces which are especially long and thin and have highly
+irregular boudaries will tend to have underestimated surface areas.
+Setting a high cell resolution (small area) will greatly reduce this
+impact, but will cause longer processing times.

-<em>r.surf.area</em> works best when the surface being evaluated extends
-to the edges of the current GRASS region and the cell resolution is small.
-Surfaces which are especially long and thin and have highly irregular
-boudaries will tend to have underestimated surface areas.  Setting a
-high cell resolution (small area) will greatly reduce this impact, but will
-cause longer processing times.
+<h2>EXAMPLES</h2>

+<div class="code"><pre>
+g.region rast=elevation
+
+r.surf.area map=elevation units=hectares
+Null value area ignored in calculation: 0.000000
+Plan area used in calculation: 20221.510000
+Surface area calculation(low, high, avg):
+        20294.310421 20320.936368 20307.623395
+Current region plan area: 20250.000000
+Estimated region Surface Area: 20336.234719
+</pre></div>
+

-<em><a href="r.sum.html">r.sum</a></em>,
-<em><a href="r.surf.idw.html">r.surf.idw</a></em>,
-<em><a href="r.surf.idw2.html">r.surf.idw2</a></em>,
-<em><a href="r.surf.fractal.html">r.surf.fractal</a></em>,
-<em><a href="r.surf.gauss.html">r.surf.gauss</a></em>,
-<em><a href="r.volume.html">r.volume</a></em>,
-<em><a href="v.to.rast.html">v.to.rast</a></em>,
-<em><a href="r.slope.aspect.html">r.slope.aspect</a></em>,
-<em><a href="g.region.html">g.region</a></em>
+<em>
+  <a href="r.surf.idw.html">r.surf.idw</a>,
+  <a href="r.surf.idw2.html">r.surf.idw2</a>,
+  <a href="r.surf.fractal.html">r.surf.fractal</a>,
+  <a href="r.surf.gauss.html">r.surf.gauss</a>,
+  <a href="r.volume.html">r.volume</a>,
+  <a href="v.to.rast.html">v.to.rast</a>,
+  <a href="r.slope.aspect.html">r.slope.aspect</a>,
+  <a href="g.region.html">g.region</a>
+</em>

<h2>AUTHOR</h2>
Bill Brown, USACERL December 21, 1994
+<br>
+Modified for floating point rasters and <tt>NULL</tt> values by Eric
+G. Miller (October 17, 2000)
+<br>
+Updated for GRASS 7, and units option by Martin Landa, Czech Technical
+University in Prague, Czech Republic (October 2011)
+
<p>
-Modified for floating point rasters and <tt>NULL</tt> values by
-Eric G. Miller (October 17, 2000)
-
-<p><i>Last changed: \$Date\$</i>
+<i>Last changed: \$Date\$</i>

```