[GRASS-SVN] r34568 - in grass/trunk/raster: . r.watershed
r.watershed2 r.watershed2/front r.watershed2/ram
r.watershed2/seg r.watershed2/shed
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Nov 27 18:00:38 EST 2008
Author: neteler
Date: 2008-11-27 18:00:37 -0500 (Thu, 27 Nov 2008)
New Revision: 34568
Added:
grass/trunk/raster/r.watershed/DEPRECATED
grass/trunk/raster/r.watershed2/
grass/trunk/raster/r.watershed2/Makefile
grass/trunk/raster/r.watershed2/README
grass/trunk/raster/r.watershed2/front/
grass/trunk/raster/r.watershed2/front/Makefile
grass/trunk/raster/r.watershed2/front/main.c
grass/trunk/raster/r.watershed2/front/r.watershed.html
grass/trunk/raster/r.watershed2/ram/
grass/trunk/raster/r.watershed2/ram/Gwater.h
grass/trunk/raster/r.watershed2/ram/Makefile
grass/trunk/raster/r.watershed2/ram/close_maps.c
grass/trunk/raster/r.watershed2/ram/close_maps2.c
grass/trunk/raster/r.watershed2/ram/def_basin.c
grass/trunk/raster/r.watershed2/ram/do_astar.c
grass/trunk/raster/r.watershed2/ram/do_astar.h
grass/trunk/raster/r.watershed2/ram/do_cum.c
grass/trunk/raster/r.watershed2/ram/find_pour.c
grass/trunk/raster/r.watershed2/ram/flag.h
grass/trunk/raster/r.watershed2/ram/flag_clr_all.c
grass/trunk/raster/r.watershed2/ram/flag_create.c
grass/trunk/raster/r.watershed2/ram/flag_destroy.c
grass/trunk/raster/r.watershed2/ram/flag_get.c
grass/trunk/raster/r.watershed2/ram/flag_set.c
grass/trunk/raster/r.watershed2/ram/flag_unset.c
grass/trunk/raster/r.watershed2/ram/haf_side.c
grass/trunk/raster/r.watershed2/ram/init_vars.c
grass/trunk/raster/r.watershed2/ram/main.c
grass/trunk/raster/r.watershed2/ram/no_stream.c
grass/trunk/raster/r.watershed2/ram/over_cells.c
grass/trunk/raster/r.watershed2/ram/ramseg.c
grass/trunk/raster/r.watershed2/ram/ramseg.h
grass/trunk/raster/r.watershed2/ram/sg_factor.c
grass/trunk/raster/r.watershed2/ram/slope_len.c
grass/trunk/raster/r.watershed2/ram/split_str.c
grass/trunk/raster/r.watershed2/ram/usage.c
grass/trunk/raster/r.watershed2/seg/
grass/trunk/raster/r.watershed2/seg/Gwater.h
grass/trunk/raster/r.watershed2/seg/Makefile
grass/trunk/raster/r.watershed2/seg/bseg_close.c
grass/trunk/raster/r.watershed2/seg/bseg_get.c
grass/trunk/raster/r.watershed2/seg/bseg_open.c
grass/trunk/raster/r.watershed2/seg/bseg_put.c
grass/trunk/raster/r.watershed2/seg/bseg_read.c
grass/trunk/raster/r.watershed2/seg/bseg_write.c
grass/trunk/raster/r.watershed2/seg/close_maps.c
grass/trunk/raster/r.watershed2/seg/close_maps2.c
grass/trunk/raster/r.watershed2/seg/cseg.h
grass/trunk/raster/r.watershed2/seg/cseg_close.c
grass/trunk/raster/r.watershed2/seg/cseg_get.c
grass/trunk/raster/r.watershed2/seg/cseg_open.c
grass/trunk/raster/r.watershed2/seg/cseg_put.c
grass/trunk/raster/r.watershed2/seg/cseg_read.c
grass/trunk/raster/r.watershed2/seg/cseg_write.c
grass/trunk/raster/r.watershed2/seg/def_basin.c
grass/trunk/raster/r.watershed2/seg/do_astar.c
grass/trunk/raster/r.watershed2/seg/do_astar.h
grass/trunk/raster/r.watershed2/seg/do_cum.c
grass/trunk/raster/r.watershed2/seg/dseg_close.c
grass/trunk/raster/r.watershed2/seg/dseg_get.c
grass/trunk/raster/r.watershed2/seg/dseg_open.c
grass/trunk/raster/r.watershed2/seg/dseg_put.c
grass/trunk/raster/r.watershed2/seg/dseg_read.c
grass/trunk/raster/r.watershed2/seg/dseg_write.c
grass/trunk/raster/r.watershed2/seg/find_pour.c
grass/trunk/raster/r.watershed2/seg/haf_side.c
grass/trunk/raster/r.watershed2/seg/init_vars.c
grass/trunk/raster/r.watershed2/seg/main.c
grass/trunk/raster/r.watershed2/seg/no_stream.c
grass/trunk/raster/r.watershed2/seg/over_cells.c
grass/trunk/raster/r.watershed2/seg/sg_factor.c
grass/trunk/raster/r.watershed2/seg/slope_len.c
grass/trunk/raster/r.watershed2/seg/split_str.c
grass/trunk/raster/r.watershed2/seg/sseg_close.c
grass/trunk/raster/r.watershed2/seg/sseg_get.c
grass/trunk/raster/r.watershed2/seg/sseg_open.c
grass/trunk/raster/r.watershed2/seg/sseg_put.c
grass/trunk/raster/r.watershed2/seg/usage.c
grass/trunk/raster/r.watershed2/shed/
grass/trunk/raster/r.watershed2/shed/Makefile
grass/trunk/raster/r.watershed2/shed/accum_down.c
grass/trunk/raster/r.watershed2/shed/basin_maps.c
grass/trunk/raster/r.watershed2/shed/com_line.c
grass/trunk/raster/r.watershed2/shed/file_in.c
grass/trunk/raster/r.watershed2/shed/free.c
grass/trunk/raster/r.watershed2/shed/insert_cat.c
grass/trunk/raster/r.watershed2/shed/intro.c
grass/trunk/raster/r.watershed2/shed/main.c
grass/trunk/raster/r.watershed2/shed/print.c
grass/trunk/raster/r.watershed2/shed/read.c
grass/trunk/raster/r.watershed2/shed/valid.c
grass/trunk/raster/r.watershed2/shed/watershed.h
Modified:
grass/trunk/raster/Makefile
Log:
moved here from GRASS Addons; deactivated r.watershed
Modified: grass/trunk/raster/Makefile
===================================================================
--- grass/trunk/raster/Makefile 2008-11-27 23:00:35 UTC (rev 34567)
+++ grass/trunk/raster/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -114,7 +114,7 @@
r.volume \
r.walk \
r.water.outlet \
- r.watershed \
+ r.watershed2 \
r.what \
r.what.color
Added: grass/trunk/raster/r.watershed/DEPRECATED
===================================================================
--- grass/trunk/raster/r.watershed/DEPRECATED (rev 0)
+++ grass/trunk/raster/r.watershed/DEPRECATED 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,2 @@
+no longer maintained.
+Continued in ../r.watershed2/
Added: grass/trunk/raster/r.watershed2/Makefile
===================================================================
--- grass/trunk/raster/r.watershed2/Makefile (rev 0)
+++ grass/trunk/raster/r.watershed2/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../..
+
+SUBDIRS = \
+ front \
+ ram \
+ seg
+# inter unused:
+# shed
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: parsubdirs
+
+clean: cleansubdirs
+
Property changes on: grass/trunk/raster/r.watershed2/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/README
===================================================================
--- grass/trunk/raster/r.watershed2/README (rev 0)
+++ grass/trunk/raster/r.watershed2/README 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,50 @@
+
+ fast version of r.watershed
+ ===========================
+
+ Implemented is a new sorting algorithm for both the ram and the
+ segmented version, resulting in significant speed improvements,
+ with identical results.
+
+ The following procedure assumes that GRASS was compiled from source
+ and that the configured source tree is available.
+ To compile and install the test version of r.watershed for GRASS,
+
+ 1) copy the directory <r.watershed2> into the directory
+ <raster> of the GRASS souorce tree
+
+ 2) edit <Makefile> in directory <raster>
+ add the line
+ r.watershed2 \
+ e.g. after
+ r.watershed \
+
+ 3) change to directory <r.watershed2>
+ run make
+
+ 4) optionally change to GRASS source directory, something like
+ grass-6.4.svn
+ run as su
+ make install
+
+ Testing
+ =======
+
+ The new test version can be invoked with r.watershed2.
+ The new sorting algorithm is implemented for both the ram and the
+ segmented version.
+ The new segmented version is still much slower than the new ram
+ version, but the new segmented version is already much faster than
+ the original ram version. A new option is available for the segmented
+ mode specifying the maximum amount of memory to be used in MB.
+ Set up a watershed analysis with the original version, run it.
+ The command will have been r.watershed <options>
+ To test, use the same input options, but different output maps
+ and run r.watershed2 <test options>
+ There should be no differences in the outputs of the original and
+ the new version.
+ If you find critical differences between the original and the fast version,
+ please let me know at
+ markus.metz.giswork at gmail.com
+
+
Added: grass/trunk/raster/r.watershed2/front/Makefile
===================================================================
--- grass/trunk/raster/r.watershed2/front/Makefile (rev 0)
+++ grass/trunk/raster/r.watershed2/front/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.watershed
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass/trunk/raster/r.watershed2/front/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/front/main.c
===================================================================
--- grass/trunk/raster/r.watershed2/front/main.c (rev 0)
+++ grass/trunk/raster/r.watershed2/front/main.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,409 @@
+
+/****************************************************************************
+ *
+ * MODULE: front end
+ * AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
+ * Brad Douglas <rez touchofmadness.com>, Hamish Bowman <hamish_nospam yahoo.com>
+ * PURPOSE: Watershed determination
+ * COPYRIGHT: (C) 1999-2006 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.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+int write_hist(char *, char *, char *, int);
+
+int main(int argc, char *argv[])
+{
+ int err, ret;
+ char command[512];
+ struct Option *opt1;
+ struct Option *opt2;
+ struct Option *opt3;
+ struct Option *opt4;
+ struct Option *opt5;
+ struct Option *opt6;
+ struct Option *opt7;
+ struct Option *opt8;
+ struct Option *opt9;
+ struct Option *opt10;
+ struct Option *opt11;
+ struct Option *opt12;
+ struct Option *opt13;
+ struct Option *opt14;
+ struct Option *opt15;
+ struct Option *opt16;
+ struct Flag *flag1;
+ struct Flag *flag2;
+ struct GModule *module;
+
+ G_gisinit(argv[0]);
+
+ /* Set description */
+ module = G_define_module();
+ module->keywords = _("raster");
+ module->description = _("Watershed basin analysis program.");
+
+ opt1 = G_define_standard_option(G_OPT_R_ELEV);
+ opt1->description =
+ _("Input map: elevation on which entire analysis is based");
+ opt1->guisection = _("Input_options");
+
+ opt2 = G_define_option();
+ opt2->key = "depression";
+ opt2->description = _("Input map: locations of real depressions");
+ opt2->required = NO;
+ opt2->type = TYPE_STRING;
+ opt2->gisprompt = "old,cell,raster";
+ opt2->guisection = _("Input_options");
+
+ opt3 = G_define_option();
+ opt3->key = "flow";
+ opt3->description = _("Input map: amount of overland flow per cell");
+ opt3->required = NO;
+ opt3->type = TYPE_STRING;
+ opt3->gisprompt = "old,cell,raster";
+ opt3->guisection = _("Input_options");
+
+ opt4 = G_define_option();
+ opt4->key = "disturbed_land";
+ opt4->description =
+ _("Input map or value: percent of disturbed land, for USLE");
+ opt4->required = NO;
+ opt4->type = TYPE_STRING;
+ opt4->gisprompt = "old,cell,raster";
+ opt4->guisection = _("Input_options");
+
+ opt5 = G_define_option();
+ opt5->key = "blocking";
+ opt5->description =
+ _("Input map: terrain blocking overland surface flow, for USLE");
+ opt5->required = NO;
+ opt5->type = TYPE_STRING;
+ opt5->gisprompt = "old,cell,raster";
+ opt5->guisection = _("Input_options");
+
+ opt6 = G_define_option();
+ opt6->key = "threshold";
+ opt6->description =
+ _("Input value: minimum size of exterior watershed basin");
+ opt6->required = NO;
+ opt6->type = TYPE_INTEGER;
+ opt6->gisprompt = "new,cell,raster";
+ opt6->guisection = _("Input_options");
+
+ opt7 = G_define_option();
+ opt7->key = "max_slope_length";
+ opt7->description =
+ _("Input value: maximum length of surface flow, for USLE");
+ opt7->required = NO;
+ opt7->type = TYPE_DOUBLE;
+ opt7->gisprompt = "new,cell,raster";
+ opt7->guisection = _("Input_options");
+
+ opt8 = G_define_option();
+ opt8->key = "accumulation";
+ opt8->description =
+ _("Output map: number of cells that drain through each cell");
+ opt8->required = NO;
+ opt8->type = TYPE_STRING;
+ opt8->gisprompt = "new,cell,raster";
+ opt8->guisection = _("Output_options");
+
+ opt9 = G_define_option();
+ opt9->key = "drainage";
+ opt9->description = _("Output map: drainage direction");
+ opt9->required = NO;
+ opt9->type = TYPE_STRING;
+ opt9->gisprompt = "new,cell,raster";
+ opt9->guisection = _("Output_options");
+
+ opt10 = G_define_option();
+ opt10->key = "basin";
+ opt10->description =
+ _("Output map: unique label for each watershed basin");
+ opt10->required = NO;
+ opt10->type = TYPE_STRING;
+ opt10->gisprompt = "new,cell,raster";
+ opt10->guisection = _("Output_options");
+
+ opt11 = G_define_option();
+ opt11->key = "stream";
+ opt11->description = _("Output map: stream segments");
+ opt11->required = NO;
+ opt11->type = TYPE_STRING;
+ opt11->gisprompt = "new,cell,raster";
+ opt11->guisection = _("Output_options");
+
+ opt12 = G_define_option();
+ opt12->key = "half_basin";
+ opt12->description =
+ _("Output map: each half-basin is given a unique value");
+ opt12->required = NO;
+ opt12->type = TYPE_STRING;
+ opt12->gisprompt = "new,cell,raster";
+ opt12->guisection = _("Output_options");
+
+ opt13 = G_define_option();
+ opt13->key = "visual";
+ opt13->description =
+ _("Output map: useful for visual display of results");
+ opt13->required = NO;
+ opt13->type = TYPE_STRING;
+ opt13->gisprompt = "new,cell,raster";
+ opt13->guisection = _("Output_options");
+
+ opt14 = G_define_option();
+ opt14->key = "length_slope";
+ opt14->description =
+ _("Output map: slope length and steepness (LS) factor for USLE");
+ opt14->required = NO;
+ opt14->type = TYPE_STRING;
+ opt14->gisprompt = "new,cell,raster";
+ opt14->guisection = _("Output_options");
+
+ opt15 = G_define_option();
+ opt15->key = "slope_steepness";
+ opt15->description = _("Output map: slope steepness (S) factor for USLE");
+ opt15->required = NO;
+ opt15->type = TYPE_STRING;
+ opt15->gisprompt = "new,cell,raster";
+ opt15->guisection = _("Output_options");
+
+ opt16 = G_define_option() ;
+ opt16->key = "memory";
+ opt16->type = TYPE_INTEGER;
+ opt16->required = NO;
+ opt16->answer = "300"; /* 300MB default value */
+ opt16->description = _("Maximum memory to be used with -m flag (in MB)");
+
+
+ flag1 = G_define_flag();
+ flag1->key = '4';
+ flag1->description =
+ _("Allow only horizontal and vertical flow of water");
+
+ flag2 = G_define_flag();
+ flag2->key = 'm';
+ flag2->description =
+ _("Enable disk swap memory option: Operation is slow");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* Check option combinations */
+
+ /* Check for some output map */
+ if ((opt8->answer == NULL)
+ && (opt9->answer == NULL)
+ && (opt10->answer == NULL)
+ && (opt11->answer == NULL)
+ && (opt12->answer == NULL)
+ && (opt13->answer == NULL)
+ && (opt14->answer == NULL)
+ && (opt15->answer == NULL)) {
+ G_fatal_error(_("Sorry, you must choose an output map."));
+ }
+
+ err = 0;
+ /* basin and basin.thresh */
+ err += (opt10->answer != NULL && opt6->answer == NULL);
+ /* stream and basin.thresh */
+ err += (opt11->answer != NULL && opt6->answer == NULL);
+ /* half.basin and basin.thresh */
+ err += (opt12->answer != NULL && opt6->answer == NULL);
+ /* slope and basin.thresh */
+ err += (opt15->answer != NULL && opt6->answer == NULL);
+
+ if (err) {
+ G_message(_("Sorry, if any of the following options are set:\n"
+ " basin, stream, half.basin, slope, or lS\n"
+ " you MUST provide a value for the basin "
+ "threshold parameter."));
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+
+ /* Build command line */
+ sprintf(command, "%s/etc/", G_gisbase());
+
+ if (flag2->answer)
+ strcat(command, "r.watershed2.seg");
+ else
+ strcat(command, "r.watershed2.ram");
+
+ if (flag1->answer)
+ strcat(command, " -4");
+
+ if (opt1->answer) {
+ strcat(command, " el=");
+ strcat(command, "\"");
+ strcat(command, opt1->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt2->answer) {
+ strcat(command, " de=");
+ strcat(command, "\"");
+ strcat(command, opt2->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt3->answer) {
+ strcat(command, " ov=");
+ strcat(command, "\"");
+ strcat(command, opt3->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt4->answer) {
+ strcat(command, " r=");
+ strcat(command, "\"");
+ strcat(command, opt4->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt5->answer) {
+ strcat(command, " ob=");
+ strcat(command, "\"");
+ strcat(command, opt5->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt6->answer) {
+ strcat(command, " t=");
+ strcat(command, opt6->answer);
+ }
+
+ if (opt7->answer) {
+ strcat(command, " ms=");
+ strcat(command, opt7->answer);
+ }
+
+ if (opt8->answer) {
+ strcat(command, " ac=");
+ strcat(command, "\"");
+ strcat(command, opt8->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt9->answer) {
+ strcat(command, " dr=");
+ strcat(command, "\"");
+ strcat(command, opt9->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt10->answer) {
+ strcat(command, " ba=");
+ strcat(command, "\"");
+ strcat(command, opt10->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt11->answer) {
+ strcat(command, " se=");
+ strcat(command, "\"");
+ strcat(command, opt11->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt12->answer) {
+ strcat(command, " ha=");
+ strcat(command, "\"");
+ strcat(command, opt12->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt13->answer) {
+ strcat(command, " di=");
+ strcat(command, "\"");
+ strcat(command, opt13->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt14->answer) {
+ strcat(command, " LS=");
+ strcat(command, "\"");
+ strcat(command, opt14->answer);
+ strcat(command, "\"");
+ }
+
+ if (opt15->answer) {
+ strcat(command, " S=");
+ strcat(command, "\"");
+ strcat(command, opt15->answer);
+ strcat(command, "\"");
+ }
+
+ if (flag2->answer && opt16->answer) {
+ strcat(command, " mb=");
+ strcat(command, "\"");
+ strcat(command, opt16->answer);
+ strcat(command, "\"");
+ }
+
+ G_debug(1, "Mode: %s", flag1->answer ? "Segmented" : "All in RAM");
+ G_debug(1, "Running: %s", command);
+
+ ret = system(command);
+
+
+ /* record map metadata/history info */
+ if (opt8->answer)
+ write_hist(opt8->answer,
+ "Watershed accumulation: overland flow that traverses each cell",
+ opt1->answer, flag1->answer);
+ if (opt9->answer)
+ write_hist(opt9->answer,
+ "Watershed drainage direction (divided by 45deg)",
+ opt1->answer, flag1->answer);
+ if (opt10->answer)
+ write_hist(opt10->answer,
+ "Watershed basins", opt1->answer, flag1->answer);
+ if (opt11->answer)
+ write_hist(opt11->answer,
+ "Watershed stream segments", opt1->answer, flag1->answer);
+ if (opt12->answer)
+ write_hist(opt12->answer,
+ "Watershed half-basins", opt1->answer, flag1->answer);
+ if (opt13->answer)
+ write_hist(opt13->answer,
+ "Watershed visualization map (filtered accumulation map)",
+ opt1->answer, flag1->answer);
+ if (opt14->answer)
+ write_hist(opt14->answer,
+ "Watershed slope length and steepness (LS) factor",
+ opt1->answer, flag1->answer);
+ if (opt15->answer)
+ write_hist(opt15->answer,
+ "Watershed slope steepness (S) factor",
+ opt1->answer, flag1->answer);
+
+ exit(ret);
+}
+
+/* record map history info */
+int write_hist(char *map_name, char *title, char *source_name, int mode)
+{
+ struct History history;
+
+ G_put_cell_title(map_name, title);
+
+ G_short_history(map_name, "raster", &history);
+ strncpy(history.datsrc_1, source_name, RECORD_LEN);
+ history.datsrc_1[RECORD_LEN - 1] = '\0'; /* strncpy() doesn't null terminate if maxfill */
+ sprintf(history.edhist[0],
+ "Processing mode: %s", mode ? "Segmented" : "All in RAM");
+ history.edlinecnt = 1;
+ G_command_history(&history);
+
+ return G_write_history(map_name, &history);
+}
Property changes on: grass/trunk/raster/r.watershed2/front/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/front/r.watershed.html
===================================================================
--- grass/trunk/raster/r.watershed2/front/r.watershed.html (rev 0)
+++ grass/trunk/raster/r.watershed2/front/r.watershed.html 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,427 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.watershed</em> generates a set of maps indicating:
+1) the location of watershed basins, and
+2) the LS and S factors of the Revised Universal Soil Loss Equation (RUSLE).
+
+<p>
+<!-- Interactive mode not activated in GRASS 6.
+<em>r.watershed</em> can be run either interactively or non-interactively.
+If the user types <tt>r.watershed</tt>
+on the command line without program arguments, the program will prompt the user
+with a verbose description of the input maps. The interactive version of
+<em>r.watershed</em> can prepare inputs to lumped-parameter hydrologic models.
+After a verbose interactive session, <em>r.watershed</em> will query the user
+for a number of
+map layers. Each map layer's values will be tabulated by watershed basin and sent
+to an output file. This output file is organized to ease data entry into a
+lumped-parameter hydrologic model program. The non-interactive version of
+<em>r.watershed</em> cannot create this file.
+
+<p>
+The user can run the program non-interactively, by specifying input map names
+on the command line. Parameter names may be specified by their
+full names, or by any initial string that distinguish them from other parameter names.
+In <em>r.watershed</em>'s case, the first two letters of each name sufficiently
+distinguishes parameter names.
+For example, the two expressions below are equivalent inputs to <em>r.watershed</em>:
+<p>
+<pre>
+ el=elev.map th=100 st=stream.map ba=basin.map
+
+ elevation=elev.map threshold=100 stream=stream.map basin=basin.map
+</pre>
+-->
+<h2>OPTIONS</h2>
+
+<dl>
+<dt><em>-m</em>
+
+<dd>Without this flag set, the entire analysis is run in memory
+maintained by the operating system. This can be limiting, but is
+relatively fast. Setting the flag causes the program to manage memory
+on disk which allows larger maps to be processes but is considerably
+slower.
+
+<dt><em>-4</em>
+
+<dd>Allow only horizontal and vertical flow of water.
+Stream and slope lengths are approximately the same as outputs from default surface
+flow (allows horizontal, vertical, and diagonal flow of water).
+This flag will also make the drainage basins look more homogeneous.
+
+<dt><em>elevation</em>
+
+<dd>Input map: Elevation on which entire analysis is based.
+
+<dt><em>depression</em>
+
+<dd>Input map: Map layer of actual depressions or sinkholes in the
+landscape that are large enough to slow and store surface runoff from
+a storm event. Any non-zero values indicate depressions.
+
+<dt><em>flow</em>
+
+<dd>Input map: amount of overland flow per cell. This map indicates the
+amount of overland flow units that each cell will contribute to the
+watershed basin model. Overland flow units represent the amount of
+overland flow each cell contributes to surface flow. If omitted, a
+value of one (1) is assumed. The algorithm is D8 flow accumulation.
+
+<dt><em>disturbed_land</em>
+
+<dd>Raster map input layer or value containing the percent of disturbed
+land (i.e., croplands, and construction sites) where the raster or input
+value of 17 equals 17%. If no map or value is given, <em>r.watershed</em>
+assumes no disturbed land. This input is used for the RUSLE calculations.
+
+<dt><em>blocking</em>
+
+<dd>Input map: terrain that will block overland surface flow. Terrain
+that will block overland surface flow and restart the slope length
+for the RUSLE. Any non-zero values indicate blocking terrain.
+
+<dt><em>threshold</em>
+
+<dd>The minimum size of an exterior watershed basin in cells, if no flow
+map is input, or overland flow units when a flow map is given.
+Warning: low threshold values will dramactically increase run time and
+generate difficult too read basin and half_basin results.
+This parameter also controls the level of detail in the <em>stream</em>
+segments map.
+
+<dt><em>max_slope_length</em>
+
+<dd>Input value indicating the maximum length of overland surface flow
+in meters. If overland flow travels greater than the maximum length,
+the program assumes the maximum length (it assumes that landscape
+characteristics not discernible in the digital elevation model exist
+that maximize the slope length). This input is used for the RUSLE calculations
+and is a sensitive parameter.
+
+<dt><em>accumulation</em>
+
+<dd>Output map: The absolute value of each cell in this output map layer is
+the amount of overland flow that traverses the cell. This value will be
+the number of upland cells plus one if no overland flow map is given. If
+the overland flow map is given, the value will be in overland flow units.
+Negative numbers indicate that those cells possibly have surface runoff
+from outside of the current geographic region. Thus, any cells with
+negative values cannot have their surface runoff and sedimentation yields
+calculated accurately.
+
+<dt><em>drainage</em>
+
+<dd>Output map: drainage direction. Provides the "aspect" for each
+cell. Multiplying positive values by 45 will give the direction in
+degrees that the surface runoff will travel from that cell. The
+value -1 indicates that the cell is a depression area (defined by
+the depression input map). Other negative values indicate that
+surface runoff is leaving the boundaries of the current geographic
+region. The absolute value of these negative cells indicates the
+direction of flow.
+
+<dt><em>basin</em>
+
+<dd>Output map: Unique label for each watershed basin. Each basin will
+be given a unique positive even integer. Areas along edges may not
+be large enough to create an exterior watershed basin. 0 values
+indicate that the cell is not part of a complete watershed basin
+in the current geographic region.
+
+<dt><em>stream</em>
+
+<dd>Output map: stream segments. Values correspond to the watershed
+basin values.
+
+<dt><em>half_basin</em>
+
+<dd>Output map: each half-basin is given a unique value. Watershed
+basins are divided into left and right sides. The right-hand side
+cell of the watershed basin (looking upstream) are given even values
+corresponding to the values in basin. The left-hand side
+cells of the watershed basin are given odd values which are one less
+than the value of the watershed basin.
+
+<dt><em>visual</em>
+
+<dd>Output map: useful for visual display of results.
+Surface runoff accumulation with the values
+modified to provide for easy display. All negative accumulation values
+are changed to zero. All positive values above the basin threshold
+are given the value of the <em>threshold</em> parameter.
+
+<dt><em>length_slope</em>
+
+<dd>Output map: slope length and steepness (LS) factor. Contains the LS
+factor for the Revised Universal Soil Loss Equation. Equations taken
+from <em>Revised Universal Soil Loss Equation for Western Rangelands</em>
+(Weltz et al. 1987).
+Since the LS factor is a small number, it is multiplied by 100 for the
+GRASS output map.
+
+<dt><em>slope_steepness</em>
+
+<dd>Output map: slope steepness (S) factor for RUSLE.
+Contains the revised S factor for the Universal Soil
+Loss Equation. Equations taken from article entitled
+<em>Revised Slope Steepness Factor for the Universal Soil
+Loss Equation</em> (McCool et al. 1987). Since the S factor
+is a small number (usually less than one), it is multiplied
+by 100 for the GRASS output map layer.
+</dd>
+</dl>
+
+
+<h2>NOTES</h2>
+
+<em>r.watershed</em> uses an algorithm designed to minimize the impact of
+DEM data errors. This algorithm works slower than <em>r.terraflow</em> but
+provides more accurate results in areas of low slope as well as DEMs
+constructed with techniques that mistake canopy tops as the ground elevation.
+Kinner et al. (2005), for example, used SRTM and IFSAR DEMs to compare
+<em>r.watershed</em> against <em>r.terraflow</em> results in Panama.
+<em>r.terraflow</em> was unable to replicate stream locations in the larger
+valleys while <em>r.watershed</em> performed much better. Thus, if forest
+canopy exists in valleys, SRTM, IFSAR, and similar data products will cause
+major errors in <em>r.terraflow</em> stream output. Under similar conditions,
+<em>r.watershed</em> will generate better <b>stream</b> and <b>half_basin</b>
+results. If watershed divides contain flat to low slope, <em>r.watershed</em>
+will generate better basin results than <em>r.terraflow</em>.
+(<em>r.terraflow</em> uses the same type of algorithm as ESRI's ArcGIS
+watershed software which fails under these conditions.) Also, if watershed
+divides contain forest canopy mixed with uncanopied areas using SRTM, IFSAR,
+and similar data products, <em>r.watershed</em> will generate better basin
+results than <em>r.terraflow</em>.
+
+<p>
+There are two versions of this program: <em>ram</em> and <em>seg</em>.
+Which is version is run depends on whether the <em>-m</em> flag is set.
+<br>
+The <em>ram</em> version uses virtual memory managed by the operating
+system to store all the data structures and is faster than the <em>seg</em>
+version;
+<em>seg</em> uses the GRASS segmentation library which manages data in disk
+files. Thus <em>seg</em> uses much less system memory (RAM) allowing other
+processes to operate on the same CPU, even when the current geographic
+region is huge.
+<br>
+Due to memory requirements of both programs, it is quite easy to run out of
+memory when working with huge map regions. If the <em>ram</em> version runs
+out of memory and the resolution size of the current geographic region
+cannot be increased, either more memory needs to be added to the computer,
+or the swap space size needs to be increased. If <em>seg</em> runs out of
+memory, additional disk space needs to be freed up for the program to run.
+
+<p>
+Both versions use the A<sup>T</sup> least-cost search algorithm to determine
+the flow of water over the landscape (see <a href="#seealso">SEE ALSO</a>
+section).
+The algorithm produces results similar to those obtained when running
+<em><a href="r.cost.html">r.cost</a></em> and
+<em><a href="r.drain.html">r.drain</a></em> on every cell on the map.
+
+<p>
+In many situations, the elevation data will be too finely detailed for
+the amount of time or memory available. Running <em>r.watershed</em> may
+require use of a coarser resolution. To make the results more closely
+resemble the finer terrain data, create a map layer containing the
+lowest elevation values at the coarser resolution. This is done by:
+1) Setting the current geographic region equal to the elevation map
+layer with <em>g.region</em>, and 2) Use the <em>r.neighbors</em> or
+<em>r.resamp.stats</em> command to find the lowest value for an area
+equal in size to the desired resolution. For example, if the resolution
+of the elevation data is 30 meters and the resolution of the geographic
+region for <em>r.watershed</em> will be 90 meters: use the minimum
+function for a 3 by 3 neighborhood. After changing to the resolution at
+which <em>r.watershed</em> will be run, <em>r.watershed</em> should be run
+using the values from the <em>neighborhood</em> output map layer that
+represents the minimum elevation within the region of the coarser cell.
+
+<p>
+The minimum size of drainage basins, defined by the <em>threshold</em>
+parameter, is only relevant for those watersheds with a single stream
+having at least the <em>threshold</em> of cells flowing into it.
+(These watersheds are called exterior basins.)
+Interior drainage basins contain stream segments below multiple tributaries.
+Interior drainage basins can be of any size because the length of
+an interior stream segment is determined by the distance between the
+tributaries flowing into it.
+
+<p>
+The <em>r.watershed</em> program does not require the user to have the
+current geographic region filled with elevation values. Areas without
+elevation data MUST be masked out, by creating a raster map (or raster
+reclassification) named <tt>MASK</tt>. Areas
+masked out will be treated as if they are off the edge of the region.
+MASKs will reduce the memory necessary to run the program. Masking out
+unimportant areas can significantly reduce processing time if the watersheds
+of interest occupy a small percentage of the overall area.
+
+<p>
+Zero data values will be treated as elevation data (not no_data).
+
+<p>
+To isolate an individual river network using the output of this module,
+a number of approaches may be considered.
+<ol>
+<li>Use a resample of the basins catchment raster map as a MASK.<br>
+ The equivalent vector map method is similar using <em>v.select</em> or
+ <em>v.overlay</em>.
+<li>Use the <em>r.cost</em> module with a point in the river as a starting
+ point.
+<li>Use the <em>v.net.iso</em> module with a node in the river as a
+ starting point.
+</ol>
+
+
+<p>
+To create <i>river mile</i> segmentation from a vectorized streams map,
+try the <em>v.net.iso</em> or <em>v.lrs.segment</em> modules.
+
+
+<h2>EXAMPLES</h2>
+<i>These examples use the Spearfish sample dataset.</i>
+<p>
+Convert <em>r.watershed</em> streams map output to a vector layer.
+<p>
+If you want a detailed stream network, set the threshold option
+small to create lots of catchment basins, as only one stream is
+presented per catchment. The r.to.vect -v flag preserves the
+catchment ID as the vector category number.
+
+<div class="code"><pre>
+ r.watershed elev=elevation.dem stream=rwater.stream
+ r.to.vect -v in=rwater.stream out=rwater_stream
+</pre></div>
+<br>
+
+<p>
+Set a nice color table for the accumulation map:
+<div class="code"><pre>
+ MAP=rwater.accum
+ r.watershed elev=elevation.dem accum=$MAP
+
+ eval `r.univar -g "$MAP"`
+ stddev_x_2=`echo $stddev | awk '{print $1 * 2}'`
+ stddev_div_2=`echo $stddev | awk '{print $1 / 2}'`
+
+ r.colors $MAP col=rules << EOF
+ 0% red
+ -$stddev_x_2 red
+ -$stddev yellow
+ -$stddev_div_2 cyan
+ -$mean_of_abs blue
+ 0 white
+ $mean_of_abs blue
+ $stddev_div_2 cyan
+ $stddev yellow
+ $stddev_x_2 red
+ 100% red
+ EOF
+</pre></div>
+<br>
+
+
+<p>
+Create a more detailed stream map using the accumulation map and convert
+it to a vector output map. The accumulation cut-off, and therefore fractal
+dimension, is arbitrary; in this example we use the map's mean number of
+upstream catchment cells (calculated in the above example by <em>r.univar</em>)
+as the cut-off value.
+<div class="code"><pre>
+ r.watershed elev=elevation.dem accum=rwater.accum
+
+ r.mapcalc 'MASK = if(!isnull(elevation.dem))'
+ r.mapcalc "rwater.course = \
+ if( abs(rwater.accum) > $mean_of_abs, \
+ abs(rwater.accum), \
+ null() )"
+ r.colors -g rwater.course col=bcyr
+ g.remove MASK
+
+ # <i>Thinning is required before converting raster lines to vector</i>
+ r.thin in=rwater.course out=rwater.course.Thin
+ r.colors -gn rwater.course.Thin color=grey
+ r.to.vect in=rwater.course.Thin out=rwater_course feature=line
+ v.db.dropcol map=rwater_course column=label
+</pre></div>
+<!-- can't set line attribute to catchment it is in as v.what.rast and
+ v.distance only work for point features. Could create endpoint node
+ points map and upload to that ?? -->
+<!-- Note value column containing accumulation cells in output vector
+ may not necessarily reference the downstream end of the line! drop it? -->
+<br>
+
+<p>
+Create watershed basins map and convert to a vector polygon map
+<div class="code"><pre>
+ r.watershed elev=elevation.dem basin=rwater.basin thresh=15000
+ r.to.vect -s in=rwater.basin out=rwater_basins feature=area
+ v.db.dropcol map=rwater_basins column=label
+ v.db.renamecol map=rwater_basins column=value,catchment
+</pre></div>
+<br>
+
+<p>
+Display output in a nice way
+<div class="code"><pre>
+ r.shaded.relief map=elevation.dem
+ d.shadedmap rel=elevation.dem.shade drape=rwater.basin bright=40
+ d.vect rwater_course color=orange
+</pre></div>
+<br>
+
+<a name="references"></a>
+<h2>REFERENCES</h2>
+
+
+Ehlschlaeger, C. (1989). <i>Using the A<sup>T</sup> Search Algorithm
+to Develop Hydrologic Models from Digital Elevation Data</i>,
+<b>Proceedings of International Geographic Information Systems (IGIS)
+Symposium '89</b>, pp 275-281 (Baltimore, MD, 18-19 March 1989).<br>
+URL: <a href="http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html">
+http://faculty.wiu.edu/CR-Ehlschlaeger2/older/IGIS/paper.html</a>
+
+<p>
+Kinner D., H. Mitasova, R. Harmon, L. Toma, R., Stallard. (2005).
+<i>GIS-based Stream Network Analysis for The Chagres River Basin,
+Republic of Panama</i>. <b>The Rio Chagres: A Multidisciplinary Profile of
+a Tropical Watershed</b>, R. Harmon (Ed.), Springer/Kluwer, p.83-95.<br>
+URL: <a href="http://skagit.meas.ncsu.edu/%7Ehelena/measwork/panama/panama.html">
+http://skagit.meas.ncsu.edu/~helena/measwork/panama/panama.html</a>
+
+<p>
+McCool et al. (1987). <i>Revised Slope Steepness Factor for the Universal
+Soil Loss Equation</i>, <b>Transactions of the ASAE</b> Vol 30(5).
+
+<p>
+Weltz M. A., K. G. Renard, J. R. Simanton (1987). <i>Revised Universal Soil
+Loss Equation for Western Rangelands</i>, <b>U.S.A./Mexico Symposium of
+Strategies for Classification and Management of Native Vegetation for
+Food Production In Arid Zones</b> (Tucson, AZ, 12-16 Oct. 1987).
+
+<a name="seealso"></a>
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a>,
+<a href="r.cost.html">r.cost</a>,
+<a href="r.drain.html">r.drain</a>,
+<a href="r.flow.html">r.flow</a>,
+<!-- <a href="r.flowmd.html">r.flowmd</a>, -->
+<a href="r.neighbors.html">r.neighbors</a>,
+<a href="r.param.scale.html">r.param.scale</a>,
+<a href="r.resamp.interp.html">r.resamp.interp</a>,
+<a href="r.terraflow.html">r.terraflow</a>,
+<a href="r.topidx.html">r.topidx</a>,
+<a href="r.water.outlet.html">r.water.outlet</a>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Original version: Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory<br>
+Markus Metz <markus.metz.giswork at gmail.com>
+<p>
+<i>Last changed: $Date: 2008-09-28 01:36:15 +0200 (Sun, 28 Sep 2008) $</i>
Added: grass/trunk/raster/r.watershed2/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed2/ram/Gwater.h (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/Gwater.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,130 @@
+#ifndef __G_WATER_H__
+#define __G_WATER_H__
+
+
+/* program to map out drainage basin structure */
+/* this one uses the A * search algorithm */
+/* written by Chuck Ehlschlaeger */
+/* last modified 03/26/91 */
+#include <math.h>
+#include <grass/gis.h>
+#include "ramseg.h"
+#include "flag.h"
+
+/* redefining G_malloc allows you to see where in */
+/* program that memory runs out */
+/* #define G_malloc malloc */
+
+#define AR_SIZE 16
+#define AR_INCR 16
+#define SHORT int
+#define NOMASK 1
+#define MIN_SLOPE .00001
+#define MIN_GRADIENT_DEGREES 1
+#define DEG_TO_RAD .017453293 /* pi / 180 */
+#define METER_TO_FOOT 3.281
+#define MAX_BYTES 2000000
+#define PAGE_BLOCK 512
+#define RITE 1
+#define LEFT 2
+#define NEITHER 0
+#define ABS(x) (((x) < 0) ? -(x) : (x))
+#define TSTSTR(a) (fprintf (stderr, "%s\n", a))
+#define TST(a) (fprintf (stderr, "%e\n", (double) (a)))
+
+#define POINT struct points
+POINT
+{
+ SHORT r, c, downr, downc;
+ int nxt;
+};
+
+#ifdef MAIN
+#define GLOBAL
+#define DRAINVAR = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }}
+#define UPDRAINVAR = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }}
+#define NEXTDRVAR = { 1,-1,0,0,-1,1,1,-1 }
+#define NEXTDCVAR = { 0,0,-1,1,1,-1,1,-1 }
+#else
+#define GLOBAL extern
+#define DRAINVAR
+#define UPDRAINVAR
+#define NEXTDRVAR
+#define NEXTDCVAR
+#endif
+
+GLOBAL struct Cell_head window;
+
+GLOBAL int *heap_index, heap_size;
+GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+GLOBAL SHORT nrows, ncols;
+GLOBAL double half_res, diag, max_length, dep_slope;
+GLOBAL int bas_thres, tot_parts;
+GLOBAL FLAG *worked, *in_list, *s_b, *swale;
+GLOBAL RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
+GLOBAL RAMSEG r_h_seg, dep_seg;
+GLOBAL RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
+GLOBAL POINT *astar_pts;
+GLOBAL CELL *dis, *alt, *wat, *asp, *bas, *haf, *r_h, *dep;
+GLOBAL CELL *ril_buf;
+GLOBAL int ril_fd;
+GLOBAL double *s_l, *s_g, *l_s;
+GLOBAL CELL one, zero;
+GLOBAL double ril_value, dzero;
+GLOBAL SHORT sides;
+GLOBAL SHORT drain[3][3] DRAINVAR;
+GLOBAL SHORT updrain[3][3] UPDRAINVAR;
+GLOBAL SHORT nextdr[8] NEXTDRVAR;
+GLOBAL SHORT nextdc[8] NEXTDCVAR;
+GLOBAL char ele_name[GNAME_MAX], *ele_mapset, pit_name[GNAME_MAX], *pit_mapset;
+GLOBAL char run_name[GNAME_MAX], *run_mapset, ob_name[GNAME_MAX], *ob_mapset;
+GLOBAL char ril_name[GNAME_MAX], *ril_mapset, dep_name[GNAME_MAX], *dep_mapset;
+GLOBAL char *this_mapset;
+GLOBAL char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
+GLOBAL char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
+GLOBAL char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+GLOBAL char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+GLOBAL char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+GLOBAL char bas_flag, seg_flag, haf_flag, er_flag;
+GLOBAL char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+GLOBAL FILE *fp;
+/* close_maps.c */
+int close_maps(void);
+/* close_maps2.c */
+int close_array_seg(void);
+/* def_basin.c */
+CELL def_basin(int, int, CELL, double, CELL);
+/* do_astar.c */
+int do_astar(void);
+int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int drop_pt(void);
+double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int replace(SHORT, SHORT, SHORT, SHORT);
+/* do_cum.c */
+int do_cum(void);
+/* find_pour.c */
+int find_pourpts(void);
+/* haf_side.c */
+int haf_basin_side(SHORT, SHORT, SHORT);
+/* init_vars.c */
+int init_vars(int, char *[]);
+int do_legal(char *);
+char *do_exist(char *);
+/* no_stream.c */
+int no_stream(int, int, CELL, double, CELL);
+/* over_cells.c */
+int overland_cells(int, int, CELL, CELL, CELL *);
+/* ramseg.c */
+int size_array(int *, int, int);
+/* sg_factor.c */
+int sg_factor(void);
+int len_slp_equ(double, double, double, int, int);
+/* slope_len.c */
+int slope_length(SHORT, SHORT, SHORT, SHORT);
+/* split_str.c */
+CELL split_stream(int, int, int [], int [], int, CELL, double, CELL);
+/* usage.c */
+void usage(char *);
+
+
+#endif /* __G_WATER_H__ */
Property changes on: grass/trunk/raster/r.watershed2/ram/Gwater.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/Makefile
===================================================================
--- grass/trunk/raster/r.watershed2/ram/Makefile (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.watershed.ram
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: etc
+
+htmletc:
+ @echo "No docs to generate."
Property changes on: grass/trunk/raster/r.watershed2/ram/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/close_maps.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/close_maps.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,152 @@
+#include "Gwater.h"
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int close_maps(void)
+{
+ struct Colors colors;
+ int value, r, c, fd;
+ CELL *buf = NULL;
+
+ if (wat_flag || asp_flag || dis_flag || ls_flag || sl_flag || sg_flag)
+ buf = G_allocate_cell_buf();
+ G_free(alt);
+ if (ls_flag || sg_flag)
+ G_free(r_h);
+
+ if (wat_flag) {
+ fd = G_open_cell_new(wat_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new accum map layer."));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ }
+
+ if (asp_flag) {
+ fd = G_open_cell_new(asp_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new aspect map layer."));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = asp[SEG_INDEX(asp_seg, r, c)];
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ G_init_colors(&colors);
+ G_make_aspect_colors(&colors, 0, 8);
+ G_write_colors(asp_name, this_mapset, &colors);
+ }
+ G_free(asp);
+
+ if (dis_flag) {
+ fd = G_open_cell_new(dis_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new accum map layer."));
+ }
+ else {
+ if (bas_thres <= 0)
+ bas_thres = 60;
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+ if (buf[c] < 0) {
+ buf[c] = 0;
+ }
+ else {
+ value = FLAG_GET(swale, r, c);
+ if (value) {
+ buf[c] = bas_thres;
+ }
+ }
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ G_init_colors(&colors);
+ G_make_rainbow_colors(&colors, 1, 120);
+ G_write_colors(dis_name, this_mapset, &colors);
+ }
+ flag_destroy(swale);
+ /*
+ G_free_colors(&colors);
+ */
+ G_free(wat);
+
+ if (ls_flag) {
+ fd = G_open_cell_new(ls_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new L map layer."));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = l_s[SEG_INDEX(l_s_seg, r, c)] + .5;
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ G_free(l_s);
+ }
+
+ if (sl_flag) {
+ fd = G_open_cell_new(sl_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new slope length map layer."));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = s_l[SEG_INDEX(s_l_seg, r, c)] + .5;
+ if (buf[c] > max_length)
+ buf[c] = max_length + .5;
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ }
+
+ if (sl_flag || ls_flag || sg_flag)
+ G_free(s_l);
+
+ if (sg_flag) {
+ fd = G_open_cell_new(sg_name);
+ if (fd < 0) {
+ G_warning(_("unable to open new S map layer."));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ buf[c] = s_g[SEG_INDEX(s_g_seg, r, c)] * 100 + .5;
+ }
+ G_put_raster_row(fd, buf, CELL_TYPE);
+ }
+ if (G_close_cell(fd) < 0)
+ G_warning(_("Close failed."));
+ }
+ G_free(s_g);
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/close_maps.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/close_maps2.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/close_maps2.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,138 @@
+#include "Gwater.h"
+#include <unistd.h>
+
+int close_array_seg(void)
+{
+ struct Colors colors;
+ int incr, max, red, green, blue, rd, gr, bl, flag;
+ int c, r, map_fd;
+ CELL *cellrow, value;
+ CELL *theseg;
+ RAMSEG thesegseg;
+
+ cellrow = G_allocate_cell_buf();
+ if (seg_flag || bas_flag || haf_flag) {
+ if (seg_flag) {
+ theseg = bas;
+ thesegseg = bas_seg;
+ }
+ else if (bas_flag) {
+ theseg = bas;
+ thesegseg = bas_seg;
+ }
+ else {
+ theseg = haf;
+ thesegseg = haf_seg;
+ }
+ max = -9;
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ if (max < theseg[SEG_INDEX(thesegseg, r, c)])
+ max = theseg[SEG_INDEX(thesegseg, r, c)];
+ }
+ }
+ G_debug(1, "%d basins created", max);
+ G_init_colors(&colors);
+ G_make_random_colors(&colors, 1, max);
+
+ if (max < 10000) {
+ G_set_color((CELL) 0, 0, 0, 0, &colors);
+ r = 1;
+ incr = 0;
+ while (incr >= 0) {
+ G_percent(r, max, 3);
+ for (gr = 130 + incr; gr <= 255; gr += 20) {
+ for (rd = 90 + incr; rd <= 255; rd += 30) {
+ for (bl = 90 + incr; bl <= 255; bl += 40) {
+ flag = 1;
+ while (flag) {
+ G_get_color(r, &red, &green, &blue, &colors);
+ if ((blue * .11 + red * .30 + green * .59) <
+ 100) {
+ G_set_color(r, rd, gr, bl, &colors);
+ flag = 0;
+ }
+ if (++r > max) {
+ gr = rd = bl = 300;
+ flag = 0;
+ incr = -1;
+ }
+ if (r % 200 == 0)
+ G_debug(5,
+ "adjusting colors: r=%d\tof %d basins",
+ r, max);
+ }
+ }
+ }
+ }
+ if (incr >= 0) {
+ incr += 15;
+ if (incr > 120)
+ incr = 7;
+ }
+ }
+ G_percent(r - 1, max, 3); /* finish it */
+ }
+ else
+ G_debug(1,
+ "Too many subbasins to reasonably check neighboring color spread");
+ }
+
+ /* stream segments map */
+ if (seg_flag) {
+ map_fd = G_open_cell_new(seg_name);
+ for (r = 0; r < nrows; r++) {
+ G_set_c_null_value(cellrow, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ value = FLAG_GET(swale, r, c);
+ if (value)
+ cellrow[c] = bas[SEG_INDEX(bas_seg, r, c)];
+ }
+ G_put_raster_row(map_fd, cellrow, CELL_TYPE);
+ }
+ G_close_cell(map_fd);
+ G_write_colors(seg_name, this_mapset, &colors);
+ }
+
+ /* basins map */
+ if (bas_flag) {
+ map_fd = G_open_cell_new(bas_name);
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ cellrow[c] = bas[SEG_INDEX(bas_seg, r, c)];
+ if (cellrow[c] == 0)
+ G_set_c_null_value(cellrow + c, 1);
+ }
+ G_put_raster_row(map_fd, cellrow, CELL_TYPE);
+ }
+ G_close_cell(map_fd);
+ G_write_colors(bas_name, this_mapset, &colors);
+ }
+
+ /* half.basins map */
+ if (haf_flag) {
+ map_fd = G_open_cell_new(haf_name);
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ cellrow[c] = haf[SEG_INDEX(haf_seg, r, c)];
+ if (cellrow[c] == 0)
+ G_set_c_null_value(cellrow + c, 1);
+ }
+ G_put_raster_row(map_fd, cellrow, CELL_TYPE);
+ }
+ G_close_cell(map_fd);
+ G_write_colors(haf_name, this_mapset, &colors);
+ }
+
+ if (seg_flag || bas_flag || haf_flag)
+ G_free_colors(&colors);
+
+ G_free(haf);
+ G_free(bas);
+ G_free(cellrow);
+ if (arm_flag)
+ fclose(fp);
+ close_maps();
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/close_maps2.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/def_basin.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/def_basin.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/def_basin.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,101 @@
+#include "Gwater.h"
+
+CELL def_basin(int row, int col, CELL basin_num,
+ double stream_length, CELL old_elev)
+{
+ int r, rr, c, cc, ct, new_r[9], new_c[9];
+ CELL downdir, direction, asp_value, value, new_elev;
+ SHORT oldupdir, riteflag, leftflag, thisdir;
+
+ for (;;) {
+ bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
+ asp_value = asp[SEG_INDEX(asp_seg, row, col)];
+ FLAG_SET(swale, row, col);
+ if (asp_value < 0)
+ asp_value = -asp_value;
+ ct = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ value = asp[SEG_INDEX(asp_seg, r, c)];
+ if (value < -1)
+ value = -value;
+ if (value == drain[rr][cc]) {
+ value = FLAG_GET(swale, r, c);
+ if (value) {
+ new_r[++ct] = r;
+ new_c[ct] = c;
+ }
+ }
+ }
+ }
+ }
+ if (ct == 0) {
+ no_stream(row, col, basin_num, stream_length, old_elev);
+ return (basin_num);
+ }
+ if (ct >= 2) {
+ basin_num = split_stream(row, col, new_r, new_c, ct,
+ basin_num, stream_length, old_elev);
+ return (basin_num);
+ }
+ oldupdir = drain[row - new_r[1] + 1][col - new_c[1] + 1];
+ downdir = asp[SEG_INDEX(asp_seg, row, col)];
+ if (downdir < 0)
+ downdir = -downdir;
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ direction = asp[SEG_INDEX(asp_seg, r, c)];
+ if (direction == drain[rr][cc]) {
+ thisdir = updrain[rr][cc];
+ switch (haf_basin_side
+ (oldupdir, (SHORT) downdir, thisdir)) {
+ case LEFT:
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_elev);
+ leftflag++;
+ break;
+ case RITE:
+ overland_cells(r, c, basin_num, basin_num,
+ &new_elev);
+ riteflag++;
+ break;
+ }
+ }
+ }
+ }
+ }
+ if (leftflag > riteflag)
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num - 1;
+ else
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
+ if (sides == 8) {
+ if (new_r[1] != row && new_c[1] != col)
+ stream_length += diag;
+ else if (new_r[1] != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += window.ew_res;
+ }
+ else { /* sides == 4 */
+
+ if (asp_value == 2 || asp_value == 6) {
+ if (new_r[1] != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += diag;
+ }
+ else { /* asp_value == 4, 8 */
+
+ if (new_c[1] != col)
+ stream_length += window.ew_res;
+ else
+ stream_length += diag;
+ }
+ }
+ row = new_r[1];
+ col = new_c[1];
+ }
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/def_basin.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/do_astar.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/do_astar.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,298 @@
+#include "Gwater.h"
+#include "do_astar.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+int sift_up(int, CELL);
+
+int do_astar(void)
+{
+ POINT *point;
+ int doer, count;
+ SHORT upr, upc, r, c, ct_dir;
+ CELL alt_val, alt_up, asp_up, wat_val;
+ CELL in_val, drain_val;
+ int index_doer, index_up;
+
+
+ G_message(_("SECTION 2: A * Search."));
+
+ count = 0;
+ first_astar = heap_index[1];
+
+ /* A * Search: search uphill, get downhill path */
+ while (first_astar != -1) {
+ G_percent(count++, do_points, 1);
+
+ /* start with point with lowest elevation, in case of equal elevation
+ * of following points, oldest point = point added earliest */
+ /* old routine: astar_pts[first_astar] (doer = first_astar) */
+ /* new routine: astar_pts[heap_index[1]] */
+
+ doer = heap_index[1];
+
+ point = &(astar_pts[doer]);
+
+ /* drop astar_pts[doer] from heap */
+ /* necessary to protect the current point (doer) from modification */
+ /* equivalent to first_astar = point->next in old code */
+ drop_pt();
+
+ /* can go, dragged on from old code */
+ first_astar = heap_index[1];
+
+ /* downhill path for flow accumulation is set here */
+ /* this path determines the order for flow accumulation calculation */
+ point->nxt = first_cum;
+ first_cum = doer;
+
+ r = point->r;
+ c = point->c;
+
+ G_debug(3, "R:%2d C:%2d, ", r, c);
+
+ index_doer = SEG_INDEX(alt_seg, r, c);
+ alt_val = alt[index_doer];
+ wat_val = wat[index_doer];
+
+ FLAG_SET(worked, r, c);
+
+ /* check all neighbours */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ index_up = SEG_INDEX(alt_seg, upr, upc);
+ /* check that r, c are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ /* check if neighbour is in the list */
+ /* if not, add as new point */
+ in_val = FLAG_GET(in_list, upr, upc);
+ if (in_val == 0) {
+ alt_up = alt[index_up];
+ /* flow direction is set here */
+ add_pt(upr, upc, r, c, alt_up, alt_val);
+ drain_val = drain[upr - r + 1][upc - c + 1];
+ asp[index_up] = drain_val;
+
+ }
+ else {
+ /* check if neighbour has not been worked on,
+ * update values for asp and wat */
+ in_val = FLAG_GET(worked, upr, upc);
+ if (in_val == 0) {
+ asp_up = asp[index_up];
+ if (asp_up < -1) {
+ asp[index_up] = drain[upr - r + 1][upc - c + 1];
+
+ if (wat_val > 0)
+ wat[index_doer] = -wat_val;
+
+ replace(upr, upc, r, c); /* alt_up used to be */
+ }
+ }
+ }
+ }
+ }
+ }
+ G_percent(count, do_points, 1); /* finish it */
+ flag_destroy(worked);
+ flag_destroy(in_list);
+ G_free(heap_index);
+
+ return 0;
+}
+
+/* new add point routine for min heap */
+int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+{
+
+ FLAG_SET(in_list, r, c);
+
+ /* add point to next free position */
+
+ heap_size++;
+
+ if (heap_size > do_points)
+ G_fatal_error(_("heapsize too large"));
+
+ heap_index[heap_size] = nxt_avail_pt;
+
+ astar_pts[nxt_avail_pt].r = r;
+ astar_pts[nxt_avail_pt].c = c;
+ astar_pts[nxt_avail_pt].downr = downr;
+ astar_pts[nxt_avail_pt].downc = downc;
+
+ nxt_avail_pt++;
+
+ /* sift up: move new point towards top of heap */
+
+ sift_up(heap_size, ele);
+
+ return 0;
+}
+
+/* new drop point routine for min heap */
+int drop_pt(void)
+{
+ int child, childr, parent;
+ CELL ele, eler;
+ int i;
+
+ if (heap_size == 1) {
+ heap_index[1] = -1;
+ heap_size = 0;
+ return 0;
+ }
+
+ /* start with root */
+ parent = 1;
+
+ /* sift down: move hole back towards bottom of heap */
+ /* sift-down routine customised for A * Search logic */
+
+ while ((child = GET_CHILD(parent)) <= heap_size) {
+ /* select child with lower ele, if equal, older child
+ * older child is older startpoint for flow path, important
+ * chances are good that the left child will be the older child,
+ * but just to make sure we test */
+ ele =
+ alt[SEG_INDEX
+ (alt_seg, astar_pts[heap_index[child]].r,
+ astar_pts[heap_index[child]].c)];
+ if (child < heap_size) {
+ childr = child + 1;
+ i = 1;
+ while (childr <= heap_size && i < 3) {
+ eler =
+ alt[SEG_INDEX
+ (alt_seg, astar_pts[heap_index[childr]].r,
+ astar_pts[heap_index[childr]].c)];
+ if (eler < ele) {
+ child = childr;
+ ele = eler;
+ /* make sure we get the oldest child */
+ }
+ else if (ele == eler &&
+ heap_index[child] > heap_index[childr]) {
+ child = childr;
+ }
+ childr++;
+ i++;
+ }
+ }
+
+ /* move hole down */
+
+ heap_index[parent] = heap_index[child];
+ parent = child;
+
+ }
+
+ /* hole is in lowest layer, move to heap end */
+ if (parent < heap_size) {
+ heap_index[parent] = heap_index[heap_size];
+
+ ele =
+ alt[SEG_INDEX
+ (alt_seg, astar_pts[heap_index[parent]].r,
+ astar_pts[heap_index[parent]].c)];
+ /* sift up last swapped point, only necessary if hole moved to heap end */
+ sift_up(parent, ele);
+
+ }
+
+ /* the actual drop */
+ heap_size--;
+
+ return 0;
+}
+
+/* standard sift-up routine for d-ary min heap */
+int sift_up(int start, CELL ele)
+{
+ int parent, parentp, child, childp;
+ CELL elep;
+
+ child = start;
+ childp = heap_index[child];
+
+ while (child > 1) {
+ parent = GET_PARENT(child);
+ parentp = heap_index[parent];
+
+ elep =
+ alt[SEG_INDEX
+ (alt_seg, astar_pts[parentp].r, astar_pts[parentp].c)];
+ /* parent ele higher */
+ if (elep > ele) {
+
+ /* push parent point down */
+ heap_index[child] = parentp;
+ child = parent;
+
+ }
+ /* same ele, but parent is younger */
+ else if (elep == ele && parentp > childp) {
+
+ /* push parent point down */
+ heap_index[child] = parentp;
+ child = parent;
+
+ }
+ else
+ /* no more sifting up, found new slot for child */
+ break;
+ }
+
+ /* set heap_index for child */
+ if (child < start) {
+ heap_index[child] = childp;
+ }
+
+ return 0;
+
+}
+
+double
+get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+{
+ double slope;
+
+ if (r == downr)
+ slope = (ele - downe) / window.ew_res;
+ else if (c == downc)
+ slope = (ele - downe) / window.ns_res;
+ else
+ slope = (ele - downe) / diag;
+
+ if (slope < MIN_SLOPE)
+ return (MIN_SLOPE);
+
+ return (slope);
+}
+
+
+/* new replace */
+int replace( /* ele was in there */
+ SHORT upr, SHORT upc, SHORT r, SHORT c)
+/* CELL ele; */
+{
+ int now, heap_run;
+
+ /* find the current neighbour point and
+ * set flow direction to focus point */
+
+ heap_run = 0;
+
+ while (heap_run <= heap_size) {
+ now = heap_index[heap_run];
+ if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
+ astar_pts[now].downr = r;
+ astar_pts[now].downc = c;
+ return 0;
+ }
+ heap_run++;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/do_astar.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed2/ram/do_astar.h (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/do_astar.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,7 @@
+#ifndef __DO_ASTAR_H__
+#define __DO_ASTAR__
+
+#define GET_PARENT(c) (int) (((c) - 2) / 3 + 1)
+#define GET_CHILD(p) (int) (((p) * 3) - 1)
+
+#endif /* __DO_ASTAR_H__ */
Property changes on: grass/trunk/raster/r.watershed2/ram/do_astar.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/do_cum.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/do_cum.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,58 @@
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int do_cum(void)
+{
+ SHORT r, c, dr, dc;
+ CELL is_swale, value, valued;
+ int killer, threshold, count;
+
+ G_message(_("SECTION 3: Accumulating Surface Flow."));
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+ while (first_cum != -1) {
+ G_percent(count++, do_points, 3);
+ killer = first_cum;
+ first_cum = astar_pts[killer].nxt;
+ if ((dr = astar_pts[killer].downr) > -1) {
+ r = astar_pts[killer].r;
+ c = astar_pts[killer].c;
+ dc = astar_pts[killer].downc;
+ value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (ABS(value) >= threshold)
+ FLAG_SET(swale, r, c);
+ valued = wat[SEG_INDEX(wat_seg, dr, dc)];
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
+ is_swale = FLAG_GET(swale, r, c);
+ if (is_swale || ABS(valued) >= threshold) {
+ FLAG_SET(swale, dr, dc);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, dr, dc);
+ }
+ }
+ }
+ G_percent(count, do_points, 3); /* finish it */
+ G_free(astar_pts);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/do_cum.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/find_pour.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/find_pour.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,40 @@
+#include "Gwater.h"
+
+int find_pourpts(void)
+{
+ int value, row, col;
+ double easting, northing, stream_length;
+ CELL old_elev, basin_num;
+
+ basin_num = 0;
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ northing = window.north - (row + .5) * window.ns_res;
+ for (col = 0; col < ncols; col++) {
+ value = FLAG_GET(swale, row, col);
+ if (value && asp[SEG_INDEX(asp_seg, row, col)] < 0) {
+ basin_num += 2;
+ if (arm_flag) {
+ easting = window.west + (col + .5) * window.ew_res;
+ fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
+ (int)basin_num, 0, row, col, easting, northing);
+ }
+ if (col == 0 || col == ncols - 1) {
+ stream_length = .5 * window.ew_res;
+ }
+ else if (row == 0 || row == nrows - 1) {
+ stream_length = .5 * window.ns_res;
+ }
+ else {
+ stream_length = 0.0;
+ }
+ old_elev = alt[SEG_INDEX(alt_seg, row, col)];
+ basin_num =
+ def_basin(row, col, basin_num, stream_length, old_elev);
+ }
+ }
+ }
+ G_percent(nrows, nrows, 1); /* finish it */
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/find_pour.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag.h
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag.h (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,73 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+
+/* flag.[ch] is a set of routines which will set up an array of bits
+** that allow the programmer to "flag" cells in a raster map.
+**
+** FLAG *
+** flag_create(nrows,ncols)
+** int nrows, ncols;
+** opens the structure flag.
+** The flag structure will be a two dimensional array of bits the
+** size of nrows by ncols. Will initalize flags to zero (unset).
+**
+** flag_destroy(flags)
+** FLAG *flags;
+** closes flags and gives the memory back to the system.
+**
+** flag_clear_all(flags)
+** FLAG *flags;
+** sets all values in flags to zero.
+**
+** flag_unset(flags, row, col)
+** FLAG *flags;
+** int row, col;
+** sets the value of (row, col) in flags to zero.
+**
+** flag_set(flags, row, col)
+** FLAG *flags;
+** int row, col;
+** will set the value of (row, col) in flags to one.
+**
+** int
+** flag_get(flags, row, col)
+** FLAG *flags;
+** int row, col;
+** returns the value in flags that is at (row, col).
+**
+** idea by Michael Shapiro
+** code by Chuck Ehlschlaeger
+** April 03, 1989
+*/
+#define FLAG struct _flagsss_
+FLAG
+{
+ int nrows, ncols, leng;
+ unsigned char **array;
+};
+
+#define FLAG_UNSET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] &= ~(1<<((col) & 7))
+
+#define FLAG_SET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+
+#define FLAG_GET(flags,row,col) \
+ (flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
+
+/* flag_clr_all.c */
+int flag_clear_all(FLAG *);
+/* flag_create.c */
+FLAG *flag_create(int, int);
+/* flag_destroy.c */
+int flag_destroy(FLAG *);
+/* flag_get.c */
+int flag_get(FLAG *, int, int);
+/* flag_set.c */
+int flag_set(FLAG *, int, int);
+/* flag_unset.c */
+int flag_unset(FLAG *, int, int);
+
+
+#endif /* __FLAG_H__ */
Property changes on: grass/trunk/raster/r.watershed2/ram/flag.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_clr_all.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_clr_all.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_clr_all.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,14 @@
+#include "flag.h"
+
+int flag_clear_all(FLAG * flags)
+{
+ register int r, c;
+
+ for (r = 0; r < flags->nrows; r++) {
+ for (c = 0; c < flags->leng; c++) {
+ flags->array[r][c] = 0;
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_clr_all.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_create.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_create.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_create.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,27 @@
+#include <grass/gis.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+ unsigned char *temp;
+ FLAG *new_flag;
+ register int i;
+
+ new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+ new_flag->nrows = nrows;
+ new_flag->ncols = ncols;
+ new_flag->leng = (ncols + 7) / 8;
+ new_flag->array =
+ (unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+ temp =
+ (unsigned char *)G_malloc(nrows * new_flag->leng *
+ sizeof(unsigned char));
+
+ for (i = 0; i < nrows; i++) {
+ new_flag->array[i] = temp;
+ temp += new_flag->leng;
+ }
+
+ return (new_flag);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_create.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_destroy.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_destroy.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_destroy.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,11 @@
+#include <grass/gis.h>
+#include "flag.h"
+
+int flag_destroy(FLAG * flags)
+{
+ G_free(flags->array[0]);
+ G_free(flags->array);
+ G_free(flags);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_destroy.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_get.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_get.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_get.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,6 @@
+#include "flag.h"
+
+int flag_get(FLAG * flags, int row, int col)
+{
+ return (flags->array[row][col >> 3] & (1 << (col & 7)));
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_get.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_set.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_set.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_set.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,8 @@
+#include "flag.h"
+
+int flag_set(FLAG * flags, int row, int col)
+{
+ flags->array[row][col >> 3] |= (1 << (col & 7));
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_set.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/flag_unset.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/flag_unset.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/flag_unset.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,8 @@
+#include "flag.h"
+
+int flag_unset(FLAG * flags, int row, int col)
+{
+ flags->array[row][col >> 3] &= ~(1 << (col & 7));
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/flag_unset.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/haf_side.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/haf_side.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/haf_side.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,18 @@
+#include "Gwater.h"
+
+int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+{
+ SHORT newup, newthis;
+
+ newup = updir - downdir;
+ if (newup < 0)
+ newup += 8;
+ newthis = thisdir - downdir;
+ if (newthis < 0)
+ newthis += 8;
+ if (newthis < newup)
+ return (LEFT);
+ if (newthis > newup)
+ return (RITE);
+ return (NEITHER);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/haf_side.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/init_vars.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/init_vars.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,409 @@
+#include <stdlib.h>
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int init_vars(int argc, char *argv[])
+{
+ SHORT r, c;
+ CELL *buf, alt_value, wat_value, asp_value;
+ int fd, index;
+ char MASK_flag;
+
+ G_gisinit(argv[0]);
+ ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
+ ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
+ zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+ nxt_avail_pt = 0;
+ /* dep_flag = 0; */
+ max_length = dzero = 0.0;
+ ril_value = -1.0;
+ /* dep_slope = 0.0; */
+ sides = 8;
+ for (r = 1; r < argc; r++) {
+ if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
+ ele_flag++;
+ else if (sscanf(argv[r], "ac=%[^\n]", wat_name) == 1)
+ wat_flag++;
+ else if (sscanf(argv[r], "dr=%[^\n]", asp_name) == 1)
+ asp_flag++;
+ else if (sscanf(argv[r], "de=%[^\n]", pit_name) == 1)
+ pit_flag++;
+ else if (sscanf(argv[r], "t=%d", &bas_thres) == 1) ;
+ else if (sscanf(argv[r], "ms=%lf", &max_length) == 1) ;
+ else if (sscanf(argv[r], "ba=%[^\n]", bas_name) == 1)
+ bas_flag++;
+ else if (sscanf(argv[r], "se=%[^\n]", seg_name) == 1)
+ seg_flag++;
+ else if (sscanf(argv[r], "ha=%[^\n]", haf_name) == 1)
+ haf_flag++;
+ else if (sscanf(argv[r], "ov=%[^\n]", run_name) == 1)
+ run_flag++;
+ else if (sscanf(argv[r], "ar=%[^\n]", arm_name) == 1)
+ arm_flag++;
+ else if (sscanf(argv[r], "di=%[^\n]", dis_name) == 1)
+ dis_flag++;
+ else if (sscanf(argv[r], "sl=%[^\n]", sl_name) == 1)
+ sl_flag++;
+ else if (sscanf(argv[r], "S=%[^\n]", sg_name) == 1)
+ sg_flag++;
+ else if (sscanf(argv[r], "LS=%[^\n]", ls_name) == 1)
+ ls_flag++;
+ else if (sscanf(argv[r], "ob=%[^\n]", ob_name) == 1)
+ ob_flag++;
+ else if (sscanf(argv[r], "r=%[^\n]", ril_name) == 1) {
+ if (sscanf(ril_name, "%lf", &ril_value) == 0) {
+ ril_value = -1.0;
+ ril_flag++;
+ }
+ }
+ /* else if (sscanf (argv[r], "sd=%[^\n]", dep_name) == 1) dep_flag++; */
+ else if (sscanf(argv[r], "-%d", &sides) == 1) {
+ if (sides != 4)
+ usage(argv[0]);
+ }
+ else
+ usage(argv[0]);
+ }
+ if ((ele_flag != 1)
+ ||
+ ((arm_flag == 1) &&
+ ((bas_thres <= 0) || ((haf_flag != 1) && (bas_flag != 1))))
+ ||
+ ((bas_thres <= 0) &&
+ ((bas_flag == 1) || (seg_flag == 1) || (haf_flag == 1) ||
+ (sl_flag == 1) || (sg_flag == 1) || (ls_flag == 1)))
+ )
+ usage(argv[0]);
+ tot_parts = 4;
+ if (ls_flag || sg_flag)
+ tot_parts++;
+ if (bas_thres > 0)
+ tot_parts++;
+ G_message(_("SECTION 1a (of %1d): Initiating Memory."), tot_parts);
+ this_mapset = G_mapset();
+ if (asp_flag)
+ do_legal(asp_name);
+ if (bas_flag)
+ do_legal(bas_name);
+ if (seg_flag)
+ do_legal(seg_name);
+ if (haf_flag)
+ do_legal(haf_name);
+ if (sl_flag)
+ do_legal(sl_name);
+ if (sg_flag)
+ do_legal(sg_name);
+ if (ls_flag)
+ do_legal(ls_name);
+ if (sl_flag || sg_flag || ls_flag)
+ er_flag = 1;
+ ele_mapset = do_exist(ele_name);
+ if (pit_flag)
+ pit_mapset = do_exist(pit_name);
+ if (ob_flag)
+ ob_mapset = do_exist(ob_name);
+ if (ril_flag)
+ ril_mapset = do_exist(ril_name);
+ /* for sd factor
+ if (dep_flag) {
+ if (sscanf (dep_name, "%lf", &dep_slope) != 1) {
+ dep_mapset = do_exist (dep_name);
+ dep_flag = -1;
+ }
+ }
+ */
+ G_get_set_window(&window);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ total_cells = nrows * ncols;
+ if (max_length <= dzero)
+ max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
+ if (window.ew_res < window.ns_res)
+ half_res = .5 * window.ew_res;
+ else
+ half_res = .5 * window.ns_res;
+ diag = sqrt(window.ew_res * window.ew_res +
+ window.ns_res * window.ns_res);
+ if (sides == 4)
+ diag *= 0.5;
+ buf = G_allocate_cell_buf();
+ alt =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
+ r_h =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+
+ fd = G_open_cell_old(ele_name, ele_mapset);
+ if (fd < 0) {
+ G_fatal_error(_("unable to open elevation map layer"));
+ }
+
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ index = SEG_INDEX(alt_seg, r, c);
+ alt[index] = r_h[index] = buf[c];
+ }
+ }
+ G_close_cell(fd);
+ wat =
+ (CELL *) G_malloc(sizeof(CELL) * size_array(&wat_seg, nrows, ncols));
+
+ if (run_flag) {
+ run_mapset = do_exist(run_name);
+ fd = G_open_cell_old(run_name, run_mapset);
+ if (fd < 0) {
+ G_fatal_error(_("unable to open runoff map layer"));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
+ }
+ }
+ G_close_cell(fd);
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++)
+ wat[SEG_INDEX(wat_seg, r, c)] = 1;
+ }
+ }
+ asp =
+ (CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
+
+ if (pit_flag) {
+ pit_mapset = do_exist(pit_name);
+ fd = G_open_cell_old(pit_name, pit_mapset);
+ if (fd < 0) {
+ G_fatal_error(_("unable to open depression map layer"));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ asp[SEG_INDEX(asp_seg, r, c)] = buf[c];
+ }
+ }
+ G_close_cell(fd);
+ }
+ swale = flag_create(nrows, ncols);
+ if (ob_flag) {
+ fd = G_open_cell_old(ob_name, ob_mapset);
+ if (fd < 0) {
+ G_fatal_error(_("unable to open blocking map layer"));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ if (buf[c])
+ FLAG_SET(swale, r, c);
+ }
+ }
+ G_close_cell(fd);
+ }
+ if (ril_flag) {
+ ril_fd = G_open_cell_old(ril_name, ril_mapset);
+ if (ril_fd < 0) {
+ G_fatal_error(_("unable to open rill map layer"));
+ }
+ }
+ in_list = flag_create(nrows, ncols);
+ worked = flag_create(nrows, ncols);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
+ if (NULL != G_find_file("cell", "MASK", G_mapset())) {
+ MASK_flag = 1;
+ if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
+ G_fatal_error(_("unable to open MASK"));
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row_nomask(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ if (!buf[c]) {
+ FLAG_SET(worked, r, c);
+ FLAG_SET(in_list, r, c);
+ do_points--;
+ }
+ }
+ }
+ G_close_cell(fd);
+ }
+ }
+ s_l =
+ (double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
+ sizeof(double));
+ /* astar_pts = (POINT *) G_malloc (nrows * ncols * sizeof (POINT)); */
+ astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
+
+ if (sg_flag) {
+ s_g =
+ (double *)G_malloc(size_array(&s_g_seg, nrows, ncols) *
+ sizeof(double));
+ }
+ if (ls_flag) {
+ l_s =
+ (double *)G_malloc(size_array(&l_s_seg, nrows, ncols) *
+ sizeof(double));
+ }
+
+ /* heap_index will track astar_pts in the binary min-heap */
+ /* heap_index is one-based */
+ heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
+
+ G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
+
+ /* heap is empty */
+ heap_size = 0;
+
+ first_astar = first_cum = -1;
+ if (MASK_flag) {
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 3);
+ for (c = 0; c < ncols; c++) {
+ if (FLAG_GET(worked, r, c)) {
+ wat[SEG_INDEX(wat_seg, r, c)] = 0;
+ }
+ else {
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ asp_value = asp[SEG_INDEX(asp_seg, r, c)];
+ if (r == 0 || c == 0 || r == nrows - 1 ||
+ c == ncols - 1 || asp_value != 0) {
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ if (r == 0)
+ asp_value = -2;
+ else if (c == 0)
+ asp_value = -4;
+ else if (r == nrows - 1)
+ asp_value = -6;
+ else if (c == ncols - 1)
+ asp_value = -8;
+ else
+ asp_value = -1;
+ asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ }
+ else if (FLAG_GET(worked, r - 1, c)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -2;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (FLAG_GET(worked, r + 1, c)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -6;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (FLAG_GET(worked, r, c - 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -4;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (FLAG_GET(worked, r, c + 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -8;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (sides == 8 && FLAG_GET(worked, r - 1, c - 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -3;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (sides == 8 && FLAG_GET(worked, r - 1, c + 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -1;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (sides == 8 && FLAG_GET(worked, r + 1, c - 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -5;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ else if (sides == 8 && FLAG_GET(worked, r + 1, c + 1)) {
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp[SEG_INDEX(asp_seg, r, c)] = -7;
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0)
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ }
+ }
+ }
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 3);
+ for (c = 0; c < ncols; c++) {
+ s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+ asp_value = asp[SEG_INDEX(asp_seg, r, c)];
+ if (r == 0 || c == 0 || r == nrows - 1 ||
+ c == ncols - 1 || asp_value != 0) {
+ wat_value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (wat_value > 0) {
+ wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
+ }
+ if (r == 0)
+ asp_value = -2;
+ else if (c == 0)
+ asp_value = -4;
+ else if (r == nrows - 1)
+ asp_value = -6;
+ else if (c == ncols - 1)
+ asp_value = -8;
+ else
+ asp_value = -1;
+ asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
+ alt_value = alt[SEG_INDEX(alt_seg, r, c)];
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ }
+ }
+ }
+ }
+
+ G_percent(r, nrows, 3); /* finish it */
+
+ return 0;
+}
+
+int do_legal(char *file_name)
+{
+ if (G_legal_filename(file_name) == -1)
+ G_fatal_error(_("map layer [%s] not legal for GRASS"), file_name);
+
+ return 0;
+}
+
+char *do_exist(char *file_name)
+{
+ char *file_mapset = G_find_cell2(file_name, "");
+
+ if (file_mapset == NULL)
+ G_fatal_error(_("[%s] map not found."), file_name);
+
+ return (file_mapset);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/init_vars.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/main.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/main.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,56 @@
+
+/****************************************************************************
+ *
+ * MODULE: ram - uses virtual memory
+ * AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
+ * Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>,
+ * Brad Douglas <rez touchofmadness.com>, Hamish Bowman <hamish_nospam yahoo.com>
+ * Markus Metz <markus.metz.giswork gmail.com>
+ * PURPOSE: Watershed determination
+ * COPYRIGHT: (C) 1999-2006 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 <stdlib.h>
+#include <unistd.h>
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#undef MAIN
+
+
+int main(int argc, char *argv[])
+{
+ init_vars(argc, argv);
+ do_astar();
+ do_cum();
+ if (sg_flag || ls_flag) {
+ sg_factor();
+ }
+ if (bas_thres <= 0) {
+ G_message(_("SECTION %d: Closing Maps."), tot_parts);
+ close_maps();
+ }
+ else {
+ if (arm_flag) {
+ fp = fopen(arm_name, "w");
+ }
+ bas =
+ (CELL *) G_calloc(sizeof(CELL),
+ size_array(&bas_seg, nrows, ncols));
+ haf =
+ (CELL *) G_calloc(sizeof(CELL),
+ size_array(&haf_seg, nrows, ncols));
+
+ G_message(_("SECTION %d: Watershed determination."), tot_parts - 1);
+ find_pourpts();
+ G_message(_("SECTION %d: Closing Maps."), tot_parts);
+ close_array_seg();
+ }
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/no_stream.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/no_stream.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,101 @@
+#include "Gwater.h"
+
+int no_stream(int row, int col, CELL basin_num,
+ double stream_length, CELL old_elev)
+{
+ int r, rr, c, cc, uprow = 0, upcol = 0;
+ double slope;
+ CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ SHORT updir, riteflag, leftflag, thisdir;
+
+ while (1) {
+ max_drain = -1;
+ for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
+ for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect == drain[rr][cc]) {
+ value = wat[SEG_INDEX(wat_seg, r, c)];
+ if (value < 0)
+ value = -value;
+ if (value > max_drain) {
+ uprow = r;
+ upcol = c;
+ max_drain = value;
+ }
+ }
+ }
+ }
+ }
+ if (max_drain > -1) {
+ updir = drain[row - uprow + 1][col - upcol + 1];
+ downdir = asp[SEG_INDEX(asp_seg, row, col)];
+ if (downdir < 0)
+ downdir = -downdir;
+ if (sides == 8) {
+ if (uprow != row && upcol != col)
+ stream_length += diag;
+ else if (uprow != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += window.ew_res;
+ }
+ else { /* sides == 4 */
+
+ asp_value = asp[SEG_INDEX(asp_seg, uprow, upcol)];
+ if (downdir == 2 || downdir == 6) {
+ if (asp_value == 2 || asp_value == 6)
+ stream_length += window.ns_res;
+ else
+ stream_length += diag;
+ }
+ else { /* downdir == 4,8 */
+
+ if (asp_value == 4 || asp_value == 8)
+ stream_length += window.ew_res;
+ else
+ stream_length += diag;
+ }
+ }
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect == drain[rr][cc]) {
+ thisdir = updrain[rr][cc];
+ if (haf_basin_side(updir,
+ (SHORT) downdir,
+ thisdir) == RITE) {
+ overland_cells(r, c, basin_num, basin_num,
+ &new_ele);
+ riteflag++;
+ }
+ else {
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_ele);
+ leftflag++;
+ }
+ }
+ }
+ }
+ }
+ if (leftflag >= riteflag)
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num - 1;
+ else
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
+ row = uprow;
+ col = upcol;
+ }
+ else {
+ if (arm_flag) {
+ hih_ele = alt[SEG_INDEX(alt_seg, row, col)];
+ slope = (hih_ele - old_elev) / stream_length;
+ if (slope < MIN_SLOPE)
+ slope = MIN_SLOPE;
+ fprintf(fp, " %f %f\n", slope, stream_length);
+ }
+ return 0;
+ }
+ }
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/no_stream.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/over_cells.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/over_cells.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/over_cells.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,39 @@
+#include "Gwater.h"
+#define BIGNEG -9999999
+
+int
+overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
+{
+ int r, rr, c, cc;
+ CELL new_ele, new_max_ele, value;
+
+ bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
+ haf[SEG_INDEX(haf_seg, row, col)] = haf_num;
+ new_max_ele = BIGNEG;
+ for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
+ for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ value = asp[SEG_INDEX(asp_seg, r, c)];
+ if (value == drain[rr][cc]) {
+ if (r != row && c != col) {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ else if (r != row) {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ else {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ }
+ }
+ }
+ }
+ if (new_max_ele == BIGNEG) {
+ *hih_ele = alt[SEG_INDEX(alt_seg, row, col)];
+ }
+ else {
+ *hih_ele = new_max_ele;
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/over_cells.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/ramseg.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/ramseg.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/ramseg.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,15 @@
+#include <stdio.h>
+#include "ramseg.h"
+
+int size_array(int *ram_seg, int nrows, int ncols)
+{
+ int size, segs_in_col;
+
+ segs_in_col = ((nrows - 1) >> RAMSEGBITS) + 1;
+ *ram_seg = ((ncols - 1) >> RAMSEGBITS) + 1;
+ size = ((((nrows - 1) >> RAMSEGBITS) + 1) << RAMSEGBITS) *
+ ((((ncols - 1) >> RAMSEGBITS) + 1) << RAMSEGBITS);
+ size -= ((segs_in_col << RAMSEGBITS) - nrows) << RAMSEGBITS;
+ size -= (*ram_seg << RAMSEGBITS) - ncols;
+ return (size);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/ramseg.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/ramseg.h
===================================================================
--- grass/trunk/raster/r.watershed2/ram/ramseg.h (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/ramseg.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,15 @@
+#ifndef __RAMSEG_H__
+#define __RAMSEG_H__
+
+
+#define RAMSEG int
+#define RAMSEGBITS 4
+#define DOUBLEBITS 8 /* 2 * ramsegbits */
+#define SEGLENLESS 15 /* 2 ^ ramsegbits - 1 */
+
+#define SEG_INDEX(s,r,c) (int) \
+ (((((r) >> RAMSEGBITS) * (s) + ((c) >> RAMSEGBITS)) << DOUBLEBITS) \
+ + (((r) & SEGLENLESS) << RAMSEGBITS) + ((c) & SEGLENLESS))
+
+
+#endif /* __RAMSEG_H__ */
Property changes on: grass/trunk/raster/r.watershed2/ram/ramseg.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/sg_factor.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/sg_factor.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,73 @@
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int sg_factor(void)
+{
+ int r, c;
+ CELL low_elev, hih_elev;
+ double height, length, S, sin_theta;
+
+ G_message(_("SECTION 4: Length Slope determination."));
+
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 3);
+ if (ril_flag) {
+ G_get_c_raster_row(ril_fd, ril_buf, r);
+ }
+ for (c = 0; c < ncols; c++) {
+ low_elev = alt[SEG_INDEX(alt_seg, r, c)];
+ hih_elev = r_h[SEG_INDEX(r_h_seg, r, c)];
+ length = s_l[SEG_INDEX(s_l_seg, r, c)];
+ height = hih_elev - low_elev;
+ if (length > max_length) {
+ height *= max_length / length;
+ length = max_length;
+ }
+ sin_theta = height / sqrt(height * height + length * length);
+ if (height / length < .09)
+ S = 10.8 * sin_theta + .03;
+ else
+ S = 16.8 * sin_theta - .50;
+ if (sg_flag)
+ s_g[SEG_INDEX(s_g_seg, r, c)] = S;
+ if (ls_flag) {
+ length *= METER_TO_FOOT;
+ len_slp_equ(length, sin_theta, S, r, c);
+ }
+ }
+ }
+ G_percent(r, nrows, 3); /* finish it */
+
+ if (ril_flag) {
+ G_free(ril_buf);
+ G_close_cell(ril_fd);
+ }
+
+ return 0;
+}
+
+int len_slp_equ(double slope_length, double sin_theta, double S, int r, int c)
+{
+ double ril, s_l_exp, /* m */
+ rill_ratio, /* Beta */
+ L;
+
+ rill_ratio = (sin_theta / 0.0896) / (3.0 * pow(sin_theta, 0.8) + 0.56);
+ if (ril_flag) {
+ ril = ril_buf[c];
+ }
+ else if (ril_value >= 0.0) {
+ ril = ril_value;
+ }
+ else
+ ril = 0.0;
+ /* rill_ratio equation from Steve Warren */
+ rill_ratio *= .5 + .005 * ril + .0001 * ril * ril;
+ s_l_exp = rill_ratio / (1 + rill_ratio);
+ L = 100 * pow((slope_length / 72.6), s_l_exp);
+ l_s[SEG_INDEX(l_s_seg, r, c)] = L * S;
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/sg_factor.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/slope_len.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/slope_len.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,51 @@
+#include "Gwater.h"
+
+int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+{
+ CELL top_alt, bot_alt, asp_value;
+ double res, top_ls, bot_ls;
+
+ if (sides == 8) {
+ if (r == dr)
+ res = window.ns_res;
+ else if (c == dc)
+ res = window.ew_res;
+ else
+ res = diag;
+ }
+ else { /* sides == 4 */
+
+ asp_value = asp[SEG_INDEX(asp_seg, dr, dc)];
+ if (r == dr) {
+ if (asp_value == 2 || asp_value == 6)
+ res = window.ns_res;
+ else /* asp_value == 4, 8, -2, -4, -6, or -8 */
+ res = diag;
+ }
+ else { /* c == dc */
+
+ if (asp_value == 4 || asp_value == 8)
+ res = window.ew_res;
+ else /* asp_value == 2, 6, -2, -4, -6, or -8 */
+ res = diag;
+ }
+ }
+ top_ls = s_l[SEG_INDEX(s_l_seg, r, c)];
+ if (top_ls == half_res)
+ top_ls = res;
+ else
+ top_ls += res;
+ s_l[SEG_INDEX(s_l_seg, r, c)] = top_ls;
+ top_alt = alt[SEG_INDEX(alt_seg, r, c)];
+ bot_alt = alt[SEG_INDEX(alt_seg, dr, dc)];
+ if (top_alt > bot_alt) {
+ bot_ls = s_l[SEG_INDEX(s_l_seg, dr, dc)];
+ if (top_ls > bot_ls) {
+ bot_ls = top_ls + res;
+ s_l[SEG_INDEX(s_l_seg, dr, dc)] = bot_ls;
+ r_h[SEG_INDEX(r_h_seg, dr, dc)] = r_h[SEG_INDEX(s_l_seg, r, c)];
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/slope_len.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/split_str.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/split_str.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/split_str.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,87 @@
+#include "Gwater.h"
+
+CELL
+split_stream(int row, int col, int new_r[], int new_c[], int ct,
+ CELL basin_num, double stream_length, CELL old_elev)
+{
+ CELL downdir, old_basin, new_elev, aspect;
+ double slope, easting, northing;
+ SHORT doit, ctr, updir, splitdir[9];
+ SHORT thisdir, leftflag, riteflag;
+ int r, c, rr, cc;
+
+ for (ctr = 1; ctr <= ct; ctr++)
+ splitdir[ctr] = drain[row - new_r[ctr] + 1][col - new_c[ctr] + 1];
+ updir = splitdir[1];
+ downdir = asp[SEG_INDEX(asp_seg, row, col)];
+ if (downdir < 0)
+ downdir = -downdir;
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ aspect = asp[SEG_INDEX(asp_seg, r, c)];
+ if (aspect == drain[rr][cc]) {
+ doit = 1;
+ thisdir = updrain[rr][cc];
+ for (ctr = 1; ctr <= ct; ctr++) {
+ if (thisdir == splitdir[ctr]) {
+ doit = 0;
+ ctr = ct;
+ }
+ }
+ if (doit) {
+ thisdir = updrain[rr][cc];
+ switch (haf_basin_side
+ (updir, (SHORT) downdir, thisdir)) {
+ case LEFT:
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_elev);
+ leftflag++;
+ break;
+ case RITE:
+ overland_cells(r, c, basin_num, basin_num,
+ &new_elev);
+ riteflag++;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (leftflag >= riteflag) {
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num - 1;
+ }
+ else {
+ haf[SEG_INDEX(haf_seg, row, col)] = basin_num;
+ }
+ old_basin = basin_num;
+ new_elev = alt[SEG_INDEX(alt_seg, row, col)];
+ if ((slope = (new_elev - old_elev) / stream_length) < MIN_SLOPE)
+ slope = MIN_SLOPE;
+ if (arm_flag)
+ fprintf(fp, " %f %f\n", slope, stream_length);
+ for (r = 1; r <= ct; r++) {
+ basin_num += 2;
+ easting = window.west + (new_c[r] + .5) * window.ew_res;
+ northing = window.north - (new_r[r] + .5) * window.ns_res;
+ if (arm_flag) {
+ fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
+ (int)basin_num, old_basin, new_r[r], new_c[r], easting,
+ northing);
+ }
+ if (new_r[r] != row && new_c[r] != col)
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, diag, new_elev);
+ else if (new_r[r] != row)
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, window.ns_res,
+ new_elev);
+ else
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, window.ew_res,
+ new_elev);
+ }
+ return (basin_num);
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/split_str.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/ram/usage.c
===================================================================
--- grass/trunk/raster/r.watershed2/ram/usage.c (rev 0)
+++ grass/trunk/raster/r.watershed2/ram/usage.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,35 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+void usage(char *me)
+{
+ G_fatal_error(_("USAGE for basin delineation:\n%s -4 el=elevation_map "
+ "t=swale_threshold [ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "ba=watershed_basin_map [se=stream_segment_map]\n\nUSAGE for "
+ "ARMSED FILE creation:\n%s [-4] el=elevation_map t=swale_threshold "
+ "[ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "[ba=watershed_basin_map] [se=stream_segment_map] "
+ "ha=half_basin_map ar=ARMSED_file_name\n\nUSAGE for slope length "
+ "determination:\n%s [-4] el=elevation_map t=swale_threshold "
+ "[dr=drain_direction_map] [de=depression_map] "
+ "[ac=accumulation_map] [di=display_map] [ms=max_slope_length] "
+ "[ob=overland_blocking_map] [S=slope_steepness_map] "
+ "LS=length_slope_map [r=rill_erosion_map] "
+ "[sd=slope_deposition value or map]"), me, me, me);
+ /*
+ G_fatal_error(_("USAGE for basin delineation:\n%s -4 el=elevation_map "
+ "t=swale_threshold [ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "ba=watershed_basin_map [se=stream_segment_map]\n\nUSAGE for "
+ "ARMSED FILE creation:\n%s [-4] el=elevation_map t=swale_threshold "
+ "[ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "[ba=watershed_basin_map] [se=stream_segment_map] "
+ "ha=half_basin_map ar=ARMSED_file_name"), me, me);
+ */
+}
Property changes on: grass/trunk/raster/r.watershed2/ram/usage.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed2/seg/Gwater.h (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/Gwater.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,130 @@
+#ifndef __G_WATER_H__
+#define __G_WATER_H__
+
+
+/* program to map out drainage basin structure */
+/* this one uses the A * search algorithm */
+/* written by Chuck Ehlschlaeger */
+/* last modified 03/26/91 */
+#include <math.h>
+#include <grass/gis.h>
+#include "cseg.h"
+
+#define AR_SIZE 16
+#define AR_INCR 16
+#define SHORT int
+#define NOMASK 1
+#define MIN_SLOPE .00001
+#define MIN_GRADIENT_DEGREES 1
+#define DEG_TO_RAD .017453293 /* pi / 180 */
+#define METER_TO_FOOT 3.281
+#define MAX_BYTES 10485760
+#define PAGE_BLOCK 1024
+#define SROW 200
+#define SCOL 200
+#define RITE 1
+#define LEFT 2
+#define NEITHER 0
+#define ABS(x) (((x) < 0) ? -(x) : (x))
+#define TSTSTR(a) (fprintf (stderr, "%s\n", a))
+#define TST(a) (fprintf (stderr, "%e\n", (double) (a)))
+
+#define POINT struct points
+POINT
+{
+ SHORT r, c, downr, downc;
+ int nxt;
+};
+
+#ifdef MAIN
+#define GLOBAL
+#define DRAINVAR = {{ 7,6,5 },{ 8,0,4 },{ 1,2,3 }}
+#define UPDRAINVAR = {{ 3,2,1 },{ 4,0,8 },{ 5,6,7 }}
+#define NEXTDRVAR = { 1,-1,0,0,-1,1,1,-1 }
+#define NEXTDCVAR = { 0,0,-1,1,1,-1,1,-1 }
+#else
+#define GLOBAL extern
+#define DRAINVAR
+#define UPDRAINVAR
+#define NEXTDRVAR
+#define NEXTDCVAR
+#endif
+
+#define HEAP struct heap_item
+HEAP
+{
+ int point;
+ CELL ele;
+};
+
+GLOBAL struct Cell_head window;
+
+GLOBAL SSEG heap_index;
+GLOBAL int heap_size;
+GLOBAL int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
+GLOBAL SHORT nrows, ncols;
+GLOBAL double half_res, diag, max_length, dep_slope;
+GLOBAL int bas_thres, tot_parts;
+GLOBAL SSEG astar_pts;
+GLOBAL BSEG worked, in_list, s_b, swale;
+GLOBAL CSEG dis, alt, wat, asp, bas, haf, r_h, dep;
+GLOBAL DSEG slp, s_l, s_g, l_s, ril;
+GLOBAL CELL one, zero;
+GLOBAL double ril_value, dzero;
+GLOBAL SHORT sides;
+GLOBAL SHORT drain[3][3] DRAINVAR;
+GLOBAL SHORT updrain[3][3] UPDRAINVAR;
+GLOBAL SHORT nextdr[8] NEXTDRVAR;
+GLOBAL SHORT nextdc[8] NEXTDCVAR;
+GLOBAL char ele_name[GNAME_MAX], *ele_mapset, pit_name[GNAME_MAX], *pit_mapset;
+GLOBAL char run_name[GNAME_MAX], *run_mapset, ob_name[GNAME_MAX], *ob_mapset;
+GLOBAL char ril_name[GNAME_MAX], *ril_mapset, dep_name[GNAME_MAX], *dep_mapset;
+
+GLOBAL char *this_mapset;
+GLOBAL char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
+GLOBAL char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
+GLOBAL char wat_name[GNAME_MAX], asp_name[GNAME_MAX], arm_name[GNAME_MAX], dis_name[GNAME_MAX];
+GLOBAL char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+GLOBAL char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+GLOBAL char bas_flag, seg_flag, haf_flag, er_flag;
+GLOBAL char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
+GLOBAL FILE *fp;
+
+/* close_maps.c */
+int close_maps(void);
+/* close_maps2.c */
+int close_array_seg(void);
+/* def_basin.c */
+CELL def_basin(int, int, CELL, double, CELL);
+/* do_astar.c */
+int do_astar(void);
+int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int drop_pt(void);
+double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int replace(SHORT, SHORT, SHORT, SHORT);
+/* do_cum.c */
+int do_cum(void);
+/* find_pour.c */
+int find_pourpts(void);
+/* haf_side.c */
+int haf_basin_side(SHORT, SHORT, SHORT);
+/* init_vars.c */
+int init_vars(int, char *[]);
+int do_legal(char *);
+char *do_exist(char *);
+/* no_stream.c */
+int no_stream(int, int, CELL, double, CELL);
+/* over_cells.c */
+int overland_cells(int, int, CELL, CELL, CELL *);
+/* sg_factor.c */
+int sg_factor(void);
+int len_slp_equ(double, double, double, int, int);
+/* slope_len.c */
+int slope_length(SHORT, SHORT, SHORT, SHORT);
+/* split_str.c */
+CELL split_stream(int, int, int [], int [], int, CELL, double, CELL);
+/* usage.c */
+void usage(char *);
+
+
+#endif /* __G_WATER_H__ */
Property changes on: grass/trunk/raster/r.watershed2/seg/Gwater.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/Makefile
===================================================================
--- grass/trunk/raster/r.watershed2/seg/Makefile (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.watershed.seg
+
+LIBES = $(GISLIB) $(SEGMENTLIB)
+DEPENDENCIES = $(GISDEP) $(SEGMENTDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: etc
+
+htmletc:
+ @echo "No docs to generate."
Property changes on: grass/trunk/raster/r.watershed2/seg/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_close.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_close.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_close.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,19 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include "cseg.h"
+
+int bseg_close(BSEG * bseg)
+{
+ segment_release(&(bseg->seg));
+ close(bseg->fd);
+ unlink(bseg->filename);
+ if (bseg->name) {
+ G_free(bseg->name);
+ bseg->name = NULL;
+ }
+ if (bseg->mapset) {
+ G_free(bseg->mapset);
+ bseg->mapset = NULL;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_close.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_get.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_get.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,15 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int bseg_get(BSEG * bseg, CELL * value, int row, int col)
+{
+ CELL x;
+
+ if (segment_get(&(bseg->seg), &x, row, col >> 3) < 0) {
+ G_warning("bseg_get(): could not read segment file");
+ return -1;
+ }
+ *value = (CELL) ((x & (1 << (col & 7))) ? 1 : 0);
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_get.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_open.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_open.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_open.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,59 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+
+int bseg_open(BSEG * bseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ bseg->filename = NULL;
+ bseg->fd = -1;
+ bseg->name = NULL;
+ bseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning("bseg_open(): unable to create segment file");
+ return -2;
+ }
+ if (0 > (errflag = segment_format(fd, G_window_rows(),
+ (G_window_cols() + 7) / 8, srows, scols,
+ sizeof(char)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("bseg_open(): could not write segment file");
+ return -1;
+ }
+ else {
+ G_warning("bseg_open(): illegal configuration parameter(s)");
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning("bseg_open(): unable to re-open segment file");
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(bseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("bseg_open(): could not read segment file");
+ return -5;
+ }
+ else {
+ G_warning("bseg_open(): out of memory");
+ return -6;
+ }
+ }
+ bseg->filename = filename;
+ bseg->fd = fd;
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_put.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_put.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,22 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int bseg_put(BSEG * bseg, CELL * value, int row, int col)
+{
+ CELL old_value;
+
+ if (segment_get(&(bseg->seg), &old_value, row, col >> 3) < 0) {
+ G_warning("bseg_put(): could not read segment file");
+ return -1;
+ }
+ if (*value)
+ old_value |= (1 << (col & 7));
+ else
+ old_value &= ~(1 << (col & 7));
+ if (segment_put(&(bseg->seg), &old_value, row, col >> 3) < 0) {
+ G_warning("bseg_put(): could not write segment file");
+ return -2;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_put.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_read.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_read.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,49 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include "cseg.h"
+
+static char *me = "bseg_read_cell";
+
+int bseg_read_cell(BSEG * bseg, char *map_name, char *mapset)
+{
+ int row, nrows;
+ int col, ncols;
+ int map_fd;
+ char msg[100];
+ CELL *buffer;
+
+ bseg->name = NULL;
+ bseg->mapset = NULL;
+
+ map_fd = G_open_cell_old(map_name, mapset);
+ if (map_fd < 0) {
+ sprintf(msg, "%s(): unable to open file [%s] in [%s], %d",
+ me, map_name, mapset, map_fd);
+ G_warning(msg);
+ return -3;
+ }
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ buffer = G_allocate_cell_buf();
+ for (row = 0; row < nrows; row++) {
+ if (G_get_c_raster_row(map_fd, buffer, row) < 0) {
+ G_free(buffer);
+ G_close_cell(map_fd);
+ sprintf(msg, "%s(): unable to read file [%s] in [%s], %d %d",
+ me, map_name, mapset, row, nrows);
+ G_warning(msg);
+ return -2;
+ }
+ for (col = ncols; col >= 0; col--) {
+ bseg_put(bseg, &(buffer[col]), row, col);
+ }
+ }
+
+ G_close_cell(map_fd);
+ G_free(buffer);
+
+ bseg->name = G_store(map_name);
+ bseg->mapset = G_store(mapset);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_read.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/bseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/bseg_write.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/bseg_write.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,38 @@
+#include <grass/gis.h>
+#include "cseg.h"
+
+static char *me = "bseg_write_cell";
+
+int bseg_write_cellfile(BSEG * bseg, char *map_name)
+{
+ int map_fd;
+ int row, nrows;
+ int col, ncols;
+ CELL *buffer;
+ CELL value;
+
+ map_fd = G_open_cell_new(map_name);
+ if (map_fd < 0) {
+ G_warning("%s(): unable to open new map layer [%s]", me, map_name);
+ return -1;
+ }
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ buffer = G_allocate_cell_buf();
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ bseg_get(bseg, &value, row, col);
+ buffer[col] = value;
+ }
+ if (G_put_raster_row(map_fd, buffer, CELL_TYPE) < 0) {
+ G_free(buffer);
+ G_unopen_cell(map_fd);
+ G_warning("%s(): unable to write new map layer [%s], row %d",
+ me, map_name, row);
+ return -2;
+ }
+ }
+ G_free(buffer);
+ G_close_cell(map_fd);
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/bseg_write.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/close_maps.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/close_maps.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,77 @@
+#include "Gwater.h"
+#include <unistd.h>
+
+int close_maps(void)
+{
+ struct Colors colors;
+ int r, c;
+ CELL is_swale, value;
+ double dvalue;
+
+ dseg_close(&slp);
+ cseg_close(&alt);
+ if (wat_flag)
+ cseg_write_cellfile(&wat, wat_name);
+ if (asp_flag) {
+ cseg_write_cellfile(&asp, asp_name);
+ G_init_colors(&colors);
+ G_make_grey_scale_colors(&colors, 1, 8);
+ G_write_colors(asp_name, this_mapset, &colors);
+ }
+ cseg_close(&asp);
+ if (dis_flag) {
+ if (bas_thres <= 0)
+ bas_thres = 60;
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ cseg_get(&wat, &value, r, c);
+ if (value < 0) {
+ value = 0;
+ cseg_put(&wat, &value, r, c);
+ }
+ else {
+ bseg_get(&swale, &is_swale, r, c);
+ if (is_swale) {
+ value = bas_thres;
+ cseg_put(&wat, &value, r, c);
+ }
+ }
+ }
+ }
+ cseg_write_cellfile(&wat, dis_name);
+ G_init_colors(&colors);
+ G_make_rainbow_colors(&colors, 1, 120);
+ G_write_colors(dis_name, this_mapset, &colors);
+ }
+ /* error in gislib.a
+ G_free_colors(&colors);
+ */
+ cseg_close(&wat);
+ if (ls_flag) {
+ dseg_write_cellfile(&l_s, ls_name);
+ dseg_close(&l_s);
+ }
+ bseg_close(&swale);
+ if (sl_flag) {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ dseg_get(&s_l, &dvalue, r, c);
+ if (dvalue > max_length)
+ dseg_put(&s_l, &max_length, r, c);
+ }
+ }
+ dseg_write_cellfile(&s_l, sl_name);
+ }
+ if (sl_flag || ls_flag || sg_flag)
+ dseg_close(&s_l);
+ if (ril_flag)
+ dseg_close(&ril);
+ if (sg_flag)
+ dseg_write_cellfile(&s_g, sg_name);
+ if (sg_flag)
+ dseg_close(&s_g);
+ if (ls_flag || sg_flag)
+ cseg_close(&r_h);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/close_maps.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/close_maps2.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/close_maps2.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,109 @@
+#include "Gwater.h"
+#include <unistd.h>
+
+int close_array_seg(void)
+{
+ struct Colors colors;
+ int incr, max, red, green, blue, rd, gr, bl, flag;
+ int c, r, map_fd;
+ CELL *cellrow, value;
+ CSEG *theseg;
+
+ if (seg_flag || bas_flag || haf_flag) {
+ if (seg_flag)
+ theseg = &bas;
+ else if (bas_flag)
+ theseg = &bas;
+ else
+ theseg = &haf;
+ max = -9;
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ cseg_get(theseg, &value, r, c);
+ if (value > max)
+ max = value;
+ }
+ }
+ G_debug(1, "%d basins created", max);
+ G_init_colors(&colors);
+ G_make_random_colors(&colors, 1, max);
+
+ if (max < 10000) {
+ G_set_color((CELL) 0, 0, 0, 0, &colors);
+ r = 1;
+ incr = 0;
+ while (incr >= 0) {
+ G_percent(r, max, 3);
+ for (gr = 130 + incr; gr <= 255; gr += 20) {
+ for (rd = 90 + incr; rd <= 255; rd += 30) {
+ for (bl = 90 + incr; bl <= 255; bl += 40) {
+ flag = 1;
+ while (flag) {
+ G_get_color(r, &red, &green, &blue, &colors);
+ if ((blue * .11 + red * .30 + green * .59) <
+ 100) {
+ G_set_color(r, rd, gr, bl, &colors);
+ flag = 0;
+ }
+ if (++r > max) {
+ gr = rd = bl = 300;
+ flag = 0;
+ incr = -1;
+ }
+ }
+ }
+ }
+ }
+ if (incr >= 0) {
+ incr += 15;
+ if (incr > 120)
+ incr = 7;
+ }
+ }
+ G_percent(r - 1, max, 3); /* finish it */
+ }
+ else
+ G_debug(1,
+ "Too many subbasins to reasonably check neighboring color spread");
+ }
+
+ /* stream segments map */
+ if (seg_flag) {
+ cellrow = (CELL *) G_malloc(ncols * sizeof(CELL));
+ map_fd = G_open_cell_new(seg_name);
+ for (r = 0; r < nrows; r++) {
+ G_set_c_null_value(cellrow, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ bseg_get(&swale, &value, r, c);
+ if (value)
+ cseg_get(&bas, &(cellrow[c]), r, c);
+ }
+ G_put_raster_row(map_fd, cellrow, CELL_TYPE);
+ }
+ G_free(cellrow);
+ G_close_cell(map_fd);
+ G_write_colors(seg_name, this_mapset, &colors);
+ }
+
+ /* basins map */
+ if (bas_flag) {
+ cseg_write_cellfile(&bas, bas_name);
+ G_write_colors(bas_name, this_mapset, &colors);
+ }
+
+ /* half.basins map */
+ if (haf_flag) {
+ cseg_write_cellfile(&haf, haf_name);
+ G_write_colors(haf_name, this_mapset, &colors);
+ }
+
+ if (seg_flag || bas_flag || haf_flag)
+ G_free_colors(&colors);
+ cseg_close(&haf);
+ cseg_close(&bas);
+ if (arm_flag)
+ fclose(fp);
+ close_maps();
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/close_maps2.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg.h
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg.h (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,98 @@
+#ifndef __CSEG_H__
+#define __CSEG_H__
+
+
+#include <grass/segment.h>
+
+#define CSEG struct _c_s_e_g_
+CSEG
+{
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+} ;
+
+#define DSEG struct _d_s_e_g_
+DSEG
+{
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+} ;
+
+#define BSEG struct _b_s_e_g_
+BSEG
+{
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+ char *name; /* raster map read into segment file */
+ char *mapset;
+} ;
+
+#define SSEG struct _s_s_e_g_
+SSEG
+{
+ SEGMENT seg; /* segment structure */
+ int fd; /* fd for reading/writing segment file */
+ char *filename; /* name of segment file */
+} ;
+/* bseg_close.c */
+int bseg_close(BSEG *);
+/* bseg_get.c */
+int bseg_get(BSEG *, CELL *, int, int);
+/* bseg_open.c */
+int bseg_open(BSEG *, int, int, int);
+/* bseg_put.c */
+int bseg_put(BSEG *, CELL *, int, int);
+/* bseg_read.c */
+int bseg_read_cell(BSEG *, char *, char *);
+/* bseg_write.c */
+int bseg_write_cellfile(BSEG *, char *);
+/* cseg_close.c */
+int cseg_close(CSEG *);
+/* cseg_get.c */
+int cseg_get(CSEG *, CELL *, int, int);
+/* cseg_open.c */
+int cseg_open(CSEG *, int, int, int);
+/* cseg_put.c */
+int cseg_put(CSEG *, CELL *, int, int);
+/* cseg_read.c */
+int cseg_read_cell(CSEG *, char *, char *);
+/* cseg_write.c */
+int cseg_write_cellfile(CSEG *, char *);
+/* dseg_close.c */
+int dseg_close(DSEG *);
+/* dseg_get.c */
+int dseg_get(DSEG *, double *, int, int);
+/* dseg_open.c */
+int dseg_open(DSEG *, int, int, int);
+/* dseg_put.c */
+int dseg_put(DSEG *, double *, int, int);
+/* dseg_read.c */
+int dseg_read_cell(DSEG *, char *, char *);
+/* dseg_write.c */
+int dseg_write_cellfile(DSEG *, char *);
+/* seg_close.c */
+int seg_close(SSEG *);
+/* seg_get.c */
+int seg_get(SSEG *, char *, int, int);
+/* seg_open.c */
+int seg_open(SSEG *, int, int, int, int, int, int);
+/* seg_put.c */
+int seg_put(SSEG *, char *, int, int);
+/* sseg_close.c */
+int seg_close(SSEG *);
+/* sseg_get.c */
+int seg_get(SSEG *, char *, int, int);
+/* sseg_open.c */
+int seg_open(SSEG *, int, int, int, int, int, int);
+/* sseg_put.c */
+int seg_put(SSEG *, char *, int, int);
+
+
+#endif /* __CSEG_H__ */
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_close.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_close.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_close.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,20 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int cseg_close(CSEG * cseg)
+{
+ segment_release(&(cseg->seg));
+ close(cseg->fd);
+ unlink(cseg->filename);
+ if (cseg->name) {
+ G_free(cseg->name);
+ cseg->name = NULL;
+ }
+ if (cseg->mapset) {
+ G_free(cseg->mapset);
+ cseg->mapset = NULL;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_close.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_get.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_get.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int cseg_get(CSEG * cseg, CELL * value, int row, int col)
+{
+ if (segment_get(&(cseg->seg), value, row, col) < 0) {
+ G_warning("cseg_get(): could not read segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_get.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_open.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_open.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_open.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,59 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include "cseg.h"
+
+
+int cseg_open(CSEG * cseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ cseg->filename = NULL;
+ cseg->fd = -1;
+ cseg->name = NULL;
+ cseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning("cseg_open(): unable to create segment file");
+ return -2;
+ }
+ if (0 >
+ (errflag =
+ segment_format(fd, G_window_rows(), G_window_cols(), srows, scols,
+ sizeof(CELL)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("cseg_open(): could not write segment file");
+ return -1;
+ }
+ else {
+ G_warning("cseg_open(): illegal configuration parameter(s)");
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning("cseg_open(): unable to re-open segment file");
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(cseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("cseg_open(): could not read segment file");
+ return -5;
+ }
+ else {
+ G_warning("cseg_open(): out of memory");
+ return -6;
+ }
+ }
+ cseg->filename = filename;
+ cseg->fd = fd;
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_put.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_put.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int cseg_put(CSEG * cseg, CELL * value, int row, int col)
+{
+ if (segment_put(&(cseg->seg), value, row, col) < 0) {
+ G_warning("cseg_put(): could not write segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_put.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_read.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_read.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,52 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+static char *me = "cseg_read_cell";
+
+int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
+{
+ int row, nrows;
+ int map_fd;
+ char msg[100];
+ CELL *buffer;
+
+ cseg->name = NULL;
+ cseg->mapset = NULL;
+
+ map_fd = G_open_cell_old(map_name, mapset);
+ if (map_fd < 0) {
+ sprintf(msg, "%s(): unable to open file [%s] in [%s], %d",
+ me, map_name, mapset, map_fd);
+ G_warning(msg);
+ return -3;
+ }
+ nrows = G_window_rows();
+ buffer = G_allocate_cell_buf();
+ for (row = 0; row < nrows; row++) {
+ if (G_get_c_raster_row(map_fd, buffer, row) < 0) {
+ G_free(buffer);
+ G_close_cell(map_fd);
+ sprintf(msg, "%s(): unable to read file [%s] in [%s], %d %d",
+ me, map_name, mapset, row, nrows);
+ G_warning(msg);
+ return -2;
+ }
+ if (segment_put_row(&(cseg->seg), buffer, row) < 0) {
+ G_free(buffer);
+ G_close_cell(map_fd);
+ sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+ me, map_name, mapset);
+ G_warning(msg);
+ return (-1);
+ }
+ }
+
+ G_close_cell(map_fd);
+ G_free(buffer);
+
+ cseg->name = G_store(map_name);
+ cseg->mapset = G_store(mapset);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_read.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/cseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/cseg_write.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/cseg_write.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,34 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+static char *me = "cseg_write_cell";
+
+int cseg_write_cellfile(CSEG * cseg, char *map_name)
+{
+ int map_fd;
+ int row, nrows;
+ CELL *buffer;
+
+ map_fd = G_open_cell_new(map_name);
+ if (map_fd < 0) {
+ G_warning("%s(): unable to open new map layer [%s]", me, map_name);
+ return -1;
+ }
+ nrows = G_window_rows();
+ buffer = G_allocate_cell_buf();
+ segment_flush(&(cseg->seg));
+ for (row = 0; row < nrows; row++) {
+ segment_get_row(&(cseg->seg), buffer, row);
+ if (G_put_raster_row(map_fd, buffer, CELL_TYPE) < 0) {
+ G_free(buffer);
+ G_unopen_cell(map_fd);
+ G_warning("%s(): unable to write new map layer [%s], row %d",
+ me, map_name, row);
+ return -2;
+ }
+ }
+ G_free(buffer);
+ G_close_cell(map_fd);
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/cseg_write.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/def_basin.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/def_basin.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/def_basin.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,105 @@
+#include "Gwater.h"
+
+CELL
+def_basin(int row, int col, CELL basin_num, double stream_length,
+ CELL old_elev)
+{
+ int r, rr, c, cc, ct, new_r[9], new_c[9];
+ CELL downdir, direction, asp_value, value, new_elev;
+ SHORT oldupdir, riteflag, leftflag, thisdir;
+
+ for (;;) {
+ cseg_put(&bas, &basin_num, row, col);
+ bseg_put(&swale, &one, row, col);
+ cseg_get(&asp, &asp_value, row, col);
+ if (asp_value < 0)
+ asp_value = -asp_value;
+ ct = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ cseg_get(&asp, &value, r, c);
+ if (value < -1)
+ value = -value;
+ if (value == drain[rr][cc]) {
+ bseg_get(&swale, &value, r, c);
+ if (value) {
+ new_r[++ct] = r;
+ new_c[ct] = c;
+ }
+ }
+ }
+ }
+ }
+ if (ct == 0) {
+ no_stream(row, col, basin_num, stream_length, old_elev);
+ return (basin_num);
+ }
+ if (ct >= 2) {
+ basin_num = split_stream(row, col, new_r, new_c, ct,
+ basin_num, stream_length, old_elev);
+ return (basin_num);
+ }
+ oldupdir = drain[row - new_r[1] + 1][col - new_c[1] + 1];
+ cseg_get(&asp, &downdir, row, col);
+ if (downdir < 0)
+ downdir = -downdir;
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ cseg_get(&asp, &direction, r, c);
+ if (direction == drain[rr][cc]) {
+ thisdir = updrain[rr][cc];
+ switch (haf_basin_side
+ (oldupdir, (SHORT) downdir, thisdir)) {
+ case LEFT:
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_elev);
+ leftflag++;
+ break;
+ case RITE:
+ overland_cells(r, c, basin_num, basin_num,
+ &new_elev);
+ riteflag++;
+ break;
+ }
+ }
+ }
+ }
+ }
+ if (leftflag > riteflag) {
+ value = basin_num - 1;
+ cseg_put(&haf, &value, row, col);
+ }
+ else {
+ cseg_put(&haf, &basin_num, row, col);
+ }
+ if (sides == 8) {
+ if (new_r[1] != row && new_c[1] != col)
+ stream_length += diag;
+ else if (new_r[1] != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += window.ew_res;
+ }
+ else { /* sides == 4 */
+
+ if (asp_value == 2 || asp_value == 6) {
+ if (new_r[1] != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += diag;
+ }
+ else { /* asp_value == 4, 8 */
+
+ if (new_c[1] != col)
+ stream_length += window.ew_res;
+ else
+ stream_length += diag;
+ }
+ }
+ row = new_r[1];
+ col = new_c[1];
+ }
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/def_basin.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/do_astar.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/do_astar.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,306 @@
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "Gwater.h"
+#include "do_astar.h"
+
+int sift_up(int, CELL);
+
+int do_astar(void)
+{
+ POINT point;
+ int doer, count;
+ SHORT upr, upc, r, c, ct_dir;
+ CELL work_val, alt_val, alt_up, asp_up, wat_val;
+ CELL in_val, drain_val;
+ HEAP heap_pos;
+
+ /* double slope; */
+
+ G_message(_("SECTION 2: A * Search."));
+
+ count = 0;
+ seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+ first_astar = heap_pos.point;
+
+ /* A * Search: downhill path through all not masked cells */
+ while (first_astar != -1) {
+ G_percent(count++, do_points, 1);
+
+ seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+ doer = heap_pos.point;
+
+ seg_get(&astar_pts, (char *)&point, 0, doer);
+
+ /* drop astar_pts[doer] from heap */
+ drop_pt();
+
+ seg_get(&heap_index, (char *)&heap_pos, 0, 1);
+ first_astar = heap_pos.point;
+
+ point.nxt = first_cum;
+ seg_put(&astar_pts, (char *)&point, 0, doer);
+
+ first_cum = doer;
+ r = point.r;
+ c = point.c;
+
+ bseg_put(&worked, &one, r, c);
+ cseg_get(&alt, &alt_val, r, c);
+
+ /* check all neighbours, breadth first search */
+ for (ct_dir = 0; ct_dir < sides; ct_dir++) {
+ /* get r, c (upr, upc) for this neighbour */
+ upr = r + nextdr[ct_dir];
+ upc = c + nextdc[ct_dir];
+ /* check that upr, upc are within region */
+ if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
+ /* check if neighbour is in the list */
+ /* if not, add as new point */
+ bseg_get(&in_list, &in_val, upr, upc);
+ if (in_val == 0) {
+ cseg_get(&alt, &alt_up, upr, upc);
+ add_pt(upr, upc, r, c, alt_up, alt_val);
+ /* flow direction is set here */
+ drain_val = drain[upr - r + 1][upc - c + 1];
+ cseg_put(&asp, &drain_val, upr, upc);
+ }
+ else {
+ /* check if neighbour has not been worked on,
+ * update values for asp, wat and slp */
+ bseg_get(&worked, &work_val, upr, upc);
+ if (!work_val) {
+ cseg_get(&asp, &asp_up, upr, upc);
+ if (asp_up < -1) {
+ drain_val = drain[upr - r + 1][upc - c + 1];
+ cseg_put(&asp, &drain_val, upr, upc);
+ cseg_get(&wat, &wat_val, r, c);
+ if (wat_val > 0)
+ wat_val = -wat_val;
+ cseg_put(&wat, &wat_val, r, c);
+ cseg_get(&alt, &alt_up, upr, upc);
+ replace(upr, upc, r, c); /* alt_up used to be */
+ /* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
+ dseg_put (&slp, &slope, upr, upc); */
+ }
+ }
+ }
+ }
+ }
+ }
+ bseg_close(&worked);
+ bseg_close(&in_list);
+ seg_close(&heap_index);
+
+ G_percent(count, do_points, 1); /* finish it */
+ return 0;
+}
+
+/* new add point routine for min heap */
+int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+{
+ POINT point;
+ HEAP heap_pos;
+
+ /* double slp_value; */
+
+ /* slp_value = get_slope(r, c, downr, downc, ele, downe);
+ dseg_put (&slp, &slp_value, r, c); */
+ bseg_put(&in_list, &one, r, c);
+
+ /* add point to next free position */
+
+ heap_size++;
+
+ if (heap_size > do_points)
+ G_fatal_error(_("heapsize too large"));
+
+ heap_pos.point = nxt_avail_pt;
+ heap_pos.ele = ele;
+ seg_put(&heap_index, (char *)&heap_pos, 0, heap_size);
+
+ point.r = r;
+ point.c = c;
+ point.downr = downr;
+ point.downc = downc;
+
+ seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
+
+ nxt_avail_pt++;
+
+ /* sift up: move new point towards top of heap */
+
+ sift_up(heap_size, ele);
+
+ return 0;
+}
+
+/* drop point routine for min heap */
+int drop_pt(void)
+{
+ int child, childr, parent;
+ int childp, childrp, parentp;
+ CELL ele, eler;
+ int i;
+ POINT point;
+ HEAP heap_pos;
+
+ if (heap_size == 1) {
+ parent = -1;
+ heap_pos.point = -1;
+ heap_pos.ele = 0;
+ seg_put(&heap_index, (char *)&heap_pos, 0, 1);
+ heap_size = 0;
+ return 0;
+ }
+
+ /* start with heap root */
+ parent = 1;
+
+ /* sift down: move hole back towards bottom of heap */
+ /* sift-down routine customised for A * Search logic */
+
+ while ((child = GET_CHILD(parent)) <= heap_size) {
+ /* select child with lower ele, if both are equal, older child
+ * older child is older startpoint for flow path, important */
+ seg_get(&heap_index, (char *)&heap_pos, 0, child);
+ childp = heap_pos.point;
+ ele = heap_pos.ele;
+ if (child < heap_size) {
+ childr = child + 1;
+ i = 1;
+ while (childr <= heap_size && i < 3) {
+ seg_get(&heap_index, (char *)&heap_pos, 0, childr);
+ childrp = heap_pos.point;
+ eler = heap_pos.ele;
+ if (eler < ele) {
+ child = childr;
+ childp = childrp;
+ ele = eler;
+ }
+ /* make sure we get the oldest child */
+ else if (ele == eler && childp > childrp) {
+ child = childr;
+ childp = childrp;
+ }
+ childr++;
+ i++;
+ }
+ }
+
+ /* move hole down */
+ heap_pos.point = childp;
+ heap_pos.ele = ele;
+ seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+ parent = child;
+
+ }
+
+ /* hole is in lowest layer, move to heap end */
+ if (parent < heap_size) {
+ seg_get(&heap_index, (char *)&heap_pos, 0, heap_size);
+ seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+
+ ele = heap_pos.ele;
+
+ /* sift up last swapped point, only necessary if hole moved to heap end */
+ sift_up(parent, ele);
+
+ }
+
+ /* the actual drop */
+ heap_size--;
+
+ return 0;
+
+}
+
+/* standard sift-up routine for d-ary min heap */
+int sift_up(int start, CELL ele)
+{
+ int parent, parentp, child, childp;
+ POINT point;
+ CELL elep;
+ HEAP heap_pos;
+
+ child = start;
+ seg_get(&heap_index, (char *)&heap_pos, 0, child);
+ childp = heap_pos.point;
+
+ while (child > 1) {
+ parent = GET_PARENT(child);
+ seg_get(&heap_index, (char *)&heap_pos, 0, parent);
+ parentp = heap_pos.point;
+ elep = heap_pos.ele;
+
+ /* parent ele higher */
+ if (elep > ele) {
+
+ /* push parent point down */
+ seg_put(&heap_index, (char *)&heap_pos, 0, child);
+ child = parent;
+
+ }
+ /* same ele, but parent is younger */
+ else if (elep == ele && parentp > childp) {
+
+ /* push parent point down */
+ seg_put(&heap_index, (char *)&heap_pos, 0, child);
+ child = parent;
+
+ }
+ else
+ /* no more sifting up, found new slot for child */
+ break;
+ }
+
+ /* no more sifting up, found new slot for child */
+ if (child < start) {
+ heap_pos.point = childp;
+ heap_pos.ele = ele;
+ seg_put(&heap_index, (char *)&heap_pos, 0, child);
+ }
+ return 0;
+}
+
+double
+get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+{
+ double slope;
+
+ if (r == downr)
+ slope = (ele - downe) / window.ew_res;
+ else if (c == downc)
+ slope = (ele - downe) / window.ns_res;
+ else
+ slope = (ele - downe) / diag;
+ if (slope < MIN_SLOPE)
+ return (MIN_SLOPE);
+ return (slope);
+}
+
+int replace( /* ele was in there */
+ SHORT upr, SHORT upc, SHORT r, SHORT c)
+/* CELL ele; */
+{
+ int now, heap_run;
+ POINT point;
+ HEAP heap_pos;
+
+ heap_run = 0;
+
+ while (heap_run <= heap_size) {
+ seg_get(&heap_index, (char *)&heap_pos, 0, heap_run);
+ now = heap_pos.point;
+ seg_get(&astar_pts, (char *)&point, 0, now);
+ if (point.r == upr && point.c == upc) {
+ point.downr = r;
+ point.downc = c;
+ seg_put(&astar_pts, (char *)&point, 0, now);
+ return 0;
+ }
+ heap_run++;;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/do_astar.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed2/seg/do_astar.h (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/do_astar.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,7 @@
+#ifndef __DO_ASTAR_H__
+#define __DO_ASTAR__
+
+#define GET_PARENT(c) ((int) (((c) - 2) / 3 + 1))
+#define GET_CHILD(p) ((int) ((p) * 3 - 1))
+
+#endif /* __DO_ASTAR_H__ */
Property changes on: grass/trunk/raster/r.watershed2/seg/do_astar.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/do_cum.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/do_cum.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,61 @@
+#include "Gwater.h"
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int do_cum(void)
+{
+ SHORT r, c, dr, dc;
+ CELL is_swale, value, valued;
+ POINT point;
+ int killer, threshold, count;
+
+ G_message(_("SECTION 3: Accumulating Surface Flow."));
+
+ count = 0;
+ if (bas_thres <= 0)
+ threshold = 60;
+ else
+ threshold = bas_thres;
+ while (first_cum != -1) {
+ G_percent(count++, do_points, 3);
+ killer = first_cum;
+ seg_get(&astar_pts, (char *)&point, 0, killer);
+ first_cum = point.nxt;
+ if ((dr = point.downr) > -1) {
+ r = point.r;
+ c = point.c;
+ dc = point.downc;
+ cseg_get(&wat, &value, r, c);
+ if (ABS(value) >= threshold)
+ bseg_put(&swale, &one, r, c);
+ cseg_get(&wat, &valued, dr, dc);
+ if (value > 0) {
+ if (valued > 0)
+ valued += value;
+ else
+ valued -= value;
+ }
+ else {
+ if (valued < 0)
+ valued += value;
+ else
+ valued = value - valued;
+ }
+ cseg_put(&wat, &valued, dr, dc);
+ bseg_get(&swale, &is_swale, r, c);
+ if (is_swale || ABS(valued) >= threshold) {
+ bseg_put(&swale, &one, dr, dc);
+ }
+ else {
+ if (er_flag && !is_swale)
+ slope_length(r, c, dr, dc);
+ }
+ }
+ }
+ seg_close(&astar_pts);
+
+ G_percent(count, do_points, 3); /* finish it */
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/do_cum.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_close.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_close.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_close.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,19 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include "cseg.h"
+
+int dseg_close(DSEG * dseg)
+{
+ segment_release(&(dseg->seg));
+ close(dseg->fd);
+ unlink(dseg->filename);
+ if (dseg->name) {
+ G_free(dseg->name);
+ dseg->name = NULL;
+ }
+ if (dseg->mapset) {
+ G_free(dseg->mapset);
+ dseg->mapset = NULL;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_close.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_get.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_get.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int dseg_get(DSEG * dseg, double *value, int row, int col)
+{
+ if (segment_get(&(dseg->seg), (CELL *) value, row, col) < 0) {
+ G_warning("dseg_get(): could not read segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_get.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_open.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_open.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_open.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,60 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+
+int dseg_open(DSEG * dseg, int srows, int scols, int nsegs_in_memory)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ dseg->filename = NULL;
+ dseg->fd = -1;
+ dseg->name = NULL;
+ dseg->mapset = NULL;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning("dseg_open(): unable to create segment file");
+ return -2;
+ }
+ if (0 >
+ (errflag =
+ segment_format(fd, G_window_rows(), G_window_cols(), srows, scols,
+ sizeof(double)))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("dseg_open(): could not write segment file");
+ return -1;
+ }
+ else {
+ G_warning("dseg_open(): illegal configuration parameter(s)");
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning("dseg_open(): unable to re-open segment file");
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(dseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("dseg_open(): could not read segment file");
+ return -5;
+ }
+ else {
+ G_warning("dseg_open(): out of memory");
+ return -6;
+ }
+ }
+ dseg->filename = filename;
+ dseg->fd = fd;
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_put.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_put.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int dseg_put(DSEG * dseg, double *value, int row, int col)
+{
+ if (segment_put(&(dseg->seg), (CELL *) value, row, col) < 0) {
+ G_warning("dseg_put(): could not write segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_put.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_read.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_read.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,61 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+static char *me = "dseg_read_cell";
+
+int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
+{
+ int row, col, nrows, ncols;
+ int map_fd;
+ char msg[100];
+ CELL *buffer;
+ double *dbuffer;
+
+ dseg->name = NULL;
+ dseg->mapset = NULL;
+
+ map_fd = G_open_cell_old(map_name, mapset);
+ if (map_fd < 0) {
+ sprintf(msg, "%s(): unable to open file [%s] in [%s], %d",
+ me, map_name, mapset, map_fd);
+ G_warning(msg);
+ return -3;
+ }
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ buffer = G_allocate_cell_buf();
+ dbuffer = (double *)G_malloc(ncols * sizeof(double));
+ for (row = 0; row < nrows; row++) {
+ if (G_get_c_raster_row(map_fd, buffer, row) < 0) {
+ G_free(buffer);
+ G_free(dbuffer);
+ G_close_cell(map_fd);
+ sprintf(msg, "%s(): unable to read file [%s] in [%s], %d %d",
+ me, map_name, mapset, row, nrows);
+ G_warning(msg);
+ return -2;
+ }
+ for (col = ncols - 1; col >= 0; col--) {
+ dbuffer[col] = (double)buffer[col];
+ }
+ if (segment_put_row(&(dseg->seg), (CELL *) dbuffer, row) < 0) {
+ G_free(buffer);
+ G_free(dbuffer);
+ G_close_cell(map_fd);
+ sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+ me, map_name, mapset);
+ G_warning(msg);
+ return (-1);
+ }
+ }
+
+ G_close_cell(map_fd);
+ G_free(buffer);
+ G_free(dbuffer);
+
+ dseg->name = G_store(map_name);
+ dseg->mapset = G_store(mapset);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_read.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/dseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/dseg_write.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/dseg_write.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,42 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+static char *me = "dseg_write_cell";
+
+int dseg_write_cellfile(DSEG * dseg, char *map_name)
+{
+ int map_fd;
+ int row, col, nrows, ncols;
+ CELL *buffer;
+ double *dbuffer;
+
+ map_fd = G_open_cell_new(map_name);
+ if (map_fd < 0) {
+ G_warning("%s(): unable to open new map layer [%s]", me, map_name);
+ return -1;
+ }
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ buffer = G_allocate_cell_buf();
+ dbuffer = (double *)G_malloc(ncols * sizeof(double));
+ segment_flush(&(dseg->seg));
+ for (row = 0; row < nrows; row++) {
+ segment_get_row(&(dseg->seg), (CELL *) dbuffer, row);
+ for (col = ncols - 1; col >= 0; col--) {
+ buffer[col] = (CELL) (dbuffer[col] + 0.5);
+ }
+ if (G_put_raster_row(map_fd, buffer, CELL_TYPE) < 0) {
+ G_free(buffer);
+ G_free(dbuffer);
+ G_unopen_cell(map_fd);
+ G_warning("%s(): unable to write new map layer [%s], row %d",
+ me, map_name, row);
+ return -2;
+ }
+ }
+ G_free(buffer);
+ G_free(dbuffer);
+ G_close_cell(map_fd);
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/dseg_write.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/find_pour.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/find_pour.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,46 @@
+#include "Gwater.h"
+
+int find_pourpts(void)
+{
+ int row, col;
+ double easting, northing, stream_length;
+ CELL old_elev, basin_num, value, is_swale;
+
+ basin_num = 0;
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
+ northing = window.north - (row + .5) * window.ns_res;
+ for (col = 0; col < ncols; col++) {
+ /* cseg_get (&wat, &value, row, col);
+ if (value < 0)
+ {
+ value = -value;
+ } */
+ cseg_get(&asp, &value, row, col);
+ bseg_get(&swale, &is_swale, row, col);
+ if (value < 0 && is_swale > 0) {
+ basin_num += 2;
+ cseg_get(&alt, &old_elev, row, col);
+ if (arm_flag) {
+ easting = window.west + (col + .5) * window.ew_res;
+ fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
+ (int)basin_num, 0, row, col, easting, northing);
+ }
+ if (col == 0 || col == ncols - 1) {
+ stream_length = .5 * window.ew_res;
+ }
+ else if (row == 0 || row == nrows - 1) {
+ stream_length = .5 * window.ns_res;
+ }
+ else {
+ stream_length = 0.0;
+ }
+ basin_num =
+ def_basin(row, col, basin_num, stream_length, old_elev);
+ }
+ }
+ }
+ G_percent(nrows, nrows, 1); /* finish it */
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/find_pour.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/haf_side.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/haf_side.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/haf_side.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,18 @@
+#include "Gwater.h"
+
+int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+{
+ SHORT newup, newthis;
+
+ newup = updir - downdir;
+ if (newup < 0)
+ newup += 8;
+ newthis = thisdir - downdir;
+ if (newthis < 0)
+ newthis += 8;
+ if (newthis < newup)
+ return (LEFT);
+ if (newthis > newup)
+ return (RITE);
+ return (NEITHER);
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/haf_side.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/init_vars.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/init_vars.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,466 @@
+#include <stdlib.h>
+#include <unistd.h>
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int init_vars(int argc, char *argv[])
+{
+ SHORT r, c;
+ int fd, num_cseg_total, num_open_segs;
+ int seg_rows, seg_cols;
+ double segs_mb;
+ char *mb_opt;
+
+ /* int page_block, num_cseg; */
+ int max_bytes;
+ CELL *buf, alt_value, wat_value, asp_value, worked_value;
+ extern FILE *fopen();
+ char MASK_flag, *do_exist();
+
+ G_gisinit(argv[0]);
+ ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
+ ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
+ zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+ nxt_avail_pt = 0;
+ /* dep_flag = 0; */
+ max_length = dzero = 0.0;
+ ril_value = -1.0;
+ /* dep_slope = 0.0; */
+ max_bytes = 0;
+ sides = 8;
+ for (r = 1; r < argc; r++) {
+ if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
+ ele_flag++;
+ else if (sscanf(argv[r], "ac=%[^\n]", wat_name) == 1)
+ wat_flag++;
+ else if (sscanf(argv[r], "dr=%[^\n]", asp_name) == 1)
+ asp_flag++;
+ else if (sscanf(argv[r], "de=%[^\n]", pit_name) == 1)
+ pit_flag++;
+ else if (sscanf(argv[r], "t=%d", &bas_thres) == 1) ;
+ else if (sscanf(argv[r], "ms=%lf", &max_length) == 1) ;
+ else if (sscanf(argv[r], "ba=%[^\n]", bas_name) == 1)
+ bas_flag++;
+ else if (sscanf(argv[r], "se=%[^\n]", seg_name) == 1)
+ seg_flag++;
+ else if (sscanf(argv[r], "ha=%[^\n]", haf_name) == 1)
+ haf_flag++;
+ else if (sscanf(argv[r], "ov=%[^\n]", run_name) == 1)
+ run_flag++;
+ else if (sscanf(argv[r], "ar=%[^\n]", arm_name) == 1)
+ arm_flag++;
+ else if (sscanf(argv[r], "di=%[^\n]", dis_name) == 1)
+ dis_flag++;
+ else if (sscanf(argv[r], "sl=%[^\n]", sl_name) == 1)
+ sl_flag++;
+ else if (sscanf(argv[r], "S=%[^\n]", sg_name) == 1)
+ sg_flag++;
+ else if (sscanf(argv[r], "LS=%[^\n]", ls_name) == 1)
+ ls_flag++;
+ else if (sscanf(argv[r], "ob=%[^\n]", ob_name) == 1)
+ ob_flag++;
+ else if (sscanf(argv[r], "mb=%[^\n]", mb_opt) == 1) {
+ if (sscanf(mb_opt, "%lf", &segs_mb) == 0) {
+ segs_mb = 300;
+ }
+ }
+ else if (sscanf(argv[r], "r=%[^\n]", ril_name) == 1) {
+ if (sscanf(ril_name, "%lf", &ril_value) == 0) {
+ ril_value = -1.0;
+ ril_flag++;
+ }
+ }
+ /* else if (sscanf (argv[r], "sd=%[^\n]", dep_name) == 1) dep_flag++; */
+ else if (sscanf(argv[r], "-%d", &sides) == 1) {
+ if (sides != 4)
+ usage(argv[0]);
+ }
+ else
+ usage(argv[0]);
+ }
+ if ((ele_flag != 1)
+ ||
+ ((arm_flag == 1) &&
+ ((bas_thres <= 0) || ((haf_flag != 1) && (bas_flag != 1))))
+ ||
+ ((bas_thres <= 0) &&
+ ((bas_flag == 1) || (seg_flag == 1) || (haf_flag == 1) ||
+ (sl_flag == 1) || (sg_flag == 1) || (ls_flag == 1)))
+ )
+ usage(argv[0]);
+ tot_parts = 4;
+ if (ls_flag || sg_flag)
+ tot_parts++;
+ if (bas_thres > 0)
+ tot_parts++;
+
+ G_message(_("SECTION 1 beginning: Initiating Variables. %d sections total."),
+ tot_parts);
+
+ this_mapset = G_mapset();
+ if (asp_flag)
+ do_legal(asp_name);
+ if (bas_flag)
+ do_legal(bas_name);
+ if (seg_flag)
+ do_legal(seg_name);
+ if (haf_flag)
+ do_legal(haf_name);
+ if (sl_flag)
+ do_legal(sl_name);
+ if (sg_flag)
+ do_legal(sg_name);
+ if (ls_flag)
+ do_legal(ls_name);
+ if (sl_flag || sg_flag || ls_flag)
+ er_flag = 1;
+ ele_mapset = do_exist(ele_name);
+ /* for sd factor
+ if (dep_flag) {
+ if (sscanf (dep_name, "%lf", &dep_slope) != 1) {
+ dep_mapset = do_exist (dep_name);
+ dep_flag = -1;
+ }
+ }
+ */
+ G_get_set_window(&window);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ if (max_length <= dzero)
+ max_length = 10 * nrows * window.ns_res + 10 * ncols * window.ew_res;
+ if (window.ew_res < window.ns_res)
+ half_res = .5 * window.ew_res;
+ else
+ half_res = .5 * window.ns_res;
+ diag = sqrt(window.ew_res * window.ew_res +
+ window.ns_res * window.ns_res);
+ if (sides == 4)
+ diag *= 0.5;
+
+ /* segment parameters: one size fits all. Fine tune? */
+ /* Segment rows and cols: 200 */
+ /* 1 segment open for all rasters: 2.86 MB */
+ /* num_open_segs = segs_mb / 2.86 */
+
+ seg_rows = SROW;
+ seg_cols = SCOL;
+
+ if (segs_mb < 3.0) {
+ segs_mb = 300;
+ G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
+ }
+
+ num_open_segs = segs_mb / 2.86;
+
+ G_debug(1, "segs MB: %.0f", segs_mb);
+ G_debug(1, "seg cols: %d", seg_cols);
+ G_debug(1, "seg rows: %d", seg_rows);
+
+ num_cseg_total = nrows / SROW + 1;
+ G_debug(1, " row segments:\t%d", num_cseg_total);
+
+ num_cseg_total = ncols / SCOL + 1;
+ G_debug(1, "column segments:\t%d", num_cseg_total);
+
+ num_cseg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
+ G_debug(1, " total segments:\t%d", num_cseg_total);
+ G_debug(1, " open segments:\t%d", num_open_segs);
+
+ /* nonsense to have more segments open than exist */
+ if (num_open_segs > nrows)
+ num_open_segs = nrows;
+ G_debug(1, " open segments after adjusting:\t%d", num_open_segs);
+
+ cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
+ cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
+ cseg_read_cell(&alt, ele_name, ele_mapset);
+ cseg_read_cell(&r_h, ele_name, ele_mapset);
+ cseg_open(&wat, seg_rows, seg_cols, num_open_segs);
+
+ if (run_flag) {
+ run_mapset = do_exist(run_name);
+ cseg_read_cell(&wat, run_name, run_mapset);
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++)
+ if (-1 == cseg_put(&wat, &one, r, c))
+ exit(EXIT_FAILURE);
+ }
+ }
+ cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+ if (pit_flag) {
+ pit_mapset = do_exist(pit_name);
+ cseg_read_cell(&asp, pit_name, pit_mapset);
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++)
+ if (-1 == cseg_put(&asp, &zero, r, c))
+ exit(EXIT_FAILURE);
+ }
+ }
+ bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
+ if (ob_flag) {
+ ob_mapset = do_exist(ob_name);
+ bseg_read_cell(&swale, ob_name, ob_mapset);
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++)
+ bseg_put(&swale, &zero, r, c);
+ }
+ }
+ if (ril_flag) {
+ ril_mapset = do_exist(ril_name);
+ dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
+ dseg_read_cell(&ril, ril_name, ril_mapset);
+ }
+ bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+ bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
+ MASK_flag = 0;
+ do_points = nrows * ncols;
+ if (NULL != G_find_file("cell", "MASK", G_mapset())) {
+ MASK_flag = 1;
+ if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
+ G_fatal_error(_("unable to open MASK"));
+ }
+ else {
+ buf = G_allocate_cell_buf();
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row_nomask(fd, buf, r);
+ for (c = 0; c < ncols; c++) {
+ if (!buf[c]) {
+ do_points--;
+ bseg_put(&worked, &one, r, c);
+ bseg_put(&in_list, &one, r, c);
+ }
+ }
+ }
+ G_close_cell(fd);
+ G_free(buf);
+ }
+ }
+ /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
+ dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
+ if (sg_flag)
+ dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
+ if (ls_flag)
+ dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
+ seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols,
+ num_open_segs, sizeof(POINT));
+
+ /* heap_index will track astar_pts in the binary min-heap */
+ /* heap_index is one-based */
+ seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
+ num_open_segs, sizeof(HEAP));
+
+ G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
+
+ /* heap is empty */
+ heap_size = 0;
+
+ first_astar = first_cum = -1;
+
+ if (MASK_flag) {
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 3);
+ for (c = 0; c < ncols; c++) {
+ bseg_get(&worked, &worked_value, r, c);
+ if (worked_value) {
+ cseg_put(&wat, &zero, r, c);
+ }
+ else {
+ dseg_put(&s_l, &half_res, r, c);
+ cseg_get(&asp, &asp_value, r, c);
+ if (r == 0 || c == 0 || r == nrows - 1 ||
+ c == ncols - 1 || asp_value != 0) {
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ if (r == 0)
+ asp_value = -2;
+ else if (c == 0)
+ asp_value = -4;
+ else if (r == nrows - 1)
+ asp_value = -6;
+ else if (c == ncols - 1)
+ asp_value = -8;
+ else
+ asp_value = -1;
+ if (-1 == cseg_put(&asp, &asp_value, r, c))
+ exit(EXIT_FAILURE);
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ }
+ else if (!bseg_get(&worked, &worked_value, r - 1, c)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -2;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (!bseg_get(&worked, &worked_value, r + 1, c)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -6;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (!bseg_get(&worked, &worked_value, r, c - 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -4;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (!bseg_get(&worked, &worked_value, r, c + 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -8;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (sides == 8 &&
+ !bseg_get(&worked, &worked_value, r - 1, c - 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -3;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (sides == 8 &&
+ !bseg_get(&worked, &worked_value, r - 1, c + 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -1;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (sides == 8 &&
+ !bseg_get(&worked, &worked_value, r + 1, c - 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -5;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else if (sides == 8 &&
+ !bseg_get(&worked, &worked_value, r + 1, c + 1)
+ && worked_value != 0) {
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ asp_value = -7;
+ cseg_put(&asp, &asp_value, r, c);
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ cseg_put(&wat, &wat_value, r, c);
+ }
+ }
+ else {
+ bseg_put(&in_list, &zero, r, c);
+ /* dseg_put(&slp, &dzero, r, c); */
+ }
+ }
+ }
+ }
+ }
+ else {
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 3);
+ for (c = 0; c < ncols; c++) {
+ bseg_put(&worked, &zero, r, c);
+ dseg_put(&s_l, &half_res, r, c);
+ cseg_get(&asp, &asp_value, r, c);
+ if (r == 0 || c == 0 || r == nrows - 1 ||
+ c == ncols - 1 || asp_value != 0) {
+ cseg_get(&wat, &wat_value, r, c);
+ if (wat_value > 0) {
+ wat_value = -wat_value;
+ if (-1 == cseg_put(&wat, &wat_value, r, c))
+ exit(EXIT_FAILURE);
+ }
+ if (r == 0)
+ asp_value = -2;
+ else if (c == 0)
+ asp_value = -4;
+ else if (r == nrows - 1)
+ asp_value = -6;
+ else if (c == ncols - 1)
+ asp_value = -8;
+ else
+ asp_value = -1;
+ if (-1 == cseg_put(&asp, &asp_value, r, c))
+ exit(EXIT_FAILURE);
+ cseg_get(&alt, &alt_value, r, c);
+ add_pt(r, c, -1, -1, alt_value, alt_value);
+ }
+ else {
+ bseg_put(&in_list, &zero, r, c);
+ /* dseg_put(&slp, &dzero, r, c); */
+ }
+ }
+ }
+ }
+ G_percent(r, nrows, 3); /* finish it */
+
+ return 0;
+}
+
+int do_legal(char *file_name)
+{
+ if (G_legal_filename(file_name) == -1)
+ G_fatal_error(_("map layer [%s] not legal for GRASS"), file_name);
+
+ return 0;
+}
+
+char *do_exist(char *file_name)
+{
+ char *file_mapset = G_find_cell2(file_name, "");
+
+ if (file_mapset == NULL)
+ G_fatal_error(_("[%s] map not found"), file_name);
+
+ return (file_mapset);
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/init_vars.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/main.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/main.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,59 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.watershed/seg - uses the GRASS segmentation library
+ * AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
+ * Markus Neteler <neteler itc.it>,
+ * Roberto Flor <flor itc.it>,
+ * Brad Douglas <rez touchofmadness.com>,
+ * Hamish Bowman <hamish_nospam yahoo.com>
+ * Markus Metz <markus.metz.giswork gmail.com>
+ * PURPOSE: Watershed determination using the GRASS segmentation lib
+ * COPYRIGHT: (C) 1999-2006 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 <stdlib.h>
+#include <unistd.h>
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#undef MAIN
+
+
+int main(int argc, char *argv[])
+{
+ extern FILE *fopen();
+
+ one = 1;
+ zero = 0;
+ dzero = 0.0;
+ init_vars(argc, argv);
+ do_astar();
+ do_cum();
+ if (sg_flag || ls_flag) {
+ sg_factor();
+ }
+
+ if (bas_thres <= 0) {
+ G_message(_("SECTION %d: Closing Maps."), tot_parts);
+ close_maps();
+ }
+ else {
+ if (arm_flag) {
+ fp = fopen(arm_name, "w");
+ }
+ cseg_open(&bas, SROW, SCOL, 4);
+ cseg_open(&haf, SROW, SCOL, 4);
+ G_message(_("SECTION %d: Watershed determination."), tot_parts - 1);
+ find_pourpts();
+ G_message(_("SECTION %d: Closing Maps."), tot_parts);
+ close_array_seg();
+ }
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/no_stream.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/no_stream.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,105 @@
+#include "Gwater.h"
+
+int
+no_stream(int row, int col, CELL basin_num, double stream_length,
+ CELL old_elev)
+{
+ int r, rr, c, cc, uprow = 0, upcol = 0;
+ double slope;
+ CELL downdir, max_drain, value, asp_value, hih_ele, new_ele, aspect;
+ SHORT updir, riteflag, leftflag, thisdir;
+
+ while (1) {
+ max_drain = -1;
+ for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
+ for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+
+ cseg_get(&asp, &aspect, r, c);
+ if (aspect == drain[rr][cc]) {
+ cseg_get(&wat, &value, r, c);
+ if (value < 0)
+ value = -value;
+ if (value > max_drain) {
+ uprow = r;
+ upcol = c;
+ max_drain = value;
+ }
+ }
+ }
+ }
+ }
+ if (max_drain > -1) {
+ updir = drain[row - uprow + 1][col - upcol + 1];
+ cseg_get(&asp, &downdir, row, col);
+ if (downdir < 0)
+ downdir = -downdir;
+ if (sides == 8) {
+ if (uprow != row && upcol != col)
+ stream_length += diag;
+ else if (uprow != row)
+ stream_length += window.ns_res;
+ else
+ stream_length += window.ew_res;
+ }
+ else { /* sides == 4 */
+
+ cseg_get(&asp, &asp_value, uprow, upcol);
+ if (downdir == 2 || downdir == 6) {
+ if (asp_value == 2 || asp_value == 6)
+ stream_length += window.ns_res;
+ else
+ stream_length += diag;
+ }
+ else { /* downdir == 4,8 */
+
+ if (asp_value == 4 || asp_value == 8)
+ stream_length += window.ew_res;
+ else
+ stream_length += diag;
+ }
+ }
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ cseg_get(&asp, &aspect, r, c);
+ if (aspect == drain[rr][cc]) {
+ thisdir = updrain[rr][cc];
+ if (haf_basin_side(updir,
+ (SHORT) downdir,
+ thisdir) == RITE) {
+ overland_cells(r, c, basin_num, basin_num,
+ &new_ele);
+ riteflag++;
+ }
+ else {
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_ele);
+ leftflag++;
+ }
+ }
+ }
+ }
+ }
+ if (leftflag >= riteflag) {
+ value = basin_num - 1;
+ cseg_put(&haf, &value, row, col);
+ }
+ else
+ cseg_put(&haf, &basin_num, row, col);
+ row = uprow;
+ col = upcol;
+ }
+ else {
+ if (arm_flag) {
+ cseg_get(&alt, &hih_ele, row, col);
+ slope = (hih_ele - old_elev) / stream_length;
+ if (slope < MIN_SLOPE)
+ slope = MIN_SLOPE;
+ fprintf(fp, " %f %f\n", slope, stream_length);
+ }
+ return 0;
+ }
+ }
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/no_stream.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/over_cells.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/over_cells.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/over_cells.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,39 @@
+#include "Gwater.h"
+#define BIGNEG -9999999
+
+int
+overland_cells(int row, int col, CELL basin_num, CELL haf_num, CELL * hih_ele)
+{
+ int r, rr, c, cc;
+ CELL new_ele, new_max_ele, value;
+
+ cseg_put(&bas, &basin_num, row, col);
+ cseg_put(&haf, &haf_num, row, col);
+ new_max_ele = BIGNEG;
+ for (r = row - 1, rr = 0; r <= row + 1; r++, rr++) {
+ for (c = col - 1, cc = 0; c <= col + 1; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ cseg_get(&asp, &value, r, c);
+ if (value == drain[rr][cc]) {
+ if (r != row && c != col) {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ else if (r != row) {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ else {
+ overland_cells(r, c, basin_num, haf_num, &new_ele);
+ }
+ }
+ }
+ }
+ }
+ if (new_max_ele == BIGNEG) {
+ cseg_get(&alt, hih_ele, row, col);
+ }
+ else {
+ *hih_ele = new_max_ele;
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/over_cells.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/sg_factor.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/sg_factor.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,67 @@
+#include "Gwater.h"
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int sg_factor(void)
+{
+ int r, c;
+ CELL downer, low_elev, hih_elev;
+ double height, length, S, sin_theta;
+
+ G_message(_("SECTION 4: Length Slope determination."));
+ for (r = nrows - 1; r >= 0; r--) {
+ G_percent(r, nrows, 3);
+ for (c = ncols - 1; c >= 0; c--) {
+ cseg_get(&alt, &low_elev, r, c);
+ cseg_get(&r_h, &hih_elev, r, c);
+ dseg_get(&s_l, &length, r, c);
+ cseg_get(&asp, &downer, r, c);
+ height = hih_elev - low_elev;
+ if (length > max_length) {
+ height *= max_length / length;
+ length = max_length;
+ }
+ sin_theta = height / sqrt(height * height + length * length);
+ if (height / length < .09)
+ S = 10.8 * sin_theta + .03;
+ else
+ S = 16.8 * sin_theta - .50;
+ if (ls_flag) {
+ length *= METER_TO_FOOT;
+ len_slp_equ(length, sin_theta, S, r, c);
+ }
+ if (sg_flag) {
+ S *= 100.0;
+ dseg_put(&s_g, &S, r, c);
+ }
+ }
+ }
+ G_percent(r, nrows, 3); /* finish it */
+
+ return 0;
+}
+
+int len_slp_equ(double slope_length, double sin_theta, double S, int r, int c)
+{
+ double rill, s_l_exp, /* m */
+ rill_ratio, /* Beta */
+ LS;
+
+ rill_ratio = (sin_theta / 0.0896) / (3.0 * pow(sin_theta, 0.8) + 0.56);
+ if (ril_flag) {
+ dseg_get(&ril, &rill, r, c);
+ }
+ else if (ril_value >= 0.0) {
+ rill = ril_value;
+ }
+ else
+ rill = 0.0;
+ /* rill_ratio equation from Steve Warren */
+ rill_ratio *= .5 + .005 * rill + .0001 * rill * rill;
+ s_l_exp = rill_ratio / (1 + rill_ratio);
+ LS = S * 100 * pow((slope_length / 72.6), s_l_exp);
+ dseg_put(&l_s, &LS, r, c);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/sg_factor.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/slope_len.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/slope_len.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,52 @@
+#include "Gwater.h"
+
+int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+{
+ CELL top_alt, bot_alt, asp_value, ridge;
+ double res, top_ls, bot_ls;
+
+ if (sides == 8) {
+ if (r == dr)
+ res = window.ns_res;
+ else if (c == dc)
+ res = window.ew_res;
+ else
+ res = diag;
+ }
+ else { /* sides == 4 */
+
+ cseg_get(&asp, &asp_value, dr, dc);
+ if (r == dr) {
+ if (asp_value == 2 || asp_value == 6)
+ res = window.ns_res;
+ else /* asp_value == 4, 8, -2, -4, -6, or -8 */
+ res = diag;
+ }
+ else { /* c == dc */
+
+ if (asp_value == 4 || asp_value == 8)
+ res = window.ew_res;
+ else /* asp_value == 2, 6, -2, -4, -6, or -8 */
+ res = diag;
+ }
+ }
+ dseg_get(&s_l, &top_ls, r, c);
+ if (top_ls == half_res)
+ top_ls = res;
+ else
+ top_ls += res;
+ dseg_put(&s_l, &top_ls, r, c);
+ cseg_get(&alt, &top_alt, r, c);
+ cseg_get(&alt, &bot_alt, dr, dc);
+ if (top_alt > bot_alt) {
+ dseg_get(&s_l, &bot_ls, dr, dc);
+ if (top_ls > bot_ls) {
+ bot_ls = top_ls + res;
+ dseg_put(&s_l, &bot_ls, dr, dc);
+ cseg_get(&r_h, &ridge, r, c);
+ cseg_put(&r_h, &ridge, dr, dc);
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/slope_len.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/split_str.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/split_str.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/split_str.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,88 @@
+#include "Gwater.h"
+
+CELL
+split_stream(int row, int col, int new_r[], int new_c[], int ct,
+ CELL basin_num, double stream_length, CELL old_elev)
+{
+ CELL downdir, old_basin, new_elev, aspect;
+ double slope, easting, northing;
+ SHORT doit, ctr, updir, splitdir[9];
+ SHORT thisdir, leftflag, riteflag;
+ int r, c, rr, cc;
+
+ for (ctr = 1; ctr <= ct; ctr++)
+ splitdir[ctr] = drain[row - new_r[ctr] + 1][col - new_c[ctr] + 1];
+ updir = splitdir[1];
+ cseg_get(&asp, &downdir, row, col);
+ if (downdir < 0)
+ downdir = -downdir;
+ riteflag = leftflag = 0;
+ for (r = row - 1, rr = 0; rr < 3; r++, rr++) {
+ for (c = col - 1, cc = 0; cc < 3; c++, cc++) {
+ if (r >= 0 && c >= 0 && r < nrows && c < ncols) {
+ cseg_get(&asp, &aspect, r, c);
+ if (aspect == drain[rr][cc]) {
+ doit = 1;
+ thisdir = updrain[rr][cc];
+ for (ctr = 1; ctr <= ct; ctr++) {
+ if (thisdir == splitdir[ctr]) {
+ doit = 0;
+ ctr = ct;
+ }
+ }
+ if (doit) {
+ thisdir = updrain[rr][cc];
+ switch (haf_basin_side
+ (updir, (SHORT) downdir, thisdir)) {
+ case LEFT:
+ overland_cells(r, c, basin_num, basin_num - 1,
+ &new_elev);
+ leftflag++;
+ break;
+ case RITE:
+ overland_cells(r, c, basin_num, basin_num,
+ &new_elev);
+ riteflag++;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (leftflag >= riteflag) {
+ old_basin = basin_num - 1;
+ cseg_put(&haf, &old_basin, row, col);
+ }
+ else {
+ cseg_put(&haf, &basin_num, row, col);
+ }
+ old_basin = basin_num;
+ cseg_get(&alt, &new_elev, row, col);
+ if ((slope = (new_elev - old_elev) / stream_length) < MIN_SLOPE)
+ slope = MIN_SLOPE;
+ if (arm_flag)
+ fprintf(fp, " %f %f\n", slope, stream_length);
+ for (r = 1; r <= ct; r++) {
+ basin_num += 2;
+ easting = window.west + (new_c[r] + .5) * window.ew_res;
+ northing = window.north - (new_r[r] + .5) * window.ns_res;
+ if (arm_flag) {
+ fprintf(fp, "%5d drains into %5d at %3d %3d %.3f %.3f",
+ (int)basin_num, old_basin, new_r[r], new_c[r], easting,
+ northing);
+ }
+ if (new_r[r] != row && new_c[r] != col)
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, diag, new_elev);
+ else if (new_r[r] != row)
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, window.ns_res,
+ new_elev);
+ else
+ basin_num =
+ def_basin(new_r[r], new_c[r], basin_num, window.ew_res,
+ new_elev);
+ }
+ return (basin_num);
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/split_str.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/sseg_close.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/sseg_close.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/sseg_close.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int seg_close(SSEG * sseg)
+{
+ segment_release(&(sseg->seg));
+ close(sseg->fd);
+ unlink(sseg->filename);
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/sseg_close.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/sseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/sseg_get.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/sseg_get.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int seg_get(SSEG * sseg, char *value, int row, int col)
+{
+ if (segment_get(&(sseg->seg), (CELL *) value, row, col) < 0) {
+ G_warning("seg_get(): could not read segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/sseg_get.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/sseg_open.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/sseg_open.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/sseg_open.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,57 @@
+#include <grass/gis.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int
+seg_open(SSEG * sseg, int nrows, int ncols, int row_in_seg, int col_in_seg,
+ int nsegs_in_memory, int size_struct)
+{
+ char *filename;
+ int errflag;
+ int fd;
+
+ sseg->filename = NULL;
+ sseg->fd = -1;
+
+ filename = G_tempfile();
+ if (-1 == (fd = creat(filename, 0666))) {
+ G_warning("seg_open(): unable to create segment file");
+ return -2;
+ }
+ if (0 > (errflag = segment_format(fd, nrows, ncols,
+ row_in_seg, col_in_seg, size_struct))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("seg_open(): could not write segment file");
+ return -1;
+ }
+ else {
+ G_warning("seg_open(): illegal configuration parameter(s)");
+ return -3;
+ }
+ }
+ close(fd);
+ if (-1 == (fd = open(filename, 2))) {
+ unlink(filename);
+ G_warning("seg_open(): unable to re-open segment file");
+ return -4;
+ }
+ if (0 > (errflag = segment_init(&(sseg->seg), fd, nsegs_in_memory))) {
+ close(fd);
+ unlink(filename);
+ if (errflag == -1) {
+ G_warning("seg_open(): could not read segment file");
+ return -5;
+ }
+ else {
+ G_warning("seg_open(): out of memory");
+ return -6;
+ }
+ }
+ sseg->filename = filename;
+ sseg->fd = fd;
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/sseg_open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/sseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/sseg_put.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/sseg_put.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+#include <grass/segment.h>
+#include "cseg.h"
+
+int seg_put(SSEG * sseg, char *value, int row, int col)
+{
+ if (segment_put(&(sseg->seg), (CELL *) value, row, col) < 0) {
+ G_warning("seg_put(): could not write segment file");
+ return -1;
+ }
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/sseg_put.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/seg/usage.c
===================================================================
--- grass/trunk/raster/r.watershed2/seg/usage.c (rev 0)
+++ grass/trunk/raster/r.watershed2/seg/usage.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,25 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+void usage(char *me)
+{
+ G_fatal_error(_("USAGE for basin delineation:\n%s -4 el=elevation_map "
+ "t=swale_threshold [ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "ba=watershed_basin_map [se=stream_segment_map]\n\nUSAGE for "
+ "ARMSED FILE creation:\n%s [-4] el=elevation_map "
+ "t=swale_threshold [ov=overland_flow_map] [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "[ba=watershed_basin_map] [se=stream_segment_map] "
+ "ha=half_basin_map ar=ARMSED_file_name\n\nUSAGE for slope "
+ "length determination:\n%s [-4] el=elevation_map "
+ "t=swale_threshold [dr=drain_direction_map] "
+ "[de=depression_map] [ac=accumulation_map] [di=display_map] "
+ "[ms=max_slope_length] [ob=overland_blocking_map] "
+ "[S=slope_steepness_map] LS=length_slope_map "
+ "[r=rill_erosion_map] [sd=slope_deposition value or map]"),
+ me, me, me);
+}
Property changes on: grass/trunk/raster/r.watershed2/seg/usage.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/Makefile
===================================================================
--- grass/trunk/raster/r.watershed2/shed/Makefile (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/Makefile 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.watershed
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: inter
Property changes on: grass/trunk/raster/r.watershed2/shed/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/accum_down.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/accum_down.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/accum_down.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,55 @@
+#include "watershed.h"
+
+int accum_down(OUTPUT * output)
+{
+ double new_length;
+ int m, b, num_basins, down_basin;
+ B_FACTS *basin_facts, *basin_info, *down_basin_info;
+ CAT *down_cat, *cat;
+ MAP *map;
+
+ num_basins = output->num_basins;
+ basin_facts = output->basin_facts;
+ for (b = num_basins - 1; b >= 0; b--) {
+ basin_facts[b].accum_length = basin_facts[b].str_length;
+ basin_facts[b].accum_slope = basin_facts[b].str_slope;
+ }
+ for (b = num_basins - 1; b >= 0; b--) {
+ basin_info = &(basin_facts[b]);
+ down_basin = basin_info->down_basin;
+ if (down_basin >= 0) {
+ down_basin_info = &(basin_facts[down_basin]);
+ down_basin_info->num_cells += basin_info->num_cells;
+ new_length =
+ basin_info->accum_length + down_basin_info->str_length;
+ if (new_length > down_basin_info->accum_length) {
+ down_basin_info->accum_length = new_length;
+ down_basin_info->accum_slope =
+ (down_basin_info->str_slope *
+ down_basin_info->str_length +
+ basin_info->accum_slope * basin_info->accum_length) /
+ new_length;
+ }
+ /* accum map layer information */
+ for (m = 0; m < output->num_maps; m++) {
+ map = &(output->maps[m]);
+ map->basins[down_basin].sum_values +=
+ map->basins[b].sum_values;
+ if (output->maps[m].do_cats != 0) {
+ cat = &(map->basins[b].first_cat);
+ down_cat = &(map->basins[down_basin].first_cat);
+ while (cat != NULL) {
+ insert_cat(down_cat, cat->cat_val, cat->num_cat);
+ cat = cat->nxt;
+ }
+ }
+ }
+ }
+ }
+ for (b = num_basins - 1; b >= 0; b--) {
+ basin_facts[b].str_length = basin_facts[b].accum_length;
+ basin_facts[b].str_slope = basin_facts[b].accum_slope;
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/accum_down.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/basin_maps.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/basin_maps.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/basin_maps.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,78 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+
+
+int basin_maps(INPUT * input, OUTPUT * output)
+{
+ char *mapset, map_layer[48];
+ int i;
+
+ G_message(_("\n\nPlease indicate which map layers you wish to use in the lumped"));
+ G_message(_("parameter hydrologic/soil erosion model. Continue inputing cell map"));
+ G_message(_("layers, one at a time, until all desired map layers are in."));
+ G_message(_("You can have %s include a list of categories in each."),
+ G_program_name());
+ G_message(_("\nHit <return> at the map prompt to continue with %s"),
+ G_program_name());
+
+ mapset = G_ask_old("", map_layer, "cell", "cell");
+ while (mapset != NULL) {
+ output->num_maps++;
+ if (output->num_maps == 1)
+ output->maps = (MAP *) G_malloc(sizeof(MAP));
+ else
+ output->maps =
+ (MAP *) G_realloc(output->maps,
+ output->num_maps * sizeof(MAP));
+ output->maps[output->num_maps - 1].mapset = mapset;
+ output->maps[output->num_maps - 1].name = G_store(map_layer);
+ output->maps[output->num_maps - 1].do_cats =
+ G_yes("Complete list of categories?", 1);
+ mapset = G_ask_old("", map_layer, "cell", "cell");
+ }
+
+ G_message(_("\nThe output from %s will be divided into watershed"),
+ G_program_name());
+ G_message(_("basins. There are two possible methods of tabulating the information:"));
+ G_message(_("1) by only including data pertaining to the basin itself, or 2) using"));
+ G_message(_("data from the basin, and all basins upstream of it."));
+
+ do {
+ G_message(_("\nWould you like the data organized:"));
+ G_message(_("1) Basin only\n2) Upstream only\n3) Both\nOR 0) to cancel program"));
+ fprintf(stderr, _("\nYour choice: "));
+ G_gets(map_layer);
+ sscanf(map_layer, "%d", &i);
+ } while (i > 3 || i < 0);
+
+ switch (i) {
+ case 0:
+ exit(EXIT_SUCCESS);
+ break;
+ case 1:
+ output->do_basin = 1;
+ output->do_accum = 0;
+ break;
+ case 2:
+ output->do_basin = 0;
+ output->do_accum = 1;
+ break;
+ case 3:
+ output->do_basin = 1;
+ output->do_accum = 1;
+ break;
+ }
+
+ if (input->fast) {
+ G_message(_("\nOK, %s should start running now using the "
+ "following form:\n%s"), RAM_NAME, input->com_line_ram);
+ }
+ else {
+ G_message(_("\nOK, %s should start running now using the "
+ "following form:\n%s"), SEG_NAME, input->com_line_seg);
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/basin_maps.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/com_line.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/com_line.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/com_line.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,456 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+#include "string.h"
+
+
+int com_line_Gwater(INPUT * input, OUTPUT * output)
+{
+ struct Cell_head *window;
+ char map_layer[48], buf[100], *prog_name, *mapset;
+ double d;
+ int i;
+
+ window = &(output->window);
+ if (0 == G_yes("Continue?", 1))
+ exit(EXIT_SUCCESS);
+ input->haf_name = (char *)G_calloc(40, sizeof(char));
+ input->accum_name = (char *)G_calloc(40, sizeof(char));
+
+ G_message(_("\nThis set of questions will organize the command line for the"));
+ G_message(_("%s program to run properly for your application."),
+ NON_NAME);
+ G_message(_("The first question is whether you want %s to run"),
+ NON_NAME);
+ G_message(_("in its fast mode or its slow mode. If you run %s"),
+ NON_NAME);
+ G_message(_("in the fast mode, the computer will finish about 10 times faster"));
+ G_message(_("than in the slow mode, but will not allow other programs to run"));
+ G_message(_("at the same time. The fast mode also places all of the data into"));
+ G_message(_("RAM, which limits the size of window that can be run. The slow"));
+ G_message(_("mode uses disk space in the same hard disk partition as where GRASS is"));
+ G_message(_("stored. Thus, if the program does not work in the slow mode, you will"));
+ G_message(_("need to remove unnecessary files from that partition. The slow mode"));
+ G_message(_("will allow other processes to run concurrently with %s.\n"),
+ NON_NAME);
+
+ sprintf(buf, "Do you want to use the fast mode of %s?", NON_NAME);
+ input->com_line_ram = input->com_line_seg = NULL;
+ input->fast = 0;
+ input->slow = 0;
+ if (G_yes(buf, 1)) {
+ input->fast = 1;
+ input->com_line_ram = (char *)G_calloc(400, sizeof(char));
+ prog_name = G_store(RAM_NAME);
+ sprintf(input->com_line_ram,
+ "%s/etc/water/%s", G_getenv("GISBASE"), RAM_NAME);
+ fprintf(stderr,
+ "\nIf there is not enough ram for the fast mode (%s) to run,\n",
+ RAM_NAME);
+ sprintf(buf, "should the slow mode (%s) be run instead?", SEG_NAME);
+ if (G_yes(buf, 1)) {
+ input->slow = 1;
+ input->com_line_seg = (char *)G_calloc(400, sizeof(char));
+ sprintf(input->com_line_seg,
+ "%s/etc/water/%s", G_getenv("GISBASE"), SEG_NAME);
+ }
+ }
+ else {
+ input->slow = 1;
+ prog_name = G_store(SEG_NAME);
+ input->com_line_seg = (char *)G_calloc(400, sizeof(char));
+ sprintf(input->com_line_seg,
+ "%s/etc/water/%s", G_getenv("GISBASE"), SEG_NAME);
+ }
+
+ G_message(_("\nIf you hit <return> by itself for the next question, this"));
+ G_message(_("program will terminate."));
+
+ mapset = G_ask_old("What is the name of the elevation map layer?",
+ map_layer, "cell", "cell");
+ if (!mapset)
+ exit(EXIT_FAILURE);
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " el=", map_layer, mapset);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " el=", map_layer, mapset);
+
+ G_message(_("\nOne of the options for %s is a `depression map'. A"),
+ prog_name);
+ G_message(_("depression map indicates all the locations in the current map window where"));
+ G_message(_("water accumulates and does not leave by the edge of the map. Lakes without"));
+ G_message(_("outlet streams and sinkholes are examples of `depressions'. If you wish to"));
+ G_message(_("have a depression map, prepare a map where non-zero values indicate the"));
+ G_message(_("locations where depressions occur.\n"));
+ G_message(_("Hit <return> by itself for the next question if there is no depression map."));
+
+ mapset = G_ask_old("What is the name of the depression map layer?",
+ map_layer, "cell", "cell");
+ if (mapset) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " de=", map_layer, mapset);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " de=", map_layer, mapset);
+ }
+
+ G_message(_("\nThe %s program will divide the elevation map into a number of"),
+ prog_name);
+ G_message(_("watershed basins. The number of watershed basins is indirectly determined"));
+ G_message(_("by the `basin threshold' value. The basin threshold is the area necessary for"));
+ G_message(_("%s to define a unique watershed basin. This area only applies to"),
+ prog_name);
+ G_message(_("`exterior drainage basins'. An exterior drainage basin does not have any"));
+ G_message(_("drainage basins flowing into it. Interior drainage basin size is determined"));
+ G_message(_("by the surface flow going into stream segments between stream interceptions."));
+ G_message(_("Thus interior drainage basins can be of any size. The %s program"),
+ prog_name);
+ G_message(_("also allows the user to relate basin size to potential overland flow"));
+ G_message(_("(i.e., areas with low infiltration capacities will need smaller areas to"));
+ G_message(_("develop stream channels than neighboring areas with high infiltration rates)."));
+ G_message(_("The user can create a map layer with potential overland flow values, and"));
+ G_message(_("%s will accumulate those values instead of area.\n"),
+ prog_name);
+ G_message(_("What unit of measure will you use for the basin threshold:"));
+
+ do {
+ G_message(_(" 1) acres, 2) meters sq., 3) miles sq., 4) hectares,"));
+ G_message(_(" 5) kilometers sq., 6) map cells, 7) overland flow units"));
+ fprintf(stderr, _("Choose 1-7 or 0 to exit this program: "));
+ G_gets(map_layer);
+ sscanf(map_layer, "%d", &i);
+ } while (i > 7 || i < 0);
+
+ if (!i)
+ exit(EXIT_SUCCESS);
+ output->type_area = (char)i;
+
+ G_message(_("\nHow large an area (or how many overland flow units) must a drainage basin"));
+ fprintf(stderr, _("be for it to be an exterior drainage basin: "));
+ G_gets(map_layer);
+ sscanf(map_layer, "%lf", &d);
+
+ switch (i) {
+ case 1:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d, ACRE_TO_METERSQ, window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d, ACRE_TO_METERSQ, window);
+ break;
+ case 2:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d, 1.0, window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d, 1.0, window);
+ break;
+ case 3:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d, MILESQ_TO_METERSQ,
+ window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d, MILESQ_TO_METERSQ,
+ window);
+ break;
+ case 4:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d, HECTACRE_TO_METERSQ,
+ window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d, HECTACRE_TO_METERSQ,
+ window);
+ break;
+ case 5:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d, KILOSQ_TO_METERSQ,
+ window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d, KILOSQ_TO_METERSQ,
+ window);
+ break;
+ case 6:
+ if (input->fast)
+ basin_com_add(&(input->com_line_ram), d,
+ (window->ns_res * window->ew_res), window);
+ if (input->slow)
+ basin_com_add(&(input->com_line_seg), d,
+ (window->ns_res * window->ew_res), window);
+ break;
+ case 7: /* needs an overland flow map */
+ G_message(_("\nIf you hit <return> by itself for the next question, this"));
+ G_message(_("program will terminate."));
+ mapset = G_ask_old("What is the name of the overland flow map layer?",
+ map_layer, "cell", "cell");
+ if (!mapset)
+ exit(EXIT_FAILURE);
+ if (input->fast) {
+ com_line_add(&(input->com_line_ram), " ov=", map_layer, mapset);
+ basin_com_add(&(input->com_line_ram), d,
+ (window->ns_res * window->ew_res), window);
+ }
+ if (input->slow) {
+ com_line_add(&(input->com_line_seg), " ov=", map_layer, mapset);
+ basin_com_add(&(input->com_line_seg), d,
+ (window->ns_res * window->ew_res), window);
+ }
+ break;
+ }
+
+ G_message(_("\n%s must create a map layer of watershed basins"),
+ prog_name);
+ G_message(_("before %s can run properly."), G_program_name());
+
+ strcpy(buf, "Please name the output watershed basin map:");
+ do {
+ mapset = G_ask_new(buf, input->haf_name, "cell", "");
+ } while (NULL == mapset);
+
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " ba=", input->haf_name, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " ba=", input->haf_name, NULL);
+
+ /*
+ This section queries the user about the armsed file input. If
+ you want to make this an option, the code below "COMMENT2" needs to be
+ modified.
+ */
+
+#ifdef ARMSED
+ G_message(_("\n%s must create a file of watershed basin relationships"),
+ prog_name);
+ G_message(_("before %s can run properly."), G_program_name());
+
+ input->ar_file_name = NULL;
+ while (input->ar_file_name == NULL) {
+ fprintf(stderr, _("\nPlease name this file:"));
+ G_gets(char_input);
+ if (1 != G_legal_filename(char_input)) {
+ G_message(_("<%s> is an illegal file name"), char_input);
+ }
+ else
+ input->ar_file_name = G_store(char_input);
+ }
+
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " ar=", input->ar_file_name,
+ NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " ar=", input->ar_file_name,
+ NULL);
+
+ /*
+ end of ARMSED comment code
+ */
+
+ /*
+ COMMENT2 This section of code tells the program where to place the statistics
+ about the watershed basin. GRASS users don't need this (w/ r.stats), but the
+ format is suppossed to be "user-friendly" to hydrologists. For the stats to be
+ created, the armsed file output needs to exist. For the stats to be an option
+ in this program: 1) it should be querried before the armsed file query, and 2)
+ make the armsed file query manditory if this option is invoked.
+ */
+
+ G_message(_("\n%s will generate a lot of output. Indicate a file"),
+ G_program_name());
+ G_message(_("name for %s to send the output to."), G_program_name());
+
+ output->file_name = NULL;
+ while (output->file_name == NULL) {
+ fprintf(stderr, _("\nPlease name this file:"));
+ G_gets(char_input);
+ if (1 != G_legal_filename(char_input)) {
+ G_message(_("<%s> is an illegal file name"), char_input);
+ }
+ else
+ output->file_name = G_store(char_input);
+ }
+
+ /*
+ end of COMMENT2
+ */
+#endif
+
+ G_message(_("\nThe accumulation map from %s must be present for"),
+ prog_name);
+ G_message(_("%s to work properly."), G_program_name());
+ strcpy(buf, "Please name the accumulation map:");
+ do {
+ mapset = G_ask_new(buf, input->accum_name, "cell", "");
+ } while (NULL == mapset);
+
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " ac=", input->accum_name, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " ac=", input->accum_name, NULL);
+
+ G_message(_("\n%s can produce several maps not necessary for"),
+ prog_name);
+ G_message(_("%s to function (stream channels, overland flow aspect, and"),
+ G_program_name());
+ G_message(_("a display version of the accumulation map). %s also has the"),
+ prog_name);
+ G_message(_("ability to generate several variables in the Revised Universal Soil Loss"));
+ G_message(_("Equation (Rusle): Slope Length (LS), and Slope Steepness (S).\n"));
+
+ sprintf(buf, "Would you like any of these maps to be created?");
+ if (G_yes(buf, 1)) {
+ mapset = G_ask_new("", map_layer, "cell", "stream channel");
+ if (mapset != NULL) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " se=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " se=", map_layer, NULL);
+ }
+ mapset = G_ask_new("", map_layer, "cell", "half basin");
+ if (mapset != NULL) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " ha=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " ha=", map_layer, NULL);
+ }
+ mapset = G_ask_new("", map_layer, "cell", "overland aspect");
+ if (mapset != NULL) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " dr=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " dr=", map_layer, NULL);
+ }
+ mapset = G_ask_new("", map_layer, "cell", "display");
+ if (mapset != NULL) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " di=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " di=", map_layer, NULL);
+ }
+ i = 0;
+ mapset = G_ask_new("", map_layer, "cell", "Slope Length");
+ if (mapset != NULL) {
+ i = 1;
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " LS=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " LS=", map_layer, NULL);
+ }
+ mapset = G_ask_new("", map_layer, "cell", "Slope Steepness");
+ if (mapset != NULL) {
+ i = 1;
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " S=", map_layer, NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " S=", map_layer, NULL);
+ }
+
+ if (i) {
+ G_message(_("\nThe Slope Length factor (LS) and Slope Steepness (S) are influenced by"));
+ G_message(_("disturbed land. %s reflects this with an optional map layer or value"),
+ prog_name);
+ G_message(_("where the value indicates the percent of disturbed (barren) land in that cell."));
+ G_message(_("Type <return> if you do not have a disturbed land map layer."));
+
+ mapset = G_ask_old("", map_layer, "cell", "disturbed land");
+ if (mapset != NULL) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " r=", map_layer,
+ NULL);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " r=", map_layer,
+ NULL);
+ }
+ else {
+ G_message(_("\nType the value indicating the percent of disturbed land. This value will"));
+ G_message(_("be used for every cell in the current region."));
+ i = -6;
+ while (i < 0 || i > 100) {
+ fprintf(stderr, _("\nInput value here [0-100]: "));
+ fgets(buf, 80, stdin);
+ sscanf(buf, "%d", &i);
+ }
+ if (input->fast)
+ com_add(&(input->com_line_ram), " r=", i);
+ if (input->slow)
+ com_add(&(input->com_line_seg), " r=", i);
+ }
+
+ /* 12345678901234567890123456789012345678901234567890123456789012345678901234567890 */
+ G_message(_("\nOverland surface flow only occurs for a set distance before swales form."));
+ G_message(_("Because of digital terrain model limitations, %s cannot pick up"),
+ prog_name);
+ G_message(_("these swales. %s allows for an input (warning: kludge factor)"),
+ prog_name);
+ G_message(_("that prevents the surface flow distance from getting too long. Normally,"));
+ G_message(_("maximum slope length is around 600 feet (about 183 meters)."));
+
+ i = -1;
+ while (i < 0) {
+ fprintf(stdout,
+ "\nInput maximum slope length here (in meters): ");
+ fgets(buf, 80, stdin);
+ sscanf(buf, "%d", &i);
+ }
+ if (input->fast)
+ com_add(&(input->com_line_ram), " ms=", i);
+ if (input->slow)
+ com_add(&(input->com_line_seg), " ms=", i);
+
+ /* 12345678901234567890123456789012345678901234567890123456789012345678901234567890 */
+ G_message(_("\nRoads, ditches, changes in ground cover, and other factors will stop"));
+ G_message(_("slope length. You may input a raster map indicating the locations of these"));
+ G_message(_("blocking factors.\n"));
+ G_message(_("Hit <return> by itself for the next question if there is no blocking map."));
+
+ mapset = G_ask_old("What is the name of the blocking map layer?",
+ map_layer, "cell", "cell");
+ if (mapset) {
+ if (input->fast)
+ com_line_add(&(input->com_line_ram), " ob=", map_layer,
+ mapset);
+ if (input->slow)
+ com_line_add(&(input->com_line_seg), " ob=", map_layer,
+ mapset);
+ }
+ }
+ }
+
+ return 0;
+}
+
+int com_line_add(char **com_line, char *prompt, char *map_layer, char *mapset)
+{
+ G_strcat(*com_line, prompt);
+ G_strcat(*com_line, "\"");
+ G_strcat(*com_line, map_layer);
+ if (mapset) {
+ G_strcat(*com_line, "@");
+ G_strcat(*com_line, mapset);
+ }
+ G_strcat(*com_line, "\"");
+
+ return 0;
+}
+
+int basin_com_add(char **com_line, double d, double modifier,
+ struct Cell_head *window)
+{
+ int i;
+ char buf[20];
+
+ i = (int)(.5 + modifier * d / window->ns_res / window->ew_res);
+ if (i < 1)
+ i = 1;
+ sprintf(buf, " t=%d", i);
+ G_strcat(*com_line, buf);
+
+ return 0;
+}
+
+int com_add(char **com_line, char *prompt, int ril_value)
+{
+ char buf[20];
+
+ G_strcat(*com_line, prompt);
+ sprintf(buf, "%d", ril_value);
+ G_strcat(*com_line, buf);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/com_line.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/file_in.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/file_in.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/file_in.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,42 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+
+
+int ar_file_in(char *file_name, OUTPUT * output)
+{
+ FILE *ar_file;
+ int bas_num, down_bas, i1, i2, bas_alloc;
+ char s1[7], s2[5], s3[3];
+ double northing, easting, str_slope, str_length;
+
+ if ((ar_file = fopen(file_name, "r")) == NULL) {
+ G_fatal_error(_("unable to open ARMSED file"));
+ }
+ output->num_basins = 0;
+ bas_alloc = INCR;
+ output->basin_facts = (B_FACTS *) G_malloc(INCR * sizeof(B_FACTS));
+ while (!feof(ar_file)) {
+ fscanf(ar_file, "%d %s %s %d %s %d %d %lf %lf %lf %lf",
+ &bas_num, s1, s2, &down_bas, s3, &i1, &i2,
+ &easting, &northing, &str_slope, &str_length);
+ if (output->num_basins >= bas_alloc) {
+ bas_alloc += INCR;
+ output->basin_facts =
+ (B_FACTS *) G_realloc(output->basin_facts,
+ bas_alloc * sizeof(B_FACTS));
+ }
+ output->basin_facts[output->num_basins].num_cells = 0;
+ output->basin_facts[output->num_basins].str_length = str_length;
+ output->basin_facts[output->num_basins].str_slope = str_slope;
+ output->basin_facts[output->num_basins].easting = easting;
+ output->basin_facts[output->num_basins].northing = northing;
+ output->basin_facts[output->num_basins].down_basin = down_bas / 2 - 1;
+ output->basin_facts[output->num_basins].valid = 1;
+ output->num_basins++;
+ }
+ fclose(ar_file);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/file_in.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/free.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/free.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/free.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,30 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include "watershed.h"
+
+int free_input(INPUT * input)
+{
+ if (input->com_line_ram)
+ G_free(input->com_line_ram);
+ if (input->com_line_seg)
+ G_free(input->com_line_seg);
+ G_free(input->haf_name);
+ G_free(input->ar_file_name);
+ G_free(input->accum_name);
+
+ return 0;
+}
+
+int free_output(OUTPUT * output)
+{
+ int c;
+
+ G_free(output->basin_facts);
+ G_free(output->file_name);
+ for (c = output->num_maps - 1; c >= 0; c--) {
+ G_free(output->maps[c].name);
+ }
+ G_free(output->maps);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/free.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/insert_cat.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/insert_cat.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/insert_cat.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,43 @@
+#include "watershed.h"
+
+int insert_cat(CAT * cat, CELL value, int num_cells)
+{
+ CAT *temp_cat;
+
+ for (;;) {
+ if (value > cat->cat_val) {
+ if (cat->nxt == NULL) {
+ cat->nxt = new_cat(value, num_cells);
+ cat = cat->nxt;
+ cat->nxt = NULL;
+ return 0;
+ }
+ else
+ cat = cat->nxt;
+ }
+ else if (value == cat->cat_val) {
+ cat->num_cat = cat->num_cat + num_cells;
+ return 0;
+ }
+ else {
+ temp_cat = new_cat(cat->cat_val, cat->num_cat);
+ temp_cat->nxt = cat->nxt;
+ cat->cat_val = value;
+ cat->num_cat = num_cells;
+ cat->nxt = temp_cat;
+ return 0;
+ }
+ }
+
+ return 0;
+}
+
+CAT *new_cat(CELL value, int num_cat)
+{
+ CAT *new;
+
+ new = (CAT *) G_malloc(sizeof(CAT));
+ new->num_cat = num_cat;
+ new->cat_val = value;
+ return (new);
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/insert_cat.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/intro.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/intro.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/intro.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,22 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+
+
+int intro(void)
+{
+ G_message(_("%s provides a text-based user-interface to the %s program."),
+ G_program_name(), NON_NAME);
+ G_message(_("%s also allows the user to prepare a report of map layers for each"),
+ G_program_name());
+ G_message(_("watershed basin determined in %s.\n"), NON_NAME);
+
+ G_message(_("%s will help the user determine which options to use for the"),
+ G_program_name());
+ G_message(_("%s program. %s will then ask for map layers that will be"),
+ NON_NAME, G_program_name());
+ G_message(_("divided by basin. %s will then run %s and create the report."),
+ G_program_name(), NON_NAME);
+
+ return (0);
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/intro.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/main.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/main.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/main.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,93 @@
+
+/****************************************************************************
+ *
+ * MODULE: shed
+ * AUTHOR(S): Charles Ehlschlaeger, CERL (original contributor)
+ * Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>,
+ * Brad Douglas <rez touchofmadness.com>, Hamish Bowman <hamish_nospam yahoo.com>
+ * PURPOSE: Watershed determination
+ * COPYRIGHT: (C) 1999-2006 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.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+#include "string.h"
+
+
+int main(int argc, char *argv[])
+{
+ INPUT input;
+ OUTPUT output;
+
+ G_gisinit(argv[0]);
+ G_set_program_name("r.watershed");
+ G_get_window(&(output.window));
+ intro();
+
+ output.num_maps = 0;
+ com_line_Gwater(&input, &output); /* develops r.watershed command line */
+ basin_maps(&input, &output); /* organizes map layers output */
+ if (input.fast || input.slow) {
+ if (input.fast) {
+ if (G_system(input.com_line_ram)) {
+ if (input.slow) {
+ G_message(_("Slow version of water analysis program starting now"));
+
+ if (G_system(input.com_line_seg)) {
+ free_input(&input);
+ free_output(&output);
+ G_fatal_error(_("<<%s>> command line failed"),
+ input.com_line_seg);
+ }
+ }
+ }
+ }
+ else if (G_system(input.com_line_seg)) {
+ free_input(&input);
+ free_output(&output);
+ G_fatal_error(_("<<%s>> command line failed"),
+ input.com_line_seg);
+ }
+ }
+
+ /*
+ ARMSED: This section exists to create the stats part.
+ input.ar_file_name could be used as a flag to determine this stuff
+ should be run.
+ */
+#ifdef ARMSED
+ ar_file_in(input.ar_file_name, &output);
+ read_basins(input.haf_name, &output);
+ valid_basins(input.accum_name, &output);
+ free_input(&input);
+ if ((output.out_file = fopen(output.file_name, "w")) == NULL) {
+ free_output(&output);
+ G_fatal_error(_("unable to open output file"));
+ }
+ if (output.do_basin) {
+ fprintf(output.out_file,
+ "\n\nThese values are accumulations within the basin itself\n");
+ fprintf(output.out_file, "They do not include sub-basins\n\n");
+ print_output(&output);
+ }
+ if (output.do_accum) {
+ accum_down(&output);
+ fprintf(output.out_file,
+ "\n\nThese values are accumulations of basins and sub-basins\n");
+ print_output(&output);
+ }
+ free_output(&output);
+
+#endif
+ /*
+ end of ARMSED comment code
+ */
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/print.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/print.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/print.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,107 @@
+#include "watershed.h"
+#include "string.h"
+
+int print_output(OUTPUT * output)
+{
+ double cell_size;
+ int b, c;
+ CAT *do_cat;
+ char *cat_name, area[32];
+
+ cell_size = output->window.ns_res * output->window.ew_res;
+ for (c = output->num_basins - 1; c >= 0; c--) {
+ if (output->basin_facts[c].valid == 1)
+ fprintf(output->out_file,
+ "\nValid Basin: %-5d flows into basin: %-5d at: E=%.1f N=%.1f\n",
+ (c + 1) * 2, (output->basin_facts[c].down_basin + 1) * 2,
+ output->basin_facts[c].easting,
+ output->basin_facts[c].northing);
+ else
+ fprintf(output->out_file,
+ "\nInvalid basin: %-5d flows into basin: %-5d at: E=%.1f N=%.1f\n",
+ (c + 1) * 2, (output->basin_facts[c].down_basin + 1) * 2,
+ output->basin_facts[c].easting,
+ output->basin_facts[c].northing);
+ fprintf(output->out_file,
+ " Str. length:%-.3f meters, %-.3f feet; Str. slope:%-.4f\n",
+ output->basin_facts[c].str_length,
+ (double)(output->basin_facts[c].str_length * METER_TO_FOOT),
+ output->basin_facts[c].str_slope);
+ switch (output->type_area) {
+ case 1:
+ fprintf(output->out_file, " Basin Area acres: %-16.4f",
+ output->basin_facts[c].num_cells * cell_size *
+ METERSQ_TO_ACRE);
+ break;
+ case 2:
+ fprintf(output->out_file, " Basin Area sq. meters: %-11.3f",
+ output->basin_facts[c].num_cells * cell_size);
+ break;
+ case 3:
+ fprintf(output->out_file, " Basin Area miles sq: %-16.5f",
+ output->basin_facts[c].num_cells * cell_size *
+ METERSQ_TO_MILESQ);
+ break;
+ case 4:
+ fprintf(output->out_file, " Basin Area hectareas: %-14.4f",
+ output->basin_facts[c].num_cells * cell_size *
+ METERSQ_TO_HECTACRE);
+ break;
+ case 5:
+ fprintf(output->out_file, " Basin Area kilometers: %-13.4f",
+ output->basin_facts[c].num_cells * cell_size *
+ METERSQ_TO_KILOSQ);
+ break;
+ case 6:
+ fprintf(output->out_file, " Basin Area in cells: %-16d",
+ output->basin_facts[c].num_cells);
+ break;
+ }
+ fprintf(output->out_file, " Area Percent Basin\n");
+ for (b = 0; b < output->num_maps; b++) {
+ fprintf(output->out_file,
+ "<< %20s >> map layer, average catagory value: %.2f\n",
+ output->maps[b].name,
+ ((double)output->maps[b].basins[c].sum_values) /
+ output->basin_facts[c].num_cells);
+ do_cat = &(output->maps[b].basins[c].first_cat);
+ while ((output->maps[b].do_cats != 0) && do_cat) {
+ cat_name =
+ G_get_cat(do_cat->cat_val, &(output->maps[b].cats));
+ switch (output->type_area) {
+ case 1:
+ sprintf(area, "%.3f acres",
+ METERSQ_TO_ACRE * cell_size * do_cat->num_cat);
+ break;
+ case 2:
+ sprintf(area, "%.2f sq. meters",
+ cell_size * do_cat->num_cat);
+ break;
+ case 3:
+ sprintf(area, "%.4f sq. miles",
+ METERSQ_TO_MILESQ * cell_size * do_cat->num_cat);
+ break;
+ case 4:
+ sprintf(area, "%.3f hectacres",
+ METERSQ_TO_HECTACRE * cell_size *
+ do_cat->num_cat);
+ break;
+ case 5:
+ sprintf(area, "%.3f sq. km.",
+ METERSQ_TO_KILOSQ * cell_size * do_cat->num_cat);
+ break;
+ case 6:
+ sprintf(area, "%6d cells", do_cat->num_cat);
+ break;
+ }
+ fprintf(output->out_file, "%3d %-43s %16s %-.4f\n",
+ do_cat->cat_val, cat_name, area,
+ ((double)do_cat->num_cat) /
+ output->basin_facts[c].num_cells);
+ do_cat = do_cat->nxt;
+ }
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/print.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/read.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/read.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/read.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,75 @@
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+
+
+int read_basins(char *haf_name, OUTPUT * output)
+{
+ int bas_fd, fd, m, r, nrows, c, ncols, tot_basins;
+ CELL v, *buf, *bas_buf, b;
+ CAT *cat;
+ MAP *map;
+ char *mapset;
+ B_FACTS *facts;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ buf = G_allocate_cell_buf();
+ bas_buf = G_allocate_cell_buf();
+ mapset = G_find_cell(haf_name, "");
+ if (!mapset) {
+ G_fatal_error(_("unable to open basin/half basin map"));
+ }
+
+ bas_fd = G_open_cell_old(haf_name, mapset);
+ facts = output->basin_facts;
+ for (r = nrows - 1; r >= 0; r--) {
+ G_get_c_raster_row(bas_fd, bas_buf, r);
+ for (c = ncols - 1; c >= 0; c--) {
+ b = bas_buf[c] / 2 - 1;
+ if (b >= 0)
+ facts[b].num_cells++;
+ }
+ }
+
+ tot_basins = output->num_basins;
+ for (m = 0; m < output->num_maps; m++) {
+ map = &(output->maps[m]);
+ G_read_cats(map->name, map->mapset, &(map->cats));
+ map->basins = (BASIN *) G_malloc(tot_basins * sizeof(BASIN));
+ for (r = 0; r < tot_basins; r++) {
+ map->basins[r].first_cat.num_cat = -1;
+ map->basins[r].first_cat.cat_val = -123456789;
+ map->basins[r].first_cat.nxt = NULL;
+ map->basins[r].sum_values = 0.0;
+ }
+ fd = G_open_cell_old(map->name, map->mapset);
+ if (fd >= 0) {
+ for (r = 0; r < nrows; r++) {
+ G_get_c_raster_row(fd, buf, r);
+ G_get_c_raster_row(bas_fd, bas_buf, r);
+ for (c = 0; c < ncols; c++) {
+ v = buf[c];
+ b = bas_buf[c] / 2 - 1;
+ if (b >= 0) {
+ map->basins[b].sum_values += v;
+ if (map->do_cats != 0) {
+ cat = &(map->basins[b].first_cat);
+ if (cat->num_cat == -1) {
+ cat->num_cat = 1;
+ cat->cat_val = v;
+ }
+ else
+ insert_cat(cat, v, (int)1);
+ }
+ }
+ }
+ }
+ G_close_cell(fd);
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/read.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/valid.c
===================================================================
--- grass/trunk/raster/r.watershed2/shed/valid.c (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/valid.c 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,62 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "watershed.h"
+
+
+int valid_basins(char *accum_name, OUTPUT * output)
+{
+ int i, r, c, fd;
+
+ /* int nrows, ncols; */
+ CELL *buf;
+ B_FACTS *basin_facts, *basin, *down_basin;
+ char *mapset;
+ struct Cell_head *window;
+
+ window = &(output->window);
+ /*
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+ */
+ if (NULL == (mapset = G_find_cell(accum_name, ""))) {
+ free_output(output);
+ G_fatal_error(_("accum file missing in valid_basins()"));
+ }
+ if (-1 == (fd = G_open_cell_old(accum_name, mapset))) {
+ free_output(output);
+ G_fatal_error(_("unable to open accum file in valid_basins()"));
+ }
+
+ buf = G_allocate_cell_buf();
+ basin_facts = output->basin_facts;
+ for (i = output->num_basins - 1; i >= 0; i--) {
+ basin = &(basin_facts[i]);
+ if (basin->valid == 0) {
+ if (basin->down_basin >= 0)
+ basin_facts[basin->down_basin].valid = 0;
+ }
+ else {
+ if (basin->down_basin >= 0)
+ down_basin = &(basin_facts[basin->down_basin]);
+ else
+ down_basin = NULL;
+ r = (int)(window->north - basin->northing) / window->ns_res - 1;
+ c = (int)(basin->easting - window->west) / window->ew_res;
+ /*
+ if (r < 0 || r >= nrows || c < 0 || c>= ncols) {
+ G_fatal_error("r:%d c:%d big error", r,c);
+ }
+ */
+ G_get_c_raster_row(fd, buf, r);
+ if (buf[c] < 0) {
+ basin->valid = 0;
+ if (down_basin != NULL)
+ down_basin->valid = 0;
+ }
+ }
+ }
+ G_free(buf);
+
+ return 0;
+}
Property changes on: grass/trunk/raster/r.watershed2/shed/valid.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass/trunk/raster/r.watershed2/shed/watershed.h
===================================================================
--- grass/trunk/raster/r.watershed2/shed/watershed.h (rev 0)
+++ grass/trunk/raster/r.watershed2/shed/watershed.h 2008-11-27 23:00:37 UTC (rev 34568)
@@ -0,0 +1,107 @@
+#ifndef __WATERSHED_H__
+#define __WATERSHED_H__
+
+
+#include <stdio.h>
+#include <grass/gis.h>
+
+#define RAM_NAME "ram"
+#define SEG_NAME "seg"
+#define NON_NAME "watershed"
+#define ACRE_TO_METERSQ 4047.0
+#define MILESQ_TO_ACRE 640.0
+#define HECTACRE_TO_METERSQ 10000.0
+#define KILOSQ_TO_ACRE 247.1
+#define KILOSQ_TO_METERSQ 1000000.0
+#define MILESQ_TO_METERSQ 2590080.0
+#define METERSQ_TO_ACRE 0.00024709661
+#define METERSQ_TO_MILESQ 0.00000038608
+#define METERSQ_TO_HECTACRE 0.0001
+#define METERSQ_TO_KILOSQ 0.000001
+#define METER_TO_FOOT 3.281
+#define INCR 32
+#define INPUT struct inputtttttt
+#define OUTPUT struct outputttttt
+#define BASIN struct mapvallllllll
+#define MAP struct mapssssss
+#define B_FACTS struct bfactsssssss
+#define CAT struct cattttttttttt
+
+INPUT {
+ char *ar_file_name, *com_line_ram, *com_line_seg, *haf_name, *accum_name;
+ char fast, slow;
+};
+
+CAT {
+ int num_cat; /* num of cells */
+ CELL cat_val;
+ CAT *nxt;
+};
+
+BASIN {
+ CAT first_cat; /* linked list of cats with num */
+ double sum_values; /* summation */
+};
+
+MAP {
+ char *name, *mapset;
+ BASIN *basins; /* array of basins */
+ struct Categories cats;
+ char do_cats;
+};
+
+OUTPUT {
+ MAP *maps; /* map layers of output stuff */
+ B_FACTS *basin_facts; /* basin information array */
+ FILE *out_file;
+ struct Cell_head window;
+ int num_maps, num_basins;
+
+ /* output file, display map name, flag for basin by percent, flag for accum of percent */
+ char *file_name, do_basin, do_accum, type_area;
+};
+
+B_FACTS {
+ double str_length, str_slope, accum_length, accum_slope, easting,
+ northing;
+ int num_cells, down_basin, valid;
+};
+
+
+/* accum_down.c */
+int accum_down(OUTPUT *);
+
+/* basin_maps.c */
+int basin_maps(INPUT *, OUTPUT *);
+
+/* com_line.c */
+int com_line_Gwater(INPUT *, OUTPUT *);
+int com_line_add(char **, char *, char *, char *);
+int basin_com_add(char **, double, double, struct Cell_head *);
+int com_add(char **, char *, int);
+
+/* file_in.c */
+int ar_file_in(char *, OUTPUT *);
+
+/* free.c */
+int free_input(INPUT *);
+int free_output(OUTPUT *);
+
+/* insert_cat.c */
+int insert_cat(CAT *, CELL, int);
+CAT *new_cat(CELL, int);
+
+/* intro.c */
+int intro(void);
+
+/* print.c */
+int print_output(OUTPUT *);
+
+/* read.c */
+int read_basins(char *, OUTPUT *);
+
+/* valid.c */
+int valid_basins(char *, OUTPUT *);
+
+
+#endif /* __WATERSHED_H__ */
Property changes on: grass/trunk/raster/r.watershed2/shed/watershed.h
___________________________________________________________________
Name: svn:eol-style
+ native
More information about the grass-commit
mailing list