Skip to content

Commit

Permalink
r.accumulate: No need to accumulate flow for subwatershed delineation (…
Browse files Browse the repository at this point in the history 
…#189)

* No need to accumulate flow for subwatershed delineation

* comments in 80 columns
  • Loading branch information
@HuidaeCho
HuidaeCho committed Jun 1, 2020
1 parent 00d8573 commit 372f4da
Show file tree
Hide file tree
Showing 3 changed files with 114 additions and 115 deletions.
40 changes: 14 additions & 26 deletions grass7/raster/r.accumulate/delineate_subwatersheds.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,20 +5,18 @@
static struct Cell_head window;
static int rows, cols;

static void trace_up(struct cell_map *, struct raster_map *, char **, int,
int, int);
static void trace_up(struct cell_map *, char **, int, int, int);

void delineate_subwatersheds(struct Map_info *Map, struct cell_map *dir_buf,
struct raster_map *accum_buf, char **done,
int *id, struct point_list *outlet_pl)
void delineate_subwatersheds(struct cell_map *dir_buf, char **done, int *id,
struct point_list *outlet_pl)
{
int i, j;
int subwshed_id;

G_get_set_window(&window);

rows = accum_buf->rows;
cols = accum_buf->cols;
rows = dir_buf->rows;
cols = dir_buf->cols;

G_message(_("Flagging outlets..."));
for (i = 0; i < outlet_pl->n; i++) {
Expand Down Expand Up @@ -60,7 +58,7 @@ void delineate_subwatersheds(struct Map_info *Map, struct cell_map *dir_buf,
done[row][col] = 0;

/* trace up flow directions */
trace_up(dir_buf, accum_buf, done, row, col, subwshed_id);
trace_up(dir_buf, done, row, col, subwshed_id);

/* flag the current outlet again */
done[row][col] = 1;
Expand All @@ -77,24 +75,18 @@ void delineate_subwatersheds(struct Map_info *Map, struct cell_map *dir_buf,
G_percent(1, 1, 1);
}

static void trace_up(struct cell_map *dir_buf, struct raster_map *accum_buf,
char **done, int row, int col, int id)
static void trace_up(struct cell_map *dir_buf, char **done, int row, int col,
int id)
{
double cur_acc;
int i, j;

/* if the current cell is outside the computational region or already
* processed, stop tracing */
if (row < 0 || row >= rows || col < 0 || col >= cols || done[row][col])
return;

/* if the current accumulation is 1 (headwater), stop tracing */
cur_acc = get(accum_buf, row, col);
if (cur_acc == 1) {
dir_buf->c[row][col] = id;
done[row][col] = 1;
return;
}
dir_buf->c[row][col] = id;
done[row][col] = 1;

for (i = -1; i <= 1; i++) {
/* skip edge cells */
Expand All @@ -106,16 +98,12 @@ static void trace_up(struct cell_map *dir_buf, struct raster_map *accum_buf,
if ((i == 0 && j == 0) || col + j < 0 || col + j >= cols)
continue;

/* if a neighbor cell flows into the current cell, store its
* accumulation in the accum array; upstream accumulation must
* always be less than the current accumulation; upstream
* accumulation can be greater than the current accumulation in a
* subaccumulation raster */
if (dir_buf->c[row + i][col + j] == dir_checks[i + 1][j + 1][0] &&
get(accum_buf, row + i, col + j) < cur_acc) {
/* if a neighbor cell flows into the current cell, add it to the
* map and trace up further */
if (dir_buf->c[row + i][col + j] == dir_checks[i + 1][j + 1][0]) {
dir_buf->c[row][col] = id;
done[row][col] = 1;
trace_up(dir_buf, accum_buf, done, row + i, col + j, id);
trace_up(dir_buf, done, row + i, col + j, id);
}
}
}
Expand Down
3 changes: 1 addition & 2 deletions grass7/raster/r.accumulate/global.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,8 +101,7 @@ void subaccumulate(struct Map_info *, struct cell_map *, struct raster_map *,
struct point_list *);

/* delineate_subwatersheds.c */
void delineate_subwatersheds(struct Map_info *, struct cell_map *,
struct raster_map *, char **, int *,
void delineate_subwatersheds(struct cell_map *, char **, int *,
struct point_list *);

/* delineate_streams.c */
Expand Down
186 changes: 99 additions & 87 deletions grass7/raster/r.accumulate/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ int main(int argc, char *argv[])
char neg, accum, conf;
char **done;
struct cell_map dir_buf;
struct raster_map weight_buf, accum_buf;
struct raster_map accum_buf;
int rows, cols, row, col;
struct Map_info Map;
struct point_list outlet_pl;
Expand Down Expand Up @@ -440,92 +440,102 @@ int main(int argc, char *argv[])
G_percent(1, 1, 1);
Rast_close(dir_fd);

/* prepare to create accumulation buffer */
accum_buf.rows = rows;
accum_buf.cols = cols;
accum_buf.map.v = (void **)G_malloc(rows * sizeof(void *));

/* optionally, read a weight map */
weight_buf.map.v = NULL;
if (weight_name) {
int weight_fd = Rast_open_old(weight_name, "");

accum_buf.type = weight_buf.type = Rast_get_map_type(weight_fd);
weight_buf.rows = rows;
weight_buf.cols = cols;
weight_buf.map.v = (void **)G_malloc(rows * sizeof(void *));
G_message(_("Reading weight map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
weight_buf.map.v[row] =
(void *)Rast_allocate_buf(weight_buf.type);
Rast_get_row(weight_fd, weight_buf.map.v[row], row,
weight_buf.type);
/* prepare to create accumulation buffer only when needed */
if (accum_name || subaccum_name || stream_name || lfp_name) {
struct raster_map weight_buf;

accum_buf.rows = rows;
accum_buf.cols = cols;
accum_buf.map.v = (void **)G_malloc(rows * sizeof(void *));

/* optionally, read a weight map */
weight_buf.map.v = NULL;
if (weight_name) {
int weight_fd = Rast_open_old(weight_name, "");

accum_buf.type = weight_buf.type = Rast_get_map_type(weight_fd);
weight_buf.rows = rows;
weight_buf.cols = cols;
weight_buf.map.v = (void **)G_malloc(rows * sizeof(void *));
G_message(_("Reading weight map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
weight_buf.map.v[row] =
(void *)Rast_allocate_buf(weight_buf.type);
Rast_get_row(weight_fd, weight_buf.map.v[row], row,
weight_buf.type);
}
G_percent(1, 1, 1);
Rast_close(weight_fd);
}
G_percent(1, 1, 1);
Rast_close(weight_fd);
}
/* create non-weighted accumulation if input accumulation is not given */
else if (!input_accum_name && !input_subaccum_name)
accum_buf.type = CELL_TYPE;

/* optionally, read an accumulation or subaccumulation map */
if (input_accum_name || input_subaccum_name) {
int accum_fd =
Rast_open_old(input_accum_name ? input_accum_name :
input_subaccum_name, "");

accum_buf.type = Rast_get_map_type(accum_fd);
G_message(input_accum_name ? _("Reading accumulation map...") :
_("Reading subaccumulation map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
accum_buf.map.v[row] = (void *)Rast_allocate_buf(accum_buf.type);
Rast_get_row(accum_fd, accum_buf.map.v[row], row, accum_buf.type);
/* create non-weighted accumulation if input accumulation is not given
*/
else if (!input_accum_name && !input_subaccum_name)
accum_buf.type = CELL_TYPE;

/* optionally, read an accumulation or subaccumulation map */
if (input_accum_name || input_subaccum_name) {
int accum_fd =
Rast_open_old(input_accum_name ? input_accum_name :
input_subaccum_name, "");

accum_buf.type = Rast_get_map_type(accum_fd);
G_message(input_accum_name ? _("Reading accumulation map...") :
_("Reading subaccumulation map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
accum_buf.map.v[row] =
(void *)Rast_allocate_buf(accum_buf.type);
Rast_get_row(accum_fd, accum_buf.map.v[row], row,
accum_buf.type);
}
G_percent(1, 1, 1);
Rast_close(accum_fd);
}
/* accumulate flows if input accumulation is not given */
else {
for (row = 0; row < rows; row++)
accum_buf.map.v[row] =
(void *)Rast_allocate_buf(accum_buf.type);
accumulate(&dir_buf, &weight_buf, &accum_buf, done, neg);
}
G_percent(1, 1, 1);
Rast_close(accum_fd);
}
/* accumulate flows if input accumulation is not given */
else {
for (row = 0; row < rows; row++)
accum_buf.map.v[row] = (void *)Rast_allocate_buf(accum_buf.type);
accumulate(&dir_buf, &weight_buf, &accum_buf, done, neg);
}

/* free buffer memory */
for (row = 0; row < rows; row++) {
G_free(done[row]);
/* free buffer memory */
for (row = 0; row < rows; row++) {
G_free(done[row]);
if (weight_name)
G_free(weight_buf.map.v[row]);
}
G_free(done);
if (weight_name)
G_free(weight_buf.map.v[row]);
}
G_free(done);
if (weight_name)
G_free(weight_buf.map.v);
G_free(weight_buf.map.v);

/* write out buffer to the accumulation map if requested */
if (accum_name) {
int accum_fd = Rast_open_new(accum_name, accum_buf.type);
struct History hist;
/* write out buffer to the accumulation map if requested */
if (accum_name) {
int accum_fd = Rast_open_new(accum_name, accum_buf.type);
struct History hist;

G_message(_("Writing accumulation map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
Rast_put_row(accum_fd, accum_buf.map.v[row], accum_buf.type);
}
G_percent(1, 1, 1);
Rast_close(accum_fd);
G_message(_("Writing accumulation map..."));
for (row = 0; row < rows; row++) {
G_percent(row, rows, 1);
Rast_put_row(accum_fd, accum_buf.map.v[row], accum_buf.type);
}
G_percent(1, 1, 1);
Rast_close(accum_fd);

/* write history */
Rast_put_cell_title(accum_name,
weight_name ? _("Weighted flow accumulation")
: (neg ?
_("Flow accumulation with likely underestimates")
: _("Flow accumulation")));
Rast_short_history(accum_name, "raster", &hist);
Rast_command_history(&hist);
Rast_write_history(accum_name, &hist);
/* write history */
Rast_put_cell_title(accum_name,
weight_name ? _("Weighted flow accumulation")
: (neg ?
_("Flow accumulation with likely underestimates")
: _("Flow accumulation")));
Rast_short_history(accum_name, "raster", &hist);
Rast_command_history(&hist);
Rast_write_history(accum_name, &hist);
}
}
else
accum_buf.map.v = NULL;

/* delineate stream networks */
if (stream_name) {
Expand All @@ -544,8 +554,7 @@ int main(int argc, char *argv[])

/* calculate subaccumulation if needed; this process overwrite accum_buf to
* save memory */
if (subaccum_name ||
(!input_subaccum_name && (subwshed_name || (lfp_name && !accum)))) {
if (subaccum_name || (lfp_name && !input_subaccum_name && !accum)) {
subaccumulate(&Map, &dir_buf, &accum_buf, &outlet_pl);

/* write out buffer to the subaccumulation map if requested */
Expand Down Expand Up @@ -590,6 +599,13 @@ int main(int argc, char *argv[])
Vect_close(&Map);
}

/* free buffer memory */
if (accum_buf.map.v) {
for (row = 0; row < rows; row++)
G_free(accum_buf.map.v[row]);
G_free(accum_buf.map.v);
}

/* delineate subwatersheds; this process overwrites dir_buf to save memory
*/
if (subwshed_name) {
Expand All @@ -600,8 +616,7 @@ int main(int argc, char *argv[])
for (row = 0; row < rows; row++)
done[row] = (char *)G_calloc(cols, 1);

delineate_subwatersheds(&Map, &dir_buf, &accum_buf, done, id,
&outlet_pl);
delineate_subwatersheds(&dir_buf, done, id, &outlet_pl);

subwshed_fd = Rast_open_c_new(subwshed_name);
G_message(_("Writing subwatershed map..."));
Expand All @@ -624,12 +639,9 @@ int main(int argc, char *argv[])
}

/* free buffer memory */
for (row = 0; row < rows; row++) {
for (row = 0; row < rows; row++)
G_free(dir_buf.c[row]);
G_free(accum_buf.map.v[row]);
}
G_free(dir_buf.c);
G_free(accum_buf.map.v);

exit(EXIT_SUCCESS);
}

0 comments on commit 372f4da

Please sign in to comment.