r.watersheds: New addon for OpenMP parallel watershed delineation based on Cho (2025)

src/raster/r.watersheds/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.watersheds
+
+LIBES = $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB)
+EXTRA_INC = $(VECT_INC) $(OPENMP_INCPATH)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(OPENMP_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

src/raster/r.watersheds/delineate.c

@@ -0,0 +1,10 @@
+#include "global.h"
+
+void delineate(struct raster_map *dir_map, struct outlet_list *outlet_l,
+               int use_lessmem)
+{
+    if (use_lessmem)
+        delineate_lessmem(dir_map, outlet_l);
+    else
+        delineate_moremem(dir_map, outlet_l);
+}

src/raster/r.watersheds/delineate_funcs.h

@@ -0,0 +1,223 @@
+#include <stdlib.h>
+#include "global.h"
+
+#define SHED(row, col) DIR(row, col)
+
+#ifdef USE_LESS_MEMORY
+#define DELINEATE         delineate_lessmem
+#define GET_DIR(row, col) (DIR(row, col) & 0x7fffffff)
+#define SET_NOTDONE(row, col)        \
+    do {                             \
+        DIR(row, col) |= 0x80000000; \
+    } while (0)
+#define SET_DONE(row, col)           \
+    do {                             \
+        DIR(row, col) &= 0x7fffffff; \
+    } while (0)
+#define IS_NOTDONE(row, col) (DIR(row, col) & 0x80000000)
+#define IS_DONE(row, col)    !IS_NOTDONE(row, col)
+#else
+#define DELINEATE         delineate_moremem
+#define GET_DIR(row, col) DIR(row, col)
+#define DONE(row, col)    done[INDEX(row, col)]
+#define SET_DONE(row, col)  \
+    do {                    \
+        DONE(row, col) = 1; \
+    } while (0)
+#define IS_NOTDONE(row, col) !DONE(row, col)
+#define IS_DONE(row, col)    DONE(row, col)
+static char *done;
+#endif
+
+#define E  1
+#define SE 2
+#define S  4
+#define SW 8
+#define W  16
+#define NW 32
+#define N  64
+#define NE 128
+
+static int dir_checks[3][3] = {{SE, S, SW}, {E, 0, W}, {NE, N, NW}};
+
+struct cell {
+    int row;
+    int col;
+};
+
+struct cell_stack {
+    struct cell *cells;
+    int n;
+    int nalloc;
+};
+
+static int nrows, ncols;
+
+static void trace_up(struct raster_map *, int, int, int, struct cell_stack *);
+static void init_up_stack(struct cell_stack *);
+static void free_up_stack(struct cell_stack *);
+static void push_up(struct cell_stack *, struct cell *);
+static struct cell pop_up(struct cell_stack *);
+
+void DELINEATE(struct raster_map *dir_map, struct outlet_list *outlet_l)
+{
+    int i, j;
+
+    nrows = dir_map->nrows;
+    ncols = dir_map->ncols;
+
+#ifdef USE_LESS_MEMORY
+#pragma omp parallel for schedule(dynamic) private(j)
+    for (i = 0; i < nrows; i++) {
+        for (j = 0; j < ncols; j++)
+            SET_NOTDONE(i, j);
+    }
+#else
+    done = calloc((size_t)nrows * ncols, sizeof *done);
+#endif
+
+#pragma omp parallel for schedule(dynamic)
+    for (i = 0; i < outlet_l->n; i++) {
+#ifndef USE_LESS_MEMORY
+        SET_DONE(outlet_l->row[i], outlet_l->col[i]);
+#endif
+        SHED(outlet_l->row[i], outlet_l->col[i]) = outlet_l->id[i];
+    }
+
+    /* loop through all outlets and delineate the subwatershed for each */
+#pragma omp parallel for schedule(dynamic)
+    for (i = 0; i < outlet_l->n; i++) {
+        struct cell_stack up_stack;
+
+        init_up_stack(&up_stack);
+
+        /* trace up flow directions */
+        trace_up(dir_map, outlet_l->row[i], outlet_l->col[i], outlet_l->id[i],
+                 &up_stack);
+
+        free_up_stack(&up_stack);
+    }
+
+#pragma omp parallel for schedule(dynamic) private(j)
+    for (i = 0; i < nrows; i++) {
+        for (j = 0; j < ncols; j++)
+            if (IS_NOTDONE(i, j))
+                Rast_set_c_null_value(&SHED(i, j), 1);
+    }
+
+#ifndef USE_LESS_MEMORY
+    free(done);
+#endif
+}
+
+static void trace_up(struct raster_map *dir_map, int row, int col, int id,
+                     struct cell_stack *up_stack)
+{
+    int i, j;
+    int nup = 0;
+    int next_row = -1, next_col = -1;
+
+    for (i = -1; i <= 1; i++) {
+        /* skip edge cells */
+        if (row + i < 0 || row + i >= nrows)
+            continue;
+
+        for (j = -1; j <= 1; j++) {
+            /* skip the current and edge cells */
+            if ((i == 0 && j == 0) || col + j < 0 || col + j >= ncols)
+                continue;
+
+            /* if a neighbor cell flows into the current cell, trace up
+             * further; we need to check if that neighbor cell has already been
+             * processed because we don't want to misinterpret a subwatershed
+             * ID as a direction; remember we're overwriting dir_map so it can
+             * have both directions and subwatershed IDs */
+            if (GET_DIR(row + i, col + j) == dir_checks[i + 1][j + 1] &&
+                IS_NOTDONE(row + i, col + j)) {
+                if (++nup == 1) {
+                    /* climb up only to this cell at this time */
+                    next_row = row + i;
+                    next_col = col + j;
+#ifndef USE_LESS_MEMORY
+                    SET_DONE(next_row, next_col);
+#endif
+                    SHED(next_row, next_col) = id;
+                }
+                else
+                    /* if we found more than one upstream cell, we don't need
+                     * to find more at this point */
+                    break;
+            }
+        }
+        /* if we found more than one upstream cell, we don't need to find more
+         * at this point */
+        if (nup > 1)
+            break;
+    }
+
+    if (!nup) {
+        /* reached a ridge cell; if there were any up cells to visit, let's go
+         * back or simply complete tracing */
+        struct cell up;
+
+        if (!up_stack->n)
+            return;
+
+        up = pop_up(up_stack);
+        next_row = up.row;
+        next_col = up.col;
+    }
+    else if (nup > 1) {
+        /* if there are more up cells to visit, let's come back later */
+        struct cell up;
+
+        up.row = row;
+        up.col = col;
+        push_up(up_stack, &up);
+    }
+
+    /* use gcc -O2 or -O3 flags for tail-call optimization
+     * (-foptimize-sibling-calls) */
+    trace_up(dir_map, next_row, next_col, id, up_stack);
+}
+
+static void init_up_stack(struct cell_stack *up_stack)
+{
+    up_stack->nalloc = up_stack->n = 0;
+    up_stack->cells = NULL;
+}
+
+static void free_up_stack(struct cell_stack *up_stack)
+{
+    if (up_stack->cells)
+        free(up_stack->cells);
+    init_up_stack(up_stack);
+}
+
+static void push_up(struct cell_stack *up_stack, struct cell *up)
+{
+    if (up_stack->n == up_stack->nalloc) {
+        up_stack->nalloc += REALLOC_INCREMENT;
+        up_stack->cells = realloc(up_stack->cells,
+                                  sizeof *up_stack->cells * up_stack->nalloc);
+    }
+    up_stack->cells[up_stack->n++] = *up;
+}
+
+static struct cell pop_up(struct cell_stack *up_stack)
+{
+    struct cell up;
+
+    if (up_stack->n > 0) {
+        up = up_stack->cells[--up_stack->n];
+        if (up_stack->n == 0)
+            free_up_stack(up_stack);
+        else if (up_stack->n == up_stack->nalloc - REALLOC_INCREMENT) {
+            up_stack->nalloc -= REALLOC_INCREMENT;
+            up_stack->cells = realloc(up_stack->cells, sizeof *up_stack->cells *
+                                                           up_stack->nalloc);
+        }
+    }
+
+    return up;
+}

src/raster/r.watersheds/delineate_lessmem.c

@@ -0,0 +1,2 @@
+#define USE_LESS_MEMORY
+#include "delineate_funcs.h"

src/raster/r.watersheds/delineate_moremem.c

@@ -0,0 +1 @@
+#include "delineate_funcs.h"

src/raster/r.watersheds/gettimeofday.c

@@ -0,0 +1,28 @@
+#ifdef _MSC_VER
+#include <stdint.h>
+#include <winsock2.h>
+
+// https://stackoverflow.com/a/26085827/16079666
+int gettimeofday(struct timeval *tp, struct timezone *tzp)
+{
+    // Note: some broken versions only have 8 trailing zero's, the correct
+    // epoch has 9 trailing zero's This magic number is the number of 100
+    // nanosecond intervals since January 1, 1601 (UTC) until 00:00:00 January
+    // 1, 1970
+    static const uint64_t EPOCH = ((uint64_t)116444736000000000ULL);
+
+    SYSTEMTIME system_time;
+    FILETIME file_time;
+    uint64_t time;
+
+    GetSystemTime(&system_time);
+    SystemTimeToFileTime(&system_time, &file_time);
+    time = ((uint64_t)file_time.dwLowDateTime);
+    time += ((uint64_t)file_time.dwHighDateTime) << 32;
+
+    tp->tv_sec = (long)((time - EPOCH) / 10000000L);
+    tp->tv_usec = (long)(system_time.wMilliseconds * 1000);
+
+    return 0;
+}
+#endif

src/raster/r.watersheds/global.h

@@ -0,0 +1,61 @@
+#ifndef _GLOBAL_H_
+#define _GLOBAL_H_
+
+#ifdef _MSC_VER
+#include <winsock2.h>
+/* gettimeofday.c */
+int gettimeofday(struct timeval *, struct timezone *);
+#else
+#include <sys/time.h>
+#endif
+#include <grass/raster.h>
+
+#define REALLOC_INCREMENT 1024
+
+#define INDEX(row, col)   ((size_t)(row) * ncols + (col))
+#define DIR(row, col)     dir_map->cells.c[INDEX(row, col)]
+
+struct raster_map {
+    RASTER_MAP_TYPE type;
+    size_t cell_size;
+    int nrows, ncols;
+    union {
+        void *v;
+        unsigned char *uint8;
+        CELL *c;
+        FCELL *f;
+        DCELL *d;
+    } cells;
+};
+
+struct outlet_list {
+    int nalloc, n;
+    int *row, *col;
+    int *id;
+};
+
+/* timeval_diff.c */
+long long timeval_diff(struct timeval *, struct timeval *, struct timeval *);
+
+/* raster_map.c */
+struct raster_map *read_direction(char *, char *);
+void write_watersheds(char *, struct raster_map *);
+void free_raster_map(struct raster_map *);
+
+/* outlet_list.c */
+void init_outlet_list(struct outlet_list *);
+void reset_outlet_list(struct outlet_list *);
+void free_outlet_list(struct outlet_list *);
+void add_outlet(struct outlet_list *, double, double, int);
+struct outlet_list *read_outlets(char *, char *, char *);
+
+/* delineate.c */
+void delineate(struct raster_map *, struct outlet_list *, int);
+
+/* delineate_lessmem.c */
+void delineate_lessmem(struct raster_map *, struct outlet_list *);
+
+/* delineate_moremem.c */
+void delineate_moremem(struct raster_map *, struct outlet_list *);
+
+#endif