[GRASS-SVN] r54470 - in grass/branches/releasebranch_6_4/raster: r.cost r.walk
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 31 04:29:35 PST 2012
Author: mmetz
Date: 2012-12-31 04:29:35 -0800 (Mon, 31 Dec 2012)
New Revision: 54470
Modified:
grass/branches/releasebranch_6_4/raster/r.cost/description.html
grass/branches/releasebranch_6_4/raster/r.cost/main.c
grass/branches/releasebranch_6_4/raster/r.cost/stash.h
grass/branches/releasebranch_6_4/raster/r.walk/description.html
grass/branches/releasebranch_6_4/raster/r.walk/main.c
grass/branches/releasebranch_6_4/raster/r.walk/stash.h
Log:
r.cost / r.walk: backport direction output
Modified: grass/branches/releasebranch_6_4/raster/r.cost/description.html
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.cost/description.html 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.cost/description.html 2012-12-31 12:29:35 UTC (rev 54470)
@@ -6,11 +6,13 @@
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. <em>r.cost</em> will produce an <b>output</b> raster map in
+that cell. <em>r.cost</em> will produce 1) an <b>output</b> 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 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
+<a href="#move">Movement Direction</a>). This program uses the current geographic region settings.
The <b>output</b> map will be of the same data format as the <b>input</b>
map, integer or floating point.
@@ -18,7 +20,8 @@
The <b>input</b> <em>name</em> is the name of a raster map whose category values
represent the surface cost. The <b>output</b> <em>name</em> is the name of the
-resultant raster map of cumulative cost.
+resultant raster map of cumulative cost. The <b>outdir</b> <em>name</em> is the
+name of the resultant raster map of movement directions (see <a href="#move">Movement Direction</a>).
<p>
<em>r.cost</em> can be run with three different methods of identifying the
@@ -203,11 +206,11 @@
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 <em>r.cost</em>
-can be used as inputs to <em><a href="r.drain.html">r.drain</a></em>, in order
-to trace the least-cost path given by this model between any given cell
-and the <em>r.cost</em> starting location(s). The two programs, when
-used together, generate least-cost paths or corridors between any two
-map locations (cells).
+can be used as inputs to <em><a href="r.drain.html">r.drain</a></em> with
+the direction flag <b>-d</b>, in order to trace the least-cost path given by this
+model between any given cell and the <em>r.cost</em> starting location(s). The
+two programs, when used together, generate least-cost paths or corridors between any
+two map locations (cells).
<h3>Shortest distance surfaces</h3>
The <em>r.cost</em> module allows for computing the shortest distance
@@ -228,6 +231,26 @@
d.rast dist_meters
</pre></div>
+<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.cost</em> computes the cumulative cost map, <em>r.drain</em>
+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>BUGS</h2>
Modified: grass/branches/releasebranch_6_4/raster/r.cost/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.cost/main.c 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.cost/main.c 2012-12-31 12:29:35 UTC (rev 54470)
@@ -8,6 +8,9 @@
* Pierre de Mouveaux <pmx audiovu com>
* Eric G. Miller <egm2 jps net>
*
+ * Updated for calculation errors and directional surface generation
+ * Colin Nielsen <colin.nielsen gmail com>
+ *
* PURPOSE: Outputs a raster map layer showing the cumulative cost
* of moving between different geographic locations on an
* input raster map layer whose cell category values
@@ -72,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;
@@ -104,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;
@@ -134,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_standard_option(G_OPT_V_INPUT);
opt7->key = "start_points";
opt7->required = NO;
@@ -227,10 +239,16 @@
"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 */
@@ -321,11 +339,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();
@@ -379,6 +409,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);
@@ -388,6 +425,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 */
G_message(_("Reading raster map <%s>..."),
@@ -494,6 +536,31 @@
G_free(fbuff);
}
+ if (dir == 1) {
+ G_message(_("Initializing directional output "));
+ {
+ double *fbuff;
+ int i;
+
+ fbuff =
+ (double *)G_malloc((unsigned int)(ncols * sizeof(double)));
+
+ 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);
+ }
+ }
+ G_percent(1, 1, 1);
+ segment_flush(&out_seg2);
+ G_free(fbuff);
+ }
+ }
+
/* Scan the start_points layer searching for starting points.
* Create a btree of starting points ordered by increasing costs.
*/
@@ -721,55 +788,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;
}
@@ -883,15 +966,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 {
}
@@ -917,9 +1009,17 @@
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 raster map <%s>..."), cum_cost_layer);
@@ -1029,19 +1129,52 @@
}
G_percent(1, 1, 1);
+ if (dir == 1) {
+ G_message(_("Writing movement direction file %s..."), move_dir_layer);
+ for (row = 0; row < nrows; row++) {
+ double *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(1, 1, 1);
+
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 */
/*
@@ -1052,7 +1185,7 @@
*/
G_done_msg(_("Peak cost value: %f."), peak);
-
+
exit(EXIT_SUCCESS);
}
Modified: grass/branches/releasebranch_6_4/raster/r.cost/stash.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.cost/stash.h 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.cost/stash.h 2012-12-31 12:29:35 UTC (rev 54470)
@@ -29,6 +29,7 @@
#define CUM_COST_LAYER 1
#define COST_LAYER 2
#define START_PT 3
+#define MOVE_DIR_LAYER 4
struct start_pt
{
@@ -49,17 +50,18 @@
{
"output", CUM_COST_LAYER}, {
"input", COST_LAYER}, {
- "coor", START_PT}
+ "coor", START_PT}, {
+ "outdir", MOVE_DIR_LAYER}
};
-char cum_cost_layer[GNAME_MAX];
+char cum_cost_layer[GNAME_MAX], move_dir_layer[GNAME_MAX];
char cost_layer[GNAME_MAX];
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[];
extern struct start_pt *head_start_pt;
extern struct start_pt *head_end_pt;
Modified: grass/branches/releasebranch_6_4/raster/r.walk/description.html
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.walk/description.html 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.walk/description.html 2012-12-31 12:29:35 UTC (rev 54470)
@@ -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/releasebranch_6_4/raster/r.walk/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.walk/main.c 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.walk/main.c 2012-12-31 12:29:35 UTC (rev 54470)
@@ -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;
}
@@ -1351,15 +1435,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 +1479,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 +1607,21 @@
}
}
+ if (dir == 1) {
+ G_message(_("Writing movement direction file %s..."), move_dir_layer);
+ for (row = 0; row < nrows; row++) {
+ double *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 +1629,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 +1665,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);
}
Modified: grass/branches/releasebranch_6_4/raster/r.walk/stash.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.walk/stash.h 2012-12-30 17:23:10 UTC (rev 54469)
+++ grass/branches/releasebranch_6_4/raster/r.walk/stash.h 2012-12-31 12:29:35 UTC (rev 54470)
@@ -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 cost_layer[64], dtm_layer[64];
+char cum_cost_layer[GNAME_MAX], move_dir_layer[GNAME_MAX];
+char cost_layer[GNAME_MAX], dtm_layer[GNAME_MAX];
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;
More information about the grass-commit
mailing list