[GRASS-SVN] r41182 - grass-addons/raster/r.stream.distance
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 24 12:32:16 EST 2010
Author: jarekj71
Date: 2010-02-24 12:32:15 -0500 (Wed, 24 Feb 2010)
New Revision: 41182
Modified:
grass-addons/raster/r.stream.distance/description.html
grass-addons/raster/r.stream.distance/distance.c
grass-addons/raster/r.stream.distance/global.h
grass-addons/raster/r.stream.distance/main.c
Log:
add new method <upstream> method of calculation
Modified: grass-addons/raster/r.stream.distance/description.html
===================================================================
--- grass-addons/raster/r.stream.distance/description.html 2010-02-24 14:24:45 UTC (rev 41181)
+++ grass-addons/raster/r.stream.distance/description.html 2010-02-24 17:32:15 UTC (rev 41182)
@@ -1,11 +1,13 @@
<h2>OPTIONS</h2>
<DL>
<DT><b>-o</b></DT>
-<DD>Calculate distance and elevation against basin outlets instead of streams. It choose only last outlets in
+<DD>Downstream method only. Calculate distance and elevation against basin outlets instead of streams. It choose only last outlets in
the network ignoring nodes.</DD>
<p>
<DT><b>-s</b></DT>
-<DD>Calculate distance and elevation against stream nodes instead of streams. It create distance and elevation parameters not for whole basins but for subbasins
+<DD>Downstream method only. Calculate distance and elevation against stream nodes instead of streams. It create distance and elevation parameters not for whole basins but for subbasins
+<DT><b>-n</b></DT>
+<DD>For upstram method only. Calculate distance and elevation to the nearest local maximum/divide. With the default option distance/elevation is calculated to the farthest possible maximum/divide
</DD>
<p>
@@ -19,9 +21,10 @@
<p>
<DT><b>elev</b></DT>
<DD>Elevation: name of input elevation map. Map can be of type CELL, FCELL or DCELL. It is not restricted to resolution of region settings as stream and dir.</DD>
+<DT><b>method</b></DT>
+<DD>It is possible to calculate distance with two method: <b>downstream</b> from any raster cell to the nearest stream cell/ junction cell or outlet or <b>upstream</b> from any cell upstream to the local maximum</DD>
-
<h2>OUTPUTS</h2>
<DT><b>elevation</b></DT>
<DD>Returns elevation above the targer (outlet, node stream) along watercoures. The map is of FCELL type</DD>
@@ -31,9 +34,13 @@
<h2>DESCRIPTION</h2>
<P>
-Module r.stream.distance is prepared to calculate distance to streams and outlets and elevation above streams and outlets. The distance and elevation is calculated along watercourses. In outlets mode it can also calculate parameters for subbasins.
+Module r.stream.distance may calculate distance using two methods: downstream and upstream.
<P>
-In streams mode (default) it calculates that parameters upstream from streams which are added as stream mask. In outlets mode there are some additional possibilities. If subbasin is off it calculate parameters only for last point of last (downstream) CELL. In subbasin mode it calculates parameters for every subbasin separately. Subbasin mode acts similar to subbasin mask. Streams file prepared to create basins and subbasins with r.stream.basins can use t to calculate distance and elevation parameters.
+The default is downstream method when it calculate distance to streams and outlets and elevation above streams and outlets. The distance and elevation is calculated along watercourses. In outlets mode it can also calculate parameters for subbasins.
+<P>
+In streams mode (default) it calculates that parameters downstream to streams which are added as stream mask. In outlets mode there are some additional possibilities. If subbasin is off it calculate parameters only for last point of last (downstream) CELL. In subbasin mode it calculates parameters for every subbasin separately. Subbasin mode acts similar to subbasin mask. Streams file prepared to create basins and subbasins with r.stream.basins can use to to calculate distance and elevation parameters.
+<P>
+With upstream method it calculate distance to the local maximum or divide. Opposite to downstream method, where every cell has one and only one downstream cell in upstream method every cell has usssualy more than one upstream cell. So it is impossible to determine nterchangeable path from any cell. The upstream method offers two alternative modes switched with -n flag: nearest local maximum/divide: means the shortest path to local maximum and default option farthest maximum/divide means the longest path. In hydrological sense nearest mode means the shortest path which particle of water must run from divide to reach particular cell, while farthest mode means the possible longest path.
<h2>NOTES</h2>
<P>
@@ -42,7 +49,7 @@
<P>
<CODE>echo 'dist_corrected = sqrt(distance^2 + elevation ^2)'|r.mapcalc<CODE>
<P>
-Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted. this means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used. Nowadays f direction map comes from r.stream.extract must be patched by direction map from r.watershed. (with r.patch). .
+Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted. this means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used.
<h2>SEE ALSO</h2>
Modified: grass-addons/raster/r.stream.distance/distance.c
===================================================================
--- grass-addons/raster/r.stream.distance/distance.c 2010-02-24 14:24:45 UTC (rev 41181)
+++ grass-addons/raster/r.stream.distance/distance.c 2010-02-24 17:32:15 UTC (rev 41182)
@@ -1,101 +1,135 @@
#include <grass/glocale.h>
#include "global.h"
-int find_outlets(void)
-{
- int d, i, j; /* d: direction, i: iteration */
- int r, c;
- int next_stream = -1, cur_stream;
- int out_max = ncols + nrows;
+int calculate_upstream(void) {
+ int r,c;
+ int next_r, next_c;
+ float easting, northing;
+ float cell_easting, cell_northing;
+ int i,j,k;
+ int done;
+ int n_inits=0;
+ float cur_dist;
+ POINT* d_inits;
+ float* tmp_elevation=NULL; /* only for elevation */
+ float tmp_dist=0;
+ float tmp_elev=0;
+ float target_elev=0;
+ int tmp_inits=0;
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+
+ for (i=1;i<9;++i) {
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1) || c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
+ continue; /* out of border */
+
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+
+ if (dirs[r + nextr[i]][c + nextc[i]] == j &&
+ distance[r][c]!=0) {/* is contributing cell */
+ distance[r][c]=-1;
+ break;
+ }
+ }
+ if(distance[r][c]==1 && dirs[r][c]>0)
+ n_inits++;
+ else if (dirs[r][c]>0)
+ distance[r][c]=-1;
+ }
+ }
- int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
-
- G_message(_("Finding nodes..."));
- outlets = (OUTLET *) G_malloc((out_max) * sizeof(OUTLET));
-
- outlets_num = 0;
-
- for (r = 0; r < nrows; ++r) {
- for (c = 0; c < ncols; ++c) {
- if (streams[r][c] > 0) {
- if (outlets_num > (out_max - 1)) {
- out_max *= 2;
- outlets =
- (OUTLET *) G_realloc(outlets,
- out_max * sizeof(OUTLET));
+ d_inits=(POINT *) G_malloc(n_inits * sizeof(POINT));
+
+ if(out_elev) {
+ tmp_elevation=(float *) G_malloc(nrows*ncols*sizeof(float));
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ tmp_elevation[r*ncols+c]=elevation[r][c];
+ }
}
-
- d = abs(dirs[r][c]); /* abs */
- if (r + nextr[d] < 0 || r + nextr[d] > (nrows - 1) ||
- c + nextc[d] < 0 || c + nextc[d] > (ncols - 1)) {
- next_stream = -1; /* border */
+ }
+
+
+ k=0;
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+
+ if(distance[r][c]==1) {
+ distance[r][c]=0;
+ if(out_elev)
+ elevation[r][c]=0;
+
+ if(dirs[r+nextr[dirs[r][c]]][c+nextc[dirs[r][c]]]<0)
+ continue;
+ if(distance[next_r][next_c]==0)
+ continue;
+
+ d_inits[k].r=r;
+ d_inits[k].c=c;
+ d_inits[k].cur_dist=0;
+
+ if(out_elev)
+ d_inits[k].target_elev=tmp_elevation[r*ncols+c];
+
+ k++;
+ }
}
- else {
- next_stream = streams[r + nextr[d]][c + nextc[d]];
- if (next_stream < 1)
- next_stream = -1;
- }
- if (d == 0)
- next_stream = -1;
- cur_stream = streams[r][c];
+ }
+n_inits=k;
- if (subs && outs) { /* in stream mode subs is ignored */
- if (cur_stream != next_stream) { /* is outlet or node! */
- outlets[outlets_num].r = r;
- outlets[outlets_num].c = c;
- outlets[outlets_num].northing =
- window.north - (r + .5) * window.ns_res;
- outlets[outlets_num].easting =
- window.west + (c + .5) * window.ew_res;
- outlets_num++;
- }
- }
- else {
- if (cur_stream != next_stream && next_stream < 0) { /* is outlet! */
- outlets[outlets_num].r = r;
- outlets[outlets_num].c = c;
- outlets[outlets_num].northing =
- window.north - (r + .5) * window.ns_res;
- outlets[outlets_num].easting =
- window.west + (c + .5) * window.ew_res;
- outlets_num++;
- }
- } /* end lasts */
- } /* end if streams */
- } /* end for */
- } /* end for */
+while (n_inits>0) {
+ k=0;
- return 0;
-}
+ for(i=0;i<n_inits;++i) {
+ r=d_inits[i].r;
+ c=d_inits[i].c;
+ next_r=r+nextr[dirs[r][c]];
+ next_c=c+nextc[dirs[r][c]];
+ tmp_dist=d_inits[i].cur_dist;
+
+ if(out_elev)
+ target_elev=d_inits[i].target_elev;
+
+ easting=window.west + (c + 0.5) * window.ew_res;
+ northing=window.north - (r + 0.5) * window.ns_res;
+ cell_easting=window.west + (next_c + 0.5) * window.ew_res;
+ cell_northing=window.north - (next_r + 0.5) * window.ns_res;
+
+ cur_dist=tmp_dist+G_distance(easting,northing,cell_easting,cell_northing);
+
+
+ if(near)
+ done=(distance[next_r][next_c]>cur_dist||distance[next_r][next_c]==-1) ? 1: 0;
+ else
+ done=(distance[next_r][next_c]<cur_dist||distance[next_r][next_c]==-1) ? 1: 0;
+
+
+ if(done) {
+ distance[next_r][next_c]=cur_dist;
+ if(out_elev)
+ elevation[next_r][next_c]=target_elev-tmp_elevation[next_r*ncols+next_c];
-int reset_distance(void)
-{
- int r, c, i;
+ if(dirs[r+nextr[dirs[r][c]]][c+nextc[dirs[r][c]]]<1)
+ continue;
- distance = (FCELL **) G_malloc(sizeof(FCELL *) * nrows);
- if (!outs) { /* stream mode */
- for (r = 0; r < nrows; ++r) {
- distance[r] = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
- for (c = 0; c < ncols; ++c) {
- distance[r][c] = (streams[r][c]) ? 0 : -1;
- } /* r */
- } /* r */
- }
- else { /* outlets mode */
- for (r = 0; r < nrows; ++r) {
- distance[r] = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
- for (c = 0; c < ncols; ++c) {
- distance[r][c] = -1;
- }
- }
- for (i = 0; i < outlets_num; ++i) {
+ d_inits[k].r=next_r;
+ d_inits[k].c=next_c;
+ d_inits[k].cur_dist=cur_dist;
- distance[outlets[i].r][outlets[i].c] = 0;
+ if(out_elev)
+ d_inits[k].target_elev=target_elev;
+
+ k++;
+ } /* end of if done */
}
- }
- return 0;
+ n_inits=k;
}
+return 0;
+}
/*
Modified: grass-addons/raster/r.stream.distance/global.h
===================================================================
--- grass-addons/raster/r.stream.distance/global.h 2010-02-24 14:24:45 UTC (rev 41181)
+++ grass-addons/raster/r.stream.distance/global.h 2010-02-24 17:32:15 UTC (rev 41182)
@@ -16,20 +16,22 @@
*/
#define SQRT2 1.414214
-#define POINT struct points
+#define POINT struct points
+#define UPSTREAM 0
+#define DOWNSTREAM 1
+#define RELATIVE 2
+
POINT {
int r, c;
float cur_dist;
float target_elev;
-float northing;
-float easting;
};
#define OUTLET struct outs
OUTLET {
int r, c;
+ float easting;
float northing;
- float easting;
};
/* functions.c */
@@ -42,9 +44,14 @@
int write_elevation(void);
int set_null_elev(void);
-/* distance */
+/* inits */
+int find_chatchment_outlets(void);
+int fill_catchments(OUTLET outlet);
int find_outlets(void);
int reset_distance(void);
+
+/* distance */
+int calculate_upstream(void);
int fill_maps(OUTLET outlet);
int fifo_insert (POINT point);
POINT fifo_return_del (void);
@@ -60,9 +67,9 @@
GLOBAL char *in_dirs, *in_streams, *in_elev; /* input dirrection and accumulation raster names*/
GLOBAL char *out_dist, *out_elev;
-GLOBAL int zeros, outs, subs; /* flags */
+GLOBAL int zeros, outs, subs, near; /* flags */
+GLOBAL int method;
-
GLOBAL CELL **dirs, **streams; /* matrix with input data*/
GLOBAL FCELL **elevation, **distance;
/* streams and elevation are used to store internal data during processing */
Modified: grass-addons/raster/r.stream.distance/main.c
===================================================================
--- grass-addons/raster/r.stream.distance/main.c 2010-02-24 14:24:45 UTC (rev 41181)
+++ grass-addons/raster/r.stream.distance/main.c 2010-02-24 17:32:15 UTC (rev 41182)
@@ -29,14 +29,22 @@
*
*
*/
+
int main(int argc, char *argv[])
{
struct GModule *module; /* GRASS module for parsing arguments */
- struct Option *in_dir_opt, *in_stm_opt, *in_elev_opt, *out_dist_opt, *out_elev_opt; /* options */
- struct Flag *out_outs, *out_sub; /* flags */
+ struct Option *in_dir_opt,
+ *in_stm_opt,
+ *in_elev_opt,
+ *in_method_opt,
+ *out_dist_opt,
+ *out_elev_opt; /* options */
+ struct Flag *out_outs, *out_sub, *out_near; /* flags */
int link_max;
+ char* method_name[]={"upstream","downstream"};
+
/* initialize GIS environment */
G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
@@ -69,6 +77,14 @@
in_elev_opt->required = NO;
in_elev_opt->gisprompt = "old,cell,raster";
in_elev_opt->description = "Name of elevation map";
+
+ in_method_opt = G_define_option();
+ in_method_opt->key = "method";
+ in_method_opt->description ="Calculation method";
+ in_method_opt->type = TYPE_STRING;
+ in_method_opt->required = YES;
+ in_method_opt->options ="upstream,downstream";
+ in_method_opt->answer = "downstream";
/* output option - at least one is reqired */
@@ -88,17 +104,20 @@
out_elev_opt->gisprompt = "new,cell,raster";
out_elev_opt->description = "Output elevation map";
- /* Define the different flags */
out_outs = G_define_flag();
out_outs->key = 'o';
out_outs->description =
- _("Calculate parameters for outlets (outlet mode) instead of streams (stream mode: default)");
+ _("Calculate parameters for outlets (outlet mode) instead of (default) streams");
- /* Define the different flags */
out_sub = G_define_flag();
out_sub->key = 's';
out_sub->description =
_("Calculate parameters for subbasins (ignored in stream mode)");
+
+ out_near = G_define_flag();
+ out_near->key = 'n';
+ out_near->description =
+ _("Calculate nearest local maximum (ignored in downstream calculation)");
if (G_parser(argc, argv)) /* parser */
exit(EXIT_FAILURE);
@@ -116,6 +135,7 @@
out_elev = out_elev_opt->answer;
outs = (out_outs->answer != 0);
subs = (out_sub->answer != 0);
+ near = (out_near->answer != 0);
if (out_dist) {
if (G_legal_filename(out_dist) < 0)
@@ -126,6 +146,12 @@
G_fatal_error(_("<%s> is an illegal file name"), out_elev);
}
+
+ if (!strcmp(in_method_opt->answer, "upstream"))
+ method = UPSTREAM;
+ else if (!strcmp(in_method_opt->answer, "downstream"))
+ method = DOWNSTREAM;
+
nrows = G_window_rows();
ncols = G_window_cols();
G_get_window(&window);
@@ -135,20 +161,34 @@
find_outlets();
reset_distance();
free_streams();
-
+ G_message(_("Calculate %s distance"),method_name[method]);
+
+ if (method==UPSTREAM) {
+ int j;
+ fifo_max = 4 * (nrows + ncols);
+ fifo_outlet = (POINT *) G_malloc((fifo_max + 1) * sizeof(POINT));
+
+ for (j = 0; j < outlets_num; ++j) {
+ fill_catchments(outlets[j]);
+ }
+ G_free(fifo_outlet);
+ calculate_upstream();
+ }
+
+ if (method==DOWNSTREAM) {
+
/* fifo queue for contributing cell */
fifo_max = 4 * (nrows + ncols);
fifo_outlet = (POINT *) G_malloc((fifo_max + 1) * sizeof(POINT));
- { /* block */
- int j;
-
- G_message(_("Calculate..."));
- for (j = 0; j < outlets_num; ++j) {
+ {
+ int j;
+ for (j = 0; j < outlets_num; ++j) {
fill_maps(outlets[j]);
+ }
+ }
+ G_free(fifo_outlet);
}
- } /* end block */
- G_free(fifo_outlet);
if (out_elev) {
set_null(elevation);
@@ -159,6 +199,5 @@
write_maps(out_dist, distance);
}
-
exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list