Index: main.c =================================================================== --- main.c (revision 36142) +++ main.c (working copy) @@ -1,3 +1,4 @@ + /**************************************************************************** * * MODULE: r.walk @@ -24,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 * @@ -43,9 +46,9 @@ * 1) Cost of traversing each grid cell as given by an elevation map and * a cost map cell (inputcost). * 2) If starting points are not specified on the command line - * then the ouput map must exist and contain the starting locations + * then the output map must exist and contain the starting locations * - * Otherwise the ouput map need not exist and the coor points + * Otherwise the output map need not exist and the coor points * from the command line are used. * *********************************************************************/ @@ -114,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]; @@ -127,16 +130,17 @@ 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 have_start_points; - int have_stop_points; - int dtm_in_fd, cost_in_fd, out_fd; + int have_stop_points, dir = 0; + 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; @@ -151,14 +155,15 @@ 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; struct start_pt *pres_stop_pt = NULL; void *ptr2; - RASTER_MAP_TYPE dtm_data_type, data_type2, cost_data_type, cum_data_type = DCELL_TYPE, cat; + RASTER_MAP_TYPE dtm_data_type, data_type2, cost_data_type, cum_data_type = + DCELL_TYPE, dir_data_type = DCELL_TYPE, cat; double peak = 0.0; int dtm_dsize, cost_dsize; @@ -200,6 +205,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; @@ -233,7 +246,6 @@ opt5 = G_define_option(); opt5->key = "max_cost"; opt5->type = TYPE_INTEGER; - opt5->key_desc = "cost"; opt5->required = NO; opt5->multiple = NO; opt5->answer = "0"; @@ -242,7 +254,6 @@ opt6 = G_define_option(); opt6->key = "null_cost"; opt6->type = TYPE_DOUBLE; - opt6->key_desc = "null cost"; opt6->required = NO; opt6->multiple = NO; opt6->description = @@ -251,7 +262,6 @@ opt9 = G_define_option(); opt9->key = "percent_memory"; opt9->type = TYPE_INTEGER; - opt9->key_desc = "percent memory"; opt9->required = NO; opt9->multiple = NO; opt9->answer = "100"; @@ -260,11 +270,11 @@ opt14 = G_define_option(); opt14->key = "nseg"; opt14->type = TYPE_INTEGER; - opt14->key_desc = "nseg"; opt14->required = NO; opt14->multiple = NO; opt14->answer = "4"; - opt14->description = _("Number of the segment to create (segment library)"); + opt14->description = + _("Number of the segment to create (segment library)"); opt10 = G_define_option(); opt10->key = "walk_coeff"; @@ -279,16 +289,16 @@ opt11 = G_define_option(); opt11->key = "lambda"; opt11->type = TYPE_DOUBLE; - opt11->key_desc = "lambda"; opt11->required = YES; opt11->multiple = NO; + opt11->answer = "0.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"; opt13->type = TYPE_DOUBLE; - opt13->key_desc = "slope_factor"; opt13->required = NO; opt13->multiple = NO; opt13->answer = "-0.2125"; @@ -312,11 +322,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 */ @@ -329,8 +345,10 @@ NS_fac = window.ns_res; DIAG_fac = (double)sqrt((double)(NS_fac * NS_fac + EW_fac * EW_fac)); - V_DIAG_fac = (double)sqrt((double)(4 * NS_fac * NS_fac + EW_fac * EW_fac)); - H_DIAG_fac = (double)sqrt((double)(NS_fac * NS_fac + 4 * EW_fac * EW_fac)); + V_DIAG_fac = + (double)sqrt((double)(4 * NS_fac * NS_fac + EW_fac * EW_fac)); + H_DIAG_fac = + (double)sqrt((double)(NS_fac * NS_fac + 4 * EW_fac * EW_fac)); G_set_d_null_value(&null_cost, 1); @@ -350,7 +368,8 @@ if (sscanf(opt5->answer, "%d", &maxcost) != 1 || maxcost < 0) G_fatal_error(_("Inappropriate maximum cost: %d"), maxcost); - if (sscanf(opt9->answer, "%d", &maxmem) != 1 || maxmem < 0 || maxmem > 100) + if (sscanf(opt9->answer, "%d", &maxmem) != 1 || maxmem < 0 || + maxmem > 100) G_fatal_error(_("Inappropriate percent memory: %d"), maxmem); /* Getting walking energy formula parameters */ @@ -359,38 +378,37 @@ G_fatal_error(_("Missing required value: got %d instead of 4"), par_number); else { - - G_message(_("Walking costs are a=%lf b=%lf c=%lf d=%lf"), a, b, c, - d); + + G_message(_("Walking costs are a=%lf b=%lf c=%lf d=%lf"), a, b, c, d); } /* Getting lambda */ if ((par_number = sscanf(opt11->answer, "%lf", &lambda)) != 1) G_fatal_error(_("Missing required value: %d"), par_number); else { - - G_message(_("Lambda is %lf"), lambda); + + G_message(_("Lambda is %lf"), lambda); } /*Getting slope_factor */ if ((par_number = sscanf(opt13->answer, "%lf", &slope_factor)) != 1) G_fatal_error(_("Missing required value: %d"), par_number); else { - - G_message(_("Slope_factor is %lf"), slope_factor); + + G_message(_("Slope_factor is %lf"), slope_factor); } if ((par_number = sscanf(opt14->answer, "%d", &nseg)) != 1) G_fatal_error(_("Missing required value: %d"), par_number); else { - - G_message(_("Nseg is %d"), nseg); + + G_message(_("Nseg is %d"), nseg); } if ((opt6->answer == NULL) || (sscanf(opt6->answer, "%lf", &null_cost) != 1)) { - - G_message(_("Null cells excluded from cost evaluation.")); + + G_message(_("Null cells excluded from cost evaluation.")); G_set_d_null_value(&null_cost, 1); } else if (keep_nulls) @@ -398,7 +416,7 @@ if (opt7->answer) { - FILE *fp; + struct Map_info *fp; struct start_pt *new_start_pt; Site *site = NULL; /* pointer to Site */ int dims, strs, dbls; @@ -406,12 +424,13 @@ search_mapset = ""; search_mapset = G_find_sites(opt7->answer, ""); if (search_mapset == NULL) - G_fatal_error(_("Unable to find starting vector %s "), opt7->answer); + G_fatal_error(_("Unable to find starting vector <%s> "), + opt7->answer); fp = G_fopen_sites_old(opt7->answer, search_mapset); - if (G_site_describe( fp, &dims, &cat, &strs, &dbls)) - G_fatal_error( "Failed to guess site file format\n"); - site = G_site_new_struct(cat, dims, strs, dbls); + if (G_site_describe(fp, &dims, &cat, &strs, &dbls)) + G_fatal_error("Failed to guess site file format\n"); + site = G_site_new_struct(cat, dims, strs, dbls); for (; (G_site_get(fp, site) != EOF);) { if (!G_site_in_region(site, &window)) @@ -444,7 +463,7 @@ } if (opt8->answer) { - FILE *fp; + struct Map_info *fp; struct start_pt *new_start_pt; Site *site = NULL; /* pointer to Site */ int dims, strs, dbls; @@ -452,11 +471,11 @@ search_mapset = ""; search_mapset = G_find_sites(opt8->answer, ""); if (search_mapset == NULL) - G_fatal_error(_("Unable to find stop vector %s"), opt8->answer); + G_fatal_error(_("Unable to find stop vector <%s>"), opt8->answer); fp = G_fopen_sites_old(opt8->answer, search_mapset); - if (G_site_describe( fp, &dims, &cat, &strs, &dbls)) - G_fatal_error( "Failed to guess site file format\n"); + if (G_site_describe(fp, &dims, &cat, &strs, &dbls)) + G_fatal_error("Failed to guess site file format\n"); site = G_site_new_struct(cat, dims, strs, dbls); for (; (G_site_get(fp, site) != EOF);) { @@ -508,7 +527,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); @@ -520,16 +545,21 @@ cost_mapset = G_find_cell2(cost_layer, search_mapset); if (dtm_mapset == NULL) - G_fatal_error(_("%s - not found"), dtm_layer); + G_fatal_error(_("Raster map <%s> not found"), dtm_layer); if (cost_mapset == NULL) - G_fatal_error(_("%s - not found"), cost_layer); + G_fatal_error(_("Raster map <%s> not found"), cost_layer); /* Check if specified output layer name is legal */ 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(); @@ -553,9 +583,9 @@ /*Reading headers from maps */ if (!dtm_head_ok) - G_fatal_error(_("Cannot read %s"), dtm_layer); + G_fatal_error(_("Unable to read %s"), dtm_layer); if (!cost_head_ok) - G_fatal_error(_("Cannot read %s"), cost_layer); + G_fatal_error(_("Unable to read %s"), cost_layer); /*Projection */ @@ -569,38 +599,35 @@ /* Parameters for map submatrices */ - - switch (dtm_data_type) { - case (CELL_TYPE): - G_message(_("DTM_Source map is: Integer cell type")); - break; - case (FCELL_TYPE): - G_message(_("DTM_Source map is: Floating point (float) cell type")); - break; - case (DCELL_TYPE): - G_message(_ - ("DTM_Source map is: Floating point (double) cell type")); - break; - } - G_message(_(" %d rows, %d cols"), dtm_cellhd.rows, dtm_cellhd.cols); - - - switch (cost_data_type) { - case (CELL_TYPE): - G_message(_("COST_Source map is: Integer cell type")); - break; - case (FCELL_TYPE): - G_message(_ - ("COST_Source map is: Floating point (float) cell type")); - break; - case (DCELL_TYPE): - G_message(_ - ("COST_Source map is: Floating point (double) cell type")); - break; - } - G_message(_(" %d rows, %d cols"), cost_cellhd.rows, cost_cellhd.cols); - + switch (dtm_data_type) { + case (CELL_TYPE): + G_message(_("DTM_Source map is: Integer cell type")); + break; + case (FCELL_TYPE): + G_message(_("DTM_Source map is: Floating point (float) cell type")); + break; + case (DCELL_TYPE): + G_message(_("DTM_Source map is: Floating point (double) cell type")); + break; + } + G_message(_(" %d rows, %d cols"), dtm_cellhd.rows, dtm_cellhd.cols); + + + + switch (cost_data_type) { + case (CELL_TYPE): + G_message(_("COST_Source map is: Integer cell type")); + break; + case (FCELL_TYPE): + G_message(_("COST_Source map is: Floating point (float) cell type")); + break; + case (DCELL_TYPE): + G_message(_("COST_Source map is: Floating point (double) cell type")); + break; + } + G_message(_(" %d rows, %d cols"), cost_cellhd.rows, cost_cellhd.cols); + if (cost_data_type != dtm_data_type) { switch (cost_data_type) { case (CELL_TYPE): @@ -624,25 +651,25 @@ /* Data type are equal, it doesn't matter */ cum_data_type = dtm_data_type; - - switch (cum_data_type) { - case (CELL_TYPE): - G_message(_("Output map is: Integer cell type")); - break; - case (FCELL_TYPE): - G_message(_("Output map is: Floating point (float) cell type")); - break; - case (DCELL_TYPE): - G_message(_("Output map is: Floating point (double) cell type")); - break; - } - G_message(_(" %d rows, %d cols"), nrows, ncols); - G_format_resolution(window.ew_res, buf, window.proj); - G_message(_(" EW resolution %s (%lf)"), buf, window.ew_res); - G_format_resolution(window.ns_res, buf, window.proj); - G_message(_(" NS resolution %s (%lf)"), buf, window.ns_res); - + switch (cum_data_type) { + case (CELL_TYPE): + G_message(_("Output map is: Integer cell type")); + break; + case (FCELL_TYPE): + G_message(_("Output map is: Floating point (float) cell type")); + break; + case (DCELL_TYPE): + G_message(_("Output map is: Floating point (double) cell type")); + break; + } + G_message(_(" %d rows, %d cols"), nrows, ncols); + G_format_resolution(window.ew_res, buf, window.proj); + G_message(_(" EW resolution %s (%lf)"), buf, window.ew_res); + G_format_resolution(window.ns_res, buf, window.proj); + G_message(_(" NS resolution %s (%lf)"), buf, window.ns_res); + + srows = nrows / nseg + 1; scols = ncols / nseg + 1; if (maxmem > 0) @@ -653,9 +680,9 @@ /* Create segmented format files for cost layer and output layer */ - - G_message(_("Creating some temporary files ...")); + G_message(_("Creating some temporary files...")); + dtm_in_fd = creat(dtm_in_file, 0600); segment_format(dtm_in_fd, nrows, ncols, srows, scols, sizeof(double)); close(dtm_in_fd); @@ -668,6 +695,12 @@ 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); @@ -679,11 +712,16 @@ 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 */ - - G_message(_("Reading %s ..."), dtm_layer); + G_message(_("Reading %s..."), dtm_layer); + start_with_raster_vals = flag4->answer; { @@ -694,10 +732,11 @@ p = 0.0; for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (G_get_raster_row(dtm_fd, dtm_cell, row, dtm_data_type) < 0) - G_fatal_error(_("Unable to read raster map <%s> row %d"), dtm_layer, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + dtm_layer, row); /* INPUT NULL VALUES: ??? */ ptr2 = dtm_cell; switch (dtm_data_type) { @@ -742,9 +781,9 @@ } } - - G_message(_("Reading %s..."), cost_layer); + G_message(_("Reading %s..."), cost_layer); + { int i; double p; @@ -752,10 +791,11 @@ cost_dsize = G_raster_size(cost_data_type); p = 0.0; for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (G_get_raster_row(cost_fd, cost_cell, row, cost_data_type) < 0) - G_fatal_error(_("Unable to read raster map <%s> row %d"), cost_layer, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + cost_layer, row); /* INPUT NULL VALUES: ??? */ ptr2 = cost_cell; switch (cost_data_type) { @@ -803,14 +843,14 @@ segment_flush(&dtm_in_seg); segment_flush(&cost_in_seg); - - G_percent(row, nrows, 2); + G_percent(row, nrows, 2); + /* Initialize output map with NULL VALUES */ /* Initialize segmented output file */ - - G_message(_("Initializing output ")); + + G_message(_("Initializing output ")); { double *fbuff; int i; @@ -818,12 +858,13 @@ 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); for (row = 0; row < nrows; row++) { - { + { G_percent(row, nrows, 2); } for (i = 0; i < ncols; i++) { @@ -832,27 +873,57 @@ } segment_flush(&out_seg); - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); 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. */ if (!have_start_points) { int dsize2; + 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)"), + ("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); if (cum_fd < 0) - G_fatal_error(_("Unable to open raster map <%s>"), cum_cost_layer); + G_fatal_error(_("Unable to open raster map <%s>"), + cum_cost_layer); data_type2 = G_get_raster_map_type(cum_fd); @@ -864,18 +935,20 @@ G_fatal_error(_ ("Memory allocation error on reading start points from raster map %s"), cum_cost_layer); - - G_message(_("Reading %s ... "), cum_cost_layer); + + G_message(_("Reading %s... "), cum_cost_layer); for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (G_get_raster_row(cum_fd, cell2, row, data_type2) < 0) - G_fatal_error(_("Unable to read raster map <%s> row %d"), cum_cost_layer, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + cum_cost_layer, row); ptr2 = cell2; for (col = 0; col < ncols; col++) { /* Did I understand that concept of cummulative cost map? - (pmx) 12 april 2000 */ if (!G_is_null_value(ptr2, data_type2)) { double cellval; + if (start_with_raster_vals == 1) { cellval = G_get_raster_value_d(ptr2, data_type2); new_cell = insert(cellval, row, col); @@ -890,9 +963,9 @@ ptr2 = G_incr_void_ptr(ptr2, dsize2); } } - - G_percent(row, nrows, 2); + G_percent(row, nrows, 2); + G_close_cell(cum_fd); G_free(cell2); @@ -904,6 +977,7 @@ */ else { struct start_pt *top_start_pt = NULL; + top_start_pt = head_start_pt; while (top_start_pt != NULL) { value_start_pt = &zero; @@ -925,10 +999,9 @@ * 3) Free the memory allocated to the present cell. */ - - /*system("date");*/ - G_message(_("Finding cost path")); - + + G_message(_("Finding cost path")); + n_processed = 0; total_cells = nrows * ncols; at_percent = 0; @@ -965,9 +1038,9 @@ if (G_is_d_null_value(&my_cost)) continue; - - G_percent(++n_processed, total_cells, 1); + G_percent(++n_processed, total_cells, 1); + /* 9 10 Order in which neighbors * 13 5 3 6 14 are visited (Knight move). * 1 2 @@ -979,66 +1052,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; } @@ -1061,7 +1150,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; @@ -1079,7 +1168,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; @@ -1097,7 +1186,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; @@ -1115,7 +1204,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; @@ -1133,7 +1222,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; @@ -1151,7 +1240,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; @@ -1169,7 +1258,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; @@ -1187,7 +1276,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; @@ -1206,7 +1295,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; @@ -1225,7 +1315,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; @@ -1244,7 +1335,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; @@ -1263,7 +1355,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; @@ -1282,7 +1375,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; @@ -1301,7 +1395,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; @@ -1320,7 +1415,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; @@ -1339,7 +1435,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; @@ -1350,15 +1447,24 @@ 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 { } @@ -1373,8 +1479,8 @@ pres_cell = get_lowest(); if (pres_cell == NULL) { - - G_message(_("End of map!")); + + G_message(_("End of map!")); goto OUT; } if (ct == pres_cell) @@ -1385,32 +1491,39 @@ 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); - + + G_message(_("Writing output raster map %s... "), cum_cost_layer); + if (keep_nulls) { - - G_message(_("Will copy input map null values into output map")); + + G_message(_("Will copy input map null values into output map")); cell2 = G_allocate_raster_buf(dtm_data_type); } if (cum_data_type == CELL_TYPE) { int *p; int *p2; - - G_message(_("Integer cell type.\nWriting...")); + + G_message(_("Integer cell type.\nWriting...")); for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (keep_nulls) { if (G_get_raster_row(dtm_fd, cell2, row, dtm_data_type) < 0) - G_fatal_error(_ - ("Unable to read raster map <%s> row %d"), dtm_layer, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + dtm_layer, row); } p = cum_cell; p2 = cell2; @@ -1437,15 +1550,15 @@ else if (cum_data_type == FCELL_TYPE) { float *p; float *p2; - - G_message(_("Float cell type.\nWriting...")); + + G_message(_("Float cell type.\nWriting...")); for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (keep_nulls) { if (G_get_raster_row(dtm_fd, cell2, row, dtm_data_type) < 0) - G_fatal_error(_ - ("Unable to read raster map <%s> row %d"), dtm_layer, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + dtm_layer, row); } p = cum_cell; p2 = cell2; @@ -1472,14 +1585,15 @@ else if (cum_data_type == DCELL_TYPE) { double *p; double *p2; - - G_message(_("Double cell type.\nWriting...")); + + G_message(_("Double cell type.\nWriting...")); for (row = 0; row < nrows; row++) { - - G_percent(row, nrows, 2); + + G_percent(row, nrows, 2); if (keep_nulls) { if (G_get_raster_row(dtm_fd, cell2, row, dtm_data_type) < 0) - G_fatal_error(_("Unable to read raster map <%s> row %d"), cell2, row); + G_fatal_error(_("Unable to read raster map <%s> row %d"), + cell2, row); } p = cum_cell; p2 = cell2; @@ -1504,23 +1618,50 @@ } } - + 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_message(_("Peak cost value: %f"), peak); + G_percent(row, nrows, 2); + + + G_message(_("Peak cost value: %f"), peak); + 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 */ @@ -1536,6 +1677,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 36142) +++ 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

Index: stash.h =================================================================== --- stash.h (revision 36142) +++ stash.h (working copy) @@ -1,3 +1,4 @@ + /***************************************************************/ /* */ /* stash.h in ~/src/Gcost */ @@ -6,6 +7,7 @@ /* the structures that are to be used for command */ /* line processing. */ /* */ + /***************************************************************/ #ifndef __R_COST_STASH_H__ @@ -15,39 +17,43 @@ #define CUM_COST_LAYER 1 #define COST_LAYER 2 #define START_PT 3 +#define MOVE_DIR_LAYER 4 -struct start_pt{ - int row; - int col; - struct start_pt *next; +struct start_pt +{ + int row; + int col; + struct start_pt *next; }; #ifdef MAIN - struct variables +struct variables +{ + char *alias; + int position; +} + +variables[] = { + { - char *alias; - int position; - } + "output", CUM_COST_LAYER}, { + "input", COST_LAYER}, { + "coor", START_PT}, { + "outdir", MOVE_DIR_LAYER} +}; - variables [] = { - - {"output",CUM_COST_LAYER}, - {"input",COST_LAYER}, - {"coor",START_PT} - }; +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; - char cum_cost_layer[64]; - char cost_layer[64],dtm_layer[64]; - struct start_pt *head_start_pt = NULL; - struct start_pt *head_end_pt = NULL; - -#else +#else - extern char cum_cost_layer[]; - extern char cost_layer[],dtm_layer[]; - extern struct start_pt *head_start_pt; - extern struct start_pt *head_end_pt; +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; #endif @@ -55,5 +61,5 @@ int time_to_stop(int, int); #endif -/****************END OF "GCOST_CMD_LINE.H"**********************/ +/****************END OF "GCOST_CMD_LINE.H"**********************/