Support lfp for subwatersheds (non-accumulated lfp by default) (#159)

grass7/raster/r.accumulate/delineate_streams.c

@@ -3,9 +3,8 @@
 #include <grass/glocale.h>
 #include "global.h"
 
-static void trace_down(struct cell_map *, struct raster_map *accum_buf,
-                       double, char, struct Cell_head *, int, int,
-                       struct point_list *);
+static void trace_down(struct cell_map *, struct raster_map *, double, char,
+		       struct Cell_head *, int, int, struct point_list *);
 
 void delineate_streams(struct Map_info *Map, struct cell_map *dir_buf,
                        struct raster_map *accum_buf, double thresh, char conf)
@@ -94,7 +93,6 @@ static void trace_down(struct cell_map *dir_buf, struct raster_map *accum_buf,
         {-1, 1}, {-1, 0}, {-1, -1}, {0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}
     };
     int rows = dir_buf->rows, cols = dir_buf->cols;
-    int i, j;
     int dir;
 
     /* if the current cell is outside the computational region, stop tracing */
@@ -109,6 +107,7 @@ static void trace_down(struct cell_map *dir_buf, struct raster_map *accum_buf,
      * conf = 0 for split streams at a confluence */
     if (!conf) {
         int nup = 0;
+	int i, j;
 
         for (i = -1; i <= 1; i++) {
             /* skip edge cells */

grass7/raster/r.accumulate/global.h

@@ -1,9 +1,15 @@
 #ifndef _GLOBAL_H_
 #define _GLOBAL_H_
 
+#include <float.h>
 #include <grass/raster.h>
 #include <grass/vector.h>
 
+#ifndef DBL_MAX
+#define DBL_MAX 1.797693E308  /* DBL_MAX approximation */
+#endif
+#define OUTLET -DBL_MAX
+
 #define NE 1
 #define N 2
 #define NW 3
@@ -90,6 +96,10 @@ void accumulate(struct cell_map *, struct raster_map *, struct raster_map *,
 void delineate_streams(struct Map_info *, struct cell_map *,
                        struct raster_map *, double, char);
 
+/* subaccumulate.c */
+void subaccumulate(struct Map_info *, struct cell_map *, struct raster_map *,
+		   struct point_list *);
+
 /* calculate_lfp.c */
 void calculate_lfp(struct Map_info *, struct cell_map *, struct raster_map *,
                    int *, char *, struct point_list *);

grass7/raster/r.accumulate/main.c

@@ -41,6 +41,7 @@ int main(int argc, char *argv[])
         struct Option *weight;
         struct Option *input_accum;
         struct Option *accum;
+        struct Option *subaccum;
         struct Option *thresh;
         struct Option *stream;
         struct Option *coords;
@@ -54,16 +55,17 @@ int main(int argc, char *argv[])
     struct
     {
         struct Flag *neg;
+	struct Flag *accum;
         struct Flag *conf;
     } flag;
     char *desc;
     char *dir_name, *weight_name, *input_accum_name, *accum_name,
-        *stream_name, *outlet_name, *lfp_name;
+	 *subaccum_name, *stream_name, *outlet_name, *lfp_name;
     int dir_fd;
     double dir_format, thresh;
     struct Range dir_range;
     CELL dir_min, dir_max;
-    char neg, conf;
+    char neg, accum, conf;
     char **done;
     struct cell_map dir_buf;
     struct raster_map weight_buf, accum_buf;
@@ -118,6 +120,13 @@ int main(int argc, char *argv[])
     opt.accum->description =
         _("Name for output weighted flow accumulation map");
 
+    opt.subaccum = G_define_standard_option(G_OPT_R_OUTPUT);
+    opt.subaccum->key = "subaccumulation";
+    opt.subaccum->required = NO;
+    opt.subaccum->type = TYPE_STRING;
+    opt.subaccum->description =
+        _("Name for output weighted flow subaccumulation map");
+
     opt.thresh = G_define_option();
     opt.thresh->key = "threshold";
     opt.thresh->type = TYPE_DOUBLE;
@@ -168,6 +177,11 @@ int main(int argc, char *argv[])
     flag.neg->label =
         _("Use negative flow accumulation for likely underestimates");
 
+    flag.accum = G_define_flag();
+    flag.accum->key = 'a';
+    flag.accum->label =
+        _("Calculate accumulated longest flow paths");
+
     flag.conf = G_define_flag();
     flag.conf->key = 'c';
     flag.conf->label = _("Delineate streams across confluences");
@@ -193,6 +207,7 @@ int main(int argc, char *argv[])
     weight_name = opt.weight->answer;
     input_accum_name = opt.input_accum->answer;
     accum_name = opt.accum->answer;
+    subaccum_name = opt.subaccum->answer;
     stream_name = opt.stream->answer;
     outlet_name = opt.outlet->answer;
     lfp_name = opt.lfp->answer;
@@ -348,6 +363,7 @@ int main(int argc, char *argv[])
 
     thresh = opt.thresh->answer ? atof(opt.thresh->answer) : 0.0;
     neg = flag.neg->answer;
+    accum = flag.accum->answer;
     conf = flag.conf->answer;
 
     rows = Rast_window_rows();
@@ -438,8 +454,11 @@ int main(int argc, char *argv[])
         struct History hist;
 
         G_message(_("Writing accumulation map..."));
-        for (row = 0; row < rows; row++)
+        for (row = 0; row < rows; row++) {
+            G_percent(row, rows, 1);
             Rast_put_row(accum_fd, accum_buf.map.v[row], accum_buf.type);
+	}
+	G_percent(1, 1, 1);
         Rast_close(accum_fd);
 
         /* write history */
@@ -468,6 +487,35 @@ int main(int argc, char *argv[])
         Vect_close(&Map);
     }
 
+    /* calculate subaccumulation if needed */
+    if (subaccum_name || (lfp_name && !accum)) {
+	subaccumulate(&Map, &dir_buf, &accum_buf, &outlet_pl);
+
+	/* write out buffer to the subaccumulation map if requested */
+	if (subaccum_name) {
+	    int subaccum_fd = Rast_open_new(subaccum_name, accum_buf.type);
+	    struct History hist;
+
+	    G_message(_("Writing subaccumulation map..."));
+	    for (row = 0; row < rows; row++) {
+		G_percent(row, rows, 1);
+		Rast_put_row(subaccum_fd, accum_buf.map.v[row], accum_buf.type);
+	    }
+	    G_percent(1, 1, 1);
+	    Rast_close(subaccum_fd);
+
+	    /* write history */
+	    Rast_put_cell_title(subaccum_name,
+				weight_name ? "Weighted flow subaccumulation" :
+				(neg ?
+				 "Flow subaccumulation with likely underestimates" :
+				 "Flow subaccumulation"));
+	    Rast_short_history(subaccum_name, "raster", &hist);
+	    Rast_command_history(&hist);
+	    Rast_write_history(subaccum_name, &hist);
+	}
+    }
+
     /* calculate the longest flow path */
     if (lfp_name) {
         if (Vect_open_new(&Map, lfp_name, 0) < 0)

grass7/raster/r.accumulate/subaccumulate.c

@@ -0,0 +1,65 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+static void trace_down(struct cell_map *, struct raster_map *, double,
+                       struct Cell_head *, int, int);
+
+void subaccumulate(struct Map_info *Map, struct cell_map *dir_buf,
+		   struct raster_map *accum_buf, struct point_list *outlet_pl)
+{
+    struct Cell_head window;
+    int rows = accum_buf->rows, cols = accum_buf->cols;
+    int i;
+
+    G_get_set_window(&window);
+
+    /* loop through all outlets and calculate subaccumulation for each */
+    G_message(_("Subaccumulating flows..."));
+    for (i = 0; i < outlet_pl->n; i++) {
+        int row = (int)Rast_northing_to_row(outlet_pl->y[i], &window);
+        int col = (int)Rast_easting_to_col(outlet_pl->x[i], &window);
+
+        G_percent(i, outlet_pl->n, 1);
+
+        /* if the outlet is outside the computational region, skip */
+        if (row < 0 || row >= rows || col < 0 || col >= cols) {
+            G_warning(_("Skip outlet (%f, %f) outside the current region"),
+                      outlet_pl->x[i], outlet_pl->y[i]);
+            continue;
+        }
+
+        /* trace down flow accumulation and calculate subaccumulation */
+        trace_down(dir_buf, accum_buf, OUTLET, &window, row, col);
+    }
+    G_percent(1, 1, 1);
+}
+
+static void trace_down(struct cell_map *dir_buf, struct raster_map *accum_buf,
+		       double up_accum, struct Cell_head *window,
+		       int row, int col)
+{
+    static int next_cells[8][2] = {
+        {-1, 1}, {-1, 0}, {-1, -1}, {0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}
+    };
+    int rows = dir_buf->rows, cols = dir_buf->cols;
+    int dir;
+
+    /* if the current cell is outside the computational region, stop tracing */
+    if (row < 0 || row >= rows || col < 0 || col >= cols)
+        return;
+
+    if (up_accum == OUTLET)
+	/* accumulation at the outlet will need to be subtracted from all
+	 * downstream accumulation cells */
+	up_accum = get(accum_buf, row, col);
+    else
+	/* calculate subaccumulation */
+	set(accum_buf, row, col, get(accum_buf, row, col) - up_accum);
+
+    /* if the current cell doesn't flow out of the computational region
+     * (negative direction from r.watershed flows out), keep tracing */
+    dir = dir_buf->c[row][col] - 1;
+    if (dir >= 0 && dir < 8)
+        trace_down(dir_buf, accum_buf, up_accum, window,
+                   row + next_cells[dir][0], col + next_cells[dir][1]);
+}