[GRASS-SVN] r36351 - grass/branches/develbranch_6/raster/r.drain
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 13 02:14:59 EDT 2009
Author: cmbarton
Date: 2009-03-13 02:14:56 -0400 (Fri, 13 Mar 2009)
New Revision: 36351
Modified:
grass/branches/develbranch_6/raster/r.drain/description.html
grass/branches/develbranch_6/raster/r.drain/main.c
Log:
Adds movement directions map option for more accurate least cost path generation.
Modified: grass/branches/develbranch_6/raster/r.drain/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.drain/description.html 2009-03-13 06:14:05 UTC (rev 36350)
+++ grass/branches/develbranch_6/raster/r.drain/description.html 2009-03-13 06:14:56 UTC (rev 36351)
@@ -2,15 +2,17 @@
<em>r.drain</em> traces a flow through a least-cost path in an elevation
-model. The <b>input</b> elevation surface (a raster map layer) might
-be a cumulative cost map generated by the
-<em><a href="r.cost.html">r.cost</a></em> program.
+model. If the <b>input</b> surface (a raster map layer) is a cumulative
+cost map generated by the <em><a href="r.walk.html">r.walk</a></em> or
+<em><a href="r.cost.html">r.cost</a></em> modules, the -d flag and a
+movement direction surface <b>"indir"</b> must be specified.
The <b>output</b> result (also a raster map layer) will show one or more
least-cost paths between each user-provided location(s) and the minima
-(low category values) in the <b>input</b> map. By default, the <b>output</b>
-will be an integer CELL map with <tt>1</tt> along the least cost path,
-and null cells elsewhere.
+(low category values) in the <b>input</b> map. If the -d flag is used the
+output least-cost paths will be found using the direction layer.
+By default, the <b>output</b> will be an integer CELL map with <tt>1</tt>
+along the least cost path, and null cells elsewhere.
<p>
With the <b>-c</b> (<em>copy</em>) flag, the input map cell values are
@@ -24,9 +26,11 @@
The <b>-c</b>, <b>-a</b>, and <b>-n</b> flags are mutually incompatible.
<p>
-The path is calculated by choosing the steeper "slope" between adjacent
-cells. The slope calculation accurately acounts for the variable scale in
-lat-lon projections.
+For an elevation surface, the path is calculated by choosing the steeper
+"slope" between adjacent cells. The slope calculation accurately acounts
+for the variable scale in lat-lon projections. For a cost surface, the path
+is calculated by following the movement direction surface back to the start
+point given in <em>r.walk</em> or <em>r.cost</em>.
<p>
The <b>coordinate</b> parameter consists of map E and N grid coordinates of
@@ -153,8 +157,20 @@
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
</pre></div>
+With the <b>-d</b> <em>(direction)</em> flag, the direction raster is used
+for the input, rather than the elevation surface. The output is then created
+according to one of the <b>-can</b> flags.
+<div class="code"><pre>
+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 <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>
-
<h2>BUGS</h2>
<p>
Modified: grass/branches/develbranch_6/raster/r.drain/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.drain/main.c 2009-03-13 06:14:05 UTC (rev 36350)
+++ grass/branches/develbranch_6/raster/r.drain/main.c 2009-03-13 06:14:56 UTC (rev 36351)
@@ -10,6 +10,9 @@
* 24 July 2004: WebValley 2004, error checking and vector points added by
* Matteo Franchi Liceo Leonardo Da Vinci Trento
* Roberto Flor ITC-irst
+ * New code added by Colin Nielsen <colin.nielsen at gmail dot com> *
+ * to use movement direction surface from r.walk and r.cost, and to
+ * output vector paths 2/2009
*
* PURPOSE: This is the main program for tracing out the path that a
* drop of water would take if released at a certain location
@@ -39,6 +42,7 @@
#include <grass/gis.h>
#include <grass/site.h>
#include <grass/glocale.h>
+#include <grass/Vect.h>
#define DEBUG
#include "tinf.h"
@@ -59,23 +63,25 @@
int main(int argc, char **argv)
{
- int fe, fd;
+ int fe, fd, dir_fd;
int i, have_points = 0;
int new_id;
int nrows, ncols, points_row[MAX_POINTS], points_col[MAX_POINTS], npoints;
int cell_open(), cell_open_new();
- int map_id;
- char map_name[GNAME_MAX], *map_mapset, new_map_name[GNAME_MAX];
- char *tempfile1, *tempfile2;
+ int map_id, dir_id;
+ char map_name[GNAME_MAX], *map_mapset, new_map_name[GNAME_MAX],
+ dir_name[GNAME_MAX], *dir_mapset;
+ char *tempfile1, *tempfile2, *tempfile3;
char *search_mapset;
struct History history;
struct Cell_head window;
- struct Option *opt1, *opt2, *coordopt, *vpointopt;
- struct Flag *flag1, *flag2, *flag3;
+ struct Option *opt1, *opt2, *coordopt, *vpointopt, *opt3, *opt4;
+ struct Flag *flag1, *flag2, *flag3, *flag4;
struct GModule *module;
- int in_type;
+ int in_type, dir_data_type;
void *in_buf;
+ void *dir_buf;
CELL *out_buf;
struct band3 bnd, bndC;
struct metrics *m = NULL;
@@ -83,10 +89,18 @@
struct point *list;
struct point *thispoint;
int ival, bsz, start_row, start_col, mode;
+ int costmode = 0;
double east, north, val;
struct point *drain(int, struct point *, int, int);
+ struct point *drain_cost(int, struct point *, int, int);
int bsort(int, struct point *);
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ struct Map_info vout;
+ int cat;
+ double x, y;
+ char vect[GNAME_MAX], *vect_mapset;
G_gisinit(argv[0]);
@@ -97,9 +111,25 @@
opt1 = G_define_standard_option(G_OPT_R_ELEV);
opt1->key = "input";
+
+ opt3 = G_define_option();
+ opt3->key = "indir";
+ opt3->type = TYPE_STRING;
+ opt3->gisprompt = "old,cell,raster";
+ opt3->description =
+ _("Name of movement direction map associated with the cost surface");
+ opt3->required = NO;
opt2 = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt4 = G_define_option();
+ opt4->key = "voutput";
+ opt4->type = TYPE_STRING;
+ opt4->gisprompt = "new,vector,vector";
+ opt4->required = NO;
+ opt4->description =
+ _("Output drain vector map (recommended for cost surface made using knight's move)");
+
coordopt = G_define_option();
coordopt->key = "coordinate";
coordopt->type = TYPE_STRING;
@@ -127,6 +157,11 @@
flag3->key = 'n';
flag3->description = _("Count cell numbers along the path");
+ flag4 = G_define_flag();
+ flag4->key = 'd';
+ flag4->description =
+ _("The input surface is a cost surface (if checked, a direction surface must also be specified");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -134,6 +169,50 @@
strcpy(map_name, opt1->answer);
strcpy(new_map_name, opt2->answer);
+ if (flag4->answer) {
+ costmode = 1;
+ G_message(_
+ ("Directional drain selected... checking for direction raster"));
+ }
+ else {
+ G_message(_("Surface/Hydrology drain selected"));
+ }
+
+ if (costmode == 1) {
+ if (!opt3->answer) {
+ G_fatal_error(_
+ ("Direction raster not specified, if direction flag is on, a direction raster must be given"));
+ }
+ strcpy(dir_name, opt3->answer);
+ dir_mapset = G_find_cell2(dir_name, "");
+ if (dir_mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), dir_name);
+ else {
+ G_message(_("Direction raster found <%s>"), dir_name);
+ }
+ dir_data_type = G_raster_map_type(dir_name, dir_mapset);
+ }
+ if (costmode == 0) {
+ if (opt3->answer) {
+ G_fatal_error(_
+ ("Direction map <%s> should not be specified for Surface/Hydrology drains"),
+ opt3->answer);
+ }
+ }
+
+ if (opt4->answer) {
+ G_message(_("Outputting a vector path"));
+ /*vect = opt4->answer; */
+ if (G_legal_filename(opt4->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), opt4->answer);
+ /*G_ask_vector_new("",vect); */
+ if (0 > Vect_open_new(&vout, opt4->answer, 0)) {
+ G_fatal_error(_("Unable to create vector map <%s>"),
+ opt4->answer);
+ }
+ Vect_hist_command(&vout);
+ }
+
/* get the name of the elevation map layer for filling */
map_mapset = G_find_cell(map_name, "");
if (!map_mapset)
@@ -160,6 +239,10 @@
G_get_window(&window);
nrows = G_window_rows();
ncols = G_window_cols();
+ if (vect) {
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+ }
/* calculate true cell resolution */
m = (struct metrics *)G_malloc(nrows * sizeof(struct metrics));
@@ -297,14 +380,30 @@
}
G_close_cell(map_id);
- G_verbose_message(_("Calculating flow directions..."));
+ if (costmode == 1) {
+ dir_buf = G_allocate_d_raster_buf();
+ dir_id = G_open_cell_old(dir_name, dir_mapset);
+ tempfile3 = G_tempfile();
+ dir_fd = open(tempfile3, O_RDWR | O_CREAT);
- /* fill one-cell pits and take a first stab at flow directions */
- filldir(fe, fd, nrows, &bnd, m);
+ for (i = 0; i < nrows; i++) {
+ G_get_d_raster_row(dir_id, dir_buf, i);
+ write(dir_fd, dir_buf, ncols * sizeof(DCELL));
+ }
+ G_close_cell(dir_id);
+ }
- /* determine flow directions for more ambiguous cases */
- resolve(fd, nrows, &bndC);
+ /* only necessary for non-dir drain */
+ if (costmode == 0) {
+ G_verbose_message(_("Calculating flow directions..."));
+ /* fill one-cell pits and take a first stab at flow directions */
+ filldir(fe, fd, nrows, &bnd, m);
+
+ /* determine flow directions for more ambiguous cases */
+ resolve(fd, nrows, &bndC);
+ }
+
/* free the buffers already used */
G_free(bndC.b[0]);
G_free(bndC.b[1]);
@@ -323,7 +422,13 @@
thispoint->row = points_row[i];
thispoint->col = points_col[i];
thispoint->next = NULL;
- thispoint = drain(fd, thispoint, nrows, ncols);
+ /* drain algorithm selection (dir or non-dir) */
+ if (costmode == 0) {
+ thispoint = drain(fd, thispoint, nrows, ncols);
+ }
+ if (costmode == 1) {
+ thispoint = drain_cost(dir_fd, thispoint, nrows, ncols);
+ }
}
/* do the output */
@@ -429,6 +534,38 @@
G_percent(1, 1, 1);
}
+ /* Output a vector path */
+ if (opt4->answer) {
+ /* Need to modify for multiple paths */
+ thispoint = list;
+ i = 1;
+ while (thispoint->next != NULL) {
+ if (thispoint->row == INT_MAX) {
+ thispoint = thispoint->next;
+ Vect_cat_set(Cats, 1, i);
+ Vect_write_line(&vout, GV_LINE, Points, Cats);
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+ i++;
+ continue;
+ }
+ if (Vect_cat_get(Cats, 1, &cat) == 0) {
+ Vect_cat_set(Cats, 1, i);
+ }
+ /* Need to convert row and col to coordinates
+ * y = cell_head.north - ((double) p->row + 0.5) * cell_head.ns_res;
+ * x = cell_head.west + ((double) p->col + 0.5) * cell_head.ew_res;
+ */
+
+ x = window.west + ((double)thispoint->col + 0.5) * window.ew_res;
+ y = window.north - ((double)thispoint->row + 0.5) * window.ns_res;
+ Vect_append_point(Points, x, y, 0.0);
+ thispoint = thispoint->next;
+ } /* End while */
+ Vect_build(&vout);
+ Vect_close(&vout);
+ }
+
/* close files and free buffers */
G_close_cell(new_id);
@@ -446,6 +583,11 @@
G_free(in_buf);
G_free(out_buf);
+ if (costmode == 1) {
+ close(dir_fd);
+ unlink(tempfile3);
+ G_free(dir_buf);
+ }
G_done_msg(" ");
exit(EXIT_SUCCESS);
@@ -509,3 +651,134 @@
return list;
}
+
+struct point *drain_cost(int dir_fd, struct point *list, int nrow, int ncol)
+{
+ /*
+ * The idea is that each cell of the direction surface has a value representing
+ * the direction towards the next cell in the path. The direction is read from
+ * the input raster, and a simple case/switch is used to determine which cell to
+ * read next. This is repeated via a while loop until a null direction is found.
+ */
+
+ int neighbour, row, col, next_row, next_col, go = 1;
+ DCELL direction;
+ DCELL *dir_buf;
+
+ dir_buf = G_allocate_d_raster_buf();
+
+ next_row = list->row;
+ next_col = list->col;
+
+ while (go) {
+ go = 0;
+ /* Directional algorithm
+ * 1) read cell direction
+ * 2) shift to cell in that direction
+ */
+ /* find the direction recorded at row,col */
+ lseek(dir_fd, (off_t) list->row * ncol * sizeof(DCELL), SEEK_SET);
+ read(dir_fd, dir_buf, ncol * sizeof(DCELL));
+ direction = *(dir_buf + list->col);
+ neighbour = direction * 10;
+ if (G_verbose() > 2)
+ G_message(_("direction read: %lf, neighbour found: %i"),
+ direction, neighbour);
+ switch (neighbour) {
+ case 1800:
+ next_row = list->row;
+ next_col = list->col + 1;
+ break;
+ case 0:
+ next_row = list->row;
+ next_col = list->col - 1;
+ break;
+ case 900:
+ next_row = list->row + 1;
+ next_col = list->col;
+ break;
+ case 2700:
+ next_row = list->row - 1;
+ next_col = list->col;
+ break;
+ case 1350:
+ next_col = list->col + 1;
+ next_row = list->row + 1;
+ break;
+ case 450:
+ next_col = list->col - 1;
+ next_row = list->row + 1;
+ break;
+ case 3150:
+ next_row = list->row - 1;
+ next_col = list->col - 1;
+ break;
+ case 2250:
+ next_row = list->row - 1;
+ next_col = list->col + 1;
+ break;
+ case 1125:
+ next_row = list->row + 2;
+ next_col = list->col + 1;
+ break;
+ case 675:
+ next_row = list->row + 2;
+ next_col = list->col - 1;
+ break;
+ case 2925:
+ next_row = list->row - 2;
+ next_col = list->col - 1;
+ break;
+ case 2475:
+ next_row = list->row - 2;
+ next_col = list->col + 1;
+ break;
+ case 1575:
+ next_row = list->row + 1;
+ next_col = list->col + 2;
+ break;
+ case 225:
+ next_row = list->row + 1;
+ next_col = list->col - 2;
+ break;
+ case 3375:
+ next_row = list->row - 1;
+ next_col = list->col - 2;
+ break;
+ case 2025:
+ next_row = list->row - 1;
+ next_col = list->col + 2;
+ break;
+ /* default:
+ break;
+ Should probably add something here for the possibility of a non-direction map
+ G_fatal_error(_("Invalid direction given (possibly not a direction map)")); */
+ } /* end switch/case */
+
+ if (next_col >= 0 && next_col < ncol && next_row >= 0 &&
+ next_row < nrow) {
+ list->next = (struct point *)G_malloc(sizeof(struct point));
+ list = list->next;
+ list->row = next_row;
+ list->col = next_col;
+ next_row = -1;
+ next_col = -1;
+ go = 1;
+ }
+
+ } /* end while */
+
+ /* allocate and fill the end-of-path flag */
+ list->next = (struct point *)G_malloc(sizeof(struct point));
+ list = list->next;
+ list->row = INT_MAX;
+
+ /* return a pointer to an empty structure */
+ list->next = (struct point *)G_malloc(sizeof(struct point));
+ list = list->next;
+ list->next = NULL;
+
+ G_free(dir_buf);
+
+ return list;
+}
More information about the grass-commit
mailing list