Index: main.c
===================================================================
--- main.c (revision 36142)
+++ main.c (working copy)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.cost
@@ -6,6 +7,9 @@
* James Westervelt - CERL
* Pierre de Mouveaux
* Eric G. Miller
+ *
+ * Updated for calculation errors and directional surface generation
+ * Colin Nielsen
*
* PURPOSE: Outputs a raster map layer showing the cumulative cost
* of moving between different geographic locations on an
@@ -30,9 +34,9 @@
* 1) Cost of traversing each grid cell as given by a cost map
* cell (input).
* 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.
*
*********************************************************************/
@@ -71,25 +75,26 @@
int main(int argc, char *argv[])
{
- void *cell, *cell2;
- SEGMENT in_seg, out_seg;
- char *cost_mapset;
- char *in_file, *out_file;
+ void *cell, *cell2, *dir_cell;
+ SEGMENT in_seg, out_seg, out_seg2;
+ char *cost_mapset, *move_dir_mapset;
+ char *in_file, *out_file, *dir_out_file;
char *search_mapset;
double *value;
extern struct Cell_head window;
double NS_fac, EW_fac, DIAG_fac, H_DIAG_fac, V_DIAG_fac;
double fcost;
double min_cost, old_min_cost;
+ double cur_dir, old_cur_dir;
double zero = 0.0;
int at_percent = 0;
int col, row, nrows, ncols;
int maxcost;
int maxmem;
double cost;
- int cost_fd, cum_fd;
- int have_stop_points = 0;
- int in_fd, out_fd;
+ int cost_fd, cum_fd, dir_fd;
+ int have_stop_points = 0, dir = 0;
+ int in_fd, out_fd, dir_out_fd;
double my_cost;
double null_cost;
int srows, scols;
@@ -103,13 +108,13 @@
struct GModule *module;
struct Flag *flag2, *flag3, *flag4;
struct Option *opt1, *opt2, *opt3, *opt4, *opt5, *opt6, *opt7, *opt8;
- struct Option *opt9, *opt10;
+ struct Option *opt9, *opt10, *opt11;
struct cost *pres_cell, *new_cell;
struct start_pt *pres_start_pt = NULL;
struct start_pt *pres_stop_pt = NULL;
void *ptr2;
- RASTER_MAP_TYPE data_type;
+ RASTER_MAP_TYPE data_type, dir_data_type = DCELL_TYPE;
struct History history;
double peak = 0.0;
int dsize;
@@ -133,6 +138,14 @@
opt1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt11 = G_define_option();
+ opt11->key = "outdir";
+ opt11->type = TYPE_STRING;
+ opt11->required = NO;
+ opt11->gisprompt = "new,cell,raster";
+ opt11->description =
+ _("Name of output raster map to contain movement directions");
+
opt7 = G_define_option();
opt7->key = "start_points";
opt7->type = TYPE_STRING;
@@ -211,25 +224,31 @@
flag4->description = _("Start with values in raster map");
/* please, remove before GRASS 7 released */
- v_flag = G_define_flag() ;
- v_flag->key = 'v' ;
- v_flag->description = _("Run verbosely") ;
+ v_flag = G_define_flag();
+ v_flag->key = 'v';
+ v_flag->description = _("Run verbosely");
/* Parse command line */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
/* please, remove before GRASS 7 released */
- if(v_flag->answer) {
- G_putenv("GRASS_VERBOSE","2");
- G_warning(_("The '-v' flag is superseded and will be removed "
- "in future. Please use '--verbose' instead."));
+ if (v_flag->answer) {
+ G_putenv("GRASS_VERBOSE", "2");
+ G_warning(_("The '-v' flag is superseded and will be removed "
+ "in future. Please use '--verbose' instead."));
}
+ /* If no outdir is specified, set flag to skip all dir */
+ if (opt11->answer != NULL)
+ dir = 1;
+
/* Initalize access to database and create temporary files */
in_file = G_tempfile();
out_file = G_tempfile();
+ if (dir == 1)
+ dir_out_file = G_tempfile();
/* Get database window parameters */
@@ -241,8 +260,10 @@
EW_fac = 1.0;
NS_fac = window.ns_res / window.ew_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);
@@ -281,16 +302,17 @@
if (sscanf(opt5->answer, "%d", &maxcost) != 1 || maxcost < 0)
G_fatal_error(_("Inappropriate maximum cost: %d"), maxcost);
- if (sscanf(opt10->answer, "%d", &maxmem) != 1 || maxmem < 0 || maxmem > 100)
+ if (sscanf(opt10->answer, "%d", &maxmem) != 1 || maxmem < 0 ||
+ maxmem > 100)
G_fatal_error(_("Inappropriate percent memory: %d"), maxmem);
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)
- G_message (_("Input null cell will be retained into output map"));
+ G_message(_("Input null cell will be retained into output map"));
if (opt7->answer) {
search_mapset = G_find_vector(opt7->answer, "");
@@ -300,7 +322,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);
}
}
@@ -318,11 +341,23 @@
if (cost_mapset == NULL)
G_fatal_error(_("Raster map <%s> not found"), cost_layer);
+ if (dir == 1) {
+ strcpy(move_dir_layer, opt11->answer);
+
+ search_mapset = "";
+ move_dir_mapset = G_find_cell2(move_dir_layer, search_mapset);
+ }
+
/* 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();
@@ -343,27 +378,30 @@
switch (data_type) {
case (CELL_TYPE):
- G_message(_("Source map is: Integer cell type"));
- break;
+ G_message(_("Source map is: Integer cell type"));
+ break;
case (FCELL_TYPE):
- G_message(_("Source map is: Floating point (float) cell type"));
- break;
+ G_message(_("Source map is: Floating point (float) cell type"));
+ break;
case (DCELL_TYPE):
- G_message(_("Source map is: Floating point (double) cell type"));
- break;
+ G_message(_("Source map is: Floating point (double) cell type"));
+ break;
}
G_message(_(" %d rows, %d cols"), nrows, ncols);
srows = scols = SEGCOLSIZE;
if (maxmem > 0)
segments_in_memory =
- 2 + maxmem * (nrows / SEGCOLSIZE) * (ncols / SEGCOLSIZE) / 100;
+ 2 + maxmem * (1 + nrows / SEGCOLSIZE) * (1 +
+ ncols / SEGCOLSIZE) /
+ 100;
else
- segments_in_memory = 4 * (nrows / SEGCOLSIZE + ncols / SEGCOLSIZE + 2);
+ segments_in_memory =
+ 4 * (nrows / SEGCOLSIZE + ncols / SEGCOLSIZE + 2);
/* Create segmented format files for cost layer and output layer */
- G_message(_("Creating some temporary files"));
+ G_message(_("Creating some temporary files..."));
in_fd = creat(in_file, 0666);
segment_format(in_fd, nrows, ncols, srows, scols, sizeof(double));
@@ -373,6 +411,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 */
in_fd = open(in_file, 2);
@@ -382,6 +427,12 @@
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"), cost_layer);
@@ -394,9 +445,10 @@
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, cell, row, 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 = cell;
@@ -448,7 +500,7 @@
/* Initialize segmented output file */
G_message(_("Initializing output "));
-
+
{
double *fbuff;
int i;
@@ -456,28 +508,59 @@
fbuff = (double *)G_malloc(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);
+ G_percent(row, nrows, 2);
for (i = 0; i < ncols; i++) {
segment_put(&out_seg, &fbuff[i], row, i);
}
}
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 start_points layer searching for starting points.
* Create a btree of starting points ordered by increasing costs.
*/
if (opt7->answer) {
#if 1
- FILE *fp;
+ struct Map_info *fp;
struct start_pt *new_start_pt;
Site *site = NULL; /* pointer to Site */
int got_one = 0;
@@ -488,9 +571,9 @@
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))
@@ -528,7 +611,7 @@
if (opt8->answer) {
#if 1
- FILE *fp;
+ struct Map_info *fp;
struct start_pt *new_start_pt;
Site *site = NULL; /* pointer to Site */
int dims, strs, dbls;
@@ -538,8 +621,8 @@
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);) {
@@ -596,16 +679,18 @@
if (!cell2)
G_fatal_error(_("Unable to allocate memory"));
- G_message(_("Reading %s"), opt9->answer);
+ G_message(_("Reading %s"), opt9->answer);
for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 2);
+ G_percent(row, nrows, 2);
if (G_get_raster_row(fd, cell2, row, data_type2) < 0)
- G_fatal_error(_("Unable to read raster map <%s> row %d"), opt9->answer, row);
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ opt9->answer, 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);
@@ -621,7 +706,7 @@
ptr2 = G_incr_void_ptr(ptr2, dsize2);
}
}
- G_percent(row, nrows, 2);
+ G_percent(row, nrows, 2);
G_close_cell(fd);
G_free(cell2);
@@ -636,6 +721,7 @@
*/
if (head_start_pt) {
struct start_pt *top_start_pt = NULL;
+
top_start_pt = head_start_pt;
while (top_start_pt != NULL) {
value = &zero;
@@ -644,7 +730,8 @@
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, top_start_pt->row, top_start_pt->col);
+ segment_put(&out_seg, value, top_start_pt->row,
+ top_start_pt->col);
top_start_pt = top_start_pt->next;
}
/* printf("--------+++++----------\n"); */
@@ -684,7 +771,7 @@
segment_get(&in_seg, &my_cost, pres_cell->row, pres_cell->col);
- 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).
@@ -697,55 +784,71 @@
case 1:
row = pres_cell->row;
col = pres_cell->col - 1;
+ cur_dir = 180.0;
break;
case 2:
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;
+ cur_dir = 270.0;
break;
case 5:
row = pres_cell->row - 1;
col = pres_cell->col - 1;
+ cur_dir = 135.0;
break;
case 6:
col = pres_cell->col + 1;
+ cur_dir = 45.0;
break;
case 7:
row = pres_cell->row + 1;
+ cur_dir = 315.0;
break;
case 8:
col = pres_cell->col - 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:
col = pres_cell->col + 1;
+ cur_dir = 67.5;
break;
case 11:
row = pres_cell->row + 2;
+ cur_dir = 292.5;
break;
case 12:
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:
col = pres_cell->col + 2;
+ cur_dir = 22.5;
break;
case 15:
row = pres_cell->row + 1;
+ cur_dir = 337.5;
break;
case 16:
col = pres_cell->col - 2;
+ cur_dir = 202.5;
break;
}
@@ -760,97 +863,113 @@
case 1:
value = &W;
segment_get(&in_seg, value, row, col);
- fcost = (double)(W + my_cost) / 2.0;
+ fcost = (double)((W / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * EW_fac;
break;
case 2:
value = &E;
segment_get(&in_seg, value, row, col);
- fcost = (double)(E + my_cost) / 2.0;
+ fcost = (double)((E / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * EW_fac;
break;
case 3:
value = &N;
segment_get(&in_seg, value, row, col);
- fcost = (double)(N + my_cost) / 2.0;
+ fcost = (double)((N / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * NS_fac;
break;
case 4:
value = &S;
segment_get(&in_seg, value, row, col);
- fcost = (double)(S + my_cost) / 2.0;
+ fcost = (double)((S / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * NS_fac;
break;
case 5:
value = &NW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(NW + my_cost) / 2.0;
+ fcost = (double)((NW / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * DIAG_fac;
break;
case 6:
value = &NE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(NE + my_cost) / 2.0;
+ fcost = (double)((NE / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * DIAG_fac;
break;
case 7:
value = &SE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(SE + my_cost) / 2.0;
+ fcost = (double)((SE / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * DIAG_fac;
break;
case 8:
value = &SW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(SW + my_cost) / 2.0;
+ fcost = (double)((SW / 2.0) + (my_cost / 2.0));
min_cost = pres_cell->min_cost + fcost * DIAG_fac;
break;
case 9:
value = &NNW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(N + NW + NNW + my_cost) / 4.0;
+ fcost =
+ (double)((N / 4.0) + (NW / 4.0) + (NNW / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * V_DIAG_fac;
break;
case 10:
value = &NNE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(N + NE + NNE + my_cost) / 4.0;
+ fcost =
+ (double)((N / 4.0) + (NE / 4.0) + (NNE / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * V_DIAG_fac;
break;
case 11:
value = &SSE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(S + SE + SSE + my_cost) / 4.0;
+ fcost =
+ (double)((S / 4.0) + (SE / 4.0) + (SSE / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * V_DIAG_fac;
break;
case 12:
value = &SSW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(S + SW + SSW + my_cost) / 4.0;
+ fcost =
+ (double)((S / 4.0) + (SW / 4.0) + (SSW / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * V_DIAG_fac;
break;
case 13:
value = &WNW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(W + NW + WNW + my_cost) / 4.0;
+ fcost =
+ (double)((W / 4.0) + (NW / 4.0) + (WNW / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * H_DIAG_fac;
break;
case 14:
value = &ENE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(E + NE + ENE + my_cost) / 4.0;
+ fcost =
+ (double)((E / 4.0) + (NE / 4.0) + (ENE / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * H_DIAG_fac;
break;
case 15:
value = &ESE;
segment_get(&in_seg, value, row, col);
- fcost = (double)(E + SE + ESE + my_cost) / 4.0;
+ fcost =
+ (double)((E / 4.0) + (SE / 4.0) + (ESE / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * H_DIAG_fac;
break;
case 16:
value = &WSW;
segment_get(&in_seg, value, row, col);
- fcost = (double)(W + SW + WSW + my_cost) / 4.0;
+ fcost =
+ (double)((W / 4.0) + (SW / 4.0) + (WSW / 4.0) +
+ (my_cost / 4.0));
min_cost = pres_cell->min_cost + fcost * H_DIAG_fac;
break;
}
@@ -859,15 +978,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 {
}
@@ -882,7 +1010,7 @@
pres_cell = get_lowest();
if (pres_cell == NULL) {
- G_message(_("End of map!"));
+ G_message(_("End of map!"));
goto OUT;
}
if (ct == pres_cell)
@@ -893,25 +1021,33 @@
cum_fd = G_open_raster_new(cum_cost_layer, 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 */
G_message(_("Writing %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(data_type);
}
- if (data_type == CELL_TYPE) {
+ if (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(cost_fd, cell2, row, data_type) < 0)
G_fatal_error(_("Error getting input null cells"));
@@ -941,9 +1077,10 @@
else if (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(cost_fd, cell2, row, data_type) < 0)
G_fatal_error(_("Error getting input null cells"));
@@ -973,9 +1110,10 @@
else if (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(cost_fd, cell2, row, data_type) < 0)
G_fatal_error(_("Error getting input null cells"));
@@ -1003,23 +1141,56 @@
}
}
+ 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);
G_message(_("Peak cost value: %f"), peak);
segment_release(&in_seg); /* release memory */
segment_release(&out_seg);
+ if (dir == 1) {
+ segment_release(&out_seg2);
+ }
G_close_cell(cost_fd);
G_close_cell(cum_fd);
+ if (dir == 1) {
+ G_close_cell(dir_fd);
+ }
close(in_fd); /* close all files */
close(out_fd);
+ if (dir == 1) {
+ close(dir_out_fd);
+ }
unlink(in_file); /* remove submatrix files */
unlink(out_file);
+ if (dir == 1) {
+ unlink(dir_out_file);
+ }
G_short_history(cum_cost_layer, "raster", &history);
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);
+ }
+
/* Create colours for output map */
/*
@@ -1049,12 +1220,13 @@
for (n = 0; *answers != NULL; answers += 2) {
if (!G_scan_easting(*answers, &east, G_projection()))
G_fatal_error(_("Illegal x coordinate <%s>"), *answers);
- if (!G_scan_northing(*(answers + 1), &north, G_projection()))
+ if (!G_scan_northing(*(answers + 1), &north, G_projection()))
G_fatal_error(_("Illegal y coordinate <%s>"), *(answers + 1));
if (east < window.west || east > window.east ||
north < window.south || north > window.north) {
- G_warning(_("Warning, ignoring point outside window: %.4f,%.4f"), east, north);
+ G_warning(_("Warning, ignoring point outside window: %.4f,%.4f"),
+ east, north);
continue;
}
else
Index: description.html
===================================================================
--- description.html (revision 36142)
+++ description.html (working copy)
@@ -1,69 +1,72 @@
-DESCRIPTION
+DESCRIPTION
-r.cost determines the cumulative cost of moving to each
-cell on a cost surface (the input raster map) from
+
r.cost determines the cumulative cost of moving to each
+cell on a cost surface (the input raster map) from
other user-specified cell(s) whose locations are specified by their
geographic coordinate(s). Each cell in the original cost surface map
will contain a category value which represents the cost of traversing
-that cell. r.cost will produce an output raster map in
+that cell. r.cost will produce 1) an output raster map in
which each cell contains the lowest total cost of traversing the
-space between each cell and the user-specified points. (Diagonal
+space between each cell and the user-specified points (diagonal
costs are multiplied by a factor that depends on the dimensions of
-the cell.) This program uses the current geographic region settings.
-The output map will be of the same data format as the input
-map, integer or floating point.
+the cell) 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). This program uses the current geographic region settings.
+The output map will be of the same data format as the input
+map, integer or floating point.
-OPTIONS
+OPTIONS
-The input name is the name of a raster map whose category values
-represent the surface cost. The output name is the name of the
-resultant raster map of cumulative cost.
+The input name is the name of a raster map whose category values
+represent the surface cost. The output name is the name of the
+resultant raster map of cumulative cost. The outdir name is the
+name of the resultant raster map of movement directions (see Movement Direction).
-
-r.cost can be run with three different methods of identifying the
+
+r.cost can be run with three different methods of identifying the
starting point(s). One or more points (geographic coordinate pairs) can be
-provided as specified coordinates on the command line, from a vector
+provided as specified coordinates on the command line, from a vector
points file, or from a raster map.
All non-NULL cells are considered to be starting points.
-Each x,y coordinate pair gives the geographic location of a
+Each x,y coordinate pair gives the geographic location of a
point from which the transportation cost should be figured. As many points as
desired can be entered by the user. These starting points can also be read
-from a vector points file through the start_sites option or from a
-raster map through the start_rast option.
-
+from a vector points file through the start_sites option or from a
+raster map through the start_rast option.
+
r.cost will stop cumulating costs when either max_cost is reached,
-or one of the stop points given with stop_coordinates is reached.
+or one of the stop points given with stop_coordinates is reached.
Alternatively, the stop points can be read from a vector points file with the
-stop_sites option. During execution, once the cumulative cost to all
+stop_sites option. During execution, once the cumulative cost to all
stopping points has been determined, processing stops.
Both sites read from a vector points file and those given on the command line
will be processed.
-
+
The null cells in the input map can be assigned a (positive floating
point) cost with the null_cost option.
-When input map null cells are given a cost with the null_cost
+When input map null cells are given a cost with the null_cost
option, the corresponding cells in the output map are no longer null
cells. By using the -n flag, the null cells of the input map are
-retained as null cells in the output map.
+retained as null cells in the output map.
-
+
As r.cost can run for a very long time, it can be useful to
use the -v verbose flag to track progress.
-
+
The Knight's move (-k flag) may be used to improve the accuracy of
the output. In the diagram below, the center location (O) represents a
grid cell from which cumulative distances are calculated. Those
neighbors marked with an X are always considered for cumulative cost
-updates. With the -k option, the neighbors marked with a K are
+updates. With the -k option, the neighbors marked with a K are
also considered.
-
+
. . . . . . . . . . . . . . .
. . . K . . K . . .
@@ -77,10 +80,10 @@
. . . K . . K . . .
. . . . . . . . . . . . . . .
-
+
Knight's move example:
-
+
Flat cost surface without (left pane) and with the knight's move (right pane).
@@ -91,37 +94,38 @@
-NULL CELLS
+NULL CELLS
-By default null cells in the input raster map are excluded from
+ By default null cells in the input raster map are excluded from
the algorithm, and thus retained in the output map.
-
-If one wants r.cost to transparently cross any region of null cells,
-the null_cost=0.0 option should be used. Then null cells just
+
+If one wants r.cost to transparently cross any region of null cells,
+the null_cost=0.0 option should be used. Then null cells just
propagate the adjacent costs. These cells can be retained as null cells in the
-output map by using the -n flag.
+output map by using the -n flag.
- NOTES
-Sometimes, when the differences among integer cell category values in the
-r.cost cumulative cost surface output are small, this
-cumulative cost output cannot accurately be used as input to r.drain
-(r.drain will output bad
+ NOTES
+
+Sometimes, when the differences among integer cell category values in the
+r.cost cumulative cost surface output are small, this
+cumulative cost output cannot accurately be used as input to r.drain
+(r.drain will output bad
results). This problem can be circumvented by making the differences
between cell category values in the cumulative cost output bigger. It
-is recommended that, if the output from r.cost is to be used
-as input to r.drain, the user
-multiply the input cost surface map to r.cost by the value
-of the map's cell resolution, before running r.cost. This
-can be done using r.mapcalc. The map
-resolution can be found using g.region.
+is recommended that, if the output from r.cost is to be used
+as input to r.drain, the user
+multiply the input cost surface map to r.cost by the value
+of the map's cell resolution, before running r.cost. This
+can be done using r.mapcalc. The map
+resolution can be found using g.region.
This problem doesn't arise with floating point maps.
- Algorithm notes
+Algorithm notes
The fundamental approach to calculating minimum travel cost is as
-follows: The user generates a raster map indicating the cost of
+follows: The user generates a raster map indicating the cost of
traversing each cell in the north-south and east-west directions.
This map, along with a set of starting points are submitted to
r.cost. The starting points are put into a list cells from which
@@ -132,13 +136,13 @@
of selecting the lowest cumulative cost cell, computing costs to the
neighbors, putting the neighbors on the list and removing the
originating cell from the list continues until the list is empty.
-
+
The most time consuming aspect of this algorithm is the management of
the list of cells for which cumulative costs have been at least
initially computed. r.cost uses a binary tree with an linked list
at each node in the tree for efficiently holding cells with identical
cumulative costs.
-
+
r.cost, like most all GRASS raster programs, is also made to be run on
maps larger that can fit in available computer memory. As the
algorithm works through the dynamic list of cells it can move almost
@@ -147,14 +151,14 @@
and from disk) as needed. This provides a virtual memory approach
optimally designed for 2-D raster maps.
The amount of map to hold in memory at one time can be controlled with the
-percent_memory option. For large maps this value will have to be set
+percent_memory option. For large maps this value will have to be set
to a lower value.
- EXAMPLES
+EXAMPLES
-Consider the following example:
-
+Consider the following example:
+
Input:
COST SURFACE
@@ -167,7 +171,7 @@
. . . . . . . . . . . . . . .
. 8 . 7 . 8 . 8 . 8 . 8 . 5 .
. . . . . . . . . . _____ . .
- . 8 . 8 . 1 . 1 . 5 | 3 | 9 .
+ . 8 . 8 . 1 . 1 . 5 | 3 | 9 .
. . . . . . . . . . |___| . .
. 8 . 1 . 1 . 2 . 5 . 3 . 9 .
. . . . . . . . . . . . . . .
@@ -175,41 +179,41 @@
Output (using -k): Output (not using -k):
CUMULATIVE COST SURFACE CUMULATIVE COST SURFACE
- . . . . . . . . . . . . . . . . . . . * * * * * . . . . . .
- . 21. 21. 20. 19. 17. 15. 14. . 22. 21* 21* 20* 17. 15. 14.
- . . . . . . . . . . . . . . . . . . . * * * * * . . . . . .
- . 20. 19. 22. 19. 15. 12. 11. . 20. 19. 22* 20* 15. 12. 11.
- . . . . . . . . . . . . . . . . . . . . . * * * * * . . . .
- . 22. 18. 17. 17. 12. 11. 9. . 22. 18. 17* 18* 13* 11. 9.
- . . . . . . . . . . . . . . . . . . . . . * * * * * . . . .
+ . . . . . . . . . . . . . . . . . . . * * * * * . . . . . .
+ . 21. 21. 20. 19. 17. 15. 14. . 22. 21* 21* 20* 17. 15. 14.
+ . . . . . . . . . . . . . . . . . . . * * * * * . . . . . .
+ . 20. 19. 22. 19. 15. 12. 11. . 20. 19. 22* 20* 15. 12. 11.
+ . . . . . . . . . . . . . . . . . . . . . * * * * * . . . .
+ . 22. 18. 17. 17. 12. 11. 9. . 22. 18. 17* 18* 13* 11. 9.
+ . . . . . . . . . . . . . . . . . . . . . * * * * * . . . .
. 21. 14. 13. 12. 8. 6. 6. . 21. 14. 13. 12. 8. 6. 6.
. . . . . . . . . . _____. . . . . . . . . . . . . . . . .
- . 16. 13. 8. 7. 4 | 0 | 6. . 16. 13. 8. 7 . 4. 0. 6.
+ . 16. 13. 8. 7. 4 | 0 | 6. . 16. 13. 8. 7 . 4. 0. 6.
. . . . . . . . . . |___|. . . . . . . . . . . . . . . . .
. 14. 9. 8. 9. 6. 3. 8. . 14. 9. 8. 9 . 6. 3. 8.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-
+
-The user-provided ending location in the above example is the boxed 3
+The user-provided ending location in the above example is the boxed 3
in the above input map. The costs in the output map represent the total
cost of moving from each box ("cell") to one or more (here,
only one) starting location(s). Cells surrounded by asterisks are
those that are different between operations using and not using the
-Knight's move (-k) option.
+Knight's move (-k) option.
- Output analysis
+Output analysis
The output map can be viewed, for example, as an elevation model in which
-the starting location(s) is/are the lowest point(s). Outputs from r.cost
-can be used as inputs to r.drain, in order
-to trace the least-cost path given by this model between any given cell
-and the r.cost starting location(s). The two programs, when
-used together, generate least-cost paths or corridors between any two
-map locations (cells).
+the starting location(s) is/are the lowest point(s). Outputs from r.cost
+can be used as inputs to r.drain with
+the direction flag -d, in order to trace the least-cost path given by this
+model between any given cell and the r.cost starting location(s). The
+two programs, when used together, generate least-cost paths or corridors between any
+two map locations (cells).
-Shortest distance surfaces
+Shortest distance surfaces
The r.cost module allows for computing the shortest distance
of each pixel from raster lines, such as determining the shortest distances
of households to the nearby road. For this cost surfaces with cost value 1 are
@@ -228,6 +232,26 @@
d.rast dist_meters
+
+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.cost computes the cumulative cost map, r.drain
+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.
BUGS
Index: stash.h
===================================================================
--- stash.h (revision 36142)
+++ stash.h (working copy)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.cost
@@ -28,38 +29,42 @@
#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];
+struct start_pt *head_start_pt = NULL;
+struct start_pt *head_end_pt = NULL;
- char cum_cost_layer[64];
- char cost_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[];
- 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[];
+extern struct start_pt *head_start_pt;
+extern struct start_pt *head_end_pt;
#endif
|