[GRASS-SVN] r49998 - grass/trunk/raster/r.viewshed
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Dec 31 07:58:14 EST 2011
Author: mmetz
Date: 2011-12-31 04:58:14 -0800 (Sat, 31 Dec 2011)
New Revision: 49998
Modified:
grass/trunk/raster/r.viewshed/eventlist.cpp
grass/trunk/raster/r.viewshed/statusstructure.cpp
Log:
r.viewshed: add latlon support
Modified: grass/trunk/raster/r.viewshed/eventlist.cpp
===================================================================
--- grass/trunk/raster/r.viewshed/eventlist.cpp 2011-12-31 03:54:27 UTC (rev 49997)
+++ grass/trunk/raster/r.viewshed/eventlist.cpp 2011-12-31 12:58:14 UTC (rev 49998)
@@ -522,13 +522,28 @@
get_square_distance_from_viewpoint(const AEvent & a, const Viewpoint & vp)
{
- double eventy, eventx;
+ double eventy, eventx, dist;
calculate_event_position(a, vp.row, vp.col, &eventy, &eventx);
- double dist = (eventx - vp.col) * (eventx - vp.col) +
- (eventy - vp.row) * (eventy - vp.row);
- /*don't take sqrt, it is expensive; suffices for comparison */
+ if (G_projection() == PROJECTION_LL) {
+ struct Cell_head window;
+
+ Rast_get_window(&window);
+
+ dist = G_distance(Rast_col_to_easting(vp.col + 0.5, &window),
+ Rast_row_to_northing(vp.row + 0.5, &window),
+ Rast_col_to_easting(eventx + 0.5, &window),
+ Rast_row_to_northing(eventy + 0.5, &window));
+
+ dist = dist * dist;
+ }
+ else {
+ /*don't take sqrt, it is expensive; suffices for comparison */
+ dist = (eventx - vp.col) * (eventx - vp.col) +
+ (eventy - vp.row) * (eventy - vp.row);
+ }
+
return dist;
}
@@ -539,13 +554,28 @@
const Viewpoint & vp)
{
- double eventy, eventx;
+ double eventy, eventx, dist;
calculate_event_position(a, vp.row, vp.col, &eventy, &eventx);
- double dist = (eventx - vp.col) * (eventx - vp.col) +
- (eventy - vp.row) * (eventy - vp.row);
- /*don't take sqrt, it is expensive; suffices for comparison */
+ if (G_projection() == PROJECTION_LL) {
+ struct Cell_head window;
+
+ Rast_get_window(&window);
+
+ dist = G_distance(Rast_col_to_easting(vp.col + 0.5, &window),
+ Rast_row_to_northing(vp.row + 0.5, &window),
+ Rast_col_to_easting(eventx + 0.5, &window),
+ Rast_row_to_northing(eventy + 0.5, &window));
+
+ dist = dist * dist;
+ }
+ else {
+ /*don't take sqrt, it is expensive; suffices for comparison */
+ dist = (eventx - vp.col) * (eventx - vp.col) +
+ (eventy - vp.row) * (eventy - vp.row);
+ }
+
print_event(a, 2);
G_debug(2, " pos= (%.3f. %.3f) sqdist=%.3f", eventx, eventy, dist);
@@ -578,7 +608,7 @@
/* ------------------------------------------------------------
- //note: this is expensive because distance is not storedin the event
+ //note: this is expensive because distance is not stored in the event
//and must be computed on the fly */
int DistanceCompare::compare(const AEvent & a, const AEvent & b)
{
@@ -594,12 +624,44 @@
double eventy, eventx;
calculate_event_position(a, globalVP.row, globalVP.col, &eventy, &eventx);
- da = (eventx - globalVP.col) * (eventx - globalVP.col) +
- (eventy - globalVP.row) * (eventy - globalVP.row);
+ if (G_projection() == PROJECTION_LL) {
+ struct Cell_head window;
+
+ Rast_get_window(&window);
+
+ da = G_distance(Rast_col_to_easting(globalVP.col + 0.5, &window),
+ Rast_row_to_northing(globalVP.row + 0.5, &window),
+ Rast_col_to_easting(eventx + 0.5, &window),
+ Rast_row_to_northing(eventy + 0.5, &window));
+
+ da = da * da;
+ }
+ else {
+ /*don't take sqrt, it is expensive; suffices for comparison */
+ da = (eventx - globalVP.col) * (eventx - globalVP.col) +
+ (eventy - globalVP.row) * (eventy - globalVP.row);
+ }
+
+
calculate_event_position(b, globalVP.row, globalVP.col, &eventy, &eventx);
- db = (eventx - globalVP.col) * (eventx - globalVP.col) +
- (eventy - globalVP.row) * (eventy - globalVP.row);
+ if (G_projection() == PROJECTION_LL) {
+ struct Cell_head window;
+
+ Rast_get_window(&window);
+
+ db = G_distance(Rast_col_to_easting(globalVP.col + 0.5, &window),
+ Rast_row_to_northing(globalVP.row + 0.5, &window),
+ Rast_col_to_easting(eventx + 0.5, &window),
+ Rast_row_to_northing(eventy + 0.5, &window));
+ db = db * db;
+ }
+ else {
+ /*don't take sqrt, it is expensive; suffices for comparison */
+ db = (eventx - globalVP.col) * (eventx - globalVP.col) +
+ (eventy - globalVP.row) * (eventy - globalVP.row);
+ }
+
if (da > db) {
return 1;
}
Modified: grass/trunk/raster/r.viewshed/statusstructure.cpp
===================================================================
--- grass/trunk/raster/r.viewshed/statusstructure.cpp 2011-12-31 03:54:27 UTC (rev 49997)
+++ grass/trunk/raster/r.viewshed/statusstructure.cpp 2011-12-31 12:58:14 UTC (rev 49998)
@@ -116,11 +116,22 @@
//sn->gradient = (sn->elev - vp->elev)/(sn->dist2vp); */
double diffElev = elev - vp->elev;
- double dx = ((double)sn->col - vp->col) * hd.ew_res;
- double dy = ((double)sn->row - vp->row) * hd.ns_res;
- sn->dist2vp = (dx * dx) + (dy * dy);
+ if (G_projection() == PROJECTION_LL) {
+ double dist = G_distance(Rast_col_to_easting(sn->col + 0.5, &(hd.window)),
+ Rast_row_to_northing(sn->row + 0.5, &(hd.window)),
+ Rast_col_to_easting(vp->col + 0.5, &(hd.window)),
+ Rast_row_to_northing(vp->row + 0.5, &(hd.window)));
+ sn->dist2vp = dist * dist;
+ }
+ else {
+ double dx = ((double)sn->col - vp->col) * hd.ew_res;
+ double dy = ((double)sn->row - vp->row) * hd.ns_res;
+
+ sn->dist2vp = (dx * dx) + (dy * dy);
+ }
+
if (diffElev == 0) {
sn->gradient[1] = 0;
return;
@@ -162,11 +173,23 @@
//sn->gradient = (sn->elev - vp->elev)/(sn->dist2vp); */
double diffElev = elev - vp->elev;
- double dx = (col - vp->col) * hd.ew_res;
- double dy = (row - vp->row) * hd.ns_res;
- double dist2vp = (dx * dx) + (dy * dy);
+ double dist2vp;
+ if (G_projection() == PROJECTION_LL) {
+ double dist = G_distance(Rast_col_to_easting(col + 0.5, &(hd.window)),
+ Rast_row_to_northing(row + 0.5, &(hd.window)),
+ Rast_col_to_easting(vp->col + 0.5, &(hd.window)),
+ Rast_row_to_northing(vp->row + 0.5, &(hd.window)));
+ dist2vp = dist * dist;
+ }
+ else {
+ double dx = (col - vp->col) * hd.ew_res;
+ double dy = (row - vp->row) * hd.ns_res;
+
+ dist2vp = (dx * dx) + (dy * dy);
+ }
+
/* PI / 2 above, - PI / 2 below */
sn->gradient[e_idx] = atan(diffElev / sqrt(dist2vp));
More information about the grass-commit
mailing list