r.stream.order moved from grass-addons to trunk

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@59460 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.stream.order/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.stream.order
+
+LIBES = $(GISLIB) $(RASTERLIB) $(SEGMENTLIB) $(VECTORLIB) $(DBMILIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP) $(SEGMENTDEP) $(VECTORDEP) $(DBMIDEP)
+
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.stream.order/io.c

@@ -0,0 +1,564 @@
+#include "io.h"
+/* all in ram functions section */
+
+int ram_create_map(MAP * map, RASTER_MAP_TYPE data_type)
+{
+
+    /* 
+     * allocates 0 filled nrows*ncols map of type void;
+     * map parameters are stored in structure;
+     * map: map to be created;
+     * map type to be created must be CELL, FCELL, DCELL;
+     * */
+
+    int r;
+
+    if (data_type < 0 || data_type > 2)
+	G_fatal_error(_("Unable to create raster map of unrecognised type"));
+
+    map->data_type = data_type;
+    map->map_name = NULL;
+    map->nrows = Rast_window_rows();
+    map->ncols = Rast_window_cols();
+    map->data_size = Rast_cell_size(data_type);
+
+    /* preparing internal map */
+    switch (map->data_type) {
+    case CELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(CELL *));
+	break;
+
+    case FCELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(FCELL *));
+	break;
+
+    case DCELL_TYPE:
+	map->map = G_calloc(map->nrows, sizeof(DCELL *));
+	break;
+    }
+
+    for (r = 0; r < map->nrows; ++r)
+	(map->map)[r] = G_calloc(map->ncols, map->data_size);
+
+    return 0;
+}
+
+int ram_read_map(MAP * map, char *input_map_name, int check_res,
+		 RASTER_MAP_TYPE check_data_type)
+{
+    /*
+     * Funciton read external map and put it in MAP structure (created with create_map)
+     * map: map to be read can be of any data type, read map is converted to target map if neccesary.
+     * input_map_name: name of the map to be read;
+     * map pointer to map stucture (created with create_map);
+     * check_res: [1]: check res correspondence between region and map [0 no check];
+     * check_data_type [CELL, FCELL, DCELL] check if reading map is of particular type, [-1] no check;
+     */
+
+    int r, c;
+    char *mapset;
+    struct Cell_head cellhd, this_window;
+    char *maptypes[] = { "CELL", "FCELL", "DCELL" };
+    int input_map_fd;
+    RASTER_MAP_TYPE input_data_type;
+    size_t input_data_size;
+    void *input_buffer = NULL;
+    void *input_pointer;
+
+    /* checking if map exist */
+    mapset = (char *)G_find_raster2(input_map_name, "");
+    if (mapset == NULL)
+	G_fatal_error(_("Raster map <%s> not found"), input_map_name);
+
+    /* checking if region and input are the same */
+    G_get_window(&this_window);
+    Rast_get_cellhd(input_map_name, mapset, &cellhd);
+    if (check_res)
+	if (this_window.ew_res != cellhd.ew_res ||
+	    this_window.ns_res != cellhd.ns_res)
+          G_fatal_error(_("Region resolution and raster map <%s> resolution differs. "
+                          "Run 'g.region rast=%s' to set proper region resolution."),
+                        input_map_name, input_map_name);
+    
+    /* checking if input map is of required type */
+    if (check_data_type != map->data_type)
+	G_debug(1,
+		"ram_open:required map type and internal map type differs: conversion forced!");
+    input_data_type = Rast_map_type(input_map_name, mapset);
+    if (check_data_type != -1)
+	if (input_data_type != check_data_type)
+	    G_fatal_error(_("Raster map <%s> is not of type '%s'"),
+			  input_map_name, maptypes[check_data_type]);
+
+    input_map_fd = Rast_open_old(input_map_name, mapset);
+    input_data_size = Rast_cell_size(input_data_type);
+
+    {				/* reading range */
+	struct Range map_range;
+	struct FPRange map_fp_range;
+	int min, max;
+
+	if (input_data_type == CELL_TYPE) {
+	    Rast_init_range(&map_range);
+	    Rast_read_range(input_map_name, mapset, &map_range);
+	    Rast_get_range_min_max(&map_range, &min, &max);
+	    map->min = (double)min;
+	    map->max = (double)max;
+	}
+	else {
+	    Rast_init_fp_range(&map_fp_range);
+	    Rast_read_fp_range(input_map_name, mapset, &map_fp_range);
+	    Rast_get_fp_range_min_max(&map_fp_range, &(map->min),
+				      &(map->max));
+	}
+    }
+    /* end opening and checking */
+
+    input_buffer = Rast_allocate_buf(input_data_type);
+
+    /* start reading */
+    G_message(_("Reading raster map <%s>..."), input_map_name);
+
+    for (r = 0; r < map->nrows; ++r) {
+	G_percent(r, map->nrows, 2);
+
+	Rast_get_row(input_map_fd, input_buffer, r, input_data_type);
+	input_pointer = input_buffer;
+
+	for (c = 0; c < map->ncols; ++c)
+	    if (!Rast_is_null_value
+		(input_pointer + c * input_data_size, input_data_type))
+		switch (map->data_type) {
+		case CELL_TYPE:
+		    ((CELL **) map->map)[r][c] =
+			Rast_get_c_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case FCELL_TYPE:
+		    ((FCELL **) map->map)[r][c] =
+			Rast_get_f_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case DCELL_TYPE:
+		    ((DCELL **) map->map)[r][c] =
+			Rast_get_d_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		default:
+		    G_fatal_error(_("Wrong internal data type"));
+		    break;
+		}
+    }				/*end for r */
+
+    G_free(input_buffer);
+    G_percent(r, map->nrows, 2);
+    Rast_close(input_map_fd);
+    return 0;
+}				/* end create floating point map */
+
+int ram_reset_map(MAP * map, int value)
+{
+    /*
+     * set all cells in the map to value
+     */
+    int r;
+
+    for (r = 0; r < map->nrows; ++r)
+	memset((map->map)[r], value, map->ncols * map->data_size);
+    return 0;
+}
+
+int ram_write_map(MAP * map, char *output_map_name,
+		  RASTER_MAP_TYPE output_data_type, int convert_to_null,
+		  double value)
+{
+    /* 
+     * write map to disk with output_map_name and output_data_type [CELL, FCELL, DCELL];
+     * if output_data_type = -1 than internal map type is used for output;
+     * if output map != -1 and types differ data_type, conversion is forced
+     * convert to null: check if convert to null a particular value in dataset;
+     */
+
+    int r, c;
+    int output_fd = 0;
+    struct History history;
+    void *row;
+
+    /* check for output format */
+    if (output_data_type == -1)
+	output_data_type = map->data_type;
+
+    if (output_data_type != map->data_type)
+	G_debug(1,
+		"ram_write:required map type and internal map type differs: conversion forced!");
+
+    G_message(_("Writing raster map <%s>..."), output_map_name);
+    output_fd = Rast_open_new(output_map_name, output_data_type);
+
+    /* writing */
+    for (r = 0; r < map->nrows; ++r) {
+	G_percent(r, map->nrows, 2);
+
+	if (convert_to_null) {
+	    row = map->map[r];
+	    switch (map->data_type) {
+	    case CELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((CELL *) row)[c] == (CELL) value)
+			Rast_set_c_null_value(row + c * (map->data_size), 1);
+		break;
+	    case FCELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((FCELL *) row)[c] == (FCELL) value)
+			Rast_set_f_null_value(row + c * (map->data_size), 1);
+		break;
+	    case DCELL_TYPE:
+		for (c = 0; c < map->ncols; ++c)
+		    if (((DCELL *) row)[c] == (DCELL) value)
+			Rast_set_d_null_value(row + c * (map->data_size), 1);
+		break;
+	    default:
+		G_debug(1, "Unable to convert to NULL at: %d %d", r, c);
+	    }
+	}
+
+	Rast_put_row(output_fd, (map->map)[r], output_data_type);
+    }
+    G_percent(r, map->nrows, 2);
+    Rast_close(output_fd);
+    Rast_short_history(output_map_name, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(output_map_name, &history);
+    /* G_message(_("<%s> Done"), output_map_name); */
+    return 0;
+}
+
+int ram_release_map(MAP *map)
+{
+    /* 
+     * free memory allocated for map, set pointer to null;
+     */
+    int r;
+
+    for (r = 0; r < map->nrows; ++r)
+	G_free((map->map)[r]);
+    G_free(map->map);
+    map = NULL;
+    return 0;
+}
+
+
+/* memory swap functions section */
+
+
+int seg_create_map(SEG * seg, int srows, int scols, int number_of_segs,
+		   RASTER_MAP_TYPE data_type)
+{
+    /* create segment  and returns pointer to it;
+     * seg must be declared first;
+     * parameters are stored in structure;
+     * seg: segment to be created;
+     * srows, scols segment size
+     * number of segs max number of segs stored in memory
+     * data_type to be created must be CELL, FCELL, DCELL;
+     */
+
+    char *filename;
+    int fd;
+    int local_number_of_segs;
+
+    seg->fd = -1;
+    seg->filename = NULL;
+    seg->map_name = NULL;
+    seg->mapset = NULL;
+    seg->data_type = data_type;
+    seg->nrows = Rast_window_rows();
+    seg->ncols = Rast_window_cols();
+
+    local_number_of_segs =
+	(seg->nrows / srows + 1) * (seg->ncols / scols + 1);
+    number_of_segs =
+	(number_of_segs >
+	 local_number_of_segs) ? local_number_of_segs : number_of_segs;
+
+    G_debug(3, "seg_creat:number of segments %d", number_of_segs);
+
+    switch (seg->data_type) {
+    case CELL_TYPE:
+	seg->data_size = sizeof(CELL);
+	break;
+    case FCELL_TYPE:
+	seg->data_size = sizeof(FCELL);
+	break;
+    case DCELL_TYPE:
+	seg->data_size = sizeof(DCELL);
+	break;
+    default:
+	G_fatal_error(_("Uunrecognisable data type"));
+    }
+
+    filename = G_tempfile();
+    fd = creat(filename, 0666);
+
+    if (0 >
+	segment_format(fd, seg->nrows, seg->ncols, srows, scols,
+		       seg->data_size)) {
+	close(fd);
+	unlink(filename);
+	G_fatal_error(_("Unable to format segment"));
+    }
+
+    close(fd);
+    if (0 > (fd = open(filename, 2))) {
+	unlink(filename);
+	G_fatal_error(_("Unable re-open file '%s'"), filename);
+    }
+
+    if (0 > (fd = segment_init(&(seg->seg), fd, number_of_segs))) {
+	unlink(filename);
+	G_fatal_error(_("Unable to init segment file or out of memory"));
+    }
+
+    seg->filename = G_store(filename);
+    seg->fd = fd;
+    return 0;
+}
+
+int seg_read_map(SEG * seg, char *input_map_name, int check_res,
+		 RASTER_MAP_TYPE check_data_type)
+{
+
+    /*
+     * Funciton read external map and put it in SEG structure (created with seg_create_map)
+     * map to be read can be of any data type, read map is converted if neccesary.
+     * input_map_name: name of the map to be read;
+     * seg: pointer to map stucture (created with create_map);
+     * check_res: [1]: check res correspondence between region and map [0 no check];
+     * check_data_type [CELL, FCELL, DCELL] check if reading map is of particular type, [-1] no check;
+     */
+
+    int input_fd;
+    int r, c;
+    char *mapset;
+    struct Cell_head cellhd, this_window;
+    char *maptypes[] = { "CELL", "FCELL", "DCELL" };
+    RASTER_MAP_TYPE input_data_type;
+    size_t input_data_size;
+    void *input_buffer = NULL;
+    void *target_buffer = NULL;
+    void *input_pointer = NULL;
+
+    /* checking if map exist */
+    mapset = (char *)G_find_raster2(input_map_name, "");
+    if (mapset == NULL)
+	G_fatal_error(_("Raster map <%s> not found"),
+		      input_map_name);
+    seg->mapset = mapset;
+
+    /* checking if region and input are the same */
+    G_get_window(&this_window);
+    Rast_get_cellhd(input_map_name, mapset, &cellhd);
+
+    /* check resolution equal anyinteger check;  equal 0 no check */
+    if (check_res)
+	if (this_window.ew_res != cellhd.ew_res ||
+	    this_window.ns_res != cellhd.ns_res)
+            G_fatal_error(_("Region resolution and raster map <%s> resolution differs. "
+                            "Run 'g.region rast=%s' to set proper region resolution."),
+                          input_map_name, input_map_name);
+    
+    if (check_data_type != seg->data_type)
+	G_debug(1,
+		"ram_open:required map type and internal map type differs: conversion forced!");
+    input_data_type = Rast_map_type(input_map_name, mapset);
+    if (check_data_type != -1)
+	if (input_data_type != check_data_type)
+	    G_fatal_error(_("Raster map <%s> is not of type '%s'"),
+			  input_map_name, maptypes[check_data_type]);
+
+    input_fd = Rast_open_old(input_map_name, mapset);
+    input_data_size = Rast_cell_size(input_data_type);
+
+    {				/* reading range */
+	struct Range map_range;
+	struct FPRange map_fp_range;
+	int min, max;
+
+	if (input_data_type == CELL_TYPE) {
+	    Rast_init_range(&map_range);
+	    Rast_read_range(input_map_name, mapset, &map_range);
+	    Rast_get_range_min_max(&map_range, &min, &max);
+	    seg->min = (double)min;
+	    seg->max = (double)max;
+	}
+	else {
+	    Rast_init_fp_range(&map_fp_range);
+	    Rast_read_fp_range(input_map_name, mapset, &map_fp_range);
+	    Rast_get_fp_range_min_max(&map_fp_range, &(seg->min),
+				      &(seg->max));
+	}
+    }
+
+    /* end opening and checking */
+
+    G_message(_("Reading raster map <%s>..."), input_map_name);
+    input_buffer = Rast_allocate_buf(input_data_type);
+
+    target_buffer = Rast_allocate_buf(seg->data_type);
+
+    for (r = 0; r < seg->nrows; ++r) {
+	G_percent(r, seg->nrows, 2);
+	Rast_get_row(input_fd, input_buffer, r, input_data_type);
+	input_pointer = input_buffer;
+	memset(target_buffer, 0, seg->ncols * seg->data_size);
+
+	for (c = 0; c < seg->ncols; ++c)
+	    if (!Rast_is_null_value
+		(input_pointer + c * input_data_size, input_data_type)) {
+		switch (seg->data_type) {
+		case CELL_TYPE:
+		    ((CELL *) target_buffer)[c] =
+			Rast_get_c_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case FCELL_TYPE:
+		    ((FCELL *) target_buffer)[c] =
+			Rast_get_f_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		case DCELL_TYPE:
+		    ((DCELL *) target_buffer)[c] =
+			Rast_get_d_value(input_pointer + c * input_data_size,
+					 input_data_type);
+		    break;
+		default:
+		    G_fatal_error(_("Wrong internal data type"));
+		    break;
+		}
+	    }
+
+	if (0 > segment_put_row(&(seg->seg), target_buffer, r)) {
+	    G_free(input_buffer);
+	    G_free(target_buffer);
+	    Rast_close(input_fd);
+	    G_fatal_error(_("Unable to segment put row %d for raster map <%s>"),
+			  r, input_map_name);
+	}
+    }				/* end for row */
+
+    G_percent(r, seg->nrows, 2);
+    Rast_close(input_fd);
+    G_free(input_buffer);
+    G_free(target_buffer);
+
+    seg->map_name = G_store(input_map_name);
+    seg->mapset = G_store(mapset);
+
+    return 0;
+}
+
+int seg_reset_map(SEG * seg, int value)
+{
+    /*
+     * set all cells in the map to value
+     */
+    int r, c;
+
+    for (r = 0; r < seg->nrows; ++r)
+	for (c = 0; c < seg->ncols; ++c)
+	    segment_put(&(seg->seg), &value, r, c);
+
+    return 0;
+}
+
+int seg_write_map(SEG * seg, char *output_map_name,
+		  RASTER_MAP_TYPE output_data_type, int convert_to_null,
+		  double value)
+{
+    /* 
+     * write seg to disk with output_map_name and output_data_type [CELL, FCELL, DCELL];
+     * if output_data_type = -1 than internal map type is used for output;
+     * if output map != -1 and types differ data_type, conversion is forced
+     * convert to null: check if convert to null a particular value in dataset;
+     */
+    int output_fd;
+    int r, c;
+    void *output_buffer;
+    void *row;
+    struct History history;
+
+    /* check for output format */
+    if (output_data_type == -1)
+	output_data_type = seg->data_type;
+
+    if (output_data_type != seg->data_type)
+	G_debug(1,
+		"ram_write:required map type and internal map type differs: conversion forced!");
+
+    G_message(_("Writing raster map <%s>..."), output_map_name);
+    output_fd = Rast_open_new(output_map_name, output_data_type);
+    output_buffer = Rast_allocate_buf(output_data_type);
+    segment_flush(&(seg->seg));
+
+    /* writing */
+    for (r = 0; r < seg->nrows; ++r) {
+
+	G_percent(r, seg->nrows, 2);
+	if (0 > segment_get_row(&(seg->seg), output_buffer, r))
+	    G_warning(_("Unable to segment read row %d for raster map <%s>"),
+		      r, output_map_name);
+
+	if (convert_to_null) {
+
+	    row = output_buffer;
+	    switch (seg->data_type) {
+	    case CELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((CELL *) output_buffer)[c] == (CELL) value)
+			Rast_set_c_null_value(row + c * (seg->data_size), 1);
+		break;
+	    case FCELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((FCELL *) output_buffer)[c] == (FCELL) value)
+			Rast_set_f_null_value(row + c * (seg->data_size), 1);
+		break;
+	    case DCELL_TYPE:
+		for (c = 0; c < seg->ncols; ++c)
+		    if (((DCELL *) output_buffer)[c] == (DCELL) value)
+			Rast_set_d_null_value(row + c * (seg->data_size), 1);
+		break;
+	    default:
+		G_warning(_("Unable to convert to NULL at: %d %d"), r,
+			  c);
+	    }
+	}
+	Rast_put_row(output_fd, output_buffer, output_data_type);
+    }
+
+    G_percent(r, seg->nrows, 2);
+    G_free(output_buffer);
+    Rast_close(output_fd);
+    Rast_short_history(output_map_name, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(output_map_name, &history);
+    /* G_message(_("%s Done"), output_map_name); */
+
+    return 0;
+}
+
+int seg_release_map(SEG * seg)
+{
+    /* 
+     * release segment close files, set pointers to null;
+     */
+    segment_release(&(seg->seg));
+    close(seg->fd);
+    unlink(seg->filename);
+
+    if (seg->map_name)
+	G_free(seg->map_name);
+    if (seg->mapset)
+	G_free(seg->mapset);
+
+    return 0;
+}

raster/r.stream.order/io.h

@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+#define NOT_IN_REGION(x) (r + nextr[(x)] < 0 || r + nextr[(x)] > (nrows - 1) || \
+                          c + nextc[(x)] < 0 || c + nextc[(x)] > (ncols - 1))
+#define NR(x) (r + nextr[(x)])
+#define NC(x) (c + nextc[(x)])
+#define INDEX(r,c) ((r) * ncols + (c))
+#define DIAG(x) (((x) + 4) > 8 ? ((x) - 4) : ((x) + 4))
+
+#define SROWS 256
+#define SCOLS 256
+
+typedef struct {
+	void **map; /* matrix of data */
+	double min, max; /* data range : may requre casting */
+	int nrows, ncols;
+	char *map_name; /* map name, unused */
+	RASTER_MAP_TYPE data_type; /* type of data */
+	size_t data_size; /* type of data */
+} MAP;
+
+typedef struct {
+	SEGMENT seg;		/* segmented data store */
+	int fd;					/* segment temporary file name descriptor */
+	char *filename; /* segment temporary file name */
+	char *map_name; /* map name converted to segment */
+	char *mapset;
+	int nrows, ncols; /* store nrows and rcols */
+	RASTER_MAP_TYPE data_type; /* data type of the map */
+	size_t data_size; /* size of cell returned by sizeof */
+	double min, max; /* data range */
+} SEG;
+
+
+/* all in ram functions */
+int ram_create_map(MAP *, RASTER_MAP_TYPE);
+int ram_read_map(MAP *, char *, int, RASTER_MAP_TYPE);
+int ram_reset_map(MAP *, int);
+int ram_write_map(MAP *, char *, RASTER_MAP_TYPE, int, double);
+int ram_release_map(MAP *);
+int ram_destory_map(MAP *);
+
+/* memory swap functions */
+int seg_create_map(SEG *, int, int, int, RASTER_MAP_TYPE);
+int seg_read_map(SEG *, char *, int, RASTER_MAP_TYPE);
+int seg_reset_map (SEG *, int);
+int seg_write_map(SEG *, char *, RASTER_MAP_TYPE, int, double);
+int seg_release_map(SEG *);

raster/r.stream.order/local_proto.h

@@ -0,0 +1,27 @@
+#include "io.h"
+#include "local_vars.h"
+
+/* stream init */
+int stream_init(int min_index_of_stream, int max_index_of_stream);
+int stream_sample_map(char* input_map_name, int number_of_streams, int what);
+
+/* stream topology */
+int ram_stream_topology(CELL** streams, CELL** dirs, int number_of_streams);
+int ram_stream_geometry(CELL** streams, CELL** dirs);
+int seg_stream_topology(SEGMENT* streams, SEGMENT* dirs, int number_of_streams);
+int seg_stream_geometry(SEGMENT* streams, SEGMENT* dirs);
+
+/* stream order */
+int strahler(int* strahler);
+int shreve(int* shreve);
+int horton(const int* strahler, int* horton, int number_of_streams);
+int hack(int* hack, int* topo_dim, int number_of_streams);
+
+/* stream raster close */
+int ram_close_raster_order(CELL** streams, int number_of_streams, int zerofill);
+int seg_close_raster_order(SEGMENT* streams, int number_of_streams, int zerofill);
+
+/* stream vector */
+int ram_create_vector(CELL** streams, CELL** dirs, char* out_vector, int number_of_streams);
+int seg_create_vector(SEGMENT* streams, SEGMENT* dirs, char* out_vector, int number_of_streams);
+int stream_add_table (int number_of_streams);

raster/r.stream.order/local_vars.h

@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+
+#ifdef MAIN
+#  define GLOBAL
+#else
+#  define GLOBAL extern
+#endif
+
+typedef enum {o_streams, o_dirs, o_elev, o_accum, input_size } inputs;
+/* eny new oredering systems must be also declared here before orders_size*/
+typedef enum {o_strahler, o_horton, o_shreve, o_hack, o_topo, orders_size } orders;
+
+typedef struct {
+	char* name;
+	int required;
+	char* description;
+} IO;
+
+typedef struct {
+    int stream;	/* topology */
+    int next_stream;
+    int trib_num;
+    int trib[5];
+    int cells_num;
+    int init; /* index to recalculate into r,c*/
+    int outlet; /* index to recalculate into r,c */
+    double length;   /* geometry */ 
+    double accum_length;
+    double distance; /* distance to outlet */
+    double stright; /* use next to calculate sinusoid */
+    double accum;
+    double init_elev;
+    double outlet_elev; /* use next to calculate drop and gradient */
+} STREAM;
+
+GLOBAL char** output_map_names;
+GLOBAL int ** all_orders;
+
+GLOBAL int nextr[9];
+GLOBAL int nextc[9];
+GLOBAL int nrows, ncols;
+GLOBAL int use_vector, use_accum;
+
+
+/* stream topo */
+GLOBAL int init_num, outlet_num;
+GLOBAL STREAM* stream_attributes;
+GLOBAL unsigned int* init_streams;
+GLOBAL unsigned int* outlet_streams;
+GLOBAL unsigned long int* init_cells;
+
+/* output vector */
+GLOBAL struct Map_info Out; 

raster/r.stream.order/main.c

@@ -0,0 +1,286 @@
+
+/****************************************************************************
+ *
+ * MODULE:		r.stream.order
+ * AUTHOR(S):		Jarek Jasiewicz jarekj amu.edu.pl
+ *							 
+ * PURPOSE:		Calculate Strahler's and more streams hierarchy
+ *			It use r.stream.extract or r.watershed output files: 
+ * 			stream, direction, accumulation and elevation. 
+ * 			The output are set of raster maps and vector file containing
+ * 			addational stream attributes.
+ *
+ * COPYRIGHT:		(C) 2009-2014 by the GRASS Development Team
+ *
+ *			This program is free software under the GNU General Public 
+ *			License (>=v2). Read the file COPYING that comes with GRASS
+ *			for details.
+ *
+ *****************************************************************************/
+#define MAIN
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int nextr[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+int nextc[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+
+
+int main(int argc, char *argv[])
+{
+
+    IO input[] = {
+        {"stream_rast", YES, _("Name of input streams mask raster map")},
+	{"direction", YES, _("Name of input flow direction raster map")},
+	{"elevation", NO, _("Name of input elevation raster map")},
+	{"accumulation", NO, _("Name of input accumulation raster map")}
+    };
+
+    /* add new basic rdering here and in local_vars.h declaration 
+     * it will be added to output without addational programing.
+     * Ordering functions are to be added in stream_order.c file.
+     * Addational operations may shall be done in common part sections
+     * derivative orders (like Scheideggers/Shreve) shall be added only 
+     * to table definition as a formula and to description file. */
+    IO output[] = {
+	{"strahler", NO, _("Name for output Strahler's stream order raster map")},
+	{"horton", NO, _("Name for output original Hortons's stream order raster map")},
+	{"shreve", NO, _("Name for output Shereve's stream magnitude raster map")},
+	{"hack", NO, _("Name for output Hack's streams or Gravelius stream hierarchy raster map")},
+	{"topo", NO, _("Name for output topological dimension of streams raster map")}
+    };
+    struct GModule *module;	/* GRASS module for parsing arguments */
+
+    struct Option *opt_input[input_size];
+    struct Option *opt_output[orders_size];
+    struct Option *opt_swapsize;
+    struct Option *opt_vector;
+    struct Flag *flag_zerofill, *flag_accum, *flag_segmentation;
+
+    int output_num = 0;
+    /* int num_seg; */		/* number of segments */ 
+    int segmentation, zerofill;
+    int i;			/* iteration vars */
+    int number_of_segs;
+    int number_of_streams;
+    char *in_streams = NULL, *in_dirs = NULL, *in_elev = NULL, *in_accum =
+	NULL;
+    char *out_vector = NULL;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->label = _("Calculates Strahler's and more streams hierarchy.");
+    module->description =
+	_("Basic module for topological analysis of drainage network.");
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("hydrology"));
+    G_add_keyword(_("stream network"));
+    G_add_keyword(_("stream order"));
+
+    for (i = 0; i < input_size; ++i) {
+	opt_input[i] = G_define_standard_option(G_OPT_R_INPUT);
+	opt_input[i]->key = input[i].name;
+	opt_input[i]->required = input[i].required;
+	opt_input[i]->description = _(input[i].description);
+        opt_input[i]->guisection = _("Input maps");
+    }
+
+    opt_vector = G_define_standard_option(G_OPT_V_OUTPUT);
+    opt_vector->key = "stream_vect";
+    opt_vector->required = NO;
+    opt_vector->description =
+	_("Name for output vector map to write stream attributes");
+    opt_vector->guisection = _("Output maps");
+
+    for (i = 0; i < orders_size; ++i) {
+	opt_output[i] = G_define_standard_option(G_OPT_R_OUTPUT);
+	opt_output[i]->key = output[i].name;
+	opt_output[i]->required = NO;
+	opt_output[i]->description = _(output[i].description);
+	opt_output[i]->guisection = _("Output maps");
+    }
+
+    opt_swapsize = G_define_option();
+    opt_swapsize->key = "memory";
+    opt_swapsize->type = TYPE_INTEGER;
+    opt_swapsize->answer = "300";
+    opt_swapsize->description = _("Max memory used in memory swap mode (MB)");
+    opt_swapsize->guisection = _("Memory settings");
+
+    flag_zerofill = G_define_flag();
+    flag_zerofill->key = 'z';
+    flag_zerofill->description =
+	_("Create zero-valued background instead of NULL");
+
+    flag_segmentation = G_define_flag();
+    flag_segmentation->key = 'm';
+    flag_segmentation->description = _("Use memory swap (operation is slow)");
+    flag_segmentation->guisection = _("Memory settings");
+
+    flag_accum = G_define_flag();
+    flag_accum->key = 'a';
+    flag_accum->description =
+	_("Use flow accumulation to trace horton and hack models");
+
+    if (G_parser(argc, argv))	/* parser */
+	exit(EXIT_FAILURE);
+
+
+    /* check output names */
+    zerofill = (flag_zerofill->answer != 0);
+    segmentation = (flag_segmentation->answer != 0);
+    use_accum = (flag_accum->answer != 0);
+    use_vector = (opt_vector->answer != NULL);
+
+    number_of_segs = (int)atof(opt_swapsize->answer);
+    number_of_segs = number_of_segs < 32 ? (int)(32 / 0.12) : number_of_segs / 0.12;
+
+    if (use_vector)
+	if (!opt_input[o_elev]->answer || !opt_input[o_accum]->answer)
+	    G_fatal_error(_("To calculate vector map both accumulation and elevation raster maps are required"));
+    if (use_accum)
+	if (!opt_input[o_accum]->answer)
+	    G_fatal_error(_("Flag -a (use accumulation) accumulation raster map is required"));
+
+    for (i = 0; i < orders_size; ++i) {
+	if (!opt_output[i]->answer)
+	    continue;
+	output_num++;
+    }				/* end for */
+
+    if (!output_num && !opt_vector->answer)
+	G_fatal_error(_("You must select one or more output orders raster maps or insert the table name"));
+
+    /* start */
+    in_streams = opt_input[o_streams]->answer;
+    in_dirs = opt_input[o_dirs]->answer;
+    in_elev = opt_input[o_elev]->answer;
+    in_accum = opt_input[o_accum]->answer;
+    out_vector = opt_vector->answer;
+
+    output_map_names = (char **)G_malloc(orders_size * sizeof(char *));
+    for (i = 0; i < orders_size; ++i)
+	output_map_names[i] = opt_output[i]->answer;
+
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    /* ALL IN RAM VERSION */
+    if (!segmentation) {
+        G_message(_("All in RAM calculation..."));
+	MAP map_streams, map_dirs;
+	CELL **streams, **dirs;
+
+	ram_create_map(&map_streams, CELL_TYPE);
+	ram_read_map(&map_streams, in_streams, 1, CELL_TYPE);
+	ram_create_map(&map_dirs, CELL_TYPE);
+	ram_read_map(&map_dirs, in_dirs, 1, CELL_TYPE);
+	stream_init((int)map_streams.min, (int)map_streams.max);
+	number_of_streams = (int)(map_streams.max + 1);
+	streams = (CELL **) map_streams.map;
+	dirs = (CELL **) map_dirs.map;
+
+	ram_stream_topology(streams, dirs, number_of_streams);
+
+	if (out_vector || output_map_names[o_horton] ||
+	    output_map_names[o_hack] || output_map_names[o_topo])
+	    ram_stream_geometry(streams, dirs);
+
+	/* common part */
+	if (use_vector) {
+	    stream_sample_map(in_elev, number_of_streams, 0);
+	    stream_sample_map(in_elev, number_of_streams, 1);
+	}
+	if (use_accum || use_vector)
+	    stream_sample_map(in_accum, number_of_streams, 2);
+
+	if (output_map_names[o_strahler] || output_map_names[o_horton] ||
+	    out_vector)
+	    strahler(all_orders[o_strahler]);
+
+	if (output_map_names[o_horton] || out_vector)
+	    horton(all_orders[o_strahler], all_orders[o_horton],
+		   number_of_streams);
+
+	if (output_map_names[o_shreve] || out_vector)
+	    shreve(all_orders[o_shreve]);
+
+	if (output_map_names[o_hack] || output_map_names[o_topo] ||
+	    out_vector)
+	    hack(all_orders[o_hack], all_orders[o_topo], number_of_streams);
+	/* end of common part */
+
+	if (out_vector)
+	    ram_create_vector(streams, dirs, out_vector, number_of_streams);
+
+	ram_close_raster_order(streams, number_of_streams, zerofill);
+	ram_release_map(&map_streams);
+	ram_release_map(&map_dirs);
+	/* end ram section */
+    }
+
+    /* SEGMENTATION VERSION */
+    if (segmentation) {
+        G_message(_("Memory swap calculation (may take some time)..."));
+
+	SEG map_streams, map_dirs;
+	SEGMENT *streams, *dirs;
+
+	seg_create_map(&map_streams, SROWS, SCOLS, number_of_segs, CELL_TYPE);
+	seg_read_map(&map_streams, in_streams, 1, CELL_TYPE);
+	seg_create_map(&map_dirs, SROWS, SCOLS, number_of_segs, CELL_TYPE);
+	seg_read_map(&map_dirs, in_dirs, 1, CELL_TYPE);
+	stream_init((int)map_streams.min, (int)map_streams.max);
+	number_of_streams = (int)(map_streams.max + 1);
+	streams = &map_streams.seg;
+	dirs = &map_dirs.seg;
+
+	seg_stream_topology(streams, dirs, number_of_streams);
+
+	if (out_vector || output_map_names[o_horton] ||
+	    output_map_names[o_hack] || output_map_names[o_topo])
+	    seg_stream_geometry(streams, dirs);
+
+	/* common part */
+	if (use_vector) {
+	    stream_sample_map(in_elev, number_of_streams, 0);
+	    stream_sample_map(in_elev, number_of_streams, 1);
+	}
+	if (use_accum || use_vector)
+	    stream_sample_map(in_accum, number_of_streams, 2);
+
+	if (output_map_names[o_strahler] || output_map_names[o_horton] ||
+	    out_vector)
+	    strahler(all_orders[o_strahler]);
+
+	if (output_map_names[o_horton] || out_vector)
+	    horton(all_orders[o_strahler], all_orders[o_horton],
+		   number_of_streams);
+
+	if (output_map_names[o_shreve] || out_vector)
+	    shreve(all_orders[o_shreve]);
+
+	if (output_map_names[o_hack] || output_map_names[o_topo] ||
+	    out_vector)
+	    hack(all_orders[o_hack], all_orders[o_topo], number_of_streams);
+	/* end of common part */
+
+	if (out_vector)
+	    seg_create_vector(streams, dirs, out_vector, number_of_streams);
+
+	seg_close_raster_order(streams, number_of_streams, zerofill);
+	seg_release_map(&map_streams);
+	seg_release_map(&map_dirs);
+	/* end segmentation section */
+    }
+
+
+    /* free */
+    G_free(stream_attributes);
+    G_free(init_streams);
+    G_free(outlet_streams);
+    G_free(init_cells);
+
+    exit(EXIT_SUCCESS);
+
+}

raster/r.stream.order/r.stream.order.html

@@ -0,0 +1,363 @@
+<h2>DESCRIPTION</h2>
+
+The <em>r.stream.order</em> calculates Strahler's and other stream
+hierarchy methods. It is a basic module for topological analysis of
+drainage networks.
+
+<h2>OPTIONS</h2>
+<dl>
+<dt><b>-z</b></dt>
+<dd>Creates zero-value background instead of NULL. For some reason
+(like map algebra calculation) zero-valued background may be
+required. This flag produces zero-filled background instead of null
+(default).</dd>
+<dt><b>-a</b></dt>
+<dd>Uses accumulation map instead of cumulated stream length to
+determine main branch at bifurcation. Works well only with SFD
+networks</dd>
+<dt><b>-m</b></dt>
+<dd>Only for very large data sets. Use segment library to optimise
+memory consumption during analysis</dd>
+<dt><b>stream_rast</b></dt>
+<dd>Name of input stream map on which ordering will be performed
+produced by <em><a href="r.watershed.html">r.watershed</a></em> or
+<em><a href="r.stream.extract.html">r.stream.extract</a></em>. Because
+streams network produced
+by <em><a href="r.watershed.html">r.watershed</a></em> and
+<em><a href="r.stream.extract.html">r.stream.extract</a></em> may
+slightly differ in detail it is required to use both stream and
+direction map produced by the same module. Stream background shall
+have NULL value or zero value. Background values of NULL are by
+default produced by <em><a href="r.watershed.html">r.watershed</a></em> or
+<em><a href="r.stream.extract.html">r.stream.extract</a></em>. If not
+0 or NULL use <em><a href="r.mapcalc.html">r.mapcalc</a></em> to set background
+values to null.
+</dd>
+<dt><b>direction</b></dt>
+<dd>Name of input direction map produced by <em><a href="r.watershed.html">r.watershed</a></em> or
+<em><a href="r.stream.extract.html">r.stream.extract</a></em>. If
+<em><a href="r.stream.extract.html">r.stream.extract</a></em> output
+map is used, it only has non-NULL values in places where streams
+occur. NULL (nodata) cells are ignored, zero and negative values are
+valid direction data if they vary from -8 to 8 (CCW from East in steps
+of 45 degrees). Direction map shall be of type CELL values. Region
+resolution and map resolution must be the same.
+Also <b>stream_rast</b> network and <b>direction</b> maps must have
+the same resolution. It is checked by default. If resolutions differ
+the module informs about it and stops. Region boundary and maps
+boundary may be differ but it may lead to unexpected results.</dd>
+<dt><b>accumulation</b></dt>
+<dd>Flow accumulation (optional, not recommended): name of flow
+accumulation file produced by <em><a href="r.watershed.html">r.watershed</a></em> or
+<em><a href="r.stream.extract.html">r.stream.extract</a></em>. This
+map is an option only if Horton's or Hack's ordering is
+performed. Normally both Horton and Hack ordering is calculated on
+cumulative stream length which is calculated internally. Flow
+accumulation can be used if user want to calculate main stream as most
+accumulated stream. Flow accumulation map shall be of DCELL type, as
+is by default produced by r.watershed or converted do DCELL with
+<em><a href="r.mapcalc.html">r.mapcalc</a></em>.</dd>
+<dt><b>elevation</b></dt>
+<dd>Used to calculate geometrical properties of the network stored in the
+table.</dd>
+</dl>
+
+<h3>OUTPUTS</h3>
+
+<p>At least one output map is required:
+
+<dl>
+<dt><b>stream_vect</b></dt>
+<dd>Vector network with table where stream network topology can be
+stored. Because <em><a href="r.stream.order.html">r.stream.order</a></em>
+is prepared to work both with <em><a href="r.watershed.html">r.watershed</a></em> or
+<em><a href="r.stream.extract.html">r.stream.extract</a></em>, it may
+be used to create correct vector
+from <em><a href="r.watershed.html">r.watershed</a></em> results.<dd>
+
+<dt><b>strahler</b></dt>
+<dd>Name of Strahler's stream order output map: see notes for
+detail. </dd>
+
+<dt><b>shreve</b></dt>
+<dd>Name of Shreve's stream magnitude output map: see notes for detail.</dd>
+
+<dt><b>horton</b></dt>
+<dd>Name of Horton's stream order output map (require accum file): see
+notes for detail.</dd>
+
+<dt><b>hack</b></dt>
+<dd>Name of Hack's main streams output map : see notes for
+detail.</dd>
+
+<dt><b>topo</b></dt>
+<dd>Name of topological dimensions streams output map: see notes for
+detail.</dd>
+</dl>
+
+<h2>NOTES</h2>
+
+Module can work only if direction map, stream map and region map has
+same settings. It is also required that stream map and direction map
+come from the same source. For lots of reason this limitation probably
+cannot be omitted. This means if stream map comes from
+<em><a href="r.stream.extract.html">r.stream.extract</a></em> also
+direction map
+from <em><a href="r.stream.extract.html">r.stream.extract</a></em>
+must be used. If stream network was generated with MFD method also MFD
+direction map must be used. Nowadays f direction map comes from
+<em><a href="r.stream.extract.html">r.stream.extract</a></em> must be
+patched by direction map
+from <em><a href="r.watershed.html">r.watershed</a></em>. (with <em><a href="r.patch.html">r.patch</a></em>).
+
+<h3>Stream ordering example</h3>
+<center>
+<img src="orders.png" border="1"><br>
+</center>
+
+<h4>Strahler's stream order</h4>
+
+Strahler's stream order is a modification of Horton's streams order
+which fixes the ambiguity of Horton's ordering. In Strahler's
+ordering the main channel is not determined; instead the ordering is
+based on the hierarchy of tributaries. The ordering follows these
+rules:
+
+<ol>
+<li>if the node has no children, its Strahler order is 1.
+<li>if the node has one and only one tributary with Strahler greatest
+order i, and all other tributaries have order less than i, then the
+order remains i.
+<li>if the node has two or more tributaries with greatest order i,
+then the Strahler order of the node is i + 1.
+</ol>
+Strahler's stream ordering starts in initial links which assigns order
+one. It proceeds downstream. At every node it verifies that there are
+at least 2 equal tributaries with maximum order. If not it continues
+with highest order, if yes it increases the node's order by 1 and
+continues downstream with new order.
+
+<h4>Advantages and disadvantages of Strahler's ordering</h4>
+
+Strahler's stream order has a good mathematical background. All
+catchments with streams in this context are directed graphs, oriented
+from the root towards the leaves. Equivalent definition of the
+Strahler number of a tree is that it is the height of the largest
+complete binary tree that can be homeomorphically embedded into the
+given tree; the Strahler number of a node in a tree is equivalent to
+the height of the largest complete binary tree that can be embedded
+below that node. The disadvantage of that methods is the lack of
+distinguishing a main channel which may interfere with the analytical
+process in highly elongated catchments
+
+<h4>Horton's stream ordering</h4>
+
+Horton's stream order applies to the stream as a whole but not to
+segments or links since the order on any channel remains unchanged
+from source till it "dies" in the higher order stream or in the outlet
+of the catchment. The main segment of the catchment gets the order of
+the whole catchment, while its tributaries get the order of their own
+subcatchments. The main difficulties of the Horton's order are
+criteria to be considered to distinguish between "true" first order
+segments and extension of higher order segments. That is the reason
+why Horton's ordering has rather historical sense and is substituted
+by the more unequivocal Strahler's ordering system. There are no
+natural algorithms to order stream network according to Horton'
+paradigm. The algorithm used in r.stream.order requires to first
+calculate Strahler's stream order (downstream) and next recalculate to
+Horton ordering (upstream). To make a decision about proper ordering
+it uses first Strahler ordering, and next, if both branches have the
+same orders it uses flow accumulation to choose the actual link. The
+algorithm starts with the outlet, where the outlet link is assigned
+the corresponding Strahler order. Next it goes upstream and determines
+links according to Strahler ordering. If the orders of tributaries
+differ, the algorithm proceeds with the channel of highest order, if
+all orders are the same, it chooses that one with higher flow length
+rate or higher catchment area if accumulation is used. When it reaches
+the initial channel it goes back to the last undetermined branch,
+assign its Strahler order as Horton order and goes upstream to the
+next initial links. In that way stream orders remain unchanged from
+the point where Horton's order have been determined to the source.
+  
+<h4>Advantages and disadvantages of Horton's ordering</h4> 
+
+The main advantages of Horton's ordering is that it produces natural
+stream ordering with main streams and its tributaries. The main
+disadvantage is that it requires prior Strahler's ordering. In some
+cases this may result in unnatural ordering, where the highest order
+will be ascribed not to the channel with higher accumulation but to
+the channel which leads to the most branched parts of the catchment.
+
+<h4>Shreve's stream magnitude</h4>
+
+That ordering method is similar to Consisted Associated Integers
+proposed by Scheidegger. It assigns magnitude of 1 for every initial
+channel. The magnitude of the following channel is the sum of
+magnitudes of its tributaries. The number of a particular link is the
+number of initials which contribute to it.
+
+<h4>Scheidegger's stream magnitude</h4>
+
+That ordering method is similar to Shreve's stream magnitude. It
+assigns magnitude of 2 for every initial channel. The magnitude of the
+following channel is the sum of magnitudes of its tributaries. The
+number of a particular link is the number of streams -1 contributing
+to it. Consisted Associated Integers (Scheidegger) is available only
+in attribute table. To achieve Consisted Associated Integers
+(Scheidegger) raster the result of Shreve's magnitude is to be
+multiplied by 2:
+
+<div class="code"><pre>
+r.mapcalc "scheidegger = shreve * 2.0"
+</pre></div>
+
+<h4>Drwal's stream hierarchy (old style)</h4>
+
+That ordering method is a compromise between Strahler ordering and
+Shreve magnitude. It assigns order of 1 for every initial channel. The
+order of the following channel is calculated according Strahler
+formula, except, that streams which do not increase order of next
+channel are not lost. To increase next channel to the higher order R+1
+are require two channels of order R, or one R and two R-1 or one R,
+and four R-2 or one R, one R-1 and two R-2 etc. The order of
+particular link show the possible value of Strahler'order if the
+network was close to idealised binary tree. Drwal's order is
+aviallable only in attribute table.To achieve Drwal's raster the
+result of Shreve's magnitude is to be recalculated according
+formula: <tt>floor(log(shreve,2))+1</tt>
+
+<div class="code"><pre>
+r.mapcalc "drwal = int(log(shreve,2.0)) + 1.0"
+</pre></div>
+
+<h4>Advantages and disadvantages of Drwal's hierarhy</h4> 
+
+The main advantages of Drwal's hierarchy is that it produces natural
+stream ordering with which takes into advantage both ordering and
+magnitude. It shows the real impact of particular links of the network
+run-off. The main disadvantage is that it minimise bifuraction ratio
+of the network.
+
+<h4>Hack's main streams or Gravelius order</h4>
+
+This method of ordering calculates main streams of main catchment and
+every subcatchments. Main stream of every catchment is set to 1, and
+consequently all its tributaries receive order 2. Their tributaries
+receive order 3 etc. The order of every stream remains constant up to
+its initial link. The route of every main stream is determined
+according to the maximum flow length value of particular streams. So
+the main stream of every subcatchment is the longest stream or stream
+with highest accumulation rate if accumulation map is used. In most
+cases the main stream is the longest watercourse of the catchment, but
+in some cases, when a catchment consists of both rounded and elongated
+subcatchments these rules may not be maintained. The algorithm assigns
+1 to every outlets stream and goes upstream according to maximum flow
+accumulation of every branch. When it reaches an initial stream it
+step back to the first unassigned confluence. It assigns order 2 to
+unordered tributaries and again goes upstream to the next initial
+stream. The process runs until all branches of all outlets are
+ordered.
+
+<h4>Advantages and disadvantages of main stream ordering</h4>
+
+The biggest advantage of that method is the possibility to compare and
+analyze topology upstream, according to main streams. Because all
+tributaries of main channel have order of 2, streams can be quickly
+and easily filtered and its proprieties and relation to main stream
+determined. The main disadvantage of that method is the problem with
+the comparison of subcatchment topology of the same
+order. Subcatchments of the same order may be both highly branched and
+widespread in the catchment area and a small subcatchment with only
+one stream.
+
+<h4>Topological dimension streams order</h4>
+
+This method of ordering calculates topological distance of every
+stream from catchment outlet.
+
+<h4>Stream network topology table description connected with vector file</h4>
+
+<ul>
+	<li><b>cat</b> integer: category;
+	<li><b>stream</b>integer: stream number, usually equal to cat;
+	<li><b>next_stream</b> integer: stream to which contribute current
+stream (downstream);
+	<li><b>prev_streams</b>; two or more contributing streams (upstream);
+	<li><b>strahler</b> integer: Strahler's stream order:
+	<li><b>horton</b> integer: Hortons's stream order:
+	<li><b>shreve</b> integer: Shreve's stream magnitude;
+	<li><b>scheidegger</b> integer: Scheidegger's Consisted Associated
+Integers;
+	<li><b>drwal</b> integer: Drwal's stream hierarchy;
+	<li><b>hack</b> integer: Hack's main streams or Gravelius order;
+	<li><b>topo</b> integer: Topological dimension streams order;
+	<li><b>length</b> double precision: stream length;
+	<li><b>cum_length</b> double precision: length of stream from source;
+	<li><b>out_dist</b> double precision: distance of current stream init
+from outlet;
+	<li><b>stright</b> double precision: length of stream as stright line;
+	<li><b>sinusiod</b> double precision: fractal dimension: stream
+length/stright stream length;
+	<li><b>elev_init</b> double precision: elevation of  stream init;
+	<li><b>elev_outlet</b> double precision: elevation of  stream outlet;
+	<li><b>drop</b> double precision: difference between stream init and
+outlet + drop outlet;
+	<li><b>out_drop</b> double precision: drop at the outlet of the stream;
+	<li><b>gradient</b> double precision: drop/length;
+</ul>
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+g.region -p -a rast=elevation
+r.watershed elevation=elevation threshold=10000 drainage=direction stream=streams
+r.stream.order stream_rast=streams direction=direction strahler=riverorder_strahler \
+  horton=riverorder_horton shreve=riverorder_shreve hack=riverorder_hack \
+  topo=river_topodim
+
+# vector river network
+r.watershed elevation=elevation threshold=10000 accumulation=accum
+r.stream.order stream_rast=streams direction=direction elevation=elevation \
+  accumulation=accum stream_vect=river_vector
+</pre></div>
+
+<h2>REFERENCES</h2>
+<ul>
+<li>Drwal, J., (1982), <i>Wyksztalecenie i organizacja sieci hydrograficznej jako
+podstawa oceny struktury odplywu na terenach m;odoglacjalnych</i>, <b>Rozprawy i
+monografie</b>, Gdansk 1982, 130 pp (in Polish)
+<li>Hack, J., (1957), <i>Studies of longitudinal stream profiles in Virginia and
+Maryland</i>, 
+<b>U.S. Geological Survey Professional Paper</b>, 294-B
+<li>Horton, R. E. (1945), <i>Erosional development of streams and their drainage
+basins: hydro-physical approach to quantitative morphology</i>,<b>Geological
+Society of America Bulletin</b> 56 (3): 275-370<br>
+Scheidegger A. E., (1966), <i>Statistical Description of River Networks</i>.
+<b>Water Resour. Res.</b>, 2(4): 785-790
+<li>Shreve, R.,  (1966),<i>Statistical Law of Stream Numbers</i>, <b>J. Geol.</b>,
+74, 17-37.
+<li>Strahler, A. N. (1952), <i>Hypsometric (area-altitude) analysis of erosional
+topology</i>,<b>Geological Society of America Bulletin</b> 63 (11): 1117-1142
+<li>Strahler, A. N. (1957), <i>Quantitative analysis of watershed
+geomorphology</i>,<b>Transactions of the American Geophysical Union</b> 8 (6):
+913-920.
+<li>Woldenberg, M. J., (1967), <i>Geography and properties of surfaces,</i>
+<b>Harvard Papers in Theoretical Geography</b>, 1: 95-189.
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.stream.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
+<p>
+<i>Last changed: $Date$</i>
+

raster/r.stream.order/stream_init.c

@@ -0,0 +1,134 @@
+#include "local_proto.h"
+int stream_init(int min_index_of_stream, int max_index_of_stream)
+{
+    int i, j;
+    int number_of_streams = max_index_of_stream;
+
+    if (max_index_of_stream == 0)
+        G_fatal_error(_("Empty stream input raster map"));
+    if (min_index_of_stream < 0)
+	G_fatal_error(_("Stream map has negative values"));
+
+    stream_attributes =
+	(STREAM *) G_malloc((number_of_streams + 1) * sizeof(STREAM));
+    for (i = 0; i < max_index_of_stream + 1; ++i) {
+	stream_attributes[i].next_stream = -1;
+	stream_attributes[i].stream = -1;
+	stream_attributes[i].trib_num = -1;
+	stream_attributes[i].trib[0] = 0;
+	stream_attributes[i].trib[1] = 0;
+	stream_attributes[i].trib[2] = 0;
+	stream_attributes[i].trib[3] = 0;
+	stream_attributes[i].trib[4] = 0;
+	stream_attributes[i].cells_num = 0;
+	stream_attributes[i].init = 0;
+	stream_attributes[i].outlet = 0;
+	stream_attributes[i].length = 0.;
+	stream_attributes[i].accum_length = 0.;
+	stream_attributes[i].distance = 0.;
+	stream_attributes[i].stright = 0.;
+	stream_attributes[i].accum = 0.;
+	stream_attributes[i].init_elev = -10000;
+	stream_attributes[i].outlet_elev = -10000;
+    }
+
+    all_orders = (int **)G_malloc(orders_size * sizeof(int *));
+    for (i = 0; i < orders_size; ++i)
+	all_orders[i] =
+	    (int *)G_malloc((number_of_streams + 1) * sizeof(int));
+
+    for (i = 0; i < orders_size; ++i)
+	for (j = 0; j < number_of_streams + 1; ++j)
+	    all_orders[i][j] = -1;
+
+    return 0;
+}
+
+int compare_inits(const void *a, const void *b)
+{
+    return (*(STREAM **) a)->init - (*(STREAM **) b)->init;
+}
+
+int compare_outlets(const void *a, const void *b)
+{
+    return (*(STREAM **) a)->outlet - (*(STREAM **) b)->outlet;
+}
+
+int stream_sample_map(char *input_map_name, int number_of_streams, int what)
+{
+
+    /* if what == 2 sample outputs only for accum map */
+    /* if what == 1 sample outputs only for elev map */
+    /* if what == 0 sample inits only for elev map) */
+
+    char *mapset;
+    int input_map_fd;
+    int i;
+    int r, c, cur_row = -1;
+    RASTER_MAP_TYPE input_data_type;
+    size_t input_data_size;
+    void *input_buffer;
+    void *input_ptr;
+
+    STREAM **pointers_to_stream;
+
+    pointers_to_stream =
+	(STREAM **) G_malloc(sizeof(STREAM *) * (number_of_streams));
+    for (i = 0; i < (number_of_streams); ++i)
+	*(pointers_to_stream + i) = stream_attributes + i;
+
+    if (what)
+	qsort(pointers_to_stream, (number_of_streams), sizeof(STREAM *),
+	      compare_outlets);
+    else
+	qsort(pointers_to_stream, (number_of_streams), sizeof(STREAM *),
+	      compare_inits);
+
+
+    mapset = (char *)G_find_raster2(input_map_name, "");
+    if (mapset == NULL)
+	G_fatal_error(_("Raster map <%s> not found"), input_map_name);
+    input_map_fd = Rast_open_old(input_map_name, mapset);
+    input_data_type = Rast_map_type(input_map_name, mapset);
+    input_data_size = Rast_cell_size(input_data_type);
+    input_buffer = Rast_allocate_buf(input_data_type);
+
+    for (i = 0; i < (number_of_streams); ++i) {
+	if (pointers_to_stream[i]->stream == -1)
+	    continue;
+
+	if (what) {
+	    r = (int)pointers_to_stream[i]->outlet / ncols;
+	    c = (int)pointers_to_stream[i]->outlet % ncols;
+	}
+	else {
+	    r = (int)pointers_to_stream[i]->init / ncols;
+	    c = (int)pointers_to_stream[i]->init % ncols;
+	}
+
+	if (r > cur_row) {
+	    Rast_get_row(input_map_fd, input_buffer, r, input_data_type);
+	    input_ptr = input_buffer;
+	}
+	switch (what) {
+	case 0:		/* inits for elev */
+	    pointers_to_stream[i]->init_elev =
+		Rast_get_d_value(input_ptr + c * input_data_size,
+				 input_data_type);
+	    break;
+	case 1:		/* outlets for elev */
+	    pointers_to_stream[i]->outlet_elev =
+		Rast_get_d_value(input_ptr + c * input_data_size,
+				 input_data_type);
+	    break;
+	case 2:		/* outlets for accum */
+	    pointers_to_stream[i]->accum =
+		Rast_get_d_value(input_ptr + c * input_data_size,
+				 input_data_type);
+	    break;
+	}
+    }
+    G_free(input_buffer);
+    G_free(pointers_to_stream);
+    return 0;
+}

raster/r.stream.order/stream_order.c

@@ -0,0 +1,279 @@
+/* 
+   All algorithms used in analysis ar not recursive. For Strahler order and Shreve magnitude starts from initial channel and  proceed downstream. Algortitms try to assgin order for branch and if it is imposible start from next initial channel, till all branches are ordered.
+   For Hortor and Hack ordering it proceed upstream and uses stack data structure to determine unordered branch. 
+   Algorithm of Hack main stram according idea of Markus Metz.
+ */
+
+#include "local_proto.h"
+int strahler(int *strahler)
+{
+
+    int i, j, done = 1;
+    int cur_stream, next_stream;
+    int max_strahler = 0, max_strahler_num;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_message(_("Calculating Strahler's stream order..."));
+
+    for (j = 0; j < init_num; ++j) {	/* main loop on inits */
+
+	cur_stream = SA[init_streams[j]].stream;
+	do {			/* we must go at least once, if stream is of first order and is outlet */
+	    max_strahler_num = 1;
+	    max_strahler = 0;
+	    next_stream = SA[cur_stream].next_stream;
+
+	    if (SA[cur_stream].trib_num == 0) {	/* assign 1 for spring stream */
+		strahler[cur_stream] = 1;
+		cur_stream = next_stream;
+		done = 1;
+	    }
+	    else {
+		done = 1;
+
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {	/* loop for determining strahler */
+		    if (strahler[SA[cur_stream].trib[i]] < 0) {
+			done = 0;
+			break;	/* strahler is not determined, break for loop */
+		    }
+		    else if (strahler[SA[cur_stream].trib[i]] > max_strahler) {
+			max_strahler = strahler[SA[cur_stream].trib[i]];
+			max_strahler_num = 1;
+		    }
+		    else if (strahler[SA[cur_stream].trib[i]] == max_strahler) {
+			++max_strahler_num;
+		    }
+		}		/* end determining strahler */
+
+		if (done == 1) {
+		    strahler[cur_stream] = (max_strahler_num > 1) ?
+			++max_strahler : max_strahler;
+		    cur_stream = next_stream;	/* if next_stream<0 we in outlet stream */
+		}
+
+	    }
+	} while (done && next_stream > 0);
+    }				/* end for of main loop */
+    return 0;
+}				/* end strahler */
+
+int shreve(int *shreve)
+{
+
+    int i, j, done = 1;
+    int cur_stream, next_stream;
+    int max_shreve = 0;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_message(_("Calculating Shreve's stream magnitude, "
+                "Scheidegger's consistent integer and Drwal's streams hierarchy (old style)..."));
+
+    for (j = 0; j < init_num; ++j) {	/* main loop on inits */
+
+	cur_stream = SA[init_streams[j]].stream;
+	do {			/* we must go at least once, if stream is of first order and is outlet */
+
+	    max_shreve = 0;
+	    next_stream = SA[cur_stream].next_stream;
+
+	    if (SA[cur_stream].trib_num == 0) {	/* assign 1 for spring stream */
+
+		shreve[cur_stream] = 1;
+		cur_stream = next_stream;
+		done = 1;
+
+	    }
+	    else {
+		done = 1;
+
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {	/* loop for determining shreve */
+		    if (shreve[SA[cur_stream].trib[i]] < 0) {
+			done = 0;
+			break;	/* shreeve is not determined, break for loop */
+		    }
+		    else {
+			max_shreve += shreve[SA[cur_stream].trib[i]];
+		    }
+		}		/* end determining shreve */
+
+		if (done == 1) {
+		    shreve[cur_stream] = max_shreve;
+		    cur_stream = next_stream;	/* if next_stream<0 we in outlet stream */
+		}
+	    }
+
+	} while (done && next_stream > 0);
+    }				/* end main loop */
+    return 0;
+}				/* end shreeve */
+
+int horton(const int *strahler, int *horton, int number_of_streams)
+{
+
+    int *stack;
+    int stack_max = number_of_streams;
+    int top, i, j;
+    int cur_stream, cur_horton;
+    int max_strahler;
+    double max_accum, accum;
+    int up_stream = 0;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_message(_("Calculating Hortons's stream order..."));
+    stack = (int *)G_malloc(stack_max * sizeof(int));
+
+    for (j = 0; j < outlet_num; ++j) {
+	cur_stream = SA[outlet_streams[j]].stream;	/* outlet: init */
+	cur_horton = strahler[cur_stream];
+	stack[0] = 0;
+	stack[1] = cur_stream;
+	top = 1;
+
+	do {			/* on every stream */
+	    max_strahler = 0;
+	    max_accum = 0;
+
+	    if (SA[cur_stream].trib_num == 0) {	/* spring: go back on stack */
+
+		horton[cur_stream] = cur_horton;
+		cur_stream = stack[--top];
+
+	    }
+	    else if (SA[cur_stream].trib_num > 1) {	/* node */
+
+		up_stream = 0;	/* calculating up_stream */
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {
+		    if (horton[SA[cur_stream].trib[i]] < 0) {
+
+			if (strahler[SA[cur_stream].trib[i]] > max_strahler) {
+			    max_strahler = strahler[SA[cur_stream].trib[i]];
+			    max_accum =
+				(use_accum) ? SA[SA[cur_stream].trib[i]].
+				accum : SA[SA[cur_stream].trib[i]].
+				accum_length;
+			    up_stream = SA[cur_stream].trib[i];
+
+			}
+			else if (strahler[SA[cur_stream].trib[i]] ==
+				 max_strahler) {
+
+			    accum =
+				(use_accum) ? SA[SA[cur_stream].trib[i]].
+				accum : SA[SA[cur_stream].trib[i]].
+				accum_length;
+
+			    if (accum > max_accum) {
+				max_accum =
+				    (use_accum) ? SA[SA[cur_stream].trib[i]].
+				    accum : SA[SA[cur_stream].trib[i]].
+				    accum_length;
+
+				up_stream = SA[cur_stream].trib[i];
+			    }
+			}
+		    }
+		}		/* end determining up_stream */
+
+		if (up_stream) {	/* at least one branch is not assigned */
+		    if (horton[cur_stream] < 0) {
+			horton[cur_stream] = cur_horton;
+		    }
+		    else {
+			cur_horton = strahler[up_stream];
+		    }
+		    cur_stream = up_stream;
+		    stack[++top] = cur_stream;
+
+		}
+		else {		/* all asigned, go downstream */
+		    cur_stream = stack[--top];
+
+		}		/* end up_stream */
+	    }			/* end spring/node */
+	} while (cur_stream);
+    }				/* end for outlets */
+    G_free(stack);
+    return 0;
+}
+
+int hack(int *hack, int *topo_dim, int number_of_streams)
+{				/* also calculate topological dimension */
+
+    int *stack;
+    int top, i, j;
+    int cur_stream, cur_hack;
+    double accum, max_accum;
+    int up_stream = 0;
+    int stack_max = number_of_streams;
+    double cur_distance = 0;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_message(_("Calculating Hack's main streams and topological dimension..."));
+    stack = (int *)G_malloc(stack_max * sizeof(int));
+
+    for (j = 0; j < outlet_num; ++j) {
+
+	cur_stream = SA[outlet_streams[j]].stream;	/* outlet: init */
+	cur_hack = 1;
+	stack[0] = 0;
+	stack[1] = cur_stream;
+	top = 1;
+
+	topo_dim[cur_stream] = top;
+	cur_distance = SA[cur_stream].distance = SA[cur_stream].length;
+	do {
+	    max_accum = 0;
+
+	    if (SA[cur_stream].trib_num == 0) {	/* spring: go back on stack */
+
+		hack[cur_stream] = cur_hack;
+		cur_stream = stack[--top];
+
+	    }
+	    else if (SA[cur_stream].trib_num > 1) {	/* node */
+		up_stream = 0;	/* calculating up_stream */
+
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {	/* determining upstream */
+		    if (hack[SA[cur_stream].trib[i]] < 0) {
+
+			accum =
+			    (use_accum) ? SA[SA[cur_stream].trib[i]].
+			    accum : SA[SA[cur_stream].trib[i]].accum_length;
+
+			if (accum > max_accum) {
+			    max_accum =
+				(use_accum) ? SA[SA[cur_stream].trib[i]].
+				accum : SA[SA[cur_stream].trib[i]].
+				accum_length;
+			    up_stream = SA[cur_stream].trib[i];
+			}
+		    }
+		}		/* end determining up_stream */
+
+		if (up_stream) {	/* at least one branch is not assigned */
+
+		    if (hack[cur_stream] < 0) {
+			hack[cur_stream] = cur_hack;
+		    }
+		    else {
+			cur_hack = hack[cur_stream];
+			++cur_hack;
+		    }
+
+		    cur_distance = SA[cur_stream].distance;
+		    cur_stream = up_stream;
+		    stack[++top] = cur_stream;
+		    SA[cur_stream].distance =
+			cur_distance + SA[cur_stream].length;
+		    topo_dim[cur_stream] = top;
+		}
+		else {		/* all asigned, go downstream */
+		    cur_stream = stack[--top];
+
+		}		/* end up_stream */
+	    }			/* end spring/node */
+	} while (cur_stream);
+    }				/* end for outlets */
+    G_free(stack);
+    return 0;
+}

raster/r.stream.order/stream_raster_close.c

@@ -0,0 +1,134 @@
+#include "local_proto.h"
+int ram_close_raster_order(CELL **streams, int number_of_streams,
+			   int zerofill)
+{
+
+    G_debug(3, "ram_close_raster_order(): number_of_streams=%d", number_of_streams);
+
+    G_message("Writing outpout raster maps...");
+    int *output_fd;
+    int r, c, i;
+    CELL *output_buffer, *streams_buffer;
+    struct History history;
+    size_t data_size;
+
+    output_fd = (int *)G_malloc(orders_size * sizeof(int));
+    for (i = 0; i < orders_size; ++i) {
+	if (output_map_names[i] == NULL)
+	    continue;
+	output_fd[i] = Rast_open_c_new(output_map_names[i]);
+    }
+
+    /* this is not very elegant but use for compatibility with seg version */
+
+    data_size = Rast_cell_size(CELL_TYPE);
+    output_buffer = Rast_allocate_c_buf();
+
+    for (r = 0; r < nrows; ++r) {
+	streams_buffer = streams[r];
+
+	for (i = 0; i < orders_size; ++i) {
+
+	    if (output_map_names[i] == NULL)
+		continue;
+
+	    if (zerofill)
+		memset(output_buffer, 0, ncols * data_size);
+	    else
+		Rast_set_c_null_value(output_buffer, ncols);
+
+	    for (c = 0; c < ncols; ++c)
+		if (streams_buffer[c])
+		    output_buffer[c] = all_orders[i][streams_buffer[c]];
+
+	    Rast_put_c_row(output_fd[i], output_buffer);
+	}			/* end i */
+    }
+
+    G_free(output_buffer);
+    for (i = 0; i < orders_size; ++i)
+	G_free(all_orders[i]);
+    G_free(all_orders);
+
+
+    for (i = 0; i < orders_size; ++i) {
+	if (output_map_names[i] == NULL)
+	    continue;
+	Rast_close(output_fd[i]);
+	Rast_short_history(output_map_names[i], "raster", &history);
+	Rast_command_history(&history);
+	Rast_write_history(output_map_names[i], &history);
+	/* G_message(_("%s Done"), output_map_names[i]); */
+    }
+
+    G_free(output_fd);
+    return 0;
+}
+
+
+int seg_close_raster_order(SEGMENT *streams, int number_of_streams,
+			   int zerofill)
+{
+
+    int *output_fd;
+    int r, c, i;
+    CELL *output_buffer, *streams_buffer;
+    struct History history;
+    size_t data_size;
+
+    G_debug(3, "seg_close_raster_order(): number_of_streams=%d", number_of_streams);
+
+    output_fd = (int *)G_malloc(orders_size * sizeof(int));
+    for (i = 0; i < orders_size; ++i) {
+	if (output_map_names[i] == NULL)
+	    continue;
+	output_fd[i] = Rast_open_c_new(output_map_names[i]);
+    }
+
+    data_size = Rast_cell_size(CELL_TYPE);
+    output_buffer = Rast_allocate_c_buf();
+    streams_buffer = Rast_allocate_c_buf();
+    segment_flush(streams);
+
+    for (r = 0; r < nrows; ++r) {
+	if (0 > segment_get_row(streams, streams_buffer, r))
+	    G_warning(_("Unable to segment read row %d for raster map <%s>"),
+		      r, output_map_names[i]);
+
+	for (i = 0; i < orders_size; ++i) {
+
+	    if (output_map_names[i] == NULL)
+		continue;
+
+	    if (zerofill)
+		memset(output_buffer, 0, ncols * data_size);
+	    else
+		Rast_set_c_null_value(output_buffer, ncols);
+
+	    for (c = 0; c < ncols; ++c)
+		if (!Rast_is_c_null_value(&streams_buffer[c]))
+		    output_buffer[c] = all_orders[i][streams_buffer[c]];
+
+	    Rast_put_c_row(output_fd[i], output_buffer);
+	}			/* end i */
+    }
+
+    G_free(output_buffer);
+    G_free(streams_buffer);
+    for (i = 0; i < orders_size; ++i)
+	G_free(all_orders[i]);
+    G_free(all_orders);
+
+    for (i = 0; i < orders_size; ++i) {
+	if (output_map_names[i] == NULL)
+	    continue;
+	Rast_close(output_fd[i]);
+	Rast_short_history(output_map_names[i], "raster", &history);
+	Rast_command_history(&history);
+	Rast_write_history(output_map_names[i], &history);
+	/* G_message(_("%s Done"), output_map_names[i]);*/
+    }
+
+    G_free(output_fd);
+    return 0;
+}

raster/r.stream.order/stream_topology.c

@@ -0,0 +1,431 @@
+#include "local_proto.h"
+
+int ram_number_of_tribs(int r, int c, CELL **streams, CELL **dirs)
+{
+
+    int trib = 0;
+    int i, j;
+
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+	j = DIAG(i);
+	if (streams[NR(i)][NC(i)] && dirs[NR(i)][NC(i)] == j)
+	    trib++;
+    }
+
+    if (trib > 5)
+	G_fatal_error(_("Error finding nodes. "
+                        "Stream and direction maps probably do not match."));
+    if (trib > 3)
+	G_warning(_("Stream network may be too dense"));
+
+    return trib;
+}
+
+int ram_stream_topology(CELL **streams, CELL **dirs, int number_of_streams)
+{
+
+    int d, i, j;		/* d: direction, i: iteration */
+    int r, c;
+    int next_r, next_c;
+    int trib_num, trib = 0;
+    int next_stream = -1, cur_stream;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    init_num = 0, outlet_num = 0;
+
+    G_message(_("Finding nodes..."));
+
+    outlet_streams = (unsigned int *)G_malloc((number_of_streams) *
+					      sizeof(unsigned int));
+    init_streams = (unsigned int *)G_malloc((number_of_streams) *
+					    sizeof(unsigned int));
+    init_cells = (unsigned
+		  long int *)G_malloc((number_of_streams) * sizeof(unsigned
+								   long int));
+    /* free at the end */
+
+    for (r = 0; r < nrows; ++r)
+	for (c = 0; c < ncols; ++c)
+	    if (streams[r][c]) {
+		trib_num = ram_number_of_tribs(r, c, streams, dirs);
+		trib = 0;
+		d = abs(dirs[r][c]);	/* r.watershed! */
+		if (d < 1 || NOT_IN_REGION(d) || !streams[NR(d)][NC(d)])
+		    next_stream = -1;
+		else
+		    next_stream = streams[NR(d)][NC(d)];
+
+		cur_stream = streams[r][c];
+
+		if (cur_stream != next_stream) {	/* junction: building topology */
+
+		    if (outlet_num > (number_of_streams - 1))
+			G_fatal_error(_("Error finding nodes. "
+                                        "Stream and direction maps probably do not match."));
+
+		    SA[cur_stream].stream = cur_stream;
+		    SA[cur_stream].next_stream = next_stream;
+
+		    if (next_stream < 0)	/* is outlet stream */
+			outlet_streams[outlet_num++] = cur_stream;
+		}
+
+		if (trib_num == 0) {	/* is init */
+		    if (init_num > (number_of_streams - 1))
+			G_fatal_error(_("Error finding nodes. "
+                                        "Stream and direction maps probably do not match."));
+
+		    SA[cur_stream].trib_num = 0;
+		    init_cells[init_num] = r * ncols + c;
+		    init_streams[init_num++] = cur_stream;	/* collecting inits */
+		}
+
+		if (trib_num > 1) {	/* adding tributuaries */
+		    SA[cur_stream].trib_num = trib_num;
+
+		    for (i = 1; i < 9; ++i) {
+			if (trib > 4)
+			    G_fatal_error(_("Error finding nodes. "
+                                            "Stream and direction maps probably do not match."));
+			if (NOT_IN_REGION(i))
+			    continue;
+			j = DIAG(i);
+			next_r = NR(i);
+			next_c = NC(i);
+			if (streams[next_r][next_c] &&
+			    dirs[next_r][next_c] == j)
+			    SA[cur_stream].trib[trib++] =
+				streams[next_r][next_c];
+		    }		/* end for i... */
+		}
+	    }			/* end if streams */
+    return 0;
+}
+
+int ram_stream_geometry(CELL **streams, CELL **dirs)
+{
+
+    int i, s, d;		/* s - streams index; d - direction */
+    int done = 1;
+    int r, c;
+    int next_r, next_c;
+    int prev_r, prev_c;
+    int cur_stream;
+    float cur_northing, cur_easting;
+    float next_northing, next_easting;
+    float init_northing, init_easting;
+    double cur_length = 0.;
+    double cur_accum_length = 0.;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+    struct Cell_head window;
+
+    G_get_window(&window);
+
+    G_message(_("Finding longest streams..."));
+    G_begin_distance_calculations();
+
+    for (s = 0; s < init_num; ++s) {	/* main loop on springs */
+	r = (int)init_cells[s] / ncols;
+	c = (int)init_cells[s] % ncols;
+	cur_stream = streams[r][c];
+	cur_length = 0;
+	done = 1;
+
+	SA[cur_stream].init = init_cells[s];	/* stored as index */
+
+	init_northing = window.north - (r + .5) * window.ns_res;
+	init_easting = window.west + (c + .5) * window.ew_res;
+
+	while (done) {
+	    cur_northing = window.north - (r + .5) * window.ns_res;
+	    cur_easting = window.west + (c + .5) * window.ew_res;
+
+	    d = abs(dirs[r][c]);
+	    next_r = NR(d);
+	    next_c = NC(d);
+
+	    if (d < 1 || NOT_IN_REGION(d) || !streams[next_r][next_c]) {
+		cur_length = (window.ns_res + window.ew_res) / 2;
+		SA[cur_stream].accum_length += cur_length;
+		SA[cur_stream].length += cur_length;
+		SA[cur_stream].stright =
+		    G_distance(cur_easting, cur_northing, init_easting,
+			       init_northing);
+		SA[cur_stream].outlet = (r * ncols + c);	/* add outlet to sorting */
+		break;
+	    }
+
+	    next_northing = window.north - (next_r + .5) * window.ns_res;
+	    next_easting = window.west + (next_c + .5) * window.ew_res;
+	    cur_length =
+		G_distance(next_easting, next_northing, cur_easting,
+			   cur_northing);
+	    SA[cur_stream].accum_length += cur_length;
+	    SA[cur_stream].length += cur_length;
+	    prev_r = r;
+	    prev_c = c;
+	    r = next_r;
+	    c = next_c;
+
+	    if (streams[next_r][next_c] != cur_stream) {
+		SA[cur_stream].stright =
+		    G_distance(next_easting, next_northing, init_easting,
+			       init_northing);
+		init_northing = cur_northing;
+		init_easting = cur_easting;
+
+		SA[cur_stream].outlet = (prev_r * ncols + prev_c);
+		cur_stream = streams[next_r][next_c];
+
+		cur_accum_length = 0;
+		SA[cur_stream].init = (r * ncols + c);
+
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {
+		    if (SA[SA[cur_stream].trib[i]].accum_length == 0) {
+			done = 0;
+			cur_accum_length = 0;
+			break;	/* do not pass accum */
+		    }
+		    if (SA[SA[cur_stream].trib[i]].accum_length >
+			cur_accum_length)
+			cur_accum_length =
+			    SA[SA[cur_stream].trib[i]].accum_length;
+		}		/* end for i */
+		SA[cur_stream].accum_length = cur_accum_length;
+	    }			/* end if */
+	}			/* end while */
+    }				/* end for s */
+    return 0;
+}
+
+int seg_number_of_tribs(int r, int c, SEGMENT *streams, SEGMENT *dirs)
+{
+
+    int trib = 0;
+    int i, j;
+    int streams_cell = 0;
+    int dirs_cell = 0;
+
+    for (i = 1; i < 9; ++i) {
+	if (NOT_IN_REGION(i))
+	    continue;
+
+	j = DIAG(i);
+
+	segment_get(streams, &streams_cell, NR(i), NC(i));
+	segment_get(dirs, &dirs_cell, NR(i), NC(i));
+
+	if (streams_cell && dirs_cell == j)
+	    trib++;
+    }
+
+    if (trib > 5)
+	G_fatal_error(_("Error finding nodes. "
+                        "Stream and direction maps probably do not match."));
+    if (trib > 3)
+	G_warning(_("Stream network may be too dense"));
+
+    return trib;
+}
+
+int seg_stream_topology(SEGMENT *streams, SEGMENT *dirs,
+			int number_of_streams)
+{
+
+    int d, i, j;		/* d: direction, i: iteration */
+    int r, c;
+    int next_r, next_c;
+    int trib_num, trib = 0;
+    int next_stream = -1, cur_stream;
+    int streams_cell, dirs_cell;
+    int next_streams_cell, trib_dirs_cell, trib_stream_cell;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    init_num = 0, outlet_num = 0;
+
+    G_message(_("Finding nodes..."));
+
+    outlet_streams = (unsigned int *)G_malloc((number_of_streams) *
+					      sizeof(unsigned int));
+    init_streams = (unsigned int *)G_malloc((number_of_streams) *
+					    sizeof(unsigned int));
+    init_cells = (unsigned
+		  long int *)G_malloc((number_of_streams) * sizeof(unsigned
+								   long int));
+
+    for (r = 0; r < nrows; ++r) {
+	G_percent(r, nrows, 2);
+	for (c = 0; c < ncols; ++c) {
+	    segment_get(streams, &streams_cell, r, c);
+	    segment_get(dirs, &dirs_cell, r, c);
+
+	    if (streams_cell) {
+		trib_num = seg_number_of_tribs(r, c, streams, dirs);
+		trib = 0;
+
+		d = abs(dirs_cell);	/* r.watershed! */
+		if (NOT_IN_REGION(d))
+		    next_stream = -1;
+		else
+		    segment_get(streams, &next_streams_cell, NR(d), NC(d));
+
+		if (d < 1 || NOT_IN_REGION(d) || !next_streams_cell)
+		    next_stream = -1;
+		else
+		    segment_get(streams, &next_stream, NR(d), NC(d));
+
+		cur_stream = streams_cell;
+
+		if (cur_stream != next_stream) {	/* junction: building topology */
+		    if (outlet_num > (number_of_streams - 1))
+			G_fatal_error(_("Error finding nodes. "
+                                        "Stream and direction maps probably do not match."));
+
+		    SA[cur_stream].stream = cur_stream;
+		    SA[cur_stream].next_stream = next_stream;
+
+		    if (next_stream < 0)	/* is outlet stream */
+			outlet_streams[outlet_num++] = cur_stream;
+		}
+
+		if (trib_num == 0) {	/* is init */
+		    if (init_num > (number_of_streams - 1))
+			G_fatal_error(_("Error finding nodes. "
+                                        "Stream and direction maps probably do not match."));
+
+		    SA[cur_stream].trib_num = 0;
+		    init_cells[init_num] = r * ncols + c;
+		    init_streams[init_num++] = cur_stream;	/* collecting inits */
+		}
+
+		if (trib_num > 1) {	/* adding tributuaries */
+		    SA[cur_stream].trib_num = trib_num;
+
+		    for (i = 1; i < 9; ++i) {
+
+			if (trib > 4)
+			    G_fatal_error(_("Error finding nodes. "
+                                            "Stream and direction maps probably do not match."));
+			if (NOT_IN_REGION(i))
+			    continue;
+			j = DIAG(i);
+			next_r = NR(i);
+			next_c = NC(i);
+			segment_get(streams, &trib_stream_cell, next_r,
+				    next_c);
+			segment_get(dirs, &trib_dirs_cell, next_r, next_c);
+
+			if (trib_stream_cell && trib_dirs_cell == j)
+			    SA[cur_stream].trib[trib++] = trib_stream_cell;
+		    }		/* end for i... */
+		}
+	    }			/* end if streams */
+	}
+    }				/* end r, c */
+    G_percent(r, nrows, 2);
+    return 0;
+}
+
+int seg_stream_geometry(SEGMENT *streams, SEGMENT *dirs)
+{
+
+    int i, s, d;		/* s - streams index; d - direction */
+    int done = 1;
+    int r, c;
+    int next_r, next_c;
+    int prev_r, prev_c;
+    int cur_stream, next_stream, dirs_cell;
+    float cur_northing, cur_easting;
+    float next_northing, next_easting;
+    float init_northing, init_easting;
+    double cur_length = 0.;
+    double cur_accum_length = 0.;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+    struct Cell_head window;
+
+    G_get_window(&window);
+
+    G_message(_("Finding longest streams..."));
+    G_begin_distance_calculations();
+
+    for (s = 0; s < init_num; ++s) {	/* main loop on springs */
+	G_percent(s, init_num, 2);
+	r = (int)init_cells[s] / ncols;
+	c = (int)init_cells[s] % ncols;
+	segment_get(streams, &cur_stream, r, c);
+	cur_length = 0;
+	done = 1;
+
+	SA[cur_stream].init = init_cells[s];	/* stored as index */
+
+	init_northing = window.north - (r + .5) * window.ns_res;
+	init_easting = window.west + (c + .5) * window.ew_res;
+
+	while (done) {
+	    cur_northing = window.north - (r + .5) * window.ns_res;
+	    cur_easting = window.west + (c + .5) * window.ew_res;
+
+	    segment_get(dirs, &dirs_cell, r, c);
+	    d = abs(dirs_cell);
+	    next_r = NR(d);
+	    next_c = NC(d);
+	    if (NOT_IN_REGION(d))
+		Rast_set_c_null_value(&next_stream, 1);
+	    else
+		segment_get(streams, &next_stream, next_r, next_c);
+
+	    if (d < 1 || NOT_IN_REGION(d) || !next_stream) {
+		cur_length = (window.ns_res + window.ew_res) / 2;
+		SA[cur_stream].accum_length += cur_length;
+		SA[cur_stream].length += cur_length;
+		SA[cur_stream].stright =
+		    G_distance(cur_easting, cur_northing, init_easting,
+			       init_northing);
+		SA[cur_stream].outlet = (r * ncols + c);	/* add outlet to sorting */
+		break;
+	    }
+
+	    next_northing = window.north - (next_r + .5) * window.ns_res;
+	    next_easting = window.west + (next_c + .5) * window.ew_res;
+	    cur_length =
+		G_distance(next_easting, next_northing, cur_easting,
+			   cur_northing);
+	    SA[cur_stream].accum_length += cur_length;
+	    SA[cur_stream].length += cur_length;
+	    prev_r = r;
+	    prev_c = c;
+	    r = next_r;
+	    c = next_c;
+
+	    if (next_stream != cur_stream) {
+		SA[cur_stream].stright =
+		    G_distance(next_easting, next_northing, init_easting,
+			       init_northing);
+		init_northing = cur_northing;
+		init_easting = cur_easting;
+
+		SA[cur_stream].outlet = (prev_r * ncols + prev_c);
+		cur_stream = next_stream;
+		cur_accum_length = 0;
+		SA[cur_stream].init = (r * ncols + c);
+
+		for (i = 0; i < SA[cur_stream].trib_num; ++i) {
+		    if (SA[SA[cur_stream].trib[i]].accum_length == 0) {
+			done = 0;
+			cur_accum_length = 0;
+			break;	/* do not pass accum */
+		    }
+		    if (SA[SA[cur_stream].trib[i]].accum_length >
+			cur_accum_length)
+			cur_accum_length =
+			    SA[SA[cur_stream].trib[i]].accum_length;
+		}		/* end for i */
+		SA[cur_stream].accum_length = cur_accum_length;
+	    }			/* end if */
+	}			/* end while */
+    }				/* end for s */
+    G_percent(1, 1, 1);
+    return 0;
+}

raster/r.stream.order/stream_vector.c

@@ -0,0 +1,362 @@
+#include "local_proto.h"
+int ram_create_vector(CELL ** streams, CELL ** dirs, char *out_vector,
+		      int number_of_streams)
+{
+
+    int i, d;
+    int r, c;
+    int next_r, next_c;
+    int add_outlet = 0;
+    int cur_stream;
+    float northing, easting;
+    struct Cell_head window;
+    struct line_pnts *Segments;
+    struct line_cats *Cats;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_get_window(&window);
+    Segments = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    Vect_open_new(&Out, out_vector, 0);
+
+    Vect_reset_line(Segments);
+    Vect_reset_cats(Cats);
+
+    for (i = 0; i < number_of_streams; ++i) {
+
+	if (SA[i].stream == -1)
+	    continue;		/* empty category */
+
+	add_outlet = 0;
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+
+	cur_stream = SA[i].stream;
+	Vect_cat_set(Cats, 1, cur_stream);
+	easting = window.west + (c + .5) * window.ew_res;
+	northing = window.north - (r + .5) * window.ns_res;
+	Vect_append_point(Segments, easting, northing, 0);
+	Vect_write_line(&Out, GV_POINT, Segments, Cats);
+	Vect_reset_line(Segments);
+	Vect_append_point(Segments, easting, northing, 0);
+
+	while (streams[r][c] == cur_stream) {
+
+	    d = abs(dirs[r][c]);
+	    next_r = NR(d);
+	    next_c = NC(d);
+
+	    easting = window.west + (next_c + .5) * window.ew_res;
+	    northing = window.north - (next_r + .5) * window.ns_res;
+	    Vect_append_point(Segments, easting, northing, 0);
+
+	    if (d < 1 || NOT_IN_REGION(d) || !streams[next_r][next_c]) {
+		add_outlet = 1;
+		break;
+	    }
+	    r = next_r;
+	    c = next_c;
+	}			/* end while */
+
+	Vect_cat_set(Cats, 1, cur_stream);
+	Vect_write_line(&Out, GV_LINE, Segments, Cats);
+	Vect_reset_line(Segments);
+	Vect_reset_cats(Cats);
+
+	if (add_outlet) {
+	    Vect_cat_set(Cats, 1, 0);
+	    Vect_reset_line(Segments);
+	    Vect_append_point(Segments, easting, northing, 0);
+	    Vect_write_line(&Out, GV_POINT, Segments, Cats);
+	    Vect_reset_line(Segments);
+	    Vect_reset_cats(Cats);
+	}
+    }
+
+    /* build vector after adding table */
+    if (0 < stream_add_table(number_of_streams))
+	G_warning(_("Unable to add attribute table to vector map <%s>"), out_vector);
+    Vect_hist_command(&Out);
+    Vect_build(&Out);
+    Vect_close(&Out);
+
+    return 0;
+}
+
+int seg_create_vector(SEGMENT * streams, SEGMENT * dirs, char *out_vector,
+		      int number_of_streams)
+{
+
+    int i, d;
+    int r, c;
+    int next_r, next_c;
+    int add_outlet;
+    int streams_cell, dirs_cell;
+    int cur_stream, next_stream;
+    float northing, easting;
+    struct Cell_head window;
+    struct line_pnts *Segments;
+    struct line_cats *Cats;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+
+    G_get_window(&window);
+    Segments = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    Vect_open_new(&Out, out_vector, 0);
+
+    Vect_reset_line(Segments);
+    Vect_reset_cats(Cats);
+
+    for (i = 0; i < number_of_streams; ++i) {
+	if (SA[i].stream == -1)
+	    continue;
+
+	add_outlet = 0;
+	r = (int)SA[i].init / ncols;
+	c = (int)SA[i].init % ncols;
+
+	cur_stream = SA[i].stream;
+	Vect_cat_set(Cats, 1, cur_stream);
+	easting = window.west + (c + .5) * window.ew_res;
+	northing = window.north - (r + .5) * window.ns_res;
+	Vect_append_point(Segments, easting, northing, 0);
+	Vect_write_line(&Out, GV_POINT, Segments, Cats);
+	Vect_reset_line(Segments);
+	Vect_append_point(Segments, easting, northing, 0);
+
+	segment_get(streams, &streams_cell, r, c);
+	while (streams_cell == cur_stream) {
+
+	    segment_get(dirs, &dirs_cell, r, c);
+	    d = abs(dirs_cell);
+	    next_r = NR(d);
+	    next_c = NC(d);
+
+	    easting = window.west + (next_c + .5) * window.ew_res;
+	    northing = window.north - (next_r + .5) * window.ns_res;
+	    Vect_append_point(Segments, easting, northing, 0);
+
+	    if (NOT_IN_REGION(d))
+		Rast_set_c_null_value(&next_stream, 1);
+	    else
+		segment_get(streams, &next_stream, next_r, next_c);
+
+	    if (d < 1 || NOT_IN_REGION(d) || !next_stream) {
+		add_outlet = 1;
+		break;
+	    }
+	    r = next_r;
+	    c = next_c;
+	    segment_get(streams, &streams_cell, r, c);
+	}			/* end while */
+
+	Vect_cat_set(Cats, 1, cur_stream);
+	Vect_write_line(&Out, GV_LINE, Segments, Cats);
+	Vect_reset_line(Segments);
+	Vect_reset_cats(Cats);
+
+	if (add_outlet) {
+	    Vect_cat_set(Cats, 1, 0);
+	    Vect_reset_line(Segments);
+	    Vect_append_point(Segments, easting, northing, 0);
+	    Vect_write_line(&Out, GV_POINT, Segments, Cats);
+	    Vect_reset_line(Segments);
+	    Vect_reset_cats(Cats);
+	}
+    }
+
+    /* build vector after adding table */
+    if (0 < stream_add_table(number_of_streams))
+	G_warning(_("Unable to add attribute table to vector map <%s>"), out_vector);
+    Vect_hist_command(&Out);
+    Vect_build(&Out);
+    Vect_close(&Out);
+    return 0;
+}
+
+int stream_add_table(int number_of_streams)
+{
+
+    int i;
+    int max_trib = 0;
+    struct field_info *Fi;
+    STREAM *SA = stream_attributes;	/* for better code readability */
+    dbDriver *driver;
+    dbHandle handle;
+    dbString table_name, db_sql, val_string;
+    char *cat_col_name = "cat";
+    char buf[1000];
+    char ins_prev_streams[50];	/* insert */
+
+    /* rest of table definition */
+    char *tab_cat_col_name = "cat integer";
+    char *tab_stream = "stream integer";
+    char *tab_next_stream = "next_stream integer";
+    char *tab_prev_streams;
+    char *tab_orders =
+	"strahler integer, horton integer, shreve integer, hack integer, topo_dim integer";
+    char *tab_scheidegger = "scheidegger integer";
+    char *tab_drwal_old = "drwal_old integer";
+    char *tab_length = "length double precision";
+    char *tab_stright = "stright double precision";
+    char *tab_sinusoid = "sinosoid double precision";
+    char *tab_cumlength = "cum_length double precision";
+    char *tab_accum = "flow_accum double precision";
+    char *tab_distance = "out_dist double precision";
+    char *tab_elev_init = "source_elev double precision";
+    char *tab_elev_outlet = "outlet_elev double precision";
+    char *tab_drop = "elev_drop double precision";
+    char *tab_out_drop = "out_drop double precision";
+    char *tab_gradient = "gradient double precision";
+
+    /* addational atrributes */
+    int scheidegger, drwal_old = -1;
+    double sinusoid = 1, elev_drop, out_drop = 0, gradient = -1;
+    char insert_orders[60];	/* must have to be increased if new orders are added */
+
+    db_init_string(&db_sql);
+    db_init_string(&val_string);
+    db_init_string(&table_name);
+    db_init_handle(&handle);
+
+    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+    driver = db_start_driver_open_database(Fi->driver,
+					   Vect_subst_var(Fi->database,
+							         &Out));
+
+    /* create table */
+    for (i = 0; i < number_of_streams; ++i)
+	if (SA[i].trib_num > max_trib)
+	    max_trib = SA[i].trib_num;
+
+    switch (max_trib) {
+    case 2:
+	tab_prev_streams = "prev_str01 integer, prev_str02 integer";
+	break;
+    case 3:
+	tab_prev_streams =
+	    "prev_str01 integer, prev_str02 integer, prev_str03 integer";
+	break;
+    case 4:
+	tab_prev_streams =
+	    "prev_str01 integer, prev_str02 integer, prev_str03 integer, prev_str04 integer";
+	break;
+    case 5:
+	tab_prev_streams =
+	    "prev_str01 integer, prev_str02 integer, prev_str03 integer, prev_str04 integer, prev_str05 integer";
+	break;
+    default:
+	G_fatal_error(_("Error with number of tributuaries"));
+	break;
+    }
+
+    sprintf(buf, "create table %s (%s, %s, %s, %s, %s,"
+            "%s, %s, %s, %s, %s,"                                                     
+            "%s, %s, %s, %s, %s,"                                         
+            "%s, %s, %s)", Fi->table, tab_cat_col_name,	/* 1 */
+	    tab_stream, tab_next_stream, tab_prev_streams, tab_orders,	/* 5 */
+	    tab_scheidegger, tab_drwal_old, tab_length, tab_stright, tab_sinusoid,	/* 10 */
+	    tab_cumlength, tab_accum, tab_distance, tab_elev_init, tab_elev_outlet,	/* 15 */
+	    tab_drop, tab_out_drop, tab_gradient	/* 18 */
+	);
+
+    db_set_string(&db_sql, buf);
+
+    if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+	db_close_database(driver);
+	db_shutdown_driver(driver);
+	G_warning("Unable to create table: '%s'", db_get_string(&db_sql));
+	return -1;
+    }
+
+    if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+	G_warning(_("Unable to create index on table <%s>"), Fi->table);
+
+    if (db_grant_on_table(driver, Fi->table,
+			  DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK) {
+	G_warning(_("Unable to grant privileges on table <%s>"), Fi->table);
+	return -1;
+    }
+    db_begin_transaction(driver);
+
+    for (i = 0; i < number_of_streams; ++i) {
+
+	if (SA[i].stream < 0)
+	    continue;
+
+	/* calc addational parameters */
+
+	scheidegger = (all_orders[o_shreve][i]) * 2;
+
+	if (all_orders[o_shreve][i] > 0)
+	    drwal_old = (int)(log(all_orders[o_shreve][i]) / log(2)) + 1;
+
+	sinusoid = -1;
+	if (SA[i].stright > 0)
+	    sinusoid = SA[i].length / SA[i].stright;
+
+	out_drop = 0;
+	if (SA[i].next_stream > 0)
+	    out_drop = SA[i].outlet_elev - SA[SA[i].next_stream].init_elev;
+
+	elev_drop = (SA[i].init_elev - SA[i].outlet_elev) + out_drop;
+	if (elev_drop < 0)
+	    elev_drop = 0;
+
+	gradient = -1;
+	if (SA[i].length > 0)
+	    gradient = elev_drop / SA[i].length;
+
+	switch (max_trib) {
+	case 2:
+	    sprintf(ins_prev_streams, "%d, %d", SA[i].trib[0], SA[i].trib[1]);
+	    break;
+	case 3:
+	    sprintf(ins_prev_streams, "%d ,%d, %d", SA[i].trib[0],
+		    SA[i].trib[1], SA[i].trib[2]);
+	    break;
+	case 4:
+	    sprintf(ins_prev_streams, "%d, %d, %d, %d", SA[i].trib[0],
+		    SA[i].trib[1], SA[i].trib[2], SA[i].trib[3]);
+	    break;
+	case 5:
+	    sprintf(ins_prev_streams, "%d, %d, %d, %d, %d", SA[i].trib[0],
+		    SA[i].trib[1], SA[i].trib[2], SA[i].trib[3],
+		    SA[i].trib[4]);
+	    break;
+	default:
+            G_fatal_error(_("Error with number of tributuaries"));
+	    break;
+	}
+
+	sprintf(insert_orders, "%d, %d, %d, %d, %d",
+		all_orders[0][i], all_orders[1][i],
+		all_orders[2][i], all_orders[3][i], all_orders[4][i]);
+
+	sprintf(buf, "insert into %s values( %d, %d, %d, %s, %s, "
+                "%d, %d, %f, %f, %f, "
+                "%f, %f, %f, %f, %f, "
+                "%f, %f, %f)", Fi->table, i,	/* 1 */
+		SA[i].stream, SA[i].next_stream, ins_prev_streams,	/* buffer created before */
+		insert_orders,	/* 5 *//* buffer created before */
+		scheidegger, drwal_old, SA[i].length, SA[i].stright, sinusoid,	/* 10 */
+		SA[i].accum_length, fabs(SA[i].accum), SA[i].distance, SA[i].init_elev, SA[i].outlet_elev,	/* 15 */
+		elev_drop, out_drop, gradient	/* 18 */
+	    );
+
+	db_set_string(&db_sql, buf);
+
+	if (db_execute_immediate(driver, &db_sql) != DB_OK) {
+	    db_close_database(driver);
+	    db_shutdown_driver(driver);
+	    G_warning(_("Unable to inset new row: '%s'"), db_get_string(&db_sql));
+	    return -1;
+	}
+    }				/* end for */
+
+    db_commit_transaction(driver);
+    db_close_database_shutdown_driver(driver);
+    Vect_map_add_dblink(&Out, 1, NULL, Fi->table,
+			cat_col_name, Fi->database, Fi->driver);
+    return 0;
+}