r.watershed: Handle large datasets greater than INT32_MAX in cells

raster/r.watershed/ram/Gwater.h

@@ -48,8 +48,8 @@ OC_STACK
 extern struct Cell_head window;
 
 extern int mfd, c_fac, abs_acc, ele_scale;
-extern int *heap_index, heap_size;
-extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+extern size_t *heap_index, heap_size;
+extern size_t first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 extern int nrows, ncols;
 extern double half_res, diag, max_length, dep_slope;
 extern int bas_thres, tot_parts;
@@ -60,7 +60,7 @@ extern FLAG *worked, *in_list, *s_b, *swale, *flat_done;
 extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
 extern RAMSEG r_h_seg, dep_seg, rtn_seg;
 extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
-extern int *astar_pts;
+extern size_t *astar_pts;
 extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
 extern char *rtn;
 extern DCELL *wat, *sca, *tanb;
@@ -108,7 +108,7 @@ int drop_pt(void);
 double get_slope(int, int, int, int, CELL, CELL);
 
 /* do_flatarea.c */
-int do_flatarea(int, CELL, CELL *, CELL *);
+int do_flatarea(size_t, CELL, CELL *, CELL *);
 
 /* do_cum.c */
 int do_cum(void);
@@ -131,8 +131,8 @@ int no_stream(int, int, CELL, double, CELL);
 int overland_cells(int, int, CELL, CELL, CELL *);
 
 /* ramseg.c */
-int size_array(int *, int, int);
-int seg_index_rc(int, int, int *, int *);
+size_t size_array(int *, int, int);
+size_t seg_index_rc(int, size_t, int *, int *);
 
 /* sg_factor.c */
 int sg_factor(void);

raster/r.watershed/ram/do_astar.c

@@ -7,11 +7,11 @@ static double get_slope2(CELL, CELL, double);
 
 int do_astar(void)
 {
-    int count;
+    size_t count;
     int upr, upc, r, c, ct_dir;
     CELL alt_val, alt_nbr[8];
     CELL is_in_list, is_worked, flat_is_done, nbr_flat_is_done;
-    int index_doer, index_up;
+    size_t index_doer, index_up;
 
     /* sides
      * |7|1|4|
@@ -202,7 +202,7 @@ int do_astar(void)
 
 /* compare two heap points */
 /* return 1 if a < b else 0 */
-static int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
+static int cmp_pnt(CELL elea, CELL eleb, size_t addeda, size_t addedb)
 {
     if (elea == eleb) {
         return (addeda < addedb);
@@ -211,9 +211,9 @@ static int cmp_pnt(CELL elea, CELL eleb, int addeda, int addedb)
 }
 
 /* standard sift-up routine for d-ary min heap */
-static int sift_up(int start, CELL ele)
+static int sift_up(size_t start, CELL ele)
 {
-    register int parent, child, child_idx, child_added;
+    register size_t parent, child, child_idx, child_added;
     CELL elep;
 
     child = start;
@@ -268,9 +268,9 @@ int add_pt(int r, int c, CELL ele)
 /* drop point routine for min heap */
 int drop_pt(void)
 {
-    register int child, childr, parent;
+    register size_t child, childr, parent;
     CELL ele, eler;
-    register int i;
+    register size_t i;
 
     if (heap_size == 1) {
         heap_index[1] = -1;

raster/r.watershed/ram/do_cum.c

@@ -105,10 +105,11 @@ int do_cum(void)
     int r_nbr, c_nbr, ct_dir, np_side, edge;
     CELL is_swale, aspect, ele_nbr;
     DCELL value, valued;
-    int killer, threshold;
+    size_t killer;
+    int threshold;
     int asp_r[9] = {0, -1, -1, -1, 0, 1, 1, 1, 0};
     int asp_c[9] = {0, 1, 0, -1, -1, -1, 0, 1, 1};
-    int this_index, down_index, nbr_index;
+    size_t this_index, down_index, nbr_index;
     double *dist_to_nbr, *contour;
     double cell_size;
 
@@ -262,9 +263,10 @@ int do_cum(void)
 int do_cum_mfd(void)
 {
     int r, c, dr, dc;
-    CELL is_swale;
+    CELL is_swale = 0;
     DCELL value, valued, tci_div, sum_contour, cell_size;
-    int killer, threshold;
+    size_t killer;
+    int threshold;
 
     /* MFD */
     int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
@@ -275,7 +277,7 @@ int do_cum_mfd(void)
     int workedon, edge, flat;
     int asp_r[9] = {0, -1, -1, -1, 0, 1, 1, 1, 0};
     int asp_c[9] = {0, 1, 0, -1, -1, -1, 0, 1, 1};
-    int this_index, down_index, nbr_index;
+    size_t this_index, down_index, nbr_index;
 
     /* drainage directions bitmask encoded CW from North
      * drainage directions are set for each current cell
@@ -499,9 +501,9 @@ int do_cum_mfd(void)
     }
     if (workedon)
         G_warning(n_("MFD: A * path already processed when distributing flow: "
-                     "%d of %d cell",
+                     "%d of %zu cell",
                      "MFD: A * path already processed when distributing flow: "
-                     "%d of %d cells",
+                     "%d of %zu cells",
                      do_points),
                   workedon, do_points);
 

raster/r.watershed/ram/do_flatarea.c

@@ -13,20 +13,21 @@
 
 #include <limits.h>
 #include <assert.h>
+#include <stdint.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include <grass/rbtree.h>
 #include "Gwater.h"
 #include "do_astar.h"
 
 struct pq_node {
-    int idx;
+    size_t idx;
     struct pq_node *next;
 };
 
 struct pq {
     struct pq_node *first, *last;
-    int size;
+    size_t size;
 };
 
 struct pq *pq_create(void)
@@ -35,18 +36,18 @@ struct pq *pq_create(void)
 
     q->first = G_malloc(sizeof(struct pq_node));
     q->first->next = NULL;
-    q->first->idx = -1;
+    q->first->idx = SIZE_MAX;
     q->last = q->first;
     q->size = 0;
 
     return q;
 }
 
 /* dummy end must always be allocated and empty */
-int pq_add(int idx, struct pq *q)
+int pq_add(size_t idx, struct pq *q)
 {
     assert(q->last);
-    assert(q->last->idx == -1);
+    assert(q->last->idx == SIZE_MAX);
 
     q->last->idx = idx;
     if (q->last->next != NULL) {
@@ -56,7 +57,7 @@ int pq_add(int idx, struct pq *q)
     struct pq_node *n = (struct pq_node *)G_malloc(sizeof(struct pq_node));
 
     n->next = NULL;
-    n->idx = -1;
+    n->idx = SIZE_MAX;
     q->last->next = n;
     q->last = q->last->next;
 
@@ -99,7 +100,8 @@ int pq_destroy(struct pq *q)
 }
 
 struct orders {
-    int index, uphill, downhill;
+    size_t index;
+    int uphill, downhill;
     char flag;
 };
 
@@ -115,11 +117,11 @@ int cmp_orders(const void *a, const void *b)
  * return 0 if nothing was modidied
  * return 1 if elevation was modified
  */
-int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
+int do_flatarea(size_t index, CELL ele, CELL *alt_org, CELL *alt_new)
 {
     int upr, upc, r, c, ct_dir;
     CELL is_in_list, is_worked, this_in_list;
-    int index_doer, index_up;
+    size_t index_doer, index_up;
     int n_flat_cells = 0, counter;
     CELL ele_nbr, min_ele_diff;
     int uphill_order, downhill_order, max_uphill_order, max_downhill_order;
@@ -149,7 +151,7 @@ int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
     G_debug(2, "get uphill start points");
     counter = 0;
     while (down_pq->size) {
-        if ((index_doer = pq_drop(down_pq)) == -1)
+        if ((index_doer = pq_drop(down_pq)) == SIZE_MAX)
             G_fatal_error("get start points: no more points in down queue");
 
         seg_index_rc(alt_seg, index_doer, &r, &c);
@@ -218,7 +220,7 @@ int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
     while (up_pq->size) {
         int is_in_down_queue = 0;
 
-        if ((index_doer = pq_drop(up_pq)) == -1)
+        if ((index_doer = pq_drop(up_pq)) == SIZE_MAX)
             G_fatal_error("uphill order: no more points in up queue");
 
         seg_index_rc(alt_seg, index_doer, &r, &c);
@@ -302,7 +304,7 @@ int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
     G_debug(2, "got downhill start points, do downhill correction");
     downhill_order = 1;
     while (down_pq->size) {
-        if ((index_doer = pq_drop(down_pq)) == -1)
+        if ((index_doer = pq_drop(down_pq)) == SIZE_MAX)
             G_fatal_error(_("downhill order: no more points in down queue"));
 
         seg_index_rc(alt_seg, index_doer, &r, &c);
@@ -378,7 +380,7 @@ int do_flatarea(int index, CELL ele, CELL *alt_org, CELL *alt_new)
 
     G_debug(2, "adjust ele");
     while (up_pq->size) {
-        if ((index_doer = pq_drop(up_pq)) == -1)
+        if ((index_doer = pq_drop(up_pq)) == SIZE_MAX)
             G_fatal_error("no more points in up queue");
 
         seg_index_rc(alt_seg, index_doer, &r, &c);

raster/r.watershed/ram/init_vars.c

@@ -16,7 +16,7 @@ int init_vars(int argc, char *argv[])
     int fd, ele_map_type;
     size_t ele_size;
     char MASK_flag;
-    int seg_idx;
+    size_t seg_idx;
 
     G_gisinit(argv[0]);
     /* input */
@@ -135,7 +135,7 @@ int init_vars(int argc, char *argv[])
     G_get_set_window(&window);
     nrows = Rast_window_rows();
     ncols = Rast_window_cols();
-    total_cells = nrows * ncols;
+    total_cells = (size_t)nrows * ncols;
     if (max_length <= d_zero)
         max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
     if (window.ew_res < window.ns_res)
@@ -184,7 +184,7 @@ int init_vars(int argc, char *argv[])
 
     /* read elevation input and mark NULL/masked cells */
     /* initialize accumulation and drainage direction */
-    do_points = nrows * ncols;
+    do_points = (size_t)nrows * ncols;
     for (r = 0; r < nrows; r++) {
         Rast_get_row(fd, elebuf, r, ele_map_type);
         ptr = elebuf;
@@ -236,7 +236,7 @@ int init_vars(int argc, char *argv[])
     }
     Rast_close(fd);
     G_free(elebuf);
-    MASK_flag = (do_points < nrows * ncols);
+    MASK_flag = (do_points < (size_t)nrows * ncols);
 
     /* read flow accumulation from input map flow: amount of overland flow per
      * cell */
@@ -331,11 +331,11 @@ int init_vars(int argc, char *argv[])
                                  sizeof(double));
     }
 
-    astar_pts = (int *)G_malloc((do_points + 1) * sizeof(int));
+    astar_pts = (size_t *)G_malloc((do_points + 1) * sizeof(size_t));
 
     /* heap_index will track astar_pts in ternary min-heap */
     /* heap_index is one-based */
-    heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
+    heap_index = (size_t *)G_malloc((do_points + 1) * sizeof(size_t));
 
     G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
 
@@ -348,7 +348,7 @@ int init_vars(int argc, char *argv[])
     }
     else
         buf = NULL;
-    first_astar = first_cum = -1;
+    first_astar = first_cum = 0;
     for (r = 0; r < nrows; r++) {
         G_percent(r, nrows, 3);
         if (pit_flag)

raster/r.watershed/ram/main.c

@@ -24,8 +24,8 @@
 struct Cell_head window;
 
 int mfd, c_fac, abs_acc, ele_scale;
-int *heap_index, heap_size;
-int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+size_t *heap_index, heap_size;
+size_t first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
 int nrows, ncols;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
@@ -36,7 +36,7 @@ FLAG *worked, *in_list, *s_b, *swale, *flat_done;
 RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
 RAMSEG r_h_seg, dep_seg, rtn_seg;
 RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
-int *astar_pts;
+size_t *astar_pts;
 CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
 char *rtn;
 DCELL *wat, *sca, *tanb;

raster/r.watershed/ram/ramseg.c

@@ -1,9 +1,10 @@
 #include <stdio.h>
 #include "ramseg.h"
 
-int size_array(int *ram_seg, int nrows, int ncols)
+size_t size_array(int *ram_seg, int nrows_in, int ncols_in)
 {
-    int size, segs_in_col;
+    size_t size, segs_in_col;
+    size_t nrows = nrows_in, ncols = ncols_in;
 
     segs_in_col = ((nrows - 1) >> RAMSEGBITS) + 1;
     *ram_seg = ((ncols - 1) >> RAMSEGBITS) + 1;
@@ -15,9 +16,9 @@ int size_array(int *ram_seg, int nrows, int ncols)
 }
 
 /* get r, c from seg_index */
-int seg_index_rc(int ramseg, int seg_index, int *r, int *c)
+size_t seg_index_rc(int ramseg, size_t seg_index, int *r, int *c)
 {
-    int seg_no, seg_remainder;
+    size_t seg_no, seg_remainder;
 
     seg_no = seg_index >> DOUBLEBITS;
     seg_remainder = seg_index - (seg_no << DOUBLEBITS);

raster/r.watershed/ram/ramseg.h

@@ -6,8 +6,10 @@
 #define DOUBLEBITS 8  /* 2 * ramsegbits       */
 #define SEGLENLESS 15 /* 2 ^ ramsegbits - 1   */
 
-#define SEG_INDEX(s, r, c)                                                    \
-    (int)(((((r) >> RAMSEGBITS) * (s) + ((c) >> RAMSEGBITS)) << DOUBLEBITS) + \
-          (((r)&SEGLENLESS) << RAMSEGBITS) + ((c)&SEGLENLESS))
+#define SEG_INDEX(s, r, c)                                                 \
+    (size_t)(                                                              \
+        ((((size_t)(r) >> RAMSEGBITS) * (s) + ((size_t)(c) >> RAMSEGBITS)) \
+         << DOUBLEBITS) +                                                  \
+        (((size_t)(r)&SEGLENLESS) << RAMSEGBITS) + ((size_t)(c)&SEGLENLESS))
 
 #endif /* __RAMSEG_H__ */

raster/r.watershed/seg/do_cum.c

@@ -293,7 +293,7 @@ int do_cum_mfd(void)
     int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side /* , max_side */;
     CELL ele, *ele_nbr;
     double prop, max_val;
-    int workedon, edge, is_swale, flat;
+    int workedon, edge, is_swale = 0, flat;
     char *flag_nbr;
     int asp_r[9] = {0, -1, -1, -1, 0, 1, 1, 1, 0};
     int asp_c[9] = {0, 1, 0, -1, -1, -1, 0, 1, 1};