[GRASS-SVN] r39247 - in grass-addons/raster: . r.stream.basins
r.stream.distance r.stream.order r.stream.stats
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Sep 19 07:15:36 EDT 2009
Author: martinl
Date: 2009-09-19 07:15:35 -0400 (Sat, 19 Sep 2009)
New Revision: 39247
Added:
grass-addons/raster/r.stream.basins/
grass-addons/raster/r.stream.basins/Makefile
grass-addons/raster/r.stream.basins/catchment.c
grass-addons/raster/r.stream.basins/description.html
grass-addons/raster/r.stream.basins/global.h
grass-addons/raster/r.stream.basins/io.c
grass-addons/raster/r.stream.basins/main.c
grass-addons/raster/r.stream.distance/
grass-addons/raster/r.stream.distance/Makefile
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/io.c
grass-addons/raster/r.stream.distance/main.c
grass-addons/raster/r.stream.order/
grass-addons/raster/r.stream.order/Makefile
grass-addons/raster/r.stream.order/description.html
grass-addons/raster/r.stream.order/global.h
grass-addons/raster/r.stream.order/io.c
grass-addons/raster/r.stream.order/main.c
grass-addons/raster/r.stream.order/order.c
grass-addons/raster/r.stream.order/orders.png
grass-addons/raster/r.stream.stats/
grass-addons/raster/r.stream.stats/Makefile
grass-addons/raster/r.stream.stats/description.html
grass-addons/raster/r.stream.stats/global.h
grass-addons/raster/r.stream.stats/io.c
grass-addons/raster/r.stream.stats/main.c
grass-addons/raster/r.stream.stats/print_stats.c
grass-addons/raster/r.stream.stats/stats.c
Log:
r.stream.* modules added (thanks to Jarek Jasiewicz)
see
http://lists.osgeo.org/pipermail/grass-dev/2009-September/046191.html
Added: grass-addons/raster/r.stream.basins/Makefile
===================================================================
--- grass-addons/raster/r.stream.basins/Makefile (rev 0)
+++ grass-addons/raster/r.stream.basins/Makefile 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,13 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.basins
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/raster/r.stream.basins/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.basins/catchment.c
===================================================================
--- grass-addons/raster/r.stream.basins/catchment.c (rev 0)
+++ grass-addons/raster/r.stream.basins/catchment.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,166 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ Link: a channel between junction
+ Ooutlet: is final cell of every segment
+ Segment: a channel which order remains unchanged in spite it pass through junction or not
+ Number of outlets shall be equal the number of segments
+ Number of junction shall be less than number of links
+ */
+
+/*
+ find outlets create table of outlets point, with r, c and value. depending of flag:
+ if flag -l is as anly last points of segment is added to table and uses the value is the category as value for whole basins
+ if flag -l = flase (default) last points of every link is added to table and nd uses the value is the category as value for subbasin
+ In both cases if flag -c is used it add out_num+1 value as value of point. That structure is next used in reset_catchments and fill_catchments fuctions
+ */
+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 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)) {
+ outlets =
+ (OUTLET *) G_realloc(outlets,
+ out_max * 2 * sizeof(OUTLET));
+ }
+
+ 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 */
+ }
+ 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];
+
+ if (lasts) {
+ if (cur_stream != next_stream && next_stream < 0) { /* is outlet! */
+ outlets[outlets_num].r = r;
+ outlets[outlets_num].c = c;
+ outlets[outlets_num].val =
+ (cats) ? outlets_num + 1 : streams[r][c];
+ outlets_num++;
+ }
+ }
+ else {
+ if (cur_stream != next_stream) { /* is outlet or node! */
+ outlets[outlets_num].r = r;
+ outlets[outlets_num].c = c;
+ outlets[outlets_num].val =
+ (cats) ? outlets_num + 1 : streams[r][c];
+ outlets_num++;
+ }
+ } /* end lasts */
+ } /* end if streams */
+ } /* end for */
+ } /* end for */
+
+ return 0;
+}
+
+/*
+ this function reset all values to 0 and adds points from outlet structure
+ The lines below from function fill_catchments makes that autlet is not added do fifo
+ queue and basin is limited by next outlet
+
+ if (catchments[r + nextr[i]][c + nextc[i]]>0)
+ continue;
+
+ It is simple trick but allow gives module its real funcionality
+ */
+int reset_catchments(void)
+{
+ int r, c, i;
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ streams[r][c] = 0;
+ }
+ }
+ for (i = 0; i < outlets_num; ++i) {
+ streams[outlets[i].r][outlets[i].c] = outlets[i].val;
+ }
+ return 0;
+}
+
+/*
+ algorithm uses fifo queue for determining basins area.
+ */
+int fill_catchments(OUTLET outlet)
+{
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ int r, c, val, i, j;
+ POINT n_cell;
+
+ tail = 0;
+ head = -1;
+ r = outlet.r;
+ c = outlet.c;
+ val = outlet.val;
+
+ streams[r][c] = val;
+
+ while (tail != head) {
+ 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; /* border */
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (dirs[r + nextr[i]][c + nextc[i]] == j) { /* countributing cell */
+
+ if (streams[r + nextr[i]][c + nextc[i]] > 0)
+ continue; /* other outlet */
+
+ streams[r + nextr[i]][c + nextc[i]] = val;
+ n_cell.r = (r + nextr[i]);
+ n_cell.c = (c + nextc[i]);
+ fifo_insert(n_cell);
+ }
+ } /* end for i... */
+
+ n_cell = fifo_return_del();
+ r = n_cell.r;
+ c = n_cell.c;
+ }
+
+ return 0;
+}
+
+/* fifo functions */
+int fifo_insert(POINT point)
+{
+ fifo_outlet[tail++] = point;
+ if (tail > fifo_max)
+ tail = 0;
+ return 0;
+}
+
+POINT fifo_return_del(void)
+{
+ if (head > fifo_max)
+ head = -1;
+ return fifo_outlet[++head];
+}
Property changes on: grass-addons/raster/r.stream.basins/catchment.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.basins/description.html
===================================================================
--- grass-addons/raster/r.stream.basins/description.html (rev 0)
+++ grass-addons/raster/r.stream.basins/description.html 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,108 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-z</b></DT>
+<DD>Creates zero-value background instead of NULL. For some reason (like map algebra calculation) zero-valued background may be required. This flag produces zero-filled background instead of null (default).</DD>
+<p>
+<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
+</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>
+<p>
+
+<DT><b>stream</b></DT>
+<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>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>
+
+<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.
+
+
+<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>
+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>
+<P>
+To delineate all basins with categories of streams:
+<P>
+<CODE>r.stream.basins dir=dirs stream=streams basins=bas_basins_elem</CODE>
+<P>
+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>
+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>
+echo '42=42
+* = NULL' > tmp #for one output
+
+echo '42 = 42
+252 = 252
+188 = 188
+* = NULL' >tmp #for multiple outputs
+
+r.reclass input=streams output=sel_streams_1 <tmp
+</PRE>
+</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>
+r.circle -b output=circle coordinate=639936.623832,216939.836449 max=200
+r.stream.basins -c dir=dirs stream=circle basins=bas_simul
+</PRE>
+</CODE>
+<P>
+To determine areas of contribution to streams of particular order use as streams the result of ordering:
+<P>
+<CODE>r.stream.basins dir=dirs stream=ord_strahler basins=bas_basin_strahler</CODE>
+<P>
+Determination of areas of potential source of pollution. The example will be done for lake marked with FULL_HYDR 8056 in North Carolina sample dataset. The lake shall be extracted and converted to binary raster map.
+
+<CODE>
+<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
+</PRE>
+</CODE>
+<P>
+See aslo tutorial: <a href="http://heretic.livenet.pl/heretic">http://heretic.livenet.pl/heretic</a>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.reclass.html">r.reclass</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2009 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Property changes on: grass-addons/raster/r.stream.basins/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.basins/global.h
===================================================================
--- grass-addons/raster/r.stream.basins/global.h (rev 0)
+++ grass-addons/raster/r.stream.basins/global.h 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,77 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+
+
+ /* define */
+
+/*define directions for code simplicity
+
+directions according to r.watershed: MUST check all directions
+|3|2|1|
+|4| |8|
+|5|6|7|
+
+*/
+
+#define POINT struct points
+POINT {
+ int r, c;
+ };
+
+#define OUTLET struct outs
+OUTLET {
+ int r, c;
+ int val;
+ };
+
+ /* functions.c */
+
+/* io.c */
+int open_raster(char *mapname);
+int create_maps(void);
+int max_link(void);
+int write_chatchment(void);
+int set_null(void);
+
+/* cachments */
+int find_outlets(void);
+int reset_catchments(void);
+int fill_catchments (OUTLET outlet);
+int fifo_insert (POINT point);
+POINT fifo_return_del (void);
+
+
+ /* variables */
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+GLOBAL struct Cell_head window;
+GLOBAL char *in_dirs, *in_streams; /* input dirrection and accumulation raster names*/
+GLOBAL char *name_catchments;
+GLOBAL int zeros, cats, lasts; /* flags */
+
+GLOBAL CELL **dirs, **streams; /* matrix with input data streams is used as output data*/
+
+GLOBAL int nrows, ncols;
+
+POINT *fifo_outlet;
+int tail, head;
+int outlets_num;
+int fifo_max;
+
+GLOBAL int out; /* number of strahler and horton outlets: index */
+OUTLET *outlets;
+
+GLOBAL struct History history; /* holds meta-data (title, comments,..) */
+
+
+
+
+
Property changes on: grass-addons/raster/r.stream.basins/global.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.basins/io.c
===================================================================
--- grass-addons/raster/r.stream.basins/io.c (rev 0)
+++ grass-addons/raster/r.stream.basins/io.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,142 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+int open_raster(char *mapname)
+{
+ int fd = 0;
+ char *mapset;
+ struct Cell_head cellhd; /* it stores region information, and header information of rasters */
+
+ mapset = G_find_cell2(mapname, ""); /* checking if map exist */
+ if (mapset == NULL) {
+ 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);
+
+ if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
+ 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);
+
+ 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"),
+ mapname, mapname);
+
+ return fd;
+}
+
+int create_maps(void)
+{
+ int r, c;
+ int in_dir_fd, in_stm_fd; /* input file descriptors: indir_fd - direction.... etc */
+
+ CELL *r_dirs, *r_streams;
+
+ 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..."));
+
+ 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);
+ G_close_cell(in_stm_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];
+
+ }
+
+ if (G_is_c_null_value(&r_streams[c])) {
+ streams[r][c] = 0;
+ }
+ else {
+ streams[r][c] = r_streams[c];
+ }
+
+ }
+
+ /* END dirs & streams & accums */
+
+
+ } /*end for r */
+
+ G_free(r_dirs);
+ G_free(r_streams);
+ G_percent(r, nrows, 2);
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+
+ return 0;
+} /* end create maps */
+
+int max_link(void)
+{
+ char *cur_mapset;
+ CELL c_min, c_max;
+ struct Range *stream_range;
+
+ G_init_range(stream_range);
+ cur_mapset = G_find_cell2(in_streams, "");
+ G_read_range(in_streams, cur_mapset, stream_range);
+ G_get_range_min_max(stream_range, &c_min, &c_max);
+ return c_max;
+} /* end_max_link */
+
+
+int write_chatchment()
+{
+ int r;
+ 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);
+
+ 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_free(streams);
+
+ return 0;
+}
+
+int set_null(void)
+{
+ int r, c;
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (streams[r][c] == 0)
+ G_set_c_null_value(&streams[r][c], 1);
+ }
+ }
+ return 0;
+}
Property changes on: grass-addons/raster/r.stream.basins/io.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.basins/main.c
===================================================================
--- grass-addons/raster/r.stream.basins/main.c (rev 0)
+++ grass-addons/raster/r.stream.basins/main.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,128 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.basins
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Calculate basins according user' input data.
+ * It uses r.stream.order or r.watershed stream map
+ * and r.watershed direction map.
+ * Stream input map shall contains streams or points outlels
+ * with or without unique own categories
+ * If input stream comes from r..stream.exteract direction map
+ * from r.stream.extract must be patched with that of r.watersheed.
+ *
+ * COPYRIGHT: (C) 2002,2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ * main function
+ *
+ *
+ */
+int main(int argc, char *argv[])
+{
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *in_dir_opt, *in_stm_opt, *out_opt; /* options */
+ struct Flag *out_back, *out_cat, *out_last; /* flags */
+
+ int link_max;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords = _("stream, order, catchments");
+ module->description = _("Calculate basins according user' input");
+
+ in_stm_opt = G_define_option(); /* input stream mask file - optional */
+ in_stm_opt->key = "stream";
+ in_stm_opt->type = TYPE_STRING;
+ in_stm_opt->required = YES; /* for now; TO DO: is planned to be optional */
+ in_stm_opt->gisprompt = "old,cell,raster";
+ in_stm_opt->description = "Name of stream mask input map";
+
+ in_dir_opt = G_define_option(); /* input directon file */
+ in_dir_opt->key = "dir";
+ in_dir_opt->type = TYPE_STRING;
+ in_dir_opt->required = YES;
+ in_dir_opt->gisprompt = "old,cell,raster";
+ in_dir_opt->description = "Name of flow direction input map";
+
+ /* output option - at least one is reqired */
+
+ out_opt = G_define_option();
+ out_opt->key = "basins";
+ out_opt->type = TYPE_STRING;
+ out_opt->required = YES;
+ out_opt->answer = NULL;
+ out_opt->gisprompt = "new,cell,raster";
+ out_opt->description = "Output basin map";
+
+
+ /* Define the different flags */
+ out_back = G_define_flag();
+ out_back->key = 'z';
+ out_back->description = _("Create zero-value background instead of NULL");
+
+ /* Define the different flags */
+ out_cat = G_define_flag();
+ out_cat->key = 'c';
+ out_cat->description =
+ _("Use unique category sequence instead of input streams");
+
+ /* Define the different flags */
+ out_last = G_define_flag();
+ out_last->key = 'l';
+ out_last->description = _("Create basins only for last stream links");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ name_catchments = out_opt->answer;
+ zeros = (out_back->answer != 0);
+ cats = (out_cat->answer != 0);
+ lasts = (out_last->answer != 0);
+
+ if (G_legal_filename(name_catchments) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), name_catchments);
+
+ G_get_window(&window);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ create_maps();
+ link_max = max_link();
+ find_outlets();
+
+ {
+ int j;
+
+ reset_catchments();
+ G_message(_("Calculate basins..."));
+ fifo_max = 4 * (nrows + ncols);
+ fifo_outlet = (POINT *) G_malloc((fifo_max + 1) * sizeof(POINT));
+ /* fifo queue for contributing cell */
+ for (j = 0; j < outlets_num; ++j) {
+ fill_catchments(outlets[j]);
+ }
+ G_free(fifo_outlet);
+ if (!zeros)
+ set_null();
+ write_chatchment();
+ }
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/raster/r.stream.basins/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/Makefile
===================================================================
--- grass-addons/raster/r.stream.distance/Makefile (rev 0)
+++ grass-addons/raster/r.stream.distance/Makefile 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,13 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.distance
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/raster/r.stream.distance/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/description.html
===================================================================
--- grass-addons/raster/r.stream.distance/description.html (rev 0)
+++ grass-addons/raster/r.stream.distance/description.html 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,61 @@
+<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
+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>
+<p>
+
+<DT><b>stream</b></DT>
+<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>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>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>
+
+
+
+<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>
+<DT><b>distance</b></DT>
+<DD>Returns distance to the targer (outlet, node stream) along watercoures. The map is of FCELL type</DD>
+</DL>
+
+<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.
+<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.
+
+<h2>NOTES</h2>
+<P>
+If there are more than one point or one stream networks and some separate points or separate streams networks are in catchment area defined by others it will results as in subbasin mode. In stream mode subbasin options is ommited. Input maps must be in CELL format (default output of r.watershed, r.stream.order and r.stream.extract)
+The distance are calculated in meters both for planimeters and Latitude-Longitude projections. The distance is calculated for flat areas not corrected by topography. Distance correction by topography may be done with following mapcalc formula:
+<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). .
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.reclass.html">r.reclass</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
Property changes on: grass-addons/raster/r.stream.distance/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/distance.c
===================================================================
--- grass-addons/raster/r.stream.distance/distance.c (rev 0)
+++ grass-addons/raster/r.stream.distance/distance.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,217 @@
+#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 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 = 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 */
+ }
+ 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];
+
+ 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 */
+
+ return 0;
+}
+
+int reset_distance(void)
+{
+ int r, c, i;
+
+ 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) {
+
+ distance[outlets[i].r][outlets[i].c] = 0;
+ }
+ }
+ return 0;
+}
+
+
+/*
+ algorithm uses fifo queue for determining basins area.
+ it can works in stream mode (calculate distance to and elevation above streams)
+ and outlet mode (calculate distance to and elevation above outlets)
+ */
+int fill_maps(OUTLET outlet)
+{
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ int r, c, val, i, j;
+ POINT n_cell;
+ float cur_dist = 0;
+ float tmp_dist = 0;
+ float target_elev; /* eleavation at stream or outlet */
+ float easting, northing;
+ float cell_easting, cell_northing;
+
+ tail = 0;
+ head = -1;
+ r = outlet.r;
+ c = outlet.c;
+
+ if (out_elev) {
+ target_elev = elevation[r][c];
+ elevation[r][c] = 0.;
+ }
+
+ while (tail != head) {
+ easting = window.west + (c + .5) * window.ew_res;
+ northing = window.north - (r + .5) * window.ns_res;
+
+ 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; /* border */
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (dirs[r + nextr[i]][c + nextc[i]] == j) { /* countributing cell */
+ /* restet distance and elevation on streams and outlets: */
+ if (outs) { /* outlet mode */
+ if (distance[r + nextr[i]][c + nextc[i]] == 0) {
+ continue; /* continue for loop, point simply is not added to the queue! */
+ }
+ else {
+ cell_northing =
+ window.north - (r + nextr[i] +
+ .5) * window.ns_res;
+ cell_easting =
+ window.west + (c + nextc[i] + .5) * window.ew_res;
+ cur_dist =
+ tmp_dist + G_distance(easting, northing,
+ cell_easting,
+ cell_northing);
+ distance[r + nextr[i]][c + nextc[i]] = cur_dist;
+ }
+ }
+ else { /* stream mode */
+ if (distance[r + nextr[i]][c + nextc[i]] == 0) {
+ cur_dist = 0;
+ if (out_elev)
+ target_elev =
+ elevation[r + nextr[i]][c + nextc[i]];
+ }
+ else {
+ cell_northing =
+ window.north - (r + nextr[i] +
+ .5) * window.ns_res;
+ cell_easting =
+ window.west + (c + nextc[i] + .5) * window.ew_res;
+ cur_dist =
+ tmp_dist + G_distance(easting, northing,
+ cell_easting,
+ cell_northing);
+ distance[r + nextr[i]][c + nextc[i]] = cur_dist;
+ }
+ } /* end stream mode */
+
+ if (out_elev) {
+ elevation[r + nextr[i]][c + nextc[i]] =
+ elevation[r + nextr[i]][c + nextc[i]] - target_elev;
+ n_cell.target_elev = target_elev;
+ }
+
+ n_cell.r = (r + nextr[i]);
+ n_cell.c = (c + nextc[i]);
+ n_cell.cur_dist = cur_dist;
+ fifo_insert(n_cell);
+ }
+ } /* end for i... */
+
+ n_cell = fifo_return_del();
+ r = n_cell.r;
+ c = n_cell.c;
+ tmp_dist = n_cell.cur_dist;
+ target_elev = n_cell.target_elev;
+
+ } /* end while */
+ return 0;
+}
+
+/* fifo functions */
+int fifo_insert(POINT point)
+{
+ fifo_outlet[tail++] = point;
+ if (tail > fifo_max)
+ tail = 0;
+ return 0;
+}
+
+POINT fifo_return_del(void)
+{
+ if (head > fifo_max)
+ head = -1;
+ return fifo_outlet[++head];
+}
Property changes on: grass-addons/raster/r.stream.distance/distance.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/global.h
===================================================================
--- grass-addons/raster/r.stream.distance/global.h (rev 0)
+++ grass-addons/raster/r.stream.distance/global.h 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,85 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+
+
+ /* define */
+
+/*define directions for code simplicity
+
+directions according to r.watershed: MUST check all directions
+|3|2|1|
+|4| |8|
+|5|6|7|
+
+*/
+#define SQRT2 1.414214
+#define POINT struct points
+POINT {
+ int r, c;
+float cur_dist;
+float target_elev;
+float northing;
+float easting;
+ };
+
+#define OUTLET struct outs
+OUTLET {
+ int r, c;
+ float northing;
+ float easting;
+ };
+
+ /* functions.c */
+
+/* io.c */
+int open_raster(char *mapname);
+int create_maps(void);
+void free_streams (void);
+int write_distance(void);
+int write_elevation(void);
+int set_null_elev(void);
+
+/* distance */
+int find_outlets(void);
+int reset_distance(void);
+int fill_maps(OUTLET outlet);
+int fifo_insert (POINT point);
+POINT fifo_return_del (void);
+
+
+ /* variables */
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+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 CELL **dirs, **streams; /* matrix with input data*/
+GLOBAL FCELL **elevation, **distance;
+/* streams and elevation are used to store internal data during processing */
+
+GLOBAL int nrows, ncols;
+
+POINT *fifo_outlet;
+int tail, head;
+int outlets_num;
+int fifo_max;
+
+GLOBAL int out; /* number of strahler and horton outlets: index */
+OUTLET *outlets;
+
+GLOBAL struct History history; /* holds meta-data (title, comments,..) */
+GLOBAL struct Cell_head window;
+
+
+
+
Property changes on: grass-addons/raster/r.stream.distance/global.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/io.c
===================================================================
--- grass-addons/raster/r.stream.distance/io.c (rev 0)
+++ grass-addons/raster/r.stream.distance/io.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,216 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+int open_raster(char *mapname)
+{
+ int fd = 0;
+ char *mapset;
+ struct Cell_head cellhd; /* it stores region information, and header information of rasters */
+
+ mapset = G_find_cell2(mapname, ""); /* checking if map exist */
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), mapname);
+
+ if (G_get_cellhd(mapname, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read file header of <%s>"), mapname);
+
+ if (mapname != in_elev) {
+ /* checking if the input maps type is CELL and region is valid */
+ if (G_raster_map_type(mapname, mapset) != CELL_TYPE)
+ G_fatal_error(_("<%s> is not of type CELL"), 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"),
+ mapname, mapname);
+ } /* end if map not elev */
+
+ if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
+ 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);
+
+ return fd;
+}
+
+int create_maps(void)
+{
+ int r, c;
+ int elev_type;
+ int in_dir_fd, in_stm_fd, in_dem_fd;
+ CELL *r_dirs, *r_streams;
+ FCELL *r_dem_f = NULL;
+ CELL *r_dem_c = NULL;
+ DCELL *r_dem_d = NULL;
+
+ 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);
+
+ if (in_elev) {
+ in_dem_fd = open_raster(in_elev);
+
+ elev_type = G_get_raster_map_type(in_dem_fd);
+ switch (elev_type) {
+ case CELL_TYPE:
+ r_dem_c = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ break;
+ case FCELL_TYPE:
+ r_dem_f = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+ break;
+ case DCELL_TYPE:
+ r_dem_d = (DCELL *) G_malloc(sizeof(DCELL) * ncols);
+ break;
+ }
+ elevation = (FCELL **) G_malloc(sizeof(FCELL *) * nrows);
+ }
+
+ G_message(_("Reading maps..."));
+
+ 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);
+ G_close_cell(in_stm_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+
+
+ if (in_elev) {
+ elevation[r] = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+
+ switch (elev_type) {
+ case CELL_TYPE: /* CELL */
+ if (G_get_c_raster_row(in_dem_fd, r_dem_c, r) < 0) {
+ G_close_cell(in_dem_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_dem_c[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_c[c];
+ }
+ } /* end for */
+ break; /* CELL */
+
+ case DCELL_TYPE: /* DCELL */
+ if (G_get_d_raster_row(in_dem_fd, r_dem_d, r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"),
+ r);
+ }
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_d_null_value(&r_dem_d[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_d[c];
+
+ }
+ } /* end for */
+ break; /* DCELL */
+
+ case FCELL_TYPE: /* FCELL */
+ if (G_get_f_raster_row(in_dem_fd, elevation[r], r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"),
+ r);
+ }
+ break; /* FCELL */
+
+ } /* end switch */
+ } /* end if elev */
+
+ 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])) {
+ streams[r][c] = 0;
+ }
+ else {
+ streams[r][c] = r_streams[c];
+ }
+ } /* end for c */
+
+ /* END dirs & streams & accums */
+
+ } /*end for r */
+ switch (elev_type) {
+ case CELL_TYPE:
+ G_free(r_dem_c);
+ break;
+ case DCELL_TYPE:
+ G_free(r_dem_d);
+ break;
+ }
+
+
+ G_free(r_streams);
+ G_free(r_dirs);
+ G_percent(r, nrows, 2);
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+
+ return 0;
+} /* end create maps */
+
+void free_streams(void)
+{
+ int r;
+
+ for (r = 0; r < nrows; ++r)
+ G_free(streams[r]);
+ G_free(streams);
+}
+
+int write_maps(char *mapname, FCELL ** map)
+{
+ int r, c;
+ int fd = 0;
+
+ if ((fd = G_open_raster_new(mapname, FCELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), mapname);
+
+ for (r = 0; r < nrows; ++r)
+ G_put_f_raster_row(fd, map[r]);
+
+ G_close_cell(fd);
+ G_short_history(mapname, "raster", &history);
+ G_write_history(mapname, &history);
+ G_message(_("%s Done"), mapname);
+ G_free(mapname);
+
+ return 0;
+}
+
+int set_null(FCELL ** map)
+{
+ int r, c;
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (distance[r][c] == -1)
+ G_set_f_null_value(&map[r][c], 1);
+ }
+ }
+ return 0;
+}
Property changes on: grass-addons/raster/r.stream.distance/io.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.distance/main.c
===================================================================
--- grass-addons/raster/r.stream.distance/main.c (rev 0)
+++ grass-addons/raster/r.stream.distance/main.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,164 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.basins
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Calculate distance and elevation over the streams and outlets
+ * according user's input data. The elevation and distance is calculated
+ * along watercourses
+ * It uses r.stream.order stream map map and r.watershed direction map.
+ * Stram input map shall contains streams or points oultels with or
+ * without unique categories
+ * If input stream comes from r.stream.exteract direction map
+ * from r.stream.extract dir map must be patched with that of r.watersheed.
+ *
+ * COPYRIGHT: (C) 2002,2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ * main function
+ *
+ *
+ */
+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 */
+
+ int link_max;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords = _("stream, order, catchments");
+ module->description =
+ _("Calculate distance to and alavation above streams \
+ and outlets according user' input. It can work in stream mode where target are streams and outlets mode \
+ where targets are outlets");
+
+ in_stm_opt = G_define_option(); /* input stream mask file - optional */
+ in_stm_opt->key = "stream";
+ in_stm_opt->type = TYPE_STRING;
+ in_stm_opt->required = YES; /* for now; TO DO: is planned to be optional */
+ in_stm_opt->gisprompt = "old,cell,raster";
+ in_stm_opt->description = "Name of streams (outlets) mask input map";
+
+ in_dir_opt = G_define_option(); /* input directon file */
+ in_dir_opt->key = "dir";
+ in_dir_opt->type = TYPE_STRING;
+ in_dir_opt->required = YES;
+ in_dir_opt->gisprompt = "old,cell,raster";
+ in_dir_opt->description = "Name of flow direction input map";
+
+ in_elev_opt = G_define_option(); /* input stream mask file - optional */
+ in_elev_opt->key = "dem";
+ in_elev_opt->type = TYPE_STRING;
+ in_elev_opt->required = NO;
+ in_elev_opt->gisprompt = "old,cell,raster";
+ in_elev_opt->description = "Name of elevation map";
+
+ /* output option - at least one is reqired */
+
+ out_dist_opt = G_define_option();
+ out_dist_opt->key = "distance";
+ out_dist_opt->type = TYPE_STRING;
+ out_dist_opt->required = NO;
+ out_dist_opt->answer = NULL;
+ out_dist_opt->gisprompt = "new,cell,raster";
+ out_dist_opt->description = "Output distance map";
+
+ out_elev_opt = G_define_option();
+ out_elev_opt->key = "elevation";
+ out_elev_opt->type = TYPE_STRING;
+ out_elev_opt->required = NO;
+ out_elev_opt->answer = NULL;
+ 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)");
+
+ /* Define the different flags */
+ out_sub = G_define_flag();
+ out_sub->key = 's';
+ out_sub->description =
+ _("Calculate parameters for subbasins (ignored in stream mode)");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ if (!out_elev_opt->answer && !out_dist_opt->answer)
+ G_fatal_error(_("You must select at least one output maps: distance and (or) elevation"));
+ if (out_elev_opt->answer && !in_elev_opt->answer)
+ G_fatal_error(_("If you select elevation output, DEM is required"));
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ in_elev = in_elev_opt->answer;
+ out_dist = out_dist_opt->answer;
+ out_elev = out_elev_opt->answer;
+ outs = (out_outs->answer != 0);
+ subs = (out_sub->answer != 0);
+
+ if (out_dist) {
+ if (G_legal_filename(out_dist) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_dist);
+ }
+ if (out_elev) {
+ if (G_legal_filename(out_elev) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_elev);
+ }
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ G_get_window(&window);
+
+ create_maps();
+ G_begin_distance_calculations();
+ find_outlets();
+ reset_distance();
+ free_streams();
+
+ /* 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) {
+ fill_maps(outlets[j]);
+ }
+ } /* end block */
+ G_free(fifo_outlet);
+
+ if (out_elev) {
+ set_null(elevation);
+ write_maps(out_elev, elevation);
+ }
+ if (out_dist) {
+ set_null(distance);
+ write_maps(out_dist, distance);
+ }
+
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/raster/r.stream.distance/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/Makefile
===================================================================
--- grass-addons/raster/r.stream.order/Makefile (rev 0)
+++ grass-addons/raster/r.stream.order/Makefile 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,12 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.order
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/raster/r.stream.order/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/description.html
===================================================================
--- grass-addons/raster/r.stream.order/description.html (rev 0)
+++ grass-addons/raster/r.stream.order/description.html 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,111 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-z</b></DT>
+<DD>Creates zero-value background instead of NULL. For some reason (like map algebra calculation) zero-valued background may be required. This flag produces zero-filled background instead of null (default).</DD>
+<p>
+<DT><b>stream</b></DT>
+<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>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 and <em>accumulation</em> maps 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>accum</b></DT>
+<DD>Flow accumulation (optional): name of flow accumulation file produced by r.watershed or r.stream.extract. This map is requied only if Horton's or Hack's ordering is performed. Flow accumulation map shall be of DCELL type, as is by default produced by r.watershed or r.stream.extract.</DD>
+<h2>OUTPUTS</h2>
+
+<P>At least one output map is required: </p>
+<DT><b>strahler</b></DT>
+<DD>Name of Strahler's stream order output map: see notes for detail. </DD>
+
+<DT><b>shreve</b></DT>
+<DD>Name of Shreve's stream magnitude output map: see notes for detail.</DD>
+
+<DT><b>horton</b></DT>
+<DD>Name of Horton's stream order output map (require accum file): see notes for detail.</DD>
+
+<DT><b>hack</b></DT>
+<DD>Name of Hack's main streams output map (require accum file): see notes for detail.</DD>
+
+</DL>
+
+<h2>DESCRIPTION</h2>
+<center>
+<h3>Stream ordering example:<h3>
+<img src=orders.png border=1><br>
+</center>
+
+<P>
+<H4>Strahler's stream order</H4>
+Strahler's stream order is a modification of Horton's streams order which fixes the ambiguity of Horton's ordering.
+In Strahler's ordering the main channel is not determined; instead the ordering is based on the hierarchy of tributaries. The
+ordering follows these rules:
+<OL>
+<li>if the node has no children, its Strahler order is 1.
+<li>if the node has one and only one tributuary with Strahler greatest order i, and all other tributuaries have order less than i, then the order remains i.
+<li>if the node has two or more tributuaries with greatest order i, then the Strahler order of the node is i + 1.
+</OL>
+Strahler's stream ordering starts in initial links which assigns order one. It proceeds downstream. At every node it verifies that there are at least 2 equal tributaries with maximum order. If not it continues with highest order, if yes it increases the node's order by 1 and continues downstream with new order.
+<BR>
+<B>Advantages and disadvantages of Strahler's ordering: </B>
+ Strahler's stream order has a good mathematical background. All catchments with streams in this context are directed graphs, oriented from the root towards the leaves. Equivalent definition of the Strahler number of a tree is that it is the height of the largest complete binary tree that can be homeomorphically embedded into the given tree; the Strahler number of a node in a tree is equivalent to the height of the largest complete binary tree that can be embedded below that node. The disadvantage of that methods is the lack of distinguishing a main channel which may interfere with the analytical process in highly elongated catchments
+
+<H4>Horton's stream order</H4>
+Horton's stream order applies to the stream as a whole but not to segments or links since the order on any channel remains unchanged from source till it "dies" in the higher order stream or in the outlet of the catchment. The main segment of the catchment gets the order of the whole catchment, while its tributaries get the order of their own subcatchments. The main difficulties of the Horton's order are criteria to be considered to distinguish between "true" first order segments and extension of higher order segments. Thqat is the reason why Horton's ordering has rather historical sense and is substituted by the more unequivocal Strahler's ordering system. There are no natural algorithms to order stream network according to Horton' paradigm. The algorithm used in r.stream.order requires to first calculate Strahler's stream order (downstream) and next recalculate to Horton ordering (upstream). To make a decision about proper ordering it uses first Strahler ordering, and next, if both branches have the same orders it uses flow accumulation to choose the actual link. The algorithm starts with the outlet, where the outlet link is assigned the corresponding Strahler order. Next it goes upstream and determines links according to Strahler ordering. If the orders of tributaries differ, the algorithm proceeds with the channel of highest order, if all orders are the same, it chooses that one with higher flow accumulation rate (higher catchment area). When it reaches the initial channel it goes back to the last undetermined branch, assign its Strahler order as Horton order and goes upstream to the next initial links. In that way stream orders remain unchanged from the point where Horton's order have been determined to the source.
+
+<BR>
+<B>Advantages and disadvantages of Horton's ordering:</B>
+The main advantages of Horton's ordering is that it produces natural stream ordering with main streams and its tributaries. The main disadvantage is that it requires prior Strahler's ordering. In some cases this may result in unnatural ordering, where the highest order will be ascribed not to the channel with higher accumulation but to the channel which leads to the most branched parts of the the catchment.
+<P>
+<H4>Shreve's stream magnitude</H4>
+That ordering method is similar to Consisted Associated Integers proposed by Scheidegger. It assigns magnitude of 1 for every initial channel. The magnitude of the following channel is the sum of magnitudes of its tributaries. The number of a particular link is the number of initials which contribute to it. To achive Consisted Associated Integers the result of Shreve's magnitude is to be multiplied by 2:
+<P>
+<code>r.mapcalc scheidegger=shreve*2</code>
+The algorithm is very similar to Strahler's algorithm, it proceeds downstream, and at every node the stream magnitude is the sum of its tributaries.
+<P>
+<H4>Hack's main streams order</H4>
+This method of ordering calculates main streams of main catchment and every subcatchments. Main stream of every catchment is set to 1, and consequently all its tributaries receive order 2. Their tributaries receive order 3 etc. The order of every stream remains constant up to its initial link. The route of every main stream is determined according to the maximum flow accumulation value of particular streams. So the main stream of every subcatchment is the stream with high accumulation. In most cases the main stream is the longest watercourse of the catchment, but in some cases, when a catchment consists of both rounded and elongated subcatchments these rules may not be maintained. The algorithm assigns 1 to every outlets stream and goes upstream according to maximum flow accumulation of every branch. When it reaches an initial stream it step back to the first unassigned confluence. It assigns order 2 to unordered tributaries and again goes upstream to the next initial stream. The process runs until all branches of all outlets are ordered.
+<BR>
+<B>Advantages and disadvantages of main stream ordering:</B>
+The biggest advantage of that method is the possibility to compare and analyze topology upstream, according to main streams. Because all tributaries of main channel have order of 2, streams can be quickly and easily filtered and its proprieties and relation to main stream determined. The main disadvantage of that method is the problem with the comparison of subcatchment topology of the same order. Subcatchments of the same order may be both highly branched and widespread in the catchment area and a small subcatchment with only one stream.
+
+<h2>NOTES</H2>
+<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>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.stream.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+</em>
+
+
+
+<h2>REFERENCES</h2>
+Hack, J., (1957), <i>Studies of longitudinal stream profiles in Virginia and Maryland</i>,
+<b>U.S. Geological Survey Professional Paper</b>, 294-B<p>
+Horton, R. E. (1945), <i>Erosional development of streams and their drainage basins: hydro-physical approach to quantitative morphology</i>,
+<b>Geological Society of America Bulletin</b> 56 (3): 275-370<BR>
+<a h
+<p>
+Shreve, R., <i>Statistical Law of Stream Numbers</i>, <b>J. Geol.</b>, 74, (1966), 17-37.<p>
+Strahler, A. N. (1952), <i>Hypsometric (area-altitude) analysis of erosional topology</i>,
+<b>Geological Society of America Bulletin</b> 63 (11): 1117–1142<p>
+Strahler, A. N. (1957), <i>Quantitative analysis of watershed geomorphology</i>,
+<b>Transactions of the American Geophysical Union</b> 8 (6): 913–920.<p>
+<h2>AUTHOR</h2>
+Jarek Jasiewicz, Markus Metz
+
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2009 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Property changes on: grass-addons/raster/r.stream.order/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: grass-addons/raster/r.stream.order/global.h
===================================================================
--- grass-addons/raster/r.stream.order/global.h (rev 0)
+++ grass-addons/raster/r.stream.order/global.h 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,86 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+
+ /* define */
+#define NODES struct node
+#define STREAMS struct stream
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+/* dirs
+ *** next ***
+
+ 3 | 2 | 1
+ --- --- ---
+ 4 | 0 | 8
+ --- --- ---
+ 5 | 6 | 7
+
+ *** prev ***
+
+ 7 | 6 | 5
+ --- --- ---
+ 8 | 0 | 4
+ --- --- ---
+ 1 | 2 | 3
+ */
+
+
+
+#define STREAM struct streem
+STREAM {
+ int stream; /* index */
+ int strahler, shreeve, horton, hack; /* orders */
+ int next_stream;
+ int trib_num;
+ int trib[5]; /* for now, I have problems with realock */
+
+ /* int *trib; */
+ double accum;
+};
+
+
+ /* functions.c */
+
+/* order.c */
+int trib_nums(int r, int c);
+int init_streams(int stream_num);
+int find_nodes(int stream_num);
+int strahler(void);
+int horton(void);
+int shreeve(void);
+int hack(void);
+
+/* io.c */
+int open_raster(char *mapname);
+int create_base_maps(void);
+int stream_number(void);
+int open_accum(void);
+int destroy_c_data(CELL ** data_name);
+int destroy_d_data(DCELL ** data_name);
+int write_maps(void);
+
+
+ /* variables */
+
+GLOBAL struct Cell_head window;
+GLOBAL char *in_dirs, *in_streams, *in_accum; /* input dirrection and accumulation raster names */
+GLOBAL char *out_strahler, *out_shreeve, *out_hack, *out_horton; /* output strahler, shreeve and class raster names */
+GLOBAL int out_zero; /* zero-value background */
+
+GLOBAL CELL **dirs, **streams; /* matrix with input data */
+GLOBAL DCELL **accum; /* accum file is of type double */
+
+GLOBAL int stack_max;
+GLOBAL int nrows, ncols;
+
+GLOBAL STREAM *s_streams; /* stream structure all parameters we have here */
+int *springs, *outlets;
+int springs_num, outlets_num;
Property changes on: grass-addons/raster/r.stream.order/global.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/io.c
===================================================================
--- grass-addons/raster/r.stream.order/io.c (rev 0)
+++ grass-addons/raster/r.stream.order/io.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,246 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+int open_raster(char *mapname)
+{
+
+ int fd = 0;
+ char *mapset;
+ struct Cell_head cellhd; /* it stores region information, and header information of rasters */
+
+ mapset = G_find_cell2(mapname, ""); /* checking if map exist */
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), mapname);
+
+ if (mapname != in_accum) {
+ 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);
+ }
+
+ if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
+ 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);
+
+ 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"),
+ mapname, mapname);
+
+ return fd;
+}
+
+
+int create_base_maps(void)
+{
+ int r, c;
+ CELL *r_dirs=NULL, *r_streams=NULL;
+ DCELL *r_accum=NULL;
+ int in_dir_fd=0, in_stm_fd=0, in_acc_fd=0; /* input file descriptors: indir_fd - direction.... etc */
+
+ in_dir_fd = open_raster(in_dirs);
+ in_stm_fd = open_raster(in_streams);
+ if (in_accum)
+ in_acc_fd = open_raster(in_accum);
+
+ dirs = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+ streams = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+ if (in_accum)
+ accum = (DCELL **) G_malloc(sizeof(DCELL *) * nrows);
+
+ r_dirs = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ r_streams = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ if (in_accum)
+ r_accum = (DCELL *) G_malloc(sizeof(DCELL) * ncols);
+
+ G_message(_("Reading maps..."));
+
+ for (r = 0; r < nrows; ++r) {
+ G_percent(r, nrows, 2);
+ /* dirs & streams */
+ dirs[r] = G_malloc(sizeof(CELL) * ncols);
+ streams[r] = 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);
+ G_close_cell(in_stm_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+
+ if (in_accum) {
+ accum[r] = (DCELL *) G_malloc(sizeof(DCELL) * ncols);
+ if (G_get_d_raster_row(in_acc_fd, r_accum, r) < 0) {
+ G_close_cell(in_acc_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];
+ }
+ if (G_is_c_null_value(&r_streams[c])) {
+ streams[r][c] = 0;
+ }
+ else {
+ streams[r][c] = r_streams[c];
+ }
+
+ if (in_accum) {
+ if (G_is_d_null_value(&r_accum[c])) {
+ accum[r][c] = 0;
+ }
+ else {
+ accum[r][c] = r_accum[c];
+ }
+ } /* end if accum */
+
+ } /* end for c */
+ }
+ G_free(r_streams);
+ G_free(r_dirs);
+ if (in_accum)
+ G_free(r_accum);
+ G_percent(r, nrows, 2);
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+ if (in_accum)
+ G_close_cell(in_acc_fd);
+ return 0;
+} /* end create base maps */
+
+int stream_number(void)
+{
+ char *cur_mapset;
+ CELL c_min, c_max;
+ struct Range *stream_range;
+ G_init_range(stream_range);
+ cur_mapset = G_find_cell2(in_streams, "");
+ G_read_range(in_streams,cur_mapset,stream_range);
+ G_get_range_min_max(stream_range, &c_min, &c_max);
+ return c_max;
+}
+
+int write_maps(void)
+{
+ int r, c;
+ CELL *strahler_buf, *shreeve_buf, *hack_buf, *horton_buf;
+ int out_str_fd, out_shr_fd, out_hck_fd, out_hrt_fd; /* output file descriptors: outstr_fd - strahler... etc */
+ struct History history; /* holds meta-data (title, comments,..) */
+
+ /* routine check if legal file names and we able to opuen files */
+ if (out_strahler) {
+ if ((out_str_fd = G_open_raster_new(out_strahler, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"),
+ out_strahler);
+ strahler_buf = G_allocate_cell_buf();
+ }
+
+ if (out_shreeve) {
+ if ((out_shr_fd = G_open_raster_new(out_shreeve, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), out_shreeve);
+ shreeve_buf = G_allocate_cell_buf();
+ }
+
+ if (out_hack) {
+ if ((out_hck_fd = G_open_raster_new(out_hack, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), out_hack);
+ hack_buf = G_allocate_cell_buf();
+ }
+
+ if (out_horton) {
+ if ((out_hrt_fd = G_open_raster_new(out_horton, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), out_horton);
+ horton_buf = G_allocate_cell_buf();
+ }
+
+ G_message(_("Writing maps..."));
+
+ for (r = 0; r < nrows; ++r) {
+ G_percent(r, nrows, 2);
+
+ if (out_strahler)
+ G_set_c_null_value(strahler_buf, ncols);
+ if (out_shreeve)
+ G_set_c_null_value(shreeve_buf, ncols);
+ if (out_hack)
+ G_set_c_null_value(hack_buf, ncols);
+ if (out_horton)
+ G_set_c_null_value(horton_buf, ncols);
+
+ for (c = 0; c < ncols; ++c) {
+ if (!out_zero) {
+ if (streams[r][c] > 0) {
+ if (out_strahler)
+ strahler_buf[c] = s_streams[streams[r][c]].strahler;
+ if (out_shreeve)
+ shreeve_buf[c] = s_streams[streams[r][c]].shreeve;
+ if (out_hack)
+ hack_buf[c] = s_streams[streams[r][c]].hack;
+ if (out_horton)
+ horton_buf[c] = s_streams[streams[r][c]].horton;
+ } /* end if streams */
+ }
+ else if (out_zero) {
+ if (out_strahler)
+ strahler_buf[c] = s_streams[streams[r][c]].strahler;
+ if (out_shreeve)
+ shreeve_buf[c] = s_streams[streams[r][c]].shreeve;
+ if (out_hack)
+ hack_buf[c] = s_streams[streams[r][c]].hack;
+ if (out_horton)
+ horton_buf[c] = s_streams[streams[r][c]].horton;
+ }
+ } /* end for cols */
+
+ if (out_strahler)
+ G_put_c_raster_row(out_str_fd, strahler_buf);
+ if (out_shreeve)
+ G_put_c_raster_row(out_shr_fd, shreeve_buf);
+ if (out_hack)
+ G_put_c_raster_row(out_hck_fd, hack_buf);
+ if (out_horton)
+ G_put_c_raster_row(out_hrt_fd, horton_buf);
+ } /* end for r */
+
+ G_percent(r, nrows, 2);
+
+ if (out_strahler) {
+ G_close_cell(out_str_fd);
+ G_free(strahler_buf);
+ G_short_history(out_strahler, "raster", &history);
+ G_write_history(out_strahler, &history);
+ G_message(_("%s Done!"),out_strahler);
+ }
+
+ if (out_shreeve) {
+ G_close_cell(out_shr_fd);
+ G_free(shreeve_buf);
+ G_short_history(out_shreeve, "raster", &history);
+ G_write_history(out_shreeve, &history);
+ G_message(_("%s Done!"),out_shreeve);
+ }
+ if (out_hack) {
+ G_close_cell(out_hck_fd);
+ G_free(hack_buf);
+ G_short_history(out_hack, "raster", &history);
+ G_write_history(out_hack, &history);
+ G_message(_("%s Done!"),out_hack);
+ }
+
+ if (out_horton) {
+ G_close_cell(out_hrt_fd);
+ G_free(horton_buf);
+ G_short_history(out_horton, "raster", &history);
+ G_write_history(out_horton, &history);
+ G_message(_("%s Done!"),out_horton);
+ }
+
+ return 0;
+
+} /* end write_maps */
Property changes on: grass-addons/raster/r.stream.order/io.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/main.c
===================================================================
--- grass-addons/raster/r.stream.order/main.c (rev 0)
+++ grass-addons/raster/r.stream.order/main.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,173 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.order
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Calculate Strahler's and Horton's streams order maps, Hack's main streams map and Shreeve's stream magnitude
+* It use r.stream.extract or r.watershed output files: stream, direction and accumulation.
+*
+* COPYRIGHT: (C) 2002,2009 by the GRASS Development Team
+*
+* This program is free software under the GNU General Public
+* License (>=v2). Read the file COPYING that comes with GRASS
+* for details.
+*
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ * main function *
+ */
+int main(int argc, char *argv[])
+{
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *in_dir_opt, *in_stm_opt, *in_acc_opt, *out_str_opt, *out_shr_opt, *out_hck_opt, *out_hrt_opt; /* options */
+ struct Flag *out_back; /* flags */
+ int stream_num;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords =
+ _("stream, order, magnitude, Horton, Strahler, Shreeve");
+ module->description =
+ _("Calculate Strahler's and Horton's stream order Hack's main streams and Shreeve's stream magnitude. It use r.watershed or r.stream.extract output files: stream, direction and optionally accumulation. Otput data can be eighter from r.watershed or r.stream.extract but not from both together");
+
+ /*input option direction is reqired, acummulation is optional */
+
+ in_stm_opt = G_define_option(); /* input stream mask file - optional */
+ in_stm_opt->key = "stream"; /* required if strahler stream order is calculated for existing stream network */
+ in_stm_opt->type = TYPE_STRING;
+ in_stm_opt->required = YES; /* for now; TO DO: is planned to be optional */
+ in_stm_opt->gisprompt = "old,cell,raster";
+ in_stm_opt->description =
+ "Name of stream mask input map (r.watershed or r.stream.extract)";
+
+ in_dir_opt = G_define_option(); /* input directon file */
+ in_dir_opt->key = "dir";
+ in_dir_opt->type = TYPE_STRING;
+ in_dir_opt->required = YES;
+ in_dir_opt->gisprompt = "old,cell,raster";
+ in_dir_opt->description =
+ "Name of direction input map (r.watershed or r.stream.extract)";
+
+ in_acc_opt = G_define_option(); /* input stream mask file - optional */
+ in_acc_opt->key = "accum"; /* required if strahler stream order is calculated for existing stream network */
+ in_acc_opt->type = TYPE_STRING;
+ in_acc_opt->required = NO; /* for now; TO DO: is planned to be optional */
+ in_acc_opt->answer = NULL;
+ in_acc_opt->gisprompt = "old,cell,raster";
+ in_acc_opt->description =
+ "Name of accumulation file (r.watershed or r.stream.extract)";
+
+ /* output option - at least one is reqired */
+
+ out_str_opt = G_define_option(); /*Strahler stream order */
+ out_str_opt->key = "strahler";
+ out_str_opt->type = TYPE_STRING;
+ out_str_opt->required = NO;
+ out_str_opt->answer = NULL;
+ out_str_opt->gisprompt = "new,cell,raster";
+ out_str_opt->description = "Name of Strahler's stream order output map";
+ out_str_opt->guisection = _("Output options");
+
+ out_shr_opt = G_define_option();
+ out_shr_opt->key = "shreve";
+ out_shr_opt->type = TYPE_STRING;
+ out_shr_opt->required = NO;
+ out_shr_opt->answer = NULL;
+ out_shr_opt->gisprompt = "new,cell,raster";
+ out_shr_opt->description = "Name of Shreve's stream magnitude output map";
+ out_shr_opt->guisection = _("Output options");
+
+ out_hrt_opt = G_define_option();
+ out_hrt_opt->key = "horton";
+ out_hrt_opt->type = TYPE_STRING;
+ out_hrt_opt->required = NO;
+ out_hrt_opt->answer = NULL;
+ out_hrt_opt->gisprompt = "new,cell,raster";
+ out_hrt_opt->description =
+ "Name of Horton's stream order output map (require accum file)";
+ out_hrt_opt->guisection = _("Output options");
+
+ out_hck_opt = G_define_option();
+ out_hck_opt->key = "hack";
+ out_hck_opt->type = TYPE_STRING;
+ out_hck_opt->required = NO;
+ out_hck_opt->answer = NULL;
+ out_hck_opt->gisprompt = "new,cell,raster";
+ out_hck_opt->description =
+ "Name of Hack's main streams output map (require accum file)";
+ out_hck_opt->guisection = _("Output options");
+
+ /* Define flags */
+ out_back = G_define_flag();
+ out_back->key = 'z';
+ out_back->description = _("Create zero-value background instead of NULL");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ G_get_window(&window);
+
+ if (!out_str_opt->answer && !out_shr_opt->answer && !out_hck_opt->answer
+ && !out_hrt_opt->answer)
+ G_fatal_error(_("You must select one or more output maps: strahler, horton, shreeve or hack"));
+ if (out_hrt_opt->answer && !in_acc_opt->answer)
+ G_fatal_error(_("If you select horton or hack, accum is neccessary"));
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ in_accum = in_acc_opt->answer;
+
+ /* stores output options to variables */
+ out_strahler = out_str_opt->answer;
+ out_shreeve = out_shr_opt->answer;
+ out_hack = out_hck_opt->answer;
+ out_horton = out_hrt_opt->answer;
+ out_zero = (out_back->answer != 0);
+
+ if (out_strahler) {
+ if (G_legal_filename(out_strahler) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_strahler);
+ }
+ if (out_shreeve) {
+ if (G_legal_filename(out_shreeve) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_shreeve);
+ }
+ if (out_hack) {
+ if (G_legal_filename(out_hack) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_hack);
+ }
+ if (out_horton) {
+ if (G_legal_filename(out_horton) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_horton);
+ }
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ create_base_maps();
+ stream_num = stream_number();
+ stack_max = stream_num; /* stack's size depends on number of streams */
+ init_streams(stream_num);
+ find_nodes(stream_num);
+ if (out_strahler || out_horton)
+ strahler();
+ if (out_shreeve)
+ shreeve();
+ if (out_horton)
+ horton();
+ if (out_hack)
+ hack();
+ write_maps();
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/raster/r.stream.order/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/order.c
===================================================================
--- grass-addons/raster/r.stream.order/order.c (rev 0)
+++ grass-addons/raster/r.stream.order/order.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,398 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int trib_nums(int r, int c)
+{ /* calcualte number of tributuaries */
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int trib = 0;
+ int i, j;
+
+ 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;
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
+ dirs[r + nextr[i]][c + nextc[i]] == j)
+ trib++;
+ }
+ if (trib > 5)
+ G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
+ if (trib > 3)
+ G_warning(_("Stream network may be too dense..."));
+ return trib;
+} /* end trib_num */
+
+int init_streams(int stream_num)
+{
+ int i;
+ s_streams = (STREAM *) G_malloc((stream_num + 1) * sizeof(STREAM));
+
+ for (i = 0; i <= stream_num; ++i) {
+ s_streams[i].next_stream = -1;
+ s_streams[i].stream = -1;
+ s_streams[i].strahler = -1;
+ s_streams[i].shreeve = -1;
+ s_streams[i].hack = -1;
+ s_streams[i].horton = -1;
+ s_streams[i].accum = 0.0;
+ s_streams[i].trib[0] = 0;
+ s_streams[i].trib[1] = 0;
+ s_streams[i].trib[2] = 0;
+ s_streams[i].trib[3] = 0;
+ s_streams[i].trib[4] = 0;
+ }
+ return 0;
+}
+
+int find_nodes(int stream_num)
+{
+
+ int d, i, j; /* d: direction, i: iteration */
+ int r, c;
+ int trib_num, trib = 0;
+ int next_stream = -1, cur_stream;
+
+ springs_num = 0, outlets_num = 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 };
+
+ G_message(_("Finding nodes..."));
+
+ outlets = (int *)G_malloc((stream_num) * sizeof(int));
+ springs = (int *)G_malloc((stream_num) * sizeof(int));
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (streams[r][c] > 0) {
+ trib_num = trib_nums(r, c);
+ trib = 0;
+ 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;
+ }
+ else {
+ next_stream = streams[r + nextr[d]][c + nextc[d]];
+ if (next_stream < 1)
+ next_stream = -1;
+ }
+ if (dirs[r][c]==0)
+ next_stream=-1;
+ cur_stream = streams[r][c];
+
+ if (cur_stream != next_stream) { /* building hierarchy */
+
+ if (outlets_num > (stream_num - 1))
+ G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
+
+ s_streams[cur_stream].stream = cur_stream;
+ s_streams[cur_stream].next_stream = next_stream;
+ if (next_stream < 0)
+ outlets[outlets_num++] = cur_stream; /* collecting outlets */
+ }
+
+ if (trib_num == 0) { /* adding springs */
+
+ if (springs_num > (stream_num - 1))
+ G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
+
+ s_streams[cur_stream].trib_num = 0;
+ springs[springs_num++] = cur_stream; /* collecting springs */
+ }
+
+ if (trib_num > 1) { /* adding tributuaries */
+ s_streams[cur_stream].trib_num = trib_num;
+ for (i = 1; i < 9; ++i) {
+
+ if (trib > 4)
+ G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
+
+ if (r + nextr[i] < 0 || r + nextr[i] > (nrows - 1) ||
+ c + nextc[i] < 0 || c + nextc[i] > (ncols - 1))
+ continue;
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 &&
+ dirs[r + nextr[i]][c + nextc[i]] == j) {
+ if (in_accum) { /* only if horton and hack */
+ s_streams[streams[r + nextr[i]]
+ [c + nextc[i]]].accum =
+ fabs(accum[r + nextr[i]][c + nextc[i]]);
+ }
+ s_streams[cur_stream].trib[trib++] =
+ streams[r + nextr[i]][c + nextc[i]];
+ }
+ } /* end for i... */
+
+ } /* if trib_num>1 */
+ } /* end if streams */
+ } /* c */
+ } /* r */
+ return 0;
+}
+
+/*
+ All algorithms used in analysis ar not recursive. For Strahler order and Shreve magnitude starts from initial channel and proceed downstream. Algortitms try to assgin order for branch and if it is imposible start from next initial channel, till all branches are ordered.
+ For Hortor and Hack ordering it proceed upstream and uses stack data structure to determine unordered branch.
+ Strahler stream order algorithm according idea of Victor Olaya' (SAGA GIS), but without recusion.
+ Algorithm of Hack main stram according idea of Markus Metz .
+ */
+
+
+int strahler(void)
+{
+
+ int i, j, done = 1;
+ int cur_stream, next_stream;
+ int max_strahler = 0, max_strahler_num;
+
+ G_message(_("Calculating Strahler's stream order ..."));
+
+ for (j = 0; j < springs_num; ++j) { /* main loop on springs */
+
+ cur_stream = s_streams[springs[j]].stream;
+ do { /* we must go at least once, if stream is of first order and is outlet */
+ max_strahler_num = 1;
+ max_strahler = 0;
+ next_stream = s_streams[cur_stream].next_stream;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* assign 1 for spring stream */
+ s_streams[cur_stream].strahler = 1;
+ cur_stream = next_stream;
+ done = 1;
+ }
+ else {
+ done = 1;
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* loop for determining strahler */
+ if (s_streams[s_streams[cur_stream].trib[i]].strahler < 0) {
+ done = 0;
+ break; /* strahler is not determined, break for loop */
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler > max_strahler) {
+ max_strahler =
+ s_streams[s_streams[cur_stream].trib[i]].strahler;
+ max_strahler_num = 1;
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler == max_strahler) {
+ ++max_strahler_num;
+ }
+ } /* end determining strahler */
+
+ if (done == 1) {
+ s_streams[cur_stream].strahler =
+ (max_strahler_num >
+ 1) ? ++max_strahler : max_strahler;
+ cur_stream = next_stream; /* if next_stream<0 we in outlet stream */
+ }
+
+ }
+ } while (done && next_stream > 0);
+ } /* end for of main loop */
+ return 0;
+} /* end strahler */
+
+int shreeve(void)
+{
+
+ int i, j, done = 1;
+ int cur_stream, next_stream;
+ int max_shreeve = 0;
+
+ G_message(_("Calculating Shreeve's stream magnitude ..."));
+
+ for (j = 0; j < springs_num; ++j) { /* main loop on springs */
+
+ cur_stream = s_streams[springs[j]].stream;
+ do { /* we must go at least once, if stream is of first order and is outlet */
+
+ max_shreeve = 0;
+ next_stream = s_streams[cur_stream].next_stream;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* assign 1 for spring stream */
+
+ s_streams[cur_stream].shreeve = 1;
+ cur_stream = next_stream;
+ done = 1;
+
+ }
+ else {
+ done = 1;
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* loop for determining strahler */
+ if (s_streams[s_streams[cur_stream].trib[i]].shreeve < 0) {
+ done = 0;
+ break; /* shreeve is not determined, break for loop */
+ }
+ else {
+ max_shreeve +=
+ s_streams[s_streams[cur_stream].trib[i]].shreeve;
+ }
+ } /* end determining strahler */
+
+ if (done == 1) {
+ s_streams[cur_stream].shreeve = max_shreeve;
+ cur_stream = next_stream; /* if next_stream<0 we in outlet stream */
+ }
+ }
+
+ } while (done && next_stream > 0);
+ } /* end main loop */
+ return 0;
+} /* end shreeve */
+
+int horton(void)
+{
+
+ int *stack;
+ int top, i, j;
+ int cur_stream, cur_horton;
+ int max_strahler;
+ double max_accum;
+ int up_stream = 0;
+
+ G_message(_("Calculating Hortons's stream order ..."));
+ stack = (int *)G_malloc(stack_max * sizeof(int));
+
+ for (j = 0; j < outlets_num; ++j) {
+ cur_stream = s_streams[outlets[j]].stream; /* outlet: init */
+ cur_horton = s_streams[cur_stream].strahler;
+ stack[0] = 0;
+ stack[1] = cur_stream;
+ top = 1;
+
+ do { /* on every stream */
+ max_strahler = 0;
+ max_accum = 0;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* spring: go back on stack */
+
+ s_streams[cur_stream].horton = cur_horton;
+ cur_stream = stack[--top];
+
+ }
+ else if (s_streams[cur_stream].trib_num > 1) { /* node */
+
+ up_stream = 0; /* calculating up_stream */
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) {
+ if (s_streams[s_streams[cur_stream].trib[i]].horton < 0) {
+
+ if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler > max_strahler) {
+ max_strahler =
+ s_streams[s_streams[cur_stream].trib[i]].
+ strahler;
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+
+ }
+ else if (s_streams[s_streams[cur_stream].trib[i]].
+ strahler == max_strahler) {
+
+ if (s_streams[s_streams[cur_stream].trib[i]].
+ accum > max_accum) {
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+ }
+ }
+ }
+ } /* end determining up_stream */
+
+ if (up_stream) { /* at least one branch is not assigned */
+ if (s_streams[cur_stream].horton < 0) {
+ s_streams[cur_stream].horton = cur_horton;
+ }
+ else {
+ cur_horton = s_streams[up_stream].strahler;
+ }
+ cur_stream = up_stream;
+ stack[++top] = cur_stream;
+
+ }
+ else { /* all asigned, go downstream */
+ cur_stream = stack[--top];
+
+ } /* end up_stream */
+ } /* end spring/node */
+ } while (cur_stream);
+ } /* end for outlets */
+ G_free(stack);
+ return 0;
+}
+
+int hack(void)
+{
+
+ int *stack;
+ int top, i, j;
+ int cur_stream, cur_hack;
+ double max_accum;
+ int up_stream = 0;
+
+ G_message(_("Calculating Hack's main streams ..."));
+ stack = (int *)G_malloc(stack_max * sizeof(int));
+
+ for (j = 0; j < outlets_num; ++j) {
+ cur_stream = s_streams[outlets[j]].stream; /* outlet: init */
+ cur_hack = 1;
+ stack[0] = 0;
+ stack[1] = cur_stream;
+ top = 1;
+
+ do {
+ max_accum = 0;
+
+ if (s_streams[cur_stream].trib_num == 0) { /* spring: go back on stack */
+
+ s_streams[cur_stream].hack = cur_hack;
+ cur_stream = stack[--top];
+
+ }
+ else if (s_streams[cur_stream].trib_num > 1) { /* node */
+ up_stream = 0; /* calculating up_stream */
+
+ for (i = 0; i < s_streams[cur_stream].trib_num; ++i) { /* determining upstream */
+ if (s_streams[s_streams[cur_stream].trib[i]].hack < 0) {
+ if (s_streams[s_streams[cur_stream].trib[i]].accum >
+ max_accum) {
+ max_accum =
+ s_streams[s_streams[cur_stream].trib[i]].
+ accum;
+ up_stream = s_streams[cur_stream].trib[i];
+ }
+ }
+ } /* end determining up_stream */
+
+ if (up_stream) { /* at least one branch is not assigned */
+ if (s_streams[cur_stream].hack < 0) {
+ s_streams[cur_stream].hack = cur_hack;
+ }
+ else {
+ cur_hack = s_streams[cur_stream].hack;
+ ++cur_hack;
+ }
+ cur_stream = up_stream;
+ stack[++top] = cur_stream;
+
+ }
+ else { /* all asigned, go downstream */
+
+ cur_stream = stack[--top];
+
+ } /* end up_stream */
+ } /* end spring/node */
+ } while (cur_stream);
+ } /* end for outlets */
+ G_free(stack);
+ return 0;
+}
Property changes on: grass-addons/raster/r.stream.order/order.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.order/orders.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.stream.order/orders.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/raster/r.stream.stats/Makefile
===================================================================
--- grass-addons/raster/r.stream.stats/Makefile (rev 0)
+++ grass-addons/raster/r.stream.stats/Makefile 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,13 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.stats
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/raster/r.stream.stats/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/description.html
===================================================================
--- grass-addons/raster/r.stream.stats/description.html (rev 0)
+++ grass-addons/raster/r.stream.stats/description.html 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,92 @@
+<h2>OPTIONS</h2>
+<DL>
+
+<DT><b>stream</b></DT>
+<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>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>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>
+
+
+
+<h2>OUTPUTS</h2>
+Output statistics are send to standard output. If there are no errors no addational messages are send to standard output. To redirect output to file use redarection operators: > or >>.
+</DL>
+
+<h2>DESCRIPTION</h2>
+<P>
+Module r.stream.stats is prepared to calculate Hotron's statistics of drainage network.
+<P>
+These statistics are calculated according formulas given by R.Horton (1945). Because Horton do not defined precisely what is stream slope, I proposed 2 different approaches: first (slope) use cell-by-cell slope calculation, second (gradient) use difference between elevation of outlet and source of every channel to its length to calculate formula. Bifurcation ratio for every order is calculated acording formula:
+<P>
+<CODE>n_streams[1]/n_stream[i+1]</CODE>
+<P>
+where i the current order and i+1 next higher order. For max order of the map number of streams is zero. Rest of the ratios are calculated in similar mode. The bifurcation and other ratios for the whole catchment (map) is calculated as mean i.e sum of all bifurcation ratio / max_order-1 (for max_order stream bifurcation ratio = 0)
+It is strongly recommended to extract stream network using basin map created with r.stream.basin. If whole stream order map is used the calculation will be performed but results may not have hydrological sense.
+
+For every order (std) means that statstic is calculated with standard deviation:
+<UL>
+<li>number of streams
+<li>total length of streams of given order
+<li>total area of basins of given order
+<li>drainage density
+<li>stream density
+
+<li>average length of streams of given order (std)
+<li>average slope (cell by cell inclination) of streams of given order (std)
+<li>average gradient (spring to outlet inclination ) of streams of given order (std)
+<li>average area of basins of given order (std)
+<li>avarage elevation difference of given order (std)
+<P>ratios:
+<li>bifuracation ratio
+<li>length ratio
+<li>sloope and gradient ratios
+<li>area ratio
+</UL>
+for the whole basin:
+<UL>
+<li>total number of streams
+<li>total length of streams
+<li>total basin area
+<li>drainage density
+<li>stream density
+<P>ratios:
+<li>bifurcation ratio (std)
+<li>length ratio (std)
+<li>slope and gradient ratios (std)
+<li>area ratio (std)
+</ul>
+
+
+<h2>NOTES</h2>
+<P>
+Module calculates statistics for all streams in input stream map.It is strongly recomended to extract only network of one basin, but it is not necessary for computation. Streams for desired basin first can be extracted with following mapcalc formula:
+
+<P>
+<CODE>echo 'sel_streams=if(basin==xxx,streams,null())'|r.mapcalc #xxx category of desired basin<CODE>
+<P>
+
+It is also possible to calculate Horton's statistics for Shreve ordering but it has limited hydrological sense. Hack main stream is not the same what so called Horton's reverse ordering.
+<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>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.reclass.html">r.reclass</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
Property changes on: grass-addons/raster/r.stream.stats/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/global.h
===================================================================
--- grass-addons/raster/r.stream.stats/global.h (rev 0)
+++ grass-addons/raster/r.stream.stats/global.h 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,120 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+
+
+ /* define */
+
+/*define directions for code simplicity
+
+directions according to r.watershed: MUST check all directions
+|3|2|1|
+|4| |8|
+|5|6|7|
+
+*/
+#define SQRT2 1.414214
+#define POINT struct points
+POINT {
+ int r, c;
+ };
+
+#define STREAM struct strs
+STREAM {
+ int index;
+ int r, c; /* outlet */
+ float elev_diff;
+ float elev_spring, elev_outlet;
+ float slope; /* cumulative */
+ float gradient;
+ float length; /* cumulative */
+ int order;
+ double basin_area; /* basin order */
+ int cell_num;
+ };
+
+
+#define STATS struct statistics
+STATS {
+ int order;
+ int stream_num;
+ float sum_length;
+ float avg_length;
+ float std_length;
+ float avg_slope;
+ float std_slope;
+ float avg_gradient;
+ float std_gradient;
+ double sum_area;
+ double avg_area;
+ double std_area;
+ float avg_elev_diff;
+ float std_elev_diff;
+ float bifur_ratio;
+ float std_bifur_ratio;
+ float length_ratio;
+ float std_length_ratio;
+ float area_ratio;
+ float std_area_ratio;
+ float slope_ratio;
+ float std_slope_ratio;
+ float gradient_ratio;
+ float std_gradient_ratio;
+ float stream_frequency;
+ float drainage_density;
+};
+
+ /* functions.c */
+
+/* io.c */
+int open_raster(char *mapname);
+int create_maps(void);
+
+/* stats */
+int init_streams (void);
+int fill_basin (int r,int c);
+int calculate_basins (void);
+int stats(void);
+int fifo_insert (POINT point);
+POINT fifo_return_del (void);
+
+/* print stats */
+int print_stats(void);
+
+ /* variables */
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+GLOBAL char *in_dirs, *in_streams, *in_elev; /* input dirrection and accumulation raster names*/
+GLOBAL char *out_file;
+GLOBAL int hack; /* flags */
+
+
+GLOBAL CELL **dirs, **streams; /* matrix with input data*/
+GLOBAL FCELL**elevation;
+/* streams and elevation are used to store internal data during processing */
+
+GLOBAL int nrows, ncols;
+
+POINT *fifo_outlet;
+int tail, head;
+int fifo_max;
+
+GLOBAL int outlets_num; /* number outlets: index for stream statistics*/
+STREAM *stat_streams;
+STATS *ord_stats;
+STATS stats_total;
+int order_max;
+
+GLOBAL struct History history; /* holds meta-data (title, comments,..) */
+GLOBAL struct Cell_head window;
+
+
+
+
Property changes on: grass-addons/raster/r.stream.stats/global.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/io.c
===================================================================
--- grass-addons/raster/r.stream.stats/io.c (rev 0)
+++ grass-addons/raster/r.stream.stats/io.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,155 @@
+#include <grass/glocale.h>
+#include "global.h"
+
+int open_raster(char *mapname)
+{
+ int fd = 0;
+ char *mapset;
+ struct Cell_head cellhd; /* it stores region information, and header information of rasters */
+
+ mapset = G_find_cell2(mapname, ""); /* checking if map exist */
+ if (mapset == NULL) {
+ G_fatal_error(_("Raster map <%s> not found"), mapname);
+ }
+
+ if (mapname != in_elev) {
+ 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);
+ }
+
+ if ((fd = G_open_cell_old(mapname, mapset)) < 0) /* file descriptor */
+ 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);
+
+ return fd;
+}
+
+int create_maps(void)
+{
+ int r, c;
+ int elev_type;
+ int in_dir_fd, in_stm_fd, in_dem_fd;
+ CELL *r_dirs, *r_streams;
+ FCELL *r_dem_f = NULL;
+ CELL *r_dem_c = NULL;
+ DCELL *r_dem_d = NULL;
+
+ in_dir_fd = open_raster(in_dirs);
+ in_stm_fd = open_raster(in_streams);
+ in_dem_fd = open_raster(in_elev);
+
+ r_dirs = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ r_streams = (CELL *) G_malloc(sizeof(CELL) * ncols);
+
+ elev_type = G_get_raster_map_type(in_dem_fd); /* type of DEM */
+ switch (elev_type) {
+ case CELL_TYPE:
+ r_dem_c = (CELL *) G_malloc(sizeof(CELL) * ncols);
+ break;
+ case FCELL_TYPE:
+ r_dem_f = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+ break;
+ case DCELL_TYPE:
+ r_dem_d = (DCELL *) G_malloc(sizeof(DCELL) * ncols);
+ break;
+ }
+
+ dirs = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+ streams = (CELL **) G_malloc(sizeof(CELL *) * nrows);
+ elevation = (FCELL **) G_malloc(sizeof(FCELL *) * nrows);
+
+
+ for (r = 0; r < nrows; ++r) {
+
+ /* 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);
+ G_close_cell(in_stm_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+
+ elevation[r] = (FCELL *) G_malloc(sizeof(FCELL) * ncols);
+
+ switch (elev_type) {
+ case CELL_TYPE: /* CELL */
+ if (G_get_c_raster_row(in_dem_fd, r_dem_c, r) < 0) {
+ G_close_cell(in_dem_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_dem_c[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_c[c];
+ }
+ } /* end for */
+ break; /* CELL */
+
+ case DCELL_TYPE: /* CELL */
+ if (G_get_d_raster_row(in_dem_fd, r_dem_d, r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+ for (c = 0; c < ncols; ++c) {
+ if (G_is_d_null_value(&r_dem_d[c])) {
+ G_set_f_null_value(&elevation[r][c], 1);
+ }
+ else {
+ elevation[r][c] = (FCELL) r_dem_d[c];
+ }
+ }
+ break; /* CELL */
+
+ case FCELL_TYPE: /* CELL */
+ if (G_get_f_raster_row(in_dem_fd, elevation[r], r) < 0) {
+ G_close_cell(in_dem_fd);
+ G_fatal_error(_("Unable to read raster maps at row <%d>"), r);
+ }
+ break; /* CELL */
+
+ } /* end switch */
+
+
+ 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])) {
+ streams[r][c] = 0;
+ }
+ else {
+ streams[r][c] = r_streams[c];
+ }
+ } /* end for c */
+ } /*end for r */
+
+ switch (elev_type) {
+ case CELL_TYPE:
+ G_free(r_dem_c);
+ break;
+ case DCELL_TYPE:
+ G_free(r_dem_d);
+ break;
+ case FCELL_TYPE:
+ break;
+ }
+
+ G_free(r_streams);
+ G_free(r_dirs);
+ G_close_cell(in_dir_fd);
+ G_close_cell(in_stm_fd);
+ G_close_cell(in_dem_fd);
+
+ return 0;
+} /* end create maps */
Property changes on: grass-addons/raster/r.stream.stats/io.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/main.c
===================================================================
--- grass-addons/raster/r.stream.stats/main.c (rev 0)
+++ grass-addons/raster/r.stream.stats/main.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,99 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.stats
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Calculate Horton's statistics according stream network and elevation map.
+ * Program calculates: Bifuarcation ratio, length ratio, area ratio,
+ * slope ratio and drainage density
+ * It uses r.stream.order stream map map, r.watershed direction map. and DEM
+ * Stram input map shall contains streams ordered according Strahler's or Horton's
+ * orders. It also can calculate Hack's laws as an option.
+ * If input stream comes from r.stream.exteract direction map
+ * from r.stream.extract dir map must be patched with that of r.watersheed.
+ *
+ * COPYRIGHT: (C) 2002,2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "global.h"
+
+/*
+ * main function
+ *
+ *
+ */
+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; /* options */
+ struct Flag *out_hack; /* flags */
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords =
+ _("stream, order, Strahler, Horton, Hack, statisctics");
+ module->description =
+ _("Calculate Horton's and optionaly Hack's statistics according user' input.");
+
+ in_stm_opt = G_define_option(); /* input stream mask file - optional */
+ in_stm_opt->key = "stream";
+ in_stm_opt->type = TYPE_STRING;
+ in_stm_opt->required = YES; /* for now; TO DO: is planned to be optional */
+ in_stm_opt->gisprompt = "old,cell,raster";
+ in_stm_opt->description = "Name of streams mask input map";
+
+ in_dir_opt = G_define_option(); /* input directon file */
+ in_dir_opt->key = "dir";
+ in_dir_opt->type = TYPE_STRING;
+ in_dir_opt->required = YES;
+ in_dir_opt->gisprompt = "old,cell,raster";
+ in_dir_opt->description = "Name of flow direction input map";
+
+ in_elev_opt = G_define_option(); /* input stream mask file - optional */
+ in_elev_opt->key = "dem";
+ in_elev_opt->type = TYPE_STRING;
+ in_elev_opt->required = YES;
+ in_elev_opt->gisprompt = "old,cell,raster";
+ in_elev_opt->description = "Name of elevation map";
+
+
+ /* Define the different flags */
+ /*
+ out_hack = G_define_flag();
+ out_hack->key = 'h';
+ out_hack->description = _("Calculate Hack's statistic instead of Horton's");
+ */
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ in_elev = in_elev_opt->answer;
+ /* hack = (out_hack->answer != 0); */
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ G_get_window(&window);
+
+ create_maps();
+ init_streams();
+ calculate_streams();
+ calculate_basins();
+ stats();
+ print_stats();
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/raster/r.stream.stats/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/print_stats.c
===================================================================
--- grass-addons/raster/r.stream.stats/print_stats.c (rev 0)
+++ grass-addons/raster/r.stream.stats/print_stats.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,91 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int print_stats(void)
+{
+
+ int i;
+
+ /* summary statistics */
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Summary:\n");
+ fprintf(stdout,
+ "Max order | Tot.N.str. | Tot.str.len. | Tot.area. | Dr.dens. | Str.freq. \n");
+ fprintf(stdout,
+ " (num) | (num) | (km) | (km2) | (km/km2) | (num/km2) \n");
+ fprintf(stdout, " %8d | %10d | %12.4f | %9.4f | %8.4f | %7.4f \n",
+ stats_total.order, stats_total.stream_num,
+ stats_total.sum_length / 1000, stats_total.sum_area / 1000000,
+ stats_total.drainage_density * 1000,
+ stats_total.stream_frequency * 1000000);
+
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Stream ratios with standard deviations:\n");
+ fprintf(stdout, " Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. \n");
+ fprintf(stdout, " %7.4f | %7.4f | %8.4f | %7.4f | %7.4f\n",
+ stats_total.bifur_ratio,
+ stats_total.length_ratio,
+ stats_total.area_ratio,
+ stats_total.slope_ratio, stats_total.gradient_ratio);
+
+ fprintf(stdout, " %7.4f | %7.4f | %8.4f | %7.4f | %7.4f\n",
+ stats_total.std_bifur_ratio,
+ stats_total.std_length_ratio,
+ stats_total.std_area_ratio,
+ stats_total.std_slope_ratio, stats_total.std_gradient_ratio);
+ fprintf(stdout, "\n");
+
+ /* base parameters */
+ fprintf(stdout,
+ "Order | Avg.len | Avg.ar | Avg.sl | Avg.grad. | Avg.el.dif\n");
+ fprintf(stdout,
+ " num | (km) | (km2) | (m/m) | (m/m) | (m) \n");
+ for (i = 1; i <= order_max; ++i) {
+ fprintf(stdout, "%5d | %7.4f | %8.4f | %7.4f | %10.4f | %7.4f\n",
+ ord_stats[i].order,
+ ord_stats[i].avg_length / 1000,
+ ord_stats[i].avg_area / 1000000,
+ ord_stats[i].avg_slope,
+ ord_stats[i].avg_gradient, ord_stats[i].avg_elev_diff);
+ }
+ fprintf(stdout, "\n");
+ /* std dev of base parameters */
+ fprintf(stdout,
+ "Order | Std.len | Std.ar | Std.sl | Std.grad. | Std.el.dif\n");
+ fprintf(stdout,
+ " num | (km) | (km2) | (m/m) | (m/m) | (m) \n");
+ for (i = 1; i <= order_max; ++i) {
+ fprintf(stdout, "%5d | %7.4f | %8.4f | %7.4f | %10.4f | %7.4f\n",
+ ord_stats[i].order,
+ ord_stats[i].std_length / 1000,
+ ord_stats[i].std_area / 1000000,
+ ord_stats[i].std_slope,
+ ord_stats[i].std_gradient, ord_stats[i].std_elev_diff);
+ }
+ /* sum statistics of orders */
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Order | N.streams | Tot.len (km) | Tot.area (km2)\n");
+ for (i = 1; i <= order_max; ++i) {
+ fprintf(stdout, "%5d | %9d | %12.4f | %7.4f\n",
+ ord_stats[i].order,
+ ord_stats[i].stream_num,
+ ord_stats[i].sum_length / 1000,
+ ord_stats[i].sum_area / 1000000);
+ }
+ /* order ratios */
+ fprintf(stdout, "\n");
+ fprintf(stdout,
+ "Order | Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. | d.dens. | str.freq.\n");
+ for (i = 1; i <= order_max; ++i) {
+ fprintf(stdout,
+ "%5d | %7.4f | %7.4f | %8.4f | %7.4f | %7.4f | %7.4f | %7.4f\n",
+ ord_stats[i].order, ord_stats[i].bifur_ratio,
+ ord_stats[i].length_ratio, ord_stats[i].area_ratio,
+ ord_stats[i].slope_ratio, ord_stats[i].gradient_ratio,
+ ord_stats[i].drainage_density * 1000,
+ ord_stats[i].stream_frequency * 1000000);
+ }
+}
Property changes on: grass-addons/raster/r.stream.stats/print_stats.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/raster/r.stream.stats/stats.c
===================================================================
--- grass-addons/raster/r.stream.stats/stats.c (rev 0)
+++ grass-addons/raster/r.stream.stats/stats.c 2009-09-19 11:15:35 UTC (rev 39247)
@@ -0,0 +1,431 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int init_streams(void)
+{
+ int d, i, j; /* d: direction, i: iteration */
+ int r, c;
+ int next_stream = -1, cur_stream;
+ int out_max = ncols + nrows;
+ POINT *outlets;
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ outlets = (POINT *) G_malloc((out_max) * sizeof(POINT));
+
+ 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 =
+ (POINT *) G_realloc(outlets, out_max * sizeof(POINT));
+ }
+ 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 */
+ }
+ 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];
+ if (cur_stream != next_stream) { /* is outlet or node! */
+ outlets[outlets_num].r = r;
+ outlets[outlets_num].c = c;
+ outlets_num++;
+ }
+ } /* end if streams */
+ } /* end for */
+ } /* end for */
+
+ stat_streams = (STREAM *) G_malloc((outlets_num) * sizeof(STREAM));
+ for (i = 0; i < outlets_num; ++i) {
+ stat_streams[i].r = outlets[i].r;
+ stat_streams[i].c = outlets[i].c;
+ stat_streams[i].index = i;
+ stat_streams[i].slope = 0.;
+ stat_streams[i].gradient = 0.;
+ stat_streams[i].length = 0.;
+ stat_streams[i].elev_diff = 0.;
+ stat_streams[i].elev_spring = 0.;
+ stat_streams[i].elev_outlet = elevation[outlets[i].r][outlets[i].c];
+ stat_streams[i].order = streams[outlets[i].r][outlets[i].c];
+ stat_streams[i].basin_area = 0.;
+ stat_streams[i].cell_num = 0;
+ }
+ G_free(outlets);
+ return 0;
+}
+
+int calculate_streams(void)
+{
+
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ int i, j, s; /* s - streams index */
+ int done = 1;
+ int r, c;
+ float cur_northing, cur_easting;
+ float next_northing, next_easting;
+ float diff_elev, cur_length;
+ int cur_stream_order;
+
+ G_begin_distance_calculations();
+
+ for (s = 0; s < outlets_num; ++s) {
+ r = stat_streams[s].r;
+ c = stat_streams[s].c;
+ done = 1;
+
+ while (done) {
+ done = 0;
+ cur_northing = window.north - (r + .5) * window.ns_res;
+ cur_easting = window.west + (c + .5) * window.ew_res;
+
+ stat_streams[s].cell_num++;
+ stat_streams[s].elev_spring = elevation[r][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; /* border */
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] ==
+ stat_streams[s].order &&
+ dirs[r + nextr[i]][c + nextc[i]] == j) {
+
+ next_northing =
+ window.north - (r + nextr[i] + .5) * window.ns_res;
+ next_easting =
+ window.west + (c + nextc[i] + .5) * window.ew_res;
+ cur_length =
+ G_distance(next_easting, next_northing, cur_easting,
+ cur_northing);
+
+ diff_elev =
+ elevation[r + nextr[i]][c + nextc[i]] -
+ elevation[r][c];
+ diff_elev = (diff_elev < 0) ? 0. : diff_elev; /* water cannot flow up */
+
+ stat_streams[s].length += cur_length;
+ stat_streams[s].slope += (diff_elev / cur_length);
+
+ r = r + nextr[i];
+ c = c + nextc[i];
+ done = 1;
+ break;
+ } /* end if */
+ } /* end for i */
+ } /* end while */
+ }
+ return 0;
+}
+
+int calculate_basins(void)
+{
+ int i;
+
+ G_begin_cell_area_calculations();
+ fifo_max = 4 * (nrows + ncols);
+ fifo_outlet = (POINT *) G_malloc((fifo_max + 1) * sizeof(POINT));
+
+ for (i = 0; i < outlets_num; ++i)
+ stat_streams[i].basin_area =
+ fill_basin(stat_streams[i].r, stat_streams[i].c);
+
+ G_free(fifo_outlet);
+}
+
+
+int fill_basin(int r, int c)
+{
+ int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+ int i, j;
+ float area;
+ POINT n_cell;
+
+ tail = 0;
+ head = -1;
+ area = G_area_of_cell_at_row(r);
+
+ while (tail != head) {
+ 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; /* border */
+ j = (i + 4) > 8 ? i - 4 : i + 4;
+ if (dirs[r + nextr[i]][c + nextc[i]] == j) { /* countributing cell */
+
+ area += G_area_of_cell_at_row(r);
+ n_cell.r = (r + nextr[i]);
+ n_cell.c = (c + nextc[i]);
+ fifo_insert(n_cell);
+ }
+ } /* end for i... */
+
+ n_cell = fifo_return_del();
+ r = n_cell.r;
+ c = n_cell.c;
+
+ }
+
+ return area;
+}
+
+/* fifo functions */
+int fifo_insert(POINT point)
+{
+ fifo_outlet[tail++] = point;
+ if (tail > fifo_max)
+ tail = 0;
+ return 0;
+}
+
+POINT fifo_return_del(void)
+{
+ if (head > fifo_max)
+ head = -1;
+ return fifo_outlet[++head];
+}
+
+int max_order(void)
+{
+ int i;
+ int max_order = 0;
+
+ for (i = 0; i < outlets_num; ++i) {
+ if (stat_streams[i].order > max_order)
+ max_order = stat_streams[i].order;
+ }
+ return max_order;
+}
+
+int stats(void)
+{
+ int i;
+ int num, ord_num;
+ float snum, ord_snum;
+ float tmp_num;
+
+ order_max = max_order();
+
+ ord_stats = (STATS *) G_malloc((order_max + 1) * sizeof(STATS));
+
+ stats_total.order = 0;
+ stats_total.stream_num = 0;
+ stats_total.sum_length = 0.;
+ stats_total.avg_length = 0.;
+ stats_total.std_length = 0.;
+ stats_total.avg_slope = 0.;
+ stats_total.std_slope = 0.;
+ stats_total.avg_gradient = 0.;
+ stats_total.std_gradient = 0.;
+ stats_total.sum_area = 0.;
+ stats_total.avg_area = 0.;
+ stats_total.std_area = 0.;
+ stats_total.avg_elev_diff = 0.;
+ stats_total.std_elev_diff = 0.;
+ stats_total.bifur_ratio = 0.;
+ stats_total.std_bifur_ratio = 0.;
+ stats_total.length_ratio = 0.;
+ stats_total.std_length_ratio = 0.;
+ stats_total.area_ratio = 0.;
+ stats_total.std_area_ratio = 0.;
+ stats_total.slope_ratio = 0.;
+ stats_total.std_slope_ratio = 0.;
+ stats_total.gradient_ratio = 0.;
+ stats_total.std_gradient_ratio = 0.;
+ stats_total.stream_frequency = 0.;
+ stats_total.drainage_density = 0.;
+
+ for (i = 0; i <= order_max; ++i) {
+ ord_stats[i].order = i;
+ ord_stats[i].stream_num = 0;
+ ord_stats[i].sum_length = 0.;
+ ord_stats[i].avg_length = 0.;
+ ord_stats[i].std_length = 0.;
+ ord_stats[i].avg_slope = 0.;
+ ord_stats[i].std_slope = 0.;
+ ord_stats[i].avg_gradient = 0.;
+ ord_stats[i].std_gradient = 0.;
+ ord_stats[i].sum_area = 0.;
+ ord_stats[i].avg_area = 0.;
+ ord_stats[i].std_area = 0.;
+ ord_stats[i].avg_elev_diff = 0.;
+ ord_stats[i].std_elev_diff = 0.;
+ ord_stats[i].bifur_ratio = 0.;
+ ord_stats[i].std_bifur_ratio = 0.;
+ ord_stats[i].length_ratio = 0.;
+ ord_stats[i].std_length_ratio = 0.;
+ ord_stats[i].area_ratio = 0.;
+ ord_stats[i].std_area_ratio = 0.;
+ ord_stats[i].slope_ratio = 0.;
+ ord_stats[i].std_slope_ratio = 0.;
+ ord_stats[i].gradient_ratio = 0.;
+ ord_stats[i].std_gradient_ratio = 0.;
+ ord_stats[i].stream_frequency = 0.;
+ ord_stats[i].drainage_density = 0.;
+ }
+ for (i = 0; i < outlets_num; ++i) { /* recalculate and unify */
+ stat_streams[i].elev_diff =
+ stat_streams[i].elev_spring - stat_streams[i].elev_outlet;
+ tmp_num =
+ ((stat_streams[i].cell_num - 1) <
+ 1) ? 1 : stat_streams[i].cell_num - 1;
+ stat_streams[i].slope /= tmp_num;
+ stat_streams[i].gradient =
+ stat_streams[i].elev_diff / stat_streams[i].length;
+
+ /* calculation */
+ ord_stats[stat_streams[i].order].stream_num++;
+ ord_stats[stat_streams[i].order].sum_length += stat_streams[i].length;
+ ord_stats[stat_streams[i].order].std_length +=
+ (stat_streams[i].length * stat_streams[i].length);
+ ord_stats[stat_streams[i].order].avg_slope += stat_streams[i].slope;
+ ord_stats[stat_streams[i].order].std_slope +=
+ (stat_streams[i].slope * stat_streams[i].slope);
+ ord_stats[stat_streams[i].order].avg_gradient +=
+ stat_streams[i].gradient;
+ ord_stats[stat_streams[i].order].std_gradient +=
+ (stat_streams[i].gradient * stat_streams[i].gradient);
+ ord_stats[stat_streams[i].order].sum_area +=
+ stat_streams[i].basin_area;
+ ord_stats[stat_streams[i].order].std_area +=
+ (stat_streams[i].basin_area * stat_streams[i].basin_area);
+ ord_stats[stat_streams[i].order].avg_elev_diff +=
+ stat_streams[i].elev_diff;
+ ord_stats[stat_streams[i].order].std_elev_diff +=
+ (stat_streams[i].elev_diff * stat_streams[i].elev_diff);
+ }
+
+ for (i = 1; i <= order_max; ++i) {
+
+ num = ord_stats[i].stream_num;
+ snum = (ord_stats[i].stream_num > 1) ?
+ ((float)ord_stats[i].stream_num) / (ord_stats[i].stream_num -
+ 1) : 0.;
+
+ ord_stats[i].avg_length = ord_stats[i].sum_length / num;
+ ord_stats[i].avg_slope = ord_stats[i].avg_slope / num;
+ ord_stats[i].avg_gradient = ord_stats[i].avg_gradient / num;
+ ord_stats[i].avg_area = ord_stats[i].sum_area / num;
+ ord_stats[i].avg_elev_diff = ord_stats[i].avg_elev_diff / num;
+
+ ord_stats[i].std_length = sqrt((ord_stats[i].std_length / num -
+ (ord_stats[i].avg_length *
+ ord_stats[i].avg_length)) * snum);
+
+ ord_stats[i].std_slope = sqrt((ord_stats[i].std_slope / num -
+ (ord_stats[i].avg_slope *
+ ord_stats[i].avg_slope)) * snum);
+
+ ord_stats[i].std_gradient = sqrt((ord_stats[i].std_gradient / num -
+ (ord_stats[i].avg_gradient *
+ ord_stats[i].avg_gradient)) *
+ snum);
+
+ ord_stats[i].std_area = sqrt((ord_stats[i].std_area / num -
+ (ord_stats[i].avg_area *
+ ord_stats[i].avg_area)) * snum);
+
+ ord_stats[i].std_elev_diff = sqrt((ord_stats[i].std_elev_diff / num -
+ (ord_stats[i].avg_elev_diff *
+ ord_stats[i].avg_elev_diff)) *
+ snum);
+
+ ord_stats[i - 1].bifur_ratio =
+ ord_stats[i - 1].stream_num / (float)ord_stats[i].stream_num;
+ ord_stats[i - 1].length_ratio =
+ ord_stats[i - 1].avg_length / ord_stats[i].avg_length;
+ ord_stats[i - 1].area_ratio =
+ ord_stats[i - 1].avg_area / ord_stats[i].avg_area;
+ ord_stats[i - 1].slope_ratio =
+ ord_stats[i - 1].avg_slope / ord_stats[i].avg_slope;
+ ord_stats[i - 1].gradient_ratio =
+ ord_stats[i - 1].avg_gradient / ord_stats[i].avg_gradient;
+ ord_stats[i].stream_frequency =
+ ord_stats[i].stream_num / ord_stats[i].sum_area;
+ ord_stats[i].drainage_density =
+ ord_stats[i].sum_length / ord_stats[i].sum_area;
+
+ /* total */
+ stats_total.stream_num += ord_stats[i].stream_num;
+ stats_total.sum_length += ord_stats[i].sum_length;
+
+ stats_total.bifur_ratio += ord_stats[i - 1].bifur_ratio;
+ stats_total.length_ratio += ord_stats[i - 1].length_ratio;
+ stats_total.area_ratio += ord_stats[i - 1].area_ratio;
+ stats_total.slope_ratio += ord_stats[i - 1].slope_ratio;
+ stats_total.gradient_ratio += ord_stats[i - 1].gradient_ratio;
+
+ stats_total.std_bifur_ratio +=
+ (ord_stats[i - 1].bifur_ratio * ord_stats[i - 1].bifur_ratio);
+ stats_total.std_length_ratio +=
+ (ord_stats[i - 1].length_ratio * ord_stats[i - 1].length_ratio);
+ stats_total.std_area_ratio +=
+ (ord_stats[i - 1].area_ratio * ord_stats[i - 1].area_ratio);
+ stats_total.std_slope_ratio +=
+ (ord_stats[i - 1].slope_ratio * ord_stats[i - 1].slope_ratio);
+ stats_total.std_gradient_ratio +=
+ (ord_stats[i - 1].gradient_ratio *
+ ord_stats[i - 1].gradient_ratio);
+
+ } /* end for ... orders */
+ ord_num = order_max - 1;
+ ord_snum = (ord_num == 1) ? 0 : (float)ord_num / (ord_num - 1);
+
+ stats_total.order = order_max;
+ stats_total.sum_area = ord_stats[order_max].sum_area;
+ stats_total.sum_length = stats_total.sum_length;
+
+ stats_total.bifur_ratio = stats_total.bifur_ratio / ord_num;
+ stats_total.length_ratio = stats_total.length_ratio / ord_num;
+ stats_total.area_ratio = stats_total.area_ratio / ord_num;
+ stats_total.slope_ratio = stats_total.slope_ratio / ord_num;
+ stats_total.gradient_ratio = stats_total.gradient_ratio / ord_num;
+
+
+ stats_total.std_bifur_ratio =
+ sqrt((stats_total.std_bifur_ratio / ord_num -
+ (stats_total.bifur_ratio * stats_total.bifur_ratio)) *
+ ord_snum);
+
+ stats_total.std_length_ratio =
+ sqrt((stats_total.std_length_ratio / ord_num -
+ (stats_total.length_ratio * stats_total.length_ratio)) *
+ ord_snum);
+
+ stats_total.std_area_ratio = sqrt((stats_total.std_area_ratio / ord_num -
+ (stats_total.area_ratio *
+ stats_total.area_ratio)) * ord_snum);
+
+ stats_total.std_slope_ratio =
+ sqrt((stats_total.std_slope_ratio / ord_num -
+ (stats_total.slope_ratio * stats_total.slope_ratio)) *
+ ord_snum);
+
+ stats_total.std_gradient_ratio =
+ sqrt((stats_total.std_gradient_ratio / ord_num -
+ (stats_total.gradient_ratio * stats_total.gradient_ratio)) *
+ ord_snum);
+
+ stats_total.stream_frequency =
+ stats_total.stream_num / stats_total.sum_area;
+ stats_total.drainage_density =
+ stats_total.sum_length / stats_total.sum_area;
+
+}
Property changes on: grass-addons/raster/r.stream.stats/stats.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list