[GRASS-SVN] r45442 - grass-addons/raster/r.viewshed
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Feb 21 01:42:30 EST 2011
Author: mmetz
Date: 2011-02-20 22:42:30 -0800 (Sun, 20 Feb 2011)
New Revision: 45442
Modified:
grass-addons/raster/r.viewshed/distribute.cc
grass-addons/raster/r.viewshed/eventlist.cc
grass-addons/raster/r.viewshed/grass.cc
grass-addons/raster/r.viewshed/grid.h
grass-addons/raster/r.viewshed/main.cc
grass-addons/raster/r.viewshed/rbbst.cc
grass-addons/raster/r.viewshed/statusstructure.cc
grass-addons/raster/r.viewshed/statusstructure.h
grass-addons/raster/r.viewshed/viewshed.cc
grass-addons/raster/r.viewshed/visibility.h
Log:
fix correction for earth curvature, fix viewing angle, add new option for target offset, show progress in GRASS
Modified: grass-addons/raster/r.viewshed/distribute.cc
===================================================================
--- grass-addons/raster/r.viewshed/distribute.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/distribute.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -738,7 +738,7 @@
viscell.row = sn.row;
viscell.col = sn.col;
- if (max <= sn.gradient) {
+ if (max <= sn.gradient_offset) {
/*the point is visible */
viscell.angle =
get_vertical_angle(*vp, sn, viewOptions.doCurv);
Modified: grass-addons/raster/r.viewshed/eventlist.cc
===================================================================
--- grass-addons/raster/r.viewshed/eventlist.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/eventlist.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -42,6 +42,11 @@
#include <stdio.h>
#include <assert.h>
+extern "C"
+{
+#include <grass/gis.h>
+}
+
#include "eventlist.h"
@@ -145,6 +150,7 @@
calculate_angle(double eventX, double eventY,
double viewpointX, double viewpointY)
{
+ double angle = atan(fabs(eventY - viewpointY) / fabs(eventX - viewpointX));
/*M_PI is defined in math.h to represent 3.14159... */
if (viewpointY == eventY && eventX > viewpointX) {
@@ -152,7 +158,7 @@
}
else if (eventX > viewpointX && eventY < viewpointY) {
/*first quadrant */
- return atan((viewpointY - eventY) / (eventX - viewpointX));
+ return angle;
}
else if (viewpointX == eventX && viewpointY > eventY) {
@@ -162,8 +168,7 @@
}
else if (eventX < viewpointX && eventY < viewpointY) {
/*second quadrant */
- return ((M_PI) / 2 +
- atan((viewpointX - eventX) / (viewpointY - eventY)));
+ return ((M_PI) - angle);
}
else if (viewpointY == eventY && eventX < viewpointX) {
@@ -173,17 +178,16 @@
}
else if (eventY > viewpointY && eventX < viewpointX) {
/*3rd quadrant */
- return (M_PI + atan((eventY - viewpointY) / (viewpointX - eventX)));
+ return (M_PI + angle);
}
else if (viewpointX == eventX && viewpointY < eventY) {
/*between 3rd and 4th quadrant */
- return (M_PI * 3.0 / 2.0);
+ return ((M_PI) + (M_PI) / 2);
}
else if (eventX > viewpointX && eventY > viewpointY) {
/*4th quadrant */
- return (M_PI * 3.0 / 2.0 +
- atan((eventX - viewpointX) / (eventY - viewpointY)));
+ return ((M_PI) * 2.0 - angle);
}
assert(eventX == viewpointX && eventY == viewpointY);
return 0;
@@ -206,7 +210,6 @@
calculate_event_position(AEvent e, dimensionType viewpointRow,
dimensionType viewpointCol, double *y, double *x)
{
-
assert(x && y);
*x = 0;
*y = 0;
@@ -413,7 +416,20 @@
/* it is not too smart to compare floats */
if ((int)maxDist == INFINITY_DISTANCE)
return 0;
+
+#ifdef _GRASS_
+ if (maxDist < G_distance(G_col_to_easting(vp.col + 0.5, &hd.window),
+ G_row_to_northing(vp.row + 0.5, &hd.window),
+ G_col_to_easting(col + 0.5, &hd.window),
+ G_row_to_northing(row + 0.5, &hd.window))) {
+ return 1;
+ }
+ else {
+ return 0;
+ }
+#else
+
double dif_x, dif_y, sqdist;
dif_x = (vp.col - col);
@@ -428,6 +444,7 @@
else {
return 0;
}
+#endif
}
Modified: grass-addons/raster/r.viewshed/grass.cc
===================================================================
--- grass-addons/raster/r.viewshed/grass.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/grass.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -63,22 +63,22 @@
*/
float adjust_for_curvature(Viewpoint vp, dimensionType row,
dimensionType col, float h,
- ViewOptions viewOptions)
+ ViewOptions viewOptions, GridHeader *hd)
{
if (!viewOptions.doCurv)
return h;
assert(viewOptions.ellps_a != 0);
- double dif_x, dif_y, sqdist;
- dif_x = (vp.col - col);
- dif_y = (vp.row - row);
- sqdist =
- (dif_x * dif_x +
- dif_y * dif_y) * viewOptions.cellsize * viewOptions.cellsize;
+ /* distance must be in meters because ellps_a is in meters */
- return h - (sqdist / (2.0 * viewOptions.ellps_a));
+ double dist = G_distance(G_col_to_easting(vp.col + 0.5, &(hd->window)),
+ G_row_to_northing(vp.row + 0.5, &(hd->window)),
+ G_col_to_easting(col + 0.5, &(hd->window)),
+ G_row_to_northing(row + 0.5, &(hd->window)));
+
+ return h - ((dist * dist) / (2.0 * viewOptions.ellps_a));
}
@@ -149,7 +149,7 @@
MemoryVisibilityGrid * visgrid)
{
- G_verbose_message(_("computing events ..."));
+ G_message(_("computing events ..."));
assert(eventList && vp && visgrid);
//GRASS should be defined
@@ -194,9 +194,10 @@
dimensionType i, j, k;
double ax, ay;
AEvent e;
+ int nrows = G_window_rows();
e.angle = -1;
- for (i = 0; i < G_window_rows(); i++) {
+ for (i = 0; i < nrows; i++) {
/*read in the raster row */
int rasterRowResult = G_get_raster_row(infd, inrast, i, data_type);
@@ -204,6 +205,8 @@
G_fatal_error(_("Coord not read from row %d of <%s>"), i,
rastName);
+ G_percent(i, nrows, 2);
+
/*fill event list with events from this row */
for (j = 0; j < G_window_cols(); j++) {
e.row = i;
@@ -226,7 +229,7 @@
}
/* adjust for curvature */
- e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions);
+ e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions, hd);
/*write it into the row of data going through the viewpoint */
if (i == vp->row) {
@@ -236,6 +239,10 @@
/* set the viewpoint, and don't insert it into eventlist */
if (i == vp->row && j == vp->col) {
set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
+ if (viewOptions.tgtElev > 0)
+ vp->target_offset = viewOptions.tgtElev;
+ else
+ vp->target_offset = 0.;
if (isnull) {
/*what to do when viewpoint is NODATA ? */
G_warning(_("Viewpoint is NODATA."));
@@ -285,6 +292,7 @@
}
}
+ G_percent(nrows, nrows, 2);
G_message(_("...done creating event list"));
G_close_cell(infd);
@@ -310,7 +318,6 @@
double **data,
IOVisibilityGrid * visgrid)
{
-
G_message(_("computing events ..."));
assert(rastName && vp && hd && visgrid);
@@ -352,11 +359,14 @@
dimensionType i, j, k;
double ax, ay;
AEvent e;
+ int nrows = G_window_rows();
e.angle = -1;
/*start scanning through the grid */
- for (i = 0; i < G_window_rows(); i++) {
+ for (i = 0; i < nrows; i++) {
+
+ G_percent(i, nrows, 2);
/*read in the raster row */
if (G_get_raster_row(infd, inrast, i, data_type) <= 0)
@@ -385,7 +395,7 @@
}
/* adjust for curvature */
- e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions);
+ e.elev = adjust_for_curvature(*vp, i, j, e.elev, viewOptions, hd);
if (data != NULL) {
@@ -398,6 +408,10 @@
/* set the viewpoint */
if (i == vp->row && j == vp->col) {
set_viewpoint_elev(vp, e.elev + viewOptions.obsElev);
+ if (viewOptions.tgtElev > 0)
+ vp->target_offset = viewOptions.tgtElev;
+ else
+ vp->target_offset = 0.;
/*what to do when viewpoint is NODATA */
if (is_nodata(hd, e.elev)) {
G_warning("Viewpoint is NODATA.");
@@ -443,6 +457,7 @@
} /* for j */
} /* for i */
+ G_percent(nrows, nrows, 2);
G_message(_("...done creating event list\n"));
G_close_cell(infd);
Modified: grass-addons/raster/r.viewshed/grid.h
===================================================================
--- grass-addons/raster/r.viewshed/grid.h 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/grid.h 2011-02-21 06:42:30 UTC (rev 45442)
@@ -47,8 +47,13 @@
#include <stdio.h>
#include <limits.h>
+#ifdef __GRASS__
+extern "C"
+{
+#include <grass/gis.h>
+}
+#endif
-
/* this should accomodate grid sizes up to 2^16-1=65,535
If this is not enough, change type and recompile */
typedef unsigned short int dimensionType;
@@ -63,6 +68,11 @@
float yllcorner; /*yllcorner refers to the southern edge of grid */
float cellsize; /*the resolution of the grid */
float nodata_value; /*the value that represents missing data */
+
+#ifdef __GRASS__
+ /* GRASS */
+ struct Cell_head window;
+#endif
} GridHeader;
Modified: grass-addons/raster/r.viewshed/main.cc
===================================================================
--- grass-addons/raster/r.viewshed/main.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/main.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -198,6 +198,7 @@
#ifdef __GRASS__
hd = read_header_from_GRASS(viewOptions.inputfname, ®ion);
assert(hd);
+ G_get_set_window(&(hd->window));
/* LT: there is no need to exit if viewpoint is outside grid,
the algorithm will work correctly in theory. But this
@@ -239,6 +240,9 @@
G_warning(_("Problems obtaining current ellipsoid parameters, using sphere (6370997.0)"));
viewOptions.ellps_a = 6370997.00;
}
+
+ G_begin_distance_calculations();
+
#else
/* in standalone mode we do not know how to adjust for curvature */
assert(viewOptions.doCurv == 0);
@@ -247,7 +251,6 @@
-
/* ************************************************************ */
/* decide whether the computation of the viewshed will take place
in-memory or in external memory */
@@ -567,6 +570,18 @@
obsElevOpt->answer = "1.75";
obsElevOpt->guisection = _("Input_options");
+ /* target elevation offset */
+ struct Option *tgtElevOpt;
+
+ tgtElevOpt = G_define_option();
+ tgtElevOpt->key = "tgt_elev";
+ tgtElevOpt->type = TYPE_DOUBLE;
+ tgtElevOpt->required = NO;
+ tgtElevOpt->key_desc = "value";
+ tgtElevOpt->description = _("Offset for target elevation above the ground");
+ tgtElevOpt->answer = "0.0";
+ tgtElevOpt->guisection = _("Input_options");
+
/* max distance */
struct Option *maxDistOpt;
@@ -605,7 +620,7 @@
#ifdef __MINGW32__
streamdirOpt->answer = G_convert_dirseps_from_host(G_store(getenv("TEMP")));
#else
- streamdirOpt->answer = G_store("/var/tmp/");
+ streamdirOpt->answer = "/var/tmp/";
#endif
streamdirOpt->description=
_("Directory to hold temporary files (they can be large)");
@@ -621,6 +636,8 @@
strcpy(viewOptions->streamdir,streamdirOpt->answer);
viewOptions->obsElev = atof(obsElevOpt->answer);
+ if(tgtElevOpt->answer)
+ viewOptions->tgtElev = atof(tgtElevOpt->answer);
viewOptions->maxDist = atof(maxDistOpt->answer);
if (viewOptions->maxDist < 0 && viewOptions->maxDist != INFINITY_DISTANCE) {
@@ -646,6 +663,8 @@
*memSizeBytes = (long long)memSizeMB;
*memSizeBytes = (*memSizeBytes) << 20;
+ G_get_set_window(window);
+
/*The algorithm runs with the viewpoint row and col, so depending
on if the row_col flag is present we either need to store the
row and col, or convert the lat-lon coordinates to row and
Modified: grass-addons/raster/r.viewshed/rbbst.cc
===================================================================
--- grass-addons/raster/r.viewshed/rbbst.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/rbbst.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -632,7 +632,7 @@
-/*find max within the max key */
+/* find max within the max key */
double find_max_value_within_key(TreeNode * root, double maxKey)
{
TreeNode *keyNode = search_for_node(root, maxKey);
Modified: grass-addons/raster/r.viewshed/statusstructure.cc
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/statusstructure.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -74,10 +74,11 @@
/*calculate and return the angle in degrees */
assert(fabs(sn.dist2vp) > 0.001);
+
if (diffElev >= 0.0)
- return (atan(sn.dist2vp / diffElev) * (180 / M_PI));
+ return (atan(sqrt(sn.dist2vp) / diffElev) * (180 / M_PI));
else
- return ((atan(fabs(diffElev) / sn.dist2vp) * (180 / M_PI)) + 90);
+ return ((atan(fabs(diffElev) / sqrt(sn.dist2vp)) * (180 / M_PI)) + 90);
}
@@ -105,19 +106,28 @@
/*given a StatusNode, fill in its dist2vp and gradient */
void calculate_dist_n_gradient(StatusNode * sn, Viewpoint * vp)
{
-
assert(sn && vp);
/*sqrt is expensive
//sn->dist2vp = sqrt((float) ( pow(sn->row - vp->row,2.0) +
// pow(sn->col - vp->col,2.0)));
//sn->gradient = (sn->elev - vp->elev)/(sn->dist2vp); */
+
+ double elevDiff = (double)sn->elev - vp->elev;
+
sn->dist2vp = (sn->row - vp->row) * (sn->row - vp->row) +
(sn->col - vp->col) * (sn->col - vp->col);
- sn->gradient =
- (sn->elev - vp->elev) * (sn->elev - vp->elev) / (sn->dist2vp);
+ sn->gradient = elevDiff * elevDiff / (sn->dist2vp);
+ sn->gradient = atan(sqrt(sn->gradient));
/*maintain sign */
if (sn->elev < vp->elev)
sn->gradient = -sn->gradient;
+
+ elevDiff += vp->target_offset;
+ sn->gradient_offset = elevDiff * elevDiff / (sn->dist2vp);
+ sn->gradient_offset = atan(sqrt(sn->gradient_offset));
+ /*maintain sign */
+ if (sn->elev + vp->target_offset < vp->elev)
+ sn->gradient_offset = -sn->gradient_offset;
return;
}
Modified: grass-addons/raster/r.viewshed/statusstructure.h
===================================================================
--- grass-addons/raster/r.viewshed/statusstructure.h 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/statusstructure.h 2011-02-21 06:42:30 UTC (rev 45442)
@@ -58,6 +58,7 @@
float elev; /*elevation of cell */
double dist2vp; /*distance to the viewpoint */
double gradient; /*gradient of the Line of Sight */
+ double gradient_offset; /*gradient of the Line of Sight with local elevation offset */
} StatusNode;
Modified: grass-addons/raster/r.viewshed/viewshed.cc
===================================================================
--- grass-addons/raster/r.viewshed/viewshed.cc 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/viewshed.cc 2011-02-21 06:42:30 UTC (rev 45442)
@@ -315,8 +315,17 @@
long nvis = 0; /*number of visible cells */
AEvent *e;
+#ifdef __GRASS__
+ G_message("Determine visibility...");
+ G_percent(0, 100, 2);
+#endif
for (size_t i = 0; i < nevents; i++) {
+#ifdef __GRASS__
+ int perc = (int)((double)i / nevents * 1000000000.);
+ if (perc > 0)
+ G_percent(perc, 1000000000, 2);
+#endif
/*get out one event at a time and process it according to its type */
e = &(eventList[i]);
@@ -360,7 +369,7 @@
find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
/*the point is visible: store its vertical angle */
- if (max <= sn.gradient) {
+ if (max <= sn.gradient_offset) {
add_result_to_inmem_visibilitygrid(visgrid, sn.row, sn.col,
get_vertical_angle(*vp, sn,
viewOptions.
@@ -457,6 +466,9 @@
/* ------------------------------ */
/*sort the events radially by angle */
+#ifdef __GRASS__
+ G_message("Sort the events...");
+#endif
rt_start(sortEventTime);
sort_event_list(&eventList);
eventList->seek(0); /*this does not seem to be ensured by sort?? */
@@ -471,8 +483,15 @@
structure */
StatusNode sn;
+#ifdef __GRASS__
+ G_message("Calculate distances...");
+#endif
rt_start(sweepTime);
for (dimensionType i = vp->col + 1; i < hd->ncols; i++) {
+
+#ifdef __GRASS__
+ G_percent(i, hd->ncols, 2);
+#endif
sn.col = i;
sn.row = vp->row;
sn.elev = data[i];
@@ -491,6 +510,7 @@
}
}
#ifdef __GRASS__
+ G_percent(hd->ncols, hd->ncols, 2);
G_free(data);
#else
free(data);
@@ -507,8 +527,17 @@
/*printf("nbEvents = %ld\n", (long) nbEvents); */
+#ifdef __GRASS__
+ G_message("Determine visibility...");
+ G_percent(0, 100, 2);
+#endif
for (off_t i = 0; i < nbEvents; i++) {
-
+
+#ifdef __GRASS__
+ int perc = (int)((double)i / nbEvents * 1000000000.);
+ if (perc > 0)
+ G_percent(perc, 1000000000, 2);
+#endif
/*get out one event at a time and process it according to its type */
ae = eventList->read_item(&e);
assert(ae == AMI_ERROR_NO_ERROR);
@@ -553,7 +582,7 @@
find_max_gradient_in_status_struct(status_struct, sn.dist2vp);
/*the point is visible */
- if (max <= sn.gradient) {
+ if (max <= sn.gradient_offset) {
viscell.angle =
get_vertical_angle(*vp, sn, viewOptions.doCurv);
assert(viscell.angle >= 0);
Modified: grass-addons/raster/r.viewshed/visibility.h
===================================================================
--- grass-addons/raster/r.viewshed/visibility.h 2011-02-20 19:33:01 UTC (rev 45441)
+++ grass-addons/raster/r.viewshed/visibility.h 2011-02-21 06:42:30 UTC (rev 45442)
@@ -63,6 +63,7 @@
{
dimensionType row, col;
float elev;
+ float target_offset;
} Viewpoint;
@@ -107,6 +108,9 @@
float obsElev;
/* observer elevation above the terrain */
+ float tgtElev;
+ /* target elevation offset above the terrain */
+
float maxDist;
/* points that are farther than this distance from the viewpoint are
not visible */
More information about the grass-commit
mailing list