[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
Added:
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
option <units> added
Added: grass/trunk/raster/r.surf.area/area.c
===================================================================
--- 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 */
+void add_null_area(DCELL * rast, struct Cell_head *region, double *area)
+{
+ 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
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/raster/r.surf.area/local_proto.h
===================================================================
--- 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
+
+void add_row_area(DCELL *, DCELL *, double, struct Cell_head *,
+ double *, double *);
+void add_null_area(DCELL *, struct Cell_head *, 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
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ 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
-
-
-void add_row_area(DCELL *, DCELL *, double, struct Cell_head *,
- double *, double *);
-void v3cross(double v1[3], double v2[3], double v3[3]);
-void v3mag(double v1[3], double *mag);
-void add_null_area(DCELL *, struct Cell_head *, double *);
-
-
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;
struct Cell_head w;
- int cellfile = -1;
+ int cellfile, units;
double minarea, maxarea, sz, nullarea;
G_gisinit(argv[0]);
@@ -86,40 +76,44 @@
G_add_keyword(_("raster"));
G_add_keyword(_("surface"));
G_add_keyword(_("statistics"));
- module->description = _("Surface area estimation for rasters.");
+ G_add_keyword(_("area estimation"));
+ 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.vscale->answer = "1.0";
+
+ 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);
- if (vscale->answer)
- sz = atof(vscale->answer);
- else
- sz = 1.0;
-
+ sz = atof(opt.vscale->answer);
+ if (opt.units->answer) {
+ units = G_units(opt.units->answer);
+ 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 */
- cellfile = Rast_open_old(surf->answer, "");
+ cellfile = Rast_open_old(opt.surf->answer, "");
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
-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 */
-void add_null_area(DCELL * rast, struct Cell_head *region, double *area)
-{
- 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>
+
<h2>SEE ALSO</h2>
-<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>
More information about the grass-commit
mailing list