diff --git a/src/raster/r.flowaccumulation/Makefile b/src/raster/r.flowaccumulation/Makefile new file mode 100644 index 0000000000..fe3c796666 --- /dev/null +++ b/src/raster/r.flowaccumulation/Makefile @@ -0,0 +1,13 @@ +MODULE_TOPDIR = ../.. + +PGM = r.flowaccumulation + +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 diff --git a/src/raster/r.flowaccumulation/accumulate.c b/src/raster/r.flowaccumulation/accumulate.c new file mode 100644 index 0000000000..fc0c9ca114 --- /dev/null +++ b/src/raster/r.flowaccumulation/accumulate.c @@ -0,0 +1,164 @@ +#include +#include "global.h" + +#define DIR_NULL 0x80000000 +#define ACCUM(row, col) accum_map->cells[(size_t)(row)*ncols + (col)] +#define FIND_UP(row, col) \ + ((row > 0 ? (col > 0 && DIR(row - 1, col - 1) == SE ? NW : 0) | \ + (DIR(row - 1, col) == S ? N : 0) | \ + (col < ncols - 1 && DIR(row - 1, col + 1) == SW ? NE : 0) \ + : 0) | \ + (col > 0 && DIR(row, col - 1) == E ? W : 0) | \ + (col < ncols - 1 && DIR(row, col + 1) == W ? E : 0) | \ + (row < nrows - 1 \ + ? (col > 0 && DIR(row + 1, col - 1) == NE ? SW : 0) | \ + (DIR(row + 1, col) == N ? S : 0) | \ + (col < ncols - 1 && DIR(row + 1, col + 1) == NW ? SE : 0) \ + : 0)) + +#ifdef USE_LESS_MEMORY +#define UP(row, col) FIND_UP(row, col) +#else +#define UP(row, col) up_cells[(size_t)(row)*ncols + (col)] +static unsigned char *up_cells; +#endif + +static int nrows, ncols; + +static void trace_down(struct raster_map *, struct raster_map *, int, int, int); +static int sum_up(struct raster_map *, int, int, int); + +void accumulate(struct raster_map *dir_map, struct raster_map *accum_map) +{ + int row, col; + + nrows = dir_map->nrows; + ncols = dir_map->ncols; + +#ifndef USE_LESS_MEMORY + up_cells = calloc((size_t)nrows * ncols, sizeof *up_cells); + +#pragma omp parallel for schedule(dynamic) private(col) + for (row = 0; row < nrows; row++) { + for (col = 0; col < ncols; col++) + if (DIR(row, col) != DIR_NULL) + UP(row, col) = FIND_UP(row, col); + } +#endif + +#pragma omp parallel for schedule(dynamic) private(col) + for (row = 0; row < nrows; row++) { + for (col = 0; col < ncols; col++) + /* if the current cell is not null and has no upstream cells, start + * tracing down */ + if (DIR(row, col) != DIR_NULL && !UP(row, col)) + trace_down(dir_map, accum_map, row, col, 1); + } + +#ifndef USE_LESS_MEMORY + free(up_cells); +#endif +} + +static void trace_down(struct raster_map *dir_map, struct raster_map *accum_map, + int row, int col, int accum) +{ + int up, accum_up = 0; + + /* accumulate the current cell itself */ + ACCUM(row, col) = accum; + + /* find the downstream cell */ + switch (DIR(row, col)) { + case NW: + row--; + col--; + break; + case N: + row--; + break; + case NE: + row--; + col++; + break; + case W: + col--; + break; + case E: + col++; + break; + case SW: + row++; + col--; + break; + case S: + row++; + break; + case SE: + row++; + col++; + break; + } + + /* if the downstream cell is null or any upstream cells of the downstream + * cell have never been visited, stop tracing down */ + if (row < 0 || row >= nrows || col < 0 || col >= ncols || + DIR(row, col) == DIR_NULL || !(up = UP(row, col)) || + !(accum_up = sum_up(accum_map, row, col, up))) + return; + + /* use gcc -O2 or -O3 flags for tail-call optimization + * (-foptimize-sibling-calls) */ + trace_down(dir_map, accum_map, row, col, accum_up + 1); +} + +/* if any upstream cells have never been visited, 0 is returned; otherwise, the + * sum of upstream accumulation is returned */ +static int sum_up(struct raster_map *accum_map, int row, int col, int up) +{ + int sum = 0, accum; + +#pragma omp flush(accum_map) + if (up & NW) { + if (!(accum = ACCUM(row - 1, col - 1))) + return 0; + sum += accum; + } + if (up & N) { + if (!(accum = ACCUM(row - 1, col))) + return 0; + sum += accum; + } + if (up & NE) { + if (!(accum = ACCUM(row - 1, col + 1))) + return 0; + sum += accum; + } + if (up & W) { + if (!(accum = ACCUM(row, col - 1))) + return 0; + sum += accum; + } + if (up & E) { + if (!(accum = ACCUM(row, col + 1))) + return 0; + sum += accum; + } + if (up & SW) { + if (!(accum = ACCUM(row + 1, col - 1))) + return 0; + sum += accum; + } + if (up & S) { + if (!(accum = ACCUM(row + 1, col))) + return 0; + sum += accum; + } + if (up & SE) { + if (!(accum = ACCUM(row + 1, col + 1))) + return 0; + sum += accum; + } + + return sum; +} diff --git a/src/raster/r.flowaccumulation/gettimeofday.c b/src/raster/r.flowaccumulation/gettimeofday.c new file mode 100644 index 0000000000..813517c011 --- /dev/null +++ b/src/raster/r.flowaccumulation/gettimeofday.c @@ -0,0 +1,28 @@ +#ifdef _MSC_VER +#include +#include + +// 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 diff --git a/src/raster/r.flowaccumulation/global.h b/src/raster/r.flowaccumulation/global.h new file mode 100644 index 0000000000..7c4539664e --- /dev/null +++ b/src/raster/r.flowaccumulation/global.h @@ -0,0 +1,35 @@ +#ifndef _GLOBAL_H_ +#define _GLOBAL_H_ + +#ifdef _MSC_VER +#include +/* gettimeofday.c */ +int gettimeofday(struct timeval *, struct timezone *); +#else +#include +#endif +#include + +#define E 1 +#define SE 2 +#define S 4 +#define SW 8 +#define W 16 +#define NW 32 +#define N 64 +#define NE 128 + +#define DIR(row, col) dir_map->cells[(size_t)(row)*ncols + (col)] + +struct raster_map { + int nrows, ncols; + CELL *cells; +}; + +/* timeval_diff.c */ +long long timeval_diff(struct timeval *, struct timeval *, struct timeval *); + +/* accumulate.c */ +void accumulate(struct raster_map *, struct raster_map *); + +#endif diff --git a/src/raster/r.flowaccumulation/main.c b/src/raster/r.flowaccumulation/main.c new file mode 100644 index 0000000000..7c7931dbb3 --- /dev/null +++ b/src/raster/r.flowaccumulation/main.c @@ -0,0 +1,253 @@ +/**************************************************************************** + * + * MODULE: r.flowaccumulation + * + * AUTHOR(S): Huidae Cho + * + * PURPOSE: Calculates flow accumulation from a flow direction raster map + * using the Memory-Efficient Flow Accumulation (MEFA) parallel + * algorithm by Cho (2023). + * + * COPYRIGHT: (C) 2023 by Huidae Cho and the GRASS Development Team + * + * This program is free software under the GNU General Public + * License (>=v2). Read the file COPYING that comes with GRASS + * for details. + * + *****************************************************************************/ + +#define _MAIN_C_ + +#include +#include +#include +#ifdef _MSC_VER +#include +#else +#include +#endif +#ifdef _OPENMP +#include +#endif +#include +#include +#include "global.h" + +#define DIR_UNKNOWN 0 +#define DIR_DEG 1 +#define DIR_DEG45 2 +#define DIR_POW2 3 + +int main(int argc, char *argv[]) +{ + struct GModule *module; + struct { + struct Option *dir; + struct Option *format; + struct Option *accum; + struct Option *nprocs; + } opt; + char *desc; + char *dir_name, *accum_name; +#ifdef _OPENMP + int nprocs; +#endif + int dir_fd, accum_fd; + unsigned char dir_format; + struct Range dir_range; + CELL dir_min, dir_max; + struct raster_map *dir_map, *accum_map; + int nrows, ncols, row, col; + struct History hist; + struct timeval first_time, start_time, end_time; + + G_gisinit(argv[0]); + + module = G_define_module(); + G_add_keyword(_("raster")); + G_add_keyword(_("hydrology")); + G_add_keyword(_("accumulation")); + module->description = + _("Calculates flow accumulation from a flow direction raster map using " + "the Memory-Efficient Flow Accumulation (MEFA) parallel algorithm by " + "Cho (2023)"); + + opt.dir = G_define_standard_option(G_OPT_R_INPUT); + opt.dir->description = _("Name of input direction raster map"); + + opt.format = G_define_option(); + opt.format->type = TYPE_STRING; + opt.format->key = "format"; + opt.format->label = _("Format of input direction raster map"); + opt.format->required = YES; + opt.format->options = "auto,degree,45degree,power2"; + opt.format->answer = "auto"; + G_asprintf(&desc, "auto;%s;degree;%s;45degree;%s;power2;%s", + _("auto-detect direction format"), _("degrees CCW from East"), + _("degrees CCW from East divided by 45 (e.g. r.watershed)"), + _("powers of 2 CW from East (e.g., r.terraflow, ArcGIS)")); + opt.format->descriptions = desc; + + opt.accum = G_define_standard_option(G_OPT_R_OUTPUT); + opt.accum->description = _("Name for output flow accumulation raster map"); + +#ifdef _OPENMP + opt.nprocs = G_define_standard_option(G_OPT_M_NPROCS); + opt.nprocs->label = opt.nprocs->description; + opt.nprocs->description = _("0 to use OMP_NUM_THREADS"); + opt.nprocs->answer = "0"; +#endif + + if (G_parser(argc, argv)) + exit(EXIT_FAILURE); + +#ifdef _OPENMP + nprocs = atoi(opt.nprocs->answer); + if (nprocs < 0) + G_fatal_error(_("<%s> must be >= 0"), opt.nprocs->key); + + if (nprocs >= 1) + omp_set_num_threads(nprocs); + nprocs = omp_get_max_threads(); + G_message(n_("Using %d thread for serial computation", + "Using %d threads for parallel computation", nprocs), + nprocs); +#endif + + dir_name = opt.dir->answer; + accum_name = opt.accum->answer; + + /* read direction raster */ + G_message(_("Reading flow direction raster <%s>..."), dir_name); + gettimeofday(&start_time, NULL); + first_time = start_time; + + dir_fd = Rast_open_old(dir_name, ""); + if (Rast_get_map_type(dir_fd) != CELL_TYPE) + G_fatal_error(_("Invalid direction map <%s>"), dir_name); + + /* adapted from r.path */ + if (Rast_read_range(dir_name, "", &dir_range) < 0) + G_fatal_error(_("Unable to read range file")); + Rast_get_range_min_max(&dir_range, &dir_min, &dir_max); + if (dir_max <= 0) + G_fatal_error(_("Invalid direction map <%s>"), dir_name); + + dir_format = DIR_UNKNOWN; + if (strcmp(opt.format->answer, "degree") == 0) { + if (dir_max > 360) + G_fatal_error(_("Directional degrees cannot be > 360")); + dir_format = DIR_DEG; + } + else if (strcmp(opt.format->answer, "45degree") == 0) { + if (dir_max > 8) + G_fatal_error(_("Directional degrees divided by 45 cannot be > 8")); + dir_format = DIR_DEG45; + } + else if (strcmp(opt.format->answer, "power2") == 0) { + if (dir_max > 128) + G_fatal_error(_("Powers of 2 cannot be > 128")); + dir_format = DIR_POW2; + } + else if (strcmp(opt.format->answer, "auto") == 0) { + if (dir_max <= 8) { + dir_format = DIR_DEG45; + G_important_message(_("Flow direction format assumed to be " + "degrees CCW from East divided by 45")); + } + else if (dir_max <= 128) { + dir_format = DIR_POW2; + G_important_message(_("Flow direction format assumed to be " + "powers of 2 CW from East")); + } + else if (dir_max <= 360) { + dir_format = DIR_DEG; + G_important_message( + _("Flow direction format assumed to be degrees CCW from East")); + } + else + G_fatal_error( + _("Unable to detect format of input direction map <%s>"), + dir_name); + } + if (dir_format == DIR_UNKNOWN) + G_fatal_error(_("Invalid direction format '%s'"), opt.format->answer); + /* end of r.path */ + + nrows = Rast_window_rows(); + ncols = Rast_window_cols(); + + dir_map = G_malloc(sizeof *dir_map); + dir_map->nrows = nrows; + dir_map->ncols = ncols; + dir_map->cells = G_malloc(sizeof(CELL) * nrows * ncols); + for (row = 0; row < nrows; row++) { + G_percent(row, nrows, 1); + Rast_get_c_row(dir_fd, dir_map->cells + ncols * row, row); + if (dir_format == DIR_DEG) { + for (col = 0; col < ncols; col++) + DIR(row, col) = pow(2, abs(DIR(row, col) / 45.)); + } + else if (dir_format == DIR_DEG45) { + for (col = 0; col < ncols; col++) + DIR(row, col) = pow(2, 8 - abs(DIR(row, col))); + } + else { + for (col = 0; col < ncols; col++) + DIR(row, col) = abs(DIR(row, col)); + } + } + G_percent(1, 1, 1); + Rast_close(dir_fd); + + gettimeofday(&end_time, NULL); + G_message(_("Input time for flow direction: %f seconds"), + timeval_diff(NULL, &end_time, &start_time) / 1e6); + + /* accumulate flow */ + G_message(_("Accumulating flows...")); + gettimeofday(&start_time, NULL); + + accum_map = G_malloc(sizeof *accum_map); + accum_map->nrows = nrows; + accum_map->ncols = ncols; + accum_map->cells = G_malloc(sizeof(CELL) * nrows * ncols); + + accumulate(dir_map, accum_map); + + G_free(dir_map->cells); + G_free(dir_map); + + gettimeofday(&end_time, NULL); + G_message(_("Compute time for flow accumulation: %f seconds"), + timeval_diff(NULL, &end_time, &start_time) / 1e6); + + /* write out buffer to accumulation raster */ + G_message(_("Writing flow accumulation raster <%s>..."), accum_name); + gettimeofday(&start_time, NULL); + + accum_fd = Rast_open_new(accum_name, CELL_TYPE); + for (row = 0; row < nrows; row++) { + G_percent(row, nrows, 1); + Rast_put_c_row(accum_fd, accum_map->cells + ncols * row); + } + G_percent(1, 1, 1); + Rast_close(accum_fd); + G_free(accum_map->cells); + G_free(accum_map); + + /* write history */ + Rast_put_cell_title(accum_name, _("Flow accumulation")); + Rast_short_history(accum_name, "raster", &hist); + Rast_command_history(&hist); + Rast_write_history(accum_name, &hist); + + gettimeofday(&end_time, NULL); + G_message(_("Output time for flow accumulation: %f seconds"), + timeval_diff(NULL, &end_time, &start_time) / 1e6); + + G_message(_("Total elapsed time: %f seconds"), + timeval_diff(NULL, &end_time, &first_time) / 1e6); + + exit(EXIT_SUCCESS); +} diff --git a/src/raster/r.flowaccumulation/r.flowaccumulation.html b/src/raster/r.flowaccumulation/r.flowaccumulation.html new file mode 100644 index 0000000000..ff9fdc688d --- /dev/null +++ b/src/raster/r.flowaccumulation/r.flowaccumulation.html @@ -0,0 +1,87 @@ +

DESCRIPTION

+ +r.flowaccumulation calculates flow accumulation from a flow direction +raster map using the Memory-Efficient Flow Accumulation (MEFA) parallel +algorithm by Cho (2023). + +

NOTES

+ +Unlike r.watershed, but just like r.accumulate, +r.flowaccumulation does not require elevation data to calculate flow +accumulation. Instead, this module only uses a flow direction raster map to +trace and accumulate the amount of flow draining through and including each +cell. + +r.flowaccumulation supports parallel computation of flow accumulation +using OpenMP while r.accumulate does not. + +

The module recognizes three different formats of flow directions: +

+degree +
+ +

Since the module does not use elevation data (i.e., slope), flow +accumulation is calculated by single flow direction (SFD) routing and may not +be comparable to the result from multiple flow direction (MFD) routing. + +

EXAMPLES

+ +These examples use the North Carolina sample dataset. + +

Calculate flow accumulation using r.watershed and +r.flowaccumulation: +

+# set computational region
+g.region -p raster=elevation
+
+# calculate positive flow accumulation and drainage directions using r.watershed
+# for comparison, use -s (SFD)
+r.watershed -sa elevation=elevation accumulation=flow_accum drainage=drain_directions
+
+# calculate flow accumulation using r.flowaccumulation
+r.flowaccumulation input=drain_directions output=flow_accum_new
+
+# copy color table
+r.colors map=flow_accum_new raster=flow_accum
+
+# check difference between flow_accum and flow_accum_new
+r.mapcalc expression="flowaccumulation_diff=if(flow_accum-flow_accum_new, flow_accum-flow_accum_new, null())"
+
+ + + +

For some reason, there are slight differences between the two output maps. +The yellow and purple cells show the difference raster map +(flowaccumulation_diff). The red arrows and numbers represent drainage +directions (drain_directions) and flow accumulation by +r.watershed (flow_accum), respectively. Note that some cells +close to headwater cells are assigned 1 even though they are located downstream +of other cells. + +

+ +

For comparison, these numbers show the new flow accumulation by +r.flowaccumulation (flow_accum_new). The same cells are +properly accumulated from the headwater cells. + +

+ +

SEE ALSO

+ + +r.accumulate, +r.watershed, +r.stream.extract, +r.stream.distance + + +

REFERENCES

+ +Huidae Cho, July 2023. Memory-Efficient Flow Accumulation Using a +Look-Around Approach and Its OpenMP Parallelization. Environmental +Modelling & Software 167, 105771. +doi:10.1016/j.envsoft.2023.105771. + +

AUTHOR

+ +Huidae Cho, New Mexico State University diff --git a/src/raster/r.flowaccumulation/r_flowaccumulation_formats.png b/src/raster/r.flowaccumulation/r_flowaccumulation_formats.png new file mode 100644 index 0000000000..50c54bc719 Binary files /dev/null and b/src/raster/r.flowaccumulation/r_flowaccumulation_formats.png differ diff --git a/src/raster/r.flowaccumulation/r_flowaccumulation_nc_comparison.png b/src/raster/r.flowaccumulation/r_flowaccumulation_nc_comparison.png new file mode 100644 index 0000000000..fc19710670 Binary files /dev/null and b/src/raster/r.flowaccumulation/r_flowaccumulation_nc_comparison.png differ diff --git a/src/raster/r.flowaccumulation/r_flowaccumulation_nc_example.png b/src/raster/r.flowaccumulation/r_flowaccumulation_nc_example.png new file mode 100644 index 0000000000..42725b3f7c Binary files /dev/null and b/src/raster/r.flowaccumulation/r_flowaccumulation_nc_example.png differ diff --git a/src/raster/r.flowaccumulation/r_flowaccumulation_r_watershed_nc_example.png b/src/raster/r.flowaccumulation/r_flowaccumulation_r_watershed_nc_example.png new file mode 100644 index 0000000000..0b01f07c86 Binary files /dev/null and b/src/raster/r.flowaccumulation/r_flowaccumulation_r_watershed_nc_example.png differ diff --git a/src/raster/r.flowaccumulation/timeval_diff.c b/src/raster/r.flowaccumulation/timeval_diff.c new file mode 100644 index 0000000000..bfdb2e26c2 --- /dev/null +++ b/src/raster/r.flowaccumulation/timeval_diff.c @@ -0,0 +1,26 @@ +#include +#ifdef _MSC_VER +#include +#else +#include +#endif + +/* https://www.linuxquestions.org/questions/programming-9/how-to-calculate-time-difference-in-milliseconds-in-c-c-711096/#post3475339 + */ +long long timeval_diff(struct timeval *difference, struct timeval *end_time, + struct timeval *start_time) +{ + struct timeval temp_diff; + + if (difference == NULL) + difference = &temp_diff; + difference->tv_sec = end_time->tv_sec - start_time->tv_sec; + difference->tv_usec = end_time->tv_usec - start_time->tv_usec; + + /* Using while instead of if below makes the code slightly more robust. */ + while (difference->tv_usec < 0) { + difference->tv_usec += 1000000; + difference->tv_sec -= 1; + } + return 1000000LL * difference->tv_sec + difference->tv_usec; +}