Set null outside direction map

grass7/raster/r.accumulate/accumulate.c

@@ -14,7 +14,10 @@ void accumulate(struct cell_map *dir_buf, struct raster_map *weight_buf,
     for (row = 0; row < rows; row++) {
         G_percent(row, rows, 1);
         for (col = 0; col < cols; col++)
-            trace_up(dir_buf, weight_buf, accum_buf, done, neg, row, col);
+            if (Rast_is_c_null_value(&dir_buf->c[row][col]))
+                set_null(accum_buf, row, col);
+            else
+                trace_up(dir_buf, weight_buf, accum_buf, done, neg, row, col);
     }
     G_percent(1, 1, 1);
 }

grass7/raster/r.accumulate/global.h

@@ -75,6 +75,8 @@ GLOBAL int dir_checks[3][3][2]
 /* raster.c */
 void set(struct raster_map *, int, int, double);
 double get(struct raster_map *, int, int);
+int is_null(struct raster_map *, int, int);
+void set_null(struct raster_map *, int, int);
 
 /* point_list.c */
 void init_point_list(struct point_list *);

grass7/raster/r.accumulate/main.c

@@ -191,8 +191,9 @@ int main(int argc, char *argv[])
     G_option_exclusive(opt.weight, opt.lfp, flag.neg, NULL);
     G_option_exclusive(opt.input_accum, opt.accum, NULL);
     G_option_exclusive(opt.weight, opt.input_accum, NULL);
-    G_option_required(opt.accum, opt.stream, opt.lfp, NULL);
+    G_option_required(opt.accum, opt.subaccum, opt.stream, opt.lfp, NULL);
     G_option_collective(opt.thresh, opt.stream, NULL);
+    G_option_requires(opt.subaccum, opt.coords, opt.outlet, NULL);
     G_option_requires(opt.lfp, opt.coords, opt.outlet, NULL);
     G_option_requires(opt.idcol, opt.id, opt.outlet_idcol, NULL);
     G_option_requires_all(opt.id, opt.idcol, opt.coords, NULL);

grass7/raster/r.accumulate/raster.c

@@ -14,8 +14,6 @@ void set(struct raster_map *buf, int row, int col, double value)
         buf->map.d[row][col] = (DCELL) value;
         break;
     }
-
-    return;
 }
 
 double get(struct raster_map *buf, int row, int col)
@@ -36,3 +34,37 @@ double get(struct raster_map *buf, int row, int col)
 
     return value;
 }
+
+int is_null(struct raster_map *buf, int row, int col)
+{
+    int is_null_value = 0;
+
+    switch (buf->type) {
+    case CELL_TYPE:
+        is_null_value = Rast_is_c_null_value(&buf->map.c[row][col]);
+        break;
+    case FCELL_TYPE:
+        is_null_value = Rast_is_f_null_value(&buf->map.f[row][col]);
+        break;
+    case DCELL_TYPE:
+        is_null_value = Rast_is_d_null_value(&buf->map.d[row][col]);
+        break;
+    }
+
+    return is_null_value;
+}
+
+void set_null(struct raster_map *buf, int row, int col)
+{
+    switch (buf->type) {
+    case CELL_TYPE:
+        Rast_set_c_null_value(&buf->map.c[row][col], 1);
+        break;
+    case FCELL_TYPE:
+        Rast_set_f_null_value(&buf->map.f[row][col], 1);
+        break;
+    case DCELL_TYPE:
+        Rast_set_d_null_value(&buf->map.d[row][col], 1);
+        break;
+    }
+}

grass7/raster/r.accumulate/subaccumulate.c

@@ -52,9 +52,16 @@ static void trace_down(struct cell_map *dir_buf, struct raster_map *accum_buf,
         /* accumulation at the outlet will need to be subtracted from all
          * downstream accumulation cells */
         up_acc = get(accum_buf, row, col);
-    else
+    else {
         /* calculate subaccumulation */
-        set(accum_buf, row, col, get(accum_buf, row, col) - up_acc);
+        double acc = get(accum_buf, row, col);
+
+        if (acc <= up_acc)
+            /* downstream is already processed */
+            return;
+
+        set(accum_buf, row, col, acc - up_acc);
+    }
 
     /* if the current cell doesn't flow out of the computational region
      * (negative direction from r.watershed flows out), keep tracing */