Module DTM
opals::IDTM

# Aim of module

Calculates a high quality Digital Terrain Model in hybrid grid structure (i.e. regular grid + structure lines).

# General description

Airborne Laser Scanning (ALS) provides a dense set of points located at arbitrary positions. The point densities are typically in the range of 1-20 points per square meter. ALS is also capable of penetrating vegetation through small gaps in the foliage and, thus, can provide the ground beneeth the canopy. One of the main applications of ALS is therefore the derivation of Digital Terrain Models (DTM). Interpolating the terrain surface is, however, the final step in a complex processing chain starting with waveform analysis, 3D point cloud derivation (Package opalsPreprocess), quality assessment (Package opalsQuality), geometric and radiometric calibration (Package opalsGeoref), derivation of 3D structure lines (Package opalsGeomorph), and filtering ground and off-terrain points (Package opalsClassify). In general, geo-statistical interpolation strategies have proven to deliver the best results. Kriging or linear prediction, respectively, are so-called BLUE (best linear unbiased estimators), but they are computationally expensive. Dense ALS point clouds therefore justify the usage of simpler interpolation techniques like moving least squares. Module DTM, in particular, uses linear or quadratic functional models (i.e. tilted planes or paraboloids). For each grid point the neigbouring laser points are used to locally fit such basic models and the grid height is calculated from the resulting geometry primitive. In addition to the point cloud, explicit line information e.g. from semi-automatic 3D structure line modelling (cf. Package opalsGeomorph) is considered within the DTM interpolation. Please note, that whereas the final goal of Module DTM is the construction of a hybrid DTM structure consisting of a regular elevation grid (i.e. height field) and vectors (interlinked structure lines), up to now the line information is specifically used for the grid point interpolation whereas the DTM is stored as a regular raster and a separate line structure file is written.

A single or multiple OPALS Datamanger (ODM) files containing, both, point and line data constitute the main input for Module DTM (inFile). The DTM raster (outFile) is stored in a GDAL (GeoData Abstraction Library) supported format. The grid spacing (gridSize) should be chosen to match the laser point density. As a rule of thumb, approximately 1-4 laser points should be used for the height estimation of a single grid point. Specification of the size of the internal processing and storage units (tileSize) is optional. Whereas larger tile sizes improve the computation performance they also require more RAM storage resources.

## Data selection for non-triangulation interpolators

The neighbourhood definition is crucial for local surface fitting. Module DTM uses a generation based point selection strategy illustrated in Figure 1.

Fig.1: Generation based point selection

Data selection works as follows: The output grid structure is sub-sampled so that each output DTM grid cell is divided into four sub-sampling cells. The data points are then sorted into these sub-sampling cells. If more than one points falls within such a cell, a representative point is calculated based on weighted centre-of-gravity calculation. Therefore a user defined weight function can be specified (weightFunc). As a practical example for this, the user might want to favor points with a higher (raw) amplitude or (calibrated) reflectance over points with low signal strength. In this case the respective point attributes (amplitude, reflectance, intensity, etc.) can directly used as weight function but, however, arbitrarily complex weight functions can be specified in generic filter syntax. Once all the points are sorted into the respective subsampling cell, the generation based point selection starts with a search radius equal to the DTM grid size. In the 1st> generation all sub-sampling grid cells with a maximum distance to the investigated grid post of 1 * DTM grid width are considered. If the number of neighbouring points is less then the requested knn-count (neighbours) the search radius is extended by one sub-sampling grid cell (i.e. half of the DTM grid spacing) and all points of the circular ring are added. The procedure is repeated until either the knn point count is greater than or equal (neighbours) or the maximum search radius (searchRadius) is reached. Please note that, although this is a raster based point selection approach, still the arbitrary positions of the representative cell points are maintained. In addition to the point data, all line data in the vicinity of the grid point are considered within the interpolation.

## Data selection for triangulation based interpolation

Whereas moving least squares imply noise reduction and gross error filtering (at least to a certain degree), triangulations use the original unfiltered data for describing the interpolation surface. Hence, such surfaces contain measurement noise, gross errors and reflect the often random distribution of the input data. Therefore, grid models derived by moving least squares typically outperform tin based interpolators by precision and visual appearance. Nevertheless, tin based interpolators are well known and can easily support structure line information which is why Module DTM provides a Constraint Delaunay Triangulation (CDT) interpolator based on the CGAL library. For completeness, it is mentioned that the data selection parameters like neighbours, weightFunc, etc. do not influence the triangulated surface and therefore, they will be ignored during processing.

Based on the multi-threaded tile-wise processing strategy Module DTM provides two different triangulation strategies which will vary in processing speed and memory consumption:

• tile-based triangulation objects using a defined overlap (searchRadius > 0)
• a global triangulation object based on all relevant data (searchRadius = 0)

The tile-based triangulation strategy is activated by a searchRadius larger than zero which implicitly serves as overlap parameter (cf. Figure 2). Tile-based triangulation is typically faster and requires less memory since each processing thread has its own small triangulation object. This strategy, however, can result in interpolation artefacts due to differently formed triangles at the grid tile boundaries. Furthermore, large data gaps can lead to empty grid cell as visible in example 2. By increasing the searchRadius parameter which will lead to a large overlap, those issues can be eliminated or at least greatly reduced. Nevertheless, especially for data with inhomogeneous data distribution, interpolation artefacts may always occur. For such data sets the global triangulation strategy resolves those issues. There, one global triangulation object is created single-threaded which is then used by all interpolation threads. The downside of global triangulation strategy is that it requires enough memory to build a triangulation that contains all input data. So for huge data set it might be impossible to use the global triangulation strategy.

Fig. 2: tile-based triangulation with overlap (left), global triangulation (right)

All remaining parameters are optional. As mentioned above, the primary aim of Module DTM is the derivation of a terrain model. Furthermore, additional feature models (feature) can be derived as by-products of the grid interpolation (i.e. without loss of performance). The following additional features are supported:

• sigmaz: std.dev. of interpolated grid height
• sigma0: std.dev. of the unit weight observation
• pcount: point count
• pdens: point density estimate
• excentricity: distance grid point - center of gravity of data points
• slope: steepest slope [%]
• slpDeg: steepest slope [deg]
• exposition: slope aspect [rad] = azimuth of steepest slope line, N=0, clockwise sense of rotation
• normalx: x-component of the surface normal unit vector
• normaly: y-component of the surface normal unit vector

By default, the entire ODM content is used as input. Specific point/line data selection is enabled by providing an appropriate filter string (filter). In addition, a user defined limit can be specified (limit). By default, the height raster and the additional feature rasters are written to separate files. To generate a single multi-band raster file set multiBand=1. Please note, that all raster files are created in GeoTiff file format and the individual rasters are kept together by a GDAL Virtual Raster master file. Future releases will extend this concept to implement a hybrid DTM structure (i.e. joint storage of both raster and vector contents).

# Parameter description

-inFileinput data manager file name(s)
Type: opals::Vector<opals::Path>
Remarks: mandatory
The path(s) to the opals datamanager(s) whose point and line data are used to derive terrain model. In case multiple ODMs are specified, the data are queried from all (potentially overlapping) ODMs.
-outFileoutput DTM grid file name
Type: opals::Path
Remarks: estimable
Path of output DTM model file. Please note that the resulting DTM is stored as a hybrid raster. The regular grid of elevations and additional features are stored in GeoTiff format and the the structure information (spot heights, break lines, form lines, border lines) are stored in ESRI Shapefile format. If multiband output is enabled, a single raster output file containing the DTM grid and all features in separate bands (extension: '.tif') and a single vector file containing the stucture information (extension: '.shp') are created. If -multiband=false, a postfix is appended for each category (e.g. _dtm, _slope, _prec, _lines).
Estimation rule: The current directory and the name (body) of the input file are used as file name basis. Additionally, the postfixes ('_dtm', '_lines', ...) and extensions ('.tif' or '.shp') are appended.
-tileSizetile (block) size
Type: uint32
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 or finer grids are to be derived, the tile size may need to be adapted. Finer grids allow for larger tile sizes (and vice versa).
-gridSizeDTM grid width x/y
Type: double
Remarks: default=1
Size of the rectangular grid cell in x- and y-direction. If only one value is passed, quadratic cells are used.
-interpolationinterpolation method
Type: opals::GridInterpolator
Remarks: default=movingPlanes
Possible values:
delaunayTriangulation ... Delaunay triangulation interpolation
movingPlanes ............ moving (tilted) plane interpolation
kriging ................. kriging interpolation

Supported interpolation methods: movingPlanes|kriging.
-neighboursminimum nr of neighbours for grid interpolation
Type: int32
Remarks: default=0
The neighbour search is performed raster based in 'generations'. First, a regular raster featuring half of the final cell size is derived (storing the CoG of all points in each cell). All (raster) points within a search radius of the final grid size are taken into account first. If there are less valid cells available than the required number of neighours, the search radius is extendended by half of the grid size. This procedure is repeated until either the required neighbour count is exceeded or the maximum search radius is reached. In case the minimum neighbour count is set to zero, all points within th maximum search radius are considered for interpolation of the grid height.
Type: double
Remarks: estimable
Only points/lines 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 * grid size.
-weightFuncpoint/line point weight formula
Type: opals::Vector<opals::String>
Remarks: optional
An algebraic formula in generic filter syntax can be specified to calculate an individual weights for different classes of observations (bulk points, individual spot height points, soft structure line points (form lines), hard structure lines (break lines), and border lines. Please note, that additional point weights according to the distance of a data point w.r.t the grid post location are automatically calculated by the program to ensure a smooth terrain surface (i.e. inverse distance weighting).
Type: opals::Vector<opals::GridFeature::Enum>
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 [%] (deprecated: use slpPerc instead)
slpPerc ........ steepest slope [%]
slpDeg ......... steepest slope [deg]
exposDeg ....... slope aspect [deg] (azimuth of steepest slope)
normalx ........ x-component of surface normal unit vector
normaly ........ y-component of surface normal unit vector

If one or more features are specified, additional grid models in the same format and structure as the basic terrain model are derived. By default, a multiband raster is created in this case. Please note, that if each feature is written to a separate file (i.e. multiband=false) and 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
-filtermodify the input using a (tree of) filter(s)
Type: opals::String
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).
-limit2D clipping window
Type: opals::GridLimit
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).
GridLimit: Defines xy-limits for output grid/raster datasets. The limit coordinates 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.
Syntax: ['center'|'corner']['round']['(' left lower right upper ')']
-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>
-multiBandCreate multi-band output raster
Type: bool
Remarks: default=0
When activated the resulting DTM grid and the selected feature rasters will be stored in a single multiband raster file instead of as multiple files. The single feature bands appear in the order of specification in the multiband raster file.
-extrapolationCheckavoid extrapolation flag
Type: bool
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.

# Examples

The data used in the following examples are located in the \$OPALS_ROOT/demo/ directory. As pre-processing step, the 3D point cloud (flyover.laz) and the 3D break lines (flyover_breaklines.shp) are imported into a single ODM file using the following command:

opalsImport -inFile flyover.laz flyover_breaklines.shp -tileSize 100

## Example 1: DTM creation using moving planes interpolator

To create a 25cm DTM grid based using 12 nearest neighbours in a maximum search distance of 1.25m and additional feature models (surface normal vector and grid height precision) use the following command:

opalsDTM -inFile flyover.odm -gridSize 0.25 -neigh 12 -searchrad 1.25 -feature normalx normaly sigmaZ

Visualizations of the resulting raster files (DTM shaded relief map, color coded sigmaZ, normalx and normaly map) are plotted in Figure 3:

Fig. 3: (a) shaded DTM relief, (b) color coded height precision map, (c+d) color coded maps of normal vector components (c=x, d=y)

Figure 4 shows a detail of the shaded DTM relief compared to a variant without considering the breaklines in the surface interpolation (Module Grid):

opalsGrid -inFile flyover.odm -interpol movingPlanes -gridSize 0.25 -neigh 12 -searchrad 1.25 -outFile flyover_mpl.tif
Fig. 4: (a) shaded relief (opalsDTM), (b) shaded relief map (opalsGrid), (c) color coded height differences between the DTMs of 8a) and (b)

## Example 2: DTM creation using triangulation interpolator

In this example the triangulation interpolator is used and the effect of the two different triangulation strategies is shown. First, tile-based triangulation strategy using a 2m overlap (searchRadius -> 1) is applied and second, the global triangulation approach is used (searchRadius -> 0)

opalsDTM -inFile flyover.odm -outf tinLocal -interp delaunayTri -gridSize 0.25 -searchradius 1
opalsDTM -inFile flyover.odm -outf tinGlobal -interp delaunayTri -gridSize 0.25 -searchradius 0

Figure 5 shows shaded reliefs of the two different triangulation strategies. As it can be seen the two final DTM model only differ at larger data gaps (cf. Figure 4a). For completeness it is mentioned that gaps in raster models can be specifically closed by applying Module FillGaps.

Fig. 5: shaded relief of tile-based (left) and global (right) triangulation interpolator

Figure 6 visualizes the differences of the different interpolators and triangulation strategies in a side-by-side figure. Whereas the two triangulation strategies only differ at data gaps (Figure 6d), the height differences between triangulation and moving plane interpolator show speckle patterns (Figure 6e) especially on vegetated and non-sealed ground. Those speckle range from 2 to 4 cm and result from a higher filtering/smoothing effect of the moving plane interpolator. When closely comparing the corresponding shaded relief (Figure 6b and 6c) this smoothing effect is directly visible.

Fig. 6: differences of interpolations: shaded relief of (a) global triangulation, (b) tiled-based triangulation and (c) moving plane interpolator; height difference of (d) global minus tiled-based triangulation and (e) global tin minus moving plane interpolator

# References

Date
05.05.2022
@ movingPlanes
moving (tilted) plane interpolation
@ normalx
x-component of surface normal unit vector
Definition: GridFeature.hpp:23
opalsImport is the executable file of Module Import
Definition: ModuleExecutables.hpp:113
@ normaly
y-component of surface normal unit vector
Definition: GridFeature.hpp:24
opalsGrid is the executable file of Module Grid
Definition: ModuleExecutables.hpp:93
opalsDTM is the executable file of Module DTM
Definition: ModuleExecutables.hpp:48
@ feature
Use a statistic feature of the boundary gap points for filling.