[GRASS-SVN] r40255 - grass-addons/raster/r.stream.basins
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 5 09:24:50 EST 2010
Author: jarekj71
Date: 2010-01-05 09:24:49 -0500 (Tue, 05 Jan 2010)
New Revision: 40255
Modified:
grass-addons/raster/r.stream.basins/Makefile
grass-addons/raster/r.stream.basins/description.html
grass-addons/raster/r.stream.basins/global.h
grass-addons/raster/r.stream.basins/io.c
Log:
update
Modified: grass-addons/raster/r.stream.basins/Makefile
===================================================================
--- grass-addons/raster/r.stream.basins/Makefile 2010-01-05 14:10:20 UTC (rev 40254)
+++ grass-addons/raster/r.stream.basins/Makefile 2010-01-05 14:24:49 UTC (rev 40255)
@@ -4,7 +4,10 @@
PGM = r.stream.basins
-LIBES = $(GISLIB)
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
DEPENDENCIES = $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/raster/r.stream.basins/description.html
===================================================================
--- grass-addons/raster/r.stream.basins/description.html 2010-01-05 14:10:20 UTC (rev 40254)
+++ grass-addons/raster/r.stream.basins/description.html 2010-01-05 14:24:49 UTC (rev 40255)
@@ -6,11 +6,11 @@
<DT><b>-c</b></DT>
<DD>By default r.stream.basins uses streams category as basin category. In some cases - for example if
stream map is product of map algebra and separete streams may not have unique values this option will create
-new category sequence for each basin
+new category sequence for each basin (do not work in vector point mode)
</DD>
<p>
<DT><b>-l</b></DT>
-<DD>By default r.strean.basins create basins for all unique streams. This option delinate basins only for last streams ignoring upstreams.
+<DD>By default r.stream.basins create basins for all unique streams. This option delinate basins only for last streams ignoring upstreams (do not work in vector point mode).
</DD>
<p>
@@ -18,22 +18,48 @@
<DD>Stream network: name of input stream map on which ordering will be performed produced by r.watershed or r.stream.extract. Because streams network produced by r.watershed and r.stream.extract may slighty differ in detail it is required to use both stream and direction map produced by the same module. Stream background shall have NULL value or zero value. Background values of NULL are by default produced by r.watershed and r.stream.extract. If not 0 or NULL use <a href="r.mapcalc.html">r.mapcalc</a> to set background values to null.
</DD>
<p>
+
+<DT><b>cats</b></DT>
+<DD>Stream categories to delineate basins for: All categories which are not in stream map are ignored. It can be used with stream network created by r.watershed, r.stream.extract or r.stream.order. For r.stream.order use category of order for which basins must be created. For example to delineate only basins for order two use cats=2. If you need unique category for every basin use -c flag.
+</DD>
+<p>
+
<DT><b>dir</b></DT>
<DD>Flow direction: name of input direction map produced by r.watershed or r.stream.extract. If r.stream.extract output map is used, it only has non-NULL values in places where streams occur. NULL (nodata) cells are ignored, zero and negative values are valid direction data if they vary from -8 to 8 (CCW from East in steps of 45 degrees). Direction map shall be of type CELL values. Region resolution and map resoultion must be the same.
Also <em>stream</em> network map must have the same resolution. It is checked by default. If resolutions differ the module informs about it and stops. Region boundary and maps boundary may be differ but it may lead to unexpected results.</DD>
+<p>
+<DT><b>coors</b></DT>
+<DD>East and north coordinates for basin outlet. It can delinate only one basin using that option. This option simply copies funcionality of <a href="r.water.outlet.html">r.water.outlet</a>.
+</DD>
+<p>
+
+<DT><b>points</b></DT>
+<DD>Vector file containing basins outlet as vector points. Only point's categories are used to prepare basins. Table attached to it is ignored. Every point shall heve his own unique category. In that mode flags -l and -c are ignored
+</DD>
+<p>
+
+
<h2>OUTPUTS</h2>
<P>The module produces one raster map with basins acording user's rules</p>
</DL>
<h2>DESCRIPTION</h2>
<P>
-Module r.stream.basins is prepared to delineate basins and subasins according user rules. Module is prepared to delineate unrestricted number of basins in one step. It requires two maps: direction and streams. In stream map we can store all information required to proper delineation, so it may be necessary to modify it before run. Module is prepared to work with output data of <em>r.watershed, r.stream.extract, r.stream.order</em> also with modification done by <em>r.recalss</em> and <em>r.mapcalc</em>. r.stream.basin can delineate basins according outlets marked by raster streams, points and polygons. If outlets are marked by points it delineate basins which cells contribute to that points, if outlets are marked by streams it delineate cells which contribute to the last (downstream) cell of the every stream. If outlets are marked by polygon it delineate cells contributing to most downstream cell of the polygon. If polygon covers more outlets than of one basins it will create collective basin for all outlets with common category.
+Module r.stream.basins is prepared to delineate basins and subasins according user rules. Module is prepared to delineate unrestricted number of basins in one step. It can delineate basins with three methods:
+<UL>
+<LI>Using coordinates: his option simply copies funcionality of <a href="r.water.outlet.html">r.water.outlet</a>.
+<LI>Using vector points: it allow to mannually point outlets with any method
+<LI>Using streams (most advanced) it allow on lots of modifications. See examples for more details.
+</UL>
+Only one methdo can be used at once. Methods cannot be mixed.
+<P>
+The most recommended method require two maps: direction and streams. In spite of in stream map we can store information required to proper delineation, we can also enumarate stream categories for which basins are to be created (cats option). Module is prepared to work with output data of <em>r.watershed, r.stream.extract, r.stream.order</em> also with modification done by <em>r.recalss</em> and <em>r.mapcalc</em>. r.stream.basin can delineate basins according outlets marked by raster streams, and polygons, vector points and numerical coordinates. If outlets are marked by points or coordinates it delineate basins which cells contribute to that points, if outlets are marked by streams it delineate cells which contribute to the last (downstream) cell of the every stream. If outlets are marked by polygon it delineate cells contributing to most downstream cell of the polygon. If polygon covers more outlets than of one basins it will create collective basin for all outlets with common category.
<h2>NOTES</h2>
<P>
-Nowadays user's input must be in raster format. Vector format is planned in the feature. To receive good results outlets markers created by user shall overlapping with streams. On the other way basins could results with very small area. Input maps must be in CELL format (default output of r.watershed, r.stream.order and r.stream.extract)<P>
+To receive good results outlets markers created by user shall overlapping with streams. On the other way basins could results with very small area. Input maps must be in CELL format (default output of r.watershed, r.stream.order and r.stream.extract)<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).
<h2>EXAMPLES</h2>
@@ -45,10 +71,23 @@
To determine major and minor basins in area, definied by outlets, ignoring subbasins use - l flag. That flag ignores all nodes and uses only real outlets (in most cases that on map border):
<P>
<CODE>r.stream.basins -l dir=dirs stream=streams basins=bas_basins_last</CODE>
+
<P>
+<CODE>r.stream.basins dir=dirs coors=639936.623832,216939.836449</CODE>
+
+<P>
To delinaeate one or more particular basins defined by given streams, stream map must be re-classed first, to leave only desired streams:
+
<CODE>
<PRE>
+r.reclass input=streams cats=42 output=sel_streams_1
+r.reclass input=streams cats=42,252,188 output=sel_streams_1
+</PRE>
+</CODE>
+<P>
+Or alternatevely:
+<CODE>
+<PRE>
echo '42=42
* = NULL' > tmp #for one output
@@ -61,6 +100,27 @@
</PRE>
</CODE>
<P>
+Do delineate basins of particular order we must use the following procedure:
+<PRE>
+<CODE>
+r.stream.basins -lc dir=dirs stream=strahler cats=2 basins=bas_basin_strahler_2
+</CODE>
+</CODE>
+<P>
+Or alternatevely:
+<PRE>
+<CODE>
+echo '2 = 2
+* = NULL' > tmp
+r.reclass input=ord_strahler output=sel_strahler_2 < tmp
+r.stream.basins -c dir=dirs
+stream=sel_strahler_2 basins=bas_basin_strahler_2
+</CODE>
+</CODE>
+
+
+
+<P>
The usage of polygons as outlets markers is very useful when exact stream course cannot be cleary determined before running analysis, but the area of its occurrence can be determined (mostly in iterative simulations) Example uses r.circle but can be substituted by any polygon created for example with v.digit:
<CODE>
<PRE>
@@ -79,12 +139,11 @@
<PRE>
v.extract -d input=lakes at PERMANENT output=lake8056 type=area layer=1 'where=FULL_HYDRO = 8056' new=-1
v.to.rast input=lake8056 output=lake8056 use=val type=area layer=1 value=1
-echo 'sel_streams_lake=if(lake8056,streams,null())'|r.mapcalc
-r.stream.basins dir=dirs stream=sel_streams_lake basins=bas_basin_lake
+r.stream.basins dir=dirs stream=lake8056 basins=bas_basin_lake
</PRE>
</CODE>
<P>
-See aslo tutorial: <a href="http://heretic.livenet.pl/heretic">http://heretic.livenet.pl/heretic</a>
+See aslo tutorial: <a href="http://grass.osgeo.org/wiki/R.stream.*">http://grass.osgeo.org/wiki/R.stream.*</a>
<h2>SEE ALSO</h2>
Modified: grass-addons/raster/r.stream.basins/global.h
===================================================================
--- grass-addons/raster/r.stream.basins/global.h 2010-01-05 14:10:20 UTC (rev 40254)
+++ grass-addons/raster/r.stream.basins/global.h 2010-01-05 14:24:49 UTC (rev 40255)
@@ -3,6 +3,7 @@
#include <string.h>
#include <math.h>
#include <grass/gis.h>
+#include <grass/Vect.h>
/* define */
@@ -26,7 +27,7 @@
int r, c;
int val;
};
-
+
/* functions.c */
/* io.c */
@@ -35,6 +36,9 @@
int max_link(void);
int write_chatchment(void);
int set_null(void);
+int process_coors (char **answers);
+int process_cats (char **answers);
+int process_vector (char *in_point);
/* cachments */
int find_outlets(void);
@@ -43,7 +47,6 @@
int fifo_insert (POINT point);
POINT fifo_return_del (void);
-
/* variables */
#ifdef MAIN
@@ -53,7 +56,7 @@
#endif
GLOBAL struct Cell_head window;
-GLOBAL char *in_dirs, *in_streams; /* input dirrection and accumulation raster names*/
+GLOBAL char *in_dirs, *in_streams, *in_point; /* input dirrection and accumulation raster names*/
GLOBAL char *name_catchments;
GLOBAL int zeros, cats, lasts; /* flags */
@@ -67,8 +70,11 @@
int fifo_max;
GLOBAL int out; /* number of strahler and horton outlets: index */
+GLOBAL int *categories;
+
OUTLET *outlets;
+
GLOBAL struct History history; /* holds meta-data (title, comments,..) */
Modified: grass-addons/raster/r.stream.basins/io.c
===================================================================
--- grass-addons/raster/r.stream.basins/io.c 2010-01-05 14:10:20 UTC (rev 40254)
+++ grass-addons/raster/r.stream.basins/io.c 2010-01-05 14:24:49 UTC (rev 40255)
@@ -9,20 +9,20 @@
mapset = G_find_cell2(mapname, ""); /* checking if map exist */
if (mapset == NULL) {
- G_fatal_error(_("Raster map <%s> not found"), mapname);
+ G_fatal_error("Raster map <%s> not found", mapname);
}
if (G_raster_map_type(mapname, mapset) != CELL_TYPE) /* checking if the input maps type is CELL */
- G_fatal_error(_("<%s> is not of type CELL."), mapname);
+ G_fatal_error("<%s> is not of type CELL.", mapname);
if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
- G_fatal_error(_("Unable to open raster map <%s>"), mapname);
+ G_fatal_error("Unable to open raster map <%s>", mapname);
if (G_get_cellhd(mapname, mapset, &cellhd) < 0)
- G_fatal_error(_("Unable to read file header of <%s>"), mapname);
+ G_fatal_error("Unable to read file header of <%s>", mapname);
if (window.ew_res != cellhd.ew_res || window.ns_res != cellhd.ns_res)
- G_fatal_error(_("Region resolution and map %s resolution differs. Run g.region rast=%s to set proper resolution"),
+ G_fatal_error("Region resolution and map %s resolution differs. Run g.region rast=%s to set proper resolution",
mapname, mapname);
return fd;
@@ -35,62 +35,66 @@
CELL *r_dirs, *r_streams;
+ G_message("Reading maps...");
+
+
+
in_dir_fd = open_raster(in_dirs);
- in_stm_fd = open_raster(in_streams);
-
r_dirs = (CELL *) G_malloc(sizeof(CELL) * ncols);
- r_streams = (CELL *) G_malloc(sizeof(CELL) * ncols);
-
-
dirs = (CELL **) G_malloc(sizeof(CELL *) * nrows);
- streams = (CELL **) G_malloc(sizeof(CELL *) * nrows);
- G_message(_("Reading maps..."));
+ r_streams = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ streams = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+
+ if (in_streams) {
+
+ in_stm_fd = open_raster(in_streams);
+
+ for (r = 0; r < nrows; ++r) {
+ G_percent(r, nrows, 2);
+ streams[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
- for (r = 0; r < nrows; ++r) {
- G_percent(r, nrows, 2);
- /* dirs & streams */
- dirs[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
- streams[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
-
-
- if (G_get_c_raster_row(in_dir_fd, r_dirs, r) < 0 ||
- G_get_c_raster_row(in_stm_fd, r_streams, r) < 0) {
- G_close_cell(in_dir_fd);
+ if (G_get_c_raster_row(in_stm_fd, r_streams, r) < 0) {
G_close_cell(in_stm_fd);
- G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
- }
+ G_fatal_error("Unable to read raster maps at row <%d>", r);
+ }
- for (c = 0; c < ncols; ++c) {
- if (G_is_c_null_value(&r_dirs[c])) {
- dirs[r][c] = 0;
- }
- else {
- dirs[r][c] = r_dirs[c];
-
- }
-
- if (G_is_c_null_value(&r_streams[c])) {
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_c_null_value(&r_streams[c]))
streams[r][c] = 0;
- }
- else {
+ else
streams[r][c] = r_streams[c];
- }
+ }
+ }
+ G_free(r_streams);
+ G_close_cell(in_stm_fd);
+ G_percent(r, nrows, 2);
+ } /* end if streams */
- }
+ for (r = 0; r < nrows; ++r) {
+ G_percent(r, nrows, 2);
- /* END dirs & streams & accums */
+ dirs[r] = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ if (!in_streams)
+ streams[r] =(CELL *) G_malloc(sizeof(CELL) * ncols);
+ if (G_get_c_raster_row(in_dir_fd, r_dirs, r) < 0) {
+ G_close_cell(in_dir_fd);
+ G_fatal_error("Unable to read raster maps at row <%d>", r);
+ }
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_c_null_value(&r_dirs[c]))
+ dirs[r][c] = 0;
+ else
+ dirs[r][c] = r_dirs[c];
+ }
} /*end for r */
G_free(r_dirs);
- G_free(r_streams);
+ G_close_cell(in_dir_fd);
G_percent(r, nrows, 2);
- G_close_cell(in_dir_fd);
- G_close_cell(in_stm_fd);
-
return 0;
} /* end create maps */
@@ -114,17 +118,16 @@
int fdo = 0;
if ((fdo = G_open_raster_new(name_catchments, CELL_TYPE)) < 0)
- G_fatal_error(_("Unable to create raster map <%s>"), name_catchments);
+ G_fatal_error("Unable to create raster map <%s>", name_catchments);
for (r = 0; r < nrows; ++r)
G_put_c_raster_row(fdo, streams[r]);
G_close_cell(fdo);
G_short_history(name_catchments, "raster", &history);
G_write_history(name_catchments, &history);
- G_message(_("%s Done!"), name_catchments);
+ G_message("%s Done!", name_catchments);
G_free(streams);
-
return 0;
}
@@ -140,3 +143,105 @@
}
return 0;
}
+
+int process_coors (char **answers) {
+
+ int n;
+ double X,Y;
+ outlets = G_malloc(sizeof(OUTLET));
+
+ if(!answers)
+ return 0;
+
+ for (n=0; answers[n] != NULL; n+=2) {
+ if(!G_scan_easting(answers[n], &X , G_projection()))
+ G_fatal_error("Wrong coordinate <%s>",answers[n]);
+ if(!G_scan_northing(answers[n+1], &Y , G_projection()))
+ G_fatal_error("Wrong coordinate <%s>",answers[n+1]);
+ }
+
+ if (X < window.west || X > window.east || Y < window.south || Y > window.north)
+ G_fatal_error("Coordinates outside window");
+
+ outlets->r = (window.north - Y) / window.ns_res;
+ outlets->c = (X - window.west) / window.ew_res;
+ outlets->val=1;
+
+ return 1; /* not an success */
+}
+
+int process_cats (char **answers) {
+
+ int i;
+ int link_max=max_link();
+ int cat;
+
+ if(!answers)
+ return;
+
+ categories = G_malloc(link_max* sizeof(int));
+
+ for (i=0;i<(link_max+1);++i)
+ categories[i]=-1;
+
+
+ for (i = 0; answers[i] != NULL; ++i) {
+ cat=atoi(answers[i]);
+ if (cat<1 || cat>link_max)
+ G_fatal_error("Stream categories must be > 0 and < maximum stream category");
+ categories[cat]=cat;
+ }
+
+ return 0; /* success */
+}
+
+
+int process_vector (char *in_point) {
+ char *mapset;
+ struct Map_info Map;
+ struct bound_box box;
+ int num_point=0;
+ int type, i, cat;
+
+ struct line_pnts* sites;
+ struct line_cats* cats;
+
+ sites = Vect_new_line_struct();
+ cats = Vect_new_cats_struct();
+
+
+ mapset = G_find_vector2(in_point, "");
+ if (mapset == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), in_point);
+
+ if (Vect_open_old(&Map, in_point, mapset) < 0)
+ G_fatal_error("Cannot open vector map <%s>", in_point);
+
+ Vect_region_box(&window, &box);
+
+ while (type=Vect_read_next_line(&Map, sites, cats) > -1) {
+ if (type != GV_POINT)
+ continue;
+ if (Vect_point_in_box(sites->x[0], sites->y[0], sites->z[0], &box))
+ num_point++;
+ }
+
+ outlets = (OUTLET *)G_malloc(num_point * sizeof(OUTLET));
+
+ for (i=0;i<num_point;++i) {
+
+ type = Vect_read_line(&Map, sites, cats,i+1);
+ if (type != GV_POINT)
+ continue;
+
+ if (!Vect_point_in_box(sites->x[0], sites->y[0], sites->z[0], &box))
+ continue;
+
+ Vect_cat_get(cats, 1, &cat);
+
+ outlets[i].r=(int)G_northing_to_row(sites->y[0], &window);
+ outlets[i].c=(int)G_easting_to_col(sites->x[0], &window);
+ outlets[i].val=cat;
+ }
+ return num_point;
+}
More information about the grass-commit
mailing list