r.watershed: Can now optionally use already existing maps to optimize performance

raster/r.watershed/front/main.c

@@ -61,6 +61,8 @@ int main(int argc, char *argv[])
     struct Option *opt17;
     struct Option *opt18;
     struct Option *opt19;
+    struct Option *opt20; /* Input accumulation map */
+    struct Option *opt21; /* Input drainage direction map */
     struct Flag *flag_sfd;
     struct Flag *flag_flow;
     struct Flag *flag_seg;
@@ -123,6 +125,31 @@ int main(int argc, char *argv[])
     opt19->required = NO;
     opt19->guisection = _("Inputs");
 
+    opt20 = G_define_standard_option(G_OPT_R_INPUT);
+    opt20->key = "accumulation_input";
+    opt20->label = _("Name of existing accumulation raster map to reuse");
+    opt20->description =
+        _("Reuses existing flow accumulation map instead of recalculating"
+          "Must be from r.watershed output with the same elevation. "
+          "Skips flow calculation for faster basin delineation with different "
+          "thresholds. "
+          "Requires drainage_input to also be specified.");
+    opt20->required = NO;
+    opt20->guisection = _("Inputs");
+
+    opt21 = G_define_standard_option(G_OPT_R_INPUT);
+    opt21->key = "drainage_input";
+    opt21->label = _("Name of existing drainage direction raster map to reuse");
+    opt21->description =
+        _("Reuses existing drainage direction map instead of recalculating"
+          "MUST be from r.watershed (not r.stream.* or r.terraflow). "
+          "Uses r.watershed conventions: values 0-8 for directions, negative "
+          "for outlets. "
+          "Using incompatible maps will produce incorrect results. "
+          "Requires accumulation_input to also be specified.");
+    opt21->required = NO;
+    opt21->guisection = _("Inputs");
+
     opt6 = G_define_option();
     opt6->key = "threshold";
     opt6->description = _("Minimum size of exterior watershed basin");
@@ -247,13 +274,29 @@ int main(int argc, char *argv[])
     G_option_requires(opt13, opt6, NULL);
     G_option_requires(opt14, opt6, NULL);
 
+    G_option_collective(opt20, opt21, NULL);
+
+    /* Cannot use input maps with certain incompatible options */
+    G_option_excludes(
+        opt20, opt2, opt3, opt4, opt5, opt19,
+        NULL); /* depression, flow, disturbed_land, blocking, retention */
+    G_option_excludes(opt21, opt2, opt3, opt4, opt5, opt19, NULL);
+
     /* Check for some output map */
     G_option_required(opt8, opt17, opt18, opt9, opt10, opt11, opt12, opt13,
                       opt14, NULL);
 
     if (G_parser(argc, argv))
         exit(EXIT_FAILURE);
 
+    /* Prevent basin delineation with input maps */
+    if ((opt20->answer || opt21->answer) &&
+        (opt10->answer || opt11->answer || opt12->answer)) {
+        G_fatal_error(
+            _("Basin delineation (basin=, stream=, half_basin=) is not "
+              "supported when using accumulation_input and drainage_input. "
+              "To generate basins, run r.watershed without input maps."));
+    }
     /* basin threshold */
     if (opt6->answer) {
         if (atoi(opt6->answer) <= 0)
@@ -287,6 +330,8 @@ int main(int argc, char *argv[])
     do_opt(opt4);
     do_opt(opt5);
     do_opt(opt19);
+    do_opt(opt20);
+    do_opt(opt21);
     do_opt(opt6);
     do_opt(opt7);
     do_opt(opt8);

raster/r.watershed/ram/Gwater.h

@@ -84,6 +84,8 @@ extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX],
 extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX],
     spi_name[GNAME_MAX];
 extern char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+extern char wat_input_name[GNAME_MAX], asp_input_name[GNAME_MAX];
+extern char wat_input_flag, asp_input_flag;
 extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag,
     rtn_flag;
 extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag,

raster/r.watershed/ram/init_vars.c

@@ -1,5 +1,6 @@
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #include "Gwater.h"
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -21,6 +22,7 @@ int init_vars(int argc, char *argv[])
     G_gisinit(argv[0]);
     /* input */
     ele_flag = pit_flag = run_flag = ril_flag = rtn_flag = 0;
+    wat_input_flag = asp_input_flag = 0;
     /* output */
     wat_flag = asp_flag = tci_flag = spi_flag = atanb_flag = 0;
     bas_flag = seg_flag = haf_flag = 0;
@@ -49,6 +51,10 @@ int init_vars(int argc, char *argv[])
             ele_flag++;
         else if (sscanf(argv[r], "accumulation=%s", wat_name) == 1)
             wat_flag++;
+        else if (sscanf(argv[r], "accumulation_input=%s", wat_input_name) == 1)
+            wat_input_flag++;
+        else if (sscanf(argv[r], "drainage_input=%s", asp_input_name) == 1)
+            asp_input_flag++;
         else if (sscanf(argv[r], "tci=%s", tci_name) == 1)
             tci_flag++;
         else if (sscanf(argv[r], "spi=%s", spi_name) == 1)
@@ -105,6 +111,19 @@ int init_vars(int argc, char *argv[])
         else
             usage(argv[0]);
     }
+    if (wat_input_flag && asp_input_flag) {
+        if (wat_input_flag != 1 || asp_input_flag != 1) {
+            G_fatal_error(_("Both accumulation_input and drainage_input must "
+                            "be specified together"));
+        }
+        /* Cannot use certain options with input maps */
+        if (pit_flag || run_flag || rtn_flag) {
+            G_fatal_error(_("Cannot use depression, flow, or retention maps "
+                            "with accumulation_input and drainage_input"));
+        }
+        G_message(_("Using existing accumulation and drainage maps - skipping "
+                    "flow calculation"));
+    }
     if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
         G_fatal_error("Convergence factor must be between 1 and 10.");
     }
@@ -120,6 +139,11 @@ int init_vars(int argc, char *argv[])
         tot_parts++;
     if (bas_thres > 0)
         tot_parts++;
+
+    /* Reduce total parts if skipping flow calculation */
+    if (wat_input_flag && asp_input_flag) {
+        tot_parts -= 2; /* Skip A* and accumulation sections */
+    }
     G_message(_("SECTION 1a (of %1d): Initiating Memory."), tot_parts);
     this_mapset = G_mapset();
     if (sl_flag || sg_flag || ls_flag)
@@ -170,6 +194,150 @@ int init_vars(int argc, char *argv[])
     in_list = flag_create(nrows, ncols);
     worked = flag_create(nrows, ncols);
 
+    if (wat_input_flag && asp_input_flag) {
+        int input_fd;
+        DCELL *wat_buf;
+        CELL *asp_buf;
+
+        G_message(_("Reading existing accumulation map: %s"), wat_input_name);
+        input_fd = Rast_open_old(wat_input_name, "");
+        wat_buf = Rast_allocate_d_buf();
+
+        for (r = 0; r < nrows; r++) {
+            G_percent(r, nrows, 2);
+            Rast_get_d_row(input_fd, wat_buf, r);
+            for (c = 0; c < ncols; c++) {
+                seg_idx = SEG_INDEX(wat_seg, r, c);
+                if (Rast_is_d_null_value(&wat_buf[c])) {
+                    Rast_set_d_null_value(&wat[seg_idx], 1);
+                    FLAG_SET(worked, r, c);
+                    FLAG_SET(in_list, r, c);
+                }
+                else {
+                    wat[seg_idx] = wat_buf[c];
+                }
+            }
+        }
+        G_percent(nrows, nrows, 2);
+        Rast_close(input_fd);
+        G_free(wat_buf);
+
+        G_message(_("Reading existing drainage direction map: %s"),
+                  asp_input_name);
+        input_fd = Rast_open_old(asp_input_name, "");
+        asp_buf = Rast_allocate_c_buf();
+
+        for (r = 0; r < nrows; r++) {
+            G_percent(r, nrows, 2);
+            Rast_get_c_row(input_fd, asp_buf, r);
+            for (c = 0; c < ncols; c++) {
+                seg_idx = SEG_INDEX(asp_seg, r, c);
+                if (Rast_is_c_null_value(&asp_buf[c])) {
+                    Rast_set_c_null_value(&asp[seg_idx], 1);
+                }
+                else {
+                    asp[seg_idx] = asp_buf[c];
+                }
+            }
+        }
+        G_percent(nrows, nrows, 2);
+        Rast_close(input_fd);
+        G_free(asp_buf);
+
+        /* Read elevation for other calculations */
+        fd = Rast_open_old(ele_name, "");
+        ele_map_type = Rast_get_map_type(fd);
+        ele_size = Rast_cell_size(ele_map_type);
+        elebuf = Rast_allocate_buf(ele_map_type);
+
+        if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
+            ele_scale = 1000;
+        if (flat_flag)
+            ele_scale = 10000;
+
+        do_points = 0; /* No A* processing needed */
+        for (r = 0; r < nrows; r++) {
+            Rast_get_row(fd, elebuf, r, ele_map_type);
+            ptr = elebuf;
+            for (c = 0; c < ncols; c++) {
+                seg_idx = SEG_INDEX(alt_seg, r, c);
+
+                if (!Rast_is_null_value(ptr, ele_map_type)) {
+                    if (ele_map_type == CELL_TYPE) {
+                        alt_value = *((CELL *)ptr);
+                        alt_value *= ele_scale;
+                    }
+                    else if (ele_map_type == FCELL_TYPE) {
+                        dvalue = *((FCELL *)ptr);
+                        dvalue *= ele_scale;
+                        alt_value = ele_round(dvalue);
+                    }
+                    else if (ele_map_type == DCELL_TYPE) {
+                        dvalue = *((DCELL *)ptr);
+                        dvalue *= ele_scale;
+                        alt_value = ele_round(dvalue);
+                    }
+                    alt[seg_idx] = alt_value;
+                    if (er_flag) {
+                        r_h[seg_idx] = alt_value;
+                    }
+                }
+                else {
+                    Rast_set_c_null_value(&alt[seg_idx], 1);
+                    if (er_flag) {
+                        Rast_set_c_null_value(&r_h[seg_idx], 1);
+                    }
+                }
+
+                if (atanb_flag) {
+                    Rast_set_d_null_value(&sca[seg_idx], 1);
+                    Rast_set_d_null_value(&tanb[seg_idx], 1);
+                }
+                ptr = G_incr_void_ptr(ptr, ele_size);
+            }
+        }
+        Rast_close(fd);
+        G_free(elebuf);
+
+        /*  Initialize RUSLE s_l array when using input maps */
+        if (er_flag) {
+            s_l = (double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
+                                     sizeof(double));
+            /* Initialize s_l with half_res for all cells */
+            for (r = 0; r < nrows; r++) {
+                for (c = 0; c < ncols; c++) {
+                    seg_idx = SEG_INDEX(s_l_seg, r, c);
+                    s_l[seg_idx] = half_res;
+                }
+            }
+
+            if (sg_flag) {
+                s_g = (double *)G_malloc(size_array(&s_g_seg, nrows, ncols) *
+                                         sizeof(double));
+            }
+            if (ls_flag) {
+                l_s = (double *)G_malloc(size_array(&l_s_seg, nrows, ncols) *
+                                         sizeof(double));
+            }
+        }
+        /* Mark swales from drainage direction if needed for basin delineation
+         */
+        if (bas_thres > 0) {
+            for (r = 0; r < nrows; r++) {
+                for (c = 0; c < ncols; c++) {
+                    seg_idx = SEG_INDEX(wat_seg, r, c);
+                    if (!Rast_is_d_null_value(&wat[seg_idx]) &&
+                        fabs(wat[seg_idx]) >= bas_thres) {
+                        FLAG_SET(swale, r, c);
+                    }
+                }
+            }
+        }
+
+        G_message(_("Skipping flow calculation - using existing maps"));
+        return 0;
+    }
+
     /* open elevation input */
     fd = Rast_open_old(ele_name, "");
 
@@ -238,8 +406,8 @@ int init_vars(int argc, char *argv[])
     G_free(elebuf);
     MASK_flag = (do_points < (size_t)nrows * ncols);
 
-    /* read flow accumulation from input map flow: amount of overland flow per
-     * cell */
+    /* read flow accumulation from input map flow: amount of overland flow
+     * per cell */
     if (run_flag) {
         dbuf = Rast_allocate_d_buf();
         fd = Rast_open_old(run_name, "");

raster/r.watershed/ram/main.c

@@ -60,6 +60,8 @@ char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX],
 char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX],
     spi_name[GNAME_MAX];
 char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+char wat_input_name[GNAME_MAX], asp_input_name[GNAME_MAX];
+char wat_input_flag = 0, asp_input_flag = 0;
 char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
 char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, rtn_flag;
 char bas_flag, seg_flag, haf_flag, er_flag, tci_flag, spi_flag, atanb_flag;
@@ -69,13 +71,21 @@ FILE *fp;
 int main(int argc, char *argv[])
 {
     init_vars(argc, argv);
-    do_astar();
-    if (mfd) {
-        do_cum_mfd();
+
+    if (!wat_input_flag || !asp_input_flag) {
+        do_astar();
+        if (mfd) {
+            do_cum_mfd();
+        }
+        else {
+            do_cum();
+        }
     }
     else {
-        do_cum();
+        G_message(
+            _("Skipping A* search and flow accumulation - using input maps"));
     }
+
     if (sg_flag || ls_flag) {
         sg_factor();
     }

raster/r.watershed/seg/Gwater.h

@@ -115,7 +115,9 @@ extern int nextdc[8];
 extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 extern char ril_name[GNAME_MAX], rtn_name[GNAME_MAX], dep_name[GNAME_MAX];
-
+/* ADDED: Input map names for reusing existing accumulation and drainage */
+extern char wat_input_name[GNAME_MAX], asp_input_name[GNAME_MAX];
+extern char wat_input_flag, asp_input_flag;
 extern const char *this_mapset;
 extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
     thr_name[8];