In this article you’ll find out how to apply two types of interpolation (TLI and IDW) simultaneously to a LIDAR point cloud, having segmented the cloud into two types of zone.
Tutorial HD LIDAR data processing with QGis
2- Download LIDAR HD data from IGN and load it into QGis
3- Tools for LIDAR data in QGis 3.32
4-Colorize a point cloud from an orthophoto
5-Colorize from an image with LAStools
6- Digital surface model (DSM) with CloudCompare and LAStools
7- digital terrain model (DTM) with CloudCompare
8- Creating a Digital Terrain Model with LAStools
9- Creating a Digital Feature Model with Open Lidar Toolbox.
10-Role of Interpolation in the context of the Digital Elevation Model (DEM)
Interpolation has a significant influence on the accuracy and visual quality of the final Digital Elevation Model (DEM). This process involves fitting a surface to point elevation data by superimposing them with a grid defined by a specific cell size. It is important to note that the points obtained as a result of this process are deduced from the original data and do not correspond directly to the actual data points.
Spatial interpolation: a major challenge
Available interpolators
Despite numerous efforts on the subject, interpolation remains a major challenge. Existing methods for generating DEMs from airborne LiDAR data encounter significant obstacles, particularly in project contexts where areas present dense non-terrestrial features or complex landscapes. The wealth of analyses aimed at assessing accuracy demonstrates the complexity of choosing the right interpolation method, as each method has its own specific advantages and disadvantages.
Dans la pratique, les méthodes fondées sur le krigeage sont réputées pour générer des MNE avec une plus grande précision, mais cela requiert davantage de ressources de calcul. Des approches telles que la pondération par la distance inverse (IDW), la méthode du plus proche voisin, et la triangulation avec interpolation linéaire (TLI) permettent d’obtenir des MNE plus précis. Notamment, la triangulation avec interpolation linéaire (TLI) et la méthode du plus proche voisin sont des méthodes simples et rapides pour créer des MNE relativement précis. Toutefois, leur performance diminue et leur sensibilité à la topographie s’accroît à mesure que la résolution spatiale du MNE augmente.
In practice, kriging-based methods are known to generate DEMs with greater accuracy, but this requires more computational resources. Approaches such as inverse distance weighting (IDW), the nearest neighbor method, and triangulation with linear interpolation (TLI) provide more accurate DEMs. In particular, triangulation with linear interpolation (TLI) and the nearest-neighbor method are simple, fast methods for creating relatively accurate DEMs. However, their performance decreases and their sensitivity to topography increases as the spatial resolution of the DEM increases.
The Spline method seems to offer a compromise between computation time and accuracy at high resolution. Although it is highly accurate at high resolution, its reliability decreases at low resolution.
The problems inherent in every interpolator
In contrast to the profusion of papers devoted to the interpolation of Digital Elevation Models (DEMs), the interpolation of archaeology-specific Digital Feature Models (DFMs) has so far only been briefly explored in a small number of studies, leading to divergent conclusions: TLI, Kriging, Spline Interpolation and the Nearest Neighbor Method have all been singled out as the most suitable interpolators.
To date, it is considered that despite the superiority of Kriging in terms of results, it is currently the Inverse Distance Weighted (IDW) interpolator that is specifically suited to archaeology, mainly due to its ease of access.
However, this does not mean that Inverse Distance Weighted (IDW) Interpolation is free of faults. In areas where sampling is grossly inadequate (where the density of measurement points is much lower than that of grid cells), IDW generates a relatively dense pattern of negative interpolation artifacts, resulting in errors in the form of exaggerated negative values. In moderately undersampled areas, the results will be visually more satisfactory, as no positive artifacts will be present. However, very small negative artifacts and a slight scaling effect (manifesting itself as fish-scale artifacts) on slopes may be observed. Although this effect is not particularly visually disturbing, it could mislead an inexperienced observer (or even an experienced one in the absence of adequate metadata) into falsely concluding that the formation is terraced.
In slightly undersampled areas, the IDW tends to present data noise resulting from the configuration of measurement points, which could be confused with geological or other features. In areas where the density of measurement points equals or exceeds that of the grid cells, IDW leads to excessive smoothing, which is not the optimal solution for interpolation adapted to the archaeological context.
Ordinary kriging stands out as a superior interpolator, yet it remains poorly accessible in most relevant software. If its accessibility were improved, it could quickly become the preferred option for archaeology-specific interpolation. On the other hand, the most widespread opinion highlights the potential of a hybrid interpolator that would merge the features of IDW and TLI. These two methods are distinguished by their high accessibility and moderate cost of implementation. While TLI is the best choice for over-sampled areas, IDW excels in under-sampled areas.
Setting up a hybrid interpolator
In the past, a key element was missing for the realization of such a hybrid interpolator. This was the lack of a suitable segmentation method to divide the data into zones better suited to each interpolator. Recently, this problem has been solved by using classification and regression tree analysis to generate a Digital Terrain Model (DTM) uncertainty prediction map. This method has been adapted, in the Open Lidar Toolbox plugin, to create an archaeology-specific confidence map for the Digital Feature Model (DEM), assigning a confidence level ranging from one to six to each grid cell. This approach makes it possible to quantify the quality of the Digital Feature Model at pixel level.
The confidence map
When evaluating the interpolators against the DEM confidence map, it was found that the Triangulation with Linear Interpolation (TLI) method performed best for areas with confidence levels five and six, and gave good results for levels three and four. Having solved the segmentation problem, a challenge remained in implementing the hybrid interpolator. By definition, the TLI and IDW interpolators produce slightly different grid models.
If these two models were simply combined using the raster calculator, artifacts in the form of steps would appear in the transition zones between the two interpolators. While this phenomenon is intrinsic to any hybrid interpolator and cannot be totally eliminated, it can be mitigated. Three distinct mitigation measures have been implemented in Open Lidar Toolbox.
Solutions for the hybridization of two interpolators
Firstly, transition zones can be reduced by a defragmentation operation. Indeed, under certain landscape and/or data collection conditions, confidence map values appear in the form of scattered fragments. In such circumstances, the transition zone grows exponentially, diminishing the positive effects of hybrid interpolation. This problem can be solved by imposing a minimum size for each interpolation zone, which would be much larger than the size of a single cell.
The zonal statistical method with a neighborhood of 11 cells was selected as the most efficient. This choice was motivated by the fact that this neighborhood guarantees the absence of zones smaller than 6 cells, thus ensuring consistency in the analyses.
It should also be noted that the differences between the TLI and IDW interpolators are not uniform. Their variations are most pronounced in confidence zones three and four, resulting in relatively small differences between the two interpolators in these confidence ranges. As a result, transitions between the two interpolators are smoother and less noticeable when confidence levels are between three and four.
In addition, steps were taken to mitigate the remaining effect of transitions between interpolators. One strategy adopted was the introduction of a buffer zone in the contact zones, where the elevation is calculated as the weighted average between the results of the TLI and IDW interpolators. Various experiments were carried out with buffers of one, three and five cells, involving weighted averages. In the majority of cases, the results were similar, although in a few extreme cases, the use of larger buffers produced more significant artifacts. As a result, the decision was taken to use a single-cell buffer to achieve an optimal compromise.
It became apparent that “doughnut-shaped” artefacts could form when confidence level 1 zones on the DFM map were directly adjacent to zones of level 4 or higher. To solve this problem, the contact zone was widened, i.e. moved further into the level four zones: a three-cell move proved to be an adequate solution for this situation.
How hybrid interpolation works
All the conditions required for hybrid interpolation have now been met. The process is as follows:
- First, the entire point cloud is interpolated using both interpolation methods.
- Next, the area is divided into segments by reclassifying the DFM confidence map into two distinct categories: “red” segments (corresponding to levels one to three for IDW) and “blue” segments (corresponding to levels four to six for TLI).
- Next, a defragmentation step is performed, followed by a displacement of three cells.
- In addition, a contact buffer zone one cell wide is defined. Within this buffer zone, the average elevation value of the results obtained by the two interpolators is used.
- The final phase of the procedure involves merging the “red” and “blue” segments, as well as the buffer zone. This is achieved using cartographic algebra techniques.
Algorithm for creating the confidence map
This is the processing model used by Open Lidar Toolbox in QGis to create the confidence map:
The confidence map is produced from three input files, two point clouds and one raster:
- the digital model raster calculated using the triangulation method (TLI)
- the density per square meter of points corresponding to low vegetation
- the density per square metre of points corresponding to the ground
Depending on the density and slope values calculated from the digital model, we’ll have six classes:
- Class 1: density of ‘soil’ points: <1 , Slope: <22.5° , Density of low vegetation points: not used.
- Class 2: Soil point density: 1-4 , Slope: 22.5°-42.5° , Low vegetation point density: not used.
- Class 3: Soil point density: 1-2 , Slope: <22.5° , Low vegetation point density: not used.
- Class 4: Soil point density: >2 , Slope: <22.5° , Low vegetation point density: not used, or Soil point density: >4 , Slope: >12.5° , Low vegetation point density: >4.
- Class 5: Soil point density: >4 , Slope: <12.5° , Low vegetation point density: >4 or Soil point density: >4 , Slope: >12.5° , Low vegetation point density: <4.
- Class 6: Ground point density: >4 , Slope: <12.5° , Low vegetation point density: <4.
IDW interpolation is used for zones ranked 1 to 3, and TLI interpolation for zones ranked 4 to 6.
Hybrid interpolation of elevation models
To apply hybrid interpolation to create a Digital Terrain Model, as we saw in the previous chapter, we use the One Step Processing tool in the Open Lidar Toolbox plugin.
To create a Digital Terrain Model (DTM)
To use the hybrid interpolation tool, we need..:
- a confidence map
- a raster interpolated with IDW
- an interpolated raster with TIN
The rasters created by the Open Lidar Toolbox plugin use the ‘soil’ and ‘buildings’ classes as terrain layers. You can modify the parameters in the plugin’s Python code, but it’s still easier to run the plugin to obtain the confidence map, interpolate only the ‘ground’ class and replace the raster created by the plugin with these new rasters.
To create the confidence map, we need a cloud of classified LIDAR points. It’s not enough to simply have the classified ground points, as we also need the lower vegetation. You can either start from a classified IGN LIDAR point cloud, or from an unclassified point cloud to which you can apply LASTools’ LASclassify, Open Lidar Toolbox’s Classify LAS/LAZ or any other classification process, depending on your needs and the type of data you’re using.
In the example below, we start with a classified LIDAR HD point cloud from IGN.
To create the confidence map, we can use the One Step Processing tool
Make sure you check The input Las/Laz file is already classified and uncheck all visualizations that are of no use.
The result is as follows:
You can use the IDW and TIN interpolation tools of your choice, but here we’ll keep the same tools as those used by the plugin:
- WhiteboxTools->Lidar Tools -> LidarIDWInterpolation
- WhiteboxTools->Lidar Tools -> LidarTINGridding
Please note that these tools don’t work with .laz files, so you’ll need to convert them to .las beforehand (Treatments->Point Cloud Conversion -> Convert format).
Enter all classes other than ‘soil’ (2) in the exclusion list and increase the Search Radius according to the ‘holes’ without soil points. If the Search Radius is too small, you’ll have holes in your IDW raster.
Now that we have the two rasters we need, we can run the Hybrid Interpolation tool:
We use the confidence map from One Step Processing, but the two rasters interpolated using exclusively ‘ground’ points.
Load the resulting DTM into QGis and, with the Elevation Profile View, you can see all three elevation models, especially in hybrid areas:
The blue line corresponds to the hybrid DTM. It’s easy to see that the results of IDW and TIN interpolations differ where terrain slopes are steeper.
It’s up to you to decide, in the field, which method is best suited to your objectives.
