← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/r.watershed/ram/do_cum.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.watershed/ram/do_cum.c
1
1
#include "Gwater.h"
2
2
#include <grass/gis.h>
 
3
#include <grass/raster.h>
3
4
#include <grass/glocale.h>
4
5
 
 
6
double get_dist(double *dist_to_nbr, double *contour)
 
7
{
 
8
    int ct_dir, r_nbr, c_nbr;
 
9
    double dx, dy, ns_res, ew_res;
 
10
 
 
11
    if (G_projection() == PROJECTION_LL) {
 
12
        double ew_dist1, ew_dist2, ew_dist3;
 
13
        double ns_dist1, ns_dist2, ns_dist3;
 
14
 
 
15
        G_begin_distance_calculations();
 
16
 
 
17
        /* EW Dist at North edge */
 
18
        ew_dist1 = G_distance(window.east, window.north,
 
19
                              window.west, window.north);
 
20
        /* EW Dist at Center */
 
21
        ew_dist2 = G_distance(window.east, (window.north + window.south) / 2.,
 
22
                              window.west, (window.north + window.south) / 2.);
 
23
        /* EW Dist at South Edge */
 
24
        ew_dist3 = G_distance(window.east, window.south,
 
25
                              window.west, window.south);
 
26
        /* NS Dist at East edge */
 
27
        ns_dist1 = G_distance(window.east, window.north,
 
28
                              window.east, window.south);
 
29
        /* NS Dist at Center */
 
30
        ns_dist2 = G_distance((window.west + window.east) / 2., window.north,
 
31
                              (window.west + window.east) / 2., window.south);
 
32
        /* NS Dist at West edge */
 
33
        ns_dist3 = G_distance(window.west, window.north,
 
34
                              window.west, window.south);
 
35
 
 
36
        ew_res = (ew_dist1 + ew_dist2 + ew_dist3) / (3 * window.cols);
 
37
        ns_res = (ns_dist1 + ns_dist2 + ns_dist3) / (3 * window.rows);
 
38
    }
 
39
    else {
 
40
        ns_res = window.ns_res;
 
41
        ew_res = window.ew_res;
 
42
    }
 
43
    
 
44
    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 
45
        /* get r, c (r_nbr, c_nbr) for neighbours */
 
46
        r_nbr = nextdr[ct_dir];
 
47
        c_nbr = nextdc[ct_dir];
 
48
        /* account for rare cases when ns_res != ew_res */
 
49
        dy = ABS(r_nbr) * ns_res;
 
50
        dx = ABS(c_nbr) * ew_res;
 
51
        if (ct_dir < 4)
 
52
            dist_to_nbr[ct_dir] = (dx + dy) * ele_scale;
 
53
        else
 
54
            dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy) * ele_scale;
 
55
    }
 
56
    /* Quinn et al. 1991:
 
57
     * ns contour: ew_res / 2
 
58
     * ew contour: ns_res / 2
 
59
     * diag contour: sqrt(ns_res * nsres / 4 + ew_res * ew_res / 4) / 2
 
60
     *             = sqrt(ns_res * nsres + ew_res * ew_res) / 4
 
61
     * if ns_res == ew_res:
 
62
     *             sqrt(2 * (res * res) / 4 = res * 0.354
 
63
     *
 
64
     * contour lengths "have been subjectively chosen", 
 
65
     * no justification why the diagonal contour is shorter
 
66
     * BUT: if the diag contour is a bit shorter than the cardinal contour,
 
67
     * this is further enhancing the correction for diagonal flow bias
 
68
     * diagonal slope is already corrected for longer distance
 
69
     * smaller slope and shorter contour length have the same effect:
 
70
     * higher TCI
 
71
     */
 
72
    if (sides == 8) {
 
73
        /* contours are sides of an octagon, irregular if ns_res != ew_res
 
74
         * ideally: arc lengths of an ellipse */
 
75
        contour[0] = contour[1] = tan(atan(ew_res / ns_res) / 2.) * ns_res;
 
76
        contour[2] = contour[3] = tan(atan(ns_res / ew_res) / 2.) * ew_res;
 
77
        G_debug(1, "ns contour: %.4f", contour[0]);
 
78
        G_debug(1, "ew contour: %.4f", contour[2]);
 
79
        contour[4] = (ew_res - contour[0]);
 
80
        contour[5] = (ns_res - contour[2]);
 
81
        contour[7] = sqrt(contour[4] * contour[4] + contour[5] * contour[5]) / 2.;
 
82
        G_debug(1, "diag contour: %.4f", contour[7]);
 
83
        contour[4] = contour[5] = contour[6] = contour[7];
 
84
    }
 
85
    else {
 
86
        /* contours are sides of a rectangle */
 
87
        contour[0] = contour[1] = ew_res;
 
88
        contour[2] = contour[3] = ns_res;
 
89
    }
 
90
    return ew_res * ns_res;
 
91
}
 
92
 
 
93
double get_slope_tci(CELL ele, CELL down_ele, double dist)
 
94
{
 
95
    if (down_ele >= ele)
 
96
        return 0.5 / dist;
 
97
    else
 
98
        return (double)(ele - down_ele) / dist;
 
99
}
5
100
 
6
101
int do_cum(void)
7
102
{
8
 
    SHORT r, c, dr, dc;
9
 
    CELL is_swale, value, valued, aspect;
10
 
    int killer, threshold, count;
11
 
    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
12
 
    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
13
 
    int this_index, down_index;
 
103
    int r, c, dr, dc;
 
104
    int r_nbr, c_nbr, ct_dir, np_side, edge;
 
105
    CELL is_swale, aspect, ele_nbr;
 
106
    DCELL value, valued;
 
107
    int killer, threshold;
 
108
    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
 
109
    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
110
    int this_index, down_index, nbr_index;
 
111
    double *dist_to_nbr, *contour;
 
112
    double tci_div, cell_size;
14
113
 
15
114
    G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
16
115
 
17
 
    count = 0;
 
116
    /* distances to neighbours, contour lengths */
 
117
    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
 
118
    contour = (double *)G_malloc(sides * sizeof(double));
 
119
 
 
120
    cell_size = get_dist(dist_to_nbr, contour);
 
121
 
18
122
    if (bas_thres <= 0)
19
123
        threshold = 60;
20
124
    else
28
132
            dr = r + asp_r[ABS(aspect)];
29
133
            dc = c + asp_c[ABS(aspect)];
30
134
        }
 
135
        /* skip user-defined depressions */
31
136
        else
32
137
            dr = dc = -1;
33
138
        if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
34
139
            down_index = SEG_INDEX(wat_seg, dr, dc);
35
140
            value = wat[this_index];
36
 
            if ((int)(ABS(value) + 0.5) >= threshold)
 
141
            if (fabs(value) >= threshold)
37
142
                FLAG_SET(swale, r, c);
38
143
            valued = wat[down_index];
 
144
 
 
145
            edge = 0;
 
146
            np_side = -1;
 
147
            for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 
148
                /* get r, c (r_nbr, c_nbr) for neighbours */
 
149
                r_nbr = r + nextdr[ct_dir];
 
150
                c_nbr = c + nextdc[ct_dir];
 
151
 
 
152
                if (dr == r_nbr && dc == c_nbr)
 
153
                    np_side = ct_dir;
 
154
 
 
155
                /* check that neighbour is within region */
 
156
                if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 
157
                    c_nbr < ncols) {
 
158
 
 
159
                    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
 
160
                    ele_nbr = alt[nbr_index];
 
161
                    if (Rast_is_c_null_value(&ele_nbr))
 
162
                        edge = 1;
 
163
                }
 
164
                else
 
165
                    edge = 1;
 
166
                if (edge)
 
167
                    break;
 
168
            }
 
169
            /* do not distribute flow along edges, this causes artifacts */
 
170
            if (edge) {
 
171
                is_swale = FLAG_GET(swale, r, c);
 
172
                if (is_swale && aspect > 0) {
 
173
                    aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
 
174
                    asp[this_index] = aspect;
 
175
                }
 
176
                continue;
 
177
            }
 
178
 
39
179
            if (value > 0) {
40
180
                if (valued > 0)
41
181
                    valued += value;
49
189
                    valued = value - valued;
50
190
            }
51
191
            wat[down_index] = valued;
52
 
            valued = ABS(valued) + 0.5;
 
192
 
 
193
            /* topographic wetness index ln(a / tan(beta)) */
 
194
            if (tci_flag) {
 
195
                tci_div = contour[np_side] * 
 
196
                       get_slope_tci(alt[this_index], alt[down_index],
 
197
                                     dist_to_nbr[np_side]);
 
198
                tci[this_index] = log((fabs(wat[this_index]) * cell_size) / tci_div);
 
199
            }
 
200
 
53
201
            is_swale = FLAG_GET(swale, r, c);
54
 
            /* update asp for depression */
55
 
            if (is_swale && pit_flag) {
56
 
                if (aspect > 0 && asp[down_index] == 0) {
57
 
                    aspect = -aspect;
58
 
                    asp[this_index] = aspect;
59
 
                }
60
 
            }
61
 
            if (is_swale || ((int)valued) >= threshold) {
 
202
            if (is_swale || fabs(valued) >= threshold) {
62
203
                FLAG_SET(swale, dr, dc);
63
204
            }
64
205
            else {
92
233
 * before depressions/obstacles and gracefull flow divergence after 
93
234
 * depressions/obstacles
94
235
 * 
 
236
 * Topographic Convergence Index (TCI)
 
237
 * tendency of water to accumulate at a given point considering 
 
238
 * the gravitational forces to move the water accumulated at that 
 
239
 * point further downstream
 
240
 *
 
241
 * after Quinn et al. (1991), modified and adapted for the modified 
 
242
 * Holmgren MFD algorithm
 
243
 * TCI: specific catchment area divided by tangens of slope
 
244
 * specific catchment area: total catchment area divided by contour line
 
245
 * TCI for D8:     A / (L * tanb)
 
246
 * TCI for MFD:    A / (SUM(L_i) * (SUM(tanb_i * weight_i) / SUM(weight_i))
 
247
 * 
 
248
 * A: total catchment area
 
249
 * L_i: contour length towards i_th cell
 
250
 * tanb_i: slope = tan(b) towards i_th cell
 
251
 * weight_i: weight for flow distribution towards i_th cell
95
252
 * ************************************/
96
253
 
97
254
int do_cum_mfd(void)
98
255
{
99
256
    int r, c, dr, dc;
100
257
    CELL is_swale;
101
 
    DCELL value, valued;
102
 
    int killer, threshold, count;
 
258
    DCELL value, valued, tci_div, sum_contour, cell_size;
 
259
    int killer, threshold;
103
260
 
104
261
    /* MFD */
105
262
    int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
106
 
    double *dist_to_nbr, *weight, sum_weight, max_weight;
 
263
    double *dist_to_nbr, *contour, *weight, sum_weight, max_weight;
107
264
    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
108
 
    double dx, dy;
109
265
    CELL ele, ele_nbr, aspect, is_worked;
110
 
    double prop, max_acc;
 
266
    double prop, max_val;
111
267
    int workedon, edge, flat;
112
 
    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
113
 
    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
268
    int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
 
269
    int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
114
270
    int this_index, down_index, nbr_index;
115
 
 
116
 
    G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
 
271
    
 
272
    G_message(_("SECTION 3a: Accumulating Surface Flow with MFD."));
117
273
    G_debug(1, "MFD convergence factor set to %d.", c_fac);
118
274
 
119
 
    /* distances to neighbours */
 
275
    /* distances to neighbours, weights, contour lengths */
120
276
    dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
121
277
    weight = (double *)G_malloc(sides * sizeof(double));
122
 
 
123
 
    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
124
 
        /* get r, c (r_nbr, c_nbr) for neighbours */
125
 
        r_nbr = nextdr[ct_dir];
126
 
        c_nbr = nextdc[ct_dir];
127
 
        /* account for rare cases when ns_res != ew_res */
128
 
        dy = ABS(r_nbr) * window.ns_res;
129
 
        dx = ABS(c_nbr) * window.ew_res;
130
 
        if (ct_dir < 4)
131
 
            dist_to_nbr[ct_dir] = dx + dy;
132
 
        else
133
 
            dist_to_nbr[ct_dir] = sqrt(dx * dx + dy * dy);
134
 
    }
 
278
    contour = (double *)G_malloc(sides * sizeof(double));
 
279
    
 
280
    cell_size = get_dist(dist_to_nbr, contour);
135
281
 
136
282
    flag_clear_all(worked);
137
283
    workedon = 0;
138
284
 
139
 
    count = 0;
140
285
    if (bas_thres <= 0)
141
286
        threshold = 60;
142
287
    else
146
291
        G_percent(killer, do_points, 1);
147
292
        this_index = astar_pts[killer];
148
293
        seg_index_rc(alt_seg, this_index, &r, &c);
 
294
        FLAG_SET(worked, r, c);
149
295
        aspect = asp[this_index];
150
296
        if (aspect) {
151
297
            dr = r + asp_r[ABS(aspect)];
156
302
        if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
157
303
            value = wat[this_index];
158
304
            down_index = SEG_INDEX(wat_seg, dr, dc);
159
 
            valued = wat[down_index];
160
 
 
161
 
            r_max = dr;
162
 
            c_max = dc;
163
305
 
164
306
            /* get weights */
165
307
            max_weight = 0;
166
308
            sum_weight = 0;
167
309
            np_side = -1;
168
310
            mfd_cells = 0;
169
 
            stream_cells = 0;
170
 
            swale_cells = 0;
171
311
            astar_not_set = 1;
172
312
            ele = alt[this_index];
173
313
            is_null = 0;
174
314
            edge = 0;
175
 
            flat = 1;
176
315
            /* this loop is needed to get the sum of weights */
177
316
            for (ct_dir = 0; ct_dir < sides; ct_dir++) {
178
317
                /* get r, c (r_nbr, c_nbr) for neighbours */
179
318
                r_nbr = r + nextdr[ct_dir];
180
319
                c_nbr = c + nextdc[ct_dir];
181
320
                weight[ct_dir] = -1;
 
321
 
 
322
                if (dr == r_nbr && dc == c_nbr)
 
323
                    np_side = ct_dir;
 
324
 
182
325
                /* check that neighbour is within region */
183
326
                if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
184
327
                    c_nbr < ncols) {
185
328
 
186
329
                    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
187
330
 
188
 
                    /* check for swale or stream cells */
189
 
                    is_swale = FLAG_GET(swale, r_nbr, c_nbr);
190
 
                    if (is_swale)
191
 
                        swale_cells++;
192
331
                    valued = wat[nbr_index];
193
332
                    ele_nbr = alt[nbr_index];
194
 
                    if ((ABS(valued) + 0.5) >= threshold &&
195
 
                        ct_dir != np_side && ele_nbr > ele)
196
 
                        stream_cells++;
197
333
 
198
334
                    is_worked = FLAG_GET(worked, r_nbr, c_nbr);
199
335
                    if (is_worked == 0) {
200
 
                        is_null = G_is_c_null_value(&ele_nbr);
 
336
                        is_null = Rast_is_c_null_value(&ele_nbr);
201
337
                        edge = is_null;
202
 
                        if (ele_nbr != ele)
203
 
                            flat = 0;
204
338
                        if (!is_null && ele_nbr <= ele) {
205
339
                            if (ele_nbr < ele) {
206
340
                                weight[ct_dir] =
225
359
                            }
226
360
                        }
227
361
                    }
228
 
                    if (dr == r_nbr && dc == c_nbr)
229
 
                        np_side = ct_dir;
230
362
                }
231
363
                else
232
364
                    edge = 1;
235
367
            }
236
368
            /* do not distribute flow along edges, this causes artifacts */
237
369
            if (edge) {
238
 
                is_swale = FLAG_GET(swale, r, c);
239
 
                if (is_swale && aspect > 0) {
240
 
                    aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
241
 
                    asp[this_index] = aspect;
242
 
                }
243
370
                continue;
244
371
            }
245
372
 
257
384
            }
258
385
 
259
386
            /* set flow accumulation for neighbours */
260
 
            max_acc = -1;
 
387
            max_val = -1;
 
388
            tci_div = sum_contour = 0.;
261
389
 
262
390
            if (mfd_cells > 1) {
263
391
                prop = 0.0;
264
392
                for (ct_dir = 0; ct_dir < sides; ct_dir++) {
265
 
                    /* get r, c (r_nbr, c_nbr) for neighbours */
266
393
                    r_nbr = r + nextdr[ct_dir];
267
394
                    c_nbr = c + nextdc[ct_dir];
268
395
 
274
401
 
275
402
                            nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
276
403
 
 
404
                            if (tci_flag) {
 
405
                                sum_contour += contour[ct_dir];
 
406
                                tci_div += get_slope_tci(ele, alt[nbr_index],
 
407
                                                         dist_to_nbr[ct_dir]) *
 
408
                                           weight[ct_dir];
 
409
                            }
 
410
 
277
411
                            weight[ct_dir] = weight[ct_dir] / sum_weight;
278
 
                            /* check everything sums up to 1.0 */
 
412
                            /* check everything adds up to 1.0 */
279
413
                            prop += weight[ct_dir];
280
414
 
281
415
                            valued = wat[nbr_index];
292
426
                                    valued = value * weight[ct_dir] - valued;
293
427
                            }
294
428
                            wat[nbr_index] = valued;
295
 
 
296
 
                            /* get main drainage direction */
297
 
                            if (ABS(valued) >= max_acc) {
298
 
                                max_acc = ABS(valued);
299
 
                                r_max = r_nbr;
300
 
                                c_max = c_nbr;
301
 
                            }
302
429
                        }
303
430
                        else if (ct_dir == np_side) {
304
431
                            /* check for consistency with A * path */
310
437
                    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
311
438
                              prop);
312
439
                }
 
440
                if (tci_flag)
 
441
                    tci_div /= sum_weight;
313
442
            }
314
443
 
315
444
            if (mfd_cells < 2) {
327
456
                        valued = value - valued;
328
457
                }
329
458
                wat[down_index] = valued;
330
 
            }
331
 
 
 
459
 
 
460
                if (tci_flag) {
 
461
                    sum_contour = contour[np_side];
 
462
                    tci_div = get_slope_tci(ele, alt[down_index],
 
463
                                            dist_to_nbr[np_side]);
 
464
                }
 
465
            }
 
466
            /* topographic wetness index ln(a / tan(beta)) */
 
467
            if (tci_flag) {
 
468
                tci[this_index] = log((fabs(wat[this_index]) * cell_size) /
 
469
                                      (sum_contour * tci_div));
 
470
            }
 
471
        }
 
472
    }
 
473
    if (workedon)
 
474
        G_warning(_n("MFD: A * path already processed when distributing flow: %d of %d cell", 
 
475
        "MFD: A * path already processed when distributing flow: %d of %d cells", 
 
476
        do_points),
 
477
                  workedon, do_points);
 
478
 
 
479
 
 
480
    G_message(_("SECTION 3b: Adjusting drainage directions."));
 
481
 
 
482
    for (killer = 1; killer <= do_points; killer++) {
 
483
        G_percent(killer, do_points, 1);
 
484
        this_index = astar_pts[killer];
 
485
        seg_index_rc(alt_seg, this_index, &r, &c);
 
486
        FLAG_UNSET(worked, r, c);
 
487
        aspect = asp[this_index];
 
488
        if (aspect) {
 
489
            dr = r + asp_r[ABS(aspect)];
 
490
            dc = c + asp_c[ABS(aspect)];
 
491
        }
 
492
        else
 
493
            dr = dc = -1;
 
494
        if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
 
495
            value = wat[this_index];
 
496
            down_index = SEG_INDEX(wat_seg, dr, dc);
 
497
 
 
498
            r_max = dr;
 
499
            c_max = dc;
 
500
 
 
501
            /* get max flow accumulation */
 
502
            max_val = -1;
 
503
            stream_cells = 0;
 
504
            swale_cells = 0;
 
505
            ele = alt[this_index];
 
506
            is_null = 0;
 
507
            edge = 0;
 
508
            flat = 1;
 
509
 
 
510
            for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 
511
                /* get r, c (r_nbr, c_nbr) for neighbours */
 
512
                r_nbr = r + nextdr[ct_dir];
 
513
                c_nbr = c + nextdc[ct_dir];
 
514
 
 
515
                /* check that neighbour is within region */
 
516
                if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 
517
                    c_nbr < ncols) {
 
518
 
 
519
                    nbr_index = SEG_INDEX(wat_seg, r_nbr, c_nbr);
 
520
 
 
521
                    /* check for swale or stream cells */
 
522
                    is_swale = FLAG_GET(swale, r_nbr, c_nbr);
 
523
                    if (is_swale)
 
524
                        swale_cells++;
 
525
                    valued = wat[nbr_index];
 
526
                    ele_nbr = alt[nbr_index];
 
527
                    edge = Rast_is_c_null_value(&ele_nbr);
 
528
                    if ((ABS(valued) + 0.5) >= threshold  &&
 
529
                        ele_nbr > ele)
 
530
                        stream_cells++;
 
531
 
 
532
                    is_worked = !(FLAG_GET(worked, r_nbr, c_nbr));
 
533
                    if (is_worked == 0) {
 
534
                        if (ele_nbr != ele)
 
535
                            flat = 0;
 
536
                        is_null = Rast_is_c_null_value(&ele_nbr);
 
537
                        edge = is_null;
 
538
                        if (!is_null && ABS(valued) > max_val) {
 
539
                            max_val = ABS(valued);
 
540
                            r_max = r_nbr;
 
541
                            c_max = c_nbr;
 
542
                        }
 
543
                    }
 
544
                }
 
545
                else
 
546
                    edge = 1;
 
547
                if (edge)
 
548
                    break;
 
549
            }
 
550
            /* do not distribute flow along edges, this causes artifacts */
 
551
            if (edge) {
 
552
                is_swale = FLAG_GET(swale, r, c);
 
553
                if (is_swale && aspect > 0) {
 
554
                    aspect = -1 * drain[r - r_nbr + 1][c - c_nbr + 1];
 
555
                    asp[this_index] = aspect;
 
556
                }
 
557
                continue;
 
558
            }
 
559
            
332
560
            /* update asp */
333
561
            if (dr != r_max || dc != c_max) {
334
562
                aspect = drain[r - r_max + 1][c - c_max + 1];
337
565
                asp[this_index] = aspect;
338
566
            }
339
567
            is_swale = FLAG_GET(swale, r, c);
340
 
            /* update asp for depression */
341
 
            if (is_swale && pit_flag) {
342
 
                if (aspect > 0 && asp[SEG_INDEX(asp_seg, r_max, c_max)] == 0) {
343
 
                    aspect = -aspect;
344
 
                    asp[this_index] = aspect;
345
 
                }
346
 
            }
347
568
            /* start new stream */
348
569
            value = ABS(value) + 0.5;
349
570
            if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
359
580
                if (er_flag && !is_swale)
360
581
                    slope_length(r, c, r_max, c_max);
361
582
            }
362
 
            FLAG_SET(worked, r, c);
363
583
        }
364
584
    }
365
 
    if (workedon)
366
 
        G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
367
 
                  workedon, do_points);
368
585
 
369
586
    G_free(astar_pts);
370
587
 
378
595
 
379
596
double mfd_pow(double base, int exp)
380
597
{
381
 
    int x;
382
 
    double result;
383
 
 
384
 
    result = base;
385
 
    if (exp == 1)
386
 
        return result;
387
 
 
388
 
    for (x = 2; x <= exp; x++) {
 
598
    int i;
 
599
    double result = base;
 
600
 
 
601
    for (i = 2; i <= exp; i++) {
389
602
        result *= base;
390
603
    }
391
604
    return result;