[GRASS-SVN] r36350 - grass/branches/develbranch_6/raster/r.walk
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 13 02:14:10 EDT 2009
Author: cmbarton
Date: 2009-03-13 02:14:05 -0400 (Fri, 13 Mar 2009)
New Revision: 36350
Modified:
grass/branches/develbranch_6/raster/r.walk/description.html
grass/branches/develbranch_6/raster/r.walk/main.c
Log:
Adds movement directions map option for more accurate least cost path generation.
Modified: grass/branches/develbranch_6/raster/r.walk/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.walk/description.html 2009-03-13 02:39:21 UTC (rev 36349)
+++ grass/branches/develbranch_6/raster/r.walk/description.html 2009-03-13 06:14:05 UTC (rev 36350)
@@ -1,8 +1,10 @@
<h2>DESCRIPTION</h2>
-<em>r.walk</em> outputs a raster map layer showing the lowest
+<em>r.walk</em> 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
+<a href="#move">Movement Direction</a>). 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
<em>r.cost</em>, that uses the same algorithm).
+<a name="move"></a>
+<H2>Movement Direction</H2>
<p>
+The movement direction surface is created to record the sequence of
+movements that created the cost accumulation surface. Without it
+<EM>r.drain</EM> would not correctly create a path from an end point
+back to the start point. The direction shown in each cell points <b>away</b>
+from the cell that came before it. The directions are recorded as
+GRASS standard directions:<div class="code"><pre>
+ 112.5 90 67.5 i.e. a cell with the value 135
+157.5 135 0 45 22.5 means the cell <b>before</b> it is
+ 180 x 0 to the south-east.
+202.5 225 270 315 337.5
+ 247.5 292.5
+</pre></div>
+<p>
Once <em>r.walk</em> 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), <em>r.drain</em>
-can be used to find the minimum cost path.
+can be used to find the minimum cost path. Make sure to use the <b>-d</b> flag
+and the movement direction raster map when running r.drain to ensure
+the path is computed according to the proper movement directions.
<h2>SEE ALSO</h2>
Modified: grass/branches/develbranch_6/raster/r.walk/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.walk/main.c 2009-03-13 02:39:21 UTC (rev 36349)
+++ grass/branches/develbranch_6/raster/r.walk/main.c 2009-03-13 06:14:05 UTC (rev 36350)
@@ -25,6 +25,8 @@
* Updated for GRASS 6.1
* Roberto Flor and Markus Neteler
* Glynn Clements <glynn gclements.plus.com>, Soeren Gebbert <soeren.gebbert gmx.de>
+ * Updated for calculation errors and directional surface generation
+ * Colin Nielsen <colin.nielsen gmail com>
* 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,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 have_start_points;
+ int cost_fd, cum_fd, dtm_fd, dir_fd;
+ int have_start_points, dir = 0;
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 +155,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 +163,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 +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;
@@ -310,11 +321,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 */
@@ -508,7 +525,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 +553,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 +693,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 +711,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 */
@@ -836,6 +876,34 @@
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.
*/
@@ -980,66 +1048,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 +1146,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 +1164,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 +1182,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 +1200,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 +1218,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 +1236,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 +1254,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 +1272,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 +1291,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 +1311,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 +1331,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 +1351,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 +1371,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 +1391,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 +1411,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 +1431,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 +1443,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 {
}
@@ -1386,9 +1487,17 @@
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 */
@@ -1506,7 +1615,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 +1638,27 @@
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 +1674,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);
}
More information about the grass-commit
mailing list