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

@@ -83,6 +83,12 @@ <h2>NOTES</h2>
 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

raster/r.viewshed/testsuite/test_r_viewshed.py

@@ -60,6 +60,16 @@ def test_limits_extreme_value(self):
         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}