Index: main.c =================================================================== --- main.c (revision 36142) +++ 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 @@ -48,6 +54,7 @@ #include #include #include +#include #define DEBUG #include "tinf.h" @@ -68,22 +75,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, *opt4; + 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; @@ -91,10 +101,18 @@ 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 *); + struct line_pnts *Points; + struct line_cats *Cats; + struct Map_info vout; + int cat; + double x, y; + char vect[GNAME_MAX], *vect_mapset; G_gisinit(argv[0]); @@ -107,27 +125,42 @@ opt1->description = _("Name of existing raster map containing elevation surface"); + 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; + opt2 = G_define_standard_option(G_OPT_R_OUTPUT); opt2->description = _("Output drain raster map"); + opt4 = G_define_option(); + opt4->key = "voutput"; + opt4->type = TYPE_STRING; + opt4->gisprompt = "new,vector,vector"; + opt4->required = NO; + opt4->description = + _ + ("Output drain vector map (recommended for cost surface made using knight's move)"); + coordopt = G_define_option(); coordopt->key = "coordinate"; coordopt->type = TYPE_STRING; coordopt->required = NO; coordopt->multiple = YES; coordopt->key_desc = "x,y"; - coordopt->description = - _("The E and N coordinates of starting point(s)"); + coordopt->description = _("The E and N coordinates of starting point(s)"); vpointopt = G_define_option(); vpointopt->key = "vector_points"; vpointopt->type = TYPE_STRING; vpointopt->required = NO; vpointopt->multiple = YES; - vpointopt->gisprompt= "old,vector,vector"; + vpointopt->gisprompt = "old,vector,vector"; vpointopt->key_desc = "name"; - vpointopt->description = - _("Vector map(s) containing starting point(s)"); + vpointopt->description = _("Vector map(s) containing starting point(s)"); flag1 = G_define_flag(); flag1->key = 'c'; @@ -141,6 +174,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); @@ -148,10 +187,54 @@ 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); + } + } + + if (opt4->answer) { + G_message(_("Outputting a vector path")); + /*vect = opt4->answer; */ + if (G_legal_filename(opt4->answer) < 0) + G_fatal_error(_("<%s> is an illegal file name"), opt4->answer); + /*G_ask_vector_new("",vect); */ + if (0 > Vect_open_new(&vout, opt4->answer, 0)) { + G_fatal_error(_("Unable to create vector map <%s>"), + opt4->answer); + } + Vect_hist_command(&vout); + } + /* get the name of the elevation map layer for filling */ map_mapset = G_find_cell(map_name, ""); if (!map_mapset) - G_fatal_error(_("Could not access %s layer"), map_name); + G_fatal_error(_("Raster map <%s> not found"), map_name); /* allocate cell buf for the map layer */ in_type = G_raster_map_type(map_name, map_mapset); @@ -159,7 +242,7 @@ /* set the pointers for multi-typed functions */ set_func_pointers(in_type); - if( (flag1->answer + flag2->answer + flag3->answer) > 1 ) + if ((flag1->answer + flag2->answer + flag3->answer) > 1) G_fatal_error(_("Specify just one of the -c, -a and -n flags")); mode = 0; @@ -174,6 +257,10 @@ G_get_window(&window); nrows = G_window_rows(); ncols = G_window_cols(); + if (vect) { + Points = Vect_new_line_struct(); + Cats = Vect_new_cats_struct(); + } /* calculate true cell resolution */ m = (struct metrics *)G_malloc(nrows * sizeof(struct metrics)); @@ -190,20 +277,23 @@ if (start_row < 0 || start_row > nrows || start_col < 0 || start_col > ncols) { - G_warning(_("Starting point %d is outside the current region"), + G_warning(_ + ("Starting point %d is outside the current region"), i + 1); continue; } points_row[npoints] = start_row; points_col[npoints] = start_col; npoints++; - if(npoints >= MAX_POINTS) G_fatal_error(_("Too many start points.")); + if (npoints >= MAX_POINTS) + G_fatal_error(_("Too many start points")); have_points = 1; } } if (vpointopt->answer) { for (i = 0; vpointopt->answers[i] != NULL; i++) { - FILE *fp; + struct Map_info *fp; + /* struct start_pt *new_start_pt; */ Site *site = NULL; /* pointer to Site */ int dims, strs, dbls; @@ -211,12 +301,12 @@ search_mapset = G_find_sites(vpointopt->answers[i], ""); if (search_mapset == NULL) - G_fatal_error(_("Vector map %s - not found"), + G_fatal_error(_("Vector map <%s> not found"), vpointopt->answers[i]); fp = G_fopen_sites_old(vpointopt->answers[i], search_mapset); - if (0 != G_site_describe( fp, &dims, &cat, &strs, &dbls)) + if (0 != G_site_describe(fp, &dims, &cat, &strs, &dbls)) G_fatal_error(_("Failed to guess site file format")); site = G_site_new_struct(cat, dims, strs, dbls); @@ -229,20 +319,23 @@ start_row = (int)G_northing_to_row(site->north, &window); /* effectively just a duplicate check to G_site_in_region() ??? */ - if (start_row < 0 || start_row > nrows || start_col < 0 || start_col > ncols) + if (start_row < 0 || start_row > nrows || start_col < 0 || + start_col > ncols) continue; points_row[npoints] = start_row; points_col[npoints] = start_col; npoints++; - if(npoints >= MAX_POINTS) G_fatal_error(_("Too many start points.")); + if (npoints >= MAX_POINTS) + G_fatal_error(_("Too many start points")); have_points = 1; } /* only catches maps out of range until something is found, not after */ - if(!have_points) { - G_warning(_("Starting vector map <%s> contains no points in the current region."), - vpointopt->answers[i]); + if (!have_points) { + G_warning(_ + ("Starting vector map <%s> contains no points in the current region"), + vpointopt->answers[i]); } } } @@ -259,6 +352,7 @@ G_begin_distance_calculations(); { double e1, n1, e2, n2; + e1 = window.east; n1 = window.north; e2 = e1 + window.ew_res; @@ -306,20 +400,35 @@ } 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]); - G_free (bndC.b[2]); + G_free(bndC.b[0]); + G_free(bndC.b[1]); + G_free(bndC.b[2]); - G_free (bnd.b[0]); - G_free (bnd.b[1]); - G_free (bnd.b[2]); + G_free(bnd.b[0]); + G_free(bnd.b[1]); + G_free(bnd.b[2]); /* determine the drainage paths */ @@ -330,7 +439,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 */ @@ -431,8 +546,50 @@ } } + /* Output a vector path */ + if (opt4->answer) { + /* Need to modify for multiple paths */ + thispoint = list; + i = 1; + while (thispoint->next != NULL) { + if (thispoint->row == INT_MAX) { + thispoint = thispoint->next; + Vect_cat_set(Cats, 1, i); + Vect_write_line(&vout, GV_LINE, Points, Cats); + Vect_reset_line(Points); + Vect_reset_cats(Cats); + i++; + continue; + } + if (Vect_cat_get(Cats, 1, &cat) == 0) { + Vect_cat_set(Cats, 1, i); + } + /* Need to convert row and col to coordinates + * y = cell_head.north - ((double) p->row + 0.5) * cell_head.ns_res; + * x = cell_head.west + ((double) p->col + 0.5) * cell_head.ew_res; + */ + + x = window.west + ((double)thispoint->col + 0.5) * window.ew_res; + y = window.north - ((double)thispoint->row + 0.5) * window.ns_res; + Vect_append_point(Points, x, y, 0.0); + thispoint = thispoint->next; + } /* End while */ + Vect_build(&vout, stderr); + Vect_close(&vout); + } + /* close files and free buffers */ G_close_cell(new_id); + + 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); + G_write_history(new_map_name, &history); + close(fe); close(fd); @@ -440,7 +597,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); } @@ -502,3 +664,134 @@ 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: + break; + Should probably add something here for the possibility of a non-direction map + G_fatal_error(_("Invalid direction given (possibly not a direction map)")); */ + } /* 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 36142) +++ 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 modules, 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 @@ -56,6 +61,7 @@ of the map. In this case, the user could try adjusting the region extents slightly with g.region to allow additional outlet paths for r.drain. +

EXAMPLES

Consider the following example: @@ -153,6 +159,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