source: grass/branches/releasebranch_7_0/raster/r.stream.extract/main.c

Last change on this file was 63708, checked in by martinl, 10 years ago

r.stream.extract: stream_rast -> stream_raster & stream_vect -> stream_vector (#2409)

(merge r63706-7 from trunk)
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 14.1 KB
Line 
1
2/****************************************************************************
3 *
4 * MODULE: r.stream.extract
5 * AUTHOR(S): Markus Metz <markus.metz.giswork gmail.com>
6 * PURPOSE: Hydrological analysis
7 * Extracts stream networks from accumulation raster with
8 * given threshold
9 * COPYRIGHT: (C) 1999-2014 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16#include <stdlib.h>
17#include <string.h>
18#include <float.h>
19#include <math.h>
20#include <grass/raster.h>
21#include <grass/glocale.h>
22#include "local_proto.h"
23
24/* global variables */
25int nrows, ncols;
26GW_LARGE_INT n_search_points, n_points, nxt_avail_pt;
27GW_LARGE_INT heap_size;
28GW_LARGE_INT n_stream_nodes, n_alloc_nodes;
29POINT *outlets;
30struct snode *stream_node;
31GW_LARGE_INT n_outlets, n_alloc_outlets;
32char drain[3][3] = { {7, 6, 5}, {8, 0, 4}, {1, 2, 3} };
33char sides;
34int c_fac;
35int ele_scale;
36int have_depressions;
37
38SSEG search_heap;
39SSEG astar_pts;
40SSEG watalt, aspflag;
41/* BSEG bitflags, asp; */
42CSEG stream;
43
44CELL *astar_order;
45
46int main(int argc, char *argv[])
47{
48 struct
49 {
50 struct Option *ele, *acc, *depression;
51 struct Option *threshold, *d8cut;
52 struct Option *mont_exp;
53 struct Option *min_stream_length;
54 struct Option *memory;
55 } input;
56 struct
57 {
58 struct Option *stream_rast;
59 struct Option *stream_vect;
60 struct Option *dir_rast;
61 } output;
62 struct GModule *module;
63 int ele_fd, acc_fd, depr_fd;
64 double threshold, d8cut, mont_exp;
65 int min_stream_length = 0, memory;
66 int seg_cols, seg_rows;
67 double seg2kb;
68 int num_open_segs, num_open_array_segs, num_seg_total;
69 double memory_divisor, heap_mem, disk_space;
70 const char *mapset;
71
72 G_gisinit(argv[0]);
73
74 /* Set description */
75 module = G_define_module();
76 G_add_keyword(_("raster"));
77 G_add_keyword(_("hydrology"));
78 G_add_keyword(_("stream network"));
79 module->description = _("Performs stream network extraction.");
80
81 input.ele = G_define_standard_option(G_OPT_R_ELEV);
82
83 input.acc = G_define_standard_option(G_OPT_R_INPUT);
84 input.acc->key = "accumulation";
85 input.acc->label = _("Name of input accumulation raster map");
86 input.acc->required = NO;
87 input.acc->description =
88 _("Stream extraction will use provided accumulation instead of calculating it anew");
89 input.acc->guisection = _("Input maps");
90
91 input.depression = G_define_standard_option(G_OPT_R_INPUT);
92 input.depression->key = "depression";
93 input.depression->label = _("Name of input raster map with real depressions");
94 input.depression->required = NO;
95 input.depression->description =
96 _("Streams will not be routed out of real depressions");
97 input.depression->guisection = _("Input maps");
98
99 input.threshold = G_define_option();
100 input.threshold->key = "threshold";
101 input.threshold->label = _("Minimum flow accumulation for streams");
102 input.threshold->description = _("Must be > 0");
103 input.threshold->required = YES;
104 input.threshold->type = TYPE_DOUBLE;
105
106 input.d8cut = G_define_option();
107 input.d8cut->key = "d8cut";
108 input.d8cut->label = _("Use SFD above this threshold");
109 input.d8cut->description =
110 _("If accumulation is larger than d8cut, SFD is used instead of MFD."
111 " Applies only if no accumulation map is given.");
112 input.d8cut->required = NO;
113 input.d8cut->answer = "infinity";
114 input.d8cut->type = TYPE_DOUBLE;
115
116 input.mont_exp = G_define_option();
117 input.mont_exp->key = "mexp";
118 input.mont_exp->type = TYPE_DOUBLE;
119 input.mont_exp->required = NO;
120 input.mont_exp->answer = "0";
121 input.mont_exp->label =
122 _("Montgomery exponent for slope, disabled with 0");
123 input.mont_exp->description =
124 _("Montgomery: accumulation is multiplied with pow(slope,mexp) and then compared with threshold");
125
126 input.min_stream_length = G_define_option();
127 input.min_stream_length->key = "stream_length";
128 input.min_stream_length->type = TYPE_INTEGER;
129 input.min_stream_length->required = NO;
130 input.min_stream_length->answer = "0";
131 input.min_stream_length->label =
132 _("Delete stream segments shorter than stream_length cells");
133 input.min_stream_length->description =
134 _("Applies only to first-order stream segments (springs/stream heads)");
135
136 input.memory = G_define_option();
137 input.memory->key = "memory";
138 input.memory->type = TYPE_INTEGER;
139 input.memory->required = NO;
140 input.memory->answer = "300";
141 input.memory->label = _("Maximum memory to be used (in MB)");
142 input.memory->description = _("Cache size for raster rows");
143
144 output.stream_rast = G_define_standard_option(G_OPT_R_OUTPUT);
145 output.stream_rast->key = "stream_raster";
146 output.stream_rast->description =
147 _("Name for output raster map with unique stream ids");
148 output.stream_rast->required = NO;
149 output.stream_rast->guisection = _("Output maps");
150
151 output.stream_vect = G_define_standard_option(G_OPT_V_OUTPUT);
152 output.stream_vect->key = "stream_vector";
153 output.stream_vect->description =
154 _("Name for output vector map with unique stream ids");
155 output.stream_vect->required = NO;
156 output.stream_vect->guisection = _("Output maps");
157
158 output.dir_rast = G_define_standard_option(G_OPT_R_OUTPUT);
159 output.dir_rast->key = "direction";
160 output.dir_rast->description =
161 _("Name for output raster map with flow direction");
162 output.dir_rast->required = NO;
163 output.dir_rast->guisection = _("Output maps");
164
165 if (G_parser(argc, argv))
166 exit(EXIT_FAILURE);
167
168 /***********************/
169 /* check options */
170 /***********************/
171
172 /* input maps exist ? */
173 if (!G_find_raster(input.ele->answer, ""))
174 G_fatal_error(_("Raster map <%s> not found"), input.ele->answer);
175
176 if (input.acc->answer) {
177 if (!G_find_raster(input.acc->answer, ""))
178 G_fatal_error(_("Raster map <%s> not found"), input.acc->answer);
179 }
180
181 if (input.depression->answer) {
182 if (!G_find_raster(input.depression->answer, ""))
183 G_fatal_error(_("Raster map <%s> not found"), input.depression->answer);
184 have_depressions = 1;
185 }
186 else
187 have_depressions = 0;
188
189 /* threshold makes sense */
190 threshold = atof(input.threshold->answer);
191 if (threshold <= 0)
192 G_fatal_error(_("Threshold must be > 0 but is %f"), threshold);
193
194 /* d8cut */
195 if (strcmp(input.d8cut->answer, "infinity") == 0) {
196 d8cut = DBL_MAX;
197 }
198 else {
199 d8cut = atof(input.d8cut->answer);
200 if (d8cut < 0)
201 G_fatal_error(_("d8cut must be positive or zero but is %f"),
202 d8cut);
203 }
204
205 /* Montgomery stream initiation */
206 if (input.mont_exp->answer) {
207 mont_exp = atof(input.mont_exp->answer);
208 if (mont_exp < 0)
209 G_fatal_error(_("Montgomery exponent must be positive or zero but is %f"),
210 mont_exp);
211 if (mont_exp > 3)
212 G_warning(_("Montgomery exponent is %f, recommended range is 0.0 - 3.0"),
213 mont_exp);
214 }
215 else
216 mont_exp = 0;
217
218 /* Minimum stream segment length */
219 if (input.min_stream_length->answer) {
220 min_stream_length = atoi(input.min_stream_length->answer);
221 if (min_stream_length < 0)
222 G_fatal_error(_("Minimum stream length must be positive or zero but is %d"),
223 min_stream_length);
224 }
225 else
226 min_stream_length = 0;
227
228 if (input.memory->answer) {
229 memory = atoi(input.memory->answer);
230 if (memory <= 0)
231 G_fatal_error(_("Memory must be positive but is %d"),
232 memory);
233 }
234 else
235 memory = 300;
236
237 /* Check for some output map */
238 if ((output.stream_rast->answer == NULL)
239 && (output.stream_vect->answer == NULL)
240 && (output.dir_rast->answer == NULL)) {
241 G_fatal_error(_("At least one output raster maps must be specified"));
242 }
243
244 /*********************/
245 /* preparation */
246 /*********************/
247
248 /* open input maps */
249 mapset = G_find_raster2(input.ele->answer, "");
250 ele_fd = Rast_open_old(input.ele->answer, mapset);
251
252 if (input.acc->answer) {
253 mapset = G_find_raster2(input.acc->answer, "");
254 acc_fd = Rast_open_old(input.acc->answer, mapset);
255 }
256 else
257 acc_fd = -1;
258
259 if (input.depression->answer) {
260 mapset = G_find_raster2(input.depression->answer, "");
261 depr_fd = Rast_open_old(input.depression->answer, mapset);
262 }
263 else
264 depr_fd = -1;
265
266 /* set global variables */
267 nrows = Rast_window_rows();
268 ncols = Rast_window_cols();
269 sides = 8; /* not a user option */
270 c_fac = 5; /* not a user option, MFD covergence factor 5 gives best results */
271
272 /* segment structures */
273 seg_rows = seg_cols = 64;
274 seg2kb = seg_rows * seg_cols / 1024.;
275
276 /* balance segment files */
277 /* elevation + accumulation: * 2 */
278 memory_divisor = sizeof(WAT_ALT) * 2;
279 disk_space = sizeof(WAT_ALT);
280 /* stream ids: / 2 */
281 memory_divisor += sizeof(int) / 2.;
282 disk_space += sizeof(int);
283
284 /* aspect and flags: * 2 */
285 memory_divisor += sizeof(ASP_FLAG) * 4;
286 disk_space += sizeof(ASP_FLAG);
287
288 /* astar_points: / 16 */
289 /* ideally only a few but large segments */
290 memory_divisor += sizeof(POINT) / 16.;
291 disk_space += sizeof(POINT);
292 /* heap points: / 4 */
293 memory_divisor += sizeof(HEAP_PNT) / 4.;
294 disk_space += sizeof(HEAP_PNT);
295
296 /* KB -> MB */
297 memory_divisor *= seg2kb / 1024.;
298 disk_space *= seg2kb / 1024.;
299
300 num_open_segs = memory / memory_divisor;
301 heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
302 (4. * 1024.);
303 num_seg_total = (ncols / seg_cols + 1) * (nrows / seg_rows + 1);
304 if (num_open_segs > num_seg_total) {
305 heap_mem += (num_open_segs - num_seg_total) * memory_divisor;
306 heap_mem -= (num_open_segs - num_seg_total) * seg2kb *
307 sizeof(HEAP_PNT) / (4. * 1024.);
308 num_open_segs = num_seg_total;
309 }
310 if (num_open_segs < 16) {
311 num_open_segs = 16;
312 heap_mem = num_open_segs * seg2kb * sizeof(HEAP_PNT) /
313 (4. * 1024.);
314 }
315 G_verbose_message(_("%.2f%% of data are kept in memory"),
316 100. * num_open_segs / num_seg_total);
317 disk_space *= num_seg_total;
318 if (disk_space < 1024.0)
319 G_verbose_message(_("Will need up to %.2f MB of disk space"), disk_space);
320 else
321 G_verbose_message(_("Will need up to %.2f GB (%.0f MB) of disk space"),
322 disk_space / 1024.0, disk_space);
323
324 /* open segment files */
325 G_verbose_message(_("Creating temporary files..."));
326 seg_open(&watalt, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
327 sizeof(WAT_ALT), 1);
328 if (num_open_segs * 2 > num_seg_total)
329 heap_mem += (num_open_segs * 2 - num_seg_total) * seg2kb *
330 sizeof(WAT_ALT) / 1024.;
331 cseg_open(&stream, seg_rows, seg_cols, num_open_segs / 2.);
332
333 seg_open(&aspflag, nrows, ncols, seg_rows, seg_cols, num_open_segs * 2,
334 sizeof(ASP_FLAG), 1);
335/*
336 bseg_open(&asp, seg_rows, seg_cols, num_open_segs);
337 bseg_open(&bitflags, seg_rows, seg_cols, num_open_segs * 4);
338*/
339
340 if (num_open_segs * 4 > num_seg_total)
341 heap_mem += (num_open_segs * 4 - num_seg_total) * seg2kb / 1024.;
342
343 /* load maps */
344 if (load_maps(ele_fd, acc_fd) < 0)
345 G_fatal_error(_("Unable to load input raster map(s)"));
346 else if (!n_points)
347 G_fatal_error(_("No non-NULL cells in input raster map(s)"));
348
349 G_debug(1, "open segments for A* points");
350 /* columns per segment */
351 seg_cols = seg_rows * seg_rows;
352 num_seg_total = n_points / seg_cols;
353 if (n_points % seg_cols > 0)
354 num_seg_total++;
355 /* no need to have more segments open than exist */
356 num_open_array_segs = num_open_segs / 16.;
357 if (num_open_array_segs > num_seg_total)
358 num_open_array_segs = num_seg_total;
359 if (num_open_array_segs < 1)
360 num_open_array_segs = 1;
361
362 G_debug(1, "segment size for A* points: %d", seg_cols);
363 seg_open(&astar_pts, 1, n_points, 1, seg_cols, num_open_array_segs,
364 sizeof(POINT), 1);
365
366 /* one-based d-ary search_heap with astar_pts */
367 G_debug(1, "open segments for A* search heap");
368
369 /* allowed memory for search heap in MB */
370 G_debug(1, "heap memory %.2f MB", heap_mem);
371 /* columns per segment */
372 /* larger is faster */
373 seg_cols = seg_rows * seg_rows * seg_rows;
374 num_seg_total = n_points / seg_cols;
375 if (n_points % seg_cols > 0)
376 num_seg_total++;
377 /* no need to have more segments open than exist */
378 num_open_array_segs = (1 << 20) * heap_mem / (seg_cols * sizeof(HEAP_PNT));
379 if (num_open_array_segs > num_seg_total)
380 num_open_array_segs = num_seg_total;
381 if (num_open_array_segs < 2)
382 num_open_array_segs = 2;
383
384 G_debug(1, "A* search heap open segments %d, total %d",
385 num_open_array_segs, num_seg_total);
386 G_debug(1, "segment size for heap points: %d", seg_cols);
387 /* the search heap will not hold more than 5% of all points at any given time ? */
388 /* chances are good that the heap will fit into one large segment */
389 seg_open(&search_heap, 1, n_points + 1, 1, seg_cols,
390 num_open_array_segs, sizeof(HEAP_PNT), 1);
391
392 /********************/
393 /* processing */
394 /********************/
395
396 /* initialize A* search */
397 if (init_search(depr_fd) < 0)
398 G_fatal_error(_("Unable to initialize search"));
399
400 /* sort elevation and get initial stream direction */
401 if (do_astar() < 0)
402 G_fatal_error(_("Unable to sort elevation raster map values"));
403 seg_close(&search_heap);
404
405 if (acc_fd < 0) {
406 /* accumulate surface flow */
407 if (do_accum(d8cut) < 0)
408 G_fatal_error(_("Unable to calculate flow accumulation"));
409 }
410
411 /* extract streams */
412 if (extract_streams(threshold, mont_exp, acc_fd < 0) < 0)
413 G_fatal_error(_("Unable to extract streams"));
414
415 seg_close(&astar_pts);
416 seg_close(&watalt);
417
418 /* thin streams */
419 if (thin_streams() < 0)
420 G_fatal_error(_("Unable to thin streams"));
421
422 /* delete short streams */
423 if (min_stream_length) {
424 if (del_streams(min_stream_length) < 0)
425 G_fatal_error(_("Unable to delete short stream segments"));
426 }
427
428 /* write output maps */
429 if (close_maps(output.stream_rast->answer, output.stream_vect->answer,
430 output.dir_rast->answer) < 0)
431 G_fatal_error(_("Unable to write output raster maps"));
432
433 cseg_close(&stream);
434 seg_close(&aspflag);
435
436 exit(EXIT_SUCCESS);
437}
Note: See TracBrowser for help on using the repository browser.