r.viewshed: limit viewshed horizontally by specifying two angles

raster/r.viewshed/eventlist.cpp

@@ -605,7 +605,20 @@ int is_point_outside_max_dist(Viewpoint vp, GridHeader hd,
     return 0;
 }
 
-
+int is_point_inside_angle(Viewpoint vp,
+                         dimensionType row, dimensionType col,
+                         float minAngle, float maxAngle)
+{
+    double ang;
+
+    ang = atan2(vp.row -row, col - vp.col) * 180 / M_PI;
+    if (ang < 0)
+        ang += 360;
+    /* all angles are in range 0-360) */
+    if (minAngle < maxAngle)
+        return minAngle <= ang && ang <= maxAngle;
+    return minAngle <= ang || ang <= maxAngle;
+}
 
 /* ------------------------------------------------------------ 
    //note: this is expensive because distance is not stored in the event

raster/r.viewshed/eventlist.h

@@ -74,7 +74,10 @@ is_point_outside_max_dist(Viewpoint vp, GridHeader hd,
 			  dimensionType row, dimensionType col,
 			  float maxDist);
 
-
+/*determines if the point at row,col is within min and max angle (in 0-360, 0 is east, CCW) */
+int is_point_inside_angle(Viewpoint vp,
+                         dimensionType row, dimensionType col,
+                         float minAngle, float maxAngle);
 
 /*sort the event list by the angle around the viewpoint) */
 void sort_event_list(AMI_STREAM < AEvent > **eventList);

raster/r.viewshed/grass.cpp

@@ -309,6 +309,11 @@ init_event_list_in_memory(AEvent * eventList, char *rastName,
 						   hd->nodata_value);
 		continue;
 	    }
+	    if (viewOptions.doDirection &&
+		    !is_point_inside_angle(*vp, i, j, viewOptions.horizontal_angle_min,
+					   viewOptions.horizontal_angle_max)) {
+		continue;
+	    }
 
 	    /* if point is outside maxDist, do NOT include it as an
 	       event */
@@ -533,7 +538,11 @@ AMI_STREAM < AEvent > *init_event_list(char *rastName, Viewpoint * vp,
 		add_result_to_io_visibilitygrid(visgrid, &visCell);
 		continue;
 	    }
-
+	    if (viewOptions.doDirection &&
+		    !is_point_inside_angle(*vp, i, j, viewOptions.horizontal_angle_min,
+					   viewOptions.horizontal_angle_max)) {
+		continue;
+	    }
 	    /* if point is outside maxDist, do NOT include it as an
 	       event */
 	    if (is_point_outside_max_dist

raster/r.viewshed/main.cpp

@@ -151,6 +151,8 @@ int main(int argc, char *argv[])
     viewOptions.doCurv = FALSE;
     viewOptions.doRefr = FALSE;
     viewOptions.refr_coef = 1.0/7.0;
+    viewOptions.horizontal_angle_min = 0;
+    viewOptions.horizontal_angle_max = 360;
 
     parse_args(argc, argv, &vpRow, &vpCol, &viewOptions, &memSizeBytes,
 	       &region);
@@ -530,6 +532,19 @@ parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
     maxDistOpt->answer = infdist;
     maxDistOpt->guisection = _("Settings");
 
+    /* angle range */
+    struct Option *direction;
+
+    direction = G_define_option();
+    direction->key = "direction_range";
+    direction->type = TYPE_DOUBLE;
+    direction->required = NO;
+    direction->key_desc = "min,max";
+    direction->options = "0-360";
+    direction->description =
+	_("Minimum and maximum horizontal angle limiting viewshed (0 is East, counterclockwise)");
+    direction->guisection = _("Settings");
+
     /* atmospheric refraction coeff. 1/7 for visual, 0.325 for radio waves, ... */
     /* in future we might calculate this based on the physics, for now we
        just fudge by the 1/7th approximation.
@@ -608,6 +623,12 @@ parse_args(int argc, char *argv[], int *vpRow, int *vpCol,
     if (viewOptions->maxDist < 0 && viewOptions->maxDist != INFINITY_DISTANCE) {
 	G_fatal_error(_("A negative max distance value is not allowed"));
     }
+    viewOptions->doDirection = 0;
+    if (direction->answer) {
+        viewOptions->horizontal_angle_min = atof(direction->answers[0]);
+        viewOptions->horizontal_angle_max = atof(direction->answers[1]);
+        viewOptions->doDirection = 1;
+    }
 
     viewOptions->doCurv = curvature->answer;
     viewOptions->doRefr = refractionFlag->answer;

raster/r.viewshed/r.viewshed.html

@@ -84,6 +84,12 @@ visibility using option <b>max_distance</b>.  The value entered is in the
 same units as the cell size of the raster.
 
 <p>
+The user can limit view horizontally by specifying a minimum and maximum directions
+using option <b>direction_range</b>.  The angles are in degrees, CCW, East is 0.
+The angles should be between 0 and 360, e.g. direction_range=0,180 (north view),
+or direction_range=270,90 (east view).
+
+<p>
 Main memory usage: By default <em>r.viewshed</em> assumes it has
 500MB of main memory, and sets up its internal data structures so that
 it does not require more than this amount of RAM.  The user can set

raster/r.viewshed/testsuite/test_r_viewshed.py

@@ -60,6 +60,16 @@ class TestViewshed(TestCase):
         self.assertRasterMinMax(map=self.viewshed, refmin=0, refmax=180,
             msg="Viewing angle above the ground must be between 0 and 180 deg")
 
+    def test_direction(self):
+        """Test direction range
+        """
+        obs_elev = '1.75'
+        self.assertModule('r.viewshed', input='elevation',
+            coordinates=(634720,216180), output=self.viewshed,
+            observer_elevation=obs_elev, direction_range=[270, 90])
+        minmax = 'null_cells=1998900\nmin=74.90344\nmax=180'
+        self.assertRasterFitsUnivar(raster=self.viewshed, reference=minmax, precision=1e-5)
+
 
 class TestViewshedAgainstReference(TestCase):
     """

raster/r.viewshed/viewshed.cpp

@@ -282,6 +282,11 @@ MemoryVisibilityGrid *viewshed_in_memory(char *inputfname, GridHeader * hd,
 	if (!is_nodata(visgrid->grid->hd, data[1][i]) &&
 	    !is_point_outside_max_dist(*vp, *hd, sn.row, sn.col,
 				       viewOptions.maxDist)) {
+	    if (viewOptions.doDirection &&
+		    !is_point_inside_angle(*vp, sn.row, sn.col,
+					   viewOptions.horizontal_angle_min,
+					   viewOptions.horizontal_angle_max))
+		continue;
 	    /*calculate Distance to VP and Gradient, store them into sn */
 	    /* need either 3 elevation values or 
 	     * 3 gradients calculated from 3 elevation values */
@@ -527,6 +532,11 @@ IOVisibilityGrid *viewshed_external(char *inputfname, GridHeader * hd,
 	if (!is_nodata(visgrid->hd, data[1][i]) &&
 	    !is_point_outside_max_dist(*vp, *hd, sn.row, sn.col,
 				       viewOptions.maxDist)) {
+	    if (viewOptions.doDirection &&
+		    !is_point_inside_angle(*vp, sn.row, sn.col,
+					   viewOptions.horizontal_angle_min,
+					   viewOptions.horizontal_angle_max))
+		continue;
 	    /*calculate Distance to VP and Gradient, store them into sn */
 	    /* need either 3 elevation values or 
 	     * 3 gradients calculated from 3 elevation values */

raster/r.viewshed/visibility.h

@@ -110,6 +110,11 @@ typedef struct viewOptions_
     /* points that are farther than this distance from the viewpoint are
        not visible  */
 
+    float horizontal_angle_min;
+    float horizontal_angle_max;
+    int doDirection;
+    /* exclude points outside of the angle range  */
+
     OutputMode outputMode;
     /* The mode the viewshed is output; 
        - in angle mode, the values recorded are   {NODATA, INVISIBLE, angle}