[GRASS-SVN] r39603 - grass/trunk/raster/r.watershed/seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Oct 21 11:51:41 EDT 2009
Author: mmetz
Date: 2009-10-21 11:51:41 -0400 (Wed, 21 Oct 2009)
New Revision: 39603
Added:
grass/trunk/raster/r.watershed/seg/flag.h
Modified:
grass/trunk/raster/r.watershed/seg/Gwater.h
grass/trunk/raster/r.watershed/seg/bseg_get.c
grass/trunk/raster/r.watershed/seg/bseg_open.c
grass/trunk/raster/r.watershed/seg/bseg_put.c
grass/trunk/raster/r.watershed/seg/bseg_read.c
grass/trunk/raster/r.watershed/seg/bseg_write.c
grass/trunk/raster/r.watershed/seg/close_maps.c
grass/trunk/raster/r.watershed/seg/close_maps2.c
grass/trunk/raster/r.watershed/seg/cseg.h
grass/trunk/raster/r.watershed/seg/cseg_write.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_astar.h
grass/trunk/raster/r.watershed/seg/do_cum.c
grass/trunk/raster/r.watershed/seg/dseg_write.c
grass/trunk/raster/r.watershed/seg/find_pour.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/sg_factor.c
Log:
improved segmented mode
Modified: grass/trunk/raster/r.watershed/seg/Gwater.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/Gwater.h 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/Gwater.h 2009-10-21 15:51:41 UTC (rev 39603)
@@ -10,6 +10,7 @@
#include <grass/gis.h>
#include <grass/raster.h>
#include "cseg.h"
+#include "flag.h"
#define AR_SIZE 16
#define AR_INCR 16
@@ -21,8 +22,8 @@
#define METER_TO_FOOT (1 / 0.3048)
#define MAX_BYTES 10485760
#define PAGE_BLOCK 1024
-#define SROW 200
-#define SCOL 200
+#define SROW 64
+#define SCOL 64
#define RITE 1
#define LEFT 2
#define NEITHER 0
@@ -30,33 +31,53 @@
#define TSTSTR(a) (fprintf (stderr, "%s\n", a))
#define TST(a) (fprintf (stderr, "%e\n", (double) (a)))
+/* flag positions */
+#define NULLFLAG 0 /* elevation is NULL */
+#define EDGEFLAG 1 /* edge cell */
+#define INLISTFLAG 2 /* in open A* list */
+#define WORKEDFLAG 3 /* in closed A* list/ accumulation done/streams done */
+#define SWALEFLAG 4 /* swale */
+#define PITFLAG 5 /* user-defined real depression */
+#define RUSLEBLOCKFLAG 6 /* is RUSLE block */
+/* #define XXXFLAG 7 */ /* last bit unused */
+
+
#define POINT struct points
POINT {
SHORT r, c; /* , downr, downc */
- int nxt;
+ char asp; /* drainage direction */
+ char guessed; /* accumulation will likely be an underestimate */
};
-#define HEAP struct heap_item
-HEAP {
- int point;
+#define HEAP_PNT struct heap_point
+HEAP_PNT {
+ int added;
CELL ele;
+ POINT pnt;
};
+#define WAT_ALT struct wat_altitude
+WAT_ALT {
+ CELL ele;
+ DCELL wat;
+};
+
extern struct Cell_head window;
extern int mfd, c_fac, abs_acc, ele_scale;
-extern SSEG heap_index;
+extern SSEG search_heap;
extern int heap_size;
extern int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
extern SHORT nrows, ncols;
extern double half_res, diag, max_length, dep_slope;
extern int bas_thres, tot_parts;
extern SSEG astar_pts;
-extern BSEG worked, in_list, s_b, swale;
+extern BSEG bitflags, s_b;
extern CSEG dis, alt, asp, bas, haf, r_h, dep;
-extern DSEG wat;
+extern SSEG watalt;
extern DSEG slp, s_l, s_g, l_s, ril;
-extern CELL one, zero;
+extern double segs_mb;
+extern char zero, one;
extern double ril_value, d_zero, d_one;
extern SHORT sides;
extern SHORT drain[3][3];
@@ -77,6 +98,28 @@
extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
extern FILE *fp;
+/* the flags:
+ * ele_flag elevation map given
+ * pit_flag pit (depression) map given
+ * run_flag initial surface runoff given
+ * dis_flag ???
+ * ob_flag blocking map for RUSLE given
+ * wat_flag write accumulation output
+ * asp_flag write direction output
+ * arm_flag unused, for interactive mode
+ * ril_flag percentage disturbed land given
+ * dep_flag ???
+ * st_flag do stream extraction
+ * bas_flag write basin output
+ * seg_flag write stream output
+ * haf_flag write half-basin output
+ * er_flag do RUSLE
+ * sb_flag ???
+ * sg_flag write RUSLE S factor
+ * sl_flag slope length, unused
+ * ls_flag write RUSLE LS factor
+ */
+
/* close_maps.c */
int close_maps(void);
@@ -88,17 +131,23 @@
/* do_astar.c */
int do_astar(void);
-int add_pt(SHORT, SHORT, CELL, CELL);
-int drop_pt(void);
-int sift_up(int, CELL);
+int add_pt(SHORT, SHORT, 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);
-int replace(SHORT, SHORT, SHORT, SHORT);
/* do_cum.c */
int do_cum(void);
int do_cum_mfd(void);
double mfd_pow(double, int);
+/* do_stream.c */
+int do_stream(void);
+
/* find_pour.c */
int find_pourpts(void);
Modified: grass/trunk/raster/r.watershed/seg/bseg_get.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_get.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/bseg_get.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -2,8 +2,17 @@
#include <grass/segment.h>
#include "cseg.h"
-int bseg_get(BSEG * bseg, CELL * value, int row, int col)
+int bseg_get(BSEG * bseg, char *value, int row, int col)
{
+ if (segment_get(&(bseg->seg), value, row, col) < 0) {
+ G_warning("cseg_get(): could not read segment file");
+ return -1;
+ }
+ return 0;
+}
+
+int bseg_get_old(BSEG * bseg, CELL *value, int row, int col)
+{
CELL x;
if (segment_get(&(bseg->seg), &x, row, col >> 3) < 0) {
Modified: grass/trunk/raster/r.watershed/seg/bseg_open.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_open.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/bseg_open.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -22,7 +22,7 @@
return -2;
}
if (0 > (errflag = segment_format(fd, G_window_rows(),
- (G_window_cols() + 7) / 8, srows, scols,
+ G_window_cols(), srows, scols,
sizeof(char)))) {
close(fd);
unlink(filename);
Modified: grass/trunk/raster/r.watershed/seg/bseg_put.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_put.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/bseg_put.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -2,8 +2,18 @@
#include <grass/segment.h>
#include "cseg.h"
-int bseg_put(BSEG * bseg, CELL * value, int row, int col)
+int bseg_put(BSEG * bseg, char *value, int row, int col)
{
+ if (segment_put(&(bseg->seg), value, row, col) < 0) {
+ G_warning("cseg_put(): could not write segment file");
+ return -1;
+ }
+ return 0;
+}
+
+
+int bseg_put_old(BSEG * bseg, CELL * value, int row, int col)
+{
CELL old_value;
if (segment_get(&(bseg->seg), &old_value, row, col >> 3) < 0) {
Modified: grass/trunk/raster/r.watershed/seg/bseg_read.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_read.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/bseg_read.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -12,6 +12,7 @@
int map_fd;
char msg[100];
CELL *buffer;
+ char cbuf;
bseg->name = NULL;
bseg->mapset = NULL;
@@ -36,7 +37,8 @@
return -2;
}
for (col = ncols; col >= 0; col--) {
- bseg_put(bseg, &(buffer[col]), row, col);
+ cbuf = (char) buffer[col];
+ bseg_put(bseg, &cbuf, row, col);
}
}
Modified: grass/trunk/raster/r.watershed/seg/bseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/bseg_write.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/bseg_write.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -10,7 +10,7 @@
int row, nrows;
int col, ncols;
CELL *buffer;
- CELL value;
+ char value;
map_fd = Rast_open_c_new(map_name);
if (map_fd < 0) {
@@ -21,6 +21,7 @@
ncols = G_window_cols();
buffer = Rast_allocate_c_buf();
for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
for (col = 0; col < ncols; col++) {
bseg_get(bseg, &value, row, col);
buffer[col] = value;
@@ -33,6 +34,7 @@
return -2;
}
}
+ G_percent(row, nrows, 1); /* finish it */
G_free(buffer);
Rast_close(map_fd);
return 0;
Modified: grass/trunk/raster/r.watershed/seg/close_maps.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/close_maps.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -6,65 +6,55 @@
{
struct Colors colors;
int r, c;
- CELL is_swale;
DCELL *dbuf = NULL;
int fd;
struct FPRange accRange;
DCELL min, max;
DCELL clr_min, clr_max;
DCELL sum, sum_sqr, stddev, lstddev, dvalue;
+ WAT_ALT wa;
- dseg_close(&slp);
+ /* bseg_close(&swale); */
cseg_close(&alt);
if (wat_flag) {
+ G_message(_("Closing accumulation map"));
sum = sum_sqr = stddev = 0.0;
dbuf = Rast_allocate_d_buf();
if (abs_acc) {
G_warning("Writing out only positive flow accumulation values.");
G_warning("Cells with a likely underestimate for flow accumulation can no longer be identified.");
-
- fd = Rast_open_new(wat_name, DCELL_TYPE);
- if (fd < 0) {
- G_warning(_("unable to open new accum map layer."));
- }
- for (r = 0; r < nrows; r++) {
- Rast_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
- for (c = 0; c < ncols; c++) {
- dseg_get(&wat, &dvalue, r, c);
- if (Rast_is_d_null_value(&dvalue) == 0 && dvalue) {
- dvalue = ABS(dvalue);
- dbuf[c] = dvalue;
- sum += dvalue;
- sum_sqr += dvalue * dvalue;
- }
- }
- Rast_put_row(fd, dbuf, DCELL_TYPE);
- }
- if (Rast_close(fd) < 0)
- G_warning(_("Close failed."));
}
- else {
- dseg_write_cellfile(&wat, wat_name);
- /* get standard deviation */
- fd = Rast_open_old(wat_name, "");
- if (fd < 0) {
- G_fatal_error(_("unable to open flow accumulation map layer"));
- }
-
- for (r = 0; r < nrows; r++) {
- Rast_get_d_row(fd, dbuf, r);
- for (c = 0; c < ncols; c++) {
- dvalue = dbuf[c];
- if (Rast_is_d_null_value(&dvalue) == 0 && dvalue) {
- dvalue = ABS(dvalue);
+ fd = Rast_open_new(wat_name, DCELL_TYPE);
+ if (fd < 0) {
+ G_warning(_("unable to open new accum map layer."));
+ }
+ for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
+ Rast_set_d_null_value(dbuf, ncols); /* reset row to all NULL */
+ for (c = 0; c < ncols; c++) {
+ /* dseg_get(&wat, &dvalue, r, c); */
+ seg_get(&watalt, (char *)&wa, r, c);
+ dvalue = wa.wat;
+ if (Rast_is_d_null_value(&dvalue) == 0 && dvalue) {
+ if (abs_acc) {
+ dvalue = fabs(dvalue);
sum += dvalue;
- sum_sqr += dvalue * dvalue;
}
+ else
+ sum += fabs(dvalue);
+
+ dbuf[c] = dvalue;
+ sum_sqr += dvalue * dvalue;
}
}
+ Rast_put_row(fd, dbuf, DCELL_TYPE);
}
+ G_percent(r, nrows, 1); /* finish it */
+ if (Rast_close(fd) < 0)
+ G_warning(_("Close failed."));
+
stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
G_debug(1, "stddev: %f", stddev);
@@ -132,55 +122,31 @@
}
Rast_write_colors(wat_name, this_mapset, &colors);
}
+ seg_close(&watalt);
if (asp_flag) {
+ G_message(_("Closing flow direction map"));
cseg_write_cellfile(&asp, asp_name);
Rast_init_colors(&colors);
Rast_make_grey_scale_colors(&colors, 1, 8);
Rast_write_colors(asp_name, this_mapset, &colors);
}
cseg_close(&asp);
- /* visual ouput no longer needed */
- if (dis_flag) {
- if (bas_thres <= 0)
- bas_thres = 60;
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++) {
- dseg_get(&wat, &dvalue, r, c);
- if (dvalue < 0) {
- dvalue = 0;
- dseg_put(&wat, &dvalue, r, c);
- }
- else {
- bseg_get(&swale, &is_swale, r, c);
- if (is_swale) {
- dvalue = bas_thres;
- dseg_put(&wat, &dvalue, r, c);
- }
- }
- }
- }
- dseg_write_cellfile(&wat, dis_name);
- Rast_init_colors(&colors);
- Rast_make_rainbow_colors(&colors, 1, 120);
- Rast_write_colors(dis_name, this_mapset, &colors);
- }
- /* error in gislib.a
- Rast_free_colors(&colors);
- */
- dseg_close(&wat);
if (ls_flag) {
+ G_message(_("Closing LS map"));
dseg_write_cellfile(&l_s, ls_name);
dseg_close(&l_s);
}
- bseg_close(&swale);
if (sl_flag) {
+ G_message(_("Closing SL map"));
for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
for (c = 0; c < ncols; c++) {
dseg_get(&s_l, &dvalue, r, c);
if (dvalue > max_length)
dseg_put(&s_l, &max_length, r, c);
}
}
+ G_percent(r, nrows, 1); /* finish it */
dseg_write_cellfile(&s_l, sl_name);
}
if (sl_flag || ls_flag || sg_flag)
Modified: grass/trunk/raster/r.watershed/seg/close_maps2.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/close_maps2.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/close_maps2.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -1,5 +1,7 @@
#include "Gwater.h"
#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
int close_array_seg(void)
{
@@ -7,6 +9,7 @@
int incr, max, red, green, blue, rd, gr, bl, flag;
int c, r, map_fd;
CELL *cellrow, value;
+ char cvalue;
CSEG *theseg;
if (seg_flag || bas_flag || haf_flag) {
@@ -72,17 +75,22 @@
/* stream segments map */
if (seg_flag) {
+ G_message(_("Closing stream segments map"));
cellrow = (CELL *) G_malloc(ncols * sizeof(CELL));
map_fd = Rast_open_c_new(seg_name);
for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
Rast_set_c_null_value(cellrow, ncols); /* reset row to all NULL */
for (c = 0; c < ncols; c++) {
- bseg_get(&swale, &value, r, c);
- if (value)
+ /* bseg_get(&swale, &cvalue, r, c); */
+ /* if (cvalue) */
+ bseg_get(&bitflags, &cvalue, r, c);
+ if (FLAG_GET(cvalue, SWALEFLAG))
cseg_get(&bas, &(cellrow[c]), r, c);
}
Rast_put_row(map_fd, cellrow, CELL_TYPE);
}
+ G_percent(nrows, nrows, 1); /* finish it */
G_free(cellrow);
Rast_close(map_fd);
Rast_write_colors(seg_name, this_mapset, &colors);
@@ -90,12 +98,14 @@
/* basins map */
if (bas_flag) {
+ G_message(_("Closing basins map"));
cseg_write_cellfile(&bas, bas_name);
Rast_write_colors(bas_name, this_mapset, &colors);
}
/* half.basins map */
if (haf_flag) {
+ G_message(_("Closing half basins map"));
cseg_write_cellfile(&haf, haf_name);
Rast_write_colors(haf_name, this_mapset, &colors);
}
@@ -106,6 +116,9 @@
cseg_close(&bas);
if (arm_flag)
fclose(fp);
+
+ bseg_close(&bitflags);
+
close_maps();
return 0;
Modified: grass/trunk/raster/r.watershed/seg/cseg.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/cseg.h 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/cseg.h 2009-10-21 15:51:41 UTC (rev 39603)
@@ -42,13 +42,13 @@
int bseg_close(BSEG *);
/* bseg_get.c */
-int bseg_get(BSEG *, CELL *, int, int);
+int bseg_get(BSEG *, char *, int, int);
/* bseg_open.c */
int bseg_open(BSEG *, int, int, int);
/* bseg_put.c */
-int bseg_put(BSEG *, CELL *, int, int);
+int bseg_put(BSEG *, char *, int, int);
/* bseg_read.c */
int bseg_read_cell(BSEG *, char *, char *);
Modified: grass/trunk/raster/r.watershed/seg/cseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/cseg_write.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/cseg_write.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -20,6 +20,7 @@
buffer = Rast_allocate_c_buf();
segment_flush(&(cseg->seg));
for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
segment_get_row(&(cseg->seg), buffer, row);
if (Rast_put_row(map_fd, buffer, CELL_TYPE) < 0) {
G_free(buffer);
@@ -29,6 +30,7 @@
return -2;
}
}
+ G_percent(row, nrows, 1); /* finish it */
G_free(buffer);
Rast_close(map_fd);
return 0;
Modified: grass/trunk/raster/r.watershed/seg/def_basin.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/def_basin.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/def_basin.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -6,11 +6,11 @@
{
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;
for (;;) {
cseg_put(&bas, &basin_num, row, col);
- bseg_put(&swale, &one, row, col);
cseg_get(&asp, &asp_value, row, col);
if (asp_value < 0)
asp_value = -asp_value;
@@ -22,8 +22,8 @@
if (value < -1)
value = -value;
if (value == drain[rr][cc]) {
- bseg_get(&swale, &value, r, c);
- if (value) {
+ bseg_get(&bitflags, &cvalue, r, c);
+ if (FLAG_GET(cvalue, SWALEFLAG)) {
new_r[++ct] = r;
new_c[ct] = c;
}
Modified: grass/trunk/raster/r.watershed/seg/do_astar.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/do_astar.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -1,5 +1,3 @@
-#include <stdlib.h>
-#include <unistd.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include "Gwater.h"
@@ -7,48 +5,39 @@
int do_astar(void)
{
- POINT point;
int doer, count;
- SHORT upr, upc, r, c, ct_dir;
- CELL work_val, alt_val, alt_up, asp_up;
- DCELL wat_val;
- CELL in_val, drain_val;
- HEAP heap_pos;
+ SHORT upr, upc, r = -1, c = -1, ct_dir;
+ CELL alt_val, alt_up;
+ CELL asp_val;
+ char flag_value;
+ HEAP_PNT heap_p;
- /* double slope; */
+ G_message(_("SECTION 2: A* Search."));
- G_message(_("SECTION 2: A * Search."));
+ if (heap_size == 0)
+ G_fatal_error(_("No seeds for A* Search"));
+
+ G_debug(1, "heap size %d, points %d", heap_size, do_points);
count = 0;
- seg_get(&heap_index, (char *)&heap_pos, 0, 1);
- first_astar = heap_pos.point;
- /* A * Search: downhill path through all not masked cells */
- while (first_astar != -1) {
+ doer = do_points - 1;
+
+ /* A * Search
+ * determine downhill path through all not masked cells and
+ * set drainage direction */
+ while (heap_size > 0) {
G_percent(count++, do_points, 1);
- seg_get(&heap_index, (char *)&heap_pos, 0, 1);
- doer = heap_pos.point;
+ /* drop next point from heap */
+ heap_p = drop_pt();
- seg_get(&astar_pts, (char *)&point, 0, doer);
+ r = heap_p.pnt.r;
+ c = heap_p.pnt.c;
+ G_debug(3, "heap size %d, r %d, c %d", heap_size, r, c);
- /* drop astar_pts[doer] from heap */
- drop_pt();
+ alt_val = heap_p.ele;
- seg_get(&heap_index, (char *)&heap_pos, 0, 1);
- first_astar = heap_pos.point;
-
- point.nxt = first_cum;
- seg_put(&astar_pts, (char *)&point, 0, doer);
-
- first_cum = doer;
- r = point.r;
- c = point.c;
-
- bseg_put(&worked, &one, r, c);
- /* cseg_get(&alt, &alt_val, r, c); */
- alt_val = heap_pos.ele;
-
/* check all neighbours, breadth first search */
for (ct_dir = 0; ct_dir < sides; ct_dir++) {
/* get r, c (upr, upc) for this neighbour */
@@ -56,215 +45,166 @@
upc = c + nextdc[ct_dir];
/* check that upr, upc are within region */
if (upr >= 0 && upr < nrows && upc >= 0 && upc < ncols) {
- /* check if neighbour is in the list */
- /* if not, add as new point */
- bseg_get(&in_list, &in_val, upr, upc);
- if (in_val == 0) {
+ /* put neighbour in search list if not yet in */
+ bseg_get(&bitflags, &flag_value, upr, upc);
+ if (!FLAG_GET(flag_value, INLISTFLAG)) {
cseg_get(&alt, &alt_up, upr, upc);
- add_pt(upr, upc, alt_up, alt_val);
/* flow direction is set here */
- drain_val = drain[upr - r + 1][upc - c + 1];
- cseg_put(&asp, &drain_val, upr, upc);
+ asp_val = drain[upr - r + 1][upc - c + 1];
+ add_pt(upr, upc, alt_up, asp_val, 0);
+ cseg_put(&asp, &asp_val, upr, upc);
}
- else {
- /* check if neighbour has not been worked on,
- * update values for asp and wat */
- bseg_get(&worked, &work_val, upr, upc);
- if (!work_val) {
- cseg_get(&asp, &asp_up, upr, upc);
- if (asp_up < 0) {
- drain_val = drain[upr - r + 1][upc - c + 1];
- cseg_put(&asp, &drain_val, upr, upc);
- dseg_get(&wat, &wat_val, r, c);
- if (wat_val > 0) {
- wat_val = -wat_val;
- dseg_put(&wat, &wat_val, r, c);
- }
- /* cseg_get(&alt, &alt_up, upr, upc); */
- /* replace(upr, upc, r, c); */ /* alt_up used to be */
- /* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
- dseg_put (&slp, &slope, upr, upc); */
- }
+ /* in list, not worked. is it edge ? */
+ else if (!FLAG_GET(flag_value, WORKEDFLAG) && FLAG_GET(flag_value, EDGEFLAG)){
+ cseg_get(&asp, &asp_val, upr, upc);
+ if (asp_val < 0) {
+ /* adjust flow direction for edge cell */
+ asp_val = drain[upr - r + 1][upc - c + 1];
+ cseg_put(&asp, &asp_val, upr, upc);
+ heap_p.pnt.guessed = 1;
}
}
}
}
+ cseg_get(&asp, &asp_val, r, c);
+ heap_p.pnt.asp = asp_val;
+ seg_put(&astar_pts, (char *)&heap_p.pnt, 0, doer);
+ doer--;
+ bseg_get(&bitflags, &flag_value, r, c);
+ FLAG_SET(flag_value, WORKEDFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
}
+ if (doer != -1)
+ G_fatal_error(_("bug in A* Search: doer %d heap size %d count %d"), doer, heap_size, count);
- if (mfd == 0)
- bseg_close(&worked);
+ seg_close(&search_heap);
- bseg_close(&in_list);
- seg_close(&heap_index);
-
G_percent(count, do_points, 1); /* finish it */
+
return 0;
}
-/* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
+/* compare two heap points */
+/* return 1 if a < b else 0 */
+int cmp_pnt(HEAP_PNT *a, HEAP_PNT *b)
{
- POINT point;
- HEAP heap_pos;
+ if (a->ele < b->ele)
+ return 1;
+ else if (a->ele == b->ele) {
+ if (a->added < b->added)
+ return 1;
+ }
+ return 0;
+}
- /* double slp_value; */
+/* add point routine for min heap */
+int add_pt(SHORT r, SHORT c, CELL ele, char asp, int is_edge)
+{
+ HEAP_PNT heap_p;
+ char flag_value;
- /* slp_value = get_slope(r, c, downr, downc, ele, downe);
- dseg_put (&slp, &slp_value, r, c); */
- bseg_put(&in_list, &one, r, c);
+ bseg_get(&bitflags, &flag_value, r, c);
+ if (is_edge) {
+ FLAG_SET(flag_value, EDGEFLAG);
+ heap_p.pnt.guessed = 1;
+ }
+ else {
+ heap_p.pnt.guessed = 0;
+ }
+ FLAG_SET(flag_value, INLISTFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
/* add point to next free position */
heap_size++;
- if (heap_size > do_points)
- G_fatal_error(_("heapsize too large"));
+ heap_p.added = nxt_avail_pt;
+ heap_p.ele = ele;
+ heap_p.pnt.r = r;
+ heap_p.pnt.c = c;
+ heap_p.pnt.asp = asp;
- heap_pos.point = nxt_avail_pt;
- heap_pos.ele = ele;
- seg_put(&heap_index, (char *)&heap_pos, 0, heap_size);
-
- point.r = r;
- point.c = c;
-/* point.downr = downr;
- point.downc = downc; */
-
- seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
-
- /* cseg_put(&pnt_index, &nxt_avail_pt, r, c); */
-
nxt_avail_pt++;
/* sift up: move new point towards top of heap */
+ sift_up(heap_size, heap_p);
- sift_up(heap_size, ele);
-
- return 0;
+ return 1;
}
/* drop point routine for min heap */
-int drop_pt(void)
+HEAP_PNT drop_pt(void)
{
int child, childr, parent;
- int childp; /* , childrp */
- CELL ele; /* , eler */
int i;
- HEAP heap_pos;
+ HEAP_PNT child_p, childr_p, last_p, root_p;
- if (heap_size == 1) {
- parent = -1;
- heap_pos.point = -1;
- heap_pos.ele = 0;
- seg_put(&heap_index, (char *)&heap_pos, 0, 1);
- heap_size = 0;
- return 0;
- }
+ seg_get(&search_heap, (char *)&last_p, 0, heap_size);
+ seg_get(&search_heap, (char *)&root_p, 0, 1);
- /* start with heap root */
+ /* sift down */
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) {
+ while ((child = GET_CHILD(parent)) < heap_size) {
/* select child with lower ele, if both are equal, older child
* older child is older startpoint for flow path, important */
- seg_get(&heap_index, (char *)&heap_pos, 0, child);
- childp = heap_pos.point;
- ele = heap_pos.ele;
+
+ seg_get(&search_heap, (char *)&child_p, 0, child);
if (child < heap_size) {
childr = child + 1;
- i = child + 3;
- while (childr <= heap_size && childr < i) {
- seg_get(&heap_index, (char *)&heap_pos, 0, childr);
- /* childrp = heap_pos.point;
- eler = heap_pos.ele; */
- if (heap_pos.ele < ele) {
+ 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;
- childp = heap_pos.point;
- ele = heap_pos.ele;
+ child_p = childr_p;
}
- /* make sure we get the oldest child */
- else if (ele == heap_pos.ele && childp > heap_pos.point) {
- child = childr;
- childp = heap_pos.point;
- }
childr++;
- /* i++; */
}
}
+ if (cmp_pnt(&last_p, &child_p)) {
+ break;
+ }
+
/* move hole down */
- heap_pos.point = childp;
- heap_pos.ele = ele;
- seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+ seg_put(&search_heap, (char *)&child_p, 0, parent);
parent = child;
-
}
- /* hole is in lowest layer, move to heap end */
- if (parent < heap_size) {
- seg_get(&heap_index, (char *)&heap_pos, 0, heap_size);
- seg_put(&heap_index, (char *)&heap_pos, 0, parent);
+ seg_put(&search_heap, (char *)&last_p, 0, parent);
- ele = heap_pos.ele;
-
- /* sift up last swapped point, only necessary if hole moved to heap end */
- sift_up(parent, ele);
-
- }
-
/* the actual drop */
heap_size--;
- return 0;
-
+ return root_p;
}
/* standard sift-up routine for d-ary min heap */
-int sift_up(int start, CELL ele)
+int sift_up(int start, HEAP_PNT child_p)
{
- int parent, child, childp; /* parentp, */
- /* CELL elep; */
- HEAP heap_pos;
+ int parent, child; /* parentp, */
+ HEAP_PNT heap_p;
child = start;
- seg_get(&heap_index, (char *)&heap_pos, 0, child);
- childp = heap_pos.point;
while (child > 1) {
parent = GET_PARENT(child);
- seg_get(&heap_index, (char *)&heap_pos, 0, parent);
- /* parentp = heap_pos.point;
- elep = heap_pos.ele; */
+ seg_get(&search_heap, (char *)&heap_p, 0, parent);
- /* parent ele higher */
- if (heap_pos.ele > ele) {
-
- /* push parent point down */
- seg_put(&heap_index, (char *)&heap_pos, 0, child);
+ /* push parent point down if child is smaller */
+ if (cmp_pnt(&child_p, &heap_p)) {
+ seg_put(&search_heap, (char *)&heap_p, 0, child);
child = parent;
-
}
- /* same ele, but parent is younger */
- else if (heap_pos.ele == ele && heap_pos.point > childp) {
-
- /* push parent point down */
- seg_put(&heap_index, (char *)&heap_pos, 0, child);
- child = parent;
-
- }
+ /* no more sifting up, found slot for child */
else
- /* no more sifting up, found new slot for child */
break;
}
- /* no more sifting up, found new slot for child */
- if (child < start) {
- heap_pos.point = childp;
- heap_pos.ele = ele;
- seg_put(&heap_index, (char *)&heap_pos, 0, child);
- }
+ /* add child to heap */
+ seg_put(&search_heap, (char *)&child_p, 0, child);
+
return 0;
}
@@ -283,24 +223,3 @@
return (MIN_SLOPE);
return (slope);
}
-
-int replace(SHORT upr, SHORT upc, SHORT r, SHORT c)
-{ /* ele was in there */
- /* CELL ele; */
- int now;
- POINT point;
-
- now = 0;
-
- /* cseg_get(&pnt_index, &now, upr, upc); */
- seg_get(&astar_pts, (char *)&point, 0, now);
- if (point.r != upr || point.c != upc) {
- G_warning("pnt_index incorrect!");
- return 1;
- }
-/* point.downr = r;
- point.downc = c; */
- seg_put(&astar_pts, (char *)&point, 0, now);
-
- return 0;
-}
Modified: grass/trunk/raster/r.watershed/seg/do_astar.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_astar.h 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/do_astar.h 2009-10-21 15:51:41 UTC (rev 39603)
@@ -1,7 +1,11 @@
#ifndef __DO_ASTAR_H__
#define __DO_ASTAR_H__
+#define GET_PARENT(c) ((int) ((((c) - 2) >> 2) + 1))
+#define GET_CHILD(p) ((int) (((p) << 2) - 2))
+
+/*
#define GET_PARENT(c) ((int) (((c) - 2) / 3 + 1))
#define GET_CHILD(p) ((int) ((p) * 3 - 1))
-
+*/
#endif /* __DO_ASTAR_H__ */
Modified: grass/trunk/raster/r.watershed/seg/do_cum.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/do_cum.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/do_cum.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -8,39 +8,48 @@
int do_cum(void)
{
SHORT r, c, dr, dc;
- CELL is_swale, asp_val, asp_val_down;
+ CELL asp_val, asp_val_down;
+ char is_swale, this_flag_value, flag_value;
DCELL value, valued;
POINT point;
- int killer, threshold, count;
+ 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 };
+ WAT_ALT wa, wadown;
G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
- count = 0;
if (bas_thres <= 0)
threshold = 60;
else
threshold = bas_thres;
- while (first_cum != -1) {
- G_percent(count++, do_points, 2);
- killer = first_cum;
+ for (killer = 0; killer < do_points; killer++) {
+ G_percent(killer, do_points, 1);
seg_get(&astar_pts, (char *)&point, 0, killer);
- first_cum = point.nxt;
r = point.r;
c = point.c;
- cseg_get(&asp, &asp_val, r, c);
+ asp_val = point.asp;
if (asp_val) {
dr = r + asp_r[ABS(asp_val)];
dc = c + asp_c[ABS(asp_val)];
}
else
dr = dc = -1;
- if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
- dseg_get(&wat, &value, r, c);
- if (ABS(value) >= threshold)
- bseg_put(&swale, &one, r, c);
- dseg_get(&wat, &valued, dr, dc);
+
+ bseg_get(&bitflags, &this_flag_value, r, c);
+ FLAG_UNSET(this_flag_value, WORKEDFLAG);
+
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {
+ /* TODO: do not distribute flow along edges, this causes artifacts */
+ seg_get(&watalt, (char *)&wa, r, c);
+ value = wa.wat;
+ is_swale = FLAG_GET(this_flag_value, SWALEFLAG);
+ if (fabs(value) >= threshold && !is_swale) {
+ is_swale = 1;
+ FLAG_SET(this_flag_value, SWALEFLAG);
+ }
+ seg_get(&watalt, (char *)&wadown, dr, dc);
+ valued = wadown.wat;
if (value > 0) {
if (valued > 0)
valued += value;
@@ -53,8 +62,8 @@
else
valued = value - valued;
}
- dseg_put(&wat, &valued, dr, dc);
- bseg_get(&swale, &is_swale, r, c);
+ wadown.wat = valued;
+ seg_put(&watalt, (char *)&wadown, dr, dc);
/* update asp for depression */
if (is_swale && pit_flag) {
cseg_get(&asp, &asp_val_down, dr, dc);
@@ -63,18 +72,23 @@
cseg_put(&asp, &asp_val, r, c);
}
}
- if (is_swale || ABS(valued) >= threshold) {
- bseg_put(&swale, &one, dr, dc);
+ if (is_swale || fabs(valued) >= threshold) {
+ bseg_get(&bitflags, &flag_value, dr, dc);
+ FLAG_SET(flag_value, SWALEFLAG);
+ bseg_put(&bitflags, &flag_value, dr, dc);
+ is_swale = 1;
}
else {
if (er_flag && !is_swale)
slope_length(r, c, dr, dc);
}
}
+ bseg_put(&bitflags, &this_flag_value, r, c);
}
+ G_percent(do_points, do_points, 1); /* finish it */
+
seg_close(&astar_pts);
- G_percent(count, do_points, 1); /* finish it */
return 0;
}
@@ -103,19 +117,20 @@
int do_cum_mfd(void)
{
int r, c, dr, dc;
- CELL is_swale;
DCELL value, valued, *wat_nbr;
POINT point;
- int killer, threshold, count;
+ WAT_ALT wa;
+ int killer, threshold;
/* MFD */
int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
double *dist_to_nbr, *weight, sum_weight, max_weight;
int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
double dx, dy;
- CELL ele, ele_nbr, asp_val, asp_val2, cvalue, *worked_nbr;
+ CELL ele, *ele_nbr, asp_val;
double prop, max_acc;
- int workedon, edge;
+ int edge;
+ char workedon, *flag_nbr, this_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 };
@@ -139,42 +154,48 @@
dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
}
- /* reset worked, takes time... */
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++) {
- bseg_put(&worked, &zero, r, c);
- }
- }
-
- worked_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+ flag_nbr = (char *)G_malloc(sides * sizeof(char));
wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+ ele_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
workedon = 0;
- count = 0;
if (bas_thres <= 0)
threshold = 60;
else
threshold = bas_thres;
- while (first_cum != -1) {
- G_percent(count++, do_points, 2);
- killer = first_cum;
+ for (killer = 0; killer < do_points; killer++) {
+ G_percent(killer, do_points, 1);
seg_get(&astar_pts, (char *)&point, 0, killer);
- first_cum = point.nxt;
r = point.r;
c = point.c;
- cseg_get(&asp, &asp_val, r, c);
+ asp_val = point.asp;
if (asp_val) {
dr = r + asp_r[ABS(asp_val)];
dc = c + asp_c[ABS(asp_val)];
}
else
dr = dc = -1;
- if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
- /* dc = point.downc; */
- dseg_get(&wat, &value, r, c);
+ seg_get(&watalt, (char *)&wa, r, c);
+ value = wa.wat;
+ if (point.guessed && value > 0) {
+ value = -value;
+ wa.wat = value;
+ seg_put(&watalt, (char *)&wa, r, c);
+ }
+
+ bseg_get(&bitflags, &this_flag_value, r, c);
+ FLAG_UNSET(this_flag_value, WORKEDFLAG);
+
+ if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {
+ /* do not distribute flow along edges, this causes artifacts */
+ if (FLAG_GET(this_flag_value, EDGEFLAG)) {
+ bseg_put(&bitflags, &this_flag_value, r, c);
+ continue;
+ }
+
r_max = dr;
c_max = dc;
@@ -186,7 +207,7 @@
stream_cells = 0;
swale_cells = 0;
astar_not_set = 1;
- cseg_get(&alt, &ele, r, c);
+ ele = wa.ele;
is_null = 0;
edge = 0;
/* this loop is needed to get the sum of weights */
@@ -195,33 +216,34 @@
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
weight[ct_dir] = -1;
- /* check that neighbour is within region */
+
+ wat_nbr[ct_dir] = 0;
+ ele_nbr[ct_dir] = 0;
+
+ /* WORKEDFLAG has been set during A* Search
+ * reversed meaning here: 0 = done, 1 = not yet done */
+ FLAG_UNSET(flag_nbr[ct_dir], WORKEDFLAG);
+
+ /* check if neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols) {
- /* check for swale or stream cells */
- bseg_get(&swale, &is_swale, r_nbr, c_nbr);
- if (is_swale)
- swale_cells++;
- dseg_get(&wat, &valued, r_nbr, c_nbr);
- wat_nbr[ct_dir] = valued;
- if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
- stream_cells++;
+ bseg_get(&bitflags, &flag_nbr[ct_dir], r_nbr, c_nbr);
+ if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
- 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 = Rast_is_c_null_value(&ele_nbr);
- edge = is_null;
- if (!is_null && ele_nbr <= ele) {
- if (ele_nbr < ele) {
+ seg_get(&watalt, (char *)&wa, r_nbr, c_nbr);
+ wat_nbr[ct_dir] = wa.wat;
+ ele_nbr[ct_dir] = wa.ele;
+
+ edge = is_null = FLAG_GET(flag_nbr[ct_dir], NULLFLAG);
+ if (!is_null && ele_nbr[ct_dir] <= ele) {
+ if (ele_nbr[ct_dir] < ele) {
weight[ct_dir] =
mfd_pow(((ele -
- ele_nbr) / dist_to_nbr[ct_dir]),
+ ele_nbr[ct_dir]) / dist_to_nbr[ct_dir]),
c_fac);
}
- if (ele_nbr == ele) {
+ if (ele_nbr[ct_dir] == ele) {
weight[ct_dir] =
mfd_pow((0.5 / dist_to_nbr[ct_dir]),
c_fac);
@@ -248,11 +270,8 @@
}
/* do not distribute flow along edges, this causes artifacts */
if (edge) {
- bseg_get(&swale, &is_swale, r, c);
- if (is_swale && asp_val > 0) {
- asp_val = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
- cseg_put(&asp, &asp_val, r, c);
- }
+ G_debug(3, "edge: should not get to here");
+ bseg_put(&bitflags, &this_flag_value, r, c);
continue;
}
@@ -279,12 +298,12 @@
r_nbr = r + nextdr[ct_dir];
c_nbr = c + nextdc[ct_dir];
- /* check that neighbour is within region */
+ /* check if neighbour is within region */
if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
c_nbr < ncols && weight[ct_dir] > -0.5) {
- /* bseg_get(&worked, &is_worked, r_nbr, c_nbr); */
- if (worked_nbr[ct_dir] == 0) {
+ if (FLAG_GET(flag_nbr[ct_dir], WORKEDFLAG)) {
+
weight[ct_dir] = weight[ct_dir] / sum_weight;
/* check everything sums up to 1.0 */
prop += weight[ct_dir];
@@ -302,7 +321,9 @@
wat_nbr[ct_dir] = value * weight[ct_dir] - wat_nbr[ct_dir];
}
valued = wat_nbr[ct_dir];
- dseg_put(&wat, &valued, r_nbr, c_nbr);
+ wa.wat = valued;
+ wa.ele = ele_nbr[ct_dir];
+ seg_put(&watalt, (char *)&wa, r_nbr, c_nbr);
/* get main drainage direction */
if (ABS(wat_nbr[ct_dir]) >= max_acc) {
@@ -321,10 +342,18 @@
G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
prop);
}
+
+ /* adjust main drainage direction to A* path if possible */
+ if (fabs(wat_nbr[np_side]) >= max_acc) {
+ max_acc = fabs(wat_nbr[ct_dir]);
+ r_max = dr;
+ c_max = dc;
+ }
}
-
- if (mfd_cells < 2) {
- dseg_get(&wat, &valued, dr, dc);
+ /* SFD-like accumulation */
+ else {
+ seg_get(&watalt, (char *)&wa, dr, dc);
+ valued = wa.wat;
if (value > 0) {
if (valued > 0)
valued += value;
@@ -337,58 +366,37 @@
else
valued = value - valued;
}
- dseg_put(&wat, &valued, dr, dc);
+ wa.wat = valued;
+ seg_put(&watalt, (char *)&wa, dr, dc);
}
- /* update asp */
- if (dr != r_max || dc != c_max) {
- asp_val2 = drain[r - r_max + 1][c - c_max + 1];
- /* cseg_get(&asp, &asp_val, r, c); */
- if (asp_val < 0)
- asp_val2 = -asp_val2;
- cseg_put(&asp, &asp_val2, r, c);
-
- }
- /* update asp for depression */
- bseg_get(&swale, &is_swale, r, c);
- if (is_swale && pit_flag) {
- cseg_get(&asp, &asp_val2, r_max, c_max);
- if (asp_val > 0 && asp_val2 == 0) {
- asp_val = -asp_val;
- cseg_put(&asp, &asp_val, r, c);
+ if (!st_flag) {
+ /* update asp */
+ if (dr != r_max || dc != c_max) {
+ if (asp_val > 0) {
+ asp_val = drain[r - r_max + 1][c - c_max + 1];
+ cseg_put(&asp, &asp_val, r, c);
+ }
}
}
- /* start new stream */
- value = ABS(value) + 0.5;
- if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
- swale_cells < 1) {
- bseg_put(&swale, &one, r, c);
- is_swale = 1;
- }
- /* continue stream */
- if (is_swale) {
- bseg_put(&swale, &one, r_max, c_max);
- }
- else {
- if (er_flag && !is_swale)
- slope_length(r, c, r_max, c_max);
- }
- bseg_put(&worked, &one, r, c);
}
+ bseg_put(&bitflags, &this_flag_value, r, c);
}
- G_percent(count, do_points, 1); /* finish it */
+ G_percent(do_points, do_points, 1); /* finish it */
+
if (workedon)
G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
workedon, do_points);
- seg_close(&astar_pts);
-
- bseg_close(&worked);
+ if (!st_flag) {
+ seg_close(&astar_pts);
+ }
G_free(dist_to_nbr);
G_free(weight);
G_free(wat_nbr);
- G_free(worked_nbr);
+ G_free(ele_nbr);
+ G_free(flag_nbr);
return 0;
}
Modified: grass/trunk/raster/r.watershed/seg/dseg_write.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/dseg_write.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/dseg_write.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -21,6 +21,7 @@
dbuffer = Rast_allocate_d_buf();
segment_flush(&(dseg->seg));
for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 1);
segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
if (Rast_put_row(map_fd, dbuffer, DCELL_TYPE) < 0) {
G_free(dbuffer);
@@ -30,6 +31,7 @@
return -2;
}
}
+ G_percent(row, nrows, 1); /* finish it */
G_free(dbuffer);
Rast_close(map_fd);
return 0;
Modified: grass/trunk/raster/r.watershed/seg/find_pour.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/find_pour.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/find_pour.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -4,11 +4,12 @@
{
int row, col;
double easting, northing, stream_length;
- CELL old_elev, basin_num, value, is_swale;
+ CELL old_elev, basin_num, value;
+ char is_swale;
basin_num = 0;
for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 3);
+ G_percent(row, nrows, 1);
northing = window.north - (row + .5) * window.ns_res;
for (col = 0; col < ncols; col++) {
/* cseg_get (&wat, &value, row, col);
@@ -17,7 +18,9 @@
value = -value;
} */
cseg_get(&asp, &value, row, col);
- bseg_get(&swale, &is_swale, row, col);
+ bseg_get(&bitflags, &is_swale, row, col);
+ is_swale = FLAG_GET(is_swale, SWALEFLAG);
+ /* bseg_get(&swale, &is_swale, row, col); */
if (value < 0 && is_swale > 0) {
basin_num += 2;
cseg_get(&alt, &old_elev, row, col);
Added: grass/trunk/raster/r.watershed/seg/flag.h
===================================================================
--- grass/trunk/raster/r.watershed/seg/flag.h (rev 0)
+++ grass/trunk/raster/r.watershed/seg/flag.h 2009-10-21 15:51:41 UTC (rev 39603)
@@ -0,0 +1,28 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+/* a set of routines that allow the programmer to "flag" cells in a
+ * raster map. A flag is of type unsigned char, i.e. 8 bits can be set.
+ *
+ * int flag_set(flag, bitno)
+ * sets the flag at position bitno to one.
+ *
+ * int flag_unset(flag, bitno)
+ * sets the flag at position bitno to zero.
+ *
+ * int flag_get(flag, bitno)
+ * checks if the flag is set at postion bitno.
+ *
+ * Examples:
+ * set flag at position 0: FLAG_SET(flag, 0)
+ * unset (clear) flag at position 7: FLAG_UNSET(flag, 7)
+ * check flag at position 5: is_set_at_5 = FLAG_GET(flag, 5)
+ */
+
+#define FLAG_SET(flag,bitno) ((flag) |= (1 << (bitno)))
+
+#define FLAG_UNSET(flag,bitno) ((flag) &= ~(1 << (bitno)))
+
+#define FLAG_GET(flag,bitno) ((flag) & (1 << (bitno)))
+
+#endif /* __FLAG_H__ */
Modified: grass/trunk/raster/r.watershed/seg/init_vars.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/init_vars.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/init_vars.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -1,6 +1,5 @@
#include <stdlib.h>
#include <string.h>
-#include <unistd.h>
#include "Gwater.h"
#include <grass/gis.h>
#include <grass/raster.h>
@@ -11,24 +10,27 @@
int init_vars(int argc, char *argv[])
{
SHORT r, c;
- int fd, num_cseg_total, num_open_segs;
- int seg_rows, seg_cols;
- double segs_mb;
+ int ele_fd, fd, wat_fd;
+ int seg_rows, seg_cols, num_cseg_total, num_open_segs, num_open_array_segs;
/* int page_block, num_cseg; */
int max_bytes;
- CELL *buf, alt_value, asp_value, worked_value, block_value;
+ CELL *buf, alt_value, asp_value, block_value;
+ char worked_value, flag_value;
DCELL wat_value;
DCELL dvalue;
+ WAT_ALT wa;
char MASK_flag;
- void *elebuf, *ptr;
- int ele_map_type;
- size_t ele_size;
+ void *elebuf, *ptr, *watbuf, *watptr;
+ int ele_map_type, wat_map_type;
+ size_t ele_size, wat_size;
+ CELL cellone = 1;
+ int ct_dir, r_nbr, c_nbr;
G_gisinit(argv[0]);
ele_flag = wat_flag = asp_flag = pit_flag = run_flag = ril_flag = 0;
- ob_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
- zero = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
+ st_flag = bas_flag = seg_flag = haf_flag = arm_flag = dis_flag = 0;
+ ob_flag = sl_flag = sg_flag = ls_flag = er_flag = bas_thres = 0;
nxt_avail_pt = 0;
/* dep_flag = 0; */
max_length = d_zero = 0.0;
@@ -42,6 +44,7 @@
abs_acc = 0;
ele_scale = 1;
segs_mb = 300;
+ /* scan options */
for (r = 1; r < argc; r++) {
if (sscanf(argv[r], "el=%[^\n]", ele_name) == 1)
ele_flag++;
@@ -93,6 +96,7 @@
else
usage(argv[0]);
}
+ /* check options */
if (mfd == 1 && (c_fac < 1 || c_fac > 10)) {
G_fatal_error("Convergence factor must be between 1 and 10.");
}
@@ -107,17 +111,26 @@
)
usage(argv[0]);
tot_parts = 4;
- if (ls_flag || sg_flag)
+ if (sl_flag || sg_flag || ls_flag)
+ er_flag = 1;
+ /* do RUSLE */
+ if (er_flag)
tot_parts++;
- if (bas_thres > 0)
+ /* 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++;
G_message(_("SECTION 1 beginning: Initiating Variables. %d sections total."),
tot_parts);
this_mapset = G_mapset();
- if (sl_flag || sg_flag || ls_flag)
- er_flag = 1;
/* for sd factor
if (dep_flag) {
if (sscanf (dep_name, "%lf", &dep_slope) != 1) {
@@ -140,8 +153,8 @@
diag *= 0.5;
/* segment parameters: one size fits all. Fine tune? */
- /* Segment rows and cols: 200 */
- /* 1 segment open for all rasters: 1.34 MB */
+ /* Segment rows and cols: 128 */
+ /* 1 segment open for all data structures: 0.122 MB */
seg_rows = SROW;
seg_cols = SCOL;
@@ -151,7 +164,7 @@
G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
}
- num_open_segs = segs_mb / 1.34;
+ num_open_segs = segs_mb / 0.122;
G_debug(1, "segs MB: %.0f", segs_mb);
G_debug(1, "region rows: %d", nrows);
@@ -174,44 +187,74 @@
num_open_segs = num_cseg_total;
G_debug(1, " open segments after adjusting:\t%d", num_open_segs);
- cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
- cseg_read_cell(&alt, ele_name, "");
if (er_flag) {
cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
cseg_read_cell(&r_h, ele_name, "");
}
/* read elevation input and mark NULL/masked cells */
- bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
- bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
- G_verbose_message("Checking for masked and NULL cells in input elevation <%s>", ele_name);
+ /* scattered access: alt, watalt, listflag, asp */
+ cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
+ seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs, sizeof(WAT_ALT));
+ bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs);
+ cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+
/* open elevation input */
- fd = Rast_open_old(ele_name, "");
- if (fd < 0) {
+ ele_fd = Rast_open_old(ele_name, "");
+ if (ele_fd < 0) {
G_fatal_error(_("unable to open elevation map layer"));
}
- ele_map_type = Rast_get_map_type(fd);
+ ele_map_type = Rast_get_map_type(ele_fd);
ele_size = Rast_cell_size(ele_map_type);
elebuf = Rast_allocate_buf(ele_map_type);
if (ele_map_type == FCELL_TYPE || ele_map_type == DCELL_TYPE)
ele_scale = 1000; /* should be enough to do the trick */
+ /* initial flow accumulation */
+ if (run_flag) {
+ wat_fd = Rast_open_old(run_name, "");
+ if (wat_fd < 0) {
+ G_fatal_error(_("unable to open accumulation map layer"));
+ }
+
+ wat_map_type = Rast_get_map_type(ele_fd);
+ wat_size = Rast_cell_size(ele_map_type);
+ watbuf = Rast_allocate_buf(ele_map_type);
+ }
+ else {
+ watbuf = watptr = NULL;
+ wat_fd = wat_size = wat_map_type = -1;
+ }
+
/* read elevation input and mark NULL/masked cells */
+ G_message("SECTION 1a: Mark masked and NULL cells");
MASK_flag = 0;
do_points = nrows * ncols;
for (r = 0; r < nrows; r++) {
- Rast_get_row(fd, elebuf, r, ele_map_type);
+ G_percent(r, nrows, 1);
+ Rast_get_row(ele_fd, elebuf, r, ele_map_type);
ptr = elebuf;
+
+ if (run_flag) {
+ Rast_get_row(wat_fd, watbuf, r, wat_map_type);
+ watptr = watbuf;
+ }
+
for (c = 0; c < ncols; c++) {
+ flag_value = 0;
+
/* check for masked and NULL cells */
if (Rast_is_null_value(ptr, ele_map_type)) {
- bseg_put(&worked, &one, r, c);
- bseg_put(&in_list, &one, r, c);
+ FLAG_SET(flag_value, NULLFLAG);
+ FLAG_SET(flag_value, INLISTFLAG);
+ FLAG_SET(flag_value, WORKEDFLAG);
Rast_set_c_null_value(&alt_value, 1);
+ /* flow accumulation */
+ Rast_set_d_null_value(&wat_value, 1);
do_points--;
}
else {
@@ -228,50 +271,58 @@
dvalue *= ele_scale;
alt_value = ele_round(dvalue);
}
+
+ /* flow accumulation */
+ if (run_flag) {
+ if (Rast_is_null_value(watptr, wat_map_type)) {
+ wat_value = 0; /* ok ? */
+ }
+ else {
+ if (wat_map_type == CELL_TYPE) {
+ wat_value = *((CELL *)ptr);
+ }
+ else if (wat_map_type == FCELL_TYPE) {
+ wat_value = *((FCELL *)ptr);
+ }
+ else if (ele_map_type == DCELL_TYPE) {
+ wat_value = *((DCELL *)ptr);
+ }
+ }
+ }
+ else {
+ wat_value = 1;
+ }
+
}
cseg_put(&alt, &alt_value, r, c);
+ wa.wat = wat_value;
+ wa.ele = alt_value;
+ seg_put(&watalt, (char *)&wa, r, c);
+
+ if (flag_value)
+ bseg_put(&bitflags, &flag_value, r, c);
+
if (er_flag) {
cseg_put(&r_h, &alt_value, r, c);
}
ptr = G_incr_void_ptr(ptr, ele_size);
+ if (run_flag) {
+ watptr = G_incr_void_ptr(watptr, wat_size);
+ }
}
}
- Rast_close(fd);
+ G_percent(nrows, nrows, 1); /* finish it */
+ Rast_close(ele_fd);
G_free(elebuf);
- if (do_points < nrows * ncols)
- MASK_flag = 1;
- /* initial flow accumulation */
- dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
if (run_flag) {
- dseg_read_cell(&wat, run_name, "");
- if (MASK_flag) {
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++) {
- bseg_get(&worked, &worked_value, r, c);
- if (worked_value)
- dseg_put(&wat, &d_zero, r, c);
- }
- }
- }
+ Rast_close(wat_fd);
+ G_free(watbuf);
}
- else {
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++)
- if (MASK_flag) {
- bseg_get(&worked, &worked_value, r, c);
- if (worked_value)
- dseg_put(&wat, &d_zero, r, c);
- else
- dseg_put(&wat, &d_one, r, c);
- }
- else {
- if (-1 == dseg_put(&wat, &d_one, r, c))
- exit(EXIT_FAILURE);
- }
- }
- }
- cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+
+ if (do_points < nrows * ncols)
+ MASK_flag = 1;
+
/* depression: drainage direction will be set to zero later */
if (pit_flag) {
fd = Rast_open_old(pit_name, "");
@@ -280,28 +331,23 @@
}
buf = Rast_allocate_c_buf();
for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
Rast_get_c_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
asp_value = buf[c];
if (!Rast_is_c_null_value(&asp_value) && asp_value) {
- cseg_put(&asp, &one, r, c);
+ cseg_put(&asp, &cellone, r, c);
+ bseg_get(&bitflags, &flag_value, r, c);
+ FLAG_SET(flag_value, PITFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
}
- else {
- cseg_put(&asp, &zero, r, c);
- }
}
}
+ G_percent(nrows, nrows, 1); /* finish it */
Rast_close(fd);
G_free(buf);
}
- else {
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++)
- if (-1 == cseg_put(&asp, &zero, r, c))
- exit(EXIT_FAILURE);
- }
- }
- bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
+
if (ob_flag) {
fd = Rast_open_old(ob_name, "");
if (fd < 0) {
@@ -309,26 +355,22 @@
}
buf = Rast_allocate_c_buf();
for (r = 0; r < nrows; r++) {
+ G_percent(r, nrows, 1);
Rast_get_c_row(fd, buf, r);
for (c = 0; c < ncols; c++) {
block_value = buf[c];
if (!Rast_is_c_null_value(&block_value) && block_value) {
- bseg_put(&swale, &one, r, c);
+ bseg_get(&bitflags, &flag_value, r, c);
+ FLAG_SET(flag_value, RUSLEBLOCKFLAG);
+ bseg_put(&bitflags, &flag_value, r, c);
}
- else {
- bseg_put(&swale, &zero, r, c);
- }
}
}
+ G_percent(nrows, nrows, 1); /* finish it */
Rast_close(fd);
G_free(buf);
}
- else {
- for (r = 0; r < nrows; r++) {
- for (c = 0; c < ncols; c++)
- bseg_put(&swale, &zero, r, c);
- }
- }
+
if (ril_flag) {
dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
dseg_read_cell(&ril, ril_name, "");
@@ -337,7 +379,6 @@
/* dseg_open(&slp, SROW, SCOL, num_open_segs); */
/* RUSLE: LS and/or S factor */
-
if (er_flag) {
dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
}
@@ -346,16 +387,48 @@
if (ls_flag)
dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
- seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols * 2,
- num_open_segs / 2, sizeof(POINT));
+ G_debug(1, "open segments for A* points");
+ /* next smaller power of 2 */
+ seg_cols = (int) pow(2, (int)(log(num_open_segs / 8) / log(2)));
+ /* n cols in segment */
+ seg_cols *= seg_rows * seg_rows;
+ /* n segments in row */
+ num_cseg_total = do_points / seg_cols;
+ if (do_points % seg_cols > 0)
+ num_cseg_total++;
+ /* no need to have more segments open than exist */
+ if ((num_open_segs * seg_rows * seg_rows) / (seg_cols * 2) > num_cseg_total)
+ num_open_array_segs = num_cseg_total;
+ else
+ num_open_array_segs = (num_open_segs * seg_rows * seg_rows) / (seg_cols * 2);
+
+ G_debug(2, "A* points open segments %d, target 4", num_open_array_segs);
+
+ seg_open(&astar_pts, 1, do_points, 1, seg_cols, num_open_array_segs,
+ sizeof(POINT));
- /* heap_index will track astar_pts in ternary min-heap */
- /* heap_index is one-based */
- seg_open(&heap_index, 1, do_points + 1, 1, seg_cols * num_open_segs * seg_rows / 10,
- 10, sizeof(HEAP));
+ /* one-based d-ary search_heap with astar_pts */
+ G_debug(1, "open segments for A* search heap");
+ /* next smaller power of 2 */
+ seg_cols = (int) pow(2, (int)(log(num_open_segs / 16) / log(2)));
+ /* n cols in segment */
+ seg_cols *= seg_rows * seg_rows;
+ /* n segments in row */
+ num_cseg_total = (do_points + 1) / seg_cols;
+ if ((do_points + 1) % seg_cols > 0)
+ num_cseg_total++;
+ /* no need to have more segments open than exist */
+ if ((num_open_segs * seg_rows * seg_rows) / (seg_cols * 2) > num_cseg_total)
+ num_open_array_segs = num_cseg_total;
+ else
+ num_open_array_segs = (num_open_segs * seg_rows * seg_rows) / (seg_cols * 2);
- G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
+ G_debug(2, "A* search heap open segments %d, target 8", num_open_array_segs);
+ seg_open(&search_heap, 1, do_points + 1, 1, seg_cols,
+ num_open_array_segs, sizeof(HEAP_PNT));
+ G_message(_("SECTION 1b: Determining Offmap Flow."));
+
/* heap is empty */
heap_size = 0;
@@ -363,29 +436,32 @@
if (MASK_flag) {
for (r = 0; r < nrows; r++) {
- G_percent(r, nrows, 2);
+ G_percent(r, nrows, 1);
for (c = 0; c < ncols; c++) {
- bseg_get(&worked, &worked_value, r, c);
- if (worked_value) {
- dseg_put(&wat, &d_zero, r, c);
- }
- else {
+ bseg_get(&bitflags, &flag_value, r, c);
+ if (!FLAG_GET(flag_value, NULLFLAG)) {
if (er_flag)
dseg_put(&s_l, &half_res, r, c);
cseg_get(&asp, &asp_value, r, c);
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
- dseg_get(&wat, &wat_value, r, c);
+ /* dseg_get(&wat, &wat_value, r, c); */
+ seg_get(&watalt, (char *)&wa, r, c);
+ wat_value = wa.wat;
if (wat_value > 0) {
wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
+ /* dseg_put(&wat, &wat_value, r, c); */
+ wa.wat = wat_value;
+ seg_put(&watalt, (char *)&wa, r, c);
}
/* set depression */
if (asp_value) {
asp_value = 0;
if (wat_value < 0) {
wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
+ /* dseg_put(&wat, &wat_value, r, c); */
+ wa.wat = wat_value;
+ seg_put(&watalt, (char *)&wa, r, c);
}
}
else if (r == 0)
@@ -398,139 +474,59 @@
asp_value = -8;
if (-1 == cseg_put(&asp, &asp_value, r, c))
exit(EXIT_FAILURE);
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
+ /* cseg_get(&alt, &alt_value, r, c); */
+ alt_value = wa.ele;
+ add_pt(r, c, alt_value, asp_value, 1);
}
- else if (!bseg_get(&worked, &worked_value, r - 1, c)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -2;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
+ else {
+ seg_get(&watalt, (char *)&wa, r, c);
+ 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];
+
+ bseg_get(&bitflags, &worked_value, r_nbr, c_nbr);
+ if (FLAG_GET(worked_value, NULLFLAG)) {
+ asp_value = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
+ add_pt(r, c, wa.ele, asp_value, 1);
+ cseg_put(&asp, &asp_value, r, c);
+ wat_value = wa.wat;
+ if (wat_value > 0) {
+ wa.wat = -wat_value;
+ seg_put(&watalt, (char *)&wa, r, c);
+ }
+ break;
+ }
}
}
- else if (!bseg_get(&worked, &worked_value, r + 1, c)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -6;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (!bseg_get(&worked, &worked_value, r, c - 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -4;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (!bseg_get(&worked, &worked_value, r, c + 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -8;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (sides == 8 &&
- !bseg_get(&worked, &worked_value, r - 1, c - 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -3;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (sides == 8 &&
- !bseg_get(&worked, &worked_value, r - 1, c + 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -1;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (sides == 8 &&
- !bseg_get(&worked, &worked_value, r + 1, c - 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -5;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else if (sides == 8 &&
- !bseg_get(&worked, &worked_value, r + 1, c + 1)
- && worked_value != 0) {
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
- asp_value = -7;
- cseg_put(&asp, &asp_value, r, c);
- dseg_get(&wat, &wat_value, r, c);
- if (wat_value > 0) {
- wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
- }
- }
- else {
- bseg_put(&in_list, &zero, r, c);
- /* dseg_put(&slp, &dzero, r, c); */
- }
- }
- }
+ } /* end non-NULL cell */
+ } /* end column */
}
}
else {
for (r = 0; r < nrows; r++) {
- G_percent(r, nrows, 2);
+ G_percent(r, nrows, 1);
for (c = 0; c < ncols; c++) {
- bseg_put(&worked, &zero, r, c);
+ /* bseg_put(&worked, &zero, r, c); */
if (er_flag)
dseg_put(&s_l, &half_res, r, c);
cseg_get(&asp, &asp_value, r, c);
if (r == 0 || c == 0 || r == nrows - 1 ||
c == ncols - 1 || asp_value != 0) {
- dseg_get(&wat, &wat_value, r, c);
+ seg_get(&watalt, (char *)&wa, r, c);
+ wat_value = wa.wat;
if (wat_value > 0) {
wat_value = -wat_value;
- if (-1 == dseg_put(&wat, &wat_value, r, c))
- exit(EXIT_FAILURE);
+ wa.wat = wat_value;
+ seg_put(&watalt, (char *)&wa, r, c);
}
/* set depression */
if (asp_value) {
asp_value = 0;
if (wat_value < 0) {
wat_value = -wat_value;
- dseg_put(&wat, &wat_value, r, c);
+ wa.wat = wat_value;
+ seg_put(&watalt, (char *)&wa, r, c);
}
}
else if (r == 0)
@@ -543,13 +539,9 @@
asp_value = -8;
if (-1 == cseg_put(&asp, &asp_value, r, c))
exit(EXIT_FAILURE);
- cseg_get(&alt, &alt_value, r, c);
- add_pt(r, c, alt_value, alt_value);
+ /* cseg_get(&alt, &alt_value, r, c); */
+ add_pt(r, c, wa.ele, asp_value, 1);
}
- else {
- bseg_put(&in_list, &zero, r, c);
- /* dseg_put(&slp, &dzero, r, c); */
- }
}
}
}
Modified: grass/trunk/raster/r.watershed/seg/main.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/main.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/main.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -26,18 +26,19 @@
struct Cell_head window;
int mfd, c_fac, abs_acc, ele_scale;
-SSEG heap_index;
+SSEG search_heap;
int heap_size;
int first_astar, first_cum, nxt_avail_pt, total_cells, do_points;
SHORT nrows, ncols;
double half_res, diag, max_length, dep_slope;
int bas_thres, tot_parts;
SSEG astar_pts;
-BSEG worked, in_list, s_b, swale;
+BSEG bitflags, s_b;
CSEG dis, alt, asp, bas, haf, r_h, dep;
-DSEG wat;
+SSEG watalt;
DSEG slp, s_l, s_g, l_s, ril;
-CELL one, zero;
+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} };
@@ -64,8 +65,10 @@
int main(int argc, char *argv[])
{
+ int num_open_segs;
+
+ zero = 0;
one = 1;
- zero = 0;
d_zero = 0.0;
d_one = 1.0;
init_vars(argc, argv);
@@ -76,11 +79,14 @@
else {
do_cum();
}
+ if (st_flag && mfd) {
+ do_stream();
+ }
if (sg_flag || ls_flag) {
sg_factor();
}
- if (bas_thres <= 0) {
+ if (!seg_flag && !bas_flag && !haf_flag) {
G_message(_("SECTION %d: Closing Maps."), tot_parts);
close_maps();
}
@@ -88,8 +94,12 @@
if (arm_flag) {
fp = fopen(arm_name, "w");
}
- cseg_open(&bas, SROW, SCOL, 4);
- cseg_open(&haf, SROW, SCOL, 4);
+ num_open_segs = segs_mb / 0.122;
+ if (num_open_segs > (ncols / SCOL + 1) * (nrows / SROW + 1)) {
+ num_open_segs = (ncols / SCOL + 1) * (nrows / SROW + 1);
+ }
+ cseg_open(&bas, SROW, SCOL, num_open_segs);
+ cseg_open(&haf, SROW, SCOL, num_open_segs);
G_message(_("SECTION %d: Watershed determination."), tot_parts - 1);
find_pourpts();
G_message(_("SECTION %d: Closing Maps."), tot_parts);
Modified: grass/trunk/raster/r.watershed/seg/no_stream.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/no_stream.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/no_stream.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -9,6 +9,7 @@
CELL downdir, asp_value, hih_ele, new_ele, aspect, value;
DCELL dvalue, max_drain; /* flow acc is now DCELL */
SHORT updir, riteflag, leftflag, thisdir;
+ WAT_ALT wa;
while (1) {
max_drain = -1;
@@ -18,10 +19,11 @@
cseg_get(&asp, &aspect, r, c);
if (aspect == drain[rr][cc]) {
- dseg_get(&wat, &dvalue, r, c);
+ seg_get(&watalt, (char *)&wa, r, c);
+ dvalue = wa.wat;
if (dvalue < 0)
dvalue = -dvalue;
- if ((dvalue - max_drain) > 5E-8f) { /* floating point comparison problem workaround */
+ if ((dvalue - max_drain) > 5E-8f) {
uprow = r;
upcol = c;
max_drain = dvalue;
Modified: grass/trunk/raster/r.watershed/seg/sg_factor.c
===================================================================
--- grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-10-21 15:49:13 UTC (rev 39602)
+++ grass/trunk/raster/r.watershed/seg/sg_factor.c 2009-10-21 15:51:41 UTC (rev 39603)
@@ -8,11 +8,16 @@
int r, c;
CELL downer, low_elev, hih_elev;
double height, length, S, sin_theta;
+ char flag_value;
- G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
+ G_message(_("SECTION 5: RUSLE LS and/or S factor determination."));
for (r = nrows - 1; r >= 0; r--) {
G_percent(nrows - r, nrows, 3);
for (c = ncols - 1; c >= 0; c--) {
+ bseg_get(&bitflags, &flag_value, r, c);
+ if (FLAG_GET(flag_value, NULLFLAG))
+ continue;
+
cseg_get(&alt, &low_elev, r, c);
cseg_get(&r_h, &hih_elev, r, c);
dseg_get(&s_l, &length, r, c);
More information about the grass-commit
mailing list