[GRASS-SVN] r40250 - grass-addons/raster/r.stream.order
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 5 08:45:16 EST 2010
Author: jarekj71
Date: 2010-01-05 08:45:10 -0500 (Tue, 05 Jan 2010)
New Revision: 40250
Modified:
grass-addons/raster/r.stream.order/Makefile
grass-addons/raster/r.stream.order/global.h
grass-addons/raster/r.stream.order/io.c
grass-addons/raster/r.stream.order/main.c
grass-addons/raster/r.stream.order/order.c
Log:
update
Modified: grass-addons/raster/r.stream.order/Makefile
===================================================================
--- grass-addons/raster/r.stream.order/Makefile 2010-01-05 13:28:07 UTC (rev 40249)
+++ grass-addons/raster/r.stream.order/Makefile 2010-01-05 13:45:10 UTC (rev 40250)
@@ -4,8 +4,10 @@
PGM = r.stream.order
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/raster/r.stream.order/global.h
===================================================================
--- grass-addons/raster/r.stream.order/global.h 2010-01-05 13:28:07 UTC (rev 40249)
+++ grass-addons/raster/r.stream.order/global.h 2010-01-05 13:45:10 UTC (rev 40250)
@@ -3,6 +3,8 @@
#include <string.h>
#include <math.h>
#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/dbmi.h>
/* define */
#define NODES struct node
@@ -40,10 +42,13 @@
int strahler, shreeve, horton, hack; /* orders */
int next_stream;
int trib_num;
- int trib[5]; /* for now, I have problems with realock */
-
+ int trib[5];
+
/* int *trib; */
double accum;
+ double length;
+ double stright;
+ double fractal;
/* float length; */
};
@@ -73,12 +78,13 @@
int destroy_c_data(CELL ** data_name);
int destroy_d_data(DCELL ** data_name);
int write_maps(void);
+int create_table (void);
/* variables */
GLOBAL struct Cell_head window;
-GLOBAL char *in_dirs, *in_streams, *in_accum; /* input dirrection and accumulation raster names */
+GLOBAL char *in_dirs, *in_streams, *in_vector, *in_table, *in_accum; /* input dirrection and accumulation raster names */
GLOBAL char *out_strahler, *out_shreeve, *out_hack, *out_horton; /* output strahler, shreeve and class raster names */
GLOBAL int out_zero; /* zero-value background */
@@ -91,4 +97,4 @@
GLOBAL STREAM *s_streams; /* stream structure all parameters we have here */
GLOBAL INITS *s_inits;
int *springs, *outlets;
-int springs_num, outlets_num;
+int springs_num, outlets_num, stream_num;
Modified: grass-addons/raster/r.stream.order/io.c
===================================================================
--- grass-addons/raster/r.stream.order/io.c 2010-01-05 13:28:07 UTC (rev 40249)
+++ grass-addons/raster/r.stream.order/io.c 2010-01-05 13:45:10 UTC (rev 40250)
@@ -125,18 +125,7 @@
G_read_range(in_streams,cur_mapset,&stream_range);
G_get_range_min_max(&stream_range, &c_min, &c_max);
return c_max;
-/*
-int r, c;
-int max=0;
- for (r = 0; r < nrows; ++r) {
- for (c = 0; c < ncols; ++c) {
- if (streams[r][c]>max)
- max=streams[r][c];
- }
- }
- return max;
-*/
}
int write_maps(void)
@@ -257,3 +246,195 @@
return 0;
} /* end write_maps */
+
+int create_table (void)
+{
+
+ int i,j;
+ int index_cat=0;
+ int max_trib=0;
+ /*
+ char *mapset;
+ struct Map_info Map;
+ char out_table[30];
+ */
+ char* out_table;
+ dbConnection conn;
+ dbDriver *driver;
+ dbHandle handle;
+ dbString table_name, db_sql, val_string;
+ char buf[1000];
+ char ins_prev_streams[50]; /* insert */
+
+ /* table definition */
+ char* tab_cat_col_name="cat integer";
+ char* tab_stream="stream integer";
+ char* tab_next_stream="next_stream integer";
+ char* tab_prev_streams;
+ char* tab_strahler="strahler integer";
+ char* tab_horton="horton integer";
+ char* tab_shreve="shreve integer";
+ char* tab_hack="hack integer";
+ char* tab_length="length double precision";
+ char* tab_cumlength="cum_length double precision";
+ char* tab_stright="stright double precision";
+ char* tab_fractal="fractal double precision";
+
+ G_message("Adding table...");
+
+ /* trib num */
+ for (i=0;i<stream_num;++i) {
+ if (s_streams[i].trib_num>max_trib)
+ max_trib=s_streams[i].trib_num;
+ }
+
+ /*
+ mapset = G_find_vector(in_vector, "");
+ if (mapset == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), in_vector);
+
+ if (Vect_open_update(&Map, in_vector, mapset) < 0)
+ G_fatal_error("Cannot open vector map <%s>", in_vector);
+
+ if(in_table)
+ sprintf(out_table, "%s",in_table);
+ else
+ sprintf(out_table, "%s_new",in_vector);
+ */
+
+ sprintf(out_table, "%s",in_table);
+
+ /* init */
+
+ db_init_string(&db_sql);
+ db_init_string(&val_string);
+ db_init_string(&table_name);
+ db_init_handle(&handle);
+
+ /* string to db */
+
+ db_get_connection(&conn);
+ driver=db_start_driver_open_database(conn.driverName, conn.databaseName);
+
+ if(db_table_exists(conn.driverName, conn.databaseName,out_table)>0)
+ G_fatal_error("table %s exists. Choose different table name or check and remove existing table",out_table);
+
+ /* creating table */
+
+ switch (max_trib) {
+ case 2:
+ tab_prev_streams="prev_str01 integer, prev_str02 integer";
+ break;
+ case 3:
+ tab_prev_streams="prev_str01 integer, prev_str02 integer, prev_str03 integer";
+ break;
+ case 4:
+ tab_prev_streams="prev_str01 integer, prev_str02 integer, prev_str03 integer prev_str04 integer";
+ break;
+ case 5:
+ tab_prev_streams="prev_str01 integer, prev_str02 integer, prev_str03 integer prev_str04 integer, prev_str05 integer";
+ break;
+ default:
+ G_fatal_error("Error with number of tributuaries");
+ break;
+ }
+
+
+ sprintf(buf,"create table %s (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)",
+ out_table,
+ tab_cat_col_name,
+ tab_stream,
+ tab_next_stream,
+ tab_prev_streams,
+ tab_strahler,
+ tab_horton,
+ tab_shreve,
+ tab_hack,
+ tab_length,
+ tab_cumlength,
+ tab_stright,
+ tab_fractal);
+
+ db_set_string(&db_sql, buf);
+
+ if(db_execute_immediate(driver,&db_sql) !=DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error("Cannot create table %s", db_get_string(&db_sql));
+ }
+
+ if (db_create_index2(driver,out_table,tab_cat_col_name) !=DB_OK)
+ G_warning("cannot create index");
+
+ if (db_grant_on_table(driver,out_table,DB_PRIV_SELECT, DB_GROUP |DB_PUBLIC) !=DB_OK)
+ G_fatal_error("cannot grant privileges on table %s", out_table);
+
+ db_begin_transaction(driver);
+
+ for (i=0;i<stream_num;++i) {
+
+ if(s_streams[i].stream<0)
+ continue;
+
+ switch (max_trib) {
+ case 2:
+ sprintf(ins_prev_streams,"%d, %d",s_streams[i].trib[0],s_streams[i].trib[1]);
+ break;
+ case 3:
+ sprintf(ins_prev_streams,"%d ,%d, %d",s_streams[i].trib[0],s_streams[i].trib[1],s_streams[i].trib[2]);
+ break;
+ case 4:
+ sprintf(ins_prev_streams,"%d, %d, %d, %d",s_streams[i].trib[0],s_streams[i].trib[1],s_streams[i].trib[2],s_streams[i].trib[3]);
+ break;
+ case 5:
+ sprintf(ins_prev_streams,"%d, %d, %d, %d, %d",s_streams[i].trib[0],s_streams[i].trib[1],s_streams[i].trib[2],s_streams[i].trib[3],s_streams[i].trib[4]);
+ break;
+ default:
+ G_fatal_error("Error with number of tributuaries");
+ break;
+ }
+
+ sprintf(buf,"insert into %s values \
+ (%d, %d, %d, %s, %d, %d, %d, %d, %f, %f , %f, %f)",
+ out_table,
+ i,
+ s_streams[i].stream,
+ s_streams[i].next_stream,
+ ins_prev_streams,
+ s_streams[i].strahler,
+ s_streams[i].horton,
+ s_streams[i].shreeve,
+ s_streams[i].hack,
+ s_streams[i].length,
+ s_streams[i].accum,
+ s_streams[i].stright,
+ s_streams[i].fractal);
+
+ db_set_string(&db_sql,buf);
+
+ if(db_execute_immediate(driver,&db_sql) !=DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error("Cannot inset new row: %s", db_get_string(&db_sql));
+ }
+ }
+
+db_commit_transaction(driver);
+db_close_database_shutdown_driver(driver);
+
+/*
+ if(Vect_map_check_dblink(&Map,1))
+ Vect_map_del_dblink(&Map,1);
+
+ Vect_map_add_dblink(&Map, 1, NULL, out_table,
+ tab_cat_col_name, conn.driverName, conn.databaseName);
+
+ Vect_hist_command(&Map);
+ Vect_build(&Map);
+ Vect_close(&Map);
+*/
+
+G_message("Table %s created. You can join it to vector created with r.stream.extract \
+ using v.db.connect", out_table);
+ return 0;
+}
Modified: grass-addons/raster/r.stream.order/main.c
===================================================================
--- grass-addons/raster/r.stream.order/main.c 2010-01-05 13:28:07 UTC (rev 40249)
+++ grass-addons/raster/r.stream.order/main.c 2010-01-05 13:45:10 UTC (rev 40250)
@@ -25,9 +25,9 @@
{
struct GModule *module; /* GRASS module for parsing arguments */
- struct Option *in_dir_opt, *in_stm_opt, *in_acc_opt, *out_str_opt, *out_shr_opt, *out_hck_opt, *out_hrt_opt; /* options */
+ struct Option *in_dir_opt, *in_stm_opt, /* *in_vect_opt,*/ *in_table_opt, *in_acc_opt, *out_str_opt, *out_shr_opt, *out_hck_opt, *out_hrt_opt; /* options */
struct Flag *out_back; /* flags */
- int stream_num, i;
+ int i;
/* initialize GIS environment */
G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
@@ -56,7 +56,22 @@
in_dir_opt->gisprompt = "old,cell,raster";
in_dir_opt->description =
"Name of direction input map (r.watershed or r.stream.extract)";
-
+ /*
+ in_vect_opt = G_define_option();
+ in_vect_opt->key = "vector";
+ in_vect_opt->type = TYPE_STRING;
+ in_vect_opt->required = NO;
+ in_vect_opt->gisprompt = "old,vector,vector";
+ in_vect_opt->description =
+ "Name of stream vector file (r.stream.extract output only)";
+ */
+ in_table_opt = G_define_option(); /* optional tabe name */
+ in_table_opt->key = "table";
+ in_table_opt->type = TYPE_STRING;
+ in_table_opt->required = NO;
+ in_table_opt->description =
+ "Name of new table to create";
+
in_acc_opt = G_define_option(); /* input stream mask file - optional */
in_acc_opt->key = "accum"; /* required if strahler stream order is calculated for existing stream network */
in_acc_opt->type = TYPE_STRING;
@@ -124,6 +139,8 @@
/* stores input options to variables */
in_dirs = in_dir_opt->answer;
in_streams = in_stm_opt->answer;
+ /* in_vector = in_vect_opt->answer; */
+ in_table = in_table_opt->answer;
in_accum = in_acc_opt->answer;
/* stores output options to variables */
@@ -135,19 +152,19 @@
if (out_strahler) {
if (G_legal_filename(out_strahler) < 0)
- G_fatal_error(_("<%s> is an illegal file name"), out_strahler);
+ G_fatal_error("<%s> is an illegal file name", out_strahler);
}
if (out_shreeve) {
if (G_legal_filename(out_shreeve) < 0)
- G_fatal_error(_("<%s> is an illegal file name"), out_shreeve);
+ G_fatal_error("<%s> is an illegal file name", out_shreeve);
}
if (out_hack) {
if (G_legal_filename(out_hack) < 0)
- G_fatal_error(_("<%s> is an illegal file name"), out_hack);
+ G_fatal_error("<%s> is an illegal file name", out_hack);
}
if (out_horton) {
if (G_legal_filename(out_horton) < 0)
- G_fatal_error(_("<%s> is an illegal file name"), out_horton);
+ G_fatal_error("<%s> is an illegal file name", out_horton);
}
nrows = G_window_rows();
@@ -159,19 +176,23 @@
stack_max = stream_num; /* stack's size depends on number of streams */
init_streams(stream_num);
find_nodes(stream_num);
- if ((out_hack || out_horton) && !in_accum)
+ if (out_hack || out_horton || in_table)
do_cum_length();
-
- if (out_strahler || out_horton)
+
+ if (out_strahler || out_horton || in_table)
strahler();
- if (out_shreeve)
+ if (out_shreeve || in_table)
shreeve();
- if (out_horton)
+ if (out_horton || in_table)
horton();
- if (out_hack)
+ if (out_hack || in_table)
hack();
- write_maps();
+
+ write_maps();
+ if(in_table)
+ create_table();
+
exit(EXIT_SUCCESS);
}
Modified: grass-addons/raster/r.stream.order/order.c
===================================================================
--- grass-addons/raster/r.stream.order/order.c 2010-01-05 13:28:07 UTC (rev 40249)
+++ grass-addons/raster/r.stream.order/order.c 2010-01-05 13:45:10 UTC (rev 40250)
@@ -32,14 +32,18 @@
{
int i;
s_streams = (STREAM *) G_malloc((stream_num + 1) * sizeof(STREAM));
- for (i = 0; i <= stream_num; ++i) {
+ for (i = 0; i <= stream_num; ++i) {
s_streams[i].next_stream = -1;
s_streams[i].stream = -1;
+ s_streams[i].trib_num = -1;
s_streams[i].strahler = -1;
s_streams[i].shreeve = -1;
s_streams[i].hack = -1;
s_streams[i].horton = -1;
s_streams[i].accum = 0.;
+ s_streams[i].length = 0.;
+ s_streams[i].stright = 0.;
+ s_streams[i].fractal = 0.;
s_streams[i].trib[0] = 0;
s_streams[i].trib[1] = 0;
s_streams[i].trib[2] = 0;
@@ -87,7 +91,7 @@
cur_stream = streams[r][c];
if (cur_stream != next_stream) { /* building hierarchy */
-
+
if (outlets_num > (stream_num - 1))
G_fatal_error(_("Error finding nodes. Stream and direction maps probably do not match..."));
@@ -111,6 +115,7 @@
if (trib_num > 1) { /* adding tributuaries */
s_streams[cur_stream].trib_num = trib_num;
+ //G_message("HERE %d %d", cur_stream, s_streams[cur_stream].trib_num);
for (i = 1; i < 9; ++i) {
if (trib > 4)
@@ -152,8 +157,10 @@
int cur_stream;
float cur_northing, cur_easting;
float next_northing, next_easting;
+ float init_northing, init_easting;//
double cur_length=0.;
double cur_accum=0.;
+
G_message(_("Finding longests streams..."));
G_begin_distance_calculations();
@@ -164,10 +171,14 @@
cur_length=0;
done=1;
+ init_northing = window.north - (r + .5) * window.ns_res; //
+ init_easting = window.west + (c + .5) * window.ew_res; //
+
while(done) {
cur_northing = window.north - (r + .5) * window.ns_res;
cur_easting = window.west + (c + .5) * window.ew_res;
+
d=dirs[r][c];
next_r=r+nextr[d];
next_c=c+nextc[d];
@@ -177,25 +188,35 @@
c + nextc[d] < 0 || c + nextc[d] > (ncols - 1) ||
streams[next_r][next_c]<1) { /* mask */
- next_northing = window.north - (r + .5) * window.ns_res;
- next_easting = window.west + (c + .5) * window.ew_res;
- cur_length=G_distance(next_easting, next_northing, cur_easting, cur_northing);
- s_streams[cur_stream].accum += cur_length;
+ cur_length=(window.ns_res+window.ew_res)/2;
+ s_streams[cur_stream].accum += cur_length;
+ s_streams[cur_stream].length += cur_length;
+ s_streams[cur_stream].stright =
+ G_distance(cur_easting, cur_northing, init_easting, init_northing);
+ s_streams[cur_stream].fractal =
+ s_streams[cur_stream].length/s_streams[cur_stream].stright;
break;
}
-
next_northing = window.north - (next_r + .5) * window.ns_res;
next_easting = window.west + (next_c + .5) * window.ew_res;
cur_length = G_distance(next_easting, next_northing, cur_easting, cur_northing);
s_streams[cur_stream].accum += cur_length;
+ s_streams[cur_stream].length += cur_length;
r=next_r;
c=next_c;
if (streams[next_r][next_c] != cur_stream) {
+ s_streams[cur_stream].stright =
+ G_distance(next_easting, next_northing, init_easting, init_northing);
+ s_streams[cur_stream].fractal =
+ s_streams[cur_stream].length/s_streams[cur_stream].stright;
+ init_northing = cur_northing;
+ init_easting = cur_easting;
+
cur_stream=streams[next_r][next_c];
cur_accum=0;
-
+
for (i = 0; i < s_streams[cur_stream].trib_num; ++i) {
if (s_streams[s_streams[cur_stream].trib[i]].accum == 0) {
done=0;
@@ -205,6 +226,7 @@
if(s_streams[s_streams[cur_stream].trib[i]].accum>cur_accum)
cur_accum=s_streams[s_streams[cur_stream].trib[i]].accum;
} /* end for i */
+ if(!in_accum)
s_streams[cur_stream].accum=cur_accum;
}/* end if */
} /* end while */
More information about the grass-commit
mailing list