[GRASS-SVN] r43754 -
grass-addons/grass7/raster/r.stream/r.stream.channel
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Oct 1 07:20:30 EDT 2010
Author: jarekj71
Date: 2010-10-01 11:20:30 +0000 (Fri, 01 Oct 2010)
New Revision: 43754
Added:
grass-addons/grass7/raster/r.stream/r.stream.channel/Makefile
grass-addons/grass7/raster/r.stream/r.stream.channel/io.c
grass-addons/grass7/raster/r.stream/r.stream.channel/io.h
grass-addons/grass7/raster/r.stream/r.stream.channel/local_proto.h
grass-addons/grass7/raster/r.stream/r.stream.channel/local_vars.h
grass-addons/grass7/raster/r.stream/r.stream.channel/main.c
grass-addons/grass7/raster/r.stream/r.stream.channel/r.stream.channel.html
grass-addons/grass7/raster/r.stream/r.stream.channel/stream_topology.c
grass-addons/grass7/raster/r.stream/r.stream.channel/stream_write.c
grass-addons/grass7/raster/r.stream/r.stream.channel/svn-commit.tmp
Log:
files
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/Makefile 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.stream.channel
+
+LIBES = $(GISLIB) $(RASTERLIB) $(SEGMENTLIB) $(VECTLIB) $(DBMILIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP) $(SEGMENTDEP) $(VECTDEP) $(DBMIDEP)
+
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/io.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/io.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/io.c 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,516 @@
+#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, c;
+
+ if(data_type < 0 || data_type > 2)
+ G_fatal_error(_("ram_creat: Cannot create 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 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(_("<%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 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(_("ram_open: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) {
+ 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 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,"ram_null:Cannot 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(_("seg_create: unrecognisabe 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(_("seg_create: cannot format segment"));
+ }
+
+ close(fd);
+ if(0 > (fd = open(filename,2))) {
+ unlink(filename);
+ G_fatal_error(_("seg_create: cannot re-open file"));
+ }
+
+ if(0>(fd = segment_init(&(seg->seg),fd,number_of_segs))) {
+ unlink(filename);
+ G_fatal_error(_("seg_create: cannot 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(_("seg_read: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 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(_("<%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 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(_("seg_read: Cannot segment put row %d for 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);
+ }
+
+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 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(_("seg_write: Cannot segment read row %d for 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(_("ram_null:Cannot 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;
+}
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/io.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/io.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/io.h 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,56 @@
+#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/segment.h>
+#include <grass/vector.h>
+#include <grass/dbmi.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(x,y) (x)*ncols+(y)
+#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_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*);
+
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/local_proto.h 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,3 @@
+#include "io.h"
+#include "local_vars.h"
+
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/local_vars.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/local_vars.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/local_vars.h 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,41 @@
+#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>
+
+#define SQRT2 1.414214
+
+
+
+typedef struct {
+ int stream_num;
+ int number_of_cells;
+ int order;
+ unsigned long int * points;
+ float * elevation;
+ double * distance;
+ unsigned int init_r;
+ unsigned int init_c;
+} STREAM;
+
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+GLOBAL int nextr[9];
+GLOBAL int nextc[9];
+
+GLOBAL int nrows, ncols;
+GLOBAL STREAM* stream_attributes;
+
+GLOBAL struct Cell_head window;
+
+
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/main.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/main.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/main.c 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,323 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.stream.channel
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Program calculate some channel properities:
+ * local elevation change, curvature along stream,
+ * distance to channel init/join, elevation below channel init,
+ * optionally distance to outlet, elevation above outlet;
+
+ * COPYRIGHT: (C) 2002,2010 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 output[] = {
+ {"local_diff",NO,"Local elevation difference"},
+ {"out_dist",NO,"Upstream distance form init"},
+ */
+
+ struct GModule *module;
+ struct Option *in_dir_opt, /* options */
+ *in_stm_opt,
+ *in_elev_opt,
+ *out_identifier_opt,
+ *out_distance_opt,
+ *out_difference_opt,
+ *out_gradient_opt,
+ *out_curvature_opt,
+ *opt_swapsize;
+
+ struct Flag *flag_segmentation,
+ *flag_local,
+ *flag_cells,
+ *flag_downstream;
+
+ char *method_name[] = { "UPSTREAM", "DOWNSTREAM" };
+ int number_of_segs;
+ int number_of_streams;
+ int outlets_num;
+ int order_max;
+ int segmentation, downstream, local, cells; /*flags */
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword("Horton statisctics");
+ module->description =
+ _("Calculate local parameters for individual streams");
+ G_add_keyword("Stream parameters");
+ G_add_keyword("Stream gradient");
+ G_add_keyword("Stream curvature");
+
+ in_stm_opt = G_define_standard_option(G_OPT_R_INPUT);
+ in_stm_opt->key = "streams";
+ in_stm_opt->description = "Name of streams mask input map";
+
+ in_dir_opt = G_define_standard_option(G_OPT_R_INPUT);
+ in_dir_opt->key = "dirs";
+ in_dir_opt->description = "Name of flow direction input map";
+
+ in_elev_opt = G_define_standard_option(G_OPT_R_INPUT);
+ in_elev_opt->key = "elevation";
+ in_elev_opt->description = "Name of elevation map";
+
+ out_identifier_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_identifier_opt->key = "identifier";
+ out_identifier_opt->required = NO;
+ out_identifier_opt->description = "Unique identifier for stream";
+
+ out_distance_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_distance_opt->key = "distance";
+ out_distance_opt->required = NO;
+ out_distance_opt->description = "Distance from init/join/outlet";
+
+ out_difference_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_difference_opt->key = "difference";
+ out_difference_opt->required = NO;
+ out_difference_opt->description = "Elevation differecne from init/join/outlet";
+
+ out_gradient_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_gradient_opt->key = "gradient";
+ out_gradient_opt->required = NO;
+ out_gradient_opt->description = "Mean gradient of stream from init/join/outlet";
+
+ out_curvature_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ out_curvature_opt->key = "curvature";
+ out_curvature_opt->required = NO;
+ out_curvature_opt->description = "Local stream curvature";
+
+ 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=_("Optional");
+
+ flag_downstream = G_define_flag();
+ flag_downstream->key = 'd';
+ flag_downstream->description = _("Calculate parameters from outlet (downstream values)");
+
+ flag_local = G_define_flag();
+ flag_local->key = 'l';
+ flag_local->description = _("Calculate local values (for current cell)");
+
+ flag_cells = G_define_flag();
+ flag_cells->key = 'c';
+ flag_cells->description = _("Calculate distance in cell (ignored local)");
+
+ flag_segmentation = G_define_flag();
+ flag_segmentation->key = 'm';
+ flag_segmentation->description = _("Use memory swap (operation is slow)");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ segmentation = (flag_segmentation->answer != 0);
+ downstream = (flag_downstream->answer != 0);
+
+ if (!out_identifier_opt->answer &&
+ !out_distance_opt->answer &&
+ !out_difference_opt->answer &&
+ !out_gradient_opt->answer &&
+ !out_curvature_opt->answer)
+ G_fatal_error(_("You must select at least one output maps"));
+
+ local = (flag_local->answer != 0);
+ cells = (flag_cells->answer != 0);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ G_get_window(&window);
+ G_begin_distance_calculations();
+
+ if(out_identifier_opt->answer)
+ if (G_legal_filename(out_identifier_opt->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal difference map name"),
+ out_identifier_opt->answer);
+
+ if(out_difference_opt->answer)
+ if (G_legal_filename(out_difference_opt->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal difference map name"),
+ out_difference_opt->answer);
+
+ if(out_distance_opt->answer)
+ if (G_legal_filename(out_distance_opt->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal distance map name"),
+ out_distance_opt->answer);
+
+ if(out_gradient_opt->answer)
+ if (G_legal_filename(out_gradient_opt->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal gradient map name"),
+ out_gradient_opt->answer);
+
+ if(out_curvature_opt->answer )
+ if (G_legal_filename(out_curvature_opt->answer ) < 0)
+ G_fatal_error(_("<%s> is an illegal gradient map name"),
+ out_curvature_opt->answer );
+
+if(!segmentation) {
+ MAP map_dirs, map_streams, map_elevation, map_output, map_identifier;
+ CELL **streams, **dirs, **identifier=NULL;
+ FCELL **elevation;
+ DCELL **output;
+
+ G_message(_("ALL IN RAM CALCULATION - DIRECTION: %s"),method_name[downstream]);
+ ram_create_map(&map_streams,CELL_TYPE);
+ ram_read_map(&map_streams,in_stm_opt->answer,1,CELL_TYPE);
+ ram_create_map(&map_dirs,CELL_TYPE);
+ ram_read_map(&map_dirs,in_dir_opt->answer,1,CELL_TYPE);
+ ram_create_map(&map_elevation,FCELL_TYPE);
+ ram_read_map(&map_elevation,in_elev_opt->answer,0,-1);
+
+ streams=(CELL**)map_streams.map;
+ dirs=(CELL**)map_dirs.map;
+ elevation=(FCELL**)map_elevation.map;
+
+ number_of_streams=ram_number_of_streams(streams, dirs)+1;
+ ram_build_streamlines(streams,dirs,elevation,number_of_streams);
+ ram_release_map(&map_streams);
+ ram_release_map(&map_dirs);
+ ram_create_map(&map_output,DCELL_TYPE);
+ output=(DCELL**)map_output.map; /* one output for all maps */
+
+ if(out_identifier_opt->answer) {
+ ram_create_map(&map_identifier,CELL_TYPE);
+ identifier=(CELL**)map_identifier.map;
+ ram_calculate_identifiers(identifier,number_of_streams,downstream);
+ ram_write_map(&map_identifier,out_identifier_opt->answer,CELL_TYPE,1,0);
+ ram_release_map(&map_identifier);
+ }
+
+ if(out_difference_opt->answer) {
+ ram_set_null_output(output);
+ if (local)
+ ram_calculate_difference(output,number_of_streams,downstream);
+ else
+ ram_calculate_drop(output,number_of_streams,downstream);
+ ram_write_map(&map_output,out_difference_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_distance_opt->answer) {
+ ram_set_null_output(output);
+ if(local && !cells)
+ ram_calculate_local_distance(output,number_of_streams,downstream);
+ if (!local && !cells)
+ ram_calculate_distance(output,number_of_streams,downstream);
+ if(cells)
+ ram_calculate_cell(output,number_of_streams,downstream);
+ ram_write_map(&map_output,out_distance_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_gradient_opt->answer) {
+ ram_set_null_output(output);
+ if(local)
+ ram_calculate_local_gradient(output,number_of_streams,downstream);
+ else
+ ram_calculate_gradient(output,number_of_streams,downstream);
+ ram_write_map(&map_output,out_gradient_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_curvature_opt->answer) {
+ ram_set_null_output(output);
+ ram_calculate_curvature(output,number_of_streams,downstream);
+ ram_write_map(&map_output,out_curvature_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ ram_release_map(&map_output);
+}
+
+
+if(segmentation) {
+ SEG map_dirs, map_streams, map_elevation, map_output, map_identifier;
+ SEGMENT *streams, *dirs, *elevation, *output, *identifier;
+
+ G_message(_("SEGMENT CALCULATION: MAY TAKE SOME TIME- DIRECTION: %s"),
+ method_name[downstream]);
+
+ number_of_segs = (int)atof(opt_swapsize->answer);
+ number_of_segs < 32 ? (int)(32/0.18) : number_of_segs/0.18;
+
+ seg_create_map(&map_streams,SROWS, SCOLS, number_of_segs, CELL_TYPE);
+ seg_read_map(&map_streams,in_stm_opt->answer,1,CELL_TYPE);
+ seg_create_map(&map_dirs,SROWS, SCOLS, number_of_segs,CELL_TYPE);
+ seg_read_map(&map_dirs,in_dir_opt->answer,1,CELL_TYPE);
+ seg_create_map(&map_elevation,SROWS, SCOLS, number_of_segs,FCELL_TYPE);
+ seg_read_map(&map_elevation,in_elev_opt->answer,0,-1);
+
+ streams=&map_streams.seg;
+ dirs=&map_dirs.seg;
+ elevation=&map_elevation.seg;
+
+ number_of_streams=seg_number_of_streams(streams, dirs)+1;
+ seg_build_streamlines(streams,dirs,elevation,number_of_streams);
+ seg_release_map(&map_streams);
+ seg_release_map(&map_dirs);
+ seg_create_map(&map_output,SROWS, SCOLS, number_of_segs, DCELL_TYPE);
+ output=&map_output.seg; /* one output for all maps */
+
+ if(out_identifier_opt->answer) {
+ seg_create_map(&map_identifier,SROWS, SCOLS, number_of_segs,CELL_TYPE);
+ identifier=&map_identifier.seg;
+ seg_calculate_identifiers(identifier,number_of_streams,downstream);
+ seg_write_map(&map_identifier,out_identifier_opt->answer,CELL_TYPE,1,0);
+ seg_release_map(&map_identifier);
+ }
+
+ if(out_difference_opt->answer) {
+ seg_set_null_output(output);
+ if (local)
+ seg_calculate_difference(output,number_of_streams,downstream);
+ else
+ seg_calculate_drop(output,number_of_streams,downstream);
+ seg_write_map(&map_output,out_difference_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_distance_opt->answer) {
+ seg_set_null_output(output);
+ if(local && !cells)
+ seg_calculate_local_distance(output,number_of_streams,downstream);
+ if (!local && !cells)
+ seg_calculate_distance(output,number_of_streams,downstream);
+ if(cells)
+ seg_calculate_cell(output,number_of_streams,downstream);
+ seg_write_map(&map_output,out_distance_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_gradient_opt->answer) {
+ seg_set_null_output(output);
+ if(local)
+ seg_calculate_local_gradient(output,number_of_streams,downstream);
+ else
+ seg_calculate_gradient(output,number_of_streams,downstream);
+ seg_write_map(&map_output,out_gradient_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ if(out_curvature_opt->answer) {
+ seg_set_null_output(output);
+ seg_calculate_curvature(output,number_of_streams,downstream);
+ seg_write_map(&map_output,out_curvature_opt->answer,DCELL_TYPE,0,0);
+ }
+
+ seg_release_map(&map_output);
+}
+ free_attributes(number_of_streams);
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/r.stream.channel.html
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/r.stream.channel.html (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/r.stream.channel.html 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,88 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-l</b></DT>
+<DD>Calculate local values. See output for detials.</DD>
+<DT><b>-c</b></DT>
+<DD>Calculate distance in cells instead of meters. See output for detials.</DD>
+<DT><b>-d</b></DT>
+<DD>Calculate downstream distance (from current cell DOWNSTREAM to outlet/join). Default is upstream (from current cell upstream to init/join.</DD>
+<DT><b>-m</b></DT>
+<DD>Only for very large data sets. Use segment library to optimize memory consumption during analysis</DD>
+<DT><b>stream</b></DT>
+<DD>Stream network: name of input stream map. Map may be ordered according one of the r.stream.order ordering system as well as unordered (with origilan stream identifiers) Because streams network produced by r.watershed and r.stream.extract may slighty differ in detail it is required to use both stream and direction map produced by the same module. Stream background shall have NULL value or zero value. Background values of NULL are by default produced by r.watershed and r.stream.extract. If not 0 or NULL use <a href="r.mapcalc.html">r.mapcalc</a> to set background values to null.
+</DD>
+<DT><b>dir</b></DT>
+<DD>Flow direction: name of input direction map produced by r.watershed or r.stream.extract. If r.stream.extract output map is used, it only has non-NULL values in places where streams occur. NULL (nodata) cells are ignored, zero and negative values are valid direction data if they vary from -8 to 8 (CCW from East in steps of 45 degrees). Direction map shall be of type CELL values. Region resolution and map resoultion must be the same.
+Also <em>stream</em> network map must have the same resolution. It is checked by default. If resolutions differ the module informs about it and stops. Region boundary and maps boundary may be differ but it may lead to unexpected results.</DD>
+<DT><b>elevation</b></DT>
+<DD>Elevation: name of input elevation map. Map can be of type CELL, FCELL or DCELL. It is not restricted to resolution of region settings as streams and dirs.</DD>
+</DL>
+<h2>OUTPUTS</h2>
+<DL>
+<DT><b>distance</b></DT>
+<DD>Upstream distance of current cell to the init/join. Flag modifications: <BR>
+<b>d:</b> downstream distance of current cell to the join/outlet;<BR>
+<b>l:</b> local distance between current cell and next cell. In most cases cell resolution and sqrt2 of cell resolution. usefull when projection is LL or NS and WE resolutions differs. Flag d ignored<BR>
+<b>c:</b> distance in cells. Map is written as double. Use mapcalc to convetrt to integer. Flags l and d ignored.<BR>
+</DD>
+<DT><b>difference</b></DT>
+<DD>Upstream elevation difference between current cell to the init/join. It we need to calculate parameters different than elevation. If we need to calculate different parameters than elevation along streams (for example precipitation or so) use neccesary map as elevation. Flag modifications: <BR>
+<b>d:</b> downstream difference of current cell to the join/outlet;<BR>
+<b>l:</b> local difference between current cell and next cell. With flag calculates difference between previous cell and current cell<BR>
+<b>c:</b> Ignored.
+</DD>
+<DT><b>gradient</b></DT>
+<DD>Upstream mean gradient between current cell and the init/join. Flag modifications: <BR>
+<b>d:</b> downstream mean gradient between current cell and the the join/outlet;<BR>
+<b>l:</b> local gradient between current cell and next cell. Flag d ignored<BR>
+<b>c:</b> Ignored.
+</DD>
+<DT><b>curvature</b></DT>
+<DD>Local stream course curvature of current cell. Calculated according formula: <i>first_derivative/(1-second_derivative<sup>2</sup>)<sup>3/2</sup></i> Flag modifications: <BR>
+<b>d:</b> ignored;<BR>
+<b>l:</b> Ignored.<BR>
+<b>c:</b> Ignored.
+</DD>
+<DT><b>identifier</b></DT>
+<DD> Integer map: In ordered stream network there are more than one segment (segment: a part of the network where order ramains unchanged) with the same order. To identify particular segments (for further analysis) every segment recive his unique identifier.</DD>
+</DL>
+<h2>DESCRIPTION</h2>
+<P>
+Module r.stream.channel is prepared to calculate some local properties of the stream network. It is suplementary module for r.stream.order, and r.stream.distance to investigate channel subsystem. For slope subsystem parameters is r.stream.slope. It may use ordered or unordered network. It calculate parameters for every segment between it init to outlet/join to the next stream. it also may calculate parameters between outlet and segment's init. It can calculate parameters for every orders but best results are for these orders where order remains unchanged from stream init to natural outlet (Hack and Horton ordering).
+<P>
+
+<h2>EGZAMPLE</h2>
+
+This example shows how to visualise the change of gradient map along the main stream of the catchment:
+<PRE>
+<CODE>
+r.watershed elevation=elevation.10m treshold=1000 stream=streams drainage=dirs
+r.stream.order streams=streams dirs=dirs hack=hack
+r.stream.channel streams=hack dirs=dirs elevation=elevation.10m identifier=stream_identifier gradient=stream_gradient distance=stream_distance
+#495 is a stream identifier. May be different in different situaltion
+r.mapcalc stgrad=if(stream_identifier==495,float(stream_gradient),null())
+r.mapcalc stdist=if(stream_identifier==495,float(stream_distance),null())
+
+#the rest data processing in R
+R
+library(spgrass6)
+r=readRAST6(c("stdist","stgrad"),NODATA=-9999)
+p=subset(r at data,!is.na(r at data$dist))
+sorted=p[order(p$stdist),]
+plot(sorted,stdist~stgrad,type="l")
+</CODE>
+</PRE>
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.stream.slope.html">r.stream.slope</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, Adam Mickiewicz University, Geoecology and Geoinformation Institute.
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/stream_topology.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/stream_topology.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/stream_topology.c 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,421 @@
+#include "local_proto.h"
+double get_distance(int r, int c, int d)
+{
+ double northing, easting, next_northing, next_easting;
+ int next_r, next_c;
+ next_r=NR(d);
+ next_c=NC(d);
+ northing = window.north - (r +.5) * window.ns_res;
+ easting = window.west + (c +.5) * window.ew_res;
+ next_northing = window.north - (next_r +.5) * window.ns_res;
+ next_easting = window.west + (next_c +.5) * window.ew_res;
+ return G_distance(easting,northing,next_easting,next_northing);
+}
+
+int ram_trib_nums(int r, int c, CELL** streams, CELL** dirs)
+{ /* calculate number of tributuaries */
+
+ int trib_num = 0;
+ int i, j;
+ int next_r, next_c;
+
+ for (i = 1; i < 9; ++i) {
+ if (NOT_IN_REGION(i))
+ continue;
+
+ j = DIAG(i);
+ next_r=NR(i);
+ next_c=NC(i);
+
+ if (streams[next_r][next_c] > 0 && dirs[next_r][next_c] == j)
+ trib_num++;
+ }
+
+ if (trib_num > 1)
+ for (i = 1; i < 9; ++i) {
+ if (NOT_IN_REGION(i))
+ continue;
+
+ j = DIAG(i);
+ next_r=NR(i);
+ next_c=NC(i);
+
+ if (streams[next_r][next_c] == streams[r][c] &&
+ dirs[next_r][next_c] == j)
+ trib_num--;
+ }
+
+ if (trib_num > 5)
+ G_fatal_error(_
+ ("Error finding inits. Stream and direction maps probably do not match..."));
+ if (trib_num > 3)
+ G_warning(_("Stream network may be too dense..."));
+
+ return trib_num;
+} /* end trib_num */
+
+
+
+int ram_number_of_streams(CELL** streams, CELL** dirs)
+{
+ int r, c;
+ int stream_num=0;
+
+ for (r = 0; r < nrows; ++r)
+ for (c = 0; c < ncols; ++c)
+ if (streams[r][c]>0)
+ if (ram_trib_nums(r, c, streams, dirs) != 1)
+ stream_num++;
+ return stream_num;
+}
+
+int ram_build_streamlines(CELL** streams, CELL** dirs, FCELL** elevation, int number_of_streams)
+{
+ int r, c, i;
+ int d, next_d;
+ int next_r, next_c;
+ int prev_r, prev_c;
+ int stream_num=1, cell_num=0;
+ int contrib_cell;
+ STREAM* SA;
+ int tmp=0;
+ stream_num = 1;
+ int v;
+ Rast_get_window(&window);
+
+ stream_attributes = (STREAM *) G_malloc(number_of_streams * sizeof(STREAM));
+ G_message("Finding inits...");
+ SA=stream_attributes;
+
+ for (r = 0; r < nrows; ++r)
+ for (c = 0; c < ncols; ++c)
+ if (streams[r][c])
+ if (ram_trib_nums(r, c, streams, dirs) != 1) { /* adding inits */
+ if(stream_num>number_of_streams)
+ G_fatal_error(_
+ ("Error finding inits. Stream and direction maps probably do not match..."));
+
+ SA[stream_num].stream_num = stream_num;
+ SA[stream_num].init_r = r;
+ SA[stream_num++].init_c = c;
+ }
+
+ for(i=1;i<stream_num;++i) {
+
+ r=SA[i].init_r;
+ c=SA[i].init_c;
+ SA[i].order=streams[r][c];
+ SA[i].number_of_cells=0;
+ do {
+
+ SA[i].number_of_cells++;
+ d=abs(dirs[r][c]);
+ if(NOT_IN_REGION(d) || d==0)
+ break;
+ r=NR(d);
+ c=NC(d);
+ } while (streams[r][c]==SA[i].order);
+
+ SA[i].number_of_cells+=2; /* add two extra points for init+ and outlet+ */
+ }
+
+ for(i=1;i<number_of_streams;++i) {
+
+ SA[i].points = (unsigned long int *)
+ G_malloc((SA[i].number_of_cells) * sizeof(unsigned long int));
+ SA[i].elevation = (float*)
+ G_malloc((SA[i].number_of_cells) * sizeof(float));
+ SA[i].distance = (double*)
+ G_malloc((SA[i].number_of_cells) * sizeof(double));
+
+ r=SA[i].init_r;
+ c=SA[i].init_c;
+ contrib_cell=ram_find_contributing_cell(r,c, dirs,elevation);
+ prev_r=NR(contrib_cell);
+ prev_c=NC(contrib_cell);
+
+ /* add one point contributing to init to calculate parameters */
+ /* what to do if there is no contributing points? */
+ SA[i].points[0] = (contrib_cell==0) ? -1 :
+ INDEX(prev_r,prev_c);
+ SA[i].elevation[0]=(contrib_cell==0) ? -99999 :
+ elevation[prev_r][prev_c];
+ d=(contrib_cell==0) ? dirs[r][c] : dirs[prev_r][prev_c];
+ SA[i].distance[0]=(contrib_cell==0) ? get_distance(r,c,d) :
+ get_distance(prev_r,prev_c,d);
+
+ SA[i].points[1] = INDEX(r,c);
+ SA[i].elevation[1]=elevation[r][c];
+ d=abs(dirs[r][c]);
+ SA[i].distance[1]=get_distance(r,c,d);
+
+ cell_num=2;
+ do {
+ d=abs(dirs[r][c]);
+
+ if(NOT_IN_REGION(d) || d==0) {
+ SA[i].points[cell_num] = -1;
+ SA[i].distance[cell_num] = SA[i].distance[cell_num-1];
+ SA[i].elevation[cell_num] =
+ 2*SA[i].elevation[cell_num-1]-SA[i].elevation[cell_num-2];
+ break;
+ }
+ r=NR(d);
+ c=NC(d);
+ SA[i].points[cell_num] = INDEX(r,c);
+ SA[i].elevation[cell_num] = elevation[r][c];
+ next_d=(abs(dirs[r][c])==0) ? d : abs(dirs[r][c]);
+ SA[i].distance[cell_num] = get_distance(r,c,next_d);
+ cell_num++;
+ if(cell_num>SA[i].number_of_cells)
+ G_fatal_error(_("To much points in stream line..."));
+ } while (streams[r][c]==SA[i].order);
+
+ if(SA[i].elevation[0]==-99999)
+ SA[i].elevation[0]=2*SA[i].elevation[1]-SA[i].elevation[2];
+ }
+
+ return 0;
+}
+
+int ram_find_contributing_cell(int r, int c, CELL** dirs, FCELL** elevation)
+{
+ int i,j=0;
+ int next_r, next_c;
+ float elev_min=9999;
+
+ for(i=1; i<9; ++i) {
+ if(NOT_IN_REGION(i))
+ continue;
+ next_r=NR(i);
+ next_c=NC(i);
+ if(dirs[next_r][next_c]==DIAG(i) && elevation[next_r][next_c]<elev_min) {
+ elev_min=elevation[next_r][next_c];
+ j=i;
+ }
+ }
+
+ return j;
+}
+
+int seg_trib_nums(int r, int c, SEGMENT* streams, SEGMENT* dirs)
+{ /* calculate number of tributuaries */
+
+ int trib_num = 0;
+ int i, j;
+ int next_r, next_c;
+ int streams_cell, streams_next_cell, dirs_next_cell;
+
+ segment_get(streams,&streams_cell,r,c);
+ for (i = 1; i < 9; ++i) {
+ if (NOT_IN_REGION(i))
+ continue;
+
+ j = DIAG(i);
+ next_r=NR(i);
+ next_c=NC(i);
+ segment_get(streams,&streams_next_cell,next_r,next_c);
+ segment_get(dirs,&dirs_next_cell,next_r,next_c);
+ if (streams_next_cell > 0 && dirs_next_cell == j)
+ trib_num++;
+ }
+
+ if (trib_num > 1)
+ for (i = 1; i < 9; ++i) {
+ if (NOT_IN_REGION(i))
+ continue;
+
+ j = DIAG(i);
+ next_r=NR(i);
+ next_c=NC(i);
+ segment_get(streams,&streams_next_cell,next_r,next_c);
+ segment_get(dirs,&dirs_next_cell,next_r,next_c);
+ if (streams_next_cell == streams_cell &&
+ dirs_next_cell == j)
+ trib_num--;
+ }
+
+ if (trib_num > 5)
+ G_fatal_error(_
+ ("Error finding inits. Stream and direction maps probably do not match..."));
+ if (trib_num > 3)
+ G_warning(_("Stream network may be too dense..."));
+
+ return trib_num;
+} /* end trib_num */
+
+int seg_number_of_streams(SEGMENT* streams, SEGMENT* dirs)
+{
+ int r, c;
+ int stream_num=0;
+ int streams_cell;
+ for (r = 0; r < nrows; ++r)
+ for (c = 0; c < ncols; ++c) {
+ segment_get(streams,&streams_cell,r,c);
+ if (streams_cell>0)
+ if (seg_trib_nums(r, c, streams, dirs) != 1)
+ stream_num++;
+ }
+ G_message("%d",stream_num);
+ return stream_num;
+}
+
+int seg_build_streamlines(SEGMENT* streams, SEGMENT* dirs, SEGMENT* elevation, int number_of_streams)
+{
+ int r, c, i;
+ int d, next_d;
+ int next_r, next_c;
+ int prev_r, prev_c;
+ int streams_cell, dirs_cell, streams_next_cell, dirs_next_cell;
+ float elevation_cell, elevation_prev_cell;
+ int stream_num=1, cell_num=0;
+ int contrib_cell;
+ STREAM* SA;
+ int tmp=0;
+ stream_num = 1;
+
+
+ stream_attributes = (STREAM *) G_malloc(number_of_streams * sizeof(STREAM));
+ G_message("Finding inits...");
+ SA=stream_attributes;
+
+ for (r = 0; r < nrows; ++r)
+ for (c = 0; c < ncols; ++c) {
+ segment_get(streams,&streams_cell,r,c);
+ if (streams_cell)
+ if (seg_trib_nums(r, c, streams, dirs) != 1) { /* adding inits */
+ if(stream_num>number_of_streams)
+ G_fatal_error(_
+ ("Error finding inits. Stream and direction maps probably do not match..."));
+
+ SA[stream_num].stream_num = stream_num;
+ SA[stream_num].init_r = r;
+ SA[stream_num++].init_c = c;
+ }
+ }
+
+ for(i=1;i<stream_num;++i) {
+
+ r=SA[i].init_r;
+ c=SA[i].init_c;
+ segment_get(streams,&(SA[i].order),r,c);
+ //segment_get(streams,&streams_cell,r,c);
+ SA[i].number_of_cells=0;
+ do {
+
+ SA[i].number_of_cells++;
+ segment_get(dirs,&dirs_cell,r,c);
+ d=abs(dirs_cell);
+ if(NOT_IN_REGION(d) || d==0)
+ break;
+ r=NR(d);
+ c=NC(d);
+ segment_get(streams,&streams_cell,r,c);
+ } while (streams_cell==SA[i].order);
+
+ SA[i].number_of_cells+=2; /* add two extra points for init+ and outlet+ */
+ }
+
+ for(i=1;i<number_of_streams;++i) {
+
+ SA[i].points = (unsigned long int *)
+ G_malloc((SA[i].number_of_cells) * sizeof(unsigned long int));
+ SA[i].elevation = (float*)
+ G_malloc((SA[i].number_of_cells) * sizeof(float));
+ SA[i].distance = (double*)
+ G_malloc((SA[i].number_of_cells) * sizeof(double));
+
+ r=SA[i].init_r;
+ c=SA[i].init_c;
+ contrib_cell=seg_find_contributing_cell(r,c, streams, dirs,elevation);
+ prev_r=NR(contrib_cell);
+ prev_c=NC(contrib_cell);
+
+ /* add one point contributing to init to calculate parameters */
+ /* what to do if there is no contributing points? */
+ SA[i].points[0] = (contrib_cell==0) ? -1 :
+ INDEX(prev_r,prev_c);
+
+ segment_get(elevation,&elevation_prev_cell,prev_r,prev_c);
+ SA[i].elevation[0]=(contrib_cell==0) ? -99999 :
+ elevation_prev_cell;
+
+ if(contrib_cell==0)
+ segment_get(dirs,&d,r,c);
+ else
+ segment_get(dirs,&d,prev_r,prev_c);
+ SA[i].distance[0]=(contrib_cell==0) ? get_distance(r,c,d) :
+ get_distance(prev_r,prev_c,d);
+
+ SA[i].points[1] = INDEX(r,c);
+ segment_get(elevation,&(SA[i].elevation[1]),r,c);
+ segment_get(dirs,&d,r,c);
+ SA[i].distance[1]=get_distance(r,c,d);
+
+ cell_num=2;
+ do {
+ segment_get(dirs,&dirs_cell,r,c);
+ d=abs(dirs_cell);
+ if(NOT_IN_REGION(d) || d==0) {
+ SA[i].points[cell_num] = -1;
+ SA[i].distance[cell_num] = SA[i].distance[cell_num-1];
+ SA[i].elevation[cell_num] =
+ 2*SA[i].elevation[cell_num-1]-SA[i].elevation[cell_num-2];
+ break;
+ }
+ r=NR(d);
+ c=NC(d);
+ SA[i].points[cell_num] = INDEX(r,c);
+ segment_get(elevation,&(SA[i].elevation[cell_num]),r,c);
+ segment_get(dirs,&next_d,r,c);
+ next_d=(abs(next_d)==0) ? d : abs(next_d);
+ SA[i].distance[cell_num] = get_distance(r,c,next_d);
+ cell_num++;
+ segment_get(streams,&streams_cell,r,c);
+ if(cell_num>SA[i].number_of_cells)
+ G_fatal_error(_("To much points in stream line..."));
+ } while (streams_cell==SA[i].order);
+
+ if(SA[i].elevation[0]==-99999)
+ SA[i].elevation[0]=2*SA[i].elevation[1]-SA[i].elevation[2];
+ }
+ return 0;
+}
+
+int seg_find_contributing_cell(int r, int c, SEGMENT* dirs, SEGMENT* elevation)
+{
+ int i,j=0;
+ int next_r, next_c;
+ float elev_min=9999;
+ int dirs_next_cell;
+ float elevation_next_cell;
+
+ for(i=1; i<9; ++i) {
+ if(NOT_IN_REGION(i))
+ continue;
+ next_r=NR(i);
+ next_c=NC(i);
+ segment_get(dirs,&dirs_next_cell,next_r,next_c);
+ segment_get(elevation,&elevation_next_cell,next_r,next_c);
+ if(dirs_next_cell==DIAG(i) && elevation_next_cell<elev_min) {
+ elev_min=elevation_next_cell;
+ j=i;
+ }
+ }
+ return j;
+}
+
+int free_attributes (int number_of_streams)
+{
+ int i;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ G_free(SA[i].points);
+ G_free(SA[i].elevation);
+ G_free(SA[i].distance);
+ }
+ G_free(stream_attributes);
+ return 0;
+}
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/stream_write.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/stream_write.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/stream_write.c 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,453 @@
+#include "local_proto.h"
+int ram_set_null_output(DCELL** output) {
+
+ int r,c;
+ for (r=0;r<nrows;++r)
+ Rast_set_d_null_value(output[r],ncols);
+ return 0;
+}
+
+int seg_set_null_output(SEGMENT* output) {
+
+ int r,c;
+ double output_cell;
+ for (r=0;r<nrows;++r)
+ for (c=0;c<ncols;++c) {
+ Rast_set_d_null_value(&output_cell,1);
+ segment_put(output,&output_cell,r,c);
+ }
+ return 0;
+}
+
+int ram_calculate_identifiers(CELL** identifier,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ identifier[r][c] = SA[i].stream_num;
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_identifiers(SEGMENT* identifier,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ segment_put(identifier,&(SA[i].stream_num),r,c);
+ }
+ }
+ return 0;
+}
+
+
+int ram_calculate_distance(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ double cum_length;
+ int i,j,k;
+ STREAM* SA;
+
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ cum_length=0;
+ if(!downstream)
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c]=cum_length;
+ }
+ else
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c]=cum_length;
+ }
+ }
+ return 0;
+}
+int seg_calculate_distance(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ double cum_length;
+ int i,j,k;
+ STREAM* SA;
+
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ cum_length=0;
+ if(!downstream)
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ segment_put(output,&cum_length,r,c);
+ }
+ else
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ segment_put(output,&cum_length,r,c);
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_cell(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j,k;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+
+ k=SA[i].number_of_cells-2;
+ for(j=1;j<SA[i].number_of_cells-1;++j,--k) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = downstream ? k : j;
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_cell(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j,k;
+ double output_cell;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+
+ k=SA[i].number_of_cells-2;
+ for(j=1;j<SA[i].number_of_cells-1;++j,--k) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = downstream ? k : j;
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_difference(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double result;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ result= downstream ?
+ SA[i].elevation[j-1]-SA[i].elevation[j] :
+ SA[i].elevation[j]-SA[i].elevation[j+1];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = result;
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_difference(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double output_cell;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ output_cell= downstream ?
+ SA[i].elevation[j-1]-SA[i].elevation[j] :
+ SA[i].elevation[j]-SA[i].elevation[j+1];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_drop(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double init;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ if(!downstream) {
+ init = SA[i].elevation[1];
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = init-SA[i].elevation[j];
+ }
+ } else {
+ init = SA[i].elevation[SA[i].number_of_cells-2];
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = SA[i].elevation[j]-init;
+ }
+ }
+ }
+ return 0;
+}
+int seg_calculate_drop(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double init;
+ double output_cell;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ if(!downstream) {
+ init = SA[i].elevation[1];
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = init-SA[i].elevation[j];
+ segment_put(output,&output_cell,r,c);
+ }
+ } else {
+ init = SA[i].elevation[SA[i].number_of_cells-2];
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = SA[i].elevation[j]-init;
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_gradient(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double init;
+ double result;
+ double cum_length;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ cum_length=0;
+ if(!downstream) {
+ init = SA[i].elevation[0];
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = (init-SA[i].elevation[j])/cum_length;
+ }
+ } else {
+ init = SA[i].elevation[SA[i].number_of_cells-1];
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = (SA[i].elevation[j]-init)/cum_length;
+ }
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_gradient(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double init;
+ double output_cell;
+ double cum_length;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ cum_length=0;
+ if(!downstream) {
+ init = SA[i].elevation[1];
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = (init-SA[i].elevation[j])/cum_length;
+ segment_put(output,&output_cell,r,c);
+
+ }
+ } else {
+ init = SA[i].elevation[SA[i].number_of_cells-2];
+ for(j=SA[i].number_of_cells-2;j>0;--j) {
+ cum_length+=SA[i].distance[j];
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = (SA[i].elevation[j]-init)/cum_length;
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_local_gradient(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double elev_diff;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ elev_diff=(SA[i].elevation[j]-SA[i].elevation[j+1])<0 ? 0 :
+ (SA[i].elevation[j]-SA[i].elevation[j+1]);
+ output[r][c] = elev_diff/SA[i].distance[j];
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_local_gradient(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double elev_diff;
+ double output_cell;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ elev_diff=(SA[i].elevation[j]-SA[i].elevation[j+1])<0 ? 0 :
+ (SA[i].elevation[j]-SA[i].elevation[j+1]);
+ output_cell = elev_diff/SA[i].distance[j];
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ return 0;
+}
+
+
+int ram_calculate_local_distance(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output[r][c] = SA[i].distance[j];
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_local_distance(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double output_cell;
+ STREAM* SA;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ output_cell = SA[i].distance[j];
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ return 0;
+}
+
+int ram_calculate_curvature(DCELL** output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ STREAM* SA;
+ double result, first_derivative, second_derivative;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ first_derivative=(SA[i].elevation[j-1]-SA[i].elevation[j+1])/
+ (SA[i].distance[j-1]+SA[i].distance[j]);
+ second_derivative=((SA[i].elevation[j-1]-SA[i].elevation[j])-
+ (SA[i].elevation[j]-SA[i].elevation[j+1]))/
+ (SA[i].distance[j-1]+SA[i].distance[j]);
+ output[r][c] =
+ first_derivative/pow((1+second_derivative*second_derivative),1.5);
+ }
+ }
+ return 0;
+}
+
+int seg_calculate_curvature(SEGMENT* output,int number_of_streams, int downstream)
+{
+ int r,c;
+ int i,j;
+ double output_cell;
+ STREAM* SA;
+ double result, first_derivative, second_derivative;
+ SA=stream_attributes;
+
+ for(i=1;i<number_of_streams;++i) {
+ for(j=1;j<SA[i].number_of_cells-1;++j) {
+ r=(int)SA[i].points[j]/ncols;
+ c=(int)SA[i].points[j]%ncols;
+ first_derivative=(SA[i].elevation[j-1]-SA[i].elevation[j+1])/
+ (SA[i].distance[j-1]+SA[i].distance[j]);
+ second_derivative=((SA[i].elevation[j-1]-SA[i].elevation[j])-
+ (SA[i].elevation[j]-SA[i].elevation[j+1]))/
+ (SA[i].distance[j-1]+SA[i].distance[j]);
+ output_cell =
+ first_derivative/pow((1+second_derivative*second_derivative),1.5);
+ segment_put(output,&output_cell,r,c);
+ }
+ }
+ return 0;
+}
+
Added: grass-addons/grass7/raster/r.stream/r.stream.channel/svn-commit.tmp
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.channel/svn-commit.tmp (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.channel/svn-commit.tmp 2010-10-01 11:20:30 UTC (rev 43754)
@@ -0,0 +1,4 @@
+files
+--Ta linia i następne zostaną zignorowane--
+
+A r.stream
More information about the grass-commit
mailing list