Table of Contents
- See Also
- opals::IGrid
Aim of module
Derives a 2.5D digital elevation models (DEM) in regular grid structure using simple interpolation techniques.
General description
Digital Surface Models (DSM) and Digital Terrain Models (DTM) are among the most important products derived from ALS point cloud data. DTM production based on ALS data is a complex task comprising (semi-)automatic modelling of structure lines, classifying the point cloud in terrain and off-terrain points and, finally, interpolating a smooth surface using high-level interpolation techniques like linear prediction or Kriging, respectively. In contrast, there are also a series of tasks within the ALS processing chain, which do not require the utmost quality concerning the geomorphological surface quality but, instead, shortest possible derivation times are crucial. This becomes especially true for the quality check of overlapping ALS strips (relative strip differences) based on distinct strip DSMs.
For this reason, Module Grid aims at deriving surface models in regular grid structure with maximum speed. To achieve high computation performance:
- the OPALS Data Manager (ODM) providing high performance nearest neighbour queries within the ALS point cloud is used as data basis
- multi-core CPUs of modern computers are used extensively
- only simple interpolation methods are provided
- only rectangular grids (rasters) are derived, i.e. no break lines are considered
The ODM serves as data basis (parameter inFile) for the grid interpolation. Thus, the import of ALS point cloud data, stored in arbitrary data formats, is a prerequisite for Module Grid. For details concerning data import, please refer to Module Import. The interpolation results are stored as a regular grid file (parameter outFile) in a GDAL supported data format (parameter outformat). The internal processing is organised in patches (tiles), thus, maximum performance can be achieved for data formats supporting tile wise storage like, GeoTiff or SCOP.RDH. Sequential data formats (e.g. ArcInfo ASCII grid, DTED ...) can else be written directly by Module Grid, but probably with lower performance. The derived grid model is strictly stored in "pixel is point" interpretation, i.e. the grid values represent the interpolated heights at the pixel center. However, user defined limits (parameter limits) for the output grid model can either refer to the centers (=default) or corners of the outermost pixels (c.f. parameter limit description). For a more detailed discussion "grid versus raster" please refer to the glossary.
The following interpolation strategies (parameter interpolation) ordered by increasing processing time are offered:
- Snap grid: This is the simplest method: For each point the appropriate grid cell is determined and the height of the point is mapped to this cell. Please note that subsequent ALS points belonging to the same grid cell overwrite previous values. The snap grid parameter is mainly useful in case of already gridded input points in order to avoid additional interpolation.
- Nearest neighbour: For each grid cell the nearest data point is queried from the ODM and the height of this nearest neighbour point is mapped to the grid cell. The quality of the derived surface is still poor, but computation performance is very high.
- Delaunay triangulation: The input points are triangulated and the grid heights are interpolated from the TIN. Please note, that the processing is tile based; for each tile a local TIN is constructed. To avoid data gaps additional points from neighbouring tiles are considered within a distance controlled by the parameter searchRadius.
- Moving average: For each grid cell n nearest points (parameter neighbours) are queried and the average height of all n points is mapped to the grid cell.
- Moving planes: For each grid cell n nearest points (parameter neighbours) are queried and a best fitting tilted plane (minimizing the vertical distances) is estimated. The height of the resulting plane at the grid point (x,y) position is mapped to the grid cell. The tilted plane interpolator allows the derivation of slope measures (nx, ny, slope, exposition) for each grid point.
- Robust moving planes: As above, with the difference that outliers are detected during the plane fitting using a robust estimation approach. Point classified as outliers do not take part of the plane fitting. To achieve reasonable results a high number of neighbouring points should be considered and an appropriate weight function should be provided
- Moving paraboloid: For each grid cell n nearest ALS points (parameter neighbours) are queried and a best fitting paraboloid (2nd order polynomial) is estimated. The height of the resulting paraboloid at the grid point (x,y,) position is mapped to the grid cell. The paraboloid interpolator allows the derivation of curvature measures (kmin,kmax,kmean,kgauss) for each grid point.
The last four methods (moving average, moving planes, moving paraboloid, robust moving planes) require the specification of the number of neighbouring points considered for the interpolation of a single grid post (parameter neighbours). The results of the neighbour queries can be restricted to a maximum search radius around the grid point (parameter searchRadius), thus, enabling to reflect areas with sparse point density in the resulting grid as void pixels.
Data selection strategies
Furthermore, two different data selection strategies are offered: parameter selMode=nearest is the default mode representing the straight forward nearest neighbour query, whereas parameter selMode=quadrant yields a quadrant oriented data selection. In the latter case, the returned points are equally distributed among the four quadrants (NW, NE, SE, SW) around the (x,y) interpolation location (c.f. figure below), as long as, the number of neighbours is a multiple of four (4, 8, 12, ...). Hence, the first, second, third, etc. generation of quadrant neighbours are returned. If the number of neighbours has a remainder to 4, these "extra points" are also selected from different quadrants based on the distance criterion. In cases with empty (or inappropriate covered) quadrants this selection method only returns fully covered generation of quadrant neighbours, plus the appropriate number of extra points. E.g. if there are no points in two quadrants this method will always return 2 points (from the covered quadrants) even with higher number of neighbours than two. Please note, that quadrant oriented data selection avoids grid extrapolation effects but slows down the performance. Alternatively, grid extrapolation effects can also be reduced by specifying a greater number of neighbouring points together with using the default selection mode (parameter selMode nearest), which often results in better computation performance.
As mentioned above, the primary aim of Module Grid is the derivation of surface grid models. In addition, additional feature models (parameter feature) can be derived simultaneously as some of them are side products of the grid interpolation and, thus, can be provided without loss of performance. The following additional features are supported:
- sigmaz: std.dev. of interpolated grid height (moving average/planes only)
- sigma0: std.dev. of the unit weight observation (moving average/planes only)
- density: point density estimate (moving average/planes only)
- excentricity: distance grid point - center of gravity of data points
- slope: steepest slope in % (moving planes only)
- slpDeg: steepest slope in degree (moving planes only)
- slpRad: steepest slope in radians (moving planes only)
- slope: steepest slope in % (moving planes only)
- exposition: slope aspect [rad] = azimuth of steepest slope line, N=0, clockwise sense of rotation (moving planes only)
- normalx: x-component of the surface normal unit vector (moving average/planes only)
- normaly: y-component of the surface normal unit vector (moving average/planes only)
- kmin: minimum curvature (moving paraboloid only)
- kmax: minimum curvature (moving paraboloid only)
- kmean: mean curvature H: H = (kmin+kmax)/2 (moving paraboloid only)
- kgauss: gaussian curvature K: kmin*kmax (moving paraboloid only)
By default, the additional features are stored as separate grid model files in identical grid structure and data format as the corresponding surface model. The file names are either constructed automatically, if only the name of the surface grid model file is defined, based on the definition of the respective global postfix parameters (postfix_sigma, postfix_pdens, postfix_excen, postfix_slope, postfix_slpDeg, postfix_slpRad, postfix_expos, postfix_nx, postfix_ny, postfix_kmin, postfix_kmax, postfix_kmean, postfix_kgauss,), or individual output file names can be specified. In the latter case, the number of entered file names must correspond to the number of derived grid models (surface model + feature models). Additional feature models are especially useful for subsequent creation of grid masks using Module Algebra. If the boolean parameter multiband is set to true
or 1
, respectively, only a single raster file with multiple bands is created. The first band always contains the z/attribute-raster followed by the feature bands in the order specified by the user.
In general, each point of the point cloud has the same influence on the run of surface. However, sometimes it might be useful to give preference to certain points (e.g. points with higher amplitude/reflectance or points closer to the grid point). Individual point weights can be achieved by specifying a point weight function formula (weightFunc) as a string in generic filter syntax. Individual point weights are reasonable/supported for all k-nn based interpolators (e.g.: movingAverage, movingPlanes, movingParaboloid, robMovingmovingPlanes). The formula is evaluated for each of the k-nn points and must return a value which is used as weight for the respective point observation. To define a point weight, the following elements can be used:
- All attributes of the k-nn point (e.g., Amplitude, EchoWidth, etc.)
- The 2D-distance between the grid point and the knn-point. Within the formula, the grid point can be accessed as n[0] and the k-nn point as n[1]. To calculate the (squared) distance, the built-in functions SqrDist2D(n[0],n[1]) and Dist2D(n[0],n[1]). Please note, that only 2D distances can be used as the grid points do not contain height information (before the interpolation).
As a general rule, the formulas should be defined so that valid weights can be computed for each k-nn point. A grid point will automatically be set to void (no data) if the weight function cannot be calculated for a single point of the k-nn set. The most frequent causes for numeric problems are divisions by zero (e.g. distance in the denominator) which can often be avoided by using a proper epsilon. For Inverse Distance Weighing, for instance, predefined formulae are available (idw1, idw2, idw:n, c.f. parameter description below for details).
Please note, that the outlier classification procedure within the moving planes estimation relies on the individual point weights and, thus, the rejection level can be controlled by assigning global or individual point weights. It is good practice to use the height variance (varZ=sqr(sigmaZ)) of a point as the basis for the weight calculation, with w = 1 / varZ.
Parameter description
Remarks: mandatory
The path(s) to the opals datamanager(s) whose point data are used to derive the surface grid model. In case multiple ODMs are specified, the data are queried from all (potentially overlapping) ODMs.
Remarks: estimable
Path of output grid file in GDAL supported format.
Estimation rule: The current directory and the name (body) of the input file are used as file name basis. Additionally, the default extension of the output format is appended.
Remarks: estimable
Use GDAL driver names like GTiff, AAIGrid, USGSDEM, SCOP... .
Estimation rule: The output format is estimated based on the extension of the output file (*.tif->GTiff, *.dem->USGSDEM, *.dtm->SCOP...).
Remarks: default=64
The internal processing is organized in rectangular tiles. The size of these computing units(number of grid lines per tile) can be specified allowing control of the memory consumption of the module. Recommended values are in the range of 64 (default) if the average point distance approximately matches the grid size. In case coarser of finer grids are to be derived, the tile size may need to be adapted. Finer grids allow for larger tile sizes (and vice versa).
Remarks: default=1
Size of the rectangular grid cell in x- and y-direction. If only one value is passed, quadratic cells are used.
Remarks: default=snap
Possible values: snap .................... simple 'snap-height-of-nearest-data-point' interpolation nearestNeighbour ........ nearest neighbour interpolation delaunayTriangulation ... Delaunay triangulation interpolation movingAverage ........... moving average interpolation (=moving horizontal plane) movingPlanes ............ moving (tilted) plane interpolation robMovingPlanes ......... robust moving (tilted) plane interpolation movingParaboloid ........ moving circular, elliptic or hyperbolic paraboloid kriging ................. kriging interpolation
The grid interpolation methods are ordered by quality. Higher-quality surfaces can be achieved using the latter methods but, in turn, the processing performance decreases.
Remarks: default=8
This parameter is only considered for moving average, planes and parabola interpolation.
Remarks: estimable
Only points within the given search radius are considered for the interpolation of a single grid post. If the search area contains too few points for successful interpolation, the respective grid point is omitted (void).Estimation rule: searchRadius = 3 * gridSize .
Remarks: default=nearest
Possible values: nearest .... pure nearest neighbour (nn) selection quadrant ... quadrant-wise nn selection, ie. nn per quadrant, then 2nd nn per quadrant, ...
Please note that quadrant wise data selection on the one hand reduces extrapolation effects, but on the other hand decreases the processing performance
Remarks: optional
An algebraic formula in generic filter syntax can be specified to calculate an individual weight for each point. Both, the points' attributes and the distance w.r.t the grid location can be used to determine the point weight. Inverse Distance Weighting (IDW) can be realized by using the built-in "Dist2D(n[0],n[1])" or "SqrDist2D(n[0],n[1])" function returning the linear or squared distance between the current grid location n[0] and the considered neighbour point n[1]. The following predefined formulae can be used :
IDW1.... inverse linear distance (i.e., 1/(1+ad))
IDW2.... inverse squared distance (i.e., 1/(1+(ad)^2)
IDW:i... inverse distance to the power of i ( i.e., 1/1+(ad)^i)
where 'd' denotes the 2D distance between the grid point and the neighbouring data point and 'a' the shape factor of the weight function. 'a' determines how fast the weight function drops to zero with increasing distance from the grid point. It is determined automatically in a way that points located at the border of the search circle (c.f. -searchRadius) obtain a small weight w=0.05.
Remarks: optional
Possible values: sigmaz ......... sigma z of grid post adjustment (i.e. std.dev. of the interpolated height) sigma0 ......... sigma 0 of grid post adjustment (i.e. std.dev. of the unit weight observation) pdens .......... point density estimate pcount ......... number of points used for interpolation of grid point excentricity ... distance between grid points - c.o.g. of data points slope .......... steepest slope [%] slpDeg ......... steepest slope [deg] slpRad ......... steepest slope [rad] exposition ..... slope aspect [rad] (azimuth of steepest slope line, N=0, clockwise) normalx ........ x-component of surface normal unit vector normaly ........ y-component of surface normal unit vector kmin ........... minimum curvature kmax ........... maximum curvature kmean .......... mean curvature: kmean=(kmin+kmax)/2 kgauss ......... gaussian curvature: kgauss = kmin*kmax kminDir ........ Azimuth [rad] of minimum curvature kmaxDir ........ Azimuth [rad] of maximum curvature absKmaxDir ..... Azimuth [rad] of maximum absolute curvature precision ...... grid post precision (derived from sigmaz, pdens and slope)
If one or more features are specified, additional grid models in the same format and structure as the basic surface model are derived. Please note that if more than one output file is specified, the number of output file names must match the number of specified features: nr. outFile = nr. feat. + 1
Remarks: optional
A filter string in EBNF syntax as described in section 'Filter syntax' can be passed to restrict the set of input points (e.g. to consider last echoes only).
Remarks: optional
If no user defined limits are specified or -limit is even skipped, the OPALS Data Manager (ODM) extents are used. Rounding is enabled by default in case the limits are derived from the ODM and can be enforced for user defined limits (keyword: round).
Remarks: default=0
When using multiband grid features will be stored in multiple bands instead of multiple files.Features will be saved in the same order as user input.For example if -feature sigma0 slope then Band1=sigma0 and Band2=slope
Remarks: default=1
Unfavorable point distributions can lead to extrapolations, where the interpolated value clearly exceeds the local data values. Using the local data range extended by a buffer factory, such extrapolation are detected and rejected if this option is enabled.
Short examples
Generate a DEM from a LAS file
Examples
The data used in the following examples can be found in the $OPALS_ROOT/demo/
directory. As the derivation of (surface) grid model relies on the ODM, the respective point cloud data has to be imported. The datasets fullwave.fwf
and fullwave2.fwf
containing full waveform ALS point data (x, y, z, gps time, amplitude, echo width, echo nr, echo qualifier) of a wooded area are used. The following pre-processing steps are required:
Example 1:
In the first example only the mandatory parameter (inFile) is specified. Thus, the command is very short and simple:
It produces the grid file fullwave_z.tif
in the default GeoTiff file format. The grid extensions are determined from the bounding box of the entire ODM point data set. The simple snap grid interpolation is used, thus a poor quality surface is derived. The grid size was not specified and, therefore, set to 1m.
Example 2:
In the second example some of the optional parameters are defined.
opalsGrid is called twice in this example. Two 0.5m grids are derived using moving planes interpolation. Whereas in the first case, the output is written in GeoTiff file format to file all.tif
, an ArcInfo ASCII grid file (all.asc
) is created in the case. In both examples the output format is estimated based on output file extension (*.tif->GeoTiff, *.asc->AAIGrid). The ASCII grid file contains a 6 line header containing the following georeferencing information of the grid file.
We see, that a grid with 361 columns and 201 lines was created. Please note that the coordinates of the lower left grid corner (XLLCORNER/YLLCORNER) do not coincide with the 0.5m grid lines, but are shifted about half a pixel. This is an indication for the 'pixel-is-point' interpretation. Since XLLCORNER/YLLCORNER denote the coordinates of the lower left pixel corner, the centers of the pixel coincide exactly with the 0.5m grid lines, for which locations the grid values have been derived.
The following figure shows the 'raw' grid file all.tif
(as displayed by IrfanView) and a color coded variant using the default color palette of Module ZColor :
Example 3:
In the following example a filter is utilised to extract the first and last echoes from the ODM separately:
This produces the following output (again color coded using the standard color palette):
It can clearly be seen, that the right grid computed from the last echoes is much smoother than the one derived from the first echoes, as the first echoes contain all the reflections from the vegetation.
Example 4:
As before, this example uses a filter to restrict the point set considered for model interpolation. In addition to the last echo filter, a generic attribute filter is used to select all points with an echo width less than 4.5ns.
Please note, that quoting the filter string is strictly necessary in this case due to the embedded blanks and special characters. It is recommended always to quote complex filter strings. For more detailed information about pre-selecting data please refer to the Filter manual.
Example 5:
Based on Example 4, this example shows how to use individual point weights based on the echo intensity (approximated as the product of amplitude and echo width). The higher the intensity the more influence the point should have on the course of the surface.
Besides attributes, the 2D distances from the k-nn points to the grid point can be used to calculate the point weigth. This is commonly referred to as Inverse Distance Weighting (IDW) and the reciprocal of the squared distances (1/d^2) is often used as point weight. Based on the build in SqrDist2D fcuntion this can be formulated as:
Please note, that the grid point is accessed as n[0] and the k-nn point as n[1]. The above formula results in a divisions by zero if a data point is planimetrically coinciding with a data point. In fact, the fullwave.fwf
dataset contains one last echo at round meter coordinates.
Consequently, the weight function gets infinite when calculated for the grid locations...
... resulting in the following warning message:
This could be avoided by using an appropriate small epsilon value to secure that the denominator is grater than zero. In the special case of Inverse Distance Weighting predefined formulae are available. The proper replacement for the above weight function based on squared distances is IDW2.
Example 6:
Example 6 demonstrates the use of the robust moving plane interpolator, Again, weight functions are used in this example to control the outlier detection. In contrast to example 5, were individual point weights were calculated provided based on the intensity or the location of a point w.r.t the grid point, a mean estimated height accuracy is applied to all points. However, the robust plane estimation is executed twice with different point weights (-weightFunc 1 or -weightFung 25) to demonstrate the effect of more conservative or more progressive outlier removal on the resulting grid (c.f. Fig. 4)
Setting the point weight value to 1 or 25, respectively, corresponds to an estimated a-priori sigmaZ of 1m/0.2m (weight=1/sigmaZ**2). In the latter case (w=25), more points are detected as outliers, as the expected height precision is smaller (Fig. 4). The respective visualizations can be produced with the following commands:
Example 7:
Example 7 shows the usage of multiband. Multiband is a boolean parameter. When set to true (1), each feature is stored in the corresponding band of a single raster file.
The example creates a single output file fwf_multiband.tif
containing four bands. The first band is the z
raster, and the second band includes the sigma0
feature, the third one slope
and the last one the exposition
feature.
- Date
- 23.06.2016