r.viewshed: add latlon support

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

raster/r.viewshed/eventlist.cpp

@@ -522,13 +522,28 @@ double
 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,12 +554,27 @@ get_square_distance_from_viewpoint_with_print(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);
+    }
 
     print_event(a, 2);
     G_debug(2, " pos= (%.3f. %.3f) sqdist=%.3f", eventx, eventy, dist);
@@ -578,7 +608,7 @@ int is_point_outside_max_dist(Viewpoint vp, GridHeader hd,
 
 
 /* ------------------------------------------------------------ 
-   //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,11 +624,43 @@ int DistanceCompare::compare(const AEvent & a, const AEvent & b)
     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;

raster/r.viewshed/statusstructure.cpp

@@ -116,10 +116,21 @@ void calculate_dist_n_gradient(StatusNode * sn, double elev,
        //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;
@@ -162,10 +173,22 @@ void calculate_event_gradient(StatusNode * sn, int e_idx,
        //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));