Grid - Perl extension for raster-based GIS work
use Grid;
or
use Grid qw(:types);
or
use Grid qw(:types :logics :db);
Grid is object-oriented to a C library for handling (possibly very large) matrix of values. In some cases Grid adds some very useful functionality to the functionality of the C library. Currently the library supports (short) integers and (4 byte, i.e. float) real numbers.
Each cell in grid is assumed to be a square. The grid point represents the center of the cell.
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
Constructors have either a specified set of parameters in specified order named parameters.
To start with opening a previously saved grid:
$gd = new Grid(filename=>"data/dem");
or simply
$gd = new Grid("data/dem");
The default format for grids is a pair of [bil|dem]/hdr files. The extension of the binary file can be set using named parameter ext. Currently only one-channel 8, 16, or 32 bit images can be read. The datatype of the grid can be specified using named parameter datatype. byte images are by default read into integer grids. If a hdr file is not found the method attempts to read a pair of img/doc (Idrisi (c)) files.
To start with a new grid:
$gd = new Grid(datatype=>$datatype,
M=>100,
N=>100);
or simply
$gd = new Grid(1,100,100);
or even more simply
$gd = new Grid(100,100);
$datatype is optional, the default is $INTEGER_GRID, $REAL_GRID is another possibility. Constants $INTEGER_GRID and $REAL_GRID are imported by :types. Opening a previously saved grid sets the name attribute of the grid.
Other constructors exist, this is a copy:
$g2 = new Grid(copy=>$g1);
or simply
$gd = new Grid($g1);
See below a note about known extensions.
to create a grid with same size use like:
$g2 = new Grid(like=>$g1);
In both copy methods the datatype of the result is the same as in the original grid. Use named parameter datatype=>DATATYPE to upgrade an integer grid to a real grid or downgrade a real grid to an integer grid.
You can also import data:
$gd = new Grid(filename=>$filename, import=>$format);
$format may currently be 'Arc/Info ascii export', 'Arc/Info interchange' or 'ppm'. Arc/Info interchange format may contain other information besides the grid, this is printed out if the option print_header is set. Currently the number of columns in the grid as specified by the interchange file is increased by 3 - why it must do this is unexplained. Datatype may also be set while importing. When importing ppm images you can specify option 'channel' as 'red', 'green', or 'blue', the default is to use the luminance value.
Importing is launched if the filename of the grid to be opened has a recognized extension.
$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 but see below.
Exporting a grid:
$gd->save(filename=>$filename, export=>$method, options=>$options);
Currently only 'arc/info ascii' and 'ppm' have built-in native support for exporting. Filename is optional. The default filename is the name of the grid. The name of the grid is set if the grid is initialized from a disk file or it can be set using the method set_name. This method also sets the name attribute.
If a recognized extension is detected in the filename or in the name of the grid then that will be used as the export method. Recognized extensions are:
asc for arc/info ascii ppm jpeg and jpg (needs ppmtojpeg program) png (needs pnmtopng program)
Options for ppmto... programs can be given as the options argument.
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 Grid(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 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). The extension is used for import and export if it is one the recognized extensions.
Setting and getting the name:
$name = $gd->get_name();
$gd->set_name("dem");
$gd->setbounds(unitdist=>10,
minX=>0,
maxX=>0,
minY=>0,
maxY=>0);
at least three parameters must be set: unitdist, 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 (Grid<->World):
($x, $y) = $gd->g2w($i, $j);
($i, $j) = $gd->w2g($x, $y);
$gd->setmask();
$gd->removemask();
The mask is used in ALL grid operations and affects ALL grids. The removemask method does not need to be called with the same grid as setmask.
$gd->set($i, $j, $x);
$x = $gd->get($i, $j);
the returned value is undef if it is a nodata cell.
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();
($minval, $maxval) = $gd->getminmax();
or
$minval = $gd->min();
$maxval = $gd->max();
If you just want to know where the min or max value resides:
@min = $gd->min();
@max = $gd->max();
@min and @max will be the cell (i,j) of the min or max value.
Methods min and max have quite another meaning if a parameter is supplied, see below.
minvalue and maxvalue are also stored into the grid hash as attributes.
($datatype, $M, $N, $unitdist, $minX, $maxX, $minY, $maxY, $nodata) =
$gd->attrib();
works also in a perlish way:
@size = ($gd->attrib())[1..2];
for this there is also a separate method:
($M, $N) = $gd->size();
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).
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. Grid method round can be used to convert a real grid to an integer grid.
$gd->round();
or
$b = $a->round();
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.
$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 Grid /:logics/;
$b = not($a);
$c = and($a, $b);
$c = or($a, $b);
$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 return the minimum and maximum values of the 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.
$a->if_then($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_then($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).
$g->binary();
This has the same effect as writing $g->ne(0);
$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 = $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 zero (or very near zero in the case of real numbered grids) cells of grid $g the distance to the nearest non-zero cell in $g.
Generating a direction_to grid:
$g->directions();
or
$d = $g->directions();
The returned grid has in the zero (or very near zero in the case of real numbered grids) cells of grid $g the directions to the nearest non-zero 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 non-zero cells is very time consuming.
$g2 = $g1->clip($i1, $j1, $i2, $j2);
$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''.
$g2 = $g1->transform(\@tr, $M, $N);
or, again just (changes the g1 instead of creating a new grid):
$g1->transform(\@tr, $M, $N);
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);
$gd->print(%options);
if no options are given simply prints the grid, if option list=>1 is set prints the nonzero cells of the grid in format:
i,j,val
and returns the points as a reference to an array of references to point arrays ($i, $j, $val). Printing can be suppressed using option quiet=>1, i.e, $g->print(quiet=>1) (quiet implicitly assumes list mode) only returns a reference to an array of points. Other options are 'wc' which changes image coordinates i,j to world coordinates x,y.
$framed_grid = $grid->frame($with);
$with is a scalar and the border cells of $grid is set to that value.
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.
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.
$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.
$g->interpolate(method=>'nn');
or
$a = $b->interpolate(method=>'nn');
This method calculates values for the cells having a zero (or near zero in the case of real valued grids) value from the cells having non-zero values. Currently only the method of using the value of nearest cell ('nn') is implemented.
$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.
Plotting is grayscale unless a colortable is created for the grid;
$gd->colortable(number_of_colors=>number_of_colors,
extra_colors=>extra_colors,
contrast=>contra,
brightness=>brightness,
offset=>offset);
This creates a new colortable (and disposes of the old one). All parameters are optional. The default number of colors is maxvalue - minvalue + 1 for integer grids, for real grids there is no default value (an error is croaked). The default number of extra colors is 0. Contrast is the contrast of the color ramp (default is 1.0). Negative values reverse the direction of the ramp. Brightness is the brightness of the color ramp. Default is 0.5, but can sensibly hold any value between 0.0 and 1.0. Values at or beyond the latter two extremes, saturate the color ramp with the colors of the respective end of the color table. Offset is for those cases when the minimum value in the grid is not zero, the default is the minimum value in the grid. (mostly from pgplot.doc)
$gd->color($i, $l, $r, $g, $b);
or
$gd->color($i, $r, $g, $b);
or
$gd->color($i, colorname);
If colorname is given, the method tries to look up its RGB values from the file /usr/X11/lib/X11/rgb.txt.
This sets one color at index i (first index is 0 or offset if you have set it) in the colortable. l is a normalized ramp-intensity level corresponding to the RGB primary color intensities r, g, and b. l can be set by the method to value $i / ($self->{COLOR_TABLE_SIZE} - 1). Colors on the ramp are linearly interpolated from neighbouring levels. Levels must be sorted in increasing order. 0.0 places a color at the beginning of the ramp. 1.0 places a color at the end of the ramp. Colors outside these limits are legal, but will not be visible if contra=1.0 and bright=0.5. r, g, and b values should be integers between 0 and 255. (mostly from pgplot.doc)
More colors can be added to the colortable using method
$gd->addcolors($newcolors);
Typically colors are retrieved into a grid from the database using the method rgb. See below.
A color can be returned from the colormap using the color method with only the index:
($r, $g, $b) = $gd->color($i);
The whole colortable is returned by the method get_colortable;
$colortable = $gd->get_colortable();
$colortable is a hash with colortable indexes as keys and references to ($r, $g, $b) arrays as values.
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);
The plain plot:
$gd->plot(%options);
device and donotsetminmax and visual can be set using options. device is a pgplot device, the default is '/xserve'. donotsetminmax skips setminmax before plotting. Visual options are labels_off and scale. labels_off skips showing of (vector data) labels. scale skips scaling of x and y-axis. Plot returns right after it has done its job, i.e. it does not remain in interactive state as view:
$gd->view(%options);
The interactive state in view includes export, redraw, scroll, zoom, restore, and area calculation. Using the method drawon adds the capability of drawing lines and rectangles on the viewed grid. The mehod returns the drawn layer as a grid:
$drawing = $gd->drawon(%options);
In the case of interactive methods (view and drawon) the device option can be used to specify the export device, '/ps' is the default.
In all plotting integer grids are first converted to real grids internally.
Opening two pgplot windows:
$w1 = &Grid::gdwindow_open;
$w2 = &Grid::gdwindow_open;
Plotting two grids to these two windows:
$g1->plot(window=>$w1);
$g2->plot(window=>$w1);
The windows are always automatically closed after plotting. Thus you need to reopen them using &Grid::gdwindow_open again if you want to use multiple windows.
An opened window can be closed explicitly by subroutine
&Grid::gdwindow_close($window);
Either plot a grid or close explicitly an opened window.
$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.
$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).
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.
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 pixels 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.
$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.
$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.
$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.
$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.
Grid.pm maintains internally a database handle which can be initialized using the method:
db_connect(\%opt, \%attr);
A table in database can be used to store information about a grid. A database connection must be established before database methods can be used:
use Grid qw(:db);
db_connect({driver=>'driver',
hostname=>'hostname',
port=>'portnumber',
database=>'database',
username=>'username',
password=>'password'});
The default driver is 'Pg' (for PostgreSQL), the default hostname is '', the default port number is 5432 (the default of PostgreSQL), the default database is the current directory name (obtained with 'pwd'), the default username is the current user (obtained with 'whoami'), and the default password is ''.
The database connection is closed with
db_close();
The subroutine
($rows, $fields) = sql($sql, %options);
is a wrapper for running SQL commands. Simple variable substitution is done before the SQL command is excuted. Tokens which resemble Perl variables, e.g. $grid are substituted with respective grid names, i.e, $grid->{NAME}. Usage examples:
use Grid qw(:db);
db_connect();
$lu = new Grid 'lu';
($rows, $fields) = sql('select * from $lu');
print "@{$fields}\n";
foreach (@{$rows}) {print "@{$_}\n"}
another:
$c = $lu->contents();
foreach (keys %{$c}) {
sql('insert into $lu (index,count) values ' . "($_, $$c{$_})");
}
Note the use of '' and ``''. These examples assume that there is a table 'lu' in the database and that there are columns index and count in that table.
If the table has integer columns r, g, and b we can use the method:
$lu->rgb(extra_colors=>extra_colors);
to create a colortable into grid lu based on the r, g, and b values in the table lu in the database. The r, g, and b should have values in the range 0..255. extra_colors is the number of extra colors in the color table and its default is 0.
&Grid::db_initvd();
NOTE: The structure of the vector database is subject to change. See the Shape.pm project.
Vector data can be retrieved from the database and attached to a grid using the method:
$g->addlayer($layer);
$layer can be either the name or the id of the layer (this is checked using a regexp match /^\d+$/). This method looks into the database for a table layers. The table layer has to have columns layer (of type text) and id (of type int). If a layer is found, the geometric types which belong to the layer are retrieved next. Currently supported geometric types are point (table points), lseg (table lsegs), and polygon (table polygons). Each of these tables must exist and have at least columns data (of respective type point, lseg, or polygon), r (of type int), g (of type int), b (of type int), id (of type text), lid (of type int). r, g, and b define the color of the geometric element, id is the, also visible, name of the element, and lid is the id of the layer into which the element belongs to. The method looks up a suitable color index for the color and if no suitable one exist, adds the color into the colortable of the grid.
Several layers can be attached into a grid. These will be organized into a list where the last added layer is the last. Methods for organizing the layers will be added later. All layers can be removed using the method:
$g->removelayers();
These methods require an external program 'voronoi'. The mehods create vector layers into the database.
$gd->voronoi($layer);
$layer is optional, default is 'voronoi'.
$gd->delaunay($layer);
$layer is optional, default is 'delaunay'.
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.
$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.
$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.
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.
$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.
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.
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 (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);
(This method checks the internal FDG flag of the object on which the method is applied, the FDG flag is set when it is created using the fdg method).
The FDG should be calculated from the DEM and have no pits or flat area cells (the FDG returned by the filldepressions is ok). If it is not given it is calculated using standard methods.
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.
Usage:
$d = $fdg->distance_to_channel($open_water_grid)
$d = $fdg->distance_to_divide()
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.
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.
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
DInfinity grids can be made but otherwise the methods to handle them do not work.
Ari Jolma, ajolma@water.hut.fi
perl(1).