STATE-WIDE CALCULATION OF TERRAIN-VISUALISATIONS AND AUTOMATIC MAP GENERATION FOR ARCHAEOLOGICAL OBJECTS

MAP GENERATION FOR ARCHAEOLOGICAL OBJECTS Frank Thiemann 1, Malte Schulze 1, Utz Böhner2 1 Institute of Cartography and Geoinformatics, Leibniz University Hannover, Germany (thiemann, schulze)@ikg.uni-hannover.de 2 Lower Saxony State Service for Cultural Heritage, Hannover, Germany utz.boehner@nld.niedersachsen.de Commission II, WG II/8


INTRODUCTION
Over 130,000 archaeological sites, from which 30,000 monuments are visible on the ground (mostly burial mounds, megalithic graves, enclosures, residential mounds, ramparts, etc.) are known in the area of Lower Saxony and documented in the archaeological information system. Until now, most objects are only documented by a rough position (about 10 m from freehand mapping) and using textual descriptions, photos and hand drawings (sketches). A manual inspection and surveying in the field of all monuments, which are mostly located in the forest, is not possible. With the availability of a state-wide airborne laser scanning dataset provided by LGLN (Lower Saxony's mapping agency) in 2018/2019 it became possible to explore these objects from distance. Many additional earthworks have been discovered from the ALS data. With automatic snapshots of the point cloud or of the digital terrain model (DTM) and with automatically generated maps it is possible to document the exact state, shape and extend of all known and visible historical man-made structures in the terrain. As part of open data Initiative of the state Lower Saxony, the catalogue of culture heritage is made also available to the public in the online portal "Denkmalatlas Niedersachsen" (denkmalatlas.niedersachsen.de). With the automatically generated visualisations, the monuments can be conveyed more clearly to the general public.

RELATED WORK
In the last two decades airborne laser scanning (ALS) became very popular for archaeological prospection and for documentation of archaeological features such as earthworks. An early project in archaeology starting in 2000 is described in Sittler (2004). LiDAR point density and precision were increasing with the rapidly advancing development of the technology. One big improvement was the full waveform analysis, which was evaluated for archaeological purposes by Doneus, M., et al. (2008). First archaeological investigations in Lower Saxony with ALS were made by H.-W. Heine (2010). In cooperation with the Institute for Cartography and Geoinformatics ten castles in the Weserbergland (Weser Hills) were mapped from ALS data. In early research projects, the classical DTM visualisations were used, namely greyscale and pseudo-coloured elevation and hillshading. To improve the visibility of small terrain structures new visualisation techniques were developed. Yokoyama, Shirasawa and Pike (2002) developed the positive and negative openness as an angular measurement for dominance or enclosure. Chiba, Kaneta and Suzuki (2008) combined openness and slope to the red relief image (RRIM), a more natural looking and easier to understand shading of the terrain. Hesse (2010) developed the local relief model (LRM), which is an improved high pass filter to separate the archaeological objects from the natural terrain.
Zaksek, Ostir and Kokalj (2011) developed the sky-view factor, which looks like a shading with diffuse light. An overview about the state-of-the-art terrain visualisation techniques for archaeological purposes can be found in Kokalj & Hesse (2017). In newest research, machine learning techniques are used to detect, segment or classify archaeological objects from high resolution DTMs or its visualisations. Kazimi et al. (2019) determine which DTM visualisation is most effective for the object detection using deep learning.

STATE-WIDE TERRAIN MODEL AND DERIVATIVES
The state-wide ALS data is given in the compressed LAZ-format and is partitioned in tiles of 1 km x 1 km. The LiDAR points are already classified into the classes ground, noise, synthetic water, underground, surface, misc. and overlap. The point-density varies, due to different years of acquisition, surveying providers and different terrain and vegetation, between four and nine ground points per square meter. The implementations were made in Python using the interface ArcPy and standard analysis functions of ArcGIS. Because ArcGIS is not supporting the LAZ format the opensource tool LasZip (Ilsenburg, 2013) and the ESRI tool LAS optimizer (EzLAS) are used to convert the data from LAZ via LAS to the proprietary ZLAS format. The ZLAS files are aggregated in a LAS-dataset. To create the DTM, we need only ground and synthetic water points. Overlapping points from neighbouring scan strips are not used, to avoid systematic patterns in the surface, due to small systematic errors of the georeferencing of the scan strips (see Figure 1).

Figure 1.
Pattern from non-perfectly registered overlapping scan strips (left) and without the overlapping points (right).
For all calculations the 1 x 1 km² tiling is kept and a rectangular buffer of at least eight pixels is added to each tile to avoid interpolation artefacts at the borders and to provide vicinity for the following filter operations. The border width depends on the point density of the scan and on the filter radii.
With linear interpolation inside a TIN, different DTMs can be calculated. Nearly full detail is maintained with 0.25 m resolution. In most cases 0.5 m resolution is enough for recognizing all relevant features of the objects, the space is reduced but the surface is smoothed. The TIN interpolation method is used, because it is most direct and missing data points will result in clearly visible big triangles. This enables the users of the DTM to recognize unreliable parts with a low point density or missing points (see Figure 2).  The raster tiles will be aggregated using raster mosaics in ArcGIS and virtual raster for GDAL and QGIS. Virtual raster may also be used in ArcGIS, but the handling in ArcGIS for large datasets is slow. From the DTM tiles, different visualisations can be calculated. Due to the partitioning with additional borders, the tiles can be processed completely in parallel. The standard visualisations are an RGB-image with hillshade from three directions (see Figure  9 bottom), a slope image (see Figure 10 top) and contour lines with an equidistance of one meter. The contour lines are also stored in raster format to keep the efficient raster tiling. Contour lines (TIFF) 0.6 4.6 Table 1. Hard disk space and processing time per square kilometre (for one kernel on an Intel ® Core™ i7 CPU with 3.2 GHz).
The selected visualisations are easy to understand also for nonexperts, can be used for a wide range of object types and implemented with ArcGIS standard functions with an acceptable runtime. It is possible to derive other visualisations from the DTM tiles like sky-view-factor using the open source software "Relief Visualisation Toolbox". After using external software, the resulting images will be clipped to the 1 km x 1 km size and be organized in a raster mosaic and virtual raster. The complete data flow is shown in Figure 4. Average disk space and processing time is listed in Table 1. We used 64 bit processing on an Intel ® Core™ i7 CPU with 3.2 GHz and 16 GB RAM. All processes can run in parallel on different kernels or The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) different machines by distributing the input files to the different processes. Reading and writing time take up a significant part of the total runtime. Using fast and separate hard drives for reading and writing and a RAM-drive for temporary files speeds up the process significantly.

Calculation of DTM representations and visualisations for individual archaeological objects
For the calculation process, two different input datatypes can be used: A LAS dataset or an already processed raster DTM. For the LAS dataset, a digital terrain model with a pixel size of 0.25 m is interpolated in a linear fashion using a Delaunay Triangulation. A bounding box for the actual object is derived and enlarged to include a certain amount of surrounding around the object. Also, for generating a square map window, the shorter side of the box is widened to form a quadratic shape and is then generally enlarged if its side length is shorter than a certain threshold. For the maps generated in this paper, a minimum quadratic side length of 40 m is given for the box. The calculation of the DTM is then limited to the extent of the derived bounding box plus an additional buffer of 10 m which is cut afterwards from the raster DTM to minimize interpolation errors on the border of the area. If the input data is an already processed raster DTM, it can be directly cut out using the bounding box without the buffer. The Process is shown in Figure 5. The DTM forms the basis for the calculation of several derivatives needed for the map generation. Slope, hillshade, RGB-hillshade, simple local relief model as well as north-south and east-west profiles are easily derived directly from the 0.25 m DTM. The whole data flow process is shown in Figure 6. To calculate the contour lines, a few more steps are needed in the process. Depending on the vegetation, classification errors and the roughness of the given terrain, the derived DTM can contain a very noisy surface. As noise in the elevation of the terrain model leads to very twisty and noisy contour lines, which is not desired in a map generation process, the DTM gets low-pass filtered by calculating the mean using a kernel with a circular structure. Then the calculation of the contour lines with an equidistance of 0.1 m takes place. The different shape of lines resulting from an unfiltered and filtered DTM can be seen in Figure 7. Depending on the remaining noise of the filtered terrain model, contour lines may contain small artefacts, shown as small circular structures. These lines lead to unnecessary increased line densities, can be too small for the map representation or might be an unpleasant visual effect. For the map-reader the elevation of these lines is often hard to identify, being too small to be labelled. These artefacts occur especially in flat surfaces, where the elevation alternates around a given contour line value.
To reduce these line structures, a filtering process based on the shape of the line is implemented: For every contour line, it is determined whether it forms a closed loop inside the map area. When a closed loop is detected, the length of the line as well as the area of the closed structure is calculated. In a second step, the ratio of the area is divided by the length of the structure and compared to a given threshold. With this, we consider not only small area but also slim, elongated structures. If it falls short to the threshold, the contour line is deleted. The result is shown in Figure 8 where the deleted artefacts are shown in a red colour. As the selected equidistance of 0.1 m leads to quite dense contour lines, which can be too dense especially in smaller scale representations, a process is implemented to select an appropriate line density by picking the lines with a derived higher equidistance. This is estimated by summing up the length of all contour lines inside the given area divided by the scale of the representation, which results in a value representing the length of all contour lines inside the map. Based on the total length of the contour lines compared to given thresholds, an equidistance of 0.1 m, 0.2 m, 0.5 m or 1 m is chosen for visualisation and the corresponding contour lines get picked accordingly. For the counting lines, by default, every fifth line is chosen. As flat surfaces of the terrain may lead to a low density of contour and counting lines, the number of counting lines inside the area is determined and compared to a given threshold. If the number of counting lines is too low, every second line is chosen. Based on the calculated DTM, its derivatives and additional cadastral data, which is extracted from a web map service, different map representations can be automatically generated.
These are derived from a basic layout, which is present in a map document (.mxd). The basic layout describes the size of the map field and the coordinate frame as well as a basic scale bar, textbox and logo. As manipulation of the coordinate frame (i.e. number of ticks) as well as the ticks for the scale bar is not directly possible using ArcPy, we set up four basic layout files for different representations scales (1 : 250, 1 : 500, 1 : 1,000 and 1 : 2,000). The matching representation scale is given by the size of the calculated bounding box. To generate a certain map representation, the desired derivatives are picked and added as new layers to the map document. For each derivate, a specified style is defined in a layer file (.lyr). The symbology of these files is applied to the newly created map layers, ensuring a predefined look for the data representation in the map. Finally additional information as scan acquisition date and identification number of the object is imported from different data sources and put in the text box together with author information to provide additional meta-information.

Derived map representations
Not every archaeological monument type can be recognised equally well in each visualisation. Therefore, we create a set of 10 maps for each site, from which the archaeologists can select a suitable one.  The classical hillshade is dependent on the direction of the light source. Structures parallel to the direction of light or in very dark or light areas may not be visible. Combining three hillshades from different directions with azimuths of 315°, 15°and 75°as the three channels of a RGB-image reduces the problem of directional dependency and is still easy to interpret (Figure 9). A slope image visualizes the change rate in elevation without using any additional parameters. It works in all terrains. The disadvantage is that the direction of slope is not directly recognizable -hills and depressions look the same and can be confused (Figure 10 top). For all three visualisations, we use percentage clip stretching for maximizing the contrast. Our hillshade uses a two-times vertical exaggeration and an altitude angle of 30°to increase the contrast for shallow objects.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) Figure 10. Barrow at a hillside visualized with slope (top) and with dense contour lines plus hillshade (bottom).
Historical terrains structures are often blurred by erosion, so the shape boundaries are hard to identify. Elevation and slope are smoothly changing which results in continuous changing of brightness or hue in the visualisation. As the human eye reception and recognition is not good recognizing absolute brightness, it is even harder to identify these boundaries in the visualisation. Using contour lines with a small equidistance like 10 cm, the slope is visualized as the density of the lines. The distance between the lines can be perceived quantitatively, so it is easier to define the border of a blurred earthwork. Concentric small closed lines also help to detect small hills or holes (Figure 10 bottom). Large regions of Lower Saxony are very flat. In these areas, the direct visualization of the elevation with min-max-stretched greyscale or pseudo-colours by using only greyscale (brightness) up and down is easy to understand. With pseudo-colours, more different elevation values can be recognized and it is easier to compare the elevation levels of different parts of the image but the colour-ramp needs to be explained. The recognition of elevation differences is depending on the brightness of colour.
In the pseudo colour map, even shallow forms get visible ( Figure  11).
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) Figure 12. Contour maps of iron age rampart Schnippenburg near Ostercappeln with SLRM (top) and RGB-hillshade and cadastral data (bottom).
Contour lines convey the absolute elevation. A shaded background should transport the smaller terrain features and support the spontaneous 3D impression. The Simple Local Relief Model highlights small structures or edges. (Figure 12 top) The filter radius controls the size of the structures to be highlighted. We are using a fixed radius of ten metres. Concave edges are dark, convex edges are bright. The grey values of the SLRM background are stretched by five times standard deviation and the brightness is increased by 30 % to ensure enough contrast between contour lines and background. In the RGB-hillshade map (Figure 12 bottom), also buildings, parcels and the border of the side can be added. The brightness of hillshade is increased by 20 %. The building polygons covers the regions where the natural terrain cannot be observed. The parcels and the shape of the site is used for revision purposes, to check if the object shape is correct and if the object lies completely in one parcel. Elevation profiles are cut through the centre of the map in northsouth and west-east directions ( Figure 13). The elevation is exaggerated to emphasize the elevation structure of an object. These profiles give a good impression of the actual shape of a specific object, as its elevation and form are directly perceptible from a side view representation. The vertical shape of a structure is harder to imagine in top down map views, where identification of an object and its extend is the main task.

RESULTS AND DISCUSSION
Small objects such as burial mounds, which often have a diameter of 15 metres and an elevation of only 80 cm, are well visible in a combination of RGB-hillshading and 10 cm contour lines. These lines clearly show the shape of the barrow, and old traces from illegal digging are easy to identify. Large objects in flat terrain, such as residential mounds, which have an extension of 50 -100 metres and often an elevation of only 1 -2 metres, are hardly recognisable in a hillshade visualisation. A combination of a stretched DTM and contour lines is the best choice for this. The DTM stretching can generally be used in the area of elevation differences smaller than 2 metres, e.g. also to make shallow former river courses visible. Otherwise, if the relief is too steep and on a slope, the archaeological structures are hardly recognisable. RGB-hillshading as well as classical hillshading are very well suited for identifying archaeological objects and structures on a larger scale as well as finding more structures in the vicinity of known archaeological monuments. Hillshading is also very well suited for the representation of barrow fields. During the evaluation of the ALS data, further barrows were discovered in almost every field. Smaller rampart structures are clearly visible in almost all visualisations. Their trace and upper edge can be identified in the sky view factor and SLRM. In the case of large fortifications, such as from iron age or younger, which are usually located in very hilly terrain, visualisations with higher equidistance contour lines in combination with RGB-hillshading show the archaeological object and its topographical position best. At bundles of paths, the surroundings must also be shown in addition to the actual object. For long linear features like land defences, it is also necessary to calculate additional areas of particular interest, such as gateways.

SUMMARY & OUTLOOK
The state-wide ALS dataset for the first time allows successfully documenting all prehistoric archaeological structures and existing damages on the still remaining monuments. In the future, the loss of monument substance can be automatically quantified and visualised by calculating volume-differences of renewed aerial surveys. Only a DTM with grid sizes of 0.25 m to 0.5 m depicts the archaeological structures with the required accuracy, an even higher point density leads to many pseudostructures and noise. The automation of the process in the proposed way is the great advantage of the method for monument conservation.
Further steps in this work are planned. At the moment, the extent of the map view is limited to a fixed size and a quadratic shape. In order to visualize elongated structures, it is planned to look for a partitioning of the objects, visualize each partition and include a small-scale overview map for the whole object. In a second step other visualisation techniques mentioned in the related work should be tested, e.g. the red relief model image or the sky-view factor, and evaluated for their ability to show small structures in all kinds of terrain. In addition, a systematic evaluation of the different visualisation methods and their parameters is needed to find an optimal set for each kind of object in a specific landscape. In order to automate this evaluation, a rating process is needed for automatic analysis of the visualisation for features like the quality of the background contrast, the smoothness of the contour lines or the optimal scale.