Ari Jolma > Geo-Raster-0.42 > Geo::Raster

Download:
Geo-Raster-0.42.tar.gz

Dependencies

Annotate this POD

CPAN RT

New  1
Open  0
View/Report Bugs
Module Version: 0.42   Source  

NAME ^

Geo::Raster - Perl extension for raster algebra

SYNOPSIS ^

    use Geo::Raster;
or
    use Geo::Raster qw(:types);
or
    use Geo::Raster qw(:types :logics :db);

DESCRIPTION ^

Geo::Raster is an object-oriented interface to libral, a C library for rasters and raster algebra. Geo::Raster makes using libral easy and adds some very useful functionality to it. libral rasters are in-memory for fast and easy processing. libral rasters can be created from GDAL rasters. GDAL provides access to rasters in many formats.

Geo::Raster also adds the required functionality to display rasters in Gtk2::Ex::Geo.

Each cell in raster/grid is assumed to be a square.

The grid point represents the center of the cell and not the area of the cell (when such distinction needs to be made). TODO: This needs more attention.

A grid is indexed like this:

                   j = 0..N-1  
               ------------------>
              .
i = 0..M-1    .
              .
              .
              V

there is also a (x,y) world coordinate system

        maxY  ^
              .
              .
     y        .
        minY  .
               ------------------>
               minX           maxX
                              x

minX is the left edge of first cell in line. maxX is the right edge of the last cell in a line. minY and maxY represent similarly the boundaries of the raster.

BASIC FUNCTIONALITY ^

Constructors

Constructors have either a specified set of parameters in specified order named parameters.

To start with opening a previously saved grid:

    $gd = new Geo::Raster(filename=>"data/dem", load=>1);

or simply

    $gd = new Geo::Raster("data/dem");

A Geo::Raster object maintains internally a gdal object. To get data from the gdal object into the grid, use the load option (load is not assumed if only filename is given) or cache method (see below).

To start with a new grid:

    $gd = new Geo::Raster(datatype=>datatype_string, M=>100, N=>100);

or simply

    $gd = new Geo::Raster(1, 100, 100);

or even more simply

    $gd = new Geo::Raster(100, 100);

$datatype is optional, the default is 'integer', 'real' is another possibility (actual types of integer and real are defined in libral). Opening a previously saved grid sets the name attribute of the grid.

Other constructors exist, this is a copy:

    $g2 = new Geo::Raster(copy=>$g1);

or simply

    $gd = new Geo::Raster($g1); 

See below a note about known extensions.

to create a grid with same size use like:

    $g2 = new Geo::Raster(like=>$g1);

In both copy methods the datatype of the result is the same as in the original grid. Use named parameter datatype to upgrade an integer grid to a real grid or downgrade a real grid to an integer grid.

Caching data into the work raster from the dataset

Read all data:

    $gd->cache();

Use an existing raster as a model (bounding box, cell_size):

    $gd->cache($like_this_grid); 

Use a bounding box:

    $gd->cache($minX,$minY,$maxX,$maxY);
    $gd->cache($minX,$minY,$maxX,$maxY,$cell_size);

If the cell_size is not given, the cell_size of the dataset is used.

If the cell_size is specified, it is used if it is larger than the cell_size of the dataset.

The given bounding box clipped to the bounding box of the dataset. The resulting bounding box of the work raster is always adjusted to pixel boundaries of the dataset.

Saving a grid:

    $gd->save("data/dem");

If no filename is given the method tries to use the name attribute of the grid.

The default saving format is a pair of hdr/bil files.

Dump and Restore

Some grids may contain very little information, then dumping and restoring is an option which saves disk space. Note, however, that dumping does not save the size or other attributes of the grid:

$g->dump($to);

Dump is in fact grid print method (see below) using options quiet and nonzeros with a redirect to $to. Note that restore is not a constructor, so you may need to create the grid before restoring it:

$g = new Geo::Raster(like=>$some_other_grid);

$g->restore($from);

$to and $from can be undef (then dump uses STDOUT and restore uses STDIN), filename, or reference to a filehandle (e.g., \*DUMP).

The name of the grid

The name of the grid is (or may be) the same as its filename (including the path) without extension. The name is set when a grid is constructed from a file (or files).

Setting and getting the name:

    $name = $gd->get_name();
    $gd->set_name("dem");

Setting the world coordinate system:

    $gd->setbounds(cell_size=>1,
                   minX=>0,
                   minY=>0, 
                   maxX=>10,
                   maxY=>10);

at least three parameters must be set: cell_size, minX and minY; minX, maxX and minY; or minX, minY and maxY. minX (or easting) is the left edge of the leftmost cell, i.e. _not_ the center of the leftmost cell.

The world coordinate system can be copied to another grid:

    $g1->copyboundsto($g2);

Conversions between coordinate systems (Cell<->World):

    ($x, $y) = $gd->g2w($i, $j);
    ($i, $j) = $gd->w2g($x, $y);

Setting and removing a mask

    $gd->setmask();

    $gd->removemask();

The mask is used in ALL grid operations made on _this_ grids.

Setting a cell value

    $gd->set($i, $j, $x);

If $x is undefined or string "nodata", the cell value is set to nodata_value.

Setting all cells to a value

    $gd->set($x);

If $x is undefined or string "nodata", the cell value is set to nodata_value.

Copying values from another grid

    $gd->set($g);

$g needs to be a similar Geo::Raster.

The return value is 0 in the case of an error.

Retrieving a cell's value

    $x = $gd->get($i, $j);

The returned value is either a number, 'nodata', or undef if the cell is not in the grid.

The same in world coordinates

    $x = $gd->wget($x, $y);

Nodata values

Both integer and floating point grids may have nodata values. The nodata value is a special integer or real value, which is stored in the grid data structure. Nodata values behave a bit like a mask, they are not used in queries and in arithmetics operands do not affect them. For example if c is nodata value then after

a = b + c

a is nodata.

and

b += c

b is nodata.

Also in comparisons the result is nodata if either of the values to be compared is nodata.

In the construct "if a then b = c" for grids. b is assigned c only if a is data and c is data.

The method "data", which can be used as in-place or as normal method which returns a value and does not affect the grid itself, returns a binary grid, which has 0 where there were nodata values and 1 where there were data.

    $gd->data();

or

    $b = $a->data();

Calculating min and max values

    ($minval, $maxval) = $gd->getminmax(); 

or

    $minval = $gd->min();
    $maxval = $gd->max();

These methods have quite another meaning if a parameter is supplied, see below.

Retrieving the attributes of a grid (deprecated)

    ($datatype, $M, $N, $cell_size, $minX, $minY, $maxX, $maxY, $nodata_value) = 
    $gd->attributes();

Use the specific methods instead:

    $datatype = $gd->datatype(); # returns a string 

    ($M, $N) = $gd->size();

    $cell_size = $gd->cell_size();

    ($minX,$minY,$maxX,$maxY) = $gd->world();

    $nodata_value = $gd->nodata_value();

Size is interpreted as a size of a zone, if a cell is given as a parameter to size method:

    $zone_size = $gd->size($i,$j);

Arithmetics

Add/subtract a scalar or another grid to a grid:

    $b = $a + $x; # is equal to $b = $a->plus($x);
    $b = $a - $x; # is equal to $b = $a->minus($x);

In-place versions

    $gd += $x; # is equal to $gd->add($x);
    $gd -= $x; # is equal to $gd->subtract($x);

Multiply/divide the grid by a scalar or by another grid:

    $b = $a * $x; # is equal to $b = $a->times($x);
    $b = $a / $x; # is equal to $b = $a->over($x);

In-place versions

    $gd *= $x; # is equal to $gd->multiply_by($x);
    $gd /= $x; # is equal to $gd->divide_by($x);

NOTE: THIS IS NOT MATRIX MULTIPLICATION: what goes on is, e.g.:

  for all i,j: b[i,j] = a[i,j] * x[i,j]

Modulus:

    $b = $a % $x; # is equal to $b = $a->modulo($x);
    $gd %= $x;    # is equal to $gd->modulus_with($x);

Power:

    $b = $a**$x;  # is equal to $b = $a->power($x);
    $gd **= $x;   # is equal to $gd->to_power_of($x);

DO NOT use void context algebraic operations like $a + 5; The effect is not what you expect and it will generate a warning if run with the -w switch.

Integer grids are silently converted to real grids if the operand is a real number or a real grid or if the operator is "/" (except in modulus, which is defined only for integer grids).

Mathematical operations

Integer grids are silently converted to real grids if these methods are applied. The only exception is abs, which is defined for integer grids:

    $b = $a->abs();
    $b = $a->acos();
    $b = $a->atan();
    $c = $a->atan2($b);
    $b = $a->ceil();
    $b = $a->cos();
    $b = $a->cosh();
    $b = $a->exp();
    $b = $a->floor();
    $b = $a->log();
    $b = $a->log10();
    $b = $a->sin();
    $b = $a->sinh();
    $b = $a->sqrt();
    $b = $a->tan();
    $b = $a->tanh();

abs, atan2, cos, exp, log, sin, sqrt are overloaded

In-place versions (use always methods for in-place versions) change the original grid:

    $a->abs();

...etc.

If $a is not a grid, the functions fall back to standard Perl math functions.

NOTE: ceil and floor are defined only for real grids and return a real grid. Geo::Raster method round can be used to convert a real grid to an integer grid.

    $gd->round();

or

    $b = $a->round();

Comparisons between grids

Comparison of grid to a scalar or to another grid:

    $g2 = $g1 op $x;

where op is "<", ">", "<=", ">=", "==", "!=", or "<=>". $x may be a scalar or another grid. The return value is always an integer grid. For in-place versions of the comparisons use the methods lt, gt, le, ge, eq, ne, and cmp.

So there are four cases of the use of comparison operations:

                    a unchanged
 1. b = a->lt(0);      yes     
 2. a->lt(0);          no      
 3. b = a < 0;         yes     
 4. b = 0 < a;         yes     

DO NOT use void context comparisons like $a < 0; The effect is not what you expect and it will generate a warning if run with the -w switch.

Logical operations

    $b = $a->not();
    $c = $a->and($b);
    $c = $a->or($b);

in-place versions (changes a)

    $a->not();
    $a->and($b);
    $a->or($b);

or

use Geo::Raster /:logics/;

    $b = not($a);
    $c = and($a, $b);
    $c = or($a, $b);

Minimum and maximum

    $g2 = $g1->min($x);
    $g2 = $g1->max($x);

again, in-place versions also work.

The effect is (for the method "min"):

g2[i,j] = min( g1[i,j] , x[i,j] ) or, if x is scalar, g2[i,j] = min( g1[i,j] , x ).

If $x is undef these methods refer to the minimum and maximum values of the grid. In scalar context the methods return the minimum value.

Cross product of grids

Cross product of two grids is defined for integer grids only.

    $c = $a->cross($b);

or in-place

    $a->cross($b);

c = a x b

If a has values a1, ..., ana (ai < aj, na distinct values) and b has values b1, ..., bnb (bi < bj, nb distinct values) then c will have nc = na * nb distinct values 1, ..., nc. The c will have value 1 where a = a1 and b = b1, 2 where a = a1 and b = b2, etc.

if ... then construct for grids

    $a->if($b, $c);

$a and $b are grids and $c can be a grid or a scalar, the effect of this subroutine is:

for all i,j if (b[i,j]) then a[i,j]=c[i,j]

if a return value is requested

    $d = $a->if($b, $c);

then d is a but if b then c

If $c is a reference to a zonal mapping hash, i.e., it has value pairs k=>v, where k is an integer, which represents a zone in b, then a is set to v on that zone. A zone mapping hash can, for example, be obtained using the zonal functions (see below).

Convert an integer image into a binary image:

    $g->binary();

This has the same effect as writing $g->ne(0);

Bufferzone

    $g2 = $g1->bufferzone($z, $w);

Creates (or converts a grid to) a binary grid g2, where all cells within distance w of a cell (measured as cell center to cell center) in g1 having value z will have value 1, all other cells in g2 will have values 0.

g1 has to be an integer grid. The return value is optional.

count, sum, mean, ...

    $count = $gd->count();
    $sum = $gd->sum();
    $mean = $gd->mean();
    $variance = $gd->variance();

Return scalars. All cells except those having nodata values are taken into account.

Similar zonal functions are below. Min and max functions are above.

Generating a proximity grid:

    $g->distances();

or

    $d = $g->distances();

The returned grid has in the nodata cells of grid $g the distance to the nearest data cell in $g.

Generating a direction_to grid:

    $g->directions();

or

    $d = $g->directions();

The returned grid has in the nodata cells of grid $g the directions to the nearest data cell in $g. Directions are given in radians and direction zero is to the direction of x-axis, Pi/2 is to the direction of y-axis.

NOTE: The generation of proximity and direction_to grids in the case of large grids (more than 100 000 cells) and few data cells is very timeconsuming.

Clipping a grid:

    $g2 = $g1->clip($i1, $j1, $i2, $j2);

or

    $g2 = $g1->clip($g3);

to clip from $g1 a piece which is overlayable with $g3.

If there is no lvalue, $g1 is clipped.

Joining two grids: (NOTE: this is from before gdal/cache)

    $g3 = $g1->join($g2);

The joining is based on the world coordinates of the grids. clip and join without assignment clip or join the original grid, so

    $a->clip($i1, $j1, $i2, $j2);
    $a->join($b);

have the effect "clip a to i1, j1, i2, j2" and "join b to a".

Transforming a grid

    $g2 = $g1->transform(\@tr, $M, $N, $pick, $value);

or, again just (changes the g1 instead of creating a new grid):

    $g1->transform(\@tr, $M, $N, $pick, $value);

g2 will be of size M,N, transformation uses eqs:

  i1 = ai + bi * i2 + ci * j2
  j1 = aj + bj * i2 + cj * j2

whose parameters are in array @tr:

    @tr = (ai, bi, ci, aj, bj, cj);

$pick and $value are optional, $pick may be "mean", "variance", "min", "max", or "count". If $pick is "count" then $value should be the value which needs to be counted -- this works only for integer grids. Result grid is the same type as input except for mean and variance which are always floats. Division by n-1 is used for calculating variance.

In the case when $pick is not defined, the value which is stored into the target grid is looked up from the source grid using the equations above, and rounding the indexes to the nearest integer value.

In the case when $pick is defined, the value which is stored into the target grid is calculated from the (possibly rectangular) area into which the i2,j2 cell maps to. NOTE: In this case the cell coordinates are assumed to denote the upper left corner of the cell. This makes it easy to keep the (x,y) of the upper left the same BUT it is different than the usual assumption that (i,j) denotes the center of the cell.

The print method

    $gd->print();

simply prints the grid to stdout.

Making an array of data in a grid

    $aref = $gd->array;

The $aref is a reference to an array of cells and values:

(i0,j0,val0,i1,j1,val1,i2,j2,val2,i3,j3,val3,...).

Framing a grid

    $framed_grid = $grid->frame($with);

$with is a scalar and the border cells of $grid is set to that value.

Getting an overview of the contents of a grid

Histogram may be calculated using the method "histogram":

    $histogram = $gd->histogram(\@bin);

which returns a reference to an array which holds the counts of cells in each bin, first bin is values <= $bin[0], second bin is values > $bin[0] and values <= $bin[1], etc.

If the parameter is an integer value, the bins are internally created by splitting the range [minval..maxval] into that many parts and the return value is a hash where the key is the center value of the bin (suitable for gnuplot histeps plotting style).

A simpler tool that resembles histogram is contents:

    $contents = $gd->contents();

which returns a reference to a hash which has, values as keys and counts as values. This works for both integer and floating point grids.

Zonal functions

All zonal functions require two grids: the operand grid and the zones grid. The operand grid may be any grid. The zones grid has to be an integer grid. The zonal functions all return a hash, where the keys are the integers from the zones grid (not nodata but 0 yes). The values in the hash are either all the values (nodata values skipped) from the zone (as a reference to an array) or some function (count, sum, min, max, mean, variance) of them. The method which returns all the zone data may of course be used to calculate whatever function but this can take a lot of memory and computing time in the case of large grids. Division by n-1 is used for calculating variance.

    $zh = $gd->zones($zones);

    $counts = $gd->zonalcount($zones);
    $sums = $gd->zonalsum($zones);
    $mins = $gd->zonalmin($zones);
    $maxs = $gd->zonalmax($zones);
    $means = $gd->zonalmean($zones);
    $variances = $gd->zonalvariance($zones);

The zones grid can be changed using the method

    $zones->growzones($grow);

or

    $new_zones = $zones->growzones($grow);

which "grows" each zone in the zones grid iteratively to areas designated by the (binary) grid grow.

Note that also zero zone is also a zone. Only nodata areas are not zoned.

Interpolation

    $g->interpolate(method=>"nn");

or

    $a = $b->interpolate(method=>"nn");

This method interpolates values for all nodata cells. Currently only the method of using the value of nearest cell ("nn") is implemented.

Filling a grid using an arbitrary function of x and y

    $g->function("<function of x and y>");

fills the grid by calculating the z value for each grid cell separately using the world coordinates. An example of a function string is '2*$x+3*$y', which creates a plane.

GRAPHICS ^

Primitives:

Drawing a line to a grid:

    $gd->line($i1, $j1, $i2, $j2, $pen);

a filled rectangle:

    $gd->rect($i1, $j1, $i2, $j2, $pen);

a circle:

    $gd->circle($i, $j, $r, $pen);

Without $pen these change into "extended" get. Then the method returns an array (i,j,value, i,j,value, ...) of all values under the line, rect or circle.

IMAGE MANIPULATION METHODS ^

Mapping values

    $img2 = $img1->map(\%map);

or

    $img->map(\%map);

or, for example, using an anonymous hash created on the fly

    $img->map({1=>5,2=>3});

Maps cell values (keys in map) in img1 to respective values in map in img2 or within img. Works only for integer grids.

Hint: Take the contents of a grid, manipulate it and then feed it to the map.

Creating a colored map

    $colored_map = $img->colored_map();

or

    $img->colored_map();

Uses the least possible number (?) of unique colors to color the map (= image which consists of areas).

applytempl

The "apply template" method is a generic method which is, e.g., used in the thinning algorithm below.

    $gd->applytempl(\@templ, $new_val);

applytempl which takes a structuring template and a new value as parameters. A structuring template is an integer array [0..8] where 0 and 1 mean a binary value and -1 is don't care. The array is the 3x3 neighborhood:

0 1 2 3 4 5 6 7 8

The cell 4 is the center of the template. If the template matches a cell's neighborhood, the cell will get the given new value after all cells are tested. The grid on which applytemp is used should be a binary grid. In void context the method changes the grid, otherwise the method returns a new grid.

thin

Thinning:

    $thinned_img = $img->thin(%options);

or

    $img->thin(%options);

This is an implementation of the algorithm in Jang, B-K., Chin, R.T. 1990. Analysis of Thinning Algorithms Using Mathematical Morphology. IEEE Trans. Pattern Analysis and Machine Intelligence. 12(6). 541-551. (Same as in Grass but done in a bit different, and more generic way, I believe). Options are algorithm=>, trimming=>, maxiterations=>, and width=>. Algorithm is by default "B" (other option is "A"), trimming is by default 0 (other option is 1), and maxiterations is by default 0 (no maximum, will iterate until no cells are deleted), if width is used, maxiterations is set to int(width/2). Trimming removes artificial branches which grow on the side of wide lines in thinnning but it also shortens a bit the real branches. The thinned grid must be a binary grid.

The thinning algorithm defines a set of structuring templates and applies them in several passes until there are no matches or until the maxiterations is reached. Trimming means certain structuring templates are applied to kill emerging short limbs which appear because of the noise in the grid.

Borderizing

    $borders_img = $img->borders(method=>simple|recursive);

or

    $img->borders(method=>simple|recursive);

The default method is "recursive" which finds all areas (8-connected cells with non-zero values) and marks their borders with respective values and leaves the rest of the area zero.

Find the areas in an image

    $areas_in_img = $img->areas($k);

or

    $img->areas($k);

The $k (3 if not given) is the number of consecutive non-zero 8-neighbors required before the cell is assumed to be a part of an area.

Connect broken lines

    $img2 = $img1->connect();

or

    $img->connect();

If two 8-neighbor opposite cells (1-5, 2-6, etc) of a cell are the same and the cell is zero, then the value of this cell is set to the same value.

Number areas with unique id in an image

    $map = $img->number_areas($connectivity);

or

    $img->number_areas($connectivity);

$map or $img contains each consecutive area in $img colored with an unique integer. $img should be a 0,1 grid, the result is a 0,2,3,... grid. The connectivity of areas is either 8 (default) or 4.

TERRAIN ANALYSIS METHODS ^

METHODS FOR DIGITAL ELEVATION MODELS (DEMS)

Slope and aspect

Generate an aspect grid from a DEM:

    $aspect = $dem->aspect();

or

    $dem->aspect(); # to convert the DEM to an aspect grid

The returned aspect grid is a real number grid (-1,0..2*Pi) where -1 denotes flat area, 0 aspect north, Pi/2 aspect east, etc.

Similar method exist for calculating a slope grid from a DEM:

    $slope = $dem->slope($z_factor);
    $dem->slope($z_factor); # to convert the DEM to a slope grid

Slope and aspect calculations are based on fitting a 9-term quadratic polynomial:

z = Ax^2y^2 + Bx^2y + Cxy^2 + Dx^2 + Ey^2 + Fxy + Gx + Hy + I

to a 3*3 square grid. See Moore et al. 1991. Hydrol. Proc. 5, 3-30.

Slope is calculated in radians.

z_factor is the unit of z dived by the unit of x and y, the default value of z_factor is 1.

Flow direction grid (FDG) from DEM

    $fdg = $dem->fdg(method=>$method);

or

    $dem->fdg(method=>$method); # to convert the DEM to a FDG 

The method is optional. The default method is D8 (deterministic eight-neighbors steepest descent) and the returned FDG is of type D8, i.e., an integer grid (-1..8) where -1 denotes flat area, 0 a pit, 1 flow direction north, 2 north-east, etc. A pit is the lowest point of a depression. Another supported method is Rho8 (stochastic eight-neighbors aspect-based) which also produces a D8 FDG but the direction is chosen between the two steepest descent directions (assumed to being next to each other) so that the expected direction is the true aspect (Fairfield and Leymarie, Water Resour. Res. 27(5) 709-717). The third method is "many" which produces a FDG, where the bits in each byte (actually a short integer) in each cell denotes the neighbors having lower elevation, i.e., value 1 (2**0 = 1) means only the cell in direction 1 is lower, value 3 (2**0+2**1 = 3) means both cells in direction 1 and in direction 2 are lower, etc.

METHODS FOR D8 FDGS

General movecell method for 8-neighborhood

    $gd->movecell(@cell, $dir);

this returns undef if the cell moves outside of the grid.

If $dir is not given we assume this grid to be a FDG and then the dir is in the grid.

Upstream cells

A method for flow direction grid to get the directions of upstream cells of a cell:

    ($up,@up) = $fdg->upstream($streams,@cell);

or

    (@up) = $fdg->upstream(@cell);

$up is the direction of upstream stream cell and @up contain directions of other upstream cells.

Draining flat areas

    $fixed_fdg = $fdg->fixflats($dem,method=>$method);

or

    $fdg->fixflats($dem,method=>$method);

The method is either "one pour point" (short "o") or "multiple pour points" (short "m"). The first method, which is the default, finds the lowest or nodata cell just outside the flat area and, if the cell is lower than the flat area or a nodatacell, drains the whole area there, or into the flat area cell (which is made a pit cell), which is next to the cell in question. This method is guaranteed to produce a FDG without flat areas. The second method drains the flat area cells iteratively into their lowest non-higher neighboring cells having flow direction resolved.

Handling the depressions in a DEM

Methods

    $dem->fill($z_limit);

and

    $dem->cut($z_limit);

raise or lower cells which are lower or higher than all its 8-neigbors. The z_limit is the minimum elevation difference, which is needed to consider a cell lower or higher than all its neighbors. $z_limit is optional, the deafult value is 0.

A depression (or a "pit") is a connected (in the FDG sense) area in the DEM, which is lower than all its neighbors. To find and look at all the depressions use this method, which returns a grid:

    $depressions = $dem->depressions($fdg, $inc_m);

The argument $fdg is optional, the default is to calculate it using the D8 method and then route flow through flat areas using the methods "multiple pour points" and "one pour point" (in this order). The depressions grid is a binary grid unless $inc_m is given and is 1.

Depressions may be removed by filling or by breaching. Filling means raising the depression cells to the elevation of the lowest lying cell just outside the depression. Breaching means lowering the elevation of the "dam" cells. The breaching is tried at the lowest cell on the rim of the depression which has the steepest descent away from the depression (if there are more than one lowest cells) and the steepest descent into the depression (if there are more than one lowest cells with identical slope out) (see Martz, L.W. and Garbrecht, J. 1998. The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrol. Process. 12, 843-855; the breaching algorithm implemented here is close to but not the same as theirs - the biggest difference being that the depression cells are not raised here). Breaching is often limited to a certain number of cells. Both of these methods change the DEM. Both methods need to be run iteratively to remove all removable depressions. Only the filling method is guaranteed to produce a depressionless DEM.

The non-iterative versions of the methods are:

    $dem->filldepressions($fdg);

and

    $dem->breach($fdg, $limit);

The $limit in breaching is optional, the default is to not limit the breaching ($limit == 0). The $fdg, which is given to these algorithms should not contain flat areas.

If the $fdg is not given it is calculated as above in the depressions method and the depressions are removed iteratively until all depressions are removed or the number of depressions does not diminish in one iteration loop.

A method, which produces a pitless FDG is

    $fdg = $dem->pitless_fdg();

This method is similar to the above methods but it does not change the DEM. It changes the path in the FDG from the bottom of the pit to the lowest pour point of the depression. The method is also iterative as the above methods. It first computes a FDG without flat areas and then applies the method $fdg->fixpits($dem) until there are not more pits or there is no change.

Routing of water

Method

    $water->route($dem, $fdg, $flow, $k, $d, $f, $r);

routes water out from a catchment. The method is recursive and routes water from each cell downslope if water from all its upslope cells have been routed downslope.

The catchment tree is traversed using the flow direction grid, which thus must contain only valid directions (no pits nor flat area cells).

The flow from cell a to a downstream cell b is calculated using eq:

    slope = r * (h(a) - h(b)) / (UNIT_DISTANCE * distance_unit(dir(a->b)))
    flow = k * (slope + d) * water(a)

    r               is the unit of z dived by the unit of x and y, e.g, 
                    if z is given in cm and UNIT_DISTANCE = 25 m, then 
                    r = 1 cm / 1 m = 0.01. $r is by default 1

    h(x)            is the elevation of x

    dir(a->b)       is the direction from a to b

    UNIT_DISTANCE   is a property of the DEM 

    distance_unit() is 1 if direction is north, east, ... and sqrt(2) if
                    direction is north-east, south-east, ...  
    
    k               is a parameter

    d               is a parameter

    water(a)        is the amount of water at cell a

Arguments:

    $water Storage at each cell [grid]
    $dem   DEM (input) [grid]
    $fdg   FDG (input) [grid]
    $flow  Amount of water leaving each cell (output) [grid]
    $k     parameter [grid]
    $d     parameter [grid]
    $f     determines if water is routed from each cell to all of its
           neighbors having the same or lower elevation ($f == 1) or
           to the cell pointed by FDG ($f == 0) (default 1)
    $r     is the unit of z dived by the unit of x and y (default 1)

Upslope area

Upslope area (as a number of cells) grid (UAG) from D8 FDG:

    $uag = $fdg->uag();

upslope area (as a cell area) directly from a depressionless DEM

    $uag = $dem->uag(fdg=>$fdg);

The FDG should be calculated from the DEM and have no pits or flat area cells (the FDG returned by the filldepressions is ok).

UAG is a real number grid.

From a 8-direction FDG one can get the upslope area (catchment) of a cell:

    $catchment = $fdg->catchment($i, $j, $m);

or to mark on an existing grid

    $fdg->catchment($catchment, $i, $j, $m);

The method marks the catchment with $m. $m is not required 1 is used by default. In an array context the returned array is ($catchment, $size) where the $size is the size of the catchment.

Distance along flow path to open water

Usage:

$d = $fdg->distance_to_channel($open_water_grid,[$steps])

Description:

Returns a new (real-numbered) grid whose cell values represent the distance to nearest channel (non-zero value in streams grid) along the flow path (defined by flow direction grid fdg). The distance is in the grid units. If steps is given and non-zero, the returned integer grid measures the distance in steps from pixel to pixel.

METHODS FOR STREAMS GRIDS

Streams grid may be obtained from the upslope-area grid by thresholding. If it is to be used with a lakes grid in the subcatchment method, it should be elaborated using methods:

    $streams->prune($fdg, $lakes, $i, $j, $l);

which removes streams shorter than $l (in grid scale), note: also streams which end in a lake may be removed. If $l is not given the method removes one pixel streams. The lakes grid is optional and can be left out.

    $streams->number_streams($fdg, $lakes, $i, $j);

Gives a unique id for each stream section in a stream-tree, which root is at (i,j). The lakes grid is optional and can be left out.

Generation of a subcatchment grid

Subcatchment grid shows all subcatchments defined by a stream network.

    $subcatchments = $streams->subcatchments($fdg, $i, $j);

or

    ($subcatchments, $topo) = 
        $streams->subcatchments($fdg, $lakes, $i, $j);

where the i,j is the outlet point of the whole catchment. $topo is the topology of the catchment as a hash of associations:

$upstream_element=>$downstream_element

BUGS ^

DInfinity grids can be made but otherwise the methods to handle them do not work.

SEE ALSO ^

gdalconst, gdal

This module should be discussed in geo-perl@list.hut.fi.

The homepage of this module is http://libral.sf.net.

AUTHOR ^

Ari Jolma, ari.jolma _at_ tkk.fi

COPYRIGHT AND LICENSE ^

Copyright (C) 1999-2006 by Ari Jolma

This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself, either Perl version 5.8.5 or, at your option, any later version of Perl 5 you may have available.

syntax highlighting: