[GRASS-SVN] r40634 - in grass/trunk/raster/r.watershed: ram seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jan 24 12:45:32 EST 2010
Author: mmetz
Date: 2010-01-24 12:45:31 -0500 (Sun, 24 Jan 2010)
New Revision: 40634
Modified:
grass/trunk/raster/r.watershed/ram/Gwater.h
grass/trunk/raster/r.watershed/ram/close_maps2.c
grass/trunk/raster/r.watershed/ram/def_basin.c
grass/trunk/raster/r.watershed/ram/do_astar.c
grass/trunk/raster/r.watershed/ram/do_cum.c
grass/trunk/raster/r.watershed/ram/haf_side.c
grass/trunk/raster/r.watershed/ram/init_vars.c
grass/trunk/raster/r.watershed/ram/main.c
grass/trunk/raster/r.watershed/ram/no_stream.c
grass/trunk/raster/r.watershed/ram/slope_len.c
grass/trunk/raster/r.watershed/ram/split_str.c
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/bseg_read.c
grass/trunk/raster/r.watershed/seg/close_maps2.c
grass/trunk/raster/r.watershed/seg/cseg_read.c
grass/trunk/raster/r.watershed/seg/def_basin.c
grass/trunk/raster/r.watershed/seg/do_astar.c
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/dseg_read.c
grass/trunk/raster/r.watershed/seg/haf_side.c
grass/trunk/raster/r.watershed/seg/init_vars.c
grass/trunk/raster/r.watershed/seg/main.c
grass/trunk/raster/r.watershed/seg/no_stream.c
grass/trunk/raster/r.watershed/seg/slope_len.c
grass/trunk/raster/r.watershed/seg/split_str.c
Log:
more debug info for diagonal flow bias, improved stream extraction, comment code
Modified: grass/trunk/raster/r.watershed/ram/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/ram/Gwater.h 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/Gwater.h 2010-01-24 17:45:31 UTC (rev 40634)
@@ -18,7 +18,6 @@
#define AR_SIZE 16
#define AR_INCR 16
-#define SHORT int
#define NOMASK 1
#define MIN_SLOPE .00001
#define MIN_GRADIENT_DEGREES 1
@@ -35,7 +34,7 @@
#define POINT struct points
POINT {
- SHORT r, c; /* , downr, downc */
+ int r, c; /* , downr, downc */
/* int nxt; */
};
@@ -49,7 +48,7 @@
extern int mfd, c_fac, abs_acc, ele_scale;
extern int *heap_index, heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-extern SHORT nrows, ncols;
+extern int nrows, ncols;
extern double half_res, diag, max_length, dep_slope;
extern int bas_thres, tot_parts;
extern CELL n_basins;
@@ -66,11 +65,11 @@
extern double *s_l, *s_g, *l_s;
extern CELL one, zero;
extern double ril_value, d_one, d_zero;
-extern SHORT sides;
-extern SHORT drain[3][3];
-extern SHORT updrain[3][3];
-extern SHORT nextdr[8];
-extern SHORT nextdc[8];
+extern int sides;
+extern int drain[3][3];
+extern int updrain[3][3];
+extern int nextdr[8];
+extern int nextdc[8];
extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
@@ -95,11 +94,11 @@
/* do_astar.c */
int do_astar(void);
-int add_pt(SHORT, SHORT, CELL, CELL);
+int add_pt(int, int, CELL, CELL);
int drop_pt(void);
int sift_up(int, CELL);
-double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
-int replace(SHORT, SHORT, SHORT, SHORT);
+double get_slope(int, int, int, int, CELL, CELL);
+int replace(int, int, int, int);
/* do_cum.c */
int do_cum(void);
@@ -110,7 +109,7 @@
int find_pourpts(void);
/* haf_side.c */
-int haf_basin_side(SHORT, SHORT, SHORT);
+int haf_basin_side(int, int, int);
/* init_vars.c */
int init_vars(int, char *[]);
@@ -130,7 +129,7 @@
int len_slp_equ(double, double, double, int, int);
/* slope_len.c */
-int slope_length(SHORT, SHORT, SHORT, SHORT);
+int slope_length(int, int, int, int);
/* split_str.c */
CELL split_stream(int, int, int[], int[], int, CELL, double, CELL);
Modified: grass/trunk/raster/r.watershed/ram/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/close_maps2.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/close_maps2.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -29,7 +29,7 @@
Rast_init_colors(&colors);
Rast_make_random_colors(&colors, 1, max);
- if (max < 10000) {
+ if (max < 1000) {
Rast_set_c_color((CELL) 0, 0, 0, 0, &colors);
r = 1;
incr = 0;
@@ -53,10 +53,6 @@
flag = 0;
incr = -1;
}
- if (r % 200 == 0)
- G_debug(5,
- "adjusting colors: r=%d\tof %d basins",
- r, max);
}
}
}
Modified: grass/trunk/raster/r.watershed/ram/def_basin.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/def_basin.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/def_basin.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -5,7 +5,7 @@
{
int r, rr, c, cc, ct, new_r[9], new_c[9];
CELL downdir, direction, asp_value, value, new_elev;
- SHORT oldupdir, riteflag, leftflag, thisdir;
+ int oldupdir, riteflag, leftflag, thisdir;
for (;;) {
bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
@@ -48,7 +48,7 @@
if (direction == drain[rr][cc]) {
thisdir = updrain[rr][cc];
switch (haf_basin_side
- (oldupdir, (SHORT) downdir, thisdir)) {
+ (oldupdir, (int) downdir, thisdir)) {
case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_elev);
Modified: grass/trunk/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_astar.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/do_astar.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -8,7 +8,7 @@
int do_astar(void)
{
int count;
- SHORT upr, upc, r, c, ct_dir;
+ int upr, upc, r, c, ct_dir;
CELL alt_val, alt_nbr[8];
CELL is_in_list, is_worked;
int index_doer, index_up;
@@ -23,10 +23,10 @@
double slope[8];
int skip_diag;
- G_message(_("SECTION 2: A * Search."));
+ G_message(_("SECTION 2: A* Search."));
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
+ /* get r, c (upr, upc) for neighbours */
upr = nextdr[ct_dir];
upc = nextdc[ct_dir];
/* account for rare cases when ns_res != ew_res */
@@ -48,7 +48,7 @@
while (heap_size > 0) {
G_percent(count++, do_points, 1);
- /* start with point with lowest elevation, in case of equal elevation
+ /* get point with lowest elevation, in case of equal elevation
* of following points, oldest point = point added earliest */
index_doer = astar_pts[1];
@@ -66,13 +66,13 @@
FLAG_SET(worked, r, c);
- /* check neighbours */
+ /* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (upr, upc) for this neighbour */
upr = r + nextdr[ct_dir];
upc = c + nextdc[ct_dir];
slope[ct_dir] = alt_nbr[ct_dir] = 0;
- /* check that r, c are within region */
+ /* check if r, c are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
index_up = SEG_INDEX(alt_seg, upr, upc);
is_in_list = FLAG_GET(in_list, upr, upc);
@@ -145,8 +145,8 @@
return 0;
}
-/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
+/* add point routine for min heap */
+int add_pt(int r, int c, CELL ele, CELL downe)
{
FLAG_SET(in_list, r, c);
@@ -165,7 +165,7 @@
return 0;
}
-/* new drop point routine for min heap */
+/* drop point routine for min heap */
int drop_pt(void)
{
register int child, childr, parent;
@@ -182,10 +182,9 @@
parent = 1;
/* sift down: move hole back towards bottom of heap */
- /* sift-down routine customised for A * Search logic */
while ((child = GET_CHILD(parent)) <= heap_size) {
- /* select child with lower ele, if equal, older child
+ /* select child with lower ele, if both are equal, older child
* older child is older startpoint for flow path, important */
ele = alt[astar_pts[child]];
if (child < heap_size) {
@@ -193,21 +192,12 @@
i = child + 3;
while (childr <= heap_size && childr < i) {
eler = alt[astar_pts[childr]];
- /* get smallest child */
if (cmp_pnt(eler, ele, heap_index[childr], heap_index[child])) {
child = childr;
ele = eler;
}
childr++;
}
- /* break if childr > last entry? that saves sifting up again
- * OTOH, this is another comparison
- * we have a max heap height of 20: log(INT_MAX)/log(n children per node)
- * that would give us in the worst case 20*2 additional comparisons with 3 children
- * the last entry will never go far up again, less than half the way
- * so the additional comparisons for going all the way down
- * and then a bit up again are likely less than 20*2 */
- /* find the error in this reasoning */
}
/* move hole down */
@@ -271,7 +261,7 @@
}
double
-get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+get_slope(int r, int c, int downr, int downc, CELL ele, CELL downe)
{
double slope;
@@ -290,15 +280,15 @@
double get_slope2(CELL ele, CELL up_ele, double dist)
{
- if (ele == up_ele)
- return 0.5 / dist;
+ if (ele >= up_ele)
+ return 0.0;
else
return (double)(up_ele - ele) / dist;
}
/* replace is unused */
int replace( /* ele was in there */
- SHORT upr, SHORT upc, SHORT r, SHORT c)
+ int upr, int upc, int r, int c)
/* CELL ele; */
{
int now, heap_run;
Modified: grass/trunk/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/do_cum.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/do_cum.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -6,11 +6,11 @@
int do_cum(void)
{
- SHORT r, c, dr, dc;
+ int r, c, dr, dc;
CELL is_swale, value, valued, aspect;
int killer, threshold, count;
- SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
int this_index, down_index;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -110,9 +110,9 @@
double dx, dy;
CELL ele, ele_nbr, aspect, is_worked;
double prop, max_acc;
- int workedon, edge;
- SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int workedon, edge, flat;
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
int this_index, down_index, nbr_index;
G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
@@ -174,6 +174,7 @@
ele = alt[this_index];
is_null = 0;
edge = 0;
+ flat = 1;
/* this loop is needed to get the sum of weights */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -184,6 +185,9 @@
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
+ if (dr == r_nbr && dc == c_nbr)
+ np_side = ct_dir;
+
nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
/* check for swale or stream cells */
@@ -191,12 +195,15 @@
if (is_swale)
swale_cells++;
valued = wat[nbr_index];
- if ((ABS(valued) + 0.5) >= threshold)
+ ele_nbr = alt[nbr_index];
+ if ((ABS(valued) + 0.5) >= threshold &&
+ ct_dir != np_side && ele_nbr > ele)
stream_cells++;
is_worked = FLAG_GET(worked, r_nbr, c_nbr);
if (is_worked == 0) {
- ele_nbr = alt[nbr_index];
+ if (ele_nbr != ele)
+ flat = 0;
is_null = Rast_is_c_null_value(&ele_nbr);
edge = is_null;
if (!is_null && ele_nbr <= ele) {
@@ -223,8 +230,6 @@
}
}
}
- if (dr == r_nbr && dc == c_nbr)
- np_side = ct_dir;
}
else
edge = 1;
@@ -260,7 +265,6 @@
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
@@ -273,7 +277,7 @@
nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
weight[ct_dir] = weight[ct_dir] / sum_weight;
- /* check everything sums up to 1.0 */
+ /* check everything adds up to 1.0 */
prop += weight[ct_dir];
valued = wat[nbr_index];
@@ -344,8 +348,8 @@
}
/* start new stream */
value = ABS(value) + 0.5;
- if (!is_swale && (int)value >= threshold && stream_cells < 3 &&
- swale_cells < 1) {
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+ swale_cells < 1 && !flat) {
FLAG_SET(swale, r, c);
is_swale = 1;
}
Modified: grass/trunk/raster/r.watershed/ram/haf_side.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/haf_side.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/haf_side.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -1,8 +1,8 @@
#include "Gwater.h"
-int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+int haf_basin_side(int updir, int downdir, int thisdir)
{
- SHORT newup, newthis;
+ int newup, newthis;
newup = updir - downdir;
if (newup < 0)
Modified: grass/trunk/raster/r.watershed/ram/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/init_vars.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/init_vars.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -9,7 +9,7 @@
int init_vars(int argc, char *argv[])
{
- SHORT r, c;
+ int r, c;
CELL *buf, alt_value, wat_value, asp_value, block_value;
DCELL dvalue;
void *elebuf, *ptr;
@@ -34,8 +34,6 @@
c_fac = 5;
abs_acc = 0;
ele_scale = 1;
- for (r = 0; r < argc; r++)
- G_debug(0, "%s", argv[r]);
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "elevation=%s", ele_name) == 1)
Modified: grass/trunk/raster/r.watershed/ram/main.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/main.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/main.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -27,7 +27,7 @@
int mfd, c_fac, abs_acc, ele_scale;
int *heap_index, heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
-SHORT nrows, ncols;
+int nrows, ncols;
double half_res, diag, max_length, dep_slope;
int bas_thres, tot_parts;
CELL n_basins;
@@ -44,11 +44,11 @@
double *s_l, *s_g, *l_s;
CELL one, zero;
double ril_value, d_one, d_zero;
-SHORT sides;
-SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
-SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
-SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
-SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+int sides;
+int drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+int updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
char run_name[GNAME_MAX], ob_name[GNAME_MAX];
char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
Modified: grass/trunk/raster/r.watershed/ram/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/no_stream.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/no_stream.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -7,7 +7,7 @@
double slope;
CELL downdir, asp_value, hih_ele, new_ele, aspect;
DCELL dvalue, max_drain; /* flow acc is now DCELL */
- SHORT updir, riteflag, leftflag, thisdir;
+ int updir, riteflag, leftflag, thisdir;
while (1) {
bas[SEG_INDEX(bas_seg, row, col)] = basin_num;
@@ -67,7 +67,7 @@
if (aspect == drain[rr][cc]) {
thisdir = updrain[rr][cc];
switch (haf_basin_side
- (updir, (SHORT) downdir, thisdir)) {
+ (updir, (int) downdir, thisdir)) {
case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_ele);
Modified: grass/trunk/raster/r.watershed/ram/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/slope_len.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/slope_len.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -1,6 +1,6 @@
#include "Gwater.h"
-int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+int slope_length(int r, int c, int dr, int dc)
{
CELL top_alt, bot_alt, asp_value;
double res, top_ls, bot_ls;
Modified: grass/trunk/raster/r.watershed/ram/split_str.c
===================================================================
--- grass/trunk/raster/r.watershed/ram/split_str.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/ram/split_str.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -6,8 +6,8 @@
{
CELL downdir, old_basin, new_elev, aspect;
double slope, easting, northing;
- SHORT doit, ctr, updir, splitdir[9];
- SHORT thisdir, leftflag, riteflag;
+ int doit, ctr, updir, splitdir[9];
+ int thisdir, leftflag, riteflag;
int r, c, rr, cc;
new_elev = 0;
@@ -35,7 +35,7 @@
if (doit) {
thisdir = updrain[rr][cc];
switch (haf_basin_side
- (updir, (SHORT) downdir, thisdir)) {
+ (updir, (int) downdir, thisdir)) {
case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_elev);
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2010-01-24 17:45:31 UTC (rev 40634)
@@ -15,7 +15,6 @@
#define AR_SIZE 16
#define AR_INCR 16
-#define SHORT int
#define NOMASK 1
#define MIN_SLOPE .00001
#define MIN_GRADIENT_DEGREES 1
@@ -45,7 +44,7 @@
#define POINT struct points
POINT {
- SHORT r, c; /* , downr, downc */
+ int r, c; /* , downr, downc */
char asp; /* drainage direction */
char guessed; /* accumulation will likely be an underestimate */
};
@@ -77,7 +76,7 @@
extern CELL n_basins;
extern OC_STACK *ocs;
extern int ocs_alloced;
-extern SHORT nrows, ncols;
+extern int nrows, ncols;
extern double half_res, diag, max_length, dep_slope;
extern int bas_thres, tot_parts;
extern SSEG astar_pts;
@@ -88,11 +87,11 @@
extern double segs_mb;
extern char zero, one;
extern double ril_value, d_zero, d_one;
-extern SHORT sides;
-extern SHORT drain[3][3];
-extern SHORT updrain[3][3];
-extern SHORT nextdr[8];
-extern SHORT nextdc[8];
+extern int sides;
+extern int drain[3][3];
+extern int updrain[3][3];
+extern int nextdr[8];
+extern int nextdc[8];
extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
@@ -140,14 +139,14 @@
/* do_astar.c */
int do_astar(void);
-int add_pt(SHORT, SHORT, CELL, char, int);
+int add_pt(int, int, CELL, char, int);
HEAP_PNT drop_pt(void);
int sift_up(int, HEAP_PNT);
int sift_up_mem(int, HEAP_PNT);
int sift_down_disk(void);
int fill_mem_heap(void);
int cmp_pnt(HEAP_PNT *a, HEAP_PNT *b);
-double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+double get_slope(int, int, int, int, CELL, CELL);
/* do_cum.c */
int do_cum(void);
@@ -161,7 +160,7 @@
int find_pourpts(void);
/* haf_side.c */
-int haf_basin_side(SHORT, SHORT, SHORT);
+int haf_basin_side(int, int, int);
/* init_vars.c */
int init_vars(int, char *[]);
@@ -178,7 +177,7 @@
int len_slp_equ(double, double, double, int, int);
/* slope_len.c */
-int slope_length(SHORT, SHORT, SHORT, SHORT);
+int slope_length(int, int, int, int);
/* split_str.c */
CELL split_stream(int, int, int[], int[], int, CELL, double, CELL);
Modified: grass/trunk/raster/r.watershed/seg/bseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_read.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/bseg_read.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -8,7 +8,6 @@
int row, nrows;
int col, ncols;
int map_fd;
- char msg[100];
CELL *buffer;
char cbuf;
Modified: grass/trunk/raster/r.watershed/seg/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps2.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/close_maps2.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -35,7 +35,7 @@
Rast_init_colors(&colors);
Rast_make_random_colors(&colors, 1, max);
- if (max < 10000) {
+ if (max < 1000) {
Rast_set_c_color((CELL) 0, 0, 0, 0, &colors);
r = 1;
incr = 0;
Modified: grass/trunk/raster/r.watershed/seg/cseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/cseg_read.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/cseg_read.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -9,7 +9,6 @@
{
int row, nrows;
int map_fd;
- char msg[100];
CELL *buffer;
cseg->name = NULL;
@@ -23,9 +22,8 @@
if (segment_put_row(&(cseg->seg), buffer, row) < 0) {
G_free(buffer);
Rast_close(map_fd);
- sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+ G_warning("%s(): unable to segment put row for [%s] in [%s]",
me, map_name, mapset);
- G_warning(msg);
return (-1);
}
}
Modified: grass/trunk/raster/r.watershed/seg/def_basin.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/def_basin.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/def_basin.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -7,7 +7,7 @@
int r, rr, c, cc, ct, new_r[9], new_c[9];
CELL downdir, direction, asp_value, value, new_elev;
char cvalue;
- SHORT oldupdir, riteflag, leftflag, thisdir;
+ int oldupdir, riteflag, leftflag, thisdir;
for (;;) {
cseg_put(&bas, &basin_num, row, col);
@@ -49,7 +49,7 @@
if (direction == drain[rr][cc]) {
thisdir = updrain[rr][cc];
switch (haf_basin_side
- (oldupdir, (SHORT) downdir, thisdir)) {
+ (oldupdir, (int) downdir, thisdir)) {
case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_elev);
Modified: grass/trunk/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/do_astar.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -8,7 +8,7 @@
int do_astar(void)
{
int doer, count;
- SHORT upr, upc, r = -1, c = -1, ct_dir;
+ int upr, upc, r = -1, c = -1, ct_dir;
CELL alt_val, alt_nbr[8], alt_up;
CELL asp_val;
char flag_value, is_in_list, is_worked;
@@ -24,11 +24,12 @@
double dx, dy, dist_to_nbr[8], ew_res, ns_res;
double slope[8];
int skip_diag;
+ int count_edge = 0, count_diag = 0, count_edge_sink = 0, count_diag_sink = 0;
G_message(_("SECTION 2: A* Search."));
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
+ /* get r, c (upr, upc) for neighbours */
upr = nextdr[ct_dir];
upc = nextdc[ct_dir];
/* account for rare cases when ns_res != ew_res */
@@ -51,13 +52,10 @@
doer = do_points - 1;
- /* A * Search
- * determine downhill path through all not masked cells and
- * set drainage direction */
+ /* A* Search: search uphill, get downhill paths */
while (heap_size > 0) {
G_percent(count++, do_points, 1);
- /* drop next point from heap */
heap_p = drop_pt();
r = heap_p.pnt.r;
@@ -72,13 +70,13 @@
upr = r + nextdr[ct_dir];
upc = c + nextdc[ct_dir];
slope[ct_dir] = alt_nbr[ct_dir] = 0;
- /* check that upr, upc are within region */
+ /* check if upr, upc are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
- /* avoid diagonal flow direction bias */
bseg_get(&bitflags, &flag_value, upr, upc);
is_in_list = FLAG_GET(flag_value, INLISTFLAG);
is_worked = FLAG_GET(flag_value, WORKEDFLAG);
skip_diag = 0;
+ /* avoid diagonal flow direction bias */
if (!is_worked) {
cseg_get(&alt, &alt_up, upr, upc);
alt_nbr[ct_dir] = alt_up;
@@ -105,17 +103,31 @@
}
}
- /* put neighbour in search list if not yet in */
+ /* add neighbour as new point if not in the list */
if (is_in_list == 0 && skip_diag == 0) {
- cseg_get(&alt, &alt_up, upr, upc);
- /* flow direction is set here */
+ /* set flow direction */
asp_val = drain[upr - r + 1][upc - c + 1];
- add_pt(upr, upc, alt_up, asp_val, 0);
+ add_pt(upr, upc, alt_nbr[ct_dir], asp_val, 0);
cseg_put(&asp, &asp_val, upr, upc);
+
+
+ if (alt_nbr[ct_dir] < alt_val) {
+ if (ct_dir < 4)
+ count_edge_sink++;
+ else
+ count_diag_sink++;
+ }
+ /* includes flat areas */
+ else {
+ if (ct_dir < 4)
+ count_edge++;
+ else
+ count_diag++;
+ }
}
- /* in list, not worked. is it edge ? */
else if (is_in_list && is_worked == 0 &&
FLAG_GET(flag_value, EDGEFLAG)) {
+ /* neighbour is edge in list, not yet worked */
cseg_get(&asp, &asp_val, upr, upc);
if (asp_val < 0) {
/* adjust flow direction for edge cell */
@@ -126,6 +138,7 @@
}
}
}
+ /* add astar points to sorted list for flow accumulation */
cseg_get(&asp, &asp_val, r, c);
heap_p.pnt.asp = asp_val;
seg_put(&astar_pts, (char *)&heap_p.pnt, 0, doer);
@@ -143,6 +156,11 @@
G_percent(count, do_points, 1); /* finish it */
+ G_debug(1, "edge direction: %d (%.2f%%)", count_edge, (double) 100. * count_edge / (count_edge + count_diag));
+ G_debug(1, "diag direction: %d (%.2f%%)", count_diag, (double) 100. * count_diag / (count_edge + count_diag));
+ G_debug(1, "edge out of depression: %d (%.2f%%)", count_edge_sink, (double) 100. * count_edge_sink / (count_edge_sink + count_diag_sink));
+ G_debug(1, "diag out of depression: %d (%.2f%%)", count_diag_sink, (double) 100. * count_diag_sink / (count_edge_sink + count_diag_sink));
+
return 0;
}
@@ -160,7 +178,7 @@
}
/* add point routine for min heap */
-int add_pt(SHORT r, SHORT c, CELL ele, char asp, int is_edge)
+int add_pt(int r, int c, CELL ele, char asp, int is_edge)
{
HEAP_PNT heap_p;
char flag_value;
@@ -204,7 +222,7 @@
seg_get(&search_heap, (char *)&last_p, 0, heap_size);
seg_get(&search_heap, (char *)&root_p, 0, 1);
- /* sift down */
+ /* sift down: move hole back towards bottom of heap */
parent = 1;
while ((child = GET_CHILD(parent)) < heap_size) {
/* select child with lower ele, if both are equal, older child
@@ -215,7 +233,6 @@
childr = child + 1;
i = child + 4;
while (childr < i && childr < heap_size) {
- /* get smallest child */
seg_get(&search_heap, (char *)&childr_p, 0, childr);
if (cmp_pnt(&childr_p, &child_p)) {
child = childr;
@@ -245,7 +262,7 @@
/* standard sift-up routine for d-ary min heap */
int sift_up(int start, HEAP_PNT child_p)
{
- int parent, child; /* parentp, */
+ int parent, child;
HEAP_PNT heap_p;
child = start;
@@ -271,7 +288,7 @@
}
double
-get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+get_slope(int r, int c, int downr, int downc, CELL ele, CELL downe)
{
double slope;
@@ -288,8 +305,8 @@
double get_slope2(CELL ele, CELL up_ele, double dist)
{
- if (ele == up_ele)
- return 0.5 / dist;
+ if (ele >= up_ele)
+ return 0.0;
else
return (double)(up_ele - ele) / dist;
}
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -7,14 +7,14 @@
int do_cum(void)
{
- SHORT r, c, dr, dc;
+ int r, c, dr, dc;
CELL asp_val, asp_val_down;
char is_swale, this_flag_value, flag_value;
DCELL value, valued;
POINT point;
int killer, threshold;
- SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
WAT_ALT wa, wadown;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
@@ -29,11 +29,11 @@
r = point.r;
c = point.c;
asp_val = point.asp;
- /* skip user-defined depressions */
if (asp_val) {
dr = r + asp_r[ABS(asp_val)];
dc = c + asp_c[ABS(asp_val)];
}
+ /* skip user-defined depressions */
else
dr = dc = -1;
@@ -130,10 +130,10 @@
double dx, dy;
CELL ele, *ele_nbr, asp_val, asp_val_down;
double prop, max_acc;
- int workedon, edge, is_swale;
+ int workedon, edge, is_swale, flat;
char *flag_nbr, this_flag_value, flag_value;
- SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
- SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
+ int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+ int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
G_debug(1, "MFD convergence factor set to %d.", c_fac);
@@ -206,6 +206,7 @@
ele = wa.ele;
is_null = 0;
edge = 0;
+ flat = 1;
/* this loop is needed to get the sum of weights */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -224,6 +225,9 @@
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
+ if (dr == r_nbr && dc == c_nbr)
+ np_side = ct_dir;
+
bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
wat_nbr[ct_dir] = wa.wat;
@@ -233,11 +237,15 @@
is_swale = FLAG_GET(flag_nbr[ct_dir], SWALEFLAG);
if (is_swale)
swale_cells++;
- if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+ if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold &&
+ ct_dir != np_side && ele_nbr[ct_dir] > ele)
stream_cells++;
if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
+ if (ele_nbr[ct_dir] != ele)
+ flat = 0;
+
edge = is_null = FLAG_GET(flag_nbr[ct_dir], NULLFLAG);
if (!is_null && ele_nbr[ct_dir] <= ele) {
if (ele_nbr[ct_dir] < ele) {
@@ -263,8 +271,6 @@
}
}
}
- if (dr == r_nbr && dc == c_nbr)
- np_side = ct_dir;
}
else
edge = 1;
@@ -302,7 +308,6 @@
if (mfd_cells > 1) {
prop = 0.0;
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
- /* get r, c (r_nbr, c_nbr) for neighbours */
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
@@ -313,7 +318,7 @@
if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
weight[ct_dir] = weight[ct_dir] / sum_weight;
- /* check everything sums up to 1.0 */
+ /* check everything adds up to 1.0 */
prop += weight[ct_dir];
if (value > 0) {
@@ -393,9 +398,8 @@
is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
/* start new stream */
value = ABS(value) + 0.5;
- /* can use stream_cells < 4 only for highres, nsres and ewres < 30 m? */
- if (!is_swale && (int)value >= threshold && stream_cells < 3 &&
- swale_cells < 1) {
+ if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
+ swale_cells < 1 && !flat) {
FLAG_SET(this_flag_value, SWALEFLAG);
is_swale = 1;
}
@@ -426,9 +430,7 @@
G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
workedon, do_points);
- if (!st_flag) {
- seg_close(&astar_pts);
- }
+ seg_close(&astar_pts);
G_free(dist_to_nbr);
G_free(weight);
Modified: grass/trunk/raster/r.watershed/seg/dseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_read.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/dseg_read.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -9,7 +9,6 @@
{
int row, nrows, ncols;
int map_fd;
- char msg[100];
double *dbuffer;
dseg->name = NULL;
@@ -24,9 +23,8 @@
if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
G_free(dbuffer);
Rast_close(map_fd);
- sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
+ G_warning("%s(): unable to segment put row for [%s] in [%s]",
me, map_name, mapset);
- G_warning(msg);
return (-1);
}
}
Modified: grass/trunk/raster/r.watershed/seg/haf_side.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/haf_side.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/haf_side.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -1,8 +1,8 @@
#include "Gwater.h"
-int haf_basin_side(SHORT updir, SHORT downdir, SHORT thisdir)
+int haf_basin_side(int updir, int downdir, int thisdir)
{
- SHORT newup, newthis;
+ int newup, newthis;
newup = updir - downdir;
if (newup < 0)
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -9,7 +9,7 @@
int init_vars(int argc, char *argv[])
{
- SHORT r, c;
+ int r, c;
int ele_fd, fd, wat_fd;
int seg_rows, seg_cols, num_cseg_total, num_open_segs, num_open_array_segs;
@@ -116,13 +116,6 @@
/* do RUSLE */
if (er_flag)
tot_parts++;
- /* do stream extraction */
- if (er_flag || seg_flag || bas_flag || haf_flag) {
- st_flag = 1;
- /* separate stream extraction only needed for MFD */
- if (mfd)
- tot_parts++;
- }
/* define basins */
if (seg_flag || bas_flag || haf_flag)
tot_parts++;
@@ -288,7 +281,7 @@
else if (wat_map_type == FCELL_TYPE) {
wat_value = *((FCELL *)ptr);
}
- else if (ele_map_type == DCELL_TYPE) {
+ else if (wat_map_type == DCELL_TYPE) {
wat_value = *((DCELL *)ptr);
}
}
@@ -324,8 +317,7 @@
G_free(watbuf);
}
- if (do_points < nrows * ncols)
- MASK_flag = 1;
+ MASK_flag = (do_points < nrows * ncols);
/* depression: drainage direction will be set to zero later */
if (pit_flag) {
@@ -388,7 +380,11 @@
G_debug(1, "open segments for A* points");
/* rounded down power of 2 */
seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
+ if (seg_cols < 2)
+ seg_cols = 2;
num_open_array_segs = num_open_segs / seg_cols;
+ if (num_open_array_segs == 0)
+ num_open_array_segs = 1;
/* n cols in segment */
seg_cols *= seg_rows * seg_rows;
/* n segments in row */
@@ -409,7 +405,11 @@
G_debug(1, "open segments for A* search heap");
/* rounded down power of 2 */
seg_cols = (int) (pow(2, (int)(log(num_open_segs / 8.0) / log(2) + 0.1)) + 0.1);
+ if (seg_cols < 2)
+ seg_cols = 2;
num_open_array_segs = num_open_segs / seg_cols;
+ if (num_open_array_segs == 0)
+ num_open_array_segs = 1;
/* n cols in segment */
seg_cols *= seg_rows * seg_rows;
/* n segments in row */
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/main.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -32,7 +32,7 @@
CELL n_basins;
OC_STACK *ocs;
int ocs_alloced;
-SHORT nrows, ncols;
+int nrows, ncols;
double half_res, diag, max_length, dep_slope;
int bas_thres, tot_parts;
SSEG astar_pts;
@@ -43,11 +43,11 @@
double segs_mb;
char zero, one;
double ril_value, d_zero, d_one;
-SHORT sides;
-SHORT drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
-SHORT updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
-SHORT nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
-SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
+int sides;
+int drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
+int updrain[3][3] = { {3, 2, 1}, {4, 0, 8}, {5, 6, 7} };
+int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
+int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
char run_name[GNAME_MAX], ob_name[GNAME_MAX];
char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
Modified: grass/trunk/raster/r.watershed/seg/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/no_stream.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/no_stream.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -8,7 +8,7 @@
double slope;
CELL downdir, asp_value, hih_ele, new_ele, aspect, value;
DCELL dvalue, max_drain; /* flow acc is now DCELL */
- SHORT updir, riteflag, leftflag, thisdir;
+ int updir, riteflag, leftflag, thisdir;
WAT_ALT wa;
while (1) {
@@ -73,7 +73,7 @@
if (aspect == drain[rr][cc]) {
thisdir = updrain[rr][cc];
switch (haf_basin_side(updir,
- (SHORT) downdir,
+ (int) downdir,
thisdir)) {
case RITE:
overland_cells(r, c, basin_num, basin_num,
Modified: grass/trunk/raster/r.watershed/seg/slope_len.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/slope_len.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/slope_len.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -1,6 +1,6 @@
#include "Gwater.h"
-int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
+int slope_length(int r, int c, int dr, int dc)
{
CELL top_alt, bot_alt, asp_value, ridge;
double res, top_ls, bot_ls;
Modified: grass/trunk/raster/r.watershed/seg/split_str.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/split_str.c 2010-01-24 17:25:42 UTC (rev 40633)
+++ grass/trunk/raster/r.watershed/seg/split_str.c 2010-01-24 17:45:31 UTC (rev 40634)
@@ -6,8 +6,8 @@
{
CELL downdir, old_basin, new_elev, aspect;
double slope, easting, northing;
- SHORT doit, ctr, updir, splitdir[9];
- SHORT thisdir, leftflag, riteflag;
+ int doit, ctr, updir, splitdir[9];
+ int thisdir, leftflag, riteflag;
int r, c, rr, cc;
new_elev = 0;
@@ -35,7 +35,7 @@
if (doit) {
thisdir = updrain[rr][cc];
switch (haf_basin_side
- (updir, (SHORT) downdir, thisdir)) {
+ (updir, (int) downdir, thisdir)) {
case LEFT:
overland_cells(r, c, basin_num, basin_num - 1,
&new_elev);
More information about the grass-commit
mailing list