[GRASS-SVN] r39696 - grass-addons/raster/r.stream.del
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 9 09:33:09 EST 2009
Author: jarekj71
Date: 2009-11-09 09:33:09 -0500 (Mon, 09 Nov 2009)
New Revision: 39696
Added:
grass-addons/raster/r.stream.del/Makefile
grass-addons/raster/r.stream.del/delete.c
grass-addons/raster/r.stream.del/description.html
grass-addons/raster/r.stream.del/global.h
grass-addons/raster/r.stream.del/io.c
grass-addons/raster/r.stream.del/main.c
Log:
new module
Added: grass-addons/raster/r.stream.del/Makefile
===================================================================
--- grass-addons/raster/r.stream.del/Makefile (rev 0)
+++ grass-addons/raster/r.stream.del/Makefile 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,13 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.del
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/raster/r.stream.del/delete.c
===================================================================
--- grass-addons/raster/r.stream.del/delete.c (rev 0)
+++ grass-addons/raster/r.stream.del/delete.c 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,151 @@
+#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..."));
+
+ return trib;
+} /* end trib_num */
+
+
+int find_springs(int max_link)
+{
+
+ int d, i, j; /* d: direction, i: iteration */
+ int r, c;
+ int trib_num, trib = 0;
+ int cur_stream;
+ int springs_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 springs..."));
+
+ springs = (SPRING *) G_malloc(max_link * sizeof(SPRING));
+
+ for (r = 0; r < nrows; ++r) {
+ for (c = 0; c < ncols; ++c) {
+ if (streams[r][c] > 0) {
+ trib_num = trib_nums(r, c);
+ if (trib_num == 0) { /* adding springs */
+
+ if (springs_num > (max_link - 1))
+ G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
+
+ springs[springs_num].r=r;
+ springs[springs_num].c=c;
+ springs[springs_num++].val = streams[r][c]; /* collecting springs */
+ }
+
+
+ } /* end if streams */
+ } /* c */
+ } /* r */
+ return springs_num;
+}
+
+int delete_join (int springs_num)
+{
+
+ int i,j, d;
+ int r,c,nr,nc;
+ 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 cell_num=1;
+ int cur_stream,tmp_stream;
+ int trib_num;
+
+G_message(_("Deleting short streams..."));
+
+for (j=0;j<springs_num;++j) {
+
+r=springs[j].r;
+c=springs[j].c;
+cell_num=1;
+
+ while (cell_num<threshold && dirs[r][c]>0) {
+ nr=r+nextr[dirs[r][c]];
+ nc=c+nextc[dirs[r][c]];
+
+ if (streams[nr][nc]==streams[r][c]) {
+ cell_num++;
+ r=nr;
+ c=nc;
+ } else {
+ break; /* if stream isn't longer due to joining in previous interation */
+ }
+
+ } /*end while */
+
+ if(cell_num<threshold) {
+
+
+ r=springs[j].r;
+ c=springs[j].c;
+ cur_stream=streams[r][c];
+ streams[r][c]=0;
+ nr=r+nextr[dirs[r][c]];
+ nc=c+nextc[dirs[r][c]];
+
+ while (dirs[r][c]>0 && streams[nr][nc]==cur_stream)
+ { /* only when there is more than one cell in stream */
+ r=nr;
+ c=nc;
+ streams[r][c]=0;
+ nr=r+nextr[dirs[r][c]];
+ nc=c+nextc[dirs[r][c]];
+ }
+
+ trib_num = trib_nums(nr, nc); // check if there are no addational tributuaries
+
+ if (trib_num != 1) {
+ continue; //still is node (>1), or border (==0) finish and go to the next spring
+ } else {
+ r=nr;
+ c=nc;
+ cur_stream=0;
+ 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 */
+ d = (i + 4) > 8 ? i - 4 : i + 4;
+ if (streams[r + nextr[i]][c + nextc[i]] > 0 && dirs[r + nextr[i]][c + nextc[i]] == d)
+ cur_stream=streams[r + nextr[i]][c + nextc[i]];
+ }
+ if(cur_stream==0)
+ G_fatal_error(_("Problem with joining streams"));
+
+ tmp_stream=streams[r][c];
+ nr=r+nextr[dirs[r][c]];
+ nc=c+nextc[dirs[r][c]];
+ streams[r][c]=cur_stream;
+
+ while (dirs[r][c]>0 && streams[nr][nc]==tmp_stream) {
+ r=nr;
+ c=nc;
+ nr=r+nextr[dirs[r][c]];
+ nc=c+nextc[dirs[r][c]];
+ streams[r][c]=cur_stream;
+ }
+ } /*end if_elese */
+ } /* end if cell_num*/
+} /*end for */
+return 0;
+} /*end function */
Added: grass-addons/raster/r.stream.del/description.html
===================================================================
--- grass-addons/raster/r.stream.del/description.html (rev 0)
+++ grass-addons/raster/r.stream.del/description.html 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,76 @@
+<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 reduction will be performed produced by <a href="r.watershed.html">r.watershed</a>. <a href="r.stream.extract.html">r.stream.extract</a> has deleting short stream procedure build-in. Because streams network produced by r.watershed with SFD and MFD may differ it is required to use both stream and direction map produced by the same operation. 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.watershedRegion 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>
+<DT><b>threshold</b></DT>
+<DD>Integer vaule indicating the minimum number of cell in first-order streams (sensu Strhahler) to leeve it in the network.</DD>
+<h2>OUTPUTS</h2>
+<P>The module produces new stream map where streams shortrer than threshold are deleted. To keep file consistent, when short stream is deleted the next stream segment is joined with previous one by asigning its category.
+
+<PRE>
+BEFORE:
+
+ |2
+___1__|____3____
+
+AFTER:
+
+___1______1____
+
+</PRE>
+</P>
+</DL>
+
+<h2>DESCRIPTION</h2>
+<P>
+Module r.stream.del is prepared to removing short streams from stream network. The short streams are undesired outcome of networks delineated from DEMS in stream ordering and basin delinating. In stream ordering very short stream may increase network order, when two short stream of first order will increase order of next streams. The second situation is connected with network delineated with more complex criteria than only accumulation threeshold like Montgomery's method or weighted method (<a href="r.stream.extract.html">r.stream.extract</a>) It may produce a scope of very short streams on vallyes slopes. Such streams may lead to wrong results in stream ordering and basin delination so removing short streams to simplify the network. It also may heve great impact on Horton's statistics: <a href="r.stream.stats.html">r.stream.stats</a>.
+
+
+<h2>NOTES</h2>
+<P>
+<a href="r.stream.extract.html">r.stream.extract</a> has option to delete shot stream, so r.stream.del is preperad only for use with r.watershed output network. So this module will be probably depreciated in the nearest future.
+Module can work only if direction map, stream map and region map has same settings. It is also required that stream map and direction map come from the same source. For lots of reason this limitation probably cannot be omitted. this means if stream map comes from r.stream.extract also direction map from r.stream.extract must be used. If stream network was generated with MFD method also MFD direction map must be used.
+
+<h2>EXAMPLES</h2>
+
+<P>
+Strahler stream ordering before and after removing short streams
+<CODE>
+<PRE>
+r.watershed elevation=elev_ned_30m at PERMANENT threshold=100 d8cut=infinity mexp=2 stream_rast=stream_mont direction=stream_dir
+r.stream.order stream=stream_mont dir=stream_dir strahler=streahler_before
+
+r.stream.del stream=stream_mont dir=stream_dir threshold=10 reduced=stream_out
+r.stream.order stream=stream_out dir=stream_dir strahler=streahler_after
+
+
+</PRE>
+</CODE>
+<P>
+
+
+<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.stream.stats.html">r.stream.stats</a>,
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
+</body>
+</html>
Added: grass-addons/raster/r.stream.del/global.h
===================================================================
--- grass-addons/raster/r.stream.del/global.h (rev 0)
+++ grass-addons/raster/r.stream.del/global.h 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,53 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+
+
+ /* define */
+
+#define SPRING struct sprs
+SPRING {
+ int r, c;
+ int val;
+ };
+
+ /* functions.c */
+
+/* io.c */
+int create_maps(void);
+int max_link(void);
+int write_map(void);
+int set_null(void);
+
+/* delete */
+int find_springs(int max_link);
+int delete_join(int springs_num);
+
+ /* 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 *out_streams;
+GLOBAL int zeros; /* flags */
+GLOBAL int threshold;
+
+GLOBAL CELL **dirs, **streams; /* matrix with input data streams is used as output data*/
+
+GLOBAL int nrows, ncols;
+
+SPRING *springs;
+
+GLOBAL struct History history; /* holds meta-data (title, comments,..) */
+
+
+
+
+
Added: grass-addons/raster/r.stream.del/io.c
===================================================================
--- grass-addons/raster/r.stream.del/io.c (rev 0)
+++ grass-addons/raster/r.stream.del/io.c 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,143 @@
+#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_map(void)
+{
+ int r;
+ int fdo = 0;
+
+ if ((fdo = G_open_raster_new(out_streams, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), out_streams);
+
+ for (r = 0; r < nrows; ++r)
+ G_put_c_raster_row(fdo, streams[r]);
+ G_close_cell(fdo);
+ G_short_history(out_streams, "raster", &history);
+ G_write_history(out_streams, &history);
+ G_message(_("%s Done!"), out_streams);
+
+ 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;
+}
Added: grass-addons/raster/r.stream.del/main.c
===================================================================
--- grass-addons/raster/r.stream.del/main.c (rev 0)
+++ grass-addons/raster/r.stream.del/main.c 2009-11-09 14:33:09 UTC (rev 39696)
@@ -0,0 +1,114 @@
+
+/****************************************************************************
+ *
+ * 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, *in_threshold, *out_opt; /* options */
+ struct Flag *out_back; /* flags */
+
+ int num_spring, 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";
+
+ in_threshold = G_define_option();
+ in_threshold->key = "threshold";
+ in_threshold->label = _("Minimum number of cell in stream");
+ in_threshold->description = _("Must be > 0");
+ in_threshold->required = YES;
+ in_threshold->type = TYPE_INTEGER;
+
+ out_opt = G_define_option();
+ out_opt->key = "reduced";
+ out_opt->type = TYPE_STRING;
+ out_opt->required = YES;
+ out_opt->answer = NULL;
+ out_opt->gisprompt = "new,cell,raster";
+ out_opt->description = "Output reduced stream map";
+
+
+ /* Define the different 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);
+
+ /* stores input options to variables */
+ in_dirs = in_dir_opt->answer;
+ in_streams = in_stm_opt->answer;
+ out_streams = out_opt->answer;
+ zeros = (out_back->answer != 0);
+
+ threshold = atoi(in_threshold->answer);
+ if (threshold <= 0)
+ G_fatal_error(_("Threshold must be > 0"));
+
+
+ if (G_legal_filename(out_streams) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out_streams);
+
+ G_get_window(&window);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ create_maps();
+ link_max = max_link();
+ num_spring=find_springs(link_max);
+ delete_join(num_spring);
+ if (!zeros)
+ set_null();
+ write_map();
+
+ exit(EXIT_SUCCESS);
+}
More information about the grass-commit
mailing list