The idea for this post began with the reading of an article about a henge in Northern Galicia (Spain). The finding of this neolithic structure is quite remarkable, since the henges were thought to be constructions belonging to the British Isles only. The discovery inspired me to develop a digital terrain model (DTM) from LiDAR data, including the henge and its surroundings, in order to find possible annex structures to the henge.

The LiDAR technology has been profusely utilized in archaeology, [2], [3], here it will be used with more modest means and reach, working with public access data and free software.

The LiDAR point cloud and the orthophoto have been donwloaded from the Spanish Centro Nacional de Información Geográfica website, and the GIS operations were conducted with QGIS (and its extension QGIS2Threejs). Quick Terrain Reader was used to peek into the LiDAR point cloud, FUSION, ALDPAT, LAStools, MCC-LiDAR and SAGA were the programs employed to do the point cloud calculations. The 3D rendering of the definitive surface was carried out with Blender.

Description of the work zone

A Roda Henge is situated in the Barreiros municipality, on the coast of Lugo (Galicia, Spain), 13 km west from Ribadeo. The structure is accesible by both A-8 highway or N-642 road.

Figure 1 shows a 3D model of the studied zone.

Figure 1. Three-dimensional model of the studied zone, made with FUSION and QGIS2Threejs.

From this original landscape visualization (or more properly, from the generating LiDAR point cloud, Figure 2.) I plan to recreate the Digital Terrain Model, which means to calculate a virtual representation of the naked land, omiting any elevation due to non-ground elements, such as trees, buildings, etc. In order to do that some filtering algorithms (filters) are required. Some of the more commonly used are included in the programs listed before.

Figure 2. Point cloud of the LAS archive of the studied zone.

The function of a filter is to select those points within the cloud that meet any desirable characteristic. In this case, belonging to the ground. Therefore returns from trees, buildings or vehicles are excluded. The henge, as part of the ground continuum, enters in the category of ground points, thus being represented in the digital terrain model.

Preliminary works

The model shown in Figure 1 has been created from a LiDAR cloud point (LIDAR-PNOA ceded by the Spanish © Instituto Geográfico Nacional). The file was cut with a 90000 m2 squared frame, centered in the henge, and non-landscape outlier points have been removed (those would be, for example, the returns from a group of birds). Using an orthophoto as background image QGIS2Threejs creates a elevation model from the Z coordinate of each point.

A very important work parameter is the return density of the point cloud. With the CATALOG command in FUSION is possible to know such characteristic. In our work zone the parameter reaches the 0.75 points per square meter (Table 1).

Table 1. Output extract of the CATALOG command.

This means that for every square meter the cloud point registers an average of 0.75 returns, or conversely, each return covers 1.33 m2. This parameter will be important to fix the cell size utilized by the filters, and will affect the definition of the generated raster file as well.

Once the starting data is refined and debugged, the filtering process begins.

Surface Filter based on the Kraus and Pfeifer algorithm (1998): Program FUSION

This algorithm is one of the most used in forested areas. It works with a robust function giving different weights to each point of the cloud and, depending on that weight, the point is included in a preliminary interpolated surface, or is discarded as non-ground. With successive iterations more points are included in the surface, each step getting more similar to the naked ground. The key is to find, for the studied zone, the optimal values of the a, b, g and w parameters that define the robust function. Determining those values, and applying a median filter afterwards give as a result the elevation model shown in Figure 3.

Figure 3. 3D raster result of applying the Kraus and Pfeifer filter (FUSION).

The filter is capable of dealing with the NW and SW forested areas, suppressing them, and giving overall good results. As a curiosity the filter cannot proccess the whole length of the viaduct that cross the highway, south to the henge.

Filter based on the Zhang et al. Progressive Morphological Filter (2003): Program ALDPAT

This filter, with a variable window size and a user-defined slope parameter, discards non-ground points with a changing threshold in successive passes. The result is shown in Figure 4.

Figure 4. 3D curved raster result of the Zhang et al. (2003) morphological filter (ALDPAT).

The results are also good but somewhat coarse, given it was not applied any smoothing filter.

Adaptive TIN Filter (Axelsson, 2000): program ALDPAT

This algorithm is also included in ALDPAT, and is also implemented in the TERRASCAN commercial software. It classifies a given point as ground or non-ground depending on some initial parameters of angles and distances to a TIN (Triangulated Irregular Network) of the ground. The Figure 5 presents the results.

Figure 5. DTM generated using the Adaptive TIN Filter of ALDPAT.

The behaviour of this filter is noticeably inferior to both the previous filters, with a incomplete result surface (some missing data in the south part of the DTM), and a worse processing of the highway side slopes. Some strange pattern-wise structures are also generated, not corresponding to actual ground.

Adaptive TIN Filter (Axelsson, 2000): program LAStools

The LAStools LiDAR data processing suite includes the program LAStools, which also works with the Axelsson adaptive algorithm. As one can see in Figure 6, the results are this time quite good.

Figure 6. LASground adaptive filter results.

This filter gives the best DTMs so far: it eliminates the forested areas, the highway side slopes are correctly represented and even the problem that the viaduct poses to the other algorithms here is dealt with success. It is worth noting that the program is very easy to use, with a fixed and very intuitive parametrization, and very adaptive to various situations (from dense forested areas to urban built-up zones).

The good results of the adaptive algorithm in this case suggests that either the ALDPAT program cannot handle the point cloud properly, or the initial parameters were not correct.

Multiscale Curvature Algorithm (MCC) (Evans and Hudak, 2007): program MCC-LiDAR

This algorithm is based in the VDF (Virtual Deforestation Algorithm, Haugerud and Harding, 2001) progressive densification filter, designed to extract the DTM in forested areas. From the user-defined scale and curve tolerance parameters, a preliminary surface is generated, testing each passing of the finter the inclusion (or not) of new cloud points. The final surface is the DTM. Figure 7 shows the results of applying the filter to the cloud point from A Roda.

Figure 7. Digital Terrain Model obtained from the MCC filter.

The processing of the forested areas is correct, and the viaduct appears partialy represented.

Morphological filter: Vosselman (2000) in SAGA

The Vosselman filter operates from a function that encloses the maximun admissible difference depending on the distance to the analyzed point. If the neighboring point falls into a threshold defined by a slope and a pass radius, then is considered to be ground and is included in the DTM. The resolution of this filter to the problem area is shown in Figure 8.

Figure 8. DTM generated from the Vosselman filter (2000) in SAGA.

Although the filter eliminates succesfully the vegetal cover, one can see certain lack of definition in the higway side slopes and the viaduct.

It seems so far that the filter with the best performance, considering as work requirements the elimination of the vegetal cover and the conservation of the bearing elements of the landscape, is the Axelsson’s Adaptive TIN Filter, calculated via LAStools.

Comparison with the Instituto Geográfico Nacional (IGN) raster

In Figure 9 a visual comparison between the LAStools calculated DTM and that of the Spanish Instituto Geográfico Nacional (DTM with 5 m mesh size, from the Plan Nacional de Ortofotografía Aérea LiDAR data) is set.

Figure 9. Instituto Geográfico Nacional DTM (top) versus LAStools DTM (bottom).

The model generated from the Axxelson’s adaptive algorithm seems to be more realistic regarding definition (because it was generated from a 1.5 m mesh size, versus the 5 m IGN mesh size) and ground representation (the LAStools model upholds the viaduct over the higway). The statistics of both rasters are presented in Table 2.

Table 2. Statistics comparison between the IGN raster and the LAStools raster. mean =mean of the heights, stdev=height standard deviation, min=minimal height, max=maximun height. Units in meters.

IGN mean stdev min max 112.12 7.93 87.31 130.82 LAStools mean stdev min max 112.38 7.872 87.3 131.04

The table presents very similar values, which should be tested in the field anyway. Considering that altimetric LiDAR estimations tend to underestimate altitudes, the higher LAStools estimate is a good indication.

Digital terrain model visualization

In this link the henge is presented in the center of a surface from the LiDAR point cloud filtered with LAStools, renderized with Blender. The frame is a 1 km side square.

https://sketchfab.com/models/c4fbe929e6a048e4b0798f20a17d3427/embed

Henge de A Roda

by atlasenxeneria

on Sketchfab

Conclusions

All the filters are able to distinguish the henge, and at the same time consider it as part of the ground continuum.

All the algorithms are able to filter those clound points corresponding to vegetal formations, and to generate the naked digital terrain model, in spite of the low resolution of the LAS files used.

The filters find difficult, in greater or lesser extent, to distinguish the abrupt changes in the landscape (big slopes, viaducts).

In a 1 km2 area no similar structures were found.

To know more

Ludwig Boltzmann Institute for Archaeological Prospection and Virtual Archaeology

Acknowledgements

To Carmen A., for her corrections and contributions.

Bibliography

Axelsson P., 2000. DEM generation from laser scanner data using adaptive TIN models. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIII, Part B4. Amsterdam.

Evans J. and Hudak A., 2007. A multiscale curvature algorithm for classifying discrete return LiDAR in forested environments. IEEE Transactions on Geoscience and Remote Sensing, Vol. 45, No. 4, April. New York.

Kraus K. and Pfeifer N., 1998. Determination of terrain models in wooded areas with airborne laser scanner data. ISPRS Journal of Photogrammetry & Remote Sensing 53, 193–203.

Vosselman G., 2000. Slope based filtering of laser altimetry data. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIII, Part B3. Amsterdam.

Zhang K., Chen S., Whitman D., Shyu M., 2003. A progressive morphological filter for removing nonground measurements from airborne LiDAR data. IEEE Transactions on Geoscience and Remote Sensing, Vol. 41, No. 4, April. New York.



This work is licensed under a Creative Commons Attribution 4.0 International License.