Index: main.c =================================================================== --- main.c (revision 34032) +++ main.c (working copy) @@ -31,6 +31,12 @@ * Roberto Flor ITC-irst * *************************************************************************/ +/************************************************************************ +* New code added by Colin Nielsen * +* to use movement direction surface from r.walk 6/2008. * +* * +*************************************************************************/ + #include #include #include @@ -68,23 +74,25 @@ int main(int argc, char **argv) { - int fe, fd; + int fe, fd, dir_fd; int i, have_points = 0; int new_id; int nrows, ncols, points_row[MAX_POINTS], points_col[MAX_POINTS], npoints; int cell_open(), cell_open_new(); - int map_id; - char map_name[GNAME_MAX], *map_mapset, new_map_name[GNAME_MAX]; - char *tempfile1, *tempfile2; + int map_id, dir_id; + char map_name[GNAME_MAX], *map_mapset, new_map_name[GNAME_MAX], + dir_name[GNAME_MAX], *dir_mapset; + char *tempfile1, *tempfile2, *tempfile3; char *search_mapset; struct History history; struct Cell_head window; - struct Option *opt1, *opt2, *coordopt, *vpointopt; - struct Flag *flag1, *flag2, *flag3; + struct Option *opt1, *opt2, *coordopt, *vpointopt, *opt3; + struct Flag *flag1, *flag2, *flag3, *flag4; struct GModule *module; - int in_type; + int in_type, dir_data_type; void *in_buf; + void *dir_buf; CELL *out_buf; struct band3 bnd, bndC; struct metrics *m = NULL; @@ -92,11 +100,12 @@ struct point *list; struct point *thispoint; int ival, bsz, start_row, start_col, mode; + int costmode = 0; double east, north, val; struct point *drain(int, struct point *, int, int); + struct point *drain_cost(int, struct point *, int, int); int bsort(int, struct point *); - G_gisinit(argv[0]); module = G_define_module(); @@ -111,6 +120,14 @@ opt2 = G_define_standard_option(G_OPT_R_OUTPUT); opt2->description = _("Output drain raster map"); + opt3 = G_define_option(); + opt3->key = "indir"; + opt3->type = TYPE_STRING; + opt3->gisprompt = "old,cell,raster"; + opt3->description = + _("Name of movement direction map associated with the cost surface"); + opt3->required = NO; + coordopt = G_define_option(); coordopt->key = "coordinate"; coordopt->type = TYPE_STRING; @@ -140,6 +157,12 @@ flag3->key = 'n'; flag3->description = _("Count cell numbers along the path"); + flag4 = G_define_flag(); + flag4->key = 'd'; + flag4->description = + _ + ("The input surface is a cost surface (if checked, a direction surface must also be specified"); + if (G_parser(argc, argv)) exit(EXIT_FAILURE); @@ -147,6 +170,37 @@ strcpy(map_name, opt1->answer); strcpy(new_map_name, opt2->answer); + if (flag4->answer) { + costmode = 1; + G_message(_ + ("Directional drain selected... checking for direction raster")); + } + else { + G_message(_("Surface/Hydrology drain selected")); + } + + if (costmode == 1) { + if (!opt3->answer) { + G_fatal_error(_ + ("Direction raster not specified, if direction flag is on, a direction raster must be given")); + } + strcpy(dir_name, opt3->answer); + dir_mapset = G_find_cell2(dir_name, ""); + if (dir_mapset == NULL) + G_fatal_error(_("Raster map <%s> not found"), dir_name); + else { + G_message(_("Direction raster found <%s>"), dir_name); + } + dir_data_type = G_raster_map_type(dir_name, dir_mapset); + } + if (costmode == 0) { + if (opt3->answer) { + G_fatal_error(_ + ("Direction map <%s> should not be specified for Surface/Hydrology drains"), + opt3->answer); + } + } + /* get the name of the elevation map layer for filling */ map_mapset = G_find_cell(map_name, ""); if (!map_mapset) @@ -310,12 +364,27 @@ } G_close_cell(map_id); - /* fill one-cell pits and take a first stab at flow directions */ - filldir(fe, fd, nrows, &bnd, m); + if (costmode == 1) { + dir_buf = G_allocate_d_raster_buf(); + dir_id = G_open_cell_old(dir_name, dir_mapset); + tempfile3 = G_tempfile(); + dir_fd = open(tempfile3, O_RDWR | O_CREAT); - /* determine flow directions for more ambiguous cases */ - resolve(fd, nrows, &bndC); + for (i = 0; i < nrows; i++) { + G_get_d_raster_row(dir_id, dir_buf, i); + write(dir_fd, dir_buf, ncols * sizeof(DCELL)); + } + G_close_cell(dir_id); + } + /* only necessary for non-dir drain */ + if (costmode == 0) { + /* fill one-cell pits and take a first stab at flow directions */ + filldir(fe, fd, nrows, &bnd, m); + /* determine flow directions for more ambiguous cases */ + resolve(fd, nrows, &bndC); + } + /* free the buffers already used */ G_free(bndC.b[0]); G_free(bndC.b[1]); @@ -334,7 +403,13 @@ thispoint->row = points_row[i]; thispoint->col = points_col[i]; thispoint->next = NULL; - thispoint = drain(fd, thispoint, nrows, ncols); + /* drain algorithm selection (dir or non-dir) */ + if (costmode == 0) { + thispoint = drain(fd, thispoint, nrows, ncols); + } + if (costmode == 1) { + thispoint = drain_cost(dir_fd, thispoint, nrows, ncols); + } } /* do the output */ @@ -438,7 +513,10 @@ /* close files and free buffers */ G_close_cell(new_id); - G_put_cell_title(new_map_name, "Surface flow trace"); + if (costmode == 0) + G_put_cell_title(new_map_name, "Surface flow trace"); + if (costmode == 1) + G_put_cell_title(new_map_name, "Path of cost accumulation "); G_short_history(new_map_name, "raster", &history); G_command_history(&history); @@ -451,7 +529,12 @@ unlink(tempfile2); G_free(in_buf); G_free(out_buf); - + if (costmode == 1) { + close(dir_fd); + unlink(tempfile3); + G_free(dir_buf); + } + G_done_msg("Path(s) found"); exit(EXIT_SUCCESS); } @@ -513,3 +596,131 @@ return list; } + +struct point *drain_cost(int dir_fd, struct point *list, int nrow, int ncol) +{ + /* + * The idea is that each cell of the direction surface has a value representing + * the direction towards the next cell in the path. The direction is read from + * the input raster, and a simple case/switch is used to determine which cell to + * read next. This is repeated via a while loop until a null direction is found. + */ + int neighbour, row, col, next_row, next_col, go = 1; + DCELL direction; + DCELL *dir_buf; + + dir_buf = G_allocate_d_raster_buf(); + + next_row = list->row; + next_col = list->col; + + while (go) { + go = 0; + /* Directional algorithm + * 1) read cell direction + * 2) shift to cell in that direction + */ + /* find the direction recorded at row,col */ + lseek(dir_fd, (off_t) list->row * ncol * sizeof(DCELL), SEEK_SET); + read(dir_fd, dir_buf, ncol * sizeof(DCELL)); + direction = *(dir_buf + list->col); + neighbour = direction * 10; + if (G_verbose() > 2) + G_message(_("direction read: %lf, neighbour found: %i"), + direction, neighbour); + switch (neighbour) { + case 1800: + next_row = list->row; + next_col = list->col + 1; + break; + case 0: + next_row = list->row; + next_col = list->col - 1; + break; + case 900: + next_row = list->row + 1; + next_col = list->col; + break; + case 2700: + next_row = list->row - 1; + next_col = list->col; + break; + case 1350: + next_col = list->col + 1; + next_row = list->row + 1; + break; + case 450: + next_col = list->col - 1; + next_row = list->row + 1; + break; + case 3150: + next_row = list->row - 1; + next_col = list->col - 1; + break; + case 2250: + next_row = list->row - 1; + next_col = list->col + 1; + break; + case 1125: + next_row = list->row + 2; + next_col = list->col + 1; + break; + case 675: + next_row = list->row + 2; + next_col = list->col - 1; + break; + case 2925: + next_row = list->row - 2; + next_col = list->col - 1; + break; + case 2475: + next_row = list->row - 2; + next_col = list->col + 1; + break; + case 1575: + next_row = list->row + 1; + next_col = list->col + 2; + break; + case 225: + next_row = list->row + 1; + next_col = list->col - 2; + break; + case 3375: + next_row = list->row - 1; + next_col = list->col - 2; + break; + case 2025: + next_row = list->row - 1; + next_col = list->col + 2; + break; + default: + G_fatal_error(_("Raster map <%s> not a direction map"), dir_name); + } /* end switch/case */ + + if (next_col >= 0 && next_col < ncol && next_row >= 0 && + next_row < nrow) { + list->next = (struct point *)G_malloc(sizeof(struct point)); + list = list->next; + list->row = next_row; + list->col = next_col; + next_row = -1; + next_col = -1; + go = 1; + } + + } /* end while */ + + /* allocate and fill the end-of-path flag */ + list->next = (struct point *)G_malloc(sizeof(struct point)); + list = list->next; + list->row = INT_MAX; + + /* return a pointer to an empty structure */ + list->next = (struct point *)G_malloc(sizeof(struct point)); + list = list->next; + list->next = NULL; + + G_free(dir_buf); + + return list; +} Index: description.html =================================================================== --- description.html (revision 34032) +++ description.html (working copy) @@ -2,15 +2,17 @@ r.drain traces a flow through a least-cost path in an elevation -model. The input elevation surface (a raster map layer) might -be a cumulative cost map generated by the -r.cost program. +model. If the input surface (a raster map layer) is a cumulative +cost map generated by the r.walk or +r.cost programs, the -d flag and a +movement direction surface "indir" must be specified. The output result (also a raster map layer) will show one or more least-cost paths between each user-provided location(s) and the minima -(low category values) in the input map. By default, the output -will be an integer CELL map with 1 along the least cost path, -and null cells elsewhere. +(low category values) in the input map. If the -d flag is used the +output least-cost paths will be found using the direction layer. +By default, the output will be an integer CELL map with 1 +along the least cost path, and null cells elsewhere.

With the -c (copy) flag, the input map cell values are @@ -24,9 +26,12 @@ The -c, -a, and -n flags are mutually incompatible.

-The path is calculated by choosing the steeper "slope" between adjacent -cells. The slope calculation accurately acounts for the variable scale in -lat-lon projections. +For an elevation surface, the path is calculated by choosing the steeper +"slope" between adjacent cells. The slope calculation accurately acounts +for the variable scale in lat-lon projections. For a cost surface, the path +is calculated by following the movement direction surface back to the start +point given in r.walk or r.cost (not implimented for +r.cost yet).

The coordinate parameter consists of map E and N grid coordinates of @@ -153,6 +158,18 @@ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +With the -d (direction) flag, the direction raster is used +for the input, rather than the elevation surface. The output is then created +according to one of the -can flags. +

+The directions are recorded as GRASS standard directions:
+       112.5 90  67.5         i.e. a cell with the value 135 
+157.5  135   0   45   22.5    means the cell before it is 
+       180   x   0            to the south-east.
+202.5  225  270  315  337.5
+       247.5     292.5
+
+

BUGS

@@ -191,4 +208,7 @@ Matteo Franchi (Liceo Leonardo Da Vinci, Trento) and Roberto Flor (ITC-irst, Trento, Italy) -

Last changed: $Date$ +

+Updated to output a movement direction layer +

+Colin Nielsen, 2008