Module Algebra

opals::IAlgebra

# Aim of module

Provides a generic raster algebra calculator featuring mathematical, statistical and conditional operators.

# General description

Module Algebra is a powerful tool as it allows the combination of multiple input grid datasets, each potentiallz containing multiple bands, by means of a user defined algebraic formula. Basically, It can be used either as a selector, as a calculator or, in all generality, as a mixture of both. Module Algebra can be applied to derive, e.g.:

• difference models: by subtracting two grid models (->calculator)
• grid masks: by evaluating one or more grid models and returning either true (transparent) or false (opaque) as result (->calculator)
• erosion hazard models: by combining the pixels of multiple grids in an algebraic way (e.g. Wischmeier formula, ->calculator)
• landscape dependent DSM models: by choosing the output pixel values from multiple input grid sources according to certain criteria (->selector)
• DTM mosaics: by combining multiple input grids (->selector)
• Orthophoto mosaics: by combining the corresponding bands of multiple input grids (->selector)
• point density models: by summing up the point densities of multiple input grids (->selector+calculator)
• ...

In general, the output grid (parameter outFile) is derived by applying an algebraic formula (parameter formula) combining the values of one ore more input grid models (parameter inFile). The input grids may contain one or more raster bands. A mixture of single- and multiband input is explicitely allowed. The formula can either be specified as a generic filter or in python syntax. In the latter case the python interpreter is used evaluate the formula.

By default, the output grid is created with the same pixel size as the first input grid, but an arbitrary pixel size can be specified using parameter gridSize. The extents of the resulting grid can either be specified by the user (parameter limit (left,lower,right,upper)) or they are calculated from of the outermost pixel centers of all input grid datasets involved. If no user defined limits are specified, either the common area (parameter limit intersect), the overall area (parameter limit union) or the extents of the i'th dataset (parameter limit dataset:i) are applied (c.f. section Parameter description for more details).

Module Algebra can deal with non-coinciding grid/raster models simultaneously. If necessary, the input grids/rasters are re-sampled according to the output grid structure prior to evaluating the formula (parameter resampling). To strictly avoid resampling, please specify resampling=never. This would lead to a program termination if any of the contributing input grid files would require resampling.

One of the core applications of Module Algebra is mosaicing. In this case multiple input datasets are spatially disjoint so that, in general, only a single raster value is available on a specific location (opposed to apllication scencarios with overlapping datasets). To speed up processing, parameter skipVoidAreas will pre-select all relevant input datasets when processing a specific local region (i.e. processing unit) and will exlude all datasets which do not overlap with this processing unit. Please note that in case skipVoidAreas is activated only the relevant subset of input rasters is available for the formula (cf details for details).

# Formula syntax

The following sections describe in detail how arbitrary formulae can be specified in Module Algebra using either Calculator or Python syntax.

## Calculator syntax

This section describes how to define the algebraic formula using the calculator syntax. Calculators are a specific form of generic filters, which allow calculations, conditional expression, and assignments in addition to pure data sub-selection. Due to the better performance, it is highly recommended to use the calculator syntax whenever possible.

Within Module Algebra, the grid location and the grid values can be used accessed in the formula in the following way:

 Variable name Description x x-coordinate (Easting) of the pixel location y y-coordinate (Northing) of the pixel location r[ds_idx][b_idx|b_name] pixel value with dataset index ds_idx and band index bd_idx (indices starting with 0), and band name b_name invalid general void indicator

The specification of the 2nd (band) index is optional and defaults to band index 0. Thus r[0][0] is equivalent to r[0]. This simplifies processing of single-band grid/raster models (e.g. DTMs). In case of multi-band rasters, the bands can either be accessed via index (starting from 0) or via the respective band name. Named bands are commonplace, e.g., for RGB files, where the bands are typically named "red", "green", and "blue". The possibility to access bands via their names is specifically useful when combining datasets from different sources. E.g., for creating an RGB mosaic from multiple RGB as well as multi-spectral datasets, it is more convenient to use band names as the index identifying a certain color band may differ from dataset to dataset.

Using the syntax above we can now combine one ore more input grid values in an algebraic (calculator) expression to compute a scalar value representing the respective output grid value. The general assignment syntax is:

[band](type)=calculator_expression

band denotes the index or name of the raster band and type its data type (uint8, int8, uint16, ..., float, double). Alpha-numeric band names can be written with or without quotes, i.e., [red] or ["red"]. Quoted numeric values are explicitly used as (string) band names (e.g., ["0815"]) whereas non-quoted numeric values are used as band index (e.g., [0]). Both band names and types are optional. A default data type can be set via the global parameter data_type_grid. 32-bit float ((float or float32, respectively) is used if the data type is neither explicitely specified in the band assignment formula nor set via data_type_grid. Multiple bands can be calculated in a single call by separating the individual assignments with semicolons:

[band_1](type_1)=calculator_expression_1; [band_2](type_2)=calculator_expression_2; ...

Please note that in case band names (and types) are omitted, the bands appear in the output file exactly in the same order as specified in the formula. In order to maintain compatibility with provious versions of Module Algebra, a single non-assignment statement as last statement of the formula string is permitted. Thus, the following statement is valid:

[band_1](type_1)=calculator_expression_1; [band_2](type_2)=calculator_expression_2; expression_3

The above example results in a 3-band output raster although only 2 assignments are present. Assuming a single input dataset with a single band named "grad" (gradient raster), an output slope raster in % can be calculated as follows:

[slope](float)=r[0][0]*100
["slope"]=r[0][0]*100
[0]=r[0][0]*100
r[0][0]*100

The prior three and the latter two statements are equivalent. Whereas the prior three statements use explicit naming of the output band, an unnamed band is produced in the latter two staments. The shortest possible syntax (omitting the band index) is:

r[0]*100

In the above formula only a single input dataset is required, accessed by r[0]. In general, the datasets are strictely accessed by their index starting form 0, whereas the bands can either be accessed by name or index as explained above. The result of the expression (i.e. a scalar value) is assigned to the corresponding output grid cell.

It is recommended to use named bands as this facilitates interpretation of the dataset contents. Band names are consequently reported in all OPALS related log files. However, for the sake of brevity, band names and types are omitted in the following examples and the double indexed multi-band syantax is only used when necessary.

Another simple example is to subtract two input grids:

[]=r[1]-r[0]

In this case, the difference is returned for all those pixels where both r[0] and r[1] exist. If either of the operands is invalid (i.e. nodata) the resulting output pixel will be invalid (nodata) as the expression cannot be evaluated correctly.

A typical application of Module Algebra is the derivation of (binary) grid masks. In this case either 0 or 1 is returned depending on a certain condition. Conditions can be specified in calculator syntax using the ternary operator:

[]=r[0]*r[1]<1.0 ? 0 : 1

In the above example, first the condition (r[0]*r[1]<1.0) is evaluated and depending on the binary result of the conditional expression 0 or 1, respectively, is returned as result of the entire formula. Please note that if one of the operands in the above condition does not exist, the entire condition evaluates to false, i.e. the  else  branch of the ternary operator is returned.

Within conditional expressions it is often desired to check both existence and value of certain pixels. Consider a formula returning the ratio between two grids:

[]=r[0]/r[1]

For this formula we might want to check, if the grid value of the second input grid (r[1]) exists and is not equal to zero (otherwise the division will fail). Generic filters allow to query the existence of grid values by internally comparing the actual grid value with the grids (inherent) nodata value. For the example at hand we can simply write:

[]=r[1] && r[1]!=0 ? r[0]/r[1] : invalid

• existence check: r[1] checks the existence of the grid value; i.e, the grid contains a valid value at the respective location.
• concatenation of logical expressions: the existence (r[1]) and not-equal-zero (r[1]!=0) checks are concatenated with a logical "and" (&&)
• invalid is returned if the (entire) conditional expression is "false", which results in an invalid (nodata) output pixels.

To support the creation of grid/raster mosaics, the generic filters provide statistical operators (min, max, mean, ...) for the entire stack of grid/raster values. Consider a set of stripwise DSMs (i.e. single-band grid models). In order to create a mosaic of all strip DSMs, we can simply write:

[]=mean(r)

This returns the mean of all n valid input grids cells. Invalid grids cells (i.e., cell value = nodata value) are not considered for the statistical calculation. Expanding the previous single-band use case to mosaicing of multi-band RGB orthophoto tiles leads to the following syntax, this time using the full assignment syntax:

[red](uint16)=mean(r[:][red]); [green](uint16)=mean(r[:][green]); [blue](uint16)=mean(r[:][blue])

In this statemens (i) the output bands are explicitely named (as recommended in general) and (ii) the data type is set to 16-bit unsigned integer. The statement r[:][red] means "use the red band of all input datasets". The colon in the first (dataset) index is a special form of array slicing as, e.g. used in MATALB or Python. The general slicing syntax is: from:to, where from and to denote the optional start index (included, default=0) and end index (excluded, default:len+1). For calculating the mean of the red bands of the first four datasets, use either:

[red](uint16)=mean(r[0:4][red])

or:

[red](uint16)=mean(r[:4][red])

By using negative indices (-1=last element, -2=one but the last element, ...), the same can be done for the last four datasets as:

[red](uint16)=mean(r[-4:][red])

Slicing is generally supported for both dataset and band indices. Thus, a simple panchromatic mosaic can be calculated via:

[pan](uint16)=mean(r[:][:])

Slicing is neither meaningful nor supported for band names (e.g. r[:][red:blue]).

In order to increase the readability of the formula, intermediate calculation results can be stored in variables, which can immediatly be used in subsequent epressions and/or assignments. In the following example, a "raster cell density" is calculated by first defining the number of involved rasters, secondly counting the number of valid cells and, finally, returning the cell density.

_n=10;_cnt=count(r);[density]=_cnt/_n

In this example, only the latter statement follows the general rule for band assignments using square brackets. The first two statements, in turn, just store intermediate results, which are not considered for storage in the output raster.

Please note that, if no algebraic expression or band assignment is specified, a logical result is returned rather than the value. Thus, for a single-band input dataset

r[0]

returns true (1) if r0 exists, otherwise false (0). Consequently, if we simply want to return the value of the first input grid, e.g. using Algebra as a raster format converter, we would have to write:

1*r[0]

or

[]=r[0]

Generic filters are a very powerful tool providing a set of arithmetic, logical, conditional operators, and other features. As stated above, a detailed documentation can be found here.

## Python syntax

Python syntax for specifying formulae is deprecated. While it was supported till OPALS version 2.3.2 for single-band rasters, Python syntax support is currently suspended and intended to be re-activated for upcoming versions including support of Python files (*.py) as input.

    Delete from here ****************


An alternative way of specifying formulae in Module Algebra is to use Python syntax. In this case an internal Python interpreter is used to evaluate the formula. Although the performance is generally worse, the Python syntax offers more flexibility and should be used whenever the formula is too complex to be written in calculator syntax. The basic syntax of Python style formulae is:

return ...

where the return value represents the output grid value. Normally, the output grid value is a function of one or more input grid values. Let's assume that our input grid contains surface slope values (tangent of the steepest slope). In order to compute a grid containing the values of steepest slope in % we need to multiply the original slope values by 100. In this case, the formula reads:

return 100 * r[0]

where r[0] represents the grid value of the input grid. Of course, the combination formula can be much more complex, incorporating multiple grids, conditional statements and the like. A simple example based on two input grid models, which are used to create a grid mask (0|1) based on a certain threshold is:

return 0 if r[0]*r[1]<1.0 else 1

The formula must be specified as a string in Python syntax. In fact, the Python interpreter is used to parse and evaluate the condition. Please refer to the Python documentation for a detailed description of the Python syntax. In the following, important information concerning the construction of formulas is summarized:

To access certain information of the considered grid models (pixel values, no data indicators, etc.), the following well-known variables are predefined:

 Variable name Description x x-coordinate (Easting) of the pixel location y y-coordinate (Northing) of the pixel location r[i] pixel value of i-th input grid model (index starting with 0) None general void indicator (c.f. below for details).

None is a constant representing a null-object in Python. It can be used to express empty or invalid objects of any type. In Module Algebra in particular, None is used to handle invalid input or output pixels. E.g., to verify if all input pixels are valid, one can write (for better readability in multi-line style):

if r[0]==None and r[1]==None and ... :
return None
else:
return ....

The very same statement can be written in more compact format, without the None keyword:

if r[0] and r[1] and ... :
return
else:
return ....

Thus, in the first line the statement  r[0]==None is equivalent to not r[0] (i.e. r[0] is invalid), and in the second line it is even allowed to skip the return value, as an empty return is equivalent to returning an empty object.

To combine the specific pixel values, a series of logical and mathematical operations are provided by Python. For a complete reference please refer to the Python language reference. The most important operations are listed below:

Boolean operations (or,and,not) ordered by ascending priority):

 Operation Result x or y if x is false, then y, else x x and y if x is false, then x, else y not x if x is false, then True, else False

Comparison operations:

 Operation Meaning < strictly less than <= less than or equal > strictly greater than >= greater than or equal == equal != not equal is object identity is not negated object identity

Built-in arithmetic operations:

 Operation Result x + y sum of x and y x - y difference of x and y x * y product of x and y x / y quotient of x and y x // y (floored) quotient of x and y x % y remainder of x / y -x x negated +x x unchanged abs(x) absolute value or magnitude of x int(x) x converted to integer long(x) x converted to long integer float(x) x converted to floating point complex(re,im) a complex number with real part re, imaginary part im. im defaults to zero. c.conjugate() conjugate of the complex number c. (Identity on real numbers) divmod(x, y) the pair (x // y, x % y) pow(x, y) x to the power y x ** y x to the power y

Apart from these basic operations, conditional statements are of crucial importance for the construction of complex formulas. Conditions can be specified either as single-line or as multi-line python statements:

Single-line conditional statement:

return ... if condition else ...

Multi-line conditional statement:

if condition1:
return ...
elif condition2:
return ...
else:
return ...

However, since the formula must be passed as a single string argument, line-breaks have to be entered via the special \n character sequence. The correct formula string corresponding to the example above is:

if condition1:\n return ...\nelif condition2:\n return ...\nelse:\n return

Since Python is an indent-oriented language (in opposition to languages using parentheses for grouping), the blank characters before the return statements are mandatory. The same applies to the \nelif... sequence, where no blank is allowed in order to ensure the correct indentation.

The formulas are not restricted to conditional statements, in fact, all control statements (for-loops, while-loops, etc.) are supported. The only preconditions are:

• the formula string must be provided in correct Python syntax
• all control paths of the formula must return a value (if a control path does not return a value, the respective pixel will be void.)
    Delete to here ****************


# Parameter description

-inFileinput grid files
Type: opals::Vector<opals::Path>
Remarks: mandatory
The grid model files are expected in GDAL supported format. The grid values of the respective models are referred to as r[0], r[1] ... in the calculation formula.
-outFileoutput grid file name
Type: opals::Path
Remarks: estimable
Path of output grid file in GDAL supported format.
Estimation rule: The current path and the name (body) of the input file are used as file name basis. Additionally, the default postfix '_algebra' and the extension corresponding to the output format are appended
-oFormatoutput file format [tif,asc,dem,dtm ...]
Type: opals::String
Remarks: estimable
Grid file format, use GDAL driver names like GTiff, AAIGrid, USGSDEM, SCOP ... .
Estimation rule: The output format is estimated based on the extension of the output file. If neither outFile nor oFormat are specified, GeoTiff is used as default.
-formulaalgebraic formula for combining input grids
Type: opals::CalculatorWithPythonSupport<DM::ICalculator::ReadAccess::coordinates | attributes | rasters, DM::ICalculator::WriteAccess::coordinates | attributes | rasters>
Remarks: mandatory
Formula to compute a single grid value as a string in in Generic Filter or Python syntax.
A formula must contain return a value (e.g. generic filter: r[0]-r[r1] or python: return r[0]-r[1])). For more details please refer to the module documentation.
CalculatorWithPythonSupport<coordinates | attributes | rasters, coordinates | attributes | rasters>: has read access to: coordinates | attributes | rasters; has write access to: coordinates | attributes | rasters; Prefix with 'python:' to indicate Python syntax.
See Calculator
-gridSizegrid width x/y
Type: opals::Vector<double>
Remarks: estimable
Size of output grid cell in x- and y-direction. Estimation rule: output grid size = grid size of 1st inFile
-noDatanodata value indicating void pixels
Type: opals::NoDataType
Remarks: default=max
Value representing an undefined value in the output grid model.
NoDataType: Defines or disables (use 'none') the nodata value for the output raster. The nodata values can be a specific value or a type independent label, such as, 'min', 'max' and 'nan'. Those labels will be translate to values based on the raster band type. See the OPALS Docu (opals::NoDataType) for more Details.
Syntax: none | min | max | nan | <value>
-limit2D clipping window
Type: opals::MultiGridLimit
Remarks: optional
If -limit is skipped, the common area (intersection) of all input datasets is used.
MultiGridLimit: Defines xy-limits for output grid/raster datasets, either explicitly or depending on the limits of input datasets. The limit coordinates explicitly specified by left/lower/right/upper either refer to pixel centers (default) or corners, and may optionally get rounded to multiples of the grid size. 'intersect' and 'union' choose the intersection/union of limits of all inputs. dataset:i chooses the extents of the i'th dataset, with i in [0;nInputs).
Syntax: ['center'|'corner']['round']['intersect'|'union'|'dataset:'i|'('left lower right upper')']
-resamplinggrid resampling method.
Type: opals::ResamplingMethod
Remarks: default=bilinear
Possible values:
nearestNeighbour ... nearest neighbour
bilinear ........... bi-linear interpolation
bicubic ............ bi-cubic interpolation
never .............. indicator avoiding resampling

Method to be used for resampling the input grids to the specified grid structure of the output grid.
-maxMemorymaximum RAM memory [MB]
Type: uint32
Remarks: default=500
The memory consumption scales with the number of input grids in order to provide good reading performance. An upper memory limit can be specified in order to avoid overflows in case of numerous input grids (e.g. large mosaic),
-skipVoidAreasskip void input datatset parts
Type: bool
Remarks: mandatory
Estimation rule: In opalsAlgebra, the output raster is processed block-wise. In general, for eachblock the pixels of all input datasets are queried, even if they only contain void pixels. The latter is typically the case in merging/mosaicing applications where the input datasets are either disjunct or only partially overlapping. To improve the run-time performance in such cases activating the skipVoidAreas flag allows ignoring blocks of all-void input pixels. Please note, that in this case the order of the input datasets is not strictly maintained which might lead to invalid results if the calculation formula explicitely uses input dataset indices (e.g. -formula r[0]?r[1]:r[2]). However, for all formulae which aggregate only the valid pixels (e.g. -formula mean(r)), activating skipVoidAreas substantially speeds up processing.

# Examples

The data used in the subsequent examples are located in the \$OPALS_ROOT/demo/ directory and the following preprocessing steps (data import and surface grid calculation) are required:

opalsImport -inFile strip31.laz
opalsImport -inFile strip21.laz
opalsImport -inFile strip11.laz
opalsImport -inFile fullwave2.fwf -iFormat fwf
opalsGrid -inFile strip31.odm -outFile strip31.tif -interpolation movingPlanes
opalsGrid -inFile strip21.odm -outFile strip21.tif -interpolation movingPlanes -feature sigmaZ excen
opalsGrid -inFile strip11.odm -outFile strip11.tif -interpolation movingPlanes
opalsGrid -inFile fullwave2.odm -outFile dsm_mp.tif -interpolation movingPlanes -feature sigma0
opalsCell -inFile fullwave2.odm -outFIle dsm_max.tif -cellSize 1 -feature max -limit center

Example 1: Roof model

In the first example, the roofs of a digital surface model are extracted by applying a simple height threshold. All other values are set to void (no data).

Generic filter syntax:

opalsAlgebra -inFile strip21.tif -outFile strip21_roofs.tif -formula "r[0]>275. ? r[0] : invalid"

Python syntax:

opalsAlgebra -inFile strip21.tif -outFile strip21_roofs_py.tif -formula "python: return r[0][0] if r[0][0] is not None and r[0][0]>275. else None"

The output is written to a separate grid file in GeoTiff format (strip21_roofs[_py].tif) and a standard color coded visualization of it is shown in the following figure:

Figure 1: roof model derived by Module Algebra

In this example, two grid models describing roughness (sigmaZ) and shadowing effects (excen), both by-products of the moving planes surface interpolation with /ref ModuleGrid, are used to create a binary grid mask (0..invalid,1..valid pixel).

Generic filter syntax:

opalsAlgebra -inFile strip21_sigmaZ.tif -inFile strip21_excen.tif -outFile strip21_mask.tif -formula " r[0]<0.1 && r[1]<0.8 ? 1 : 0"

Python syntax:

opalsAlgebra -inFile strip21_sigmaZ.tif -inFile strip21_excen.tif -outFile strip21_mask_py.tif -formula "python: return 1 if r[0][0] is not None and r[0][0]<0.1 and r[1][0]<0.8 else 0"

The very same condition can be written more compact as a single (binary) return statement:

Generic filter syntax:

-formula "r[0]<0.1 && r[1]<0.8"

Python syntax:

-formula "python: return r[0][0] is not None and r[0][0]<0.1 and r[1][0]<0.8"

Example 3: DTM Mosaic I

A simple way of creating a DTM mosaic based on multiple input grids (e.g. each grid covering a single map sheet) is by selecting the first valid grid value from the list of input grids for each grid position using Python's for .. in .. : statement:

Python syntax:

opalsAlgebra -inFile strip31.tif strip21.tif strip11.tif -limit union -outFile mosaic_py.tif -formula "python:for x in r:\n if x[0] is not None and x[0]:\n return x[0]"

This commands takes three input grids (strip31.tif strip21.tif strip11.tif), calculates the xy-extents of the output grid as the union of all input datasets (-limit union) and, for each grid post, returns the first grid point with a valid height (-formula "... x:\n return x") within a loop over all input grids (-formula "for x in r[:]: ...").

To calculate a DTM mosaic based on the mean value of all valid pixels the following command in generic filter sysntax can be used:

opalsAlgebra -inFile strip31.tif strip21.tif strip11.tif -limit union -outFile mosaic.tif -formula "mean(r)"

The same result can be achieved using the wild card syntax for specifying the input files:

opalsAlgebra -inFile strip?1.tif -limit union -outFile mosaic_wc.tif -formula "mean(r)"

Example 4: DTM Mosaic II

In the following example a (global) point density model is created based on strip-wise point density models. The strip-wise point density models can be created with the following command:

opalsCell -inFile strip31.odm -feature pdens -cellSize 5
opalsCell -inFile strip21.odm -feature pdens -cellSize 5
opalsCell -inFile strip11.odm -feature pdens -cellSize 5

Now, we can calculate the global point density model as:

Generic filter syntax:

opalsAlgebra -inFile strip31_z_pdens.tif strip21_z_pdens.tif strip11_z_pdens.tif -outFile pdens.tif -limit union -formula "count(r)>0 ? sum(r) : invalid"

Python syntax:

opalsAlgebra -inFile strip31_z_pdens.tif strip21_z_pdens.tif strip11_z_pdens.tif -outFile pdens_py.tif -limit union -formula "python:sum=0\nfor x in r:\n if x[0] and x[0] is not None:\n sum+=x[0]\nreturn sum if sum>0 else None"

Contrary to the previous example, all valid point densities values are summed up and the sum is returned (if the sum is greater 0).

Example 5: Landscape dependent DSM

Finally, opalsAlgebra is used for the derivation of a landscape dependent surface model as described in Hollaus et al (2010). Hereby, the single DSM pixels are either selected from a raster map containing the highest elevation within a cell (dsm_max.tif) or from a grid model calculated via moving planes interpolation (dsm_mp.tif) depending on a roughness map (dsm_mp_sigma0.tif, by-product of moving planes interpolation).

Generic filter syntax:

opalsAlgebra -inFile dsm_max.tif dsm_mp.tif dsm_mp_sigma0.tif -outFile DSM.tif -formula " r[2] < 0.2 || !r[0] ? r[1] : r[0]"

Python syntax:

opalsAlgebra -inFile dsm_max.tif dsm_mp.tif dsm_mp_sigma0.tif -outFile DSM_py.tif -formula "python:return r[1][0] if (r[2][0] is not None and r[2][0] < 0.2) or not r[0][0] or r[0][0] is None else r[0][0]"

Figure 2 shows hill shadings of all three DSM variants derived with the following commands:

Figure 2: Hill shadings; DSM max (left), DSM moving planes (center), landspace dependent DSM (right)

Example 6: Split RGB orthophoto in 4 parts (tiles)

In this example an orthophoto with RGB bands is splitted to sub images using limit parameter.

Generic filter syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_1.tif -limit "(600865 5340530 600957 5340575)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "[red](uint8)=r[0][red];[green](uint8)=r[0][green];[blue](uint8)=r[0][blue]" -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_2.tif -limit "(600865 5340575 600957 5340620)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "[red](uint8)=r[0][red];[green](uint8)=r[0][green];[blue](uint8)=r[0][blue]" -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_3.tif -limit "(600865 5340620 600957 5340665)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "[red](uint8)=r[0][red];[green](uint8)=r[0][green];[blue](uint8)=r[0][blue]" -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_4.tif -limit "(600865 5340665 600957 5340710)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "[red](uint8)=r[0][red];[green](uint8)=r[0][green];[blue](uint8)=r[0][blue]" -noData none

Python syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_1_py.tif -limit "(600865 5340530 600957 5340575)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "python: return (('red',r[0]['red']), ('green',r[0]['green']), ('blue',r[0]['blue']))" -data_type_grid uint8 -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_2_py.tif -limit "(600865 5340575 600957 5340620)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "python: return (('red',r[0]['red']), ('green',r[0]['green']), ('blue',r[0]['blue']))" -data_type_grid uint8 -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_3_py.tif -limit "(600865 5340620 600957 5340665)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "python: return (('red',r[0]['red']), ('green',r[0]['green']), ('blue',r[0]['blue']))" -data_type_grid uint8 -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_rgb_4_py.tif -limit "(600865 5340665 600957 5340710)" -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "python: return (('red',r[0]['red']), ('green',r[0]['green']), ('blue',r[0]['blue']))" -data_type_grid uint8 -noData none

Each sub-raster is stored with band names and data types which are defined in formula.

Example 7: Split RGB orthophoto in red/green/blue color channels

In the following example each band of an RGB orthophoto is seperately saved as single-band rasters.

Generic filter syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_r.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "[red] (uint8)=r[0][0]" -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_g.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "[green](uint8)=r[0][1]" -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_b.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "[blue] (uint8)=r[0][2]" -noData none

Python syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_r_py.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "python:return ('red',r[0][0])" -data_type_grid uint8 -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_g_py.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "python:return ('green',r[0][1])" -data_type_grid uint8 -noData none
opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_b_py.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=MINISBLACK -formula "python:return ('blue',r[0][2])" -data_type_grid uint8 -noData none

Example 8: Mosaicing

Following command generates an RGB mosaic image from 4 raster files.

Generic filter syntax:

opalsAlgebra -inFile vie_cityHall_rgb_?.tif -outFile vie_cityHall_rgb_mosaic.tif -formula "[red](uint8)=mean(r[:][red]); [green](uint8)=mean(r[:][green]); [blue](uint8)=mean(r[:][blue])" -limit union -noData none

Python syntax:

opalsAlgebra -inFile vie_cityHall_rgb_?.tif -outFile vie_cityHall_rgb_mosaic_py.tif -formula "python:return (('red',mean(r[:,'red'])), ('green',mean(r[:,'green'])), ('blue',mean(r[:,'blue'])))" -data_type_grid uint8 -limit union -noData none

Example 9: Merge colour channels to single mult/band RGBI raster

In this example selected bands from input rasters are organized and stored to a multiband single raster image.

Generic filter syntax:

opalsAlgebra -inFile vie_cityHall_r.tif vie_cityHall_g.tif vie_cityHall_b.tif vie_cityHall_nir.tif -outFile vie_cityHall_rgbi.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "[red](uint8)=r[0];[green](uint8)=r[1];[blue](uint8)=r[2];[nir](uint8)=r[3]" -noData none

Python syntax:

opalsAlgebra -inFile vie_cityHall_r.tif vie_cityHall_g.tif vie_cityHall_b.tif vie_cityHall_nir.tif -outFile vie_cityHall_rgbi_py.tif -create_option COMPRESS=JPEG JPEG_QUALITY=90 PHOTOMETRIC=RGB -formula "python:return (('red',r[0][0]), ('green',r[1][0]), ('blue',r[2][0]), ('nir',r[3][0]))" -data_type_grid uint8 -noData none

Example 10: Calculate indices

The following commands calculate VARI(Visible Atmospherically Resistant Index) and NDVI(Normalized Difference Vegetation Index) indices of input rasters. Then calculated indices are stored as a band in output raster.

Generic filter syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_vari.tif -formula "[vari]=(r[0][green]-r[0][red])/(r[0][green]+r[0][red]-r[0][blue])"
opalsZColor -inFile vie_cityHall_vari.tif -palFile greyPal.xml -scalePal -5,5 -interval 0.2
opalsAlgebra -inFile vie_cityHall_rgbi.tif -outFile vie_cityHall_ndvi.tif -formula "[ndvi](float)=(r[0][nir]-r[0][red])/(r[0][nir]+r[0][red])"
opalsZColor -inFile vie_cityHall_ndvi.tif -palFile greyPal.xml -scalePal -1,1 -interval 0.1

Python syntax:

opalsAlgebra -inFile vie_cityHall_rgb.tif -outFile vie_cityHall_vari_py.tif -formula "python:return ('vari',(r[0]['green']-r[0]['red'])/(r[0]['green']+r[0]['red']-r[0]['blue']))"
opalsAlgebra -inFile vie_cityHall_rgbi.tif -outFile vie_cityHall_ndvi_py.tif -formula "python:return ('ndvi',(r[0]['nir']-r[0]['red'])/(r[0]['nir']+r[0]['red']))" -data_type_grid float32

# References

M. Hollaus, G. Mandlburger, N. Pfeifer, W. Mücke: Land cover dependent derivation of digital surface models from Airborne laser Scanning data; in: "IAPRS Volume XXXVIII Part 3A", (2010), ISSN: 168299750; 221 - 226.

Date
26.05.2016
opalsZColor is the executable file of Module ZColor
Definition: ModuleExecutables.hpp:248
@ pdens
point density estimate
Definition: GridFeature.hpp:13
@ uint16
16 bit unsigned integer
@ movingPlanes
moving (tilted) plane interpolation
@ sigma0
sigma 0 of grid post adjustment (i.e. std.dev. of the unit weight observation)
Definition: GridFeature.hpp:12
@ band
raster band number/index in case of multi-layer rasters
opalsAlgebra is the executable file of Module Algebra
Definition: ModuleExecutables.hpp:13
opalsImport is the executable file of Module Import
Definition: ModuleExecutables.hpp:113
@ strip
strip ID of an image (opalsStripAdjust)
@ slope
steepest slope [%]
Definition: GridFeature.hpp:16
@ scalePal
scale factor to be applied to the given palette (opalsZcolor)
@ mean
Mean.
This is the fake namespace of all opals Python scripts.
Definition: doc/temp_doxy/python/__init__.py:1
opalsGrid is the executable file of Module Grid
Definition: ModuleExecutables.hpp:93
opalsCell is the executable file of Module Cell
Definition: ModuleExecutables.hpp:23
@ uint8
8 bit unsigned integer
@ count
always last element
@ float32
32 bit floating point
@ condition
@ feature
Use a statistic feature of the boundary gap points for filling.
@ none
Suppress all logging output.
@ formula
formula string for albegraic grid computations (opalsAlgebra)
@ data_type_grid
default output grid/raster data type
@ create_option
dataset create options
@ palFile
palette file (opalsZcolor)