[GRASS-SVN] r40635 - in grass/branches/develbranch_6/raster/r.watershed: ram seg

[svn_grass at osgeo.org]

Author: mmetz
Date: 2010-01-24 13:33:05 -0500 (Sun, 24 Jan 2010)
New Revision: 40635

Modified:
   grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
   grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c
   grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c
   grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c
Log:
avoid diagonal flow bias, improved stream extraction

Modified: grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c	2010-01-24 17:45:31 UTC (rev 40634)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_astar.c	2010-01-24 18:33:05 UTC (rev 40635)
@@ -295,8 +295,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/branches/develbranch_6/raster/r.watershed/ram/do_cum.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c	2010-01-24 17:45:31 UTC (rev 40634)
+++ grass/branches/develbranch_6/raster/r.watershed/ram/do_cum.c	2010-01-24 18:33:05 UTC (rev 40635)
@@ -108,7 +108,7 @@
     double dx, dy;
     CELL ele, ele_nbr, aspect, is_worked;
     double prop, max_acc;
-    int workedon, edge;
+    int workedon, edge, flat;
     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 this_index, down_index, nbr_index;
@@ -172,6 +172,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 */
@@ -189,14 +190,17 @@
 		    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];
 			is_null = G_is_c_null_value(&ele_nbr);
 			edge = is_null;
+			if (ele_nbr != ele)
+			    flat = 0;
 			if (!is_null && ele_nbr <= ele) {
 			    if (ele_nbr < ele) {
 				weight[ct_dir] =
@@ -343,7 +347,7 @@
 	    /* start new stream */
 	    value = ABS(value) + 0.5;
 	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
-		swale_cells < 1) {
+		swale_cells < 1 && !flat) {
 		FLAG_SET(swale, r, c);
 		is_swale = 1;
 	    }

Modified: grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c	2010-01-24 17:45:31 UTC (rev 40634)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_astar.c	2010-01-24 18:33:05 UTC (rev 40635)
@@ -113,8 +113,7 @@
 		}
 		/* put neighbour in search list if not yet in */
 		if (in_val == 0 && skip_diag == 0) {
-		    cseg_get(&alt, &alt_up, upr, upc);
-		    add_pt(upr, upc, alt_up, alt_val);
+		    add_pt(upr, upc, alt_nbr[ct_dir], alt_val);
 		    /* flow direction is set here */
 		    drain_val = drain[upr - r + 1][upc - c + 1];
 		    cseg_put(&asp, &drain_val, upr, upc);
@@ -316,8 +315,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/branches/develbranch_6/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c	2010-01-24 17:45:31 UTC (rev 40634)
+++ grass/branches/develbranch_6/raster/r.watershed/seg/do_cum.c	2010-01-24 18:33:05 UTC (rev 40635)
@@ -114,7 +114,7 @@
     double dx, dy;
     CELL ele, ele_nbr, asp_val, asp_val2, cvalue, *worked_nbr;
     double prop, max_acc;
-    int workedon, edge;
+    int workedon, edge, flat;
     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 };
 
@@ -188,6 +188,7 @@
 	    cseg_get(&alt, &ele, r, c);
 	    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 */
@@ -205,15 +206,18 @@
 			swale_cells++;
 		    dseg_get(&wat, &valued, r_nbr, c_nbr);
 		    wat_nbr[ct_dir] = valued;
-		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+		    cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
+		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold &&
+		        ct_dir != np_side && ele_nbr > ele)
 			stream_cells++;
 
 		    bseg_get(&worked, &cvalue, r_nbr, c_nbr);
 		    worked_nbr[ct_dir] = cvalue;
 		    if (worked_nbr[ct_dir] == 0) {
-			cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
 			is_null = G_is_c_null_value(&ele_nbr);
 			edge = is_null;
+			if (ele_nbr != ele)
+			    flat = 0;
 			if (!is_null && ele_nbr <= ele) {
 			    if (ele_nbr < ele) {
 				weight[ct_dir] =
@@ -361,7 +365,7 @@
 	    /* start new stream */
 	    value = ABS(value) + 0.5;
 	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
-		swale_cells < 1) {
+		swale_cells < 1 && !flat) {
 		bseg_put(&swale, &one, r, c);
 		is_swale = 1;
 	    }