Geoinformatica  0.90
TerrainAnalysis.pm
Go to the documentation of this file.
1 #** @file TerrainAnalysis.pm
2 # @brief Adds terrain analysis methods into Geo::Raster
3 # In this file there are methods mainly for digital elevation model
4 # (DEM) rasters, for flow direction (FDG) rasters, for streams
5 # rasters.
6 #*
7 
8 package Geo::Raster;
9 
10 use strict;
11 use Scalar::Util 'blessed';
12 use File::Basename; # for fileparse
14 
15 #** @method @fit_surface($z_factor)
16 #
17 # @brief Fit a 9-term quadratic polynomial to the 3*3 neighborhood of
18 # each cell in a DEM.
19 #
20 # The 9-term quadratic polynomial:
21 #
22 # z = Ax^2y^2 + Bx^2y + Cxy^2 + Dx^2 + Ey^2 + Fxy + Gx + Hy + I
23 #
24 # @see Moore et al. 1991. Hydrol. Proc. 5, 3-30.
25 # @param[in] z_factor is the unit of z divided by the unit of x and y, the
26 # default value of z_factor is 1.
27 # @return 9 rasters in a list, one for each parameter.
28 #*
29 sub fit_surface {
30  my($dem, $z_factor) = @_;
31  $z_factor = 1 unless $z_factor;
32  my $a = ral_dem_fit_surface($dem->{GRID}, $z_factor);
33  my @ret;
34  if ($a and @$a) {
35  for my $a (@$a) {
36  my $b = Geo::Raster->new($a);
37  push @ret, $b;
38  }
39  }
40  return @ret;
41 }
42 
43 #** @method Geo::Raster aspect()
44 #
45 # @brief Estimate aspects from a DEM.
46 #
47 # The aspect is computed from a 9-term quadratic polynomial fitted
48 # independently on each cell of the DEM. The aspect is stored in
49 # radians increasing clockwise starting from zero in north. For flat
50 # cells the aspect is undefined, denoted with the value -1.
51 # @return an aspect raster. In void context converts the DEM.
52 #*
53 sub aspect {
54  my $self = shift;
55  if (defined wantarray) {
56  return Geo::Raster->new(ral_dem_aspect($self->{GRID}));
57  } else {
58  $self->_new_grid(ral_dem_aspect($self->{GRID}));
59  }
60 }
61 
62 #** @method Geo::Raster slope(scalar z_factor)
63 #
64 # @brief Estimate the slope from a DEM.
65 #
66 # The slope is computed from a 9-term quadratic polynomial fitted
67 # independently on each cell of the DEM. The slope is stored in
68 # radians.
69 # @param[in] z_factor The unit of z divided by the unit of x and
70 # y. Default is 1.
71 # @return a slope raster. In void context converts the DEM.
72 #*
73 sub slope {
74  my $self = shift;
75  my $z_factor = shift;
76  $z_factor = 1 unless $z_factor;
77  if (defined wantarray) {
78  return Geo::Raster->new(ral_dem_slope($self->{GRID}, $z_factor));
79  } else {
80  $self->_new_grid(ral_dem_slope($self->{GRID}, $z_factor));
81  }
82 }
83 
84 #** @method Geo::Raster fdg(%params)
85 #
86 # @brief Compute a flow direction raster (FDG) from a DEM.
87 #
88 # @param[in] params A hash of named parameters:
89 # - <I>method</I>=>string (optional). The method for computing the
90 # FDG. Currently supported methods are 'D8' and 'Rho8' and 'many'. D8
91 # selects the neighbor with steepest descent. Rho8 selects randomly,
92 # using the descent as a weight, a lower cell. In the 'many' method
93 # all lower cells are coded with the 8 bits of a byte. The directions
94 # are from 1 (up) to 8 (up left)<br>:
95 # 8 1 2<br>
96 # 7 X 3<br>
97 # 6 5 4<br
98 # Flat cells (no lower neighbors) are marked -1. Pit cells (all
99 # neighbor cells are higher) are marked with 0. The default is D8.
100 # - <I>drain_all</I>=>boolean (optional). Whether to use iteration to
101 # produce a FDG with drainage resolved for all cell. The iteration
102 # algorithm first applies the drain_flat_areas method using the option
103 # 'multiple pour points', then the same method with option 'one pour
104 # point'. Next the drainage of all depressions is iteratively resolved
105 # using the drain_depressions method. The default is false.
106 # - <I>quiet</I>=>boolean (optional). Whether to report the result
107 # (number of cells or areas drained).
108 # @exception - Unsupported method
109 # @exception - No progress in the iteration
110 # @return a FDG. In void context converts the DEM.
111 #*
112 sub fdg {
113  my($dem, %opt) = @_;
114  $opt{method} = 'D8' unless $opt{method};
115  my $method;
116  if ($opt{method} eq 'D8') {
117  $method = 1;
118  } elsif ($opt{method} eq 'Rho8') {
119  $method = 2;
120  } elsif ($opt{method} eq 'many') {
121  $method = 3;
122  } else {
123  croak "fdg: $opt{method}: unsupported method";
124  }
125  my $fdg = Geo::Raster->new(ral_dem_fdg($dem->{GRID}, $method));
126 
127  if ($opt{drain_all}) {
128  my $step = 1;
129  my $pits_last_time = -1;
130  my $flats_last_time = -1;
131  while (1) {
132  $fdg->drain_flat_areas($dem, method=>'m', quiet=>1);
133  $fdg->drain_flat_areas($dem, method=>'o', quiet=>1);
134  my $c = $fdg->contents();
135  my $pits = $$c{0} || 0;
136  my $flats = $$c{-1} || 0;
137  print "drain_all: iteration step $step: $pits pits and $flats flat cells\n" unless $opt{quiet};
138  last if ($pits == 0 and $flats == 0);
139  my $n = ral_fdg_drain_depressions($fdg->{GRID}, $dem->{GRID});
140  print "drain_all: iteration step $step: $n depressions fixed\n" unless $opt{quiet};
141  croak "there is no progress" if ($pits_last_time == $pits and $flats_last_time == $flats);
142  $pits_last_time = $pits;
143  $flats_last_time = $flats;
144  $step++;
145  }
146  }
147 
148  if (defined wantarray) {
149  return $fdg;
150  } else {
151  $dem->_new_grid($fdg->{GRID});
152  delete $fdg->{GRID};
153  }
154 }
155 
156 sub many2ds {
157  my($fdg) = @_;
158  my %map;
159  for my $i (1..255) {
160  my $c = 0;
161  for my $j (0..7) {
162  $c++ if $i & 1 << $j;
163  }
164  $map{$i} = $c;
165  }
166  $fdg->map(\%map);
167 }
168 
169 #** @method @movecell(@cell, $dir)
170 #
171 # @brief Return the cell in the given direction.
172 # @param[in] cell The current cell
173 # @param[in] dir (optional) Direction of the movement. If not given,
174 # the direction is taken from the cell, thus assuming this raster is a
175 # FDG.
176 # @return the cell in the given direction or undef if the cell is not
177 # within the raster.
178 #*
179 sub movecell {
180  my($fdg, $i, $j, $dir) = @_;
181  $dir = $fdg->get($i, $j) unless $dir;
182  SWITCH: {
183  if ($dir == 1) { $i--; last SWITCH; }
184  if ($dir == 2) { $i--; $j++; last SWITCH; }
185  if ($dir == 3) { $j++; last SWITCH; }
186  if ($dir == 4) { $i++; $j++; last SWITCH; }
187  if ($dir == 5) { $i++; last SWITCH; }
188  if ($dir == 6) { $i++; $j--; last SWITCH; }
189  if ($dir == 7) { $j--; last SWITCH; }
190  if ($dir == 8) { $i--; $j--; last SWITCH; }
191  croak "movecell: $dir: bad direction";
192  }
193  if ($fdg) {
194  return if ($i < 0 or $j < 0 or $i >= ral_grid_get_height($fdg->{GRID}) or $j >= ral_grid_get_width($fdg->{GRID}));
195  }
196  return ($i, $j);
197 }
198 
199 *downstream = *movecell;
200 
201 ## @fn $dirsum($dir, $add)
202 #
203 # @brief Increases the given direction with <i>add</i> steps clockwise.
204 # @param[in] dir Direction
205 # @param[in] add Steps to increase dir.
206 # @return the new direction.
207 #*
208 sub dirsum {
209  my($dir, $add) = @_;
210  $dir += $add;
211  $dir -= 8 if $dir > 8;
212  return $dir;
213 }
214 
215 #** @method Geo::Raster drain_flat_areas(Geo::Raster dem, hash params)
216 #
217 # @brief Resolve the flow direction for flat areas in a FDG.
218 #
219 # The method uses either "one pour point" (short "o") or "multiple
220 # pour points" (short "m") technique for resolving the drainage of
221 # flat area cells of a FDG.
222 #
223 # In the one pour point technique (the default) all cells on the flat
224 # area are drained to one pour point cell. The pour point cell is
225 # selected either from the border of the flat area or just outside the
226 # flat area depending whether the outside cell is not higher than the
227 # inside cell. If the pour point is on the border it is converted into
228 # a pit cell. Thus this technique is guaranteed to produce a flatless
229 # FDG but it may increase the number of pits.
230 #
231 # In the multiple pour points approach all flat cells are iteratively
232 # drained to any non-higher neighbor whose drainage is resolved. Thus
233 # this technique is not guaranteed to produce a flatless FDG.
234 #
235 # @param[in] dem The DEM.
236 # @param[in] params Named parameters:
237 # - <I>method</I>=>string (optional) Either "one pour point" (short
238 # "o") or "multiple pour points" (short "m"). Default is one pour
239 # point.
240 # - <I>quiet</I>=>boolean (optional) Whether to report the result
241 # (number of cells or areas drained).
242 # @exception - No DEM supplied
243 # @exception - Unsupported method
244 # @return In void context the method changes this flow direction raster,
245 # otherwise the method returns a new FDG.
246 #*
247 sub drain_flat_areas {
248  my($fdg, $dem, %opt) = @_;
249  croak "drain_flat_areas: no DEM supplied" unless $dem and ref($dem);
250  $fdg = Geo::Raster->new($fdg) if defined wantarray;
251  $opt{method} = 'one pour point' unless $opt{method};
252  if ($opt{method} =~ /^m/) {
253  my $n = ral_fdg_drain_flat_areas1($fdg->{GRID}, $dem->{GRID});
254  print "drain_flat_areas (multiple pour points): $n flat cells drained\n" unless $opt{quiet};
255  } elsif ($opt{method} =~ /^o/) {
256  my $n = ral_fdg_drain_flat_areas2($fdg->{GRID}, $dem->{GRID});
257  print "drain_flat_areas (one pour point): $n flat areas drained\n" unless $opt{quiet};
258  } else {
259  croak "drain_flat_areas: $opt{method}: unknown method";
260  }
261  return $fdg if defined wantarray;
262 }
263 
264 #** @method Geo::Raster drain_depressions(Geo::Raster dem)
265 #
266 # @brief Scan FDG once and drain the depressions that are found.
267 #
268 # This method scans the given FDG once and drains all depressions that
269 # are found in the FDG by reversing the flowpath from the lowest pour
270 # point of the depression to the pit cell. The DEM remains unchanged.
271 # @param[in] dem The DEM. The DEM is not changed in the method.
272 # @return In a void context the method changes this FDG, otherwise the
273 # method returns a new FDG.
274 #*
275 sub drain_depressions {
276  my($fdg, $dem) = @_;
277  $fdg = Geo::Raster->new($fdg) if defined wantarray;
278  ral_fdg_drain_depressions($fdg->{GRID}, $dem->{GRID});
279  return $fdg if defined wantarray;
280 }
281 
282 #** @method @outlet(@cell)
283 #
284 # @brief Return the outlet of a catchment in a FDG.
285 #
286 # @param[in] cell A cell on the catchment.
287 # @return the outlet cell of the catchment.
288 #*
289 sub outlet {
290  my($fdg, $streams, @cell) = @_;
291  my $cell = ral_fdg_outlet($fdg->{GRID}, $streams->{GRID}, @cell);
292  return @{$cell};
293 }
294 
295 #** @method Geo::Raster ucg()
296 #
297 # @brief Compute an upslope cell raster (UCG) from a FDG.
298 #
299 # In a UCG the upslope cells of a cell are coded with the 8 bits of a
300 # byte. The directions are from 1 (up) to 8 (up left).
301 #
302 # @return a UCG. In void context converts the FDG.
303 #*
304 sub ucg {
305  my($dem) = @_;
306  my $ucg = ral_dem_ucg($dem->{GRID});
307  if (defined wantarray) {
308  return Geo::Raster->new($ucg);
309  } else {
310  $dem->_new_grid($ucg);
311  }
312 }
313 
314 #** @method @upstream(Geo::Raster streams, array cell)
315 #
316 # @brief Return the direction(s) to the upslope cells of a cell in a FDG.
317 #
318 # Example of getting directions to the upstream stream cells:
319 # @code
320 # @up = $fdg->upstream($streams, $row, $column);
321 # @endcode
322 #
323 # @param[in] streams (optional) A streams raster.
324 # @param[in] cell The cell.
325 # @return The directions to the upstream cells. If streams raster is
326 # given, only directions to stream cells are returned. The directions
327 # are coded as usual: 1 is up, 2 is up right, etc.
328 #*
329 sub upstream {
330  my $fdg = shift;
331  my $streams;
332  my @cell;
333  if (@_ > 2) {
334  ($streams,@cell) = @_;
335  } else {
336  @cell = @_;
337  }
338  my @up;
339  my $d;
340  for $d (1..8) {
341  my @test = $fdg->movecell(@cell, $d);
342  next unless @test;
343  my $u = $fdg->get(@test);
344  next if $streams and !($streams->get(@test));
345  if ($u == ($d - 4 <= 0 ? $d + 4 : $d - 4)) {
346  push @up, $d;
347  }
348  }
349  return @up;
350 }
351 
352 #** @method Geo::Raster raise_pits(%params)
353 #
354 # @brief Raise each pit cell to the level of its lowest neighbor in a
355 # DEM.
356 #
357 # @param[in] params Named parameters:
358 # - <I>z_limit</I>=>number (optional) A threshold value for how much
359 # lower that its neighbors a cell may be before it is raised. Default
360 # is 0.
361 # - <I>quiet</I>=>boolean (optional) Whether the method should report
362 # the number of cells raised.
363 # @return In void context the method changes this DEM, otherwise the
364 # method returns a new DEM.
365 #*
366 sub raise_pits {
367  my($dem, %opt) = @_;
368  $opt{z_limit} = 0 unless defined($opt{z_limit});
369  $dem = Geo::Raster->new($dem) if defined wantarray;
370  my $n = ral_dem_raise_pits($dem->{GRID}, $opt{z_limit});
371  print "raise_pits: $n pit cells raised\n" unless $opt{quiet};
372  return $dem if defined wantarray;
373 }
374 
375 #** @method lower_peaks(%params)
376 #
377 # @brief Lower each peak cell to the level of its highest neighbor in
378 # a DEM.
379 #
380 # @param[in] params Named parameters:
381 # - <I>z_limit</I>=>number (optional) A threshold value for how much
382 # higher than its neighbors a cell may be before it is
383 # lowered. Default is 0.
384 # - <I>quiet</I>=>boolean (optional) Whether the method should report
385 # the number of cells raised.
386 # @return In void context the method changes this DEM, otherwise
387 # the method returns a new DEM.
388 #*
389 sub lower_peaks {
390  my($dem, %opt) = @_;
391  $opt{z_limit} = 0 unless defined($opt{z_limit});
392  $dem = Geo::Raster->new($dem) if defined wantarray;
393  my $n = ral_dem_lower_peaks($dem->{GRID}, $opt{z_limit});
394  print "lower_peaks: $n peak cells lowered\n" unless $opt{quiet};
395  return $dem if defined wantarray;
396 }
397 
398 #** @method Geo::Raster depressions($inc_m)
399 #
400 # @brief Return depressions defined by a FDG.
401 #
402 # @param[in] inc_m (optional) A boolean value indicating whether each
403 # depression is marked with a unique integer. Default is false.
404 # @return Depressions raster.
405 #*
406 sub depressions {
407  my($fdg, $inc_m) = @_;
408  $inc_m = 0 unless defined($inc_m) and $inc_m;
409  return Geo::Raster->new(ral_fdg_depressions($fdg->{GRID}, $inc_m));
410 }
411 
412 #** @method scalar fill_depressions(%params)
413 #
414 # @brief Fill the depressions in a DEM.
415 #
416 # The depressions in the DEM are filled up to the lowest cell just
417 # outside the depression. The depressions are obtained from a FDG.
418 # @param[in] params Named parameters:
419 # - <i>iterative</i>=>boolean (optional) Whether to run the depression
420 # filling algorithm iteratively until a pitless and flatless FDG is
421 # obtained from the DEM using the D8 method. In an iteration step a
422 # FDG is first computed from the DEM using the D8 method. Then the
423 # flat areas are removed from the FDG using first the multiple pour
424 # point method and then the one pour point method. The depression
425 # removing algorithm is then run once (one scan of the whole raster)
426 # on the DEM using this FDG. If the flatless FDG that is computed from
427 # this changed DEM is also pitless, then the iteration stops. If
428 # iterative is true, the DEM is changed and a pitless and flatless FDG
429 # is returned. Default is true.
430 # - <i>FDG</i>=>raster (optional unless iterative is false) The FDG
431 # for the depression filling algorithm. If FDG is given, the
432 # depression filling algorithm is run only once, i.e., one scan of the
433 # FDG is performed. Default is undefined.
434 # - <i>quiet</i>=>boolean (optional) Whether to report the progressing
435 # of the iteration.
436 # @exception - FDG is not given but iterative is false.
437 # @exception - There is no progress in the iteration.
438 # @return a DEM from which some depressions are removed (if the
439 # context is non-void and iterative is false), the number of filled
440 # depressions (if FDG is given), or a pitless and flatless FDG.
441 #*
442 sub fill_depressions {
443  my($dem, %opt) = @_;
444  $opt{iterative} = 1 unless exists $opt{iterative} and CORE::not $opt{iterative};
445  $opt{fdg} = $opt{FDG} if exists $opt{FDG};
446  if (CORE::not $opt{iterative}) {
447  croak "fill_depressions: FDG needed if not iterative" unless $opt{fdg};
448  $dem = Geo::Raster->new($dem) if defined wantarray;
449  ral_dem_fill_depressions($dem->{GRID}, $opt{fdg}->{GRID});
450  return $dem if defined wantarray;
451  return;
452  }
453  if ($opt{fdg}) {
454  $dem = Geo::Raster->new($dem) if defined wantarray;
455  ral_dem_fill_depressions($dem->{GRID}, $opt{fdg}->{GRID});
456  return $dem if defined wantarray;
457  return;
458  } else {
459  my $step = 1;
460  my $pits_last_time = -1;
461  my $flats_last_time = -1;
462  while (1) {
463  my $fdg = $dem->fdg(method=>'D8', quiet=>1);
464  $fdg->drain_flat_areas($dem, method=>'m', quiet=>1);
465  $fdg->drain_flat_areas($dem, method=>'o', quiet=>1);
466  my $c = $fdg->contents();
467  my $pits = $$c{0} || 0;
468  my $flats = $$c{-1} || 0;
469  print "fill_depressions: iteration step $step: $pits pits and $flats flat cells\n" unless $opt{quiet};
470  return $fdg if ($pits == 0 and $flats == 0);
471  croak "there is no progress" if ($pits_last_time == $pits and $flats_last_time == $flats);
472  my $n = ral_dem_fill_depressions($dem->{GRID}, $fdg->{GRID});
473  print "fill_depressions: iteration step $step: $n depressions filled\n" unless $opt{quiet};
474  $pits_last_time = $pits;
475  $flats_last_time = $flats;
476  $step++;
477  Gtk2->main_iteration while Gtk2->events_pending;
478  }
479  }
480 }
481 
482 #** @method scalar breach(%params)
483 #
484 # @brief Breach the depressions in a DEM.
485 #
486 # Breaching is a depression removal method, which lowers the elevation
487 # of cells, which form a dam. Breaching is tried at the lowest cell on
488 # the rim of the depression which has the steepest descent away from
489 # the depression (if there are more than one lowest cells) and the
490 # steepest descent into the depression (if there are more than one
491 # lowest cells with identical slope out)
492 #
493 # The breaching algorithm implemented here is close to but not the
494 # same as in Martz and Garbrecht (1998). The biggest difference is
495 # that the depression cells are not raised in this implementation.
496 #
497 # @param[in] params Named parameters:
498 # - <I>limit</I> (optional) Maximum amount of cells (width of the dam)
499 # to be breached. Default is to not limit the breaching ($limit == 0).
500 # - <i>iterative</i>=>boolean (optional) Whether to run the breach
501 # algorithm iteratively until a pitless and flatless FDG is obtained
502 # from the DEM using the D8 method or there is no progress in the
503 # iteration. In an iteration step a FDG is first computed from the DEM
504 # using the D8 method. Then the flat areas are removed from the FDG
505 # using first the multiple pour point method and then the one pour
506 # point method. The breaching algorithm is then run once (one scan of
507 # the whole raster) on the DEM using this FDG. If the flatless FDG
508 # that is computed from this changed DEM is also pitless or there is
509 # no progress in the iteration, then the iteration stops. If iterative
510 # is true, the DEM is changed and a FDG is returned. Default is true.
511 # - <i>FDG</i>=>raster (optional unless iterative is false) The FDG
512 # for the breaching algorithm. If FDG is given it must not contain
513 # flat areas. The algorithm is run only once, i.e., one scan of the
514 # FDG is performed. Default is undefined.
515 # - <i>quiet</i>=>boolean (optional) Whether to not to report the
516 # progressing of the iteration.
517 # @return a DEM from which some depressions are removed (if the
518 # context is non-void and iterative is false), nothing (if FDG is
519 # given), or a pitless and flatless FDG.
520 # @exception - FDG is not given but iterative is false.
521 # @see Martz, L.W. and Garbrecht, J. 1998. The treatment of flat areas
522 # and depressions in automated drainage analysis of raster digital
523 # elevation models. Hydrol. Process. 12, 843-855
524 #*
525 sub breach {
526  my($dem, %opt) = @_;
527  $opt{fdg} = $opt{FDG} if exists $opt{FDG};
528  $opt{limit} = 0 unless defined($opt{limit});
529  $opt{iterative} = 1 unless exists $opt{iterative} and CORE::not $opt{iterative};
530  if (CORE::not $opt{iterative}) {
531  croak "breach: FDG needed if not iterative" unless $opt{fdg};
532  $dem = Geo::Raster->new($dem) if defined wantarray;
533  ral_dem_breach($dem->{GRID}, $opt{fdg}->{GRID}, $opt{limit});
534  return $dem if defined wantarray;
535  return;
536  }
537  if ($opt{fdg}) {
538  $dem = Geo::Raster->new($dem) if defined wantarray;
539  return ral_dem_breach($dem->{GRID}, $opt{fdg}->{GRID}, $opt{limit});
540  } else {
541  my $step = 1;
542  my $pits_last_time = -1;
543  my $flats_last_time = -1;
544  while (1) {
545  my $fdg = $dem->fdg(method=>'D8', quiet=>1);
546  $fdg->drain_flat_areas($dem, method=>'m', quiet=>1);
547  $fdg->drain_flat_areas($dem, method=>'o', quiet=>1);
548  my $c = $fdg->contents();
549  my $pits = $$c{0} || 0;
550  my $flats = $$c{-1} || 0;
551  print "breach: iteration step $step: $pits pits and $flats flat cells\n" unless $opt{quiet};
552  return $fdg if ($pits == 0 and $flats == 0);
553  return $fdg if ($pits_last_time == $pits and $flats_last_time == $flats);
554  my $n = ral_dem_breach($dem->{GRID}, $fdg->{GRID}, $opt{limit});
555  print "breach: iteration step $step: $n depressions filled\n" unless $opt{quiet};
556  $pits_last_time = $pits;
557  $flats_last_time = $flats;
558  $step++;
559  }
560  }
561 }
562 
563 #** @method Geo::Raster path(@cell, Geo::Raster stop)
564 #
565 # @brief Return the flow path from the given FDG cell onwards.
566 #
567 # The end of the path is where flow direction is not specified, where it goes
568 # out of the raster, or where the stop raster has a positive value.
569 #
570 # @param[in] cell The origin of the path.
571 # @param[in] stop (optional) Raster denoting end cells for paths.
572 # @return raster, where the path cells have the value 1 and otherwise
573 # no data value. In void context changes the FDG.
574 #*
575 sub path {
576  my($fdg, $i, $j, $stop) = @_;
577  my $g = ral_fdg_path($fdg->{GRID}, $i, $j, $stop ? $stop->{GRID} : undef);
578  if (defined wantarray) {
579  return Geo::Raster->new($g);
580  } else {
581  $fdg->_new_grid($g);
582  }
583 }
584 
585 #** @method Geo::Raster path_length(Geo::Raster stop, Geo::Raster op)
586 #
587 # @brief Compute a path length raster from a FDG.
588 #
589 # The path is assumed to go from a center point of a cell to another
590 # center point. The length is not recorded if op is nodata. The length
591 # is calculated in the raster units. A path ends at the border of the
592 # FDG, at a FDG cell with undefined direction, or at a cell in the
593 # stop raster with positive value.
594 # @param[in] stop (optional) Raster denoting end cells for paths.
595 # @param[in] op (optional) Raster denoting cells which are included in
596 # the length computation.
597 # @return a flow path length raster. In void context changes the FDG.
598 #*
599 sub path_length {
600  my($fdg, $stop, $op) = @_;
601  my $g = ral_fdg_path_length($fdg->{GRID}, $stop ? $stop->{GRID} : undef, $op ? $op->{GRID} : undef);
602  if (defined wantarray) {
603  return new Geo::Raster $g;
604  } else {
605  $fdg->_new_grid($g);
606  }
607 }
608 
609 #** @method Geo::Raster path_sum(Geo::Raster stop, Geo::Raster op)
610 #
611 # @brief Compute a cost-to-go raster from a FDG.
612 #
613 # This FDG method returns a raster, where the value of each cell is
614 # the weighted (with length) sum along the path from that cell to the
615 # end of the path. The path is assumed to go from a center point of a
616 # cell to another center point. The length is not recorded if op is
617 # nodata. The length is calculated in the raster units.
618 # @param[in] stop (optional) Raster denoting end cells for paths.
619 # @param[in] op Weights (cost) for the summing.
620 # @return a cost-to-go raster. In void context changes the FDG.
621 #*
622 sub path_sum {
623  my($fdg, $stop, $op) = @_;
624  my $g = ral_fdg_path_sum($fdg->{GRID}, $stop ? $stop->{GRID} : undef, $op->{GRID});
625  if (defined wantarray) {
626  return new Geo::Raster $g;
627  } else {
628  $fdg->_new_grid($g);
629  }
630 }
631 
632 #** @method Geo::Raster upslope_count(Geo::Raster mask, $include_self)
633 #
634 # @brief Compute the count of the upslope cells in a FDG.
635 #
636 # @param[in] mask (optional) Can be used to mask out cells from the
637 # count. Nodata cells of the mask are not included in the count.
638 # @param[in] include_self (optional) Boolean value, which specifies
639 # whether the cell itself is included in its upslope area. Default is
640 # true.
641 # @return the upslope count raster. In void context changes the FDG.
642 # @note DO NOT call if the FDG contains loops.
643 # @note This is the method for computing an upslope area raster (UAG)
644 # from a flow direction raster (FDG).
645 #*
646 sub upslope_count {
647  my($fdg, $a, $b) = @_;
648  my $op = ref $a ? $a : $b;
649  my $include_self = ref $a ? $b : $a;
650  $include_self = 1 unless defined $include_self;
651  my $g = $op ?
652  ral_fdg_upslope_count($fdg->{GRID}, $op->{GRID}, $include_self) :
653  ral_fdg_upslope_count_without_op($fdg->{GRID}, $include_self);
654  if (defined wantarray) {
655  return new Geo::Raster $g;
656  } else {
657  $fdg->_new_grid($g);
658  }
659 }
660 
661 #** @method Geo::Raster upslope_sum(Geo::Raster a, $include_self)
662 #
663 # @brief Compute the sum of the values of the upslope cells in a
664 # raster (FDG method).
665 #
666 # @param[in] a The operand raster.
667 # @param[in] include_self (optional). Boolean value specifying whether
668 # the cell itself is included in its upslope area. Default is true.
669 # @return In void context changes this flow direction raster,
670 # otherwise returns a new FDG.
671 # @note DO NOT call if the FDG contains loops.
672 #*
673 sub upslope_sum {
674  my($fdg, $a, $b) = @_;
675  croak "usage: \$fdg->upslope_sum(\$op); where \$op is a raster whose values are to be summed" unless $a and $a->{GRID};
676  $b = 1 unless defined $b;
677  my $g = ral_fdg_upslope_sum($fdg->{GRID}, $a->{GRID}, $b);
678  if (defined wantarray) {
679  return Geo::Raster->new($g);
680  } else {
681  $fdg->_new_grid($g);
682  }
683 }
684 
685 #** @method Geo::Raster kill_extra_outlets(Geo::Raster lakes, Geo::Raster uag)
686 #
687 # @brief Checks and possibly correct the sanity of the flow paths in a
688 # terrain with lakes (FDG method).
689 #
690 # The FDG is modified so that each lake has only one outlet. The lakes
691 # raster is typically a reclassified land cover raster.
692 # @param[in] lakes an integer raster of the lakes in the terrain (each
693 # lake may have its own non-zero id)
694 # @param[in] uag (optional) Upslope area raster. Computed using the
695 # upslope_count method unless given.
696 # @return In void context changes this flow direction raster,
697 # otherwise returns a new FDG.
698 #*
699 sub kill_extra_outlets {
700  my ($fdg, $lakes, $uag) = @_;
701  $fdg = new Geo::Raster $fdg if defined wantarray;
702  $uag = $fdg->upslope_count unless $uag;
703  ral_fdg_kill_extra_outlets($fdg->{GRID}, $lakes->{GRID}, $uag->{GRID});
704  return $fdg if defined wantarray;
705 }
706 
707 #** @method Geo::Raster catchment(Geo::Raster catchment, @cell, $m)
708 #
709 # @brief Return the catchment area of the given cell (FDG method).
710 #
711 # @param[in,out] catchment (optional) The raster in which to return
712 # the catchment of the given cell.
713 # @param[in] cell The cell.
714 # @param[in] m (optional). Number used to mark the cells belonging to
715 # the catchment. Default is 1.
716 # @return Depending on the context returns the catchment raster or a
717 # list containing the catchment raster and the number of cells in the
718 # catchment.
719 #*
720 sub catchment {
721  my $fdg = shift;
722  my $i = shift;
723  my ($M, $N) = $fdg->size();
724  my ($j, $m, $catchment);
725  if (blessed($i) and $i->isa('Geo::Raster')) {
726  $catchment = $i;
727  $i = shift;
728  $j = shift;
729  $m = shift;
730  } else {
731  $catchment = new Geo::Raster(like=>$fdg);
732  $j = shift;
733  $m = shift;
734  }
735  $m = 1 unless defined($m);
736  my $size = ral_fdg_catchment($fdg->{GRID}, $catchment->{GRID}, $i, $j, $m);
737  return wantarray ? ($catchment, $size) : $catchment;
738 }
739 
740 #** @method Geo::Raster prune(Geo::Raster fdg, Geo::Raster lakes, $min_length, @cell)
741 #
742 # @brief Delete streams that are shorter than min_length in a streams raster.
743 #
744 # Example of removing streams shorter than min_lenght:
745 # @code
746 # $streams->prune($fdg, $lakes, $min_lenght, @cell);
747 # @endcode
748 # @param[in] fdg Flow direction raster.
749 # @param[in] lakes (optional) Lakes raster.
750 # @param[in] min_length (optional) Minimum length (in raster scale!)
751 # of streams to be not deleted. Default is 1.5*cell_size.
752 # @param[in] cell (optional) The root cell (a catchment outlet) from
753 # which the method begins to remove too short streams. If not given
754 # then all streams are pruned.
755 # @return In void context changes this streams raster, otherwise
756 # returns a new streams raster.
757 #*
758 sub prune {
759  my $streams = shift;
760  my $fdg = shift;
761  my $lakes;
762  my $min_length = shift;
763  if (blessed($min_length) and $min_length->isa('Geo::Raster')) {
764  $lakes = $min_length;
765  $min_length = shift;
766  }
767  my $i = shift;
768  my $j = shift;
769  $min_length = 1.5*$streams->cell_size unless defined($min_length);
770  $streams = Geo::Raster->new($streams) if defined wantarray;
771  $i = -1 unless defined $i;
772  if ($lakes) {
773  ral_streams_prune($streams->{GRID}, $fdg->{GRID}, $lakes->{GRID}, $i, $j, $min_length);
774  } else {
775  ral_streams_prune_without_lakes($streams->{GRID}, $fdg->{GRID}, $i, $j, $min_length);
776  }
777  return $streams if defined wantarray;
778 }
779 
780 #** @method Geo::Raster number_streams(Geo::Raster fdg, Geo::Raster lakes, @cell, $id)
781 #
782 # @brief Number streams in a streams raster with unique id.
783 #
784 # @param[in] fdg Flow direction raster.
785 # @param[in] lakes (optional) Lakes raster.
786 # @param[in] cell (optional) The root cell (a catchment outlet) from
787 # which the method begins to treat the streams. If not given
788 # then all stream trees are treated.
789 # @param[in] id (optional) Number for the first found stream, the next
790 # streams will get higher unique numbers. Default is 1.
791 # @return In void context changes this streams raster, otherwise
792 # returns a new streams raster.
793 #*
794 sub number_streams {
795  my $streams = shift;
796  my $fdg = shift;
797  my $lakes;
798  my $i = shift;
799  if (blessed($i) and $i->isa('Geo::Raster')) {
800  $lakes = $i;
801  $i = shift;
802  }
803  my $j = shift;
804  my $sid = shift;
805  $sid = 1 unless defined($sid);
806  $streams = Geo::Raster->new($streams) if defined wantarray;
807  $i = -1 unless defined $i;
808  ral_streams_number($streams->{GRID}, $fdg->{GRID}, $i, $j, $sid);
809  if ($lakes) {
810  $sid = $streams->max() + 1;
811  ral_streams_break($streams->{GRID}, $fdg->{GRID}, $lakes->{GRID}, $sid);
812  }
813  return $streams if defined wantarray;
814 }
815 
816 #** @method Geo::Raster subcatchments(Geo::Raster fdg, Geo::Raster lakes, @cell, $head)
817 #
818 # @brief Divide catchments into subcatchments defined by a streams
819 # raster.
820 #
821 # Example of usage:
822 # @code
823 # $subcatchments = $streams->subcatchments($fdg, @cell);
824 # @endcode
825 # or
826 # @code
827 # ($subcatchments, $topo) = $streams->subcatchments($fdg, $lakes, @cell);
828 # @endcode
829 #
830 # @param[in] fdg The FDG from which the streams raster has been computed.
831 # @param[in] lakes (optional) Lakes raster.
832 # @param[in] cell (optional) The outlet cell of the catchment.
833 # @param[in] head (optional) Boolean value denoting whether the
834 # algorithm should divide the catchment of a headstream to the
835 # catchment of the first cell of the headstream and the catchment
836 # draining to the rest of the headstream. Default is false.
837 # @return Returns a subcatchments raster or a subcatchments raster and
838 # topology, topology is a reference to a hash of
839 # $upstream_element=>$downstream_element associations.
840 #*
841 sub subcatchments {
842  my $streams = shift;
843  my $fdg = shift;
844  my $lakes;
845  my $i = shift;
846  if (blessed($i) and $i->isa('Geo::Raster')) {
847  $lakes = $i;
848  $i = undef;
849  }
850  my $j;
851  my $headwaters;
852  if (@_ == 3) {
853  ($i, $j, $headwaters) = @_;
854  } else {
855  ($headwaters) = @_;
856  }
857  $headwaters = 0 unless defined $headwaters;
858  unless ($lakes) {
859  $lakes = Geo::Raster->new(like=>$streams);
860  $lakes->set(0);
861  }
862  my $subs = Geo::Raster->new(like=>$streams);
863  unless (defined $i) {
864  $i = -1;
865  $j = -1;
866  }
867  my $r = ral_catchment_create($subs->{GRID},
868  $streams->{GRID},
869  $fdg->{GRID},
870  $lakes->{GRID}, $i, $j, $headwaters);
871 
872  # drainage structure:
873  # head -> stream (if exist)<
874  # sub -> lake or stream
875  # lake -> stream
876  # stream -> lake or stream
877 
878  my %ds;
879 
880  for my $key (keys %{$r}) {
881  ($i, $j) = split (/,/, $key);
882  my($i_down, $j_down) = split (/,/, $r->{$key});
883  my $sub = $subs->get($i, $j);
884  my $stream = $streams->get($i, $j);
885  next unless defined $stream;
886  my $lake = $lakes->get($i, $j);
887  my $sub_down = $subs->get($i_down, $j_down);
888  my $stream_down = $streams->get($i_down, $j_down);
889  next unless defined $stream_down;
890  my $lake_down = $lakes->get($i_down, $j_down);
891  if (!defined($lake) or $lake <= 0) {
892  if (defined($lake_down) and $lake_down > 0) {
893  $ds{"sub $sub $i $j"} = "stream $stream";
894  } elsif ($stream != $stream_down or ($i == $i_down and $j == $j_down)) {
895  $ds{"sub $sub $i $j"} = "stream $stream";
896  $ds{"stream $stream $i $j"} = "stream $stream_down";
897  } else {
898  $ds{"head $sub $i $j"} = "stream $stream";
899  }
900  } else {
901  $ds{"sub $sub $i $j"} = "lake $lake";
902  $ds{"lake $lake $i $j"} = "stream $stream_down";
903  }
904  if (defined($lake_down) and $lake_down > 0) {
905  $ds{"stream $stream $i $j"} = "lake $lake_down";
906  }
907  }
908  return wantarray ? ($subs, \%ds) : $subs;
909 }
910 
911 #** @method vectorize_catchment(hashref topology, Geo::Raster streams, Geo::Raster lakes, %params)
912 #
913 # @brief Save the subcatchment structure as a vector layer (a subcatchments raster method).
914 #
915 # @param[in] topology Reference to an hash having as keys = type id i j, and as
916 # values = type id.
917 # @param[in] fdg Flow direction raster.
918 # @param[in] streams Streams raster with unique ids for segments.
919 # @param[in] lakes (required if topology contains lakes) Lakes raster.
920 # @param[in] params Parameters (driver, data_source) for the new Geo::Vectors.
921 # @return Two Geo::Vector objects in an array: ($subcatchments, $streams).
922 #*
923 sub vectorize_catchment {
924  my $self = shift;
925  my $topology = shift;
926  my $fdg = shift;
927  my $streams = shift;
928  my $lakes = shift if ref($_[0]);
929  my %params = @_;
930 
931  my $cell_size = $self->cell_size();
932 
933  my ($minX, $minY, $maxX, $maxY) = $self->world();
934 
935  my %schema = ( Fields => [
936  { Name => 'element', Type => 'Integer' },
937  { Name => 'type', Type => 'String' },
938  { Name => 'down', Type => 'Integer' },
939  { Name => 'type_down', Type => 'String' },
940  ]);
941  my $catchments = $self->polygonize(%params,
942  create => 'subcatchments',
943  schema => \%schema,
944  pixel_value_field => 'element');
945 
946 
947  my %subs;
948  my %streams;
949  my %lakes;
950 
951  # topology, keys = type id i j, values = type id
952  for my $key (keys %$topology) {
953  if ($key =~ /^(\w+) (\d+) \d+ \d+/) {
954  my $type = $1;
955  my $element = $2;
956  if ($topology->{$key} =~ /^(\w+) (\d+)/) {
957  if ($type eq 'stream') {
958  $streams{$element}{type_down} = $1;
959  $streams{$element}{down} = $2;
960  } elsif ($type eq 'lake') {
961  $lakes{$element}{type_down} = $1;
962  $lakes{$element}{down} = $2;
963  } else {
964  #print STDERR "subs $type $1 $2\n";
965  $subs{$element}{type} = $type;
966  $subs{$element}{type_down} = $1;
967  $subs{$element}{down} = $2;
968  }
969  }
970  }
971  }
972 
973  for my $catchment (@{$catchments->features}) {
974 
975  # add subcatchment and lake polygons
976 
977  my $element = $catchment->GetField('element');
978 
979  #print STDERR "$subs{$element} $subs{$element}{type}\n";
980 
981  next unless $subs{$element};
982 
983  $catchment->SetField('type', $subs{$element}{type});
984  $catchment->SetField('down', $subs{$element}{down});
985  $catchment->SetField('type_down', $subs{$element}{type_down});
986  $catchments->feature($catchment->GetFID, $catchment);
987  }
988 
989  %schema = ( Fields => [
990  { Name => 'element', Type => 'Integer' },
991  { Name => 'down', Type => 'Integer' },
992  { Name => 'type_down', Type => 'String' },
993  ]);
994  my $tree = Geo::Vector->new(%params,
995  create => 'stream_segments',
996  schema => \%schema);
997 
998  my %done;
999  my($M, $N) = $streams->size;
1000  for my $i (0..$M-1) {
1001  for my $j (0..$N-1) {
1002  my $s = $streams->get($i, $j);
1003  next unless defined($s);
1004  next if $done{$s};
1005  next unless $streams{$s};
1006  print "segment $s\n";
1007  Gtk2->main_iteration while Gtk2->events_pending;
1008  my $segment = $streams->segment($fdg, $i, $j);
1009  $tree->feature({ Geometry => $segment,
1010  element => $s,
1011  down => $streams{$s}{down},
1012  type_down => $streams{$s}{type_down}
1013  });
1014  $done{$s} = 1;
1015  }
1016  }
1017  return($catchments, $tree);
1018 }
1019 
1020 sub segment {
1021  my($self, $fdg, $i, $j) = @_;
1022  my $s = $self->get($i, $j);
1023  # upstream until not s, Note: assuming a clean segment
1024  while (1) {
1025  $i--;
1026  my $x = $fdg->get($i, $j);
1027  my $s2 = $self->get($i, $j);
1028  next if defined($x) and $x == 5 and defined($s2) and $s2 == $s;
1029  $j++;
1030  $x = $fdg->get($i, $j);
1031  $s2 = $self->get($i, $j);
1032  next if defined($x) and $x == 6 and defined($s2) and $s2 == $s;
1033  $i++;
1034  $x = $fdg->get($i, $j);
1035  $s2 = $self->get($i, $j);
1036  next if defined($x) and $x == 7 and defined($s2) and $s2 == $s;
1037  $i++;
1038  $x = $fdg->get($i, $j);
1039  $s2 = $self->get($i, $j);
1040  next if defined($x) and $x == 8 and defined($s2) and $s2 == $s;
1041  $j--;
1042  $x = $fdg->get($i, $j);
1043  $s2 = $self->get($i, $j);
1044  next if defined($x) and $x == 1 and defined($s2) and $s2 == $s;
1045  $j--;
1046  $x = $fdg->get($i, $j);
1047  $s2 = $self->get($i, $j);
1048  next if defined($x) and $x == 2 and defined($s2) and $s2 == $s;
1049  $i--;
1050  $x = $fdg->get($i, $j);
1051  $s2 = $self->get($i, $j);
1052  next if defined($x) and $x == 3 and defined($s2) and $s2 == $s;
1053  $i--;
1054  $x = $fdg->get($i, $j);
1055  $s2 = $self->get($i, $j);
1056  next if defined($x) and $x == 4 and defined($s2) and $s2 == $s;
1057  $i++;
1058  $j++;
1059  last;
1060  }
1061  # downstream until not s, create the segment
1062  my $segment = Geo::OGC::LineString->new;
1063  while (1) {
1064  $segment->AddPoint(Geo::OGC::Point->new($self->g2w($i, $j)));
1065  my $x = $self->get($i, $j);
1066  last if !defined($x) or $x != $s;
1067  ($i, $j) = $fdg->downstream($i, $j);
1068  last unless defined $i;
1069  }
1070  return $segment;
1071 }
1072 
1073 #** @method route(Geo::Raster dem, Geo::Raster fdg, Geo::Raster k, $r)
1074 #
1075 # @brief Route water downstream (a water state raster method).
1076 #
1077 # Water in each cell of the self raster is routed downstream one time step.
1078 #
1079 # @param[in] dem The DEM.
1080 # @param[in] fdg [optional] The flow directions raster.
1081 # @param[in] k The flow coefficient.
1082 # @param[in] r [optional]. Unit of z divided by the unit of x and y in
1083 # the DEM. Default is 1.
1084 # @todo IN DEVELOPMENT DO NOT USE
1085 # @return The change in water state.
1086 #*
1087 sub route {
1088  my $water = shift;
1089  my $dem = shift;
1090  my $fdg = shift if @_ > 1 and ref($_[1]);
1091  my $k = shift;
1092  my $r = shift;
1093  $r = 1 unless defined $r;
1094  croak ("usage: $water->route($dem, [$fdg,] $k, $r)") unless $k and ref($k);
1095  if ($fdg) {
1096  return Geo::Raster->new(ral_water_route($water->{GRID}, $dem->{GRID}, $fdg->{GRID}, $k->{GRID}, $r));
1097  } else {
1098  return Geo::Raster->new(ral_water_route2($water->{GRID}, $dem->{GRID}, $k->{GRID}, $r));
1099  }
1100 }
1101 
1102 #** @method void vectorize_streams(Geo::Raster fdg, @cell)
1103 #
1104 # @brief Create an OGR layer from a streams raster.
1105 #
1106 # @param fdg The FDG from which the streams raster has been computed.
1107 # @param cell The outlet cell of the catchment.
1108 # @todo IN DEVELOPMENT DO NOT USE
1109 #*
1110 sub vectorize_streams {
1111  my ($self, $fdg, $i, $j, $datasource, $layer) = @_;
1112  ral_streams_vectorize($self->{GRID}, $fdg->{GRID}, $i, $j);
1113 }
1114 
1115 sub compare_dem_derived_ws_attribs {
1116  my ($self, $uag, $dem, $filename, $iname, $ielev, $idarea) = @_;
1117  #my ($self, $filename) = @_;
1118  (my $fileBaseName, my $dirName, my $fileExtension) = fileparse($filename,('\.shp'));
1119  #my $jep = $fileBaseName.$fileExtension;
1120  #print STDERR "$dirName and $jep\n";
1121  ral_compare_dem_derived_ws_attribs($self->{GRID}, $uag->{GRID}, $dem->{GRID}, $dirName, $fileBaseName, $iname, $ielev, $idarea);
1122 }
1123 
1124 1;
int RAL_CALL ral_dem_fill_depressions(ral_grid *dem, ral_grid *fdg)
int RAL_CALL ral_dem_lower_peaks(ral_grid *dem, double z_limit)
ral_grid_handle RAL_CALL ral_fdg_path_sum(ral_grid *fdg, ral_grid *stop, ral_grid *op)
ral_grid_handle RAL_CALL ral_fdg_upslope_sum(ral_grid *fdg, ral_grid *op, int include_self)
ral_catchment_handle RAL_CALL ral_catchment_create(ral_grid *subs, ral_grid *streams, ral_grid *fdg, ral_grid *lakes, ral_cell outlet, int headwaters)
A class for geospatial rasters.
Definition: Algorithms.pm:9
public Geo::Raster new()
Create a new raster.
int RAL_CALL ral_dem_raise_pits(ral_grid *dem, double z_limit)
ral_grid_handle RAL_CALL ral_fdg_depressions(ral_grid *fdg, int inc_m)
public method new()
ral_grid_handle RAL_CALL ral_fdg_path_length(ral_grid *fdg, ral_grid *stop, ral_grid *op)
int RAL_CALL ral_fdg_drain_flat_areas2(ral_grid *fdg, ral_grid *dem)
public hashref contents()
Returns the histogram of an integer raster in a hash.
ral_grid_handle RAL_CALL ral_fdg_upslope_count(ral_grid *fdg, ral_grid *op, int include_self)
private method _new_grid()
int RAL_CALL ral_fdg_drain_depressions(ral_grid *fdg, ral_grid *dem)
ral_grid_handle RAL_CALL ral_fdg_path(ral_grid *fdg, ral_cell c, ral_grid *stop)
public Geo::Raster upslope_count(Geo::Raster mask, scalar include_self)
Compute the count of the upslope cells in a FDG.
int RAL_CALL ral_fdg_kill_extra_outlets(ral_grid *fdg, ral_grid *lakes, ral_grid *uag)
int RAL_CALL ral_streams_prune(ral_grid *streams, ral_grid *fdg, ral_grid *lakes, ral_cell c, double min_l)
ral_grid_handle RAL_CALL ral_water_route(ral_grid *water, ral_grid *dem, ral_grid *fdg, ral_grid *k, double r)
ral_grid_handle RAL_CALL ral_dem_slope(ral_grid *dem, double z_factor)
int RAL_CALL ral_streams_number(ral_grid *streams, ral_grid *fdg, ral_cell c, int sid0)
ral_grid_handle RAL_CALL ral_dem_ucg(ral_grid *dem)
int RAL_CALL ral_dem_breach(ral_grid *dem, ral_grid *fdg, int limit)
ral_cell RAL_CALL ral_fdg_outlet(ral_grid *fdg, ral_grid *streams, ral_cell c)
ral_grid_handle RAL_CALL ral_dem_fdg(ral_grid *dem, int method)
int RAL_CALL ral_streams_break(ral_grid *streams, ral_grid *fdg, ral_grid *lakes, int nsid)
A geospatial layer that consists of Geo::OGR::Features.
Definition: Vector.pm:176
int RAL_CALL ral_dem_fit_surface(ral_grid *dem, double z_factor, ral_grid ***params)
int RAL_CALL ral_fdg_drain_flat_areas1(ral_grid *fdg, ral_grid *dem)
public Geo::Raster drain_flat_areas(Geo::Raster dem, hash params)
Resolve the flow direction for flat areas in a FDG.
A root class for geospatial features.
Definition: Geometry.pm:69
public void set(array cell, scalar value)
Set the value of a cell.