Index: main.c =================================================================== --- main.c (revision 34032) +++ main.c (working copy) @@ -25,6 +25,8 @@ * Updated for GRASS 6.1 * Roberto Flor and Markus Neteler * Glynn Clements , Soeren Gebbert + * Updated for calculation errors and directional surface generation + * Colin Nielsen * PURPOSE: anisotropic movements on cost surfaces * COPYRIGHT: (C) 1999-2006 by the GRASS Development Team * @@ -115,12 +117,12 @@ /* *************************************************************** */ int main(int argc, char *argv[]) { - void *dtm_cell, *cost_cell, *cum_cell, *cell2 = NULL; - SEGMENT dtm_in_seg, cost_in_seg, out_seg; + void *dtm_cell, *cost_cell, *cum_cell, *dir_cell, *cell2 = NULL; + SEGMENT dtm_in_seg, cost_in_seg, out_seg, out_seg2; char *dtm_mapset, *cost_mapset; - char *cum_cost_mapset; + char *cum_cost_mapset, *move_dir_mapset; char *current_mapset; - char *dtm_in_file, *cost_in_file, *out_file; + char *dtm_in_file, *cost_in_file, *out_file, *dir_out_file; char *search_mapset; double *dtm_value, *cost_value, *value_start_pt; char buf[400]; @@ -128,16 +130,18 @@ double NS_fac, EW_fac, DIAG_fac, H_DIAG_fac, V_DIAG_fac; double fcost_dtm, fcost_cost; double min_cost, old_min_cost; + double cur_dir, old_cur_dir; double zero = 0.0; int at_percent = 0; int col = 0, row = 0, nrows = 0, ncols = 0; int maxcost, par_number; int maxmem; int nseg; - int cost_fd, cum_fd, dtm_fd; + int cost_fd, cum_fd, dtm_fd, dir_fd; + int dir = 0; int have_start_points; int have_stop_points; - int dtm_in_fd, cost_in_fd, out_fd; + int dtm_in_fd, cost_in_fd, out_fd, dir_out_fd; double my_dtm, my_cost; double null_cost; double a, b, c, d, lambda, slope_factor; @@ -152,7 +156,7 @@ struct GModule *module; struct Flag *flag2, *flag3, *flag4; struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8; - struct Option *opt9, *opt10, *opt11, *opt12, *opt13, *opt14; + struct Option *opt9, *opt10, *opt11, *opt12, *opt13, *opt14, *opt15; struct cost *pres_cell, *new_cell; struct History history; struct start_pt *pres_start_pt = NULL; @@ -160,7 +164,7 @@ void *ptr2; RASTER_MAP_TYPE dtm_data_type, data_type2, cost_data_type, cum_data_type = - DCELL_TYPE, cat; + DCELL_TYPE, dir_data_type = DCELL_TYPE, cat; double peak = 0.0; int dtm_dsize, cost_dsize; @@ -202,6 +206,14 @@ opt1->gisprompt = "new,cell,raster"; opt1->description = _("Name of raster map to contain results"); + opt15 = G_define_option(); + opt15->key = "outdir"; + opt15->type = TYPE_STRING; + opt15->required = NO; + opt15->gisprompt = "new,cell,raster"; + opt15->description = + _("Name of output raster map to contain movement directions"); + opt7 = G_define_option(); opt7->key = "start_points"; opt7->type = TYPE_STRING; @@ -282,7 +294,8 @@ opt11->multiple = NO; opt11->answer = "1.0"; opt11->description = - _("Lambda coefficients for combining walking energy and friction cost"); + _ + ("Lambda coefficients for combining walking energy and friction cost"); opt13 = G_define_option(); opt13->key = "slope_factor"; @@ -310,11 +323,17 @@ if (G_parser(argc, argv)) exit(EXIT_FAILURE); + /* If no outdir is specified, set flag to skip all dir */ + if (opt15->answer != NULL) + dir = 1; + /* Initalize access to database and create temporary files */ dtm_in_file = G_tempfile(); cost_in_file = G_tempfile(); out_file = G_tempfile(); + if (dir == 1) + dir_out_file = G_tempfile(); /* Get database window parameters */ @@ -492,7 +511,8 @@ if (!G_is_d_null_value(&null_cost)) { if (null_cost < 0.0) { - G_warning(_("Warning: assigning negative cost to null cell. Null cells excluded.")); + G_warning(_ + ("Warning: assigning negative cost to null cell. Null cells excluded.")); G_set_d_null_value(&null_cost, 1); } } @@ -508,7 +528,13 @@ search_mapset = ""; cum_cost_mapset = G_find_cell2(cum_cost_layer, search_mapset); + if (dir == 1) { + strcpy(move_dir_layer, opt15->answer); + search_mapset = ""; + move_dir_mapset = G_find_cell2(move_dir_layer, search_mapset); + } + /* Check if dtm layer exists in data base */ strcpy(dtm_layer, opt2->answer); @@ -530,6 +556,11 @@ if (G_legal_filename(cum_cost_layer) < 0) G_fatal_error(_("<%s> is an illegal file name"), cum_cost_layer); + if (dir == 1) { + if (G_legal_filename(move_dir_layer) < 0) + G_fatal_error(_("<%s> is an illegal file name"), move_dir_layer); + } + /* Find number of rows and columns in window */ nrows = G_window_rows(); @@ -665,6 +696,13 @@ segment_format(out_fd, nrows, ncols, srows, scols, sizeof(double)); close(out_fd); + if (dir == 1) { + dir_out_fd = creat(dir_out_file, 0600); + segment_format(dir_out_fd, nrows, ncols, srows, scols, + sizeof(double)); + close(dir_out_fd); + } + /* Open initialize and segment all files */ dtm_in_fd = open(dtm_in_file, 2); @@ -676,6 +714,11 @@ out_fd = open(out_file, 2); segment_init(&out_seg, out_fd, segments_in_memory); + if (dir == 1) { + dir_out_fd = open(dir_out_file, 2); + segment_init(&out_seg2, dir_out_fd, segments_in_memory); + } + /* Write the cost layer in the segmented file */ @@ -817,7 +860,8 @@ fbuff = (double *)G_malloc((unsigned int)(ncols * sizeof(double))); if (fbuff == NULL) - G_fatal_error(_("Unable to allocate memory for segment fbuff == NULL")); + G_fatal_error(_ + ("Unable to allocate memory for segment fbuff == NULL")); G_set_d_null_value(fbuff, ncols); @@ -836,6 +880,35 @@ G_free(fbuff); } + if (dir == 1) { + G_message(_("Initializing directional output ")); + { + double *fbuff; + int i; + + fbuff = + (double *)G_malloc((unsigned int)(ncols * sizeof(double))); + + if (fbuff == NULL) + G_fatal_error(_ + ("Unable to allocate memory for segment fbuff == NULL")); + + G_set_d_null_value(fbuff, ncols); + + for (row = 0; row < nrows; row++) { + { + G_percent(row, nrows, 2); + } + for (i = 0; i < ncols; i++) { + segment_put(&out_seg2, &fbuff[i], row, i); + } + } + segment_flush(&out_seg2); + G_percent(row, nrows, 2); + G_free(fbuff); + } + } + /* Scan the existing cum_cost_layer searching for starting points. * Create a btree of starting points ordered by increasing costs. */ @@ -846,7 +919,8 @@ cum_cost_mapset = G_find_cell2(cum_cost_layer, search_mapset); if (cum_cost_mapset == NULL) - G_fatal_error(_("Raster output map <%s> not found (no start_points given)"), + G_fatal_error(_ + ("Raster output map <%s> not found (no start_points given)"), cum_cost_layer); cum_fd = G_open_cell_old(cum_cost_layer, cum_cost_mapset); @@ -861,7 +935,8 @@ cell2 = G_allocate_raster_buf(data_type2); if (cell2 == NULL) - G_fatal_error(_("Memory allocation error on reading start points from raster map %s"), + G_fatal_error(_ + ("Memory allocation error on reading start points from raster map %s"), cum_cost_layer); G_message(_("Reading %s... "), cum_cost_layer); @@ -911,7 +986,8 @@ value_start_pt = &zero; if (top_start_pt->row < 0 || top_start_pt->row >= nrows || top_start_pt->col < 0 || top_start_pt->col >= ncols) - G_fatal_error(_("Specified starting location outside database window")); + G_fatal_error(_ + ("Specified starting location outside database window")); new_cell = insert(zero, top_start_pt->row, top_start_pt->col); segment_put(&out_seg, value_start_pt, top_start_pt->row, top_start_pt->col); @@ -927,7 +1003,6 @@ */ - /*system("date"); */ G_message(_("Finding cost path")); n_processed = 0; @@ -980,66 +1055,82 @@ case 1: row = pres_cell->row; col = pres_cell->col - 1; + cur_dir = 180.0; break; case 2: row = pres_cell->row; col = pres_cell->col + 1; + cur_dir = 0.0; break; case 3: row = pres_cell->row - 1; col = pres_cell->col; + cur_dir = 90.0; break; case 4: row = pres_cell->row + 1; col = pres_cell->col; + cur_dir = 270.0; break; case 5: row = pres_cell->row - 1; col = pres_cell->col - 1; + cur_dir = 135.0; break; case 6: row = pres_cell->row - 1; col = pres_cell->col + 1; + cur_dir = 45.0; break; case 7: col = pres_cell->col + 1; row = pres_cell->row + 1; + cur_dir = 315.0; break; case 8: col = pres_cell->col - 1; row = pres_cell->row + 1; + cur_dir = 225.0; break; case 9: row = pres_cell->row - 2; col = pres_cell->col - 1; + cur_dir = 112.5; break; case 10: row = pres_cell->row - 2; col = pres_cell->col + 1; + cur_dir = 67.5; break; case 11: row = pres_cell->row + 2; col = pres_cell->col + 1; + cur_dir = 292.5; break; case 12: row = pres_cell->row + 2; col = pres_cell->col - 1; + cur_dir = 247.5; break; case 13: row = pres_cell->row - 1; col = pres_cell->col - 2; + cur_dir = 157.5; break; case 14: row = pres_cell->row - 1; col = pres_cell->col + 2; + cur_dir = 22.5; break; case 15: row = pres_cell->row + 1; col = pres_cell->col + 2; + cur_dir = 337.5; break; case 16: row = pres_cell->row + 1; col = pres_cell->col - 2; + cur_dir = 202.5; break; } @@ -1062,7 +1153,7 @@ fcost_dtm = (double)((double)(W_dtm - my_dtm) * d); else fcost_dtm = (double)((double)(W_dtm - my_dtm) * c); - fcost_cost = ((double)(W_cost + my_cost) / 2.0); + fcost_cost = (double)((W_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (EW_fac * a) + lambda * fcost_cost * EW_fac; @@ -1080,7 +1171,7 @@ fcost_dtm = (double)(E_dtm - my_dtm) * d; else fcost_dtm = (double)(E_dtm - my_dtm) * c; - fcost_cost = (double)(E_cost + my_cost) / 2.0; + fcost_cost = (double)((E_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (EW_fac * a) + lambda * fcost_cost * EW_fac; @@ -1098,7 +1189,7 @@ fcost_dtm = (double)(N_dtm - my_dtm) * d; else fcost_dtm = (double)(N_dtm - my_dtm) * c; - fcost_cost = (double)(N_cost + my_cost) / 2.0; + fcost_cost = (double)((N_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (NS_fac * a) + lambda * fcost_cost * NS_fac; @@ -1116,7 +1207,7 @@ fcost_dtm = (double)(S_dtm - my_dtm) * d; else fcost_dtm = (double)(S_dtm - my_dtm) * c; - fcost_cost = (double)(S_cost + my_cost) / 2.0; + fcost_cost = (double)((S_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (NS_fac * a) + lambda * fcost_cost * NS_fac; @@ -1134,7 +1225,7 @@ fcost_dtm = (double)(NW_dtm - my_dtm) * d; else fcost_dtm = (double)(NW_dtm - my_dtm) * c; - fcost_cost = (double)(NW_cost + my_cost) / 2.0; + fcost_cost = (double)((NW_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) + lambda * fcost_cost * DIAG_fac; @@ -1152,7 +1243,7 @@ fcost_dtm = (double)(NE_dtm - my_dtm) * d; else fcost_dtm = (double)(NE_dtm - my_dtm) * c; - fcost_cost = (double)(NE_cost + my_cost) / 2.0; + fcost_cost = (double)((NE_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) + lambda * fcost_cost * DIAG_fac; @@ -1170,7 +1261,7 @@ fcost_dtm = (double)(SE_dtm - my_dtm) * d; else fcost_dtm = (double)(SE_dtm - my_dtm) * c; - fcost_cost = (double)(SE_cost + my_cost) / 2.0; + fcost_cost = (double)((SE_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) + lambda * fcost_cost * DIAG_fac; @@ -1188,7 +1279,7 @@ fcost_dtm = (double)(SW_dtm - my_dtm) * d; else fcost_dtm = (double)(SW_dtm - my_dtm) * c; - fcost_cost = (double)(SW_cost + my_cost) / 2.0; + fcost_cost = (double)((SW_cost / 2.0) + (my_cost / 2.0)); min_cost = pres_cell->min_cost + fcost_dtm + (DIAG_fac * a) + lambda * fcost_cost * DIAG_fac; @@ -1207,7 +1298,8 @@ else fcost_dtm = (double)(NNW_dtm - my_dtm) * c; fcost_cost = - (double)(N_cost + NW_cost + NNW_cost + my_cost) / 4.0; + (double)((N_cost / 4.0) + (NW_cost / 4.0) + + (NNW_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) + lambda * fcost_cost * V_DIAG_fac; @@ -1226,7 +1318,8 @@ else fcost_dtm = (double)(NNE_dtm - my_dtm) * c; fcost_cost = - (double)(N_cost + NE_cost + NNE_cost + my_cost) / 4.0; + (double)((N_cost / 4.0) + (NE_cost / 4.0) + + (NNE_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) + lambda * fcost_cost * V_DIAG_fac; @@ -1245,7 +1338,8 @@ else fcost_dtm = (double)(SSE_dtm - my_dtm) * c; fcost_cost = - (double)(S_cost + SE_cost + SSE_cost + my_cost) / 4.0; + (double)((S_cost / 4.0) + (SE_cost / 4.0) + + (SSE_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) + lambda * fcost_cost * V_DIAG_fac; @@ -1264,7 +1358,8 @@ else fcost_dtm = (double)(SSW_dtm - my_dtm) * c; fcost_cost = - (double)(S_cost + SW_cost + SSW_cost + my_cost) / 4.0; + (double)((S_cost / 4.0) + (SW_cost / 4.0) + + (SSW_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (V_DIAG_fac * a) + lambda * fcost_cost * V_DIAG_fac; @@ -1283,7 +1378,8 @@ else fcost_dtm = (double)(WNW_dtm - my_dtm) * c; fcost_cost = - (double)(W_cost + NW_cost + WNW_cost + my_cost) / 4.0; + (double)((W_cost / 4.0) + (NW_cost / 4.0) + + (WNW_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) + lambda * fcost_cost * H_DIAG_fac; @@ -1302,7 +1398,8 @@ else fcost_dtm = (double)(ENE_dtm - my_dtm) * c; fcost_cost = - (double)(E_cost + NE_cost + ENE_cost + my_cost) / 4.0; + (double)((E_cost / 4.0) + (NE_cost / 4.0) + + (ENE_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) + lambda * fcost_cost * H_DIAG_fac; @@ -1321,7 +1418,8 @@ else fcost_dtm = (double)(ESE_dtm - my_dtm) * c; fcost_cost = - (double)(E_cost + SE_cost + ESE_cost + my_cost) / 4.0; + (double)((E_cost / 4.0) + (SE_cost / 4.0) + + (ESE_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) + lambda * fcost_cost * H_DIAG_fac; @@ -1340,7 +1438,8 @@ else fcost_dtm = (double)(WSW_dtm - my_dtm) * c; fcost_cost = - (double)(W_cost + SW_cost + WSW_cost + my_cost) / 4.0; + (double)((W_cost / 4.0) + (SW_cost / 4.0) + + (WSW_cost / 4.0) + (my_cost / 4.0)); min_cost = pres_cell->min_cost + fcost_dtm + (H_DIAG_fac * a) + lambda * fcost_cost * H_DIAG_fac; @@ -1351,15 +1450,21 @@ continue; segment_get(&out_seg, &old_min_cost, row, col); + if (dir == 1) + segment_get(&out_seg2, &old_cur_dir, row, col); if (G_is_d_null_value(&old_min_cost)) { segment_put(&out_seg, &min_cost, row, col); new_cell = insert(min_cost, row, col); + if (dir == 1) + segment_put(&out_seg2, &cur_dir, row, col); } else { if (old_min_cost > min_cost) { segment_put(&out_seg, &min_cost, row, col); new_cell = insert(min_cost, row, col); + if (dir == 1) + segment_put(&out_seg2, &cur_dir, row, col); } else { } @@ -1386,13 +1491,19 @@ cum_fd = G_open_raster_new(cum_cost_layer, cum_data_type); cum_cell = G_allocate_raster_buf(cum_data_type); + if (dir == 1) { + dir_fd = G_open_raster_new(move_dir_layer, dir_data_type); + dir_cell = G_allocate_raster_buf(dir_data_type); + } + /* Write pending updates by segment_put() to output map */ segment_flush(&out_seg); + if (dir == 1) + segment_flush(&out_seg2); /* Copy segmented map to output map */ - /* system("date"); */ G_message(_("Writing output raster map %s... "), cum_cost_layer); if (keep_nulls) { @@ -1506,7 +1617,22 @@ } } + double *p; + if (dir == 1) { + G_message(_("Writing movement direction file %s..."), move_dir_layer); + for (row = 0; row < nrows; row++) { + p = dir_cell; + for (col = 0; col < ncols; col++) { + segment_get(&out_seg2, &cur_dir, row, col); + *(p + col) = cur_dir; + } + G_put_raster_row(dir_fd, dir_cell, dir_data_type); + G_percent(row, nrows, 2); + } + } + + G_percent(row, nrows, 2); @@ -1514,15 +1640,23 @@ segment_release(&dtm_in_seg); /* release memory */ segment_release(&out_seg); + if (dir == 1) + segment_release(&out_seg2); G_close_cell(dtm_fd); G_close_cell(cost_fd); G_close_cell(cum_fd); + if (dir == 1) + G_close_cell(dir_fd); close(dtm_in_fd); /* close all files */ close(out_fd); close(cost_in_fd); + if (dir == 1) + close(dir_out_fd); unlink(dtm_in_file); /* remove submatrix files */ unlink(cost_in_file); unlink(out_file); + if (dir == 1) + unlink(dir_out_file); /* Create colours for output map */ @@ -1538,6 +1672,12 @@ G_command_history(&history); G_write_history(cum_cost_layer, &history); + if (dir == 1) { + G_short_history(move_dir_layer, "raster", &history); + G_command_history(&history); + G_write_history(move_dir_layer, &history); + } + exit(EXIT_SUCCESS); } Index: description.html =================================================================== --- description.html (revision 34032) +++ description.html (working copy) @@ -1,8 +1,10 @@

DESCRIPTION

-r.walk outputs a raster map layer showing the lowest +r.walk outputs 1) a raster map layer showing the lowest cumulative cost of moving between each cell and the user-specified -starting points. It uses an input elevation raster map layer whose +starting points and 2) a second raster map layer showing the movement +direction to the next cell on the path back to the start point (see +Movement Direction). It uses an input elevation raster map layer whose cell category values represent elevation, combined with a second input raster map layer whose cell values represent friction costs. @@ -69,11 +71,28 @@ The minimum cumulative costs are computed using Dijkstra's algorithm, that find an optimum solution (for more details see r.cost, that uses the same algorithm). + +

Movement Direction

+The movement direction surface is created to record the sequence of +movements that created the cost accumulation surface. Without it +r.drain would not correctly create a path from an end point +back to the start point. The direction shown in each cell points away +from the cell that came before it. 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
+
+

Once r.walk computes the cumulative cost map as a linear combination of friction cost (from friction map) and the altitude and distance covered (from the digital elevation model), r.drain -can be used to find the minimum cost path. +can be used to find the minimum cost path. Make sure to use the -d flag +and the movement direction raster map when running r.drain to ensure +the path is computed according to the proper movement directions.

SEE ALSO

@@ -102,15 +121,15 @@ Based on r.cost written by :

-Antony Awaida,
-Intelligent Engineering
-Systems Laboratory,
-M.I.T.
-
-James Westervelt,
+Antony Awaida,
+Intelligent Engineering
+Systems Laboratory,
+M.I.T.
+
+James Westervelt,
U.S.Army Construction Engineering Research Laboratory -

Updated for Grass 5
+

Updated for Grass 5
Pierre de Mouveaux (pmx@audiovu.com)

@@ -130,4 +149,7 @@

Roberto Flor and Markus Neteler -

Last changed: $Date$ +

+Updated to output a movement direction layer +

+Colin Nielsen, 2008 Index: stash.h =================================================================== --- stash.h (revision 34032) +++ stash.h (working copy) @@ -17,6 +17,7 @@ #define CUM_COST_LAYER 1 #define COST_LAYER 2 #define START_PT 3 +#define MOVE_DIR_LAYER 4 struct start_pt { @@ -38,17 +39,18 @@ { "output", CUM_COST_LAYER}, { "input", COST_LAYER}, { - "coor", START_PT} + "coor", START_PT}, { + "outdir", MOVE_DIR_LAYER} }; -char cum_cost_layer[64]; +char cum_cost_layer[64], move_dir_layer[64]; char cost_layer[64], dtm_layer[64]; struct start_pt *head_start_pt = NULL; struct start_pt *head_end_pt = NULL; #else -extern char cum_cost_layer[]; +extern char cum_cost_layer[], move_dir_layer[]; extern char cost_layer[], dtm_layer[]; extern struct start_pt *head_start_pt; extern struct start_pt *head_end_pt;