[GRASS-SVN] r45419 - in grass-addons/grass7/raster/r.stream: .
r.stream.snap
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 17 14:19:07 EST 2011
Author: jarekj71
Date: 2011-02-17 11:19:07 -0800 (Thu, 17 Feb 2011)
New Revision: 45419
Added:
grass-addons/grass7/raster/r.stream/r.stream.snap/
grass-addons/grass7/raster/r.stream/r.stream.snap/Makefile
grass-addons/grass7/raster/r.stream/r.stream.snap/io.c
grass-addons/grass7/raster/r.stream/r.stream.snap/io.h
grass-addons/grass7/raster/r.stream/r.stream.snap/local_proto.h
grass-addons/grass7/raster/r.stream/r.stream.snap/local_vars.h
grass-addons/grass7/raster/r.stream/r.stream.snap/main.c
grass-addons/grass7/raster/r.stream/r.stream.snap/points_io.c
grass-addons/grass7/raster/r.stream/r.stream.snap/r.stream.snap.html
grass-addons/grass7/raster/r.stream/r.stream.snap/snap.c
Log:
Initial import
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/Makefile 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.stream.snap
+
+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.snap/io.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/io.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/io.c 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,523 @@
+#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) {
+ /*
+ * 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 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.snap/io.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/io.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/io.h 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,54 @@
+#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>
+
+#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_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.snap/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/local_proto.h 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,3 @@
+#include "io.h"
+#include "local_vars.h"
+
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/local_vars.h
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/local_vars.h (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/local_vars.h 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,40 @@
+#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
+
+#define SQR(x) ((x) * (x))
+
+typedef struct {
+ int r,c;
+ int di,dj; /* shift */
+ int cat;
+ double accum;
+ int stream;
+ int status; /* 1=skipped,2=unresolved,3=snapped,4=correct */
+} OUTLET;
+
+GLOBAL int nextr[9];
+GLOBAL int nextc[9];
+
+GLOBAL OUTLET* points;
+GLOBAL int nrows, ncols;
+GLOBAL float** distance_mask;
+
+
+
+
+
+
+
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/main.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/main.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/main.c 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,146 @@
+/****************************************************************************
+ *
+ * MODULE: r.stream.snap
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ *
+ * PURPOSE: Snap ponints features to nearest pour points. Usefull both to
+ * snap outlets well as sources. Use two parameters:
+ * maximum snap distance and minimum accumulation value to snap
+ *
+ *
+ * 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[])
+{
+
+ struct GModule *module;
+ struct Option *in_points_opt,
+ *in_points_cat_opt,
+ *out_points_opt,
+ *in_stream_opt,
+ *in_accum_opt,
+ *opt_accum_treshold,
+ *opt_distance_treshold,
+ *opt_swapsize;
+
+
+ int snap_to_streams;
+ int i;
+ SEG map_dirs, map_streams, map_accum;
+ SEGMENT *streams=NULL, *accum=NULL;
+ int number_of_segs;
+ int number_of_points;
+ int radius;
+ double accum_treshold;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->description = _("Delineate basins according user' input. \
+ Input can be stream network, point file with outlets or outlet coordinates");
+ G_add_keyword("basins creation");
+
+ in_points_opt = G_define_standard_option(G_OPT_V_INPUT);
+ in_points_opt->description = _("Name of input vector points map");
+
+ out_points_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ out_points_opt->description = _("Name of output vector points map");
+
+ in_stream_opt = G_define_standard_option(G_OPT_R_INPUT);
+ in_stream_opt->key = "streams";
+ in_stream_opt->required = NO;
+ in_stream_opt->description = _("Name of stream map");
+
+ in_accum_opt = G_define_standard_option(G_OPT_R_INPUT);
+ in_accum_opt->key = "accum";
+ in_accum_opt->required = NO;
+ in_accum_opt->description = _("Name of accumulation map");
+
+ opt_accum_treshold = G_define_option();
+ opt_accum_treshold->key="accumtres";
+ opt_accum_treshold->type = TYPE_DOUBLE;
+ opt_accum_treshold->answer="-1";
+ opt_accum_treshold->description =_("Minimum accumulation streshold to snap");
+
+ opt_distance_treshold = G_define_option();
+ opt_distance_treshold->key="radius";
+ opt_distance_treshold->answer="1";
+ opt_distance_treshold->type = TYPE_INTEGER;
+ opt_distance_treshold->description =_("Maximum distance to snap (in cells)");
+
+ opt_swapsize = G_define_option();
+ opt_swapsize->key="memory";
+ opt_swapsize->type = TYPE_INTEGER;
+ opt_swapsize->answer = "300";
+ opt_swapsize->required = NO;
+ opt_swapsize->description =_("Max memory used (MB)");
+
+ if (G_parser(argc, argv)) /* parser */
+ exit(EXIT_FAILURE);
+
+ number_of_segs = (int)atof(opt_swapsize->answer);
+ number_of_segs < 32 ? (int)(32/0.12) : number_of_segs/0.12;
+
+ radius=atoi(opt_distance_treshold->answer);
+ accum_treshold=atof(opt_accum_treshold->answer);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ if (G_legal_filename(out_points_opt->answer) < 0)
+ G_fatal_error(_("<%s> is an illegal basin name"), out_points_opt->answer);
+
+ if(!in_stream_opt->answer && !in_accum_opt->answer)
+ G_fatal_error(_("At least one map of accumulation or streams is required"));
+
+ if(!in_accum_opt->answer)
+ accum_treshold=-1;
+
+
+ /* SEGMENT VERSION ONLY*/
+
+ if(in_stream_opt->answer) {
+ seg_create_map(&map_streams,SROWS, SCOLS, number_of_segs, CELL_TYPE);
+ seg_read_map(&map_streams,in_stream_opt->answer,1,CELL_TYPE);
+ streams = &map_streams.seg;
+ }
+
+ if(in_accum_opt->answer) {
+ seg_create_map(&map_accum,SROWS, SCOLS, number_of_segs, DCELL_TYPE);
+ seg_read_map(&map_accum,in_accum_opt->answer,0,-1);
+ accum = &map_accum.seg;
+ }
+
+ create_distance_mask(radius);
+ number_of_points=read_points(in_points_opt->answer, streams, accum);
+
+ for(i=0;i<number_of_points;++i)
+ snap_point(&points[i],radius,streams,accum,accum_treshold);
+
+ write_points(out_points_opt->answer,number_of_points);
+
+
+ for(i=0;i<number_of_points;++i)
+ G_message("AFTER %d %d %d %d",
+ points[i].r,points[i].c,points[i].di,points[i].dj);
+
+
+ if(in_stream_opt->answer)
+ seg_release_map (&map_streams);
+ if(in_accum_opt->answer)
+ seg_release_map (&map_accum);
+
+exit(EXIT_SUCCESS);
+ }
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/points_io.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/points_io.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/points_io.c 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,117 @@
+#include "local_proto.h"
+
+int read_points(char *in_point, SEGMENT* streams, SEGMENT* accum)
+{
+ char *mapset;
+ struct Cell_head window;
+ struct Map_info Map;
+ struct bound_box box;
+ int num_point = 0;
+ int total_points = 0;
+ int type, i,j, cat;
+ int r,c;
+ struct line_pnts *sites;
+ struct line_cats *cats;
+ double absaccum;
+ sites = Vect_new_line_struct();
+ cats = Vect_new_cats_struct();
+
+ mapset = G_find_vector2(in_point, "");
+ if (mapset == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), in_point);
+
+ if (Vect_open_old(&Map, in_point, mapset) < 0)
+ G_fatal_error("Cannot open vector map <%s>", in_point);
+
+ G_get_window(&window);
+ Vect_region_box(&window, &box);
+
+ while ((type = Vect_read_next_line(&Map, sites, cats)) > -1) {
+ if (type != GV_POINT)
+ continue;
+ if (Vect_point_in_box(sites->x[0], sites->y[0], sites->z[0], &box))
+ num_point++;
+ }
+
+ points = (OUTLET *) G_malloc(num_point * sizeof(OUTLET));
+ total_points=Vect_get_num_lines(&Map);
+ i=0;
+
+ for (j = 0; j < total_points; ++j) {
+
+ type = Vect_read_line(&Map, sites, cats,j+1);
+
+ if (type != GV_POINT)
+ continue;
+
+ if (!Vect_point_in_box(sites->x[0], sites->y[0], sites->z[0], &box))
+ continue;
+
+ Vect_cat_get(cats, 1, &cat);
+
+ points[i].r = (int)Rast_northing_to_row(sites->y[0], &window);
+ points[i].c = (int)Rast_easting_to_col(sites->x[0], &window);
+ points[i].di=0;
+ points[i].dj=0;
+ points[i].cat = cat;
+ if(streams)
+ segment_get(streams,&points[i].stream,points[i].r,points[i].c);
+ else
+ points[i].stream=0;
+ if(accum) {
+ segment_get(accum,&absaccum,points[i].r,points[i].c);
+ points[i].accum=fabs(absaccum);
+ } else {
+ points[i].accum=0;
+ points[i].status = 4; /* default status is 'correct' */
+ }
+ //dodać skip category
+
+ i++;
+ }
+ return num_point;
+}
+
+int write_points(char* out_vector, int number_of_points) {
+
+ int i;
+ int r,c;
+ int cat_layer_1, cat_layer_2;
+ float northing, easting;
+ struct Cell_head window;
+ struct Map_info Out;
+ struct line_pnts *Segments;
+ struct line_cats *Cats;
+
+ 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_points;++i) {
+
+ r = points[i].r+points[i].di;
+ c = points[i].c+points[i].dj;
+
+ cat_layer_1=points[i].cat;
+ cat_layer_2=points[i].status;
+ Vect_cat_set(Cats,1,cat_layer_1);
+ Vect_cat_set(Cats,2,cat_layer_2);
+ 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_reset_cats(Cats);
+ }
+
+/* build vector */
+Vect_hist_command(&Out);
+Vect_build(&Out);
+Vect_close(&Out);
+
+return 0;
+}
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/r.stream.snap.html
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/r.stream.snap.html (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/r.stream.snap.html 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,60 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>streams</b></DT>
+<DD>Stream network created by r.stream.extract or r.watershed. If used, points are snapped to the nearest streamline point which accumulation is grater than treshold. If accumulation is not used point is snapped to the nearest stream. </DD>
+
+<DT><b>accum</b></DT>
+<DD>Accumulation map created with r.waterhed and used to genearte stream network with r.stream.extract. If stream network is not in use, point is adaptively snapped to the to point where value is greater than mean values of accumulation grater than given treshold in a searcyh radius. See description for details.
+</DD>
+<DT><b>radius</b></DT>
+<DD>Search radius (in cells). If there are no streams in search radius, point is not snapped. If there are no cells with accumulation grater than accumtreshold point also is not snapped.
+</DD>
+
+<DT><b>accumtres</b></DT>
+<DD>Minimum value of accumulation which cell must have to snap point. This option is added to snap stream inits to the stream tubes and to distinguish between local tributuaries and main streams.
+</DD>
+
+<DT><b>input</b></DT>
+<DD>Vector file containing outlets or inits as vector points. Only point's categories are used Table attached to it is ignored. Every point shall heve his own unique category.
+</DD>
+</DL>
+
+
+<h2>OUTPUTS</h2>
+<P>Vector file containing outlets or inits after snapping. On layer 1 oryginal categories are preserved, on layer 2 there are four categories which mean:</p>
+<ol>
+<li>skipped (not in use yet)
+<li>unresolved (points remain unsnapped due to lack of streams in search radius
+<li>snapped (points snapped to streamlines)
+<li>correct (points whcich remain on its original position, wchich was originally corected
+</ol>
+
+
+
+<h2>DESCRIPTION</h2>
+<P>
+Module r.stream.snap is a suplementary module for r.stream.extract and r.stream.basins to correct position of outlets or stream init points as they do not lie on the streamliens.
+<BR>
+For outlets situation is clear. Points are snapped to nearest point whcich lies on the streamline. In situaltion where there can be a small tributuary nearer than main stream accumulation treshold shall be high enough to force program ignoring this tributuary and snap to the main stream. If there is no accumulation map points will be snapped to nearest stream line whcich in particular situation may be wrong. Bacuase r.stream is perepared to work with MFD accum maps, both stream network and accum map is neccesary to resolve the problem.
+<BR>
+While it is assumed accum map is a MFD map, if stream network is not supplied, snap point is calculated in different way: treshold is used to select only these points in search radius wchich are accum value is greater than treshold. Next mean value of these points is calculated and its value is taken as a new treshold. This procedure guarantee that points are snapped to the center of stream tube. While for inits small tresholds are in use, it is probable than points were snapped to the streamtube border instead of its center.
+</p>
+<p>
+It is strongly recommended, to use both stream network (even pre-genearted with small accum treshold) and accumulation map, than accum or stream map only.
+</p>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.watershed.html">r.watershed</a>,
+<a href="r.stream.extract.html">r.stream.extract</a>,
+<a href="r.stream.basins.html">r.stream.basins</a>,
+<a href="r.stream.order.html">r.stream.order</a>,
+<a href="r.stream.stats.html">r.stream.stats</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.reclass.html">r.reclass</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz, Adam Mickiewicz University, Geoecology and Geoinformation Institute.
Added: grass-addons/grass7/raster/r.stream/r.stream.snap/snap.c
===================================================================
--- grass-addons/grass7/raster/r.stream/r.stream.snap/snap.c (rev 0)
+++ grass-addons/grass7/raster/r.stream/r.stream.snap/snap.c 2011-02-17 19:19:07 UTC (rev 45419)
@@ -0,0 +1,140 @@
+#include "local_proto.h"
+#define SQRT(x) ((x) * (x))
+static int num;
+
+int create_distance_mask(int radius) {
+
+ int i, j;
+ int window=2*radius+1;
+
+ distance_mask = G_malloc(window * sizeof(float *));
+
+ for (i = 0; i < window; i++)
+ distance_mask[i] = G_malloc(window *sizeof(float));
+
+ for (i = 0; i < window; i++)
+ for (j = 0; j < window; j++)
+ distance_mask[i][j] =
+ (SQR(i - radius) + SQR(j - radius)<= SQR(radius)) ?
+ sqrt(SQR(i - radius) + SQR(j - radius)) : 0;
+
+ return 0;
+}
+
+int snap_point(OUTLET *point, int radius, SEGMENT* streams, SEGMENT* accum, double accum_treshold) {
+
+ int i,j,di=-1,dj=-1;
+ int status=3;
+ int number_of_cells=radius*radius;
+ int teststream=0;
+ float cur_distance=radius;
+ float distance=0;
+ double absaccum=0;
+ double sumaccum=0;
+ double maxaccum=0;
+ int naccum=-1;
+
+ if(point->stream>0 && point->accum>accum_treshold)
+ return 0; /* point lies on line(or skipped) and has proper treshold */
+
+if(streams) {
+ /* stream version: assume ve have stream network and points
+ * are snapped to stream points where accum is greater than treshold
+ * or to nearest stream point in accum is not supplied */
+
+ for(i=-radius; i<=radius; ++i)
+ for(j=-radius; j<=radius; ++j) {
+
+ if(point->r+i<0 || point->r+i>=nrows ||
+ point->c+j<0 || point->c+j>=ncols)
+ continue;
+
+ if(!distance_mask[i+radius][j+radius])
+ continue;
+
+ segment_get(streams,&teststream,point->r+i,point->c+j);
+ distance=distance_mask[i+radius][j+radius];
+
+ if(teststream) { /* is stream line */
+
+ if(accum) {
+ segment_get(accum,&absaccum,point->r+i,point->c+j);
+ absaccum=fabs(absaccum);
+ }
+
+ if(absaccum>=accum_treshold) /* if no accum absaccum always =-1 */
+ if(cur_distance>distance) {
+ cur_distance=distance;
+ di=i;
+ dj=j;
+ }
+ }
+ }
+} /* end of streams version */
+
+if(!streams) {
+ /* no stream version: problem for MFD. the snap point is found
+ * in different manner. It is not only point where accum exceed the
+ * treshold (may be far for potential streamline) but must exceed the
+ * mean value of accums in searach area taken in cells where treshold is exceeded */
+
+ for(i=-radius; i<=radius; ++i)
+ for(j=-radius; j<=radius; ++j) {
+
+ if(point->r+i<0 || point->r+i>=nrows ||
+ point->c+j<0 || point->c+j>=ncols)
+ continue;
+
+ if(!distance_mask[i+radius][j+radius])
+ continue;
+
+ segment_get(accum,&absaccum,point->r+i,point->c+j);
+ absaccum=fabs(absaccum);
+
+ if(absaccum>maxaccum)
+ maxaccum=absaccum;
+
+ if(absaccum>accum_treshold) {
+ sumaccum+=absaccum;
+ naccum++;
+ }
+ }
+
+ if(sumaccum>0)
+ /* accum_treshold=(sumaccum/naccum+maxaccum)/2 */;
+ accum_treshold=sumaccum/naccum;
+
+ for(i=-radius; i<=radius; ++i)
+ for(j=-radius; j<=radius; ++j) {
+
+ if(point->r+i<0 || point->r+i>=nrows || point->c+j<0 || point->c+j>=ncols)
+ continue;
+
+ if(!distance_mask[i+radius][j+radius])
+ continue;
+
+ segment_get(accum,&absaccum,point->r+i,point->c+j);
+ absaccum=fabs(absaccum);
+
+ if(accum_treshold >0 && absaccum>accum_treshold)
+ if(cur_distance>distance) {
+ cur_distance=distance;
+ di=i;
+ dj=j;
+ }
+ }
+}/* end of non-streams version */
+ if(di==-1&&dj==-1){
+ G_warning(_("cannot snap point with cat %d, in a given radius. \
+ Increase search radius"),point->cat);
+ di=0;
+ dj=0;
+ status=2;
+ }
+ point->di=di;
+ point->dj=dj;
+ point->status=status;
+ return 0;
+}
+
+
More information about the grass-commit
mailing list