Prostor kraj čas 14 AIRBORNE LASER SCANNING RASTER DATA VISUALIZATION A Guide to Good Practice Žiga Kokalj and Ralf Hesse PROSTOR, KRAJ, ČAS PROSTOR, KRAJ, ČAS Airborne laser scanning raster data visualization: A Guide to Good Practice Žiga Kokalj Ralf Hesse Založba ZRC Ljubljana, 2017 PROSTOR, KRAJ, ČAS 14 AIRBORNE LASER SCANNING RASTER DATA VISUALIZATION A Guide to Good Practice Žiga Kokalj and Ralf Hesse Edited by: Nataša Gregorič Bon Photography: Žiga Kokalj, Jovan Cukut, Forestry Comission Scotland, Boštjan Laharnar, Lawrence Goldman, Brian Lockett, Frank Numrich, Michael Wilson Issued by: Institute of Anthropological and Spatial Studies, ZRC SAZU For the Issuer: Ivan Šprajc Publisher: Založba ZRC For the Publisher: Oto Luthar Reviewed by: Kevin Barton, Robert Shaw CIP – Kataložni zapis o publikaciji Narodna in univerzitetna knjižnica, Ljubljana 528.8.044.6:911.5(0.034.2) KOKAJ, Žiga Airborne laser scanning raster data visualization [Elektronski vir] : a guide to good practice / Žiga Kokalj, Ralf Hesse ; [photography Žiga Kokalj ... et al.]. - 1st e-ed. - El. knjiga. - Ljubljana : Založba ZRC, 2017. - (Prostor, kraj, čas, ISSN 2335-4208 ; 14) ISBN 978-961-254-984-8 (pdf) 1. HESSE, Ralf 288463616 First e-edition. Ljubljana 2017 © 2017, authors, Institute of Anthropological and Spatial Studies ZRC SAZU, Založba ZRC. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted, in any formor by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior permission of the publisher. The authors acknowledge the financial support from the Slovenian Research Agency (research core funding No. P6-0079 and project No. J6-7085). Front Cover Image: An anisotropic sky-view factor image of the Rhône valley. SRTM Data © USGS. 2014. Shuttle Radar Topography Mission, 1 Arc-Second Global data, Sioux Falls, South Dakota, USA. Endorsed by Archaeolandscapes International. International Abstract This guide provides an insight into a range of visualization techniques for high-resolution digital elevation models (DEMs). It is provided in the context of investigation and interpretation of various types of historical and modern, cultural and natural small-scale relief features and landscape structures. It also provides concise guidance for selecting the best techniques when looking at a specific type of landscape and/or looking for particular kinds of forms. The three main sections – descriptions of visualization techniques, guidance for selection of the techniques, and visualization tools – accompany examples of visualizations, exemplar archaeological and geomorphological case studies, a glossary of terms, and a list of references and recommendations for further reading. The structure facilitates people of different academic background and level of expertise to understand different visualizations, how to read them, how to manipulate the settings in a calculation, and choose the best suited for the purpose of the intended investigation. Key Words lidar, airborne laser scanning, visualization, interpretation, settings, techniques, tools, digital elevation model Izvleček Monografija nudi vpogled v nabor tehnik prikaza visokoločljivih modelov višin. Napisana je v kontekstu preučevanja in interpretacije različnih tipov zgodovinskih in modernih, kulturnih in naravnih majhnih reliefnih oblik. Daje jedrnate napotke za izbiro najboljših tehnik prikaza določenih tipov pokrajine in izrazitih oblik. Tri glavna poglavja – opis tehnik prikazovanja digitalnih modelov višin, napotki za njihovo izbiro in orodja za izračun prikazov –, spremljajo izbrani primeri tipičnih arheoloških in geomorfoloških študij, slovarček pojmov ter seznam literature in priporočenega branja. Posameznikom z različnih znanstvenih področij in z različnim predznanjem o tematiki je struktura v pomoč pri razumevanju različnih tehnik prikazov, kako jih brati, kako izbrati prave nastavitve pri njihovem izračunu in kako prepoznati najbolj primerne za namen zasnovane raziskave. Ključne besede lidar, aerolasersko skeniranje, vizualizacija, prikaz prostorskih podatkov, interpretacija, nastavitve, tehnike, orodja, digitalni model višin Contents 1 Introduction to the Guide 10 1.1 Aim and scope of this guide 12 1.2 How to use the guide 12 1.3 Applications of lidar visualizations 13 1.4 Free data access 13 2 Description of techniques 14 2.1 Analytical hillshading and hillshading from multiple directions 16 2.2 Slope 19 2.3 Elevation differentiation 20 2.4 Trend removal and local relief model 20 2.5 Sky-view factor and anisotropic sky-view factor 22 2.6 Openness 24 2.7 Local dominance 25 2.8 Cumulative visibility 25 2.9 Accessibility 26 2.10 Multi-scale integral invariants (MSII) 26 2.11 Laplacian-of-Gaussian (LoG) 27 2.12 Visualizing uncertainty 27 3 Guidance for selection of techniques 30 3.1 Preparing the images for detection and interpretation 32 3.2 Reading and presenting the images 34 3.2.1 Choosing appropriate visualisations 34 3.2.2 Personal preferences and intercomparability 35 3.2.3 Perception 35 3.3 Visualization of datasets other than lidar 37 3.4 What’s in a name? 39 4 Tools 46 4.1 Relief Visualization Toolbox (RVT) 48 4.2 Lidar Visualisation Toolbox (LiVT) 48 5 Case studies 50 5.1 WWII anti-glider defences at Culbin, Scotland 52 5.2 Hollow ways at Volčji Potok, Slovenia 54 5.3 Prehistoric and Roman settlement above Knežak, Slovenia 56 5.4 Barrow cemetery in Pivola, Slovenia 58 5.5 Alluvial fan at Craig canyon, Saline Valley, California, USA 60 5.6 Geological features in the Julian Alps, Slovenia 62 5.7 The Granite Dells, USA 64 5.8 Cinder cones and lava flows on Mauna Loa, Hawaii, USA 66 5.9 Blind valley Odolina, Slovenia 68 5.10 Lake floor morphology, Lake Constance, Germany, Switzerland, and Austria 70 5.11 High-rise buildings in New York, USA 72 Glossary of terms and abbreviations 74 Bibliography and further reading 78 List of figures 83 List of tables 88 10 1 Introduction to the Guide 11 12 1. Introduction to the Guide This guide aims to help specialists and interested public to produce or ask for such lidar products that will facilitate ‘reading and exploring’ 1 the landscape for meaningful information. 1|1 Aim and scope of How to use the guide this guide 1|2 introduction to a variety of processing methods that help create meaningful images for observation, exploration, and interpretation of rasterized lidar The aim of this Guide to good practice The guide is aimed at everyone who data. A special chapter describes how in airborne laser scanning (ALS) raster is interested in visual exploration of to recognize and visualize uncertainty data visualization (hereafter called ‘the rasterized elevation data (DEM) which can in data. guide’) is to help specialists and interested be the product of airborne lidar or other • The guidance is meant to lead the public to produce or ask for such lidar techniques. Scientists, professionals and keen explorer through the main steps products that will facilitate ‘reading and the public are all fascinated by the details of the process from preparing the exploring’ the landscape for meaningful lidar data offer. While automatic feature data to answering all the relevant information.This guide provides an insight detection is gaining importance, it is still questions regarding the selection into a range of visualization techniques for indispensable that a human observer asks, of techniques that are best suited high-resolution digital elevation models thinks about, and answers the questions. for the purpose of the intended (DEMs) with their specifics, advantages Visual exploration and inspection is the investigation. and weaknesses. It is provided in the only way to comprehend the existence, • The tools section briefly explains context of investigation and interpretation intricate relations, and context of small- and compares the main settings and of various types of historical and modern, scale structures in lidar data. The workflows of the freely available RVT cultural and natural small-scale relief guide is structured so that people of and LiVT toolboxes for calculation features and landscape structures. It also different academic background and level of visualization techniques. Anyone provides concise guidance for selecting of expertise can understand different wishing to visualize their digital the best techniques when looking at a visualizations, how to read them, how to elevation model (DEM) with one or specific type of landscape and/or looking manipulate the settings in a calculation, several of the techniques described in for particular kinds of small-scale and which techniques to use to help the chapter Description of techniques structures. It does not, however, provide answer the questions. Throughout the should be able to do so with these insight into other, equally important, areas guide, reference is made to more detailed tools, regardless of their current of lidar data acquisition, processing, and explanations, comparisons, observations, knowledge in image processing. management. While the examples mainly and further studies. relate to archaeological relief features on scales from several metres to a few The guide is comprised of three main hundred metres, the techniques presented sections accompanied by examples in this guide are equally suitable for the of visualizations of small-scale visualisation of many other features on features, exemplar archaeological and scales from microscopic to continental. geomorphological case studies, a glossary Furthermore, while the focus is on DEMs of terms, and a list of references and based on airborne lidar, the techniques recommendations for further reading. are applicable for any DEM-like raster data set, e.g. based on Synthetic Aperture • Descriptions of visualization Radar or Structure-from-Motion (SfM). techniques provide a brief 1. Introduction to the Guide 13 1|3 Applications of lidar Free data access visualizations 1|4 files in a TIFF or other similar geospatial format. More and more classified point clouds are also shared. Table 1 lists some of the possible sources where practice Digital elevation models are raster Getting lidar data for research purposes is experimental data (and more) can be datasets containing elevation values. becoming increasingly simple as a growing accessed. Some portals are easier to use Because these numerical datasets number of science groups, governmental than others, with the more user friendly are not readable as such, visualisation agencies, city councils, counties, regions, being those of Slovenia, England, and techniques are necessary to convert the and even countries publish their data with Wales. DEM into human-readable greyscale a free and unlimited use license. Datasets or colour images. The obvious use of come in various formats and states of different visualizations is for visual processing. Gridded digital surface and examination of data. They facilitate the terrain models are usually available, either interpretation of terrain data and are used as ASCII text files or already rasterized for visual validation of the geomorphologic quality of the DEMs e.g. in geography, Table 1: Some useful resources to go search for free lidar datasets. Host sites will inevitably geomorphology, cartography, archaeology, change and new datasets are becoming available regularly, so use your favourite search engine. hydrology, glaciology, forestry, and disaster management. Often we use them for a list of sources Terrain Data · openterrain Wiki · GitHub discovering new features, especially in https://github.com/openterrain/openterrain/wiki/Terrain-Data archaeology, but also to better delineate Denmark Agency for Data Supply and Effeciency http://download.kortforsyningen.dk and more precisely position already England Environment Agency known objects, such as riverbanks, levees, http://environment.data.gov.uk/ds/survey/#/survey cultural terraces, parcel divisions, stone Finland National Land Survey of Finland open map data download walls, or areas of erosion. However, ftp://tiedostot.kartat.kapsi.fi/laser/etrs-tm35finn2000/mara_2m techniques, particularly the ones directly Netherlands Dutch National SDI (PDOK) related to physical quantities, are also https://www.pdok.nl/nl/producten/downloaden-van-data-pdok used for different other purposes. Sky-http://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_gefilterd.xml https://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_05m_int.xml view factor, for example, is used in Norway Høydedata meteorology, because it is fundamental https://hoydedata.no/LaserInnsyn/ for modelling solar insolation, and it can Spain Plan Nacional de Observación del Territorio be used even in engineering applications http://pnoa.ign.es/coberturalidar such as predicting the availability of the Slovenia Ministry of the Environment and Spatial Planning. GPS signal in urban areas. Techniques http://gis.arso.gov.si/evode/profile.aspx?id=atlas_voda_Lidar@Arso that remove the general topography, Sweden Lantmäteriet http://www.lantmateriet.se/en/Maps-and-geographic-information/Elevation-data-typical examples are local relief modelling USA OpenTopography and openness, are also well suited for http://www.opentopography.org automatic feature extraction. USGS http://earthexplorer.usgs.gov While this guide focuses on raster ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/NED/LPC/projects elevation data, it is worth mentioning NOAA - Coastal Topographic Lidar the emergence and fast adaptation https://coast.noaa.gov/dataregistry/search/collection/info/coastallidar of pseudo and real 3D techniques for NASA - G-LiHT: Goddard's LiDAR, Hyperspectral & Thermal Imager http://gliht.gsfc.nasa.gov/ext/maps/index.html presentation of lidar point clouds or Search also for data of individual states as many provide them integrated generated 3D models. These offer a whole into their own distribution systems. new world of opportunities, not only for Wales Natural Resources Wales visual examination, but even more so for http://lle.gov.wales/Catalogue/Item/LidarCompositeDataset/?lang=en scientific exploration. 14 2 Description of techniques 15 16 2. Description of techniques The explorational value of different visualizations varies significantly with respect to characteristics of structures observed; their size, shape, 2 orientation, concavity or convexity, degree of prominence, and edge type. By transforming the DEM into greyscale 90° are dark (Figure 1A). Under very low Analytical or colour images, different types of 2|1 light source angles (below 10°) features of hillshading and visualizations help us make sense of extreme subtleness can be made visible, lidar data. Some visualization techniques hillshading from but this is useful for exploratory purposes produce results that can be related to multiple directions in areas with only slight variations in physical quantities (e.g. sky-view factor – general topography (Table 2, Figure 2). SVF, and local relief model – LRM) while Relief shading (also known as hillshading others only have a presentational value or shaded relief) provides the most Areas facing directly towards or away (e.g. principal components of hillshading ‘natural’, i.e. intuitively readable, visual from illumination are commonly from multiple directions). Some work impression of all techniques. It is a saturated (homogeneously bright or dark, better on an almost flat terrain (elevation description of how the relief surface respectively), and no or very little detail differentiation), or topography with gentle reflects incoming illumination based on can be perceived there. A single light beam slopes (LRM), while others on rugged or physical laws or empirical experience. also fails to reveal linear structures that mixed terrain (SVF). Their explorational There are numerous analytical hillshading lie parallel to it, which can be problematic value varies significantly with respect to techniques (Horn 1981; Blinn 1977; Batson in some applications, especially in characteristics of structures observed; et al. 1975; Minnaert 1961), although only archaeology (Devereux et al. 2008). their size, shape, orientation, concavity the method developed by Yoëli (1965) or convexity, degree of prominence, has become a standard feature in most Producing multiple relief shading outputs and edge type. Not all methods can be GIS software. Therefore, when analytical by illuminating a surface from multiple recommended for interpretation when hillshading (or hillshading in general) is directions enhances visualization of relative or absolute comparisons need discussed it is the method developed by topography (Figure 1). Hill-shaded images to be judged. This is because most can Yoëli that is implemented. are sometimes used to guide ground change the observable properties of surveys, but comparing multiple images in structures, depending on the settings used Standard analytical hillshading is easy to the field is extremely inconvenient. A step in the calculation. The DEM can be either compute and straightforward to interpret towards an improved understanding of the a subset of the digital terrain model (DTM) even by non-experts and without training. results is combining multiple shadings by (representing the elevations of bare- It has a basic assumption that the relief considering only the mean (Hobbs 1995), earth, in the text we regard it as such) or a is a Lambertian surface illuminated by the maximum, or the range of values, for digital surface model (DSM) (representing direct light from a fictive light source at each pixel. In order to display the areas elevations of terrain or objects on it, an infinitive distance (the light beam has with a low range of values more clearly, whichever the higher), or anything a constant azimuth and elevation angle the result can be square rooted. between the two. See the glossary of for the entire area). The computed grey terms (p. 76) for further details. value is proportional to the cosine of A common example is also a combination the illumination incidence angle on the of standard hillshading (315° azimuth) relief surface – the angle between the with vertical illumination (Imhof 1982; surface normal and the light beam. Areas Hobbs 1995). Hill-shades from three perpendicular to the light beam are the different directions can be used to create most illuminated, while areas with an an RGB colour composite image by incidence angle equal or greater than assigning the three greyscale images 2. Description of techniques 17 to the red, green and blue colour channel (Devereux et al. 2008; Hobbs 1999) and this often works well indeed (Figure 3A). Because images created by illumination from several angles are highly correlated (the same scene is viewed), it is possible to ‘summarize’ information by a mathematical transformation with principal component analysis (PCA) (Devereux et al. 2008). The first three components computed from multiple (e.g. 16) directions usually contain a high percentage (typically over 99 %) of the variability in the original dataset. They can thus be expected to provide a basis for the visualization with minimal loss of small-scale features. The PCA – especially the combination of the first and second principal components (Figure 3B), or the RGB colour composite image of the first three (Figure 3C) – simplifies the interpretation of the multiple shading data. However, it does not provide consistent results with different datasets. Figure 1: Angle dependence of analytical hill-shading: 315° azimuth illumination (A), and 225° azimuth (B), both with 35° Sun elevation. Note the difference of the relief perception and the structures that can be observed. Overgrown cultural terraces near Koboli in Slovenia as evidenced by 1 m spatial resolution terrain model. Data © ARSO, Slovenia. See also Figure 23. 18 2. Description of techniques Figure 2: Very low light source angles expose features of extreme Figure 3: An RGB image of hill-shadings from three directions (315°, 0°, subtleness: a standard 45° Sun elevation (A) and (B), and low light 5° and 45° azimuth with 35° Sun elevation) (A), a composite of the first two Sun elevation (C), all with 45° azimuth. However, this only works in components (B) and an RGB composite of the first three components areas with very gentle relief morphology, such as this example of the of a principal component analysis of analytical hill-shading from 16 Site A embanked enclosure(s) in Brú na Bóinne World Heritage Site directions with 35° Sun elevation (C). A late Roman camp at St. Helena, in Ireland. 1 m spatial resolution lidar data used with permission of west of Kobarid, Slovenia. 0.5 m spatial resolution lidar data © Walks of Meath County Council and the Discovery Programme. Local histogram Peace in the Soča river Foundation. saturation is used to present (B) and (C). The first to show the difference this makes when compared with normally presented shaded relief (A), and the second because the image is otherwise too dark to expose any details. A hint of a small circular enclosure can be observed to the NW of the bigger one on (C). 2. Description of techniques 19 Table 2: Typical settings for calculation and visualization of hillshading. visualization parameter general very flat terrain steep or complex terrain Sun azimuth [°] 315 315 315 Sun elevation angle [°] 35 5 > 45 recommended histogram stretch linear stretch, 2 % cut-off linear stretch, 2 % cut-off linear stretch, 2 % cut-off Table 3: Typical settings for calculation and visualization of principal components analysis of hillshading images from multiple directions. visualization parameter general very flat terrain steep or complex terrain Sun elevation angle [°] 35 5 45 number of directions 16 16 16 recommended histogram stretch linear stretch, 2 % cut-off linear stretch, 2 % cut-off linear stretch, 2 % cut-off 2|2 Slope Slope (gradient) is the first derivative of a DEM and is aspect independent. It represents the maximum rate of change between each cell and its neighbours and can be calculated either as percentage of slope or degree of slope. Challis et al. (2011) found slope the best visualization technique for most circumstances among the methods they analysed. If presented in an inverted greyscale (steep slopes are darker), slope retains a very plastic representation of morphology (Table 4). However, additional information is needed to distinguish between positive/convex (e.g. banks) and negative/concave (e.g. ditches) features since slopes of the same gradient (regardless of rising or falling) are presented with the same colour (Figure 4). Figure 4: Slope image of Žerovinšček Iron Age hillfort near Bločice, Slovenia (A). The elevation profile (B) refers to the P -P line in (A). Note that the structure the profile line crosses can be 1 2 quite easily – and wrongly – interpreted as convex instead of concave. 1 m spatial resolution lidar data © ZRC SAZU, Slovenia. Table 4: Typical settings for visualization of slope. visualization parameter general very flat terrain steep or complex terrain recommended histogram stretch* linear stretch, 0-50° linear stretch, 0-15° linear stretch, 0-60° *Inverted grayscale (white to black) presentation works best 20 2. Description of techniques 2|3 Elevation Trend removal and A differentiation 2|4 local relief model Elevation differentiation (also referred Archaeological features are generally of to as colour coding, colour cast, or a much smaller scale than the landforms constrained colour ramps method) on which they lie. It is therefore necessary controls the range of values that are to adjust visualization techniques presented over a given range of greyscale appropriately, for example defining a tones or colours. The technique applied is 2 small search radius for sky-view factor a histogram stretch, whereby the range of calculation or set a suitable range for elevation values of interest is stretched to elevation differentiation, although this the whole range of greyscale (or colour) is not possible with all techniques. A values of the resulting image (Figure 5, 1 procedure that separates local small- Table 5). This can strongly enhance the 263 m 536 m scale features from large-scale landscape contrast between pixels or areas of B forms is called trend removal. When different elevation. 271 m 278 m working with a DEM, the trend (i.e. the larger landscape forms) is represented by There are many variants of histogram Figure 5: (A 1) A histogram of the whole a smoothed (generalized) version of that stretch that can be applied to enhance data range – elevations between 263 m DEM. The smoothing can be accomplished contrast and emphasize detail. Besides and 536 m. Because we are only interested by applying a low pass convolution filter. the simple linear stretch between in a very narrow range of values between A Gaussian low pass filter produces a 271 m and 278 m – marked by a dark arrow lower and upper cut-off values, these smoother transition between features, on (A) –, it makes sense to stretch only this include nonlinear enhancements such range to the whole ‘colour’ palette (B). (A 2) but is computationally more demanding as logarithmic stretch, square root A histogram of stretched values of the area (Reitberger et al. 2008; Wagner et al. enhancement, exponential stretch, and shown on Figure 5. Instead of everything 2008). Trend removal is then accomplished histogram equalization. However, if the being concentrated in a ‘dark corner’ of the by subtracting the smoothed DEM from the histogram, the values are much more evenly preservation of the relative differences original DEM. The resulting difference map distributed and differences in colour shade between the values (elevations) is can be more easily perceived. contains only the local deviations from the important, the basic method known as overall landscape forms. linear stretch with saturation (cut-off of extreme values in the upper and lower parts of the histogram, also known as histogram clipping) is preferable. The histogram is linearly stretched to fill the whole range of values between the defined minimum and maximum value. For more information on how to use histogram stretch for visualization see section 3.1 Preparing the images for detection and interpretation. Elevation differentiation is very useful for visualization of features of interest in flat landscapes and is very easy to interpret, especially when an appropriate colour ramp is used (Figure 6). It is also the only visualization technique presented that retains the original information about relief elevation. It is therefore easy to assess factors such as the depth of ditches or height of tumuli. However, even with Figure 6: Histogram stretch to a narrow range of values. Past riverbeds of Nadiža. 0.5 m spatial resolution lidar data, West of Kobarid in NW Slovenia. 0.5 m spatial resolution lidar data © Walks slight variations in the general morphology of Peace in the Soča river Foundation. of terrain, the technique becomes less useful, because archaeological Table 5: Typical settings for displaying elevation differentiation. earthworks are obscured by the variation in topography and because intensive visualization parameter general manipulation of the histogram is required. recommended histogram stretch linear stretch with minimum and maximum cut-off For the same reasons the technique completely fails in rugged terrain. 2. Description of techniques 21 Because small-scale features are smoothed rather than eliminated by the low-pass filter, the derived trend removal map is biased towards small features, i.e. the local relative relief elevations are progressively underestimated as the spatial extent of the features increases. Hesse (2010) has therefore proposed to refine the process by introducing a ‘purged DEM’. This is produced by creating zero contours in the trend removal map (i.e. find all points for which the values of smoothed DEM and the original DEM are equal). DEM elevation data are assigned to all points along these contours. A new approximation of the generalized DEM is then interpolated from the points to create a purged DEM. Finally, the difference map between the original DEM and the purged DEM is a local relief model (LRM). The LRM derived using this approach results in a less biased representation of small- scale topographic features and reflects more truthfully the relative heights and depths of these features with respect to the surrounding landscape (Figure 7). The level of smoothing is defined by a kernel size of the low pass filter, where a smaller kernel exposes smaller features and vice versa. The precise kernel size should reflect the size of the small-scale landforms while a generally safe bet is a kernel size of about 25 m. The method works best on terrain with gradual slopes, while it produces artefacts such as artificial banks and ditches where relief is diverse and/or changing abruptly. Figure 7: A hill-shaded image (A) and a histogram stretch (B) of a local relief model, and the LRM presented with a purposely designed colour ramp overlaid on a hill-shaded image (Hesse 2010) (C). 1 m spatial resolution lidar data of Žerovinšček Iron Age hillfort near Bločice, Slovenia © ZRC SAZU. 22 2. Description of techniques 2|5 Various software solutions provide et al. 2012). This brings back some of the Sky-view factor and different algorithms to compute the ‘plasticity’ of hill shading and gives better anisotropic sky-view SVF. The difference between results, details on very flat areas. factor especially for visualization purposes, can be enormous. For example, the code given Both, SVF and anisotropic SVF are very Sky-view factor (SVF) can be used as by ZRC SAZU (2010) that is implemented good general visualization techniques, an alternative method of relief mapping in LiVT and RVT (see Chapter 4 Tools) has because they enhance visibility of simple in order to overcome the directional no saturations, while the implementations and complex small-scale features problems of hillshading (Kokalj et al. in SAGA GIS (Conrad et al. 2015) usually whatever their orientation and shape, on 2011; Zakšek et al. 2011). SVF is a give very saturated areas with low SVF. most types of terrain. geophysical parameter (if we do not This means no details can be perceived in manipulate elevation data by vertical valleys. Omitting the closest neighbour pixels exaggeration) that represents the portion in the SVF calculation greatly improves of the sky visible from a certain point Anisotropic (directional) SVF considers an the visibility of archaeological features (considering a hemisphere centred above unevenly bright sky. The brightness can where the data are noisy due to strip each pixel and ignoring any direction depend on the azimuth and solar distance misalignment or overambitious resolution below the mathematical horizon). Diffuse from the imaginary light source (Zakšek setting. solar insolation rasters can be used to visualize archaeological features as well (Challis et al. 2011), but require additional calculations and the results are more generalized. In contrast to shading techniques based on directional illumination, features visualized by SVF (or by openness) do not contain any horizontal displacements. SVF has a range between 0 and 1. Values close to 1 indicate that almost the entire hemisphere above the pixel is visible, which is the case in exposed features (planes, ridges, and peaks), while values close to 0 are present in deep sinks and lower parts of deep valleys where almost no sky is visible. Like hillshading images, SVF images are intuitively readable. While hillshading simulates directional illumination of a surface (‘sun’), SVF simulates diffuse illumination (‘overcast sky’). The computation of SVF is influenced by the search radius of the horizon – the larger the search radius, the more generalized the results. In contrast, a small search radius can be used to visualize and classify local morphological forms. For example, a 10 km search radius can be used in meteorological studies, while a 10 m search radius is suitable for discrimination of archaeological features. Locally flat terrain, ridges and earthworks (e.g. building walls, cultivation ridges, burial mounds) which receive more illumination are highlighted and appear in light to white colours on a SVF image, while depressions (e.g. trenches, moats, ploughing furrows, mining pits) are dark Figure 8: Sky-view factor image (10 m search radius in 16 directions) (A) and anisotropic sky-view factor image (B) of Žerovinšček Iron Age hillfort. Many details can be perceived despite the because they receive less illumination. variable relief morphology. 1 m spatial resolution lidar data © ZRC SAZU. 2. Description of techniques 23 Table 6: Typical settings for trend removal and Local Relief Modelling. visualization parameter general very flat terrain steep or complex terrain filter radius [m] 10 5 25 recommended histogram linear stretch, linear stretch, linear stretch, stretch* minimum -1.0 m, minimum -0.5 m, minimum -2.0 m, maximum 1.0 m maximum 0.5 m maximum 2.0 m *Depends also on the height and depth of the observed features. Figure 9: Sky-view factor image (10 m search radius in 16 directions) of ridge and furrow near Neudingen, Germany (A). Omitting the nearest pixels from the calculation can reduce the noise in the data so that features of interest can be perceived more clearly (B).1 m spatial resolution lidar data © LGL in Baden- Württemberg. Table 7: Typical settings for calculation and visualization of sky-view factor. visualization parameter general very flat terrain steep or complex terrain radius [m] 10 10 10 number of directions 16 16 16 recommended histogram linear stretch, linear stretch, linear stretch, stretch minimum 0.65, minimum 0.9, minimum 0.55, maximum 1.0 maximum 1.0 maximum 1.0 24 2. Description of techniques 2|6 Openness Openness is another proxy for diffuse relief illumination. The method is based on estimating the mean horizon elevation angle within a defined search radius (Yokoyama et al. 2002). Positive openness equals the mean zenith angle of all determined horizons, while the negative openness is based on nadir angles. Openness does not limit the estimation of each zenith angle by the mathematical horizon (as SVF does). In other words, openness considers the whole sphere and not only the celestial hemisphere. As a result, the maximum value of openness can be greater than 90°. In addition, a plane (a long slope or a horizontal plane) without any obstacles will always have an openness value of 90° irrespective of its slope. Therefore, interpreting openness results is sometimes difficult because a slope is visualized in the same manner as a horizontal plane. However, because openness considers the whole sphere for calculation, the result is a much ‘flatter’ image, devoid of general topography – a kind of trend-removed image. It has the same valuable properties for visualization as the SVF with the exception that the visual impression of the general topography is lost. Interpretation is therefore a bit trickier, but openness has a big advantage for automatic feature detection because ‘signatures’ of features are more homogeneous because they are Figure 10: Positive (A) and negative (B) openness image (10 m search radius in 16 directions) the same irrespectively of their location on of a late Antiquity settlement of Tonovcov grad, Slovenia. The very complex terrain seems flattened. Tops of protruding features are well delineated on positive openness image, while a plane or slope (Doneus 2013). negative openness delineates the bottoms of hollows and lower edges of cliffs. 0.5 m spatial resolution lidar data © Walks of Peace in the Soča river Foundation. Negative openness is not the inverse of positive openness and it provides Table 8: Typical settings for calculation and visualization of positive and negative openness. additional information. While positive visualization parameter general very flat terrain steep or complex terrain openness highlights topographic radius [m] 10 10 10 convexities, e.g. ridges between hollow ways and rims of bomb craters, negative number of directions 16 16 16 openness emphasizes the lowest parts of recommended histogram linear stretch, linear stretch, linear stretch, concavities, e.g. the actual hollow ways, stretch for positive openness minimum 65°, minimum 85°, minimum 55°, maximum 95° maximum 91° maximum 95° the lowest parts of gorges and the lower edges of cliffs. For consistent readability, recommended histogram linear stretch, linear stretch, linear stretch, stretch for negative openness* minimum 60°, minimum 75°, minimum 45°, it is recommended that negative openness maximum 95° maximum 95° maximum 95° is displayed with inverted greyscale *Inverted grayscale (white to black) presentation works best. (i.e. darker for higher values). Thereby, concave features are always presented by dark tones. 25 2|7 Local dominance Local dominance visualisation of a DEM is based on computing, for every pixel of the DEM, how dominant an observer standing on that point would be for a local surrounding area (Hesse 2016). Dominance as used here is the average steepness of the angle at which the observer looks down at the surrounding land surface. It is higher for points on local elevations as well as on slopes and lower for points in local depressions. Local dominance is computed for pixels within a specified maximum radius and a specified observer height above the Figure 11: Local dominance image (10-20 m search radius) of former field boundaries, roads, and ridge and furrow near Hügelsheim, Germany. Ridge and furrow is only preserved in areas surface. To reduce the noisy appearance that are today covered by forest. Some other features on the image include bomb craters, earth of the resulting image due to small-scale covered bunkers, and trenches. 1 m resolution lidar data © LGL in Baden-Württemberg. surface roughness, a minimum radius can be specified. Pixel brightness is Table 9: Typical settings for calculation and visualization of local dominance. derived from the local dominance values visualization parameter general very flat terrain steep or complex terrain by applying an appropriate greyscale search radius [m] 10-20 10-20 10 histogram stretch. observer height [m] 1.7 1.7 16 This visualisation is well suited for very recommended histogram linear stretch, linear stretch, linear stretch, stretch minimum 0.5, minimum 0.5, minimum 55°, subtle positive relief features such as maximum 1.8 maximum 3.0 maximum 95° former field boundaries or strongly eroded burial mounds, but also delivers very good results for topographic depressions such as dolines, mining traces, or hollow ways. 2|8 Cumulative visibility A viewshed is the area visible from a given vantage point. The viewshed area depends on the topographic position of the vantage point and the surrounding topography, but also on the height of the observer standing on the vantage point, the height of the objects that should be visible to the observer, and the radius under consideration. Figure 12: Cumulative visibility image of burial mounds (high values) and dolines (low values) on Cumulative visibility, on the other hand, the Swabian Alb. 1 m resolution lidar data © LGL in Baden-Württemberg. specifies the size of the area from which Table 10: Typical settings for calculation and visualization of cumulative visibility. a point in the DEM (or an object on that visualization parameter general very flat terrain steep or complex terrain point) is visible to observers of a certain height. DEM visualisation by cumulative radius [m] 1-100 1-100 1-100 visibility is based on computing, for each angular resolution [°] 10 10 10 pixel of the DEM, the size of the area observer height [m] 1.7 1.7 1.7 within a specified radius from which an target height [m] 0.0 0.0 0.0 object is visible (Hesse 2016). Because the recommended histogram linear stretch, linear stretch, linear stretch, surrounding topography plays a dominant stretch minimum 15, minimum 0, minimum 10, role for intervisibility, the resulting raster maximum 55 maximum 25 maximum 65 map can also be a suitable technique to 26 2. Description of techniques visualise that topography. Besides this, such a visualisation can be used as a tool for analysing for example the locations of archaeological sites. The resulting raster map contains percentage values (0…100) of the size of the cumulative visibility area relative to the entire area within the specified radius. 2|9 Accessibility DEM data can be visualised by computing surface accessibility. This means that an algorithm determines, for every pixel of the DEM, the maximum radius of a sphere Figure 13: Accessibility image (20 m maximum radius) of very narrow agricultural terraces with that could be placed on the surface at this vineyards near Jeruzalem, Slovenia. 0.5 m resolution lidar data © ARSO, Slovenia. position without being impeded by the heights of surrounding pixels (Miller 1994). Table 11: Typical settings for calculation and visualization of accessibility. To reduce computation time, the algorithm visualization parameter general only takes into account surrounding pixels within a pre-defined radius. Computation radius [m] 10 time can be further reduced by taking into number of directions 16 account only surrounding pixels along a recommended histogram stretch linear stretch, minimum 0, maximum = radius small number of radial lines rather than all pixels. The range of values in the resulting accessibility raster map corresponds to the range of sphere radii. Greyscale or colour mapping is used to display the results as an image. Accessibility can be used for visualising negative relief features (e.g. pits, hollow ways) and features on slopes (e.g. agricultural terraces). Subtle relief features on more or less horizontal surfaces show up only poorly or not at all. 2|10 Multi-scale integral invariants (MSII) Multi-scale integral invariants is a Figure 14: Multi-scale integral invariants image of a late Roman camp at St. Helena, west of visualisation technique which has Kobarid, Slovenia. 0.5 m resolution lidar data © Walks of Peace in the Soča river Foundation. previously been applied to enhance readability and approach automatic Table 12: Typical settings for calculation and visualization of multi-scale integral invariants. interpretation of cuneiform tablets (Mara et al. 2010), but which is equally valuable visualization parameter general very flat terrain steep or complex terrain for the visualisation of elevation models. number of scales 8 8 8 To compute MSII, n spheres of different minimum radius [m] 1 1 1 diameters are centred on each pixel in maximum radius [m] 11 11 11 the DEM. The percentage of the volume of these spheres, which lies above/below the recommended histogram linear stretch, linear stretch, linear stretch, stretch minimum 1.2, minimum 1.3, minimum 1.2, DEM surface is computed. The result is maximum 1.8 maximum 1.5 maximum 2.5 2. Description of techniques 27 a set of n values (volume fractions above the DEM surface) for each DEM pixel. These sets of n values are interpreted as n-dimensional vectors. By computing the distance of these n-dimensional vectors from a reference vector, the data can be reduced to a raster map containing a single value for each pixel. Low values (low vector distance) indicate high similarity with the reference vector and vice versa. Using an appropriate greyscale histogram stretch, this raster map can be displayed as an image. The reference vector can, for example, be determined by extracting the vector values for a specific relief feature or a point within a cuneiform character or simply by choosing the origin of the n-dimensional coordinate system (i.e. zero). Figure 15: Image of a late Roman camp at St. Helena, west of Kobarid, Slovenia as evident on a 0.5 m resolution lidar data filtered by a Laplacian-of-Gaussian convolution filter with a radius of MSII is almost equally suitable for very 3 pixels (i.e. 1.5 by 1.5 m). diverse terrains from plains to mountains. Because it is a multi-scale approach, it Table 13: Typical settings for calculation and visualization of Laplacian-of-Gaussian. is able to visualize relief features within visualization parameter general a wide range of sizes, i.e. it can clearly laplacian filter radius [px] 3 show very small features while at the low pass filter radius [px] 25 same time preserving at least some impression of the landscape forms. It recommended histogram stretch* linear stretch, minimum -0.05, maximum 0.05 can be quite susceptible to data noise. *Inverted grayscale (white to black) presentation works best. This can be avoided by setting a larger minimum radius; however, this in turn partly compromises the depiction of small and interpretation processes. For example, such algorithms in particular when the details. if the interpreter knows the original goal is visual interpretation: the artefacts scanning density, the method of point are very recognizable and thus help to cloud filtering, and the method of digital avoid misinterpretation of ‘pretty’ surfaces 2|11 elevation model generation, they can that are based on insufficient data. Laplacian-of- judge and take decisions about the various Gaussian (LoG) artefacts that can be found in the data. If point density maps, data gap masks, or actual point cloud data are unavailable, The discrete Laplacian filter computes The occurrence of different artefacts forest masks or topographical maps can the second derivative of elevation, i.e. or unnaturally smooth terrain can be be used as a rough guide. In absence the change of slope. It is a measure of effectively predicted by mapping the of these, very simplified vegetation convexity and can therefore be valuable for density of laser points representing the density maps can be approximated with the visualisation of edges. The Laplacian ground and density of the vegetation vegetation indices (e.g. Normalised filter is often applied to an image that canopy above them. It can be seen from Difference Vegetation Index) derived from has first been smoothed with something Figure 16C that the density of ground freely available high-resolution satellite approximating a Gaussian smoothing points is in places insufficient (data gaps data (Landsat-8 or Sentinel-2, available filter in order to reduce its sensitivity to in red) to accurately map the ground under at earthexplorer.usgs.gov and scihub. noise and this combination is known as forest at 0.5 m resolution, despite the copernicus.eu respectively). Laplacian-of-Gaussian (LoG). high scanning density. In such cases the rasterization algorithm has to interpolate There are several types of artefacts that (infer from neighbouring points) the are commonly found in lidar-based DEMs. 2|12 representation of the ground. The It is useful to know that ‘fish scales’ Visualizing appearance of these interpolated surfaces sometimes found in forest datasets uncertainty varies from algorithm to algorithm, but the (Figure 18A), are a result of a direct point most important is influence on the general cloud rasterization (i.e. without help of a Information about how raw data has been smoothness and preservation (clarity) of Triangulated Irregular Network – TIN), acquired and processed, and about the edges (Figure 17). The artefacts produced and that wave-like features (resembling, method and settings for its presentation, by some interpolation algorithms may be at first glance, tightly spaced ridge has a great impact on the feature detection seen as unattractive and inappropriate for and furrow) are a consequence of poor display. However, we recommend to use registration of scan lines (Figure 19A). 28 2. Description of techniques Figure 16: Lidar data density maps of Castle Montfaucon east of Besancon, France. Density of all points (A), vegetation points (B) and ground points (C) calculated per every pixel (0.5 m) but plotted per m2. Density of ground points calculated and plotted per m2 (D). Figure 17: Anisotropic sky-view factor image of elevation models produced with different rasterization (interpolation) algorithms. Note the difference in interpolation of the castle area and the south-eastern slope leading to it. Natural Neighbours (NN) (A) generates a very smooth terrain, Inverse Distance Weighted (IDW) (B) introduces steps, TIN with Repetitive Interpolation (REIN) (C) introduces triangles, and splines (D) create a range of artefacts. Black stars, sometimes seen in a sky-view directions. Additional artefacts are shown factor image, can be linked to dedicated in Figures 18C, 19B and 19C. Artifacts and processing, where the point-cloud filtering anomalies are not only aesthetic but also process has been optimized to leave the affect the accuracy of analyses based on archaeology as intact as possible. Eight or the terrain surface. sixteen-pointed black stars may be formed around very small (in area) ‘bumps’ that are the unfiltered remains of conifer trees (Figure 18B) or other uprights (Figure 42). This occurs where conifers are too dense for a laser pulse to reach the ground and filtering is set so as not to over smooth the derived elevation model. Stars are generated by the fact that sky-view factor is usually calculated in eight or sixteen 2. Description of techniques 29 Figure 18: (A) ‘Fish scales’ are a result of direct point cloud Figure 19: Wave-like features are a consequence of poor registration of rasterization. 0.5 m spatial resolution lidar data of an area around scan lines (A). Charcoal burning platforms can be seen. 0.5 m spatial Besancon, France, © University of Franche-Comte. resolution lidar data of an area around Besancon, France, © University (B) Black stars typically form around small protruding features on a of Franche-Comte. Severe strip misalignment and over-ambitious sky-view factor (SVF) image when calculating it with 8 or 16 directions. resolution setting result in scrub-like artefacts and fictional steps 0.5 m spatial resolution lidar data of World War I trenches near Kobarid in terrain (B). 1 m spatial resolution lidar data of an area close to in Slovenia. © Walks of Peace in the Soča river Foundation. Tauberbischofsheim © LGL in Baden-Württemberg. Poor raw data (C) 1 m spatial resolution lidar data full of false depressions (C). They preprocessing resulted in a whole range of artefacts imprinted in are a result of poor data processing that classified many below-ground terrain NW of Volarje, Slovenia (C). 1 m spatial resolution lidar data lidar points as true ground. Poštela burial mounds can be seen in the © ARSO, Slovenia. All images display SVF calculated in 8 directions upper left corner of (C). © ARSO, Slovenia. All images display SVF with a 10 m search radius. calculated in 8 directions with a 10 m search radius. 30 3 Guidance for selection of techniques 31 32 3. Guidance f 3. Guidanc or sel e f ection of t or sel echniques ection of t As the complexity of research quiestions increases so does the need for visualizations to convey not only what is there in screaming colour, but 3 also the correct shape, size, relative elevation, degree of preservation, and the context of the immediate environment. There are usually several important To display features more clearly – i.e. with Preparing the types of questions you want answered 3|1 a higher contrast – a histogram stretch images for detection when using lidar data for archaeological is usually required. A histogram stretch prospection. Is there anything out there? Is and interpretation transforms a narrow range of input values there anything beside it? What is it? Why is (e.g. elevation values from 523.2 to 542.7 m it there? Which one is older/younger? The for a hypothetical elevation differentiation) complexity is getting higher with each one As mentioned in the introduction, DEMs onto the whole range of output values of these questions and so is the knowledge contain numerical elevation data and have (e.g. from 0=black to 255=white). There you need to choose the appropriate to be transformed into human-readable are many types of histogram stretch, visualization. As the complexity increases images for visual interpretation. Likewise, with the most frequently used being a so does the need for visualizations to most of the visualisation techniques linear minimum-maximum stretch (with convey not only what is there in screaming described in the previous chapters or without clipping of values below and/ colour, but also the correct shape, size, produce raster maps containing numerical or above a certain value or a certain relative elevation, degree of preservation, values. To display these raster maps as percentage), histogram equalization, and the context of the immediate images, greyscale or colour cast as well as standard deviation stretch, and custom environment. contrast stretch have to be applied. Such a histogram stretch, based on a user defined greyscale cast is for example achieved by frequency curve. Only the linear stretches There are two fundamental questions assigning the colour black (pixel value of preserve the relative distribution of when creating visualizations. How many zero) to the lowest and white (pixel value values between the displayed minimum do you need to accomplish your task, of 255) to the highest numerical value and maximum, while the others do not. and which ones are the best suited to do found in the computed raster file. However, The effect of non-linear stretch is that it as quickly and accurately as possible? in many cases extreme values occur (in we cannot precisely compare different Primarily, the choice of visualization particular along the edges of DEMs) which features or areas. For visual interpretation, depends on the characteristics of the can result in a very low contrast images. however, they can nevertheless be useful sought-after relief features (e.g. size, Therefore, extreme low and high values as they can ‘pull apart’ very similar values shape, convexity/concavity) and the overall usually have to be cut-off/clipped (i.e. and thereby enable a visual differentiation. landscape forms (e.g. smooth, rolling, saturated). rugged). Further factors that can influence Some software allows the user to the choice of visualization can range from While the use of a greyscale cast is automatically adjust histogram stretch computation time to personal preferences. generally the best choice, a colour cast in real time, based only on the values is preferred in some cases, in particular currently displayed on-screen. This means when a clear distinction between positive that the displayed range of values changes and negative relief features is required. when we move around the landscape, For trend removal and local relief model, which is not very useful for comparison the use of blue tones for negative and red/ of features, but works very well for their yellow tones for positive values has been detection. found to be useful. 3. Guidance f 3. Guidanc or sel e f ection of t or sel echniques ection of t 33 Figure 20: An original image where only a minimum-maximum histogram stretch has been applied – a whole range of values is displayed with relative differences between values preserved (A). Images displayed with different techniques of histogram stretch: histogram equalisation (B), standard deviation 2.5 σ (C), 0.5 % saturation of minimum and maximum (D). A histogram stretch with saturation of representation of landscapes or features the red, green and blue components of the minimum for a SVF or a positive openness and/or to reduce the number of images resulting colour image (Figure 3A and C). image renders narrow valley floors, very to be analysed. Such combinations are steep slopes, and the vicinity of high usually applied to greyscale or colour Combining results of various visualizations objects (such as buildings) black, which images rather than to the numerical in a meaningful and deliberate way builds can be useful when delineating concave output of the visualisation algorithms. upon their strengths. Some examples are features on flat or undulating terrain. If Computing the greyscale or colour presented in the Case studies chapters. this is undesirable because features of average of two or more images results in We can not only play with different interest can be ‘hidden’ in such areas, a a translucent overlay of one image over histogram stretches and colour tables, but standard deviation stretch or a standard the other; a weighted average can be used also with the degree of transparency and minimum-maximum stretch provide to adjust the translucency. Alternatively, the mathematical operations to combine a visualization with no or minimal images can be multiplied with each other the layers. Play was deliberately left out saturation, but less contrast. (and the range of values subsequently of brackets, as the process of creating an re-scaled to the 0 to 255 range). Three expressive visualization is indeed as much Two or more visualisations can be greyscale images can be combined art as it is work. combined into a single image. This can be into a RGB composite image where the a useful approach to enhance the visual different input images are displayed as 34 3. Guidance for selection of techniques 3|2 low contrast in areas facing towards soon as the overall topography deviates Reading and (homogeneously bright) or away from from nearly horizontal. Trend removal, presenting the (homogeneously dark) the light source, as local relief model, and local dominance images well as optical illusions (inverted relief). are interchangeable to a certain extent. While very low illumination elevation The advantage of simple trend removal 3|2|1 Choosing angles (< 10°) can and should be used to over local relief model is that it is a much highlight low relief features in areas of low simpler and faster algorithm. On the other appropriate slopes and flat terrain, higher illumination hand, local relief model produces more visualisations elevation angles (> 35°) are required in realistic relative elevation values of relief steeper topography. To investigate features anomalies. Local dominance retains a The choice of a visualisation technique on moderate to steep slopes, shaded relief (limited) visual impression of the overall depends on a number of factors; the should be used with (almost) vertical landscape forms as it produces higher most important are the topography of the illumination to minimize saturated bright/ values on slopes than on horizontal landscape, and morphology and size of dark areas on slopes facing towards/ planes. This entails that on steeper slopes, the sought after features. This is because away from the illumination. In such local dominance will require a different different techniques can preferentially cases, shaded relief images become histogram stretch than in areas with emphasize small-scale depressions similar to Slope images, which can be a gentle topography. or elevations, low relief features on useful alternative in moderate to steep horizontal or sloping planes, or structures topography. Laplacian-of-Gaussian is a very useful on slopes. The type of visualisation often technique to highlight edges and can be also depends on the task at hand and In areas of moderate to steep topography, used on its own for this purpose. It works the current stage in an overall workflow sky-view factor works best to highlight well as an overlay over other visualisations (e.g. general overview, feature mapping, surface depressions and features on to strongly enhance edges and thereby to analysis of details). In almost all cases, slopes. Depending on the range of slopes emphasize relief features. When combined a single visualisation technique will not in a given area under study, different in a (weighted) greyscale average with be sufficient to extract the full amount of histogram stretches may be necessary local dominance or sky-view factor, it information from the data. This entails to avoid bright saturation in gentle helps to overcome saturation on steeper that for a given area, several techniques topography and dark saturation on steep slopes that arises when histogram stretch have to be applied. In the following, an slopes. In areas with flat or very gentle for local dominance or sky-view factor is attempt is made to reduce the complexity topography, sky-view factor is generally adjusted so as to be suitable for gentle to of multiple visualisation techniques and to limited to the presentation of negative moderate topography. provide guidance regarding their selection. relief features (pits, ditches, quarries, erosion areas, dolines…) and becomes very Positive and negative openness are very It is recommended to always begin by sensitive to DEM noise. A good general useful to highlight positive and negative looking at shaded relief overview images rule is to use a histogram stretch of 0.65 relief features, respectively. As openness of the area under investigation, because to 1.0 for diverse terrain and 0.9 to 1.0 for removes the visual impression of overall they provide the most ‘natural’ visual very flat terrain. landscape forms, it is not affected by appearance of the topography and can saturation due to gentle or steep slopes thus help decide which other techniques In areas with flat or very gentle to and may be used in a varied topography. could work well (Table 14). moderate topography, local dominance, Because of the ability to differentially trend removal, and local relief model highlight positive and negative relief In areas of gentle to moderately steep are very helpful to highlight very low features, it is particularly suitable for topography, shaded relief can be relief features such as former field targeted detection of these features. successfully used to investigate relief boundaries or levelled burial mounds. details. However, care has to be taken to On very flat horizontal planes (coastal Like openness, multi-scale integral apply several illumination directions and plains or broad river valleys), elevation invariants may be used in varied to avoid the drawbacks of shaded relief differentiation (greyscale or colour coding topography because it is relatively little such as poor representation of linear of the DEM) can be a very simple and affected by saturation due to overall features parallel to illumination azimuth, effective alternative; however, it fails as topography. Because it is a multi-scale Table 14: Matrix for the suitability of visualisation techniques for selected archaeological relief features in different topographic settings. Begin with a shaded relief overview image then try/add techniques from left to right flat terrain Shaded relief Trend removal / LRM Local dominance Openness or MSII (sun elevation < 10°) (filter radius ~ 20 m) (radius 10-20 m) (radius 10 m) gentle slopes Shaded relief Sky-view factor Trend removal / LRM Local dominance Openness or MSII (sun elevation ~ 30°) (radius ~ 10 m) (filter radius ~ 20 m) (radius 10-20 m) (radius 10 m) moderate slopes Shaded relief SVF (& LoG) Trend removal / LRM LD (&LoG) Openness or MSII (sun elevation ~ 45°) (radius ~ 10 m) (filter radius ~ 20 m) (radius 10-20 m) (radius 10 m) complex Shaded relief SVF (& LoG) LD (&LoG) Openness or MSII topography (sun elevation > 45°) (radius ~ 10 m) (radius 10-20 m) (radius 10 m) 3. Guidance for selection of techniques 35 approach it results in good representation Table 15: Suitability of visualisation techniques for representing selected archaeological of features within a range of scales; topographical features. however, very wide ranges of scale are mining former field burial terraces hollow ridge and computationally intensive and can result in pits boundaries mounds ways furrow reduced contrast. Further visualisation techniques may shaded relief - - + o o - be employed on a case-by-case basis; slope - o o + + ++ however, most if not all tasks can be principal components analysis - - + o + ++ trend removal and LRM ++ + ++ + + ++ performed with the selection described sky-view factor ++ + o ++ ++ ++ above. openness ++ + + + ++ ++ local dominance ++ ++ ++ + ++ ++ 3|2|2 cumulative visibility - - + o + o accessibility - o - o o - Personal multi-scale integral invariants + + o + + + preferences and Laplacian-of-Gaussian + + ++ + + ++ intercomparability - not suitable; o indistinct; + suitable; ++ very suitable While landscape topography and feature morphology are important factors the visual detection of features (i.e. they to happen before the image can be read by that limit the suitability of any given help to see something is there), others a human observer. Thus, understanding visualisation technique, user preferences are more suitable for the interpretation of the entire process chain is necessary to also play a role in the selection. While these features (i.e. they help recognize and be able to correctly interpret the images some users prefer to look at visualisations interpret what it is). and to be able to purposefully modify/ that retain as much as possible a ‘natural’ adjust the data processing parameters. appearance of the terrain (e.g. shaded Training and experience as well as new Producing varied visualisations from a relief or sky-view factor), others prefer methods may over time change the DEM is only the prelude to the equally visualisations that provide a greater level preferences of a given user. While the important process of reading and of abstraction from terrain to image adaptation and development of personal interpreting the images. Because it (e.g. techniques that reduce or remove preferences can generally be expected to happens in our brain (where, contrary the visual impression of the overall improve the quality of the interpretation to the software we have used up to this landscape topography such as trend and the rate of detection, it also entails step, we do not have much control over removal or openness). Other personal a limited intercomparability of mapping algorithms and parameters), this process preferences may relate to the display as results between different users and even of reading and interpreting is commonly either greyscale or colour images, the between areas mapped by the same taken for granted. However, despite the combination of several visualisations (or user at different times of their skills limited control, which we have over our several variants of one visualisation with development. As intercomparability of own brain, it is important to at least different parameter settings) as RGB results is very important in scientific work, develop an awareness of how we see. composite or the type and strength of efforts should be made to improve or at histogram stretch that is applied. When least quantify it. Ideally, repeated mapping Perhaps counterintuitively, perception is producing colour images, the relatively by different persons or paired mapping an active rather than a passive process. high prevalence of various degrees would have the potential to greatly reduce The high-resolution fovea in the human of colour vision deficiency should be intercomparability issues. However, these eye covers only a very small portion of considered (e.g. in particular transitions options are usually prohibited by workload the field of view. The fact necessitates from green to brown or from blue to purple and lack of staff. eye movements, which in turn require a may be not readable for many people). complex feedback between the eye and Colour transitions should always be brain. These eye movements amount to transitions in both hue and brightness, not scan paths by which a scene or an image Perception in hue alone. 3|2|3 is scanned (Yarbus 1967). Furthermore, Yarbus (1967) noted that the scan paths In all cases, it is very helpful to ensure depend on the question or task related consistency with visualizing the same DEM visualisation is a long process chain. to the images. Because the scan path technique, e.g. inverted greyscale for We start with numerical elevation data, interactions between our brain-eye system negative openness and for Laplacian-of- transform these by diverse visualization and the image are subconscious, they Gaussian to retain dark representation of algorithms into raster maps of numerical entail that high-resolution coverage of a low/concave and bright representation of values, visualize the numerical data with given image is likely incomplete unless high/convex relief elements. Preferences greyscale/colour mapping and contrast/ the observer makes a cognitive effort to for specific visualisation techniques also histogram stretch, and present these look at all parts of the image. This effort depend on the task at hand. While some images on computer displays or print can be facilitated by applying a (visible or techniques are particularly suitable for them out to the final physical image that imaginary) grid to the image. is perceived by the human eye. All this has 36 3. Guidance for selection of techniques Figure 21: Visualization techniques in different topographic settings. (A) Plough headlands on a flat plain near Endingen am Kaiserstuhl. 1 m lidar data © LGL in Baden-Württemberg. (B) Three different types of World War I trenches with shelters on gentle NE slopes of Črni hribi, near Renče, Slovenia. 1 m lidar data © ARSO, Slovenia. (C) Charcoal burning platforms in the hills of the Black forest. 1 m lidar data © LGL in Baden-Württemberg. (D) A late Roman campo on a rocky outcrop with a church of St. Helena, west of Kobarid, Slovenia. 0.5 m lidar data © Walks of Peace in the Soča river Foundation. 3. Guidance for selection of techniques 37 Yet even if we do consciously force Visualization of ourselves to look at every portion of an 3|3 datasets other than image, we will commonly only see those features which we know and recognize and lidar which we are looking for. Adams (1982) has wonderfully described this as a ‘Somebody The described visualization techniques and Else’s Problem field’ (SEP): guidance for their selection are applicable to raster elevation data of different sources “An SEP[...] is something we can’t and resolutions, applied at various scales, see, or don’t see, or our brain and for observation of a range of entities doesn’t let us see, because we think (landscape, site object). Raster elevation that it’s somebody else’s problem models can be derived from aerial laser [...] The brain just edits it out, it’s like scanning (ALS) with a spatial resolution of a blind spot. If you look at it directly 10 cm at best, but more frequently in the you won’t see it unless you know range of 0.5 m to 1 m. ALS data are most precisely what it is. Your only hope frequently used for a site or regional scale is to catch it by surprise out of the projects. Much higher resolution DEMs can corner of your eye [...] It relies on be derived from, for example, terrestrial people’s natural predisposition not laser scanning, Structure-from-Motion to see anything they don’t want to, (SfM) modelling (also referred to as close weren’t expecting, or can’t explain” range photogrammetry), or structured (Adams 1982). light scanning. The resolution of such DEMs varies from a few centimetres to The subconscious pre-interpretation of sub-millimetre and they are best for what we see, e.g. the fact that we tend observation of smaller areas or individual to close gaps or group objects based on objects. similarity or proximity to form a ‘gestalt’ (Boeree 2009), can allow us to recognize Coarser DEMs, for example traditional relief features which are only partially national and world-wide datasets preserved. However, the very same effect produced with surveying, photogrammetry, can also be misleading and can result in or interferometry, have a spatial misinterpretations. The human eye-brain resolution ranging from 5 m to 30 m. system is a remarkably efficient and Free sources of global datasets are SRTM adaptable pattern recognition system, Global and ASTER Global DEM, both but it needs to be trained to detect these with approximately 30 m resolution and patterns. Its efficiency can also be a available at earthexplorer.usgs.gov. These hindrance, because it might see pattern are appropriate for national, continental, where there is none. A possible approach and global environmental studies. to counteracting such limitations may, for example, be paired mapping – i.e. two persons have to agree on what they see. Besides understanding the described visualization algorithms and consciously counteracting the limitations of our visual perception, training and experience in reading DEM visualizations Figure 21 continued. are therefore key for correct and (in times of overwhelming amounts of data and very limited manpower) rapid interpretation. Yet, having looked very purposefully at the same lidar image again and again, you can sometimes still discover something new ‘out of the corner of your eye’. 38 3. Guidance for selection of techniques Figure 22: Different scales of observation. Sky-view factor image of a detail from a richly decorated kerbstone from Newgrange, Ireland, produced with Structure-from-Motion of hand held photographs (A). Scattered bones (B) and looting traces (C) on a desert floor in Peru. The elevation model was computed from images taken with a camera on a pole and is presented here with a combination of local dominance and shaded relief. Grand Canyon (D) as evidenced by 30 m SRTM data presented with a combination of sky-view factor and shaded relief. Central Massif, the Jura and the Alps (E) presented with a combination of anisotropic sky-view factor and shaded relief. 3. Guidance for selection of techniques 39 3|4 What’s in a name? Table 16: Metadata required when presenting DEM visualisations. visualisation technique mandatory parameters ancillary parameters shaded relief illumination azimuth illumination elevation Metadata about data processing has to vertical exaggeration factor follow the lidar dataset to the visualization histogram stretch intended for detection and interpretation slope histogram stretch and on to the final product – a thematic (minimum/maximum slope) map. Lidar Base Specification (Heidemann trend removal and LRM low pass filter radius histogram stretch / colour code 2014) provides specifications to acquire type of low pass filter and procure lidar data and is a good openness positive/negative number of search directions source of information on what metadata greyscale / inverted greyscale histogram stretch search radius should be stored. It specifies that metadata deliverables should include sky-view factor search radius number of search directions histogram stretch (Heidemann 2014: 13): local dominance search radius observer height • collection report detailing mission histogram stretch planning and flight logs, cumulative visibility search radius observer/target height angular resolution • survey report detailing the collection of control and reference points used accessibility search radius number of search directions for calibration and QA/QC including control and calibration points, MSII reference vector (if not zero) number of scales minimum and maximum radius • processing report detailing histogram stretch calibration, classification, and product Laplacian-of-Gaussian filter radius greyscale / inverted greyscale generation procedures, histogram stretch • a QA/QC report, detailing procedures for analysis, accuracy assessment and validation of the point data, bare- earth surface, and other optional deliverables, • georeferenced, digital spatial representation of the precise extents of each delivered dataset, • product metadata files for the overall project, each scanning mission, and each deliverable product group. It also gives a descriptive template and a completed example in Appendix 3, ‘Lidar Metadata Example’ and Appendix 4, ‘Lidar Metadata Template.’ Metadata specifications, however, usually only cover the hardware and software used to process or create the dataset, with additional explanations possibly given in the data quality sections (e.g. inclusion or omissions of features). Processing parameters for filtering and visualization are rarely given. If not explicitly demanded by the financer, few technicians and Figure 23: Is it a burial mound or a hole in the ground, a quarry or a hill? Someone observing the image needs to know the direction of the light source, because it has a definitive impact on the perception of the landscape. 315° azimuth illumination (A) and 135° azimuth (B). Terrain features seem inverted on (B). 40 3. Guidance for selection of techniques scientist input all the required fields, let alone the optional fields into the metadata scheme. Nevertheless, entering the following records can assist the interpretation process at several stages: • data acquisition: lidar sensor make and model, nominal scanning density, nominal swath overlap, date of data collection; • description of post-collection processing: method(s) used, parameter settings, description of the processing goal (e.g. production of a terrain model, removing just the vegetation, production of a surface model), elevation model spatial resolution; • visualization: method(s) used, parameter settings (see below and Table 16); • interpretation process: aims of interpretation (e.g. identifying locations of individual burial mounds and marking their circumference), reliability of the results (qualitative if quantitative evaluation is not possible, e.g. low to high, description of each class is recommended). A minimum metadata requirement for visualization is the applied algorithm and the settings of its principal parameters. Providing all details may often not be practical when preparing images for publications. Figure 24: Overgrown remains of an abandoned village Novi Breg (Naubacher), Slovenia, as represented by local relief model technique with a filter radius of 25 m (A). Remains of houses Table 16 lists the mandatory and ancillary can be seen as bright rectangles. Big black blobs are dolines. It is straightforward to replicate parameters. The necessity of providing the such an image with the minimum data that is provided. However, when reproduction or relative ancillary parameters depends largely on comparison across areas is not required, combinations of visualizations can give images that the purpose of presenting the image: if the are easier to read, and details about visualizations may be omitted (B). 1 m resolution lidar data © ARSO, Slovenia. aim is solely to provide a visual illustration, many parameters may be omitted. main parameters in the name of the output remembering the necessities when However, if the aim is a quantitative file, as well as to store all the settings in a producing images for publications or analysis of specific features, a discussion processing log. trying to replicate the process with other of feature details, or a comparison of data. different visualisations, more parameters From a filename Novi_Breg_DTM_1m_ have to be provided for the reader to LRM_FSc_Ilin_FR10_MR25.tif it is understand the image and to be able to quite obvious that the local relief model reproduce the process of image creation. visualization was applied to a 1 m Sometimes a seemingly trivial parameter, resolution digital terrain model (Figure such as shaded relief illumination 24A). The settings applied were: azimuth, can have a huge impact on the perception of the landscape and has to be • FSc – filter shape: circular known to the observer (Figure 24). • Ilin – interpolation method: linear • FR – filter radius [px]: 10 It is very easy to repeatedly calculate • MR – maximum range [px]: 25 visualizations with varied settings. A good practice is to store the applied Such filenames may seem strange at visualization technique and settings of the a first, but are the best safeguard for 3. Guidance for selection of techniques 41 Figure 25: Hillshading algorithms in ArcMap (ArcMap 2012) and RVT use the same method (Yoëli 1965), but the results are quite different as evident from the image of skyscrapers in Lower Manhattan. The hillshading image calculated with ArcGIS (A) seems to be less clear (as if using a DEM of a lower resolution) than the one calculated with RVT (B). Note the distinctly visible elevation steps on (B), especially the steps on the facade of the lower right skyscraper (B), which are completely missing on (A). Both (A) and (B) have saturated high and low values (completely white and black areas). However, RVT calculates hillshading with a full range of values therefore details in black areas are revealed to a great extent on the original unsaturated image (C). The same is true for the difference in calculation of sky-view factor with RVT and SAGA GIS. All figures are calculated with 35° Sun elevation and 315° azimuth using the same DEM. 1 m spatial resolution lidar data © USGS. 42 Figure 26: Visualization techniques with different relief features. Continues on next page. 43 Figure 26 continued. 44 3. Guidance for selection of techniques Figure 26 continued. 3. Guidance for selection of techniques 45 Figure 26 continued. 46 4 Tools 47 48 4. Tools 4 Two toolboxes with implemented methods that give scientifically sound results, and are documented and supported by research papers. Researchers can only very recently benefit execution of the program also generates Relief Visualization from free software for calculation of 4|1 a processing log file per input file, thus Toolbox (RVT) advanced visualization techniques. Two automatically assembling metadata. such examples are Relief Visualization Toolbox (RVT) and Lidar Visualization Relief Visualization Toolbox (available The toolbox supports elevation raster file Toolbox (LiVT), both free, easy-to-use at http://iaps.zrc-sazu.si/en/rvt) is a data conversion. It is possible to convert all applications to create visualisations from standalone toolbox that does not require frequently used single band raster formats high-resolution digital elevation data external software to run. It provides a into GeoTIFF, ASCII gridded XYZ, Erdas derived from aerial laser scanning (lidar) narrowed range of methods and their Imagine file and ENVI file formats. RVT or other sources, e.g. structure-from- settings are limited to the most important. can thus be used to convert DEM files in motion photogrammetry. The second is The selected techniques have been proven your favourite format to formats supported oriented towards more advanced users effective for detection of small-scale by LiVT, e.g. ENVI, and to convert results with some knowledge in data processing features and default values are set to do of LiVT back to a more widely supported and geographic information systems, this task. Some techniques (e.g. sky-view format, e.g. GeoTIFF. while the first is tailored also for the factor and openness) tile large datasets. very beginners in relief interpretation. RVT can also mosaic tiled raster data Implemented methods give scientifically RVT (v 1.3) can process all GDAL (GDAL and visualizations with selected settings sound results, and are documented and Development Team 2014) supported raster can be executed on a predefined list of supported by research papers. Both formats (e.g. GeoTIFF, generic binary file, files without running the Graphical User toolboxes should work on any ordinary Erdas Imagine file, ENVI file, Arc/Info ASCII Interface (GUI). For more information, office computer. With large files, there Grid, ASCII gridded XYZ, JPEG2000…). It see the RVT Manual (Instructions for use) could be issues with RAM and TIFF size can process multiple files from various available from RVT web page. limits. Should you want to process very folders and all techniques in one go. It large files (e.g. larger than 40 km2 at 1 m is possible to run it without opening the Lidar Visualization resolution), you will most probably have to graphical user interface, on a list of files, tile them. and with predefined settings. 4|2 Toolbox (LiVT) Other raster processing software Visualizations output a pair of GeoTIFFs solutions, such as QGIS, SAGA GIS, GRASS, per each selected visualization. One file Lidar Visualization Toolbox (available at or ArcGIS, have also started to provide gives a precise calculated result, and http://sourceforge.net/projects/livt) can access to at least some of the described the other a simplified result (‘a picture’), compute all of the techniques described visualizations techniques, but are not optimized for viewing in non-GIS software, in Chapter Description of techniques and discussed here, as we have not thoroughly e.g. by Windows Photo Viewer or by has settings for intricate manipulation of tested the implemented methods. Testing Preview for Mac users. All output files the parameters for each algorithm. For is essential because software solutions are written into the folder of the input example, when computing a Local Relief can give diverse results despite using the file. Output file names for visualizations Model options allow setting the filter same method and small differences in are composed of the input file name, shape (circular or square), filter radius, implementation can have a huge impact on and suffixes describing the selected interpolation method (inverse distance, the representation of features (Figure 25). method and processing parameters. Each nearest neighbour, average, linear, or 4. Tools 49 bilinear), and a maximum range. It is limited in terms of file types it supports and can presently only process and outputs generic BIL (band interleaved by line) files and ENVI raster files. LiVT also includes a tool to rasterize ASCII XYZ point clouds (e.g. last return lidar data) with four different interpolation methods. It runs on Windows and requires the Microsoft .NET Framework version 4. It does not include a viewer; additional software (GIS) will therefore be necessary to display the processing results. All algorithms output floating point binary raster maps. For viewing them in GIS, greyscale or colour coding will have to be applied. In addition, appropriate contrast/ histogram stretches may be necessary to produce optimal results. File size limits vary from algorithm to algorithm and range from 30 to 144 million pixels per file. Computation times depend, in addition to the chosen algorithm, also on the selected parameters and of course the processing power of the CPU. LiVT only uses one processor core, but several instances of LiVT can work on separate cores in a multi-processor PC. Simpler algorithms are generally much faster Figure 27: Relief Visualization Toolbox (RVT) offers a range of only the best than more complex ones. Where filter or techniques with just the essential options. It is extremely easy to use and can search radius can be specified, doubling process multiple files and all techniques in one go. the radius will approximately double the computation time for the same input file in a radial algorithm (one that considers pixels along a limited number of radial lines) but will quadruple the computation time in an algorithm which considers all pixels within the radius. Figure 28: Lidar Visualization Toolbox (LiVT) supports calculation of various techniques with intricate options to manipulate the settings. 50 5 Case studies 51 52 5. Case studies The following chapters focus on some typical studies where lidar visualizations have helped archaeologists and geomorphologists better 5 understand their areas of interest. Archaeologists are highly involved in the et al. 2010), coastal dunes (Woolard flowing, the first material to be laid down interpretation of microrelief features and Colby 2002), glacial forms (Smith is usually coarse, with fine sands and silts from lidar data because precise relief et al. 2006)), toward the edges. models allow for a detailed mapping • modelling of processes (e.g. glacial, and measurements of overgrown periglacial (Smith et al. 2006), fluvial Alluvial fans are of practical and economic archaeological structures (dams, (French 2003), and aeolian processes importance to society, particularly in arid ramparts, ditches, pits, quarries, remains (Sankey et al. 2010), soil erosion and and semi-arid areas where they may be of houses, etc.) (Kershaw 2003; Devereux deposition (Chiverrell et al. 2008)), the principal groundwater source for et al. 2005), fossil fields and cultivation and furthermore irrigation farming and the sustenance of terraces (Sittler 2004), former land division • natural hazards risk assessment life. The use of desert fans as permanent (e.g. Roman centuriatio), abandoned (for e.g. landslides (Metternicht et sources of water is limited, however, quarries and mining pits, burial tumuli al. 2005), rock falls (Lan et al. 2010), because periodic rainfall or snowmelt and ancient roads (e.g. Roman, medieval), glacial debris (Conway et al. 2010)). provides only a very slow rate of recharge. and other remains of former cultural landscapes in specific environments Creating a settlement on an alluvial fan where other techniques or field surveying can be dangerous, because they are prone Alluvial fan at Craig do not give satisfactory results (Challis 5|1 to flooding. Rushing water, mud, and canyon, Saline 2006; Challis et al. 2008; Crutchley debris can threaten communities many 2009a; Crutchley 2009b), even in extreme Valley, California, kilometres away from the apex of the conditions, such as mapping of features USA alluvial fan. Lidar data reveals details of under a dense canopy of a tropical forest the intricate drainage system formed on (Fernandez-Diaz et al. 2014). The alluvial fan that stretches out of the fan and is often used for hydrological Craig canyon at the edge of the Inyo studies of groundwater sources and flood Aerial laser scanning data is also used Mountains onto the floor of the Saline modelling. in geomorphology, either directly in a Valley in California, USA, is a typical form of elevation data for detecting the fan-shaped deposit of unconsolidated surface discontinuities (e.g. breaklines, water-transported material (alluvium). lineaments) and forms (e.g. fans, lava Alluvial fans typically form at the base flows), or indirectly by classification of topographic features where there is of surface features relevant for a marked break in slope. The apex (the geomorphological processes. High- narrow part) of each fan is just within a resolution relief models enable: canyon mouth that serves as the outlet for a mountain drainage system. Sediment • mapping of geomorphological from erosion within the mountains is features (e.g. drainage network and moved by these drainage systems to channel heads (Passalacqua et al. the adjacent basin, forming the wide 2010), fallout deposits, epiclastic triangle – the fan’s apron. Since the rivers sediments and lava flows (Fornaciai that deposit alluvial fans tend to be fast 5. Case studies 53 Figure 29: Craig canyon alluvial fan (California, USA) as evident on 0.5 m spatial resolution lidar data (California, USA). Coloured stretched elevations over SVF and shaded relief provide a detailed view of the intricate drainage system of the fan’s apron. The typical shape of the fan is also easily recognizable. © Plate Boundary Observatory by NCALM, USA. Table 17: Lidar scanning parameters of Craig canyon alluvial fan (California, USA)1. parameter value scanner type Optech GEMINI platform fixed wing date between 18th April and 24th April 2008, exact date unknown average last and only returns per m2 on a 3.6 combined dataset spatial resolution of the final elevation model [m] 0.5 1 Lidar data acquired by Plate Boundary Observatory by NCALM, USA. PBO is operated by UNAVCO for EarthScope and supported by the National Science Foundation (No. EAR-0350028 and EAR-0732947). Figure 30: Alluvial fan’s apron at the foot of Craig Canyon in the Saline Valley, USA. Photo by Brian Lockett, Air-and-Space.com. 54 5. Case studies 5|2 were taken by people traveling along them, impossible to walk on, and have gone Hollow ways at Volčji leading in some places to formation of unexplored for decades. Their changing Potok, Slovenia river-like braided channels, branching and direction and position on slopes make converging (Edgeworth 2011: 109). them difficult to visualize on lidar data. Historical hollow ways or sunken paths A combination of several techniques is are routes that have been eroded down to “Instead of fragments of transport necessary to do this properly (Figure 32). bedrock by the flow of people, animals, networks connecting fixed points and carts, so that they are recessed in the landscape, hollow ways beneath the level of the surrounding are rather ‘messy’ landscapes of landscape. They form in the soft stone, movement, interwoven swarms not on the hard rock. Regular trampling of different scrapes and traces by people or livestock suppresses the of movement, almost biological growth of vegetation on trails and reduces shapes, ‘organic’ entanglement of the water infiltration rate. This results lines that emerges from growth in increased surface runoff along trails, and differentiation through rhythms especially on steep slopes. During the dry of human and animal movement, season, trampling displaces surface soil, change of seasons, water dynamics… providing a source of sediment during They are not about getting the rainy season. The trails become a somewhere, from point A to point conduit for surface runoff and a source B, but about being in the landscape, Table 18: Lidar scanning parameters of Žiški vrh and Volčji Potok, Slovenia. of sediment, resulting in increased living the daily life. People inhabit erosion rates (George et al. 2004). As the the landscape along these paths. parameter value roads deepened, they became natural They are lines along which past scanner type Riegl LMS-Q780 waterways. Rain drained into and down landscapes were created, in a platform helicopter them; storms turned them into temporary messy way, from things and features date July 2014 to rivers, sluicing away the loose rock debris encountered along the path.”(Mlekuž March 2015 and cutting the roads still further below 2013: 41). average last and only 3.0 the meadows and the fields.Water erosion returns per m2 on a speeded the hollowing-out process and Most have through time become combined dataset made some lanes muddy and impassable. unpassable, so overgrown by nettles, spatial resolution of the 0.5 When this happened, alternative routes brambles, and scrub that they are final elevation model [m] Figure 31: Parallel hollow ways in a forest NW of Volčji Potok, Slovenia. Photo by Žiga Kokalj. 5. Case studies 55 Figure 32: Hollow ways and headstreams below a hillfort on Žiški vrh, northern Slovenia, as evident on a combination of sky-view factor (0.65-1.0, 25 % opacity, multiply), positive openness (70-93, 50 % opacity, overlay), slope (0-45°, 50 % opacity, luminosity) and hillshading. 0.5 m spatial resolution lidar data © ARSO, Slovenia. Figure 33: Hollow ways developed on a narrow passageway along the river. Volčji Potok, northern Slovenia. 0.5 m spatial resolution lidar data © ARSO, Slovenia. 56 5. Case studies indicate loyalty of these communities to but much of the data point cloud contains Prehistoric and the Roman state (Laharnar 2012: 243– double entries and negative outlier points 5|3 Roman settlement 244). that are classified as ground (erroneous above Knežak, points lying below true ground). The Slovenia Countries, states, and regions increasingly data was also processed with a specific offer their lidar datasets for free and open requirement to classify buildings. Experts Gradišče nad Knežakom hillfort from the use. There is thus a growing opportunity interested in true representation of Iron and Roman Age was once a central for development of products and services the ground and small-scale features settlement in the Pivka area, Slovenia. based on high-resolution elevation specifically, can gain much by even just It sits on the top of a ridge and is today models. However, data users still need reprocessing the data with different totally covered by impenetrable forest. Its at least the basic knowledge about how software (see Figures 35 and 36). shape is oval-triangular, demarcated by these data were acquired, processed, and well-preserved ramparts on the western, can be developed or improved according to southern, and eastern side, and steep their requirements. The Slovenian national slopes on its northern edge. On the lidar dataset is of a generally high quality, north-western boundary the settlement is fortified by a tower. The plateau is slowly Table 19: Lidar scanning parameters of Gradišče nad Knežakom, Slovenia. reducing in height towards the south in low parameter value terraces. Small depressions on them are most probably traces of partially-buried scanner type Riegl LMS-Q780 structures (Laharnar 2012: 66–67). platform helicopter date between February 11th and November 21st Absence of Roman military finds and 2014, exact date unknown traces of military engagement from average last and only returns per m2 on a 4.2 the Late Republican period, and the combined dataset continuous occupation of old settlements spatial resolution of the final elevation model [m] 0.5 Figure 34: Gradišče nad Knežakom hillfort in the centre of the image is today totally covered by impenetrable forest. Photo by Boštjan Laharnar. 5. Case studies 57 Figure 35: Iron Age settlement above Knežak and its surroundings as evident on the classified Slovenian national lidar data (GKOT, © ARSO, Slovenia) which has been processed by gLidar software, rasterized to a 0.5 m grid. Many small bumps and non-existent holes are visible (see also Figure 18C). The image is a combination of sky-view factor (0.7-1.0, 70 % opacity, multiply), slope (0-55°, 100 % opacity, overlay) and hillshading. Figure 36: The same are as in Figure 35, but processed with different software (LASTools). Even with minimal effort (default settings), the resulting image gives a clearer picture of the archaeological and natural topography, despite losing some of the ‘sharpness’. Clearance cairns and field boundaries can be seen on the lower part of the image. 0.5 m spatial resolution lidar data © ARSO, Slovenia. 58 5. Case studies 5|4 are positioned deliberately, so as Based on this study they have found Barrow cemetery in to change the visual structure of that the landscape around Poštela had Pivola, Slovenia the landscape, to dominate the a critical role as a medium, because it foreground or short distance view, became a cultural landscape, a polygon The burial site at Pivola, on the Drava plain being immediate, close and engaging for expressing ideas and messages. at the foothills of the Pohorje Mountains to all senses. When inside this Respecting, relating to, or changing the in Slovenia, contains dozens of mounds. It group, a viewer would find himself existing spatial order was a powerful is part of a large area where dwellers of in a well-bounded visual envelope message, which reproduced or subverted the fortified hillfort at Poštela were buried and dominated by an immediate the existing social identities and power in the Early Iron Age (8th to 6th century presence of barrows (Figure 40). relations (Mlekuž and Črešnar 2014: 209). BC). As a consequence of ideological They are less striking in the middle and religious changes, burial in mounds distance, but still represent an became a new practice at the beginning of important compositional element of the Early Iron Age, with warriors becoming the landscape” (Mlekuž and Črešnar the most prominent social stratum (Gulič 2014: 205). and Črešnar 2012). In their study based on lidar data Mlekuž and Črešnar (2014) Table 20: Lidar scanning parameters of Pivola burial mounds, Slovenia. have investigated the prominence and parameter value deliberateness of positioning different groups of burial mounds around Poštela scanner type Riegl LMS-Q780 hillfort into the landscape. platform helicopter date between 12th March and 20th October 2014, “The Pivola group is situated in exact date unknown a compact visual envelope of the average last and only returns per m2 on a 4.7 valley, in a shallow depression combined dataset delimited by natural low ridges to spatial resolution of the final elevation model [m] 0.5 the north and south. The barrows Figure 37: One of the burial mounds at Pivola. Except in late winter, dense brambles make it very difficult to move around. Photo by Žiga Kokalj. 5. Case studies 59 Figure 38: Burial mounds at Pivola as displayed by a Local Relief Model transparently overlaid on a sky-view factor image and shaded relief. The state of preservation of the mounds differs significantly. The ones in the forest are well preserved (1), mounds in the Botanical Garden were damaged when the location served as a military base (2), and the mounds in the fields have been nearly completely levelled (3). 0.5 m spatial resolution lidar data © ARSO, Slovenia. Figure 39: Visual structure of landscape as represented by Higuchi’s foreground range total viewshed around Pivola barrow group. 1 m spatial resolution lidar data © ARSO, Slovenia. Reproduced with permission by D. Mlekuž. 60 5. Case studies 5|5 The general surviving layout of the anti- model (nDSM), which represents heights WWII anti-glider glider defences can now be appreciated of features relative to the ground. Features defences at Culbin, – particularly the ranked lines intended to with a height range between 0.2 m to 4 m Scotland dissuade gliders from landing along the were selected, selecting (by hand) those beach. The points depicting the location features likely to be posts (i.e. discrete Peter Crow and Matt Ritchie; Forest of the posts have not been checked on the points away from obvious patches of tall Research and Forestry Commission ground – they represent features detected vegetation). Posts shorter than 0.2 m were Scotland by lidar standing 0.2 m or higher within not plotted as it was difficult to distinguish The tidal lagoon at Culbin, on the coast a flat tidal environment. The data was them from any vegetation present of the Moray Firth in Scotland (NH 969 manipulated by cropping the lidar models (Forestry Commission Scotland 2012: 625), is dotted with regular lines of poles to remove the surrounding forest before 34–36). Sky-view factor is particularly that were built in 1940 in order to deter subtracting the digital terrain model (DTM) useful for visualizing the posts and other enemy gliders from landing. The poles from the digital surface model (DSM) ‘pointed’ features because they appear as were erected by digging a hole, inserting to produce a normalized digital surface distinct black stars on an image. an old herring barrel, dropping in the pole and securing it using stones and sand. Hundreds of poles remain standing within Table 21: Lidar scanning parameters of anti-glider defences at Culbin, Scotland. a tidal lagoon, protected by sand dunes and forest. However, many have been parameter value reduced to little more than stumps. In an scanner type Optec ALTM 3033 attempt to record those remaining before platform fixed wing they are lost to the elements, a lidar survey date 7th October 2009, low tide of the coastal area was commissioned by Forestry Commission Scotland, with average last and only returns per m2 on a 4, with a 0.9 m footprint combined dataset survey parameters specifically designed to ensure that each upstanding post would be spatial resolution of the final elevation model [m] 0.5 m recorded. Figure 40: Lines of poles at Culbin Sands. Photo © Forestry Commission Scotland. 5. Case studies 61 Figure 41: Changes in elevation with added intensity information, Culbin Sands. Lidar intensity is usually recorded at the time of scanning and can provide additional discriminating information about the landscape. However, poles are too narrow to be readily recognised on the elevation image. 0.5 m spatial resolution lidar data © Forestry Commission Scotland. Figure 42: A close-up view of Culbin Sands as seen on a sky-view factor image calculated in 8 directions. Poles are easily identifiable as black stars. 0.5 m spatial resolution lidar data © Forestry Commission Scotland. 62 5. Case studies 5|6 Valley and the Triglav Lakes Valley. The The ridge between mountains Planja Geological features area is important especially because (1964 m) and Bavški Grintavec (2333 m) in the Julian Alps, of its numerous corrosion and glacial testifies to the powerful tectonic processes Slovenia erosion forms at various developmental in the Cenozoic that have shifted the stages. High-mountain dipping limestone limestone strata almost vertically. They pavements (podi) are areas of wrinkled, now lie almost 2 km above the valley The Julian Alps are the biggest and glacially transformed karstic plateaus floor and the underlying rock is still not highest Alpine mountain chain in Slovenia, with a rocky surface, usually above the exposed (Planina 1954: 192). The location rising to 2864 m at Mount Triglav. They tree line. They are full of other intricate between Vrh Brda and Vrh Rut is known are mostly built by 200-million-year-old surface structures, such as clints, grikes, as Ribežni (‘slicers’) (Figure 43) because of layers of limestone and dolomite that and karren tables, but also subterranean the exposed strata, which are beautifully are up to 2 km deep. The landscape is features like vertical caves. Karren tables shown on an appropriate lidar visualization characterized by an extraordinary diversity stand out among the glacial karstic shapes (Figure 45). of geomorphology, especially glacial, in Velika vrata because the area is their fluvial, and karstic forms of outstanding locus typicus in Slovenia. beauty that preserves a cultural backdrop. Endemic, rare, and threatened plant and animal species are protected in the Triglav Table 22: Lidar scanning parameters of Velika vrata and Vrh Brda, Slovenia. national park that covers most of the area. Velika vrata Vrh Brda scanner type Riegl LMS-Q780 Riegl LMS-Q780 Lidar has become an indispensable tool in geomorphological and geological platform helicopter helicopter studies that traditionally rely heavily on date July 2014 to January 2015 July 2014 to January 2015 visual interpretation of the landscape. One average last and only returns per m2 4.7 5.0 example of geomorphological beauty is on a combined dataset Velika vrata (big door) area (Figure 44), a spatial resolution of the final elevation 0.5 1.0 high mountain pass between the Trenta model [m] Figure 43: A view across the ridge towards Vrh Brda and Bavški Grintavec in the background. Photo by Jovan Cukut. 5. Case studies 63 Figure 44: Glacially transformed karstic plateau at Velika vrata, Slovenia. Flat areas of clints are presented in light colours while grikes, crevices, and steep slopes are dark. The image is a combination of sky-view factor (0.5-1.0, 70 % opacity, multiply), slope (0-65°, 80 % opacity, overlay) and RGB image of hillshading from three directions (R: 315°, G: 15°, B: 75°), all computed on a digital surface model. 0.5 m spatial resolution lidar data © ARSO, Slovenia. Figure 45: Dachstein limestone strata and scree fields exposed at Vrh Brda (1952 m). The top of the ridge is highlighted in white, while colours represent different orientation of the slopes. The visualization is a combination of sky-view factor (0.55-1.0, 30 % opacity, multiply), positive openness (65-95, 50 % opacity, overlay), slope (0-65°, 80 % opacity, luminosity) and RGB image of hillshading from three directions (R: 315°, G: 15°, B: 75°). 1 m lidar data © ARSO, Slovenia. 64 5. Case studies 5|7 speed up the geomorphic lifecycle of The Granite Dells, precariously balanced rocks while high soil Arizona, USA production rates completely decompose the rocks in the subsurface. This is linked The Granite Dells is a geological feature to the density of joints; high joint densities north of Prescott, Arizona, USA. The Dells create small boulders that completely consist of exposed bedrock and large decompose prior to exposure while low granite boulders that have eroded into joint densities create relatively large and rounded bumpy and unusual shapes giving stable boulders. the rocks a rippled appearance (Figure 46). Dells’ granite formed 1.4 billion years ago at a depth of 2-3 km and has since been eroded away. Weathering along joints produced rounded shapes and other Table 23: Lidar scanning parameters of Granite Dells, Arizona, USA2. unusual rock formations that characterize the Granite Dells. The process is known parameter value as spheroidal weathering and is common scanner type Optech Gemini in granitic terrains. Also characteristic of platform fixed wing the area are precariously balanced rocks (PBRs) that are present in a range of sizes, date November 9th 2009 masses, geometric configurations, and average last and only returns per m2 on a 6.7 precariousness. Haddad and Arrowsmith combined dataset (2011: 145) have found that the emergence spatial resolution of the final elevation model [m] 0.5 of precariously balanced rocks is affected 2 Lidar data acquired by the NCALM at the University of Houston and the University of California, Berkeley, by changes in the rates of soil production on behalf of David E. Haddad (Arizona State University) as part of an NCALM graduate student Seed Grant: Geologic and Geomorphic Characterization of Precariously Balanced Rocks. NCALM is funded by the National and transport; high soil transport rates Science Foundation’s Division for Earth Sciences, Instrumentation, and Facilities Program. Figure 46: Granite Dells at Watson Lake. Photo by Michael Wilson. 5. Case studies 65 Figure 47: Digital terrain model (DTM) of Granite Dells, Arizona, USA. Much of the detail in the landscape is lost despite dedicated processing. A combination of anisotropic sky- view factor (0.2-1.0, 40 % opacity, multiply), slope (0-80°, 60 % opacity, overlay) and RGB image of hillshading from three directions (R: 337.5°, G: 0°, B: 22.5°). 0.5 m spatial resolution lidar data © NCALM, USA. Figure 48: Digital surface model (DSM) of the same area as in Figure 47. Choosing a DSM over a DTM is wise in open landscapes without extensive tree and bush cover, because we can produce a raster image considering a higher number of lidar points, thus getting a ‘sharper’ result. In case of this figure the DSM makes visible trees, houses, and cars, but more importantly also the joints in granite that are characteristic of the Dells. 0.5 m spatial resolution lidar data © NCALM, USA. 66 5. Case studies 5|8 or rough and jagged surfaces (older have found that the spatial and volumetric Cinder cones and flows). Detailed relief models permit distributions of lava reflect the effusion lava flows on Mauna mapping and morphometric analysis rate and interactions with topography. Loa, Hawaii, USA of flow components, such as channels, The main channel serves to transport, surface folds, cracks, blocks, and rather than store, lava and the pre-existing surface roughness, as well as along-flow topography exerts a primary control on the The world’s largest volcano, Mauna Loa on variations in flow type. Differencing of pre- (3D) spatial distribution of individual lava Hawai’i Big Island, has come to epitomize eruptive and post-eruptive DEMs allows flows. shield-stage volcanism. Its volume is analysis of flow thickness variations, in the range 65,000-80,000 km3 of lava which can be related to the dynamics of flows, vent deposits, and intrusions that lava emplacement. Dietterich et al. (2015) extend from the ocean floor to the 4169 m have studied lidar data of the Hawai’i and high summit (Sherrod et al. 2007). Lava eruptions from Mauna Loa are silica-poor, Table 24: Lidar scanning parameters of Mauna Loa, Hawaii, USA3. and very fluid; eruptions tend to be non- explosive and the volcano has relatively parameter value shallow slopes. Its most recent eruption scanner type Optech Gemini was in 1984. platform fixed wing Cashman et al. (2013) report that lidar date June 21st – 27th 2009 data are revolutionizing both the visual average last and only returns per m2 on a 6.0 and quantitative analysis of lava flows. combined dataset Visualizations allow accurate mapping of spatial resolution of the final elevation model [m] 0.5 flow boundaries, particularly in vegetated 3 Lidar data acquisition and processing completed by the National Center for Airborne Laser Mapping areas and areas that are difficult to (NCALM). NCALM funding is provided by NSF’s Division of Earth Sciences, Instrumentation and Facilities Program (EAR-1043051). Grant NSF EAR-0739153 and a NASA subcontract to the Jet Propulsion Laboratory access due to temperature (active flows) (Award #1290138) also supported this data collection. Figure 49: View of Mauna Loa from Mauna Kea with cinder cones in the foreground. Photo by Lawrence Goldman. 5. Case studies 67 Figure 50: Lava flows on the Mauna Loa shield volcano are difficult to observe from the ground because they are either hot, have rough and jagged surfaces, or are covered by vegetation. Many vents and cinder cones are visible on this hill-shaded image of 1 m spatial resolution lidar data © NCALM. Figure 51: This visualization provides a very plastic representation of the surface and reveals an intricate system of cinder cones, vents, cracks, folds, and channels of a lava field on the north-eastern slope of Mauna Loa. The visualization is a combination of positive openness (65-95, 50 % opacity, darken), sky- view factor (0.65-1.0, 50 % opacity, multiply), slope (0-50°, 100 % opacity, overlay) and hillshading. 0.5 m spatial resolution lidar data © NCALM. 68 5. Case studies 5|9 bottom is covered by the sediments, gravel Blind valley Odolina, and sands, cut by some younger, up to Slovenia 25 m deep alluvial ponors and sinkholes, and the riverbed of the brook. During the Odolina at the foothills of Brkini hills normal water levels the brook sinks in is Slovenia’s most characteristic blind the riverbed immediately after reaching valley. It is a special valley with corrosion the limestone, while during higher water widened bottom, a typical contact karst levels it flows a kilometre further into morphological form shaped by a surface a 117 m deep ponor cave composed of river when it reaches the karst. Brsnica potholes and shorter channels (Figure 54) sinking stream drains a 4.3 km2 large (Mihevc 1994: 102–103). water basin and forms a normal fluvial landscape on the flysch Brkini hills (Figure 53). When it reaches limestone on Matarsko podolje, it developed Odolina, about a kilometre long and up to Table 25: Lidar scanning parameters of Odolina blind valley, Slovenia. 300 m wide valley, with sharply outlined, parameter value high, and rocky slopes, ending in an scanner type Riegl LMS-Q780 amphitheatre. It was named after a small hamlet with a castle located at its northern platform helicopter part. Close to the contact of flysch and date February 2014, exact date unknown limestone the valley is 150 m deep, and average last and only returns per m2 on a 9.8 on the southern end it is carved some combined dataset 60 m deep into the karst plain. The valley spatial resolution of the final elevation model [m] 1.0 Figure 52: Blind valley Brezovica. A view to the north. Ločica stream that formed the valley sinks underground in the foreground of the photograph. Photo by Žiga Kokalj. 5. Case studies 69 Figure 53: Odolina (1) and Brezovica (2) blind valleys. Contact borderline between flysch and limestone runs NW to SE and is very clearly recognizable. The brooks form a normal fluvial relief on flysch and form corrosion widened valleys on limestone before disappearing underground. Dolines, sinkholes, and other typical karst forms are also easy to identify. 1 m spatial resolution lidar data © ARSO, Slovenia. Figure 54: A close-up view of the end of blind valley Odolina. On such a location and without prior knowledge about the area someone would expect a stream source. The abyss where Brsnica disappears underground at high water, carrying with it tons of sediment and wood, is marked with an arrow. 1 m spatial resolution lidar data © ARSO, Slovenia. 70 5. Case studies 5|10 Creative Commons Attribution 3.0 cannot have been more than a few metres Lake floor Unported license (IGKB - Internationale lower than today. This implies that large morphology, Gewässerschutzkommission für den icebergs with a thickness of more than Lake Constance, Bodensee 2015) and documents the 110 m were calving from the retreating Germany, morphology of the lake bottom in Rhine Glacier while it still occupied the unprecedented detail. While bathymetric eastern portion of today’s lake basin. Switzerland, and models are commonly presented as depth Austria contours and/or by colour coding, the high-resolution data now available allow Due to increasing erosion and the application of enhanced visualization conservation concerns, a large number techniques. This opens new perspectives of lakeshore settlements dating from the for the analysis of the sub-aqueous Neolithic to the Bronze Age are subject of geomorphology regarding e.g. wave archaeological research. Many of these erosion or mass movements. settlements at Lake Constance are now inscribed in the multinational UNESCO An unexpected detail that now becomes World Heritage Site called ‘prehistoric pile visible are numerous iceberg scour marks dwellings around the Alps.’ (Sacchetti et al. 2012) which are most prominently seen on an up to 1.5 km wide In the years 2013 and 2014, the topography and, almost horizontal (approx. 1° slope) of the lake floor was surveyed using bench in the southern half of the lake. This bathymetric lidar for the shallow near- bench lies at a depth of 110 to 130 m below shore areas and sidescan sonar for present-day lake surface. Because of the deeper portions of the lake. The geological limitations at the lake outflow, resulting digital model, with a resolution lake water level during the Late Glacial of 3 m, is freely available under the may have been somewhat higher but Figure 55: With a surface area of approximately 570 km2 and a maximum depth of 252 m, Lake Constance is one of the largest and deepest lakes in Europe. It was formed by glacial erosion during the last Ice Ages. Large drumlin fields, mainly north of the lake, document the direction of ice flow during the last glaciation. When the glaciers retreated, the basin filled with water and became Lake Constance. Photo by Frank Numrich. 5. Case studies 71 Figure 56: Iceberg scour marks (curvilinear features) and mass movement deposits (lobate features at the foot of the steep slope) in Lake Constance. Visualisation: Laplacian-of-Gaussian overlaid by colour- coded bathymetry. 3 m spatial resolution bathymetric model © IGKB (2015). Figure 57: Bathymetric profile (P -P ) from [m] 1 2 water surface data shown in figure 56 reveals the flat 400 subaqueous bench on which most iceberg scour marks occur. 3 m spatial resolution bathymetric model © IGKB (2015). 300 P1 200 P max. depth 2 100 0 500 1000 1500 2000 2500 m 72 5. Case studies 5|11 Urban planners use lidar data for inquiry is – i.e. the fewer nearby objects cast a High-rise about resources and environment, and shadow over it, the brighter it is displayed. buildings in New for integrating new developments on An added cast shadow from the southeast York, USA vacant land with the existing fabric of direction simulates shadows cast by a people, structures, and infrastructure. morning Sun in the northern hemisphere. The magnificent height of the buildings The elevation angle of the imaginary New York City is made up of five boroughs is usually lost in 2D representations and Sun has to be fine-tuned to the height sitting where the Hudson River meets maps because it is difficult to reproduce and density of structures as well as the the Atlantic Ocean: The Bronx, Brooklyn, with analytical hillshading and colour prevailing topography. The overcast sky Manhattan, Queens, and Staten Island. coding (Figure 59). As an alternative, we illumination model with a cast shadow can When people think of New York City, used an artistic approach of conveying be combined with hillshading, or some Manhattan – the island in its core –, depth with a cast shadow (Figure 60). An other visualization technique to add detail, is often the first place they picture. A overcast sky illumination model nicely and the elevation model to add colour. densely populated borough is among renders outlines of buildings and provides the world’s major commercial, financial, a beautiful background image. The less and cultural centres. It is home to big- locally obstructed a part of the scene name attractions, world-class museums, restaurants, and concert halls. Amid its Table 26: Lidar scanning parameters of Manhattan, New York, USA. iconic sites are many of the skyscrapers, parameter value such as the Empire State Building, the Chrysler Building, the Woolworth Building, scanner type Leica ALS70 Rockefeller Center, and One World Trade platform fixed wing Center. date August 5th – 15th 2013 average last and only returns per m2 on a 6.6 combined dataset spatial resolution of the final elevation model [m] 1.0 Figure 58: Brooklyn bridge across the East River with skyscrapers of Lower Manhattan in the background. Photo by Žiga Kokalj. 5. Case studies 73 Figure 59: A combination of colour coded elevations and unsaturated shaded relief. This would be a typical two dimensional representation of a city scene. While heights of buildings can be determined to a certain extent, all the values higher than 110 m had been clipped to get a better distribution of colours. Even though the image presents Lower Manhattan with many buildings exceeding 50 m in height, their magnificence is not evident. 1 m spatial resolution lidar data © USGS. Figure 60: This image is less colourful than Figure 59 but the shadows give a more grandeur feeling of the high-rise buildings. A combination of an overcast Sky Illumination Model (100 m shadow modelling distance) and a cast shadow (35° Sun elevation and 135° azimuth). 1 m spatial resolution lidar data © USGS. 74 Glossary of terms and abbreviations 75 76 Glossary of terms and abbreviations ALS Aerial laser scanning is a surveying technique that DTM A digital terrain model is an ordered set of uses a lidar scanner on an airborne platform, e.g. sampled data points that represent the spatial a helicopter, a fixed wing aircraft, or an unmanned distribution of various types of information on the aerial system. terrain (both topographic and non-topographic), e.g. elevation, slope, slope form, rivers, ridge lines, azimuth An azimuth is an angular measurement in a break lines, etc. It usually represents the elevation spherical coordinate system. The vector from of ‘bare earth’, i.e. the shape of terrain without any an observer (origin) to a point of interest is objects on it. projected perpendicularly onto a reference plane; the angle between the projected vector and a elevation The angle of the light source above the horizon. reference vector on the reference plane is called angle The elevation angle (also altitude) is expressed in the azimuth. The reference plane for an azimuth positive degrees, with 0° at the horizon and 90° is typically true north, measured as a 0° azimuth. directly overhead. Moving clockwise on a 360 degree circle, east has azimuth 90°, south 180°, and west 270°. fovea A small, central pit in the eye retina responsible for sharp central vision, which is necessary in cuneiform Clay tablets with wedge-shaped marks that humans for activities where visual detail is of tablets represent one of the earliest systems of writing primary importance. invented by Sumerians. GIS A geographic information system is a system DEM A digital elevation model is a subset of a designed to capture, store, manipulate, analyse, digital terrain model and its most fundamental manage, and present spatial or geographical data. component. It usually refers to the elevation data organized in the form of a matrix. It usually histogram A histogram is a graphical representation of represents the elevation of ‘bare earth’. the distribution of numerical data. It therefore presents the frequency of elevations by vertical DSM A digital surface model contains elevations of rectangles, with widths equal to class interval (e.g. natural terrain features including objects on it, i.e. one meter) and heights equal to the frequency of vegetation and cultural features such as buildings, elevations in that interval. bridges, and power lines. histogram A method in image processing of contrast saturation adjustment using the image's histogram where the low and/or high values are clipped. No details can therefore be perceived below the lower clipping value (percentage) and above the higher clipping value (percentage), however the rest of the values become more contrasted. 77 histogram Histogram stretch is a method in image processing RGB Additive red, green, and blue colour model, where stretch of contrast adjustment using the image's red, green, and blue light are added together in histogram. various ways to reproduce a broad array of colours. LiVT Lidar Visualisation Toolbox is software for RVT Relief Visualization Toolbox is software for computing visualizations from raster elevation computing visualizations from raster elevation models. models. LoG Laplacian-of-Gaussian is a filtering approach SAR Synthetic Aperture Radar, a microwave radar that combines a Gaussian smoothing filter with a remote sensing technology that allows the creation Laplacian edge enhancing filter. of DSMs. LRM Local relief model is a result of a procedure that separates local small-scale features from large- SfM Structure-from-Motion. A multi-image digital scale landscape forms. photogrammetry approach to create three- dimensional models. MSII Multi-scale integral invariants is a visualization technique that has low noise sensitivity and shaded A visualization technique based on shading the is especially useful for revealing very subtle relief digital terrain model with a constant Sun azimuth topography. and elevation angle. nDSM Normalised digital surface model represents SVF Sky-view factor is a geophysical parameter that heights of features relative to the ground. represents the portion of the sky visible from a certain point. noise Noise in the digital elevation model degrades the interpretability of the data. It can result from TIN Triangulated Irregular Network is a digital data measurement or production errors, and can structure used in a geographic information system usually be observed as an irregularity in the data. (GIS) for the representation of a surface. The TIN model represents a surface as a set of contiguous, openness Openness is a visualization method based on non-overlapping triangles. Within each triangle the estimating the mean horizon elevation angle within surface is represented by a plane. The triangles a defined search radius. are made from a set of points called mass points. PCA Principal component analysis is a is a viewshed The geographical area that is visible from a mathematical procedure that convert a set of location. It includes all surrounding points that observations of possibly correlated variables into a are in line-of-sight with that location and excludes set of values of linearly uncorrelated variables. points that are beyond the horizon or obstructed by terrain and other features (e.g. buildings, trees). px A pixel or picture element is a physical point in a raster image. Pixels are normally arranged in a regular two-dimensional grid, and are usually, but not necessarily, square. Each pixel is a sample of an original image; more samples typically provide more accurate representations of the original – we say the image have a higher resolution. In geographical space the pixel carries spatial information; it defines the spatial resolution of an image. 78 Bibliography and lists of figures and tables 79 80 Bibliography and further reading Adams, D. 1982. Life, the Universe and Everything. London: Pan Conrad, O., B. Bechtel, M. Bock, H. Dietrich, E. Fischer, L. Gerlitz, Books. J. Wehberg, V. Wichmann and J. Böhner. 2015. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. ArcMap. 2012 (version 10.1). Redlands (CA): Esri Inc. Geoscientific Model Development 8 (7): 1991–2007. Batson, R. M., E. Edwards and E. M. Eliason. 1975. Computer Conway, S. J., A. Decaulne, M. R. Balme, J. B. Murray and M. C. Generated Shaded Relief Images. Journal of Research of Towner. 2010. A New Approach to Estimating Hazard the U.S. Geological Survey 3 (4): 401–408. Posed by Debris Flows in the Westfjords of Iceland. Blinn, J. F. 1977. Models of Light Reflection for Computer Geomorphology 114 (4): 556–572. Synthesized Pictures. SIGGRAPH Computer Graphics 11 Crutchley, S. 2009a. Using LiDAR in Archaeological Contexts: The (2): 192–198. English Heritage Experience and Lessons Learned. In Boeree, G. 2009. General Psychology: Perception and Interaction. Laser Scanning for the Environmental Sciences, edited http://webspace.ship.edu/cgboer/genpsyperception. by G. Heritage and A. Large, 180–200. Chichester: Wiley-html. Blackwell. Cashman, K. V., S. A. Soule, B. H. Mackey, N. I. Deligne, N. D. ———. 2009b. Ancient and Modern: Combining Different Remote Deardorff and H. R. Dietterich. 2013. How Lava Flows: Sensing Techniques to Interpret Historic Landscapes. New Insights from Applications of Lidar Technologies to ICT and Remote Sensing for Cultural Resource Lava Flow Studies. Geosphere 9 (6): 1664–1680. Management and Documentation 10, Supplement 1: e65–e71. Challis, K. 2006. Airborne Laser Altimetry in Alluviated Landscapes. Archaeological Prospection 13 (2): 103–127. Devereux, B. J., G. S. Amable and P. Crow. 2008. Visualisation of LiDAR Terrain Models for Archaeological Feature Challis, K., P. Forlin and M. Kincey. 2011. A Generic Toolkit for the Detection. Antiquity 82 (316): 470–479. Visualization of Archaeological Features on Airborne Lidar Elevation Data. Archaeological Prospection 18 (4): Devereux, B. J., G. S. Amable, P. Crow and A. D. Cliff. 2005. 279–289. The Potential of Airborne Lidar for Detection of Archaeological Features under Woodland Canopies. Challis, K., Ž. Kokalj, M. Kincey, D. Moscrop and A. J. Howard. Antiquity 79 (305): 648–660. 2008. Airborne Lidar and Historic Environment Records. Antiquity 82 (318): 1055–1064. Dietterich, H. R., S. A. Soule, K. V. Cashman and B. H. Mackey. 2015. Lava Flows in 3D. In Hawaiian Volcanoes: From Chiverrell, R. C., G. S. P. Thomas and G. C. Foster. 2008. Source to Surface, edited by E. Carey, V. Cayol, M. Poland Sediment–landform Assemblages and Digital Elevation and D. Weis, 483–505. Hoboken, NJ: John Wiley & Sons, Data: Testing an Improved Methodology for the Inc. Assessment of Sand and Gravel Aggregate Resources in North-Western Britain. Engineering Geology 99 (1–2): Doneus, M. 2013. Openness as Visualization Technique for 40–50. Interpretative Mapping of Airborne LiDAR Derived Digital Terrain Models. Remote Sensing 5: 6427–6442. 81 Edgeworth, M. 2011. Fluid Pasts: Archaeology of Flow. London: Hobbs, K. F. 1995. The Rendering of Relief Images from Digital Bristol Classical Press. Contour Data. The Cartographic Journal 32 (2): 111–116. Fernandez-Diaz, J. C., W. E. Carter, R. L. Shrestha and C. ———. 1999. An Investigation of RGB Multi-Band Shading for L. Glennie. 2014. Now You See It… Now You Don’t: Relief Visualisation. International Journal of Applied Understanding Airborne Mapping LiDAR Collection and Earth Observation and Geoinformation 1 (3–4): 181–186. Data Product Generation for Archaeological Research in Mesoamerica. Remote Sensing 6 (10): 9951–10001. Horn, B. K. P. 1981. Hill Shading and the Reflectance Map. Proceedings of the IEEE 69 (1): 14–47. Forestry Commission Scotland. 2012. Archaeological Measured Survey on Scotland’s National Forest Estate. Inverness: IGKB - Internationale Gewässerschutzkommission für den Forestry Commission Scotland. http://scotland.forestry. Bodensee. 2015. IGKB-Tiefenschärfe-Bodensee gov.uk/supporting/strategy-policy-guidance/historic- Digitale Geländemodelle Mit 10 m Und 3 m Auflösung. environment. doi: 10.1594/PANGAEA.855987 https://doi.pangaea. de/10.1594/PANGAEA.855987. Fornaciai, A., M. Bisson, P. Landi, F. Mazzarini and M. T. Pareschi. 2010. A LiDAR Survey of Stromboli Volcano (Italy): Digital Imhof, E. 1982. Cartographic Relief Presentation. Berlin, New Elevation Model-Based Geomorphology and Intensity York: Walter de Gruyter. Analysis. International Journal of Remote Sensing 31 Kershaw, A. 2003. Hadrian’s Wall National Mapping Programme (12): 3177–3194. - a World Heritage Site from the Air. Archaeological French, J. R. 2003. Airborne LiDAR in Support of Prospection 10 (2): 159–161. Geomorphological and Hydraulic Modelling. Earth Kokalj, Ž., K. Zakšek and K. Oštir. 2011. Application of Sky-View Surface Processes and Landforms 28 (3): 321–335. Factor for the Visualization of Historic Landscape GDAL Development Team. 2014. GDAL - Geospatial Data Features in Lidar-Derived Relief Models. Antiquity 85 Abstraction Library. Windows (version 1.11.0). Open (327): 263–273. Source Geospatial Foundation. Laharnar, B. 2012. Notranjska med prazgodovino in antiko. George, M. R., R. E. Larsen, N. K. McDougald, K. W. Tate, J. Ljubljana: Univerza v Ljubljani, Filozofska fakulteta. Gerlach John D. and K. O. Fulgham. 2004. Cattle Grazing Lan, H., C. D. Martin, C. Zhou and C. H. Lim. 2010. Rockfall Has Varying Impacts on Stream-Channel Erosion in Oak Hazard Analysis Using LiDAR and Spatial Modeling. Woodlands. California Agriculture 58 (3). http://www. Geomorphology 118 (1–2): 213–223. escholarship.org/uc/item/72s9c4v2. Mara, H., S. Krömker, S. Jakob and B. Breuckmann. 2010. Gulič, A. and M. Črešnar, ed. 2012. Arheološka pot po Mariboru GigaMesh and Gilgamesh: 3D Multiscale Integral z okolico. Odsek I: Zgornje Radvanje – Spodnje Hoče. Invariant Cuneiform Character Extraction. In Vodnik po najdiščih / Archaeological Trail of Maribor Proceedings of the 11th International Conference on and its Surroundings. Section I: Zgornje Radvanje – Virtual Reality, Archaeology and Cultural Heritage, Spodnje Hoče. Guide to Sites. Ljubljana, Maribor: Zavod 131–138. Paris, France: Eurographics Association. za varstvo kulturne dediščine Slovenije; Društvo mladih raziskovalcev MRM Maribor. Metternicht, G., L. Hurni and R. Gogu. 2005. Remote Sensing of Landslides: An Analysis of the Potential Contribution Haddad, D. E. and J. R. Arrowsmith. 2011. Geologic and to Geo-Spatial Systems for Hazard Assessment in Geomorphic Characterization of Precariously Balanced Mountainous Environments. Remote Sensing of Rocks. Contributet report CR-11-B. Tucson: Arizona Environment 98 (2–3): 284–303. Geological Survey. Mihevc, A. 1994. Contact Karst of Brkini Hills. Acta Carsologica 23: Heidemann, H. K. 2014. Lidar Base Specification. In Book 11, 99–109. Collection and Delineation of Spatial Data, Chap. B4, ver. 1.2, 67. U.S. Geological Survey Techniques and Methods. Miller, G. 1994. Efficient Algorithms for Local and Global http://dx.doi.org/10.3133/tm11B4. Accessibility Shading. In Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Hesse, R. 2010. LiDAR-Derived Local Relief Models - a New Techniques, 319–326. ACM. Tool for Archaeological Prospection. Archaeological Prospection 17 (2): 67–72. Minnaert, M. 1961. Photometry of the Moon. In Planets and Satellites, edited by G. P. Kuiper and B. M. Middlehurst, ———. 2016. Visualisierung Hochauflösender Digitaler 213–248. Chicago: The University of Chicago Press. Geländemodelle Mit LiVT. In Computeranwendungen Und Quantitative Methoden in Der Archäologie. 4. Mlekuž, D. 2013. Roads to Nowhere? Disentangling Meshworks of Workshop Der AG CAA 2013, edited by U. Lieberwirth Holloways. In Aerial Archaeology and Remote Sensing and I. Herzog, Edition Topoi, 109–128. Berlin Studies of from the Baltic to the Adriatic, edited by Z. Czajlik and A. the Ancient World. Berlin: Topoi. Bödocs, 37–41. Budapest: Eötvös Loránd University. 82 Mlekuž, D. and M. Črešnar. 2014. Landscape and Identity Politics Yoëli, P. 1965. Analytische Schattierung. Ein Kartographischer of the Poštela Hillfort. In Studia Praehistorica in Entwurf. Kartographische Nachrichten 15 (5): 141–148. Honorem Janez Dular, edited by S. Tecco Hvala, 197–211. Opera Instituti Archaeologici Sloveniae 30. Založba ZRC. Yokoyama, R., M. Shlrasawa and R. J. Pike. 2002. Visualizing Topography by Openness: A New Application of Image Passalacqua, P., P. Tarolli and E. Foufoula-Georgiou. 2010. Processing to Digital Elevation Models. Photogrammetric Testing Space-Scale Methodologies for Automatic Engineering and Remote Sensing 68: 251–266. Geomorphic Feature Extraction from LiDAR in a Complex Mountainous Landscape. Water Resources Research 46 Zakšek, K., K. Oštir and Ž. Kokalj. 2011. Sky-View Factor as a (11): W11535. Relief Visualization Technique. Remote Sensing 3: 398–415. Planina, J. 1954. Soča (Slovenia). A Monograph of a Village and Its Surroundings (in Slovenian). Acta Geographica 2: Zakšek, K., K. Oštir, P. Pehani, Ž. Kokalj and E. Polert. 2012. Hill 187–250. Shading Based on Anisotropic Diffuse Illumination. In Symposium GIS Ostrava 2012, 1–10. Ostrava: Technical Reitberger, J., P. Krzystek and U. Stilla. 2008. Analysis of Full University of Ostrava. Waveform LIDAR Data for the Classification of Deciduous and Coniferous Trees. International Journal of Remote ZRC SAZU. 2010. IAPS ZRC SAZU | Institute of Anthropological and Sensing 29 (5): 1407–1431. Spatial Studies ZRC SAZU. http://iaps.zrc-sazu.si/en/svf 2010. Sacchetti, F., S. Benetti, C. Ó Cofaigh and A. Georgiopoulou. 2012. Geophysical Evidence of Deep-Keeled Icebergs on the Rockall Bank, Northeast Atlantic Ocean. Geomorphology 159–160: 63–72. Sankey, J. B., N. F. Glenn, M. J. Germino, A. I. N. Gironella and G. D. Thackray. 2010. Relationships of Aeolian Erosion and Deposition with LiDAR-Derived Landscape Surface Roughness Following Wildfire. Geomorphology 119 (1–2): 135–145. Sherrod, D. R., J. M. Sinton, S. E. Watkins and K. M. Brunt. 2007. Geologic Map of the State of Hawai‘i. Open-File Report 2007–1089. Reston: U.S. Geological Survey. http://pubs. usgs.gov/of/2007/1089/. Sittler, B. 2004. Revealing Historical Landscapes by Using Airborne Laser-Scanning - A 3D-Modell of Ridge and Furrow in Forests near Rastatt (Germany). International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 36 (8/W2): 258–261. Smith, M. J., J. Rose and S. Booth. 2006. Geomorphological Mapping of Glacial Landforms from Remotely Sensed Data: An Evaluation of the Principal Data Sources and an Assessment of Their Quality. Geomorphology 76 (1–2): 148–165. Wagner, W., M. Hollaus, C. Briese and V. Ducic. 2008. 3D Vegetation Mapping Using Small-footprint Full- waveform Airborne Laser Scanners. International Journal of Remote Sensing 29 (5): 1433–1452. Woolard, J. W. and J. D. Colby. 2002. Spatial Characterization, Resolution, and Volumetric Change of Coastal Dunes Using Airborne LIDAR: Cape Hatteras, North Carolina. 29th Binghamton Geomorphology Symposium: Coastal Geomorphology 48 (1–3): 269–287. Yarbus, A. L. 1967. Eye Movements and Vision. Translated from the Russian edition (Moscow, 1965) by Basil Haigh. New York: Plenum Press. 83 List of figures Figure 1: Angle dependence of analytical hill-shading: Figure 4: Slope image of Žerovinšček Iron Age hillfort 315° azimuth illumination (A), andv 225° near Bločice, Slovenia (A). The elevation profile azimuth (B), both with 35° Sun elevation. Note (B) refers to the CD line in (A). Note that the the difference of the relief perception and the structure the profile line crosses can be quite structures that can be observed. Overgrown easily – and wrongly – interpreted as convex cultural terraces near Koboli in Slovenia as instead of concave. 1 m spatial resolution lidar evidenced by 1 m spatial resolution terrain data. 19 model. See also Figure 23. 17 Figure 5: (A 1) A histogram of the whole data range – Figure 2: Very low light source angles expose features elevations between 263 m and 536 m. As we of extreme subtleness: a standard 45° Sun are interested only in a very narrow range of elevation (A) and (B), and low light 5° Sun values between 271 m and 278 m – marked by elevation (C), all with 45° azimuth. However, a dark arrow on (A) –, it makes sense to stretch this only works in areas with very gentle relief only this range to the whole ‘colour’ palette (B). morphology, such as this example of the Site A (A 2) A histogram of stretched values of the area embanked enclosure(s) in Brú na Bóinne World shown on Figure 5. Instead of everything being Heritage Site in Ireland. 1 m spatial resolution concentrated in a ‘dark corner’ of the histogram, lidar data used with permission of Meath County the values are much more evenly distributed and Council and the Discovery Programme. Local differences in colour shade can be more easily histogram saturation is used to present (B) and perceived. 20 (C). The first to show the difference this makes when compared with normally presented shaded Figure 6: Histogram stretch to a narrow range of relief (A), and the second because the image is values. Past riverbeds of Nadiža. 0.5 m spatial otherwise too dark to expose any details. 18 resolution lidar data, West of Kobarid in NW Slovenia. 20 Figure 3: An RGB image of hill-shadings from three directions (315°, 0°, and 45° azimuth with 35° Figure 7: A hill-shaded image (A) and a histogram Sun elevation) (A), a composite of the first two stretch (B) of a local relief model, and the LRM components (B) and RGB composite of the first presented with a purposely designed colour three components of a principal component ramp overlaid on a hill-shaded image (Hesse analysis of analytical hill-shading from 16 2010) (C). 1 m spatial resolution lidar data directions with 35° Sun elevation (C). A late Žerovinšček Iron Age hillfort near Bločice, Roman camp at St. Helena, west of Kobarid, Slovenia © ZRC SAZU. 21 Slovenia. 0.5 m spatial resolution lidar data © Walks of Peace in the Soča river Foundation. 18 Figure 8: Sky-view factor image (10 m search radius in 16 directions) (A) and anisotropic sky-view factor image (B) of Žerovinšček Iron Age hillfort. Many 84 details can be perceived despite the variable Figure 17: Anisotropic sky-view factor image of elevation relief morphology. 22 models produced with different rasterization (interpolation) algorithms. Note the difference in Figure 9: Sky-view factor image (10 m search radius in 16 interpolation of the castle area and the south- directions) of ridge and furrow near Neudingen, eastern slope leading to it. Natural Neighbours Germany (A). Omitting the nearest pixels from (NN) (A) generates a very smooth terrain, the calculation can reduce the noise in the data Inverse Distance Weighted (IDW) (B) introduces so that features can be perceived more clearly steps, TIN with Repetitive Interpolation (REIN) (B).1 m spatial resolution lidar data © LGL in (C) introduces triangles, and splines (D) create a Baden-Württemberg. 23 range of artefacts. 28 Figure 10: Positive (A) and negative (B) openness image Figure 18: ‘Fish scales’ are a result of direct point cloud (10 m search radius in 16 directions) of a late rasterization (A). 0.5 m spatial resolution lidar Antiquity settlement of Tonovcov grad, Slovenia. data of an area around Besancon, France, The very complex terrain seems flattened. © University of Franche-Comte. Black stars Tops of protruding features are well delineated typically form around small protruding on positive openness image, while negative features on a sky-view factor (SVF) image when openness delineates the bottoms of hollows and calculating it with 8 or 16 directions (B). 0.5 m lower edges of cliffs. 24 spatial resolution lidar data of World War I trenches near Kobarid in Slovenia. © Walks of Figure 11: Local dominance image (10-20 m search Peace in the Soča river Foundation. 1 m spatial radius) of former field boundaries, roads, and resolution lidar data full of false depressions (C). ridge and furrow near Hügelsheim, Germany. They are a result of poor data processing that Ridge and furrow is only preserved in areas that classified many below-ground lidar points as are today covered by forest. Some other features true ground. Poštela burial mounds can be seen on the image include bomb craters, earth in the upper left corner of (C). © ARSO, Slovenia. covered bunkers, and trenches. 1 m resolution All images display SVF calculated in 8 directions lidar data © LGL in Baden-Württemberg. 25 with a 10 m search radius. 29 Figure 12: Cumulative visibility image of burial mounds Figure 19: Wave-like features are a consequence of (high values) and dolines (low values) on the poor registration of scan lines (A). Charcoal Swabian Alb. 1 m resolution lidar data © LGL in burning platforms can be seen in (A). 0.5 m Baden-Württemberg. 25 spatial resolution lidar data of an area around Besancon, France, © University of Franche- Figure 13: Accessibility image (20 m maximum radius) of Comte. Severe strip misalignment and over-very narrow agricultural terraces with vineyards ambitious resolution setting result in scrub- near Jeruzalem, Slovenia. 0.5 m resolution lidar like artefacts and fictional steps in terrain data © ARSO, Slovenia. 26 (B). 1 m spatial resolution lidar data of an area close to Tauberbischofsheim © LGL in Figure 14: Multiscale integral invariants image of a late Baden-Württemberg. Both images display SVF Roman camp at St. Helena, west of Kobarid, calculated in 8 directions with a 10 m search Slovenia. 0.5 m resolution lidar data © Walks of radius. 29 Peace in the Soča river Foundation. 26 Figure 20: An original image where only a minimum- Figure 15: Image of a late Roman camp at St. Helena, maximum histogram stretch has been applied west of Kobarid, Slovenia as evident on a 0.5 m – a whole range of values is displayed with resolution lidar data filtered by a Laplacian-of- relative differences between values preserved Gaussian convolution filter with a radius of 3 (A). Images displayed with different techniques pixels (i.e. 1.5 by 1.5 m). 27 of histogram stretch: histogram equalisation (B), standard deviation 2.5σ (C), 0.5 % saturation of Figure 16: Lidar data density maps of Castle Montfaucon minimum and maximum (D). 33 east of Besancon, France. Density of all points (A), vegetation points (B) and ground points (C) Figure 21: Visualization techniques in different calculated per every pixel (0.5 m) but plotted topographic settings. (A) Plough headlands on per m2. Density of ground points calculated and a flat plain near Endingen am Kaiserstuhl. 1 m plotted per m2 (D). 0.5 m resolution lidar data lidar data © LGL in Baden-Württemberg. (B) © University of Franche-Comte. 28 Three different types of World War I trenches with shelters on gentle NE slopes of Črni hribi, near Renče, Slovenia. 1 m lidar data © ARSO, 85 Slovenia. (C) Charcoal burning platforms in the values therefore details in black areas are hills of the Black forest. 1 m lidar data © LGL in revealed to a great extent on the original Baden-Württemberg. (D) A late Roman campo unsaturated image (C). The same is true for the on a rocky outcrop with a church of St. Helena, difference in calculation of sky-view factor with west of Kobarid, Slovenia. 0.5 m lidar data RVT and SAGA GIS. All figures are calculated © Walks of Peace in the Soča river Foundation. 36 with 35° Sun elevation and 315° azimuth using the same DEM. 1 m spatial resolution lidar data Figure 22: Different scales of observation. Sky-view © USGS. 41 factor image of a detail from a richly decorated kerbstone from Newgrange, Ireland, produced Figure 26: Visualization techniques with different relief with Structure-from-Motion of hand held features. 42 photographs (A). Scattered bones (B) and looting traces (C) on a desert floor in Peru. The Figure 27: Relief Visualization Toolbox (RVT) offers a elevation model was computed from images range of only the best techniques with just the taken with a camera on a pole and is presented essential options. It is extremely easy to use and here with a combination of local dominance and can process multiple files and all techniques in shaded relief. Grand Canyon (D) as evidenced by one go. 49 30 m SRTM data presented with a combination of sky-view factor and shaded relief. Central Figure 28: Lidar Visualization Toolbox (LiVT) supports Massif, the Jura and the Alps (E) presented with calculation of various techniques with intricate a combination of anisotropic sky-view factor and options to manipulate the settings. 49 shaded relief. 38 Figure 29: Craig canyon alluvial fan (California, USA) as Figure 23: Is it a burial mound or a hole in the ground, evident on 0.5 m spatial resolution lidar data a quarry or a hill? Someone observing the (California, USA). Coloured stretched elevations image needs to know the direction of the light over SVF and shaded relief. © Plate Boundary source, because it has a definitive impact on Observatory by NCALM, USA. 53 the perception of the landscape. 315° azimuth illumination (A), and 135° azimuth (B). Terrain Figure 30: Alluvial fan’s apron at the foot of Craig Canyon features seem inverted on (B). 39 in the Saline Valley, USA. Photo by Brian Lockett, Air-and-Space.com. 53 Figure 24: Overgrown remains of an abandoned village Novi Breg (Naubacher), Slovenia, as represented Figure 31: Parallel hollow ways in a forest NW of Volčji by local relief model technique with a filter Potok, Slovenia. Photo by Žiga Kokalj. 54 radius of 25 m (A). Remains of houses can be seen as bright rectangles. Big black blobs Figure 32: Hollow ways and headstreams below a hillfort are dolines. It is straightforward to replicate on Žiški vrh, northern Slovenia, as evident on a such an image with the minimum data that combination of sky-view factor (0.65-1.0, 25 % is provided. However, when reproduction or opacity, multiply), positive openness (70-93, relative comparison across areas is not required, 50 % opacity, overlay), slope (0-45°, 50 % opacity, combinations of visualizations can give images luminosity) and hillshading. 0.5 m spatial that are easier to read, and details about resolution lidar data © ARSO, Slovenia. 55 visualizations may be omitted (B). 1 m resolution lidar data © ARSO, Slovenia. 40 Figure 33: Hollow ways developed on a narrow passageway along the river. Volčji Potok, Figure 25: Hillshading algorithms in ArcMap (ArcMap northern Slovenia. 0.5 m spatial resolution lidar 2012) and RVT use the same method (Yoëli data © ARSO, Slovenia. 55 1965), but the results are quite different as evident from the image of skyscrapers in Lower Figure 34: Gradišče nad Knežakom hillfort in the Manhattan. The hillshading image calculated centre of the image is today totally covered by with ArcGIS (A) seems to be less clear (as if impenetrable forest. Photo by Boštjan Laharnar. 56 using a DEM of a lower resolution) than the one calculated with RVT (B). Note the distinctly Figure 35: Iron Age settlement above Knežak and its visible elevation steps on (B), especially the surroundings as evident on the classified steps on the facade of the lower right skyscraper Slovenian national lidar data (GKOT, © ARSO, (B), which are completely missing on (A). Both Slovenia) which has been processed by gLidar (A) and (B) have saturated high and low values software, rasterized to a 0.5 m grid. Many small (completely white and black areas). However, bumps and non-existent holes are visible (see RVT calculates hillshading with a full range of also Figure 18C). The image is a combination of 86 sky-view factor (0.7-1.0, 70 % opacity, multiply), slope (0-65°, 80 % opacity, overlay) and RGB slope (0-55°, 100 % opacity, overlay) and image of hillshading from three directions (R: hillshading. 57 315°, G: 15°, B: 75°), all computed on a digital surface model. 0.5 m spatial resolution lidar Figure 36: The same are as in Figure 35, but processed data © ARSO, Slovenia. 63 with different software (LASTools). Even with minimal effort (default settings), the Figure 45: Dachstein limestone strata and scree fields resulting image gives a clearer picture of the exposed at Vrh Brda (1952 m). The top of the archaeological and natural topography, despite ridge is highlighted in white, while colours losing some of the ‘sharpness’. Clearance cairns represent different orientation of the slopes. and field boundaries can be seen on the lower The visualization is a combination of sky-view part of the image. 0.5 m spatial resolution lidar factor (0.55-1.0, 30 % opacity, multiply), positive data © ARSO, Slovenia. 57 openness (65-95, 50 % opacity, overlay), slope (0-65°, 80 % opacity, luminosity) and RGB image Figure 37: One of the burial mounds at Pivola, well of hillshading from three directions (R: 315°, G: protected by brambles. Photo by Žiga Kokalj. 58 15°, B: 75°). 1 m lidar data © ARSO, Slovenia. 63 Figure 38: Burial mounds at Pivola as displayed by a Local Figure 46: Granite Dells at Watson Lake. Photo by Michael Relief Model transparently overlaid on a sky-Wilson. 64 view factor image and shaded relief. The state of preservation of the mounds differs significantly. Figure 47: Digital terrain model (DTM) of Granite The ones in the forest are well preserved (1), Dells, Arizona, USA. Much of the detail in the mounds in the Botanical Garden were damaged landscape is lost despite dedicated processing. A when the location served as a military base (2), combination of anisotropic sky-view factor (0.2- and the mounds in the fields have been nearly 1.0, 40 % opacity, multiply), slope (0-80°, 60 % completely levelled (3). 0.5 m spatial resolution opacity, overlay) and RGB image of hillshading lidar data © ARSO, Slovenia. 59 from three directions (R: 337.5°, G: 0°, B: 22.5°). 0.5 m spatial resolution lidar data © NCALM, Figure 39: Visual structure of landscape as represented USA. 65 by Higuchi’s foreground range total viewshed around Pivola barrow group. 1 m spatial Figure 48: Digital surface model (DSM) of the same area resolution lidar data © ARSO, Slovenia. as in Figure 47. Visible are trees, houses, and Reproduced with permission by D. Mlekuž. 59 cars, but more importantly also the joints in granite that are characteristic of the Dells. 0.5 m Figure 40: Lines of poles at Culbin Sands. Photo spatial resolution lidar data © NCALM, USA. 65 © Forestry Commission Scotland. 60 Figure 49: View of Mauna Loa from Mauna Kea with cinder Figure 41: Changes in elevation with added intensity cones in the foreground. Photo by Lawrence information, Culbin Sands. Lidar intensity is Goldman. 66 usually recorded at the time of scanning and can provide additional discriminating information Figure 50: Lava flows on the Mauna Loa shield volcano about the landscape. 0.5 m spatial resolution are difficult to observe on the ground because lidar data © Forestry Commission Scotland. 61 they are either hot, have rough and jagged surfaces, or are covered by vegetation. Many Figure 42: A close-up view of Culbin Sands as seen on a vents and cinder cones are visible on this hill-sky-view factor image calculated in 8 directions. shaded image of 1 m spatial resolution lidar Poles are easily identifiable as black stars. data © NCALM. 67 0.5 m spatial resolution lidar data © Forestry Commission Scotland. 61 Figure 51: This visualization provides a very plastic representation of the surface and reveals an Figure 43: A view across the ridge towards Vrh Brda and intricate system of cinder cones, vents, cracks, Bavški Grintavec in the background. Photo by folds, and channels of a lava field on the north- Jovan Cukut. 62 eastern slope of Mauna Loa. The visualization is a combination of positive openness (65-95, 50 % Figure 44: Glacially transformed karstic plateau at Velika opacity, darken), sky-view factor (0.65-1.0, 50 % vrata, Slovenia. Flat areas of clints are presented opacity, multiply), slope (0-50°, 100 % opacity, in light colours while grikes, crevices, and steep overlay) and hillshading. 0.5 m spatial resolution slopes are dark. The image is a combination of lidar data © NCALM. 67 sky-view factor (0.5-1.0, 70 % opacity, multiply), 87 Figure 52: Blind valley Brezovica. A view to the north. image presents Lower Manhattan with many Ločica stream that formed the valley sinks buildings exceeding 50 m in height, their underground in the foreground of the magnificence is not evident. 1 m spatial photograph. Photo by Žiga Kokalj. 68 resolution lidar data © USGS. 73 Figure 53: Odolina (1) and Brezovica (2) blind valleys. Figure 60: This image is less colourful than Figure 59 Contact borderline between flysch and limestone but the shadows give a more grandeur feeling runs NW to SE and is very clearly recognizable. of the high-rise buildings. A combination of an The brooks form a normal fluvial relief on flysch overcast Sky Illumination Model (100 m shadow and form corrosion widened valleys on limestone modelling distance) and a cast shadow (35° before disappearing underground. Dolines, Sun elevation and 135° azimuth). 1 m spatial sinkholes, and other typical karst features are resolution lidar data © USGS. 73 also easy to identify. 1 m spatial resolution lidar data © ARSO, Slovenia. 69 Figure 54: A close-up view of the end of blind valley Chapter 1. Charcoal burning platforms in the Black Odolina. On such a location and without prior forest, Germany, as evident on a sky-view factor knowledge about the area someone would image (10 m search radius in 16 directions). expect a stream source. The abyss where 1 m resolution lidar data © LGL in Baden- Brsnica disappears underground at high water, Württemberg. 10 carrying with it tons of sediment and wood, is marked with an arrow. 1 m spatial resolution Chapter 2. Anisotropic sky-view factor image (5 m search lidar data © ARSO, Slovenia. 69 radius in 16 directions) of Ajdovščina above Rodik hillfort, Slovenia. 0.5 m resolution lidar Figure 55: With a surface area of approximately 570 km2 data © ARSO, Slovenia. 14 and a maximum depth of 252 m, Lake Constance is one of the largest and deepest lakes in Chapter 3. Local dominance image (10-20 m search Europe. It was formed by glacial erosion during radius) of plough headlands at Endingen am the last Ice Ages. Large drumlin fields, mainly Kaiserstuhl, Germany. 1 m resolution lidar data north of the lake, document the direction of © LGL in Baden-Württemberg. 30 ice flow during the last glaciation. When the glaciers retreated, the basin filled with water Chapter 4. House mounds, palaces, platforms, walls, and became Lake Constance. Photo by Frank quarries, a possible water reservoir, and a field Numrich. 70 system of an unknown Maya settlement in the Calakmul Biosphere Reserve, Campeche, Figure 56: Iceberg scour marks (curvilinear features) and Mexico. A combination of visualizations (sky-mass movement deposits (lobate features at view factor, openness, slope, hillshading). 0.5 m the foot of the steep slope) in Lake Constance. resolution lidar data © ZRC SAZU, Slovenia. 46 Visualisation: Laplacian-of-Gaussian overlaid by colour-coded bathymetry. 3 m spatial resolution Chapter 5. A part of Maiden Castle Iron Age hillfort, bathymetric model © IGKB (2015). 71 United Kingdom, as seen on a sky-view factor image (5 m search radius in 16 directions). Figure 57: Bathymetric profile (A-B) from data shown 0.5 m resolution surface model lidar data in figure 56 reveals the flat subaqueous bench © Environment Agency, UK. 50 on which most iceberg scour marks occur. 3 m spatial resolution bathymetric model © IGKB Glossary of terms. Eye of the Sahara or Guelb er Richat (2015). 71 is a deeply eroded geologic dome in Mauritania. A combination of visualizations (sky-view factor, Figure 58: Brooklyn bridge across the East River with slope, hillshading). 30 m SRTM data © USGS. 74 skyscrapers of Lower Manhattan in the background. Photo by Žiga Kokalj. 72 Bibliography and lists. Sand dunes at Guaiuba, Santa Catarina, Brazil, as seen on a sky-view factor Figure 59: A combination of colour coded elevations image (10 m search radius in 16 directions). 1 m and unsaturated shaded relief. This would be resolution surface model lidar data © FAPESP a typical two dimensional representation of a grant 2009/17675-5. Grohmann, C.H. and A.O. city scene. While heights of buildings can be Sawakuchi. 2013. Influence of cell size on determined to a certain extent, all the values volume calculation using digital terrain models: higher than 110 m had been clipped to get a a case of coastal dune fields. Geomorphology better distribution of colours. Even though the 180–181: 130–136. 78 88 List of tables Table 13: Typical settings for calculation and visualization of Laplacian-of-Gaussian. 27 Table 14: Matrix for the suitability of visualisation techniques for selected archaeological relief Table 1: Some useful resources to go search for free lidar features in different topographic settings. 34 datasets. Host sites will inevitably change and new datasets are becoming available regularly, Table 15: Suitability of visualisation techniques so use your favourite search engine. 13 for representing selected archaeological topographical features. 35 Table 2: Typical settings for calculation and visualization of hillshading. 19 Table 16: Metadata required when presenting DEM visualisations. 39 Table 3: Typical settings for calculation and visualization of principal components analysis of hillshading Table 17: Lidar scanning parameters of Craig canyon images from multiple directions. 19 alluvial fan, California, USA. 53 Table 4: Typical settings for visualization of slope. 19 Table 18: Lidar scanning parameters of Žiški vrh and Volčji Potok, Slovenia. 54 Table 5: Typical settings for displaying elevation differentiation. 20 Table 19: Lidar scanning parameters of Gradišče nad Knežakom, Slovenia. 56 Table 6: Typical settings for trend removal and Local Relief Modelling. 23 Table 20: Lidar scanning parameters of Pivola burial mounds, Slovenia. 58 Table 7: Typical settings for calculation and visualization of sky-view factor. 23 Table 21: Lidar scanning parameters of Culbin Sands, Scotland. 60 Table 8: Typical settings for calculation and visualization of positive and negative openness. 24 Table 22: Lidar scanning parameters of Velika vrata and Vrh Brda, Slovenia. 62 Table 9: Typical settings for calculation and visualization of local dominance. 25 Table 23: Lidar scanning parameters of Granite Dells, Arizona, USA. 64 Table 10: Typical settings for calculation and visualization of cumulative visibility. 25 Table 24: Lidar scanning parameters of Mauna Loa, Hawaii, USA. 66 Table 11: Typical settings for calculation and visualization of accessibility. 26 Table 25: Lidar scanning parameters of Odolina blind valley, Slovenia. 68 Table 12: Typical settings for calculation and visualization of multi-scale integral invariants. 26 Table 26: Lidar scanning parameters of Manhattan, New York, USA. 72 The series Prostor, kraj, čas (Space, Place, Time) publishes shorter, thematically Prostor rounded studies from diverse aspects of exploration of space and time, that are based on geographical information systems and remote sensing, as well as their kraj social and cultural constructs: how people of various periods and landscapes think čas 14 about, live, feel, use, and change space and time. Series editors: Nataša Gregorič Bon and Žiga Kokalj, ZRC SAZU AIRBORNE LASER SCANNING RASTER DATA VISUALIZATION A Guide to Good Practice Žiga Kokalj and Ralf Hesse Žiga Kokalj is a research fellow at the Research Centre of the Slovenian Academy of Sciences and Arts and at the Centre of Excellence for Space Sciences and Technologies. He has general research interest in application of optical imagery and aerial laser scanning, spatial analysis, and modelling of natural processes. Ralf Hesse works at the State Office for Cultural Heritage Baden-Württemberg where he uses airborne lidar data for large area archaeological prospection. His research focuses on the visualisation and analysis of digital elevation models in the fields of archaeology and physical geography. International The book is freely available in e-form (pdf) at: http:\\zalozba.zrc-sazu.si http://zalozba.zrc-sazu.si/p/P14 ZRC Publishing