Rapid mapping application of vegetated terraces based on high resolution airborne LiDAR

Abstract The aim of this work is to define a methodology for terraced areas survey and rapid mapping in a complex environment, like Ligurian (Northwestern Italy) one, where a remarkable percentage of its surface is estimated as terraced and where the canopy coverage makes their recognition very hard. Methodology steps are the definition of LiDAR survey parameters, morphometric filtering and GIS processing for final mapping. Each phase is oriented to provide a reliable terrace mapping, also practicable in canopy-covered areas due to a particular attention to land cover influence. The work considers a case study (Rupinaro basin) close to Cinque Terre, with a mixed land cover (terraces, forest and urbanized areas). The methodology provided encouraging results detecting 448 ha of terraces, 95% of them located under canopy cover. This finding pointed out that terraces mapping cannot rely only on photo-interpretation, as canopies will hamper their detection. Mapping of these areas, frequently characterized by abandonment, is crucial while identifying potential trigger factors for slope instabilities. This case study highlighted the importance of a carefully planned production chain, that should start from LiDAR survey parameter choice, providing the best input for analysis algorithm and providing the correct identification of terraces.


Introduction
Agricultural terraces have been used in the Mediterranean area since the Neolithic by different populations (Price and Nixon 2005). In Italy, the presence of man-made terraces is widely documented in several Regions like Liguria, Piemonte, Lombardia, Campania and Veneto. Some of them, featuring peculiar traits, are included in UNESCO Heritage List, too (R€ ossler 2006). Among these regions, Liguria is characterized by a remarkable percentage of such terraced areas but, up to now, an exhaustive study of the distribution of terraced areas and their mapping does not exist, nor an official survey at regional government level (Bailly and Levavasseur 2012). There are only few estimations, based on the extrapolation from the analysis of a study area. Regione Liguria (2008) in cooperation with Genova University (POLIS Department), estimated a 30% of regional surface covered by terraced areas; according to other estimations, terraces cover the 20% of the Region (Brandolini et al. 2008). The described works are mainly based on aerial image plotting (Brandolini et al. 2003;Agnoletti et al. 2011; Van der Sluis et al. 2014) and thus limited to the visible, i.e. free of coverage, portion of the analysed territory. The employment of a survey tool capable of coping with covered areas, e.g. vegetation, is a key factor in terraced areas detection. Airborne LiDAR could assure the adequate capability but it needs a correct setting in order to fully exploit its potentials (Chauve et al. 2007;Wagner et al. 2008;Mallet and Bretar 2009). Particularly by the employment of high resolution sensors it is possible to carry out a detailed survey and to deal with excessive shading by vegetation or other obstacles (Diaz-Varela et al. 2014;Stout and Belmont 2014;Passalacqua et al. 2015;Tarolli et al. 2015). Tarolli et al. (2014) mentioned this capability for terraces detection purposes in forest covered areas. In Liguria, and in other Italian regions, terraces mapping methodologies were already applied but limited to small surfaces, e.g. 1 km 2 (e.g. Sofia et al. 2014) or even smaller. Moreover, they are mainly located in vineyard or other low covered areas (among others Tarolli et al. 2015). The main focus of the above-mentioned approaches is on the anthropic terraces detection algorithm and not on the whole work-flow. In fact, the employed LiDAR dataset is often an already available one, sometimes distributed by public agencies (Cazorzi et al. 2013). Frequently carried out for other purposes, featuring an inadequate resolution and reused for terraces detection (Ninfo 2008).
In order to take advantage of LiDAR potential, a careful planning of the survey is mandatory. The aim is to cope, in particular, with forested areas. Denser canopy cover, in fact, can affect the survey with larger errors then other types of cover (Hodgson et al. 2005) and strongly reduce the number of points classified as ground. The quality of survey data will then influence further processing phases and the capability of the selected algorithm to detect terraces features properly. A high point density will provide several advantages as the increase in the quantity of pulses reaching ground, even through vegetation cover. The filtering of high-density point cloud will provide a still significant amount of ground pulses. This higher quantity will also allow a better description of terraces features (e.g. walls, risers, benches). The execution of survey flight at lower altitude will contemporary assure a high quality orthoimage. It will be helpful for visual inspection and validation of detected terraces, in vegetation free areas. In addition to sensor settings, the choice of acquisition days is crucial. Dense foliage will intercept laser pulses preventing the correct measure of terrain. Moreover, snow cover could further obstruct the measurement (Deems and Painter 2006). Therefore, survey should be carried out in a period between snow melting and the onset of spring foliage growth. In this paper, a methodology for the acquisition and processing of LiDAR dataset aimed to the expeditious identification of terraced areas is presented and discussed. In the following sections, we present the results of the developed methodology on the Rupinaro basin. An area strongly characterized by the presence of agricultural terraces. The paper also tries to demonstrate that only using a dedicated methodology that starts from LiDAR acquisition strategy and an accurate processing is possible to identify an high number of terraces. On the contrary a reduction of density of ground point cloud could generate a strong loss of terraces identification capacity of the used algorithms.

Methods
We proposed a work-flow based on LiDAR survey and GIS-based postprocessing, briefly summarized in Table 1. It does not consider only the post processing of LiDAR data but explores all its production and processing phases. In our experience, in fact, the acquisition phase is one of the key point for an adequate restitution of the studied processes and/or morphological features.
According to the previously listed aims and subjects, several requirements were identified in order to maximize LiDAR survey performance: i. adequate point density; ii. maximum ground visibility; iii. sufficient orthophoto resolution.
Starting from a carefully planned and executed LiDAR survey; the resulting point cloud was processed by using FUSION software, provided by USDA Forest Service -Pacific Northwest Research Station (McGaughey 2013). It is an opens source software providing a visual and command line user interface. The latter is aimed to the use of its command in batch mode. The integration in a scripting language is, in fact, recommended to apply processing to a tiled dataset; in the present study R was employed for this phase (R Development Core Team 2011). In each point cloud ground pulses (ASPRS Standard LiDAR Point Class ¼2, ASPRS 2013) were extracted and converted in (.asc) ASCII raster DTMs.
Terraces detection was carried out by adapting the morphometric queries proposed by Scott and Pinter (2003). Each raster tile was processed in order to compute slope angle (SA) by using a 8 Â 8 size, square shaped, moving window and relief (RE) with the employment of a 3 Â 3 diamond-shaped kernel. Then, the two variable SA and RE were filtered according to different thresholds with the purpose of generating the best representation of terraced areas. The whole set of filtered rasters was mosaicked in a unique one and further processed. The aim of this phase was to exclude from the terraces map roads, buildings and other non-terraced flat areas, by the employment of a thematic filter (geofilter). The geofilter is generated by merging different land use and map layers depicting anthropogenic land features like houses, roads and other man-made flat areas that should be excluded from the terraces layer. Land use map layers may also be used with the purpose of pointing out terraced sectors under canopy cover or other land uses, therefore adding thematic details to terraces map.

Case study
The study area (44.336 N; 9.313 E) is located in Liguria region, northwestern Italy ( Figure 1). It belongs to the municipality of Chiavari (province of Genova) close to the world-renowned area of Cinque Terre. The watershed of Rupinaro River defines its boundaries. Land use incorporates several mixed classes including forest, olive orchards, cultivated areas as well as urban areas and industrial zones ( Figure 5, inset map). Surface covers approximately 1140 ha; Elevation ranges from sea level to Mount Anchetta (547 m). Aspect is mainly concentrated in South -Southwest faces. Slope never exceed 45-50 except for limited areas. Figure 2 summarizes main watershed morphometric features.
From the geological point of view, the study area falls within the flysch succession unit, represented by Antola Unit, which overlies the Internal Ligurian units of the Northern Apennine, represented by the Gottero unit, in the northern sector of the area. Marly limestones represent the main lithology of the area (Marroni et al. 2002).

Application of the proposed methodology
In early 2015, CNR IRPI planned a LiDAR survey of the area (Giordan et al. 2017), including Rupinaro Basin, as it was one of the most involved in the 2014 floods that hit the Province of Genova (Faccini et al. 2015). The principal aim of the survey was the acquisition of high resolution DTM and orthophoto for the identification of geohydrological instabilities activated by the flood and the detailed description of investigated basins morphology. In the whole area (400 km 2 ), we mapped 1641 landslides with a mean density of 4 landslide/km 2 (Giordan et al. 2016). The flight was carried out from January to the beginning of April, according to the presented methodology specifications. Its main characteristic are listed in Table 2. The choice of the survey interval is aimed to assure the acquisition of a foliage-free dataset and in the inland,  at higher elevation, to avoid the presence of ground covered snow. Acquired data were then processed, in order to obtain the output, georeferenced, raw point cloud, in ASPRS.las format (Chen 2007;Samberg 2007;Zhang and Gao 2008). The point cloud was then clipped to the watershed borders and contemporarily filtered to obtain only ground pulses. Therefore, it was converted to ASCII raster and processed with morphometric queries obtaining a first raw terraces layer using the R implementation of Fusion Software. The result of the second step is a map of areas that satisfy the morphometric queries proposed by Scott and Pinter (2003) To avoid the possibility that selection process includes also flat areas non-terraced, a dedicated filter has been applied. The third step is the use of a geofilter that has been generated starting from the following datasets downloaded from the geoportal of the Liguria region (Regione Liguria 2016): i. the land use (LU, 1:10,000 scale); ii. the technical map (TM, 1:5000 scale); iii. vector layers (shapefile format).
By the use of thematic codes stored in their attribute table (Bossard et al. 2000;Regione Liguria 2010); Regione Liguria (2010) layers were queried and feature extracted according to Table 3.
The selected feature were merged in a unique shapefile. A 10 m buffer further expanded the resulting output in order to provide a geographic tolerance to the filter. A dissolve operator also, optimized it in a unique geometric object. The described geofilter was then rasterized and used to remove from the raw terraces map, pixels belonging to its spatial domain. Remaining terraced pixel were counted to quantify terraces surface in the study area. In order to investigate the influence of canopy cover on terraces mapping two more filters were extracted from LU. Forest areas and olive orchards polygons, in fact, were used to further subdivide terraces layer and quantify relative coverages. At the end of these procedures a terraces map layout was finally produced. The processing was carried out by the employment of Quantum GIS (QGIS Development Team 2009), R and its GIS oriented libraries, i.e. 'raster' (Hijmans 2014). The results of the algorithm detection were then validated, by a two  Attribute table query  Built up area  Polygon  TM  Attribute table query  Built up area  Polygon  LU  Attribute table query  Factory  Polygon  TM  Attribute table query  Landfill/dump  Polygon  TM  Attribute table query  Main road  Polygon  TM  Attribute table query  Other roads  Polyline  TM  Attribute table query þ2.5 m buffer  Quarry  Polygon  TM  Attribute table query  Railroad  Polygon  TM  Attribute table query  Railroad service area  Polygon  TM  Attribute table query  Reservoir  Polygon  TM  Attribute table query  Secondary roads  Polygon  TM  Attribute table query  Water bodies  Polygon  LU  Attribute table query level approach. The first step was the evaluation by comparing the obtained map with other map layer; as the study area (like the entire Region) lacks of terraces map, the comparison layer was manually produced. Then the consecutive step was the analysis of the result in term of ground truth by comparing detected terraced areas with land surveys. In order to evaluate the quality of the output raster, in arbitrarily selected sub-basin, the terraced areas were manually plotted in GIS environment. The procedure was carried out by using, as base layer, the hillshade and the orthoimage, respectively ( Figure 3). Additionally, the study area was subdivided by a 100 Â 100 m sampling grid with the aim of randomly selecting ground validation areas to check the raster result and the DTM reliability. In the six selected sectors, GPS survey of the terraces were carried out (Figure 4(a)). By employing a Leica 1200 receiver equipped with a GSM modem and VRS real-time correction (Euler et al. 2001) terraces edges were surveyed with centimetric accuracy (Figure 4(b)). Survey points were then downloaded and converted in shapefile format and compared to the terraced area layer. Moreover in order to analyse algorithm performance with lower resolution DTMs a further comparison was carried out. With the purpose of simulating a 'generic' LiDAR product, characterized by a more common point cloud densities, by employing FUSION function, the original.las format raw point cloud was decimated. Three thinned point cloud were produced: RES1, RES2 and RES5 featured by, respectively, 1, 2 and 5 points/m 2 density. These point cloud were then processed according to the previously described procedure and results were compared by quantifying percentages of terraces extracted from each dataset.

Results
The proposed methodology allowed the speditive detection of 448 ha (Table 4) of agricultural terraces featuring the 40% of analysed area, an amount remarkably greater than previous estimates (Brandolini et al. 2008;Regione Liguria 2008) based only on aerial image plotting. If we consider only the relief sector of the investigated basin, terraces ( Figure 5) cover up to the 46% of the surface, and if we exclude roads,  houses and all the elements identified, by the geofilter, the percentage raises to 54%, and represents an evidence of how the human activity modified, during centuries, this environment. According to the obtained results and by the comparison with the shaded relief map (Figure 5), the proposed methodology performed quite well when the procedure is employed as rapid mapping tool. Even if not necessary, a visual check by the comparison with the shaded relief could also evidence residual sectors of flat areas, erroneously identified as terraced pixels or erroneously classified areas or sliver polygons (SrrMara et al. 2010). The map obtained by the proposed methodology make this final visual check fast and easy to manage reducing the operator intervention to minor corrections. These processing will provide a terraces layer suitable not only for analysis purposes but also for the employment in a cartographic representation (e.g. municipality map, regional technical map). The comparison among the entire dataset and decimated ones highlighted the reduction of detection capability by, approximately, an order of magnitude. Decimation of the original point cloud showed that the use of a low quality, or non-specifically planned, survey severely scales down the algorithm capability to detect terraced areas ( Figure 6).

Results validation
To evaluate the correctness of algorithm output, a two steps validation procedure was set up. The aim of the validation test is the definition of the goodness of algorithm results and the quality of the used DTM. For this reason, the proposed test has been organized as follow: (i) manual mapping of terraces using a shaded relief derived from DTM, (ii) field mapping of terraces for a cross-check validation of DTM. The first step was carried out selecting five sub-basins as test areas and mapping terraced sectors manually. In the end, the total amount of terraced areas has been compared with the results obtained by the use of the proposed algorithm. The use of manual recognition is the only possible solution over large areas because a map of terraces is not available for this area and, generally, for Liguria region (Bailly and Levavasseur 2012). In the test area, which has an extension of 143.47 ha, the manual mapping approach gave a result of 57%. The use of automated procedure gave a result of 59%, which can be considered a very similar result. This means that the proposed algorithm is able to recognize, automatically, all the terraced areas that are described by the DTM. We also attempted the same approach employing the available orthophoto. We plotted a map of terraces using the orthoimage, but the dense vegetation coverage remarkably reduced the possibility to recognize the presence of terraces. On the same test sectors, the use of orthophoto as base layer, allowed us to identify the presence of terraces only on the 23% of the territory. This is an additional confirmation of the inadequateness of orthoimages as base layer for terraces mapping in densely vegetated sectors. The second part of the validation test is focused on the quality of the DTM. As anticipated before, the quality of the DTM is one of the most crucial elements that influence the final results. For this reason, a field activity was conducted to identify and map terraces using GPS RTK. The measurement session was performed in different areas with a different canopy cover. In Figure 7(a), comparison between shaded relief, orthophoto, and GPS RTK is presented. The obtained result is positive because the resolution of the available DTM was able to detect and describe the presence of terraces, also in densely vegetated areas. Additionally, the comparison with field survey provided encouraging results as 82% of the measured points intersected the raster layer. During checks, several areas featuring terraces and extremely dense vegetation were observed. They resulted inaccessible, and the cover hindered an adequate GPS survey. In any case, the algorithm proved its capability also in this kind of situation by recognizing underlying terraces.

Discussion
Terraces detection provided encouraging results, showing a remarkable quantity of terraced areas even in a complex territory like the investigated one. In fact, thanks to the cross check between LiDAR ground data and land use maps, up to the 95% of these terraces were located under canopy cover. Forest and olive orchards play the main role (Figure 8(a)). These two classes are highly representative of Rupinaros land use (Figure 8(b)) and of regional landscape. These land features complicate the visual recognition of terraces. Previously described approaches would have only detected, presumably, terraces located in open areas like cultivated land or shrubland, that represents only the 5% of the basin. Additionally small areas with scarce tree cover would have allowed partial recognition, summing up, an estimated, 10% of the study area as terraced. This issue underlines the importance of the base layer to be employed as it will influence the whole mapping process. The LiDAR approach provides a better terraces detection as laser pulses are able to reach the ground even under canopy cover thus supplying a valuable starting dataset. In order to exploit this capability, the LiDAR survey needs to be carefully planned by selecting adequate acquisition parameters such as the period of flight, its elevation and the point density. Only a specifically accomplished survey can cope with these needs, the reuse of an already executed one could provide only partial results.
The matter of input data density was, already discussed by several authors (e.g. Liu et al. 2007;Yilmaz and Uysal 2016) when dealing with morphological description of landscape. The availability of sensors able to provide High Resolution Topography (HRT) is of great help in this field of research (Passalacqua et al. 2015). Airborne LiDAR is one of these techniques and is largely adopted for land survey and geomorphological feature detection (Roering et al. 2013). However, the huge amount of data retrieved by a high-density survey gives several issues in its storage and usage (e.g. Oryspayev et al. 2012). Besides data management concerns, the economic aspect is equally relevant with prices reaching values of hundreds Euro per square kilometer (Lovell et al. 2005;Johansen et al. 2010;Jakubowski et al. 2013). In large areas, in order to balance the survey requirements with budget and data handling constraints, the recommended point density is set on low values, i.e. 1-2 points/m 2 (Glennie et al. 2013). The analyses carried out by Anderson et al. (2006) define even lower density thresholds for the desired DEM resolution. In US, statewide surveys were carried out at 0.06-0.11 points/m 2 (Cunningham et al. 2002;Veneziano et al. 2002). These parameter choice is noticeable in several research surveys (Takahashi et al. 2010;Jakubowski et al. 2013) like the 2011 measurements of Vernazza and Monterosso territories (unpublished data) carried out by the authors of the present paper in the aftermath of the flood occurred in those areas (Cevasco et al. 2013). The described preflight analysis pointed out the need of a high-density survey. The detailed analysis of terrains morphometric feature require an equally detailed threedimensional dataset obtainable only by a specifically designed DTM, particularly in the presented case study, where the presence of forest stand make the survey of terraces highly susceptible to input data quality.
Concerning terraces detection, land cover typology and features can have a profound impact on the density of ground returns and on survey quality in vegetated areas (Adams et al. 2011);Hodgson et al. (2003) quantified that only 10% of the LiDAR pulses may reach the ground, with a canopy cover of 80-90%.
The identification of terraces is a key element, in particular in areas where these man-made elements are now abandoned and could be considered an element at risk. The analysis of the abandon consequences of these areas started from an appropriate analysis and identification of terraces.
In fact, land cover is a key factor in terrace stability assessment; e.g. in the study area olive orchards class also include an 'abandoned' subclass. The documented presence of abandoned sectors points out several issues. Undetected terraces may suffer lack of maintenance, especially if located in abandoned areas (Stanchi et al. 2012;Tarolli et al. 2014). Damaged and/or unstable terraces could be trigger factors for slope instabilities (Shrestha et al. 2004).

Conclusions
In this paper, the method for LiDAR dataset acquisition and rapid terraced areas recognition and mapping is presented and discussed. The developed methodology was tested on the Rupinaro basin where the 40% of the surface is recognized as terraced. The proposed method also underlines the importance not only of the DTM postprocessing, but also of the correct management of the LiDAR acquisition phase. The comparison of full resolution dataset against the decimated ones highlighted notable differences between the use of the original or the decimated surveys. This should stress the importance of the careful planning of the LiDAR survey in order to assure an adequate point density. Likewise the choice of the correct season ensure the most favourable conditions in term of foliage thus allowing the maximum number of pulses to reach bare soil.