Interpolation Routines Assessment in ALS-Derived Digital Elevation Models for Forestry Applications

Airborne Laser Scanning (ALS) is capable of estimating a variety of forest parameters using different metrics extracted from the normalized heights of the point cloud using a Digital Elevation Model (DEM). In this study, six interpolation routines were tested over a range of land cover and terrain roughness in order to generate a collection of DEMs with spatial resolution of 1 and 2 m. The accuracy of the DEMs was assessed twice, first using a test sample extracted from the ALS point cloud, second using a set of 55 ground control points collected with a high precision Global Positioning System (GPS). The effects of terrain slope, land cover, ground point density and pulse penetration on the interpolation error were examined stratifying the study area with these variables. In addition, a Classification and Regression Tree (CART) analysis allowed the development of a prediction uncertainty map to identify in which areas DEMs and Airborne Light Detection and Ranging (LiDAR) derived products may be of low quality. The Triangulated Irregular Network (TIN) to raster interpolation method produced the best result in the validation process with the training data set while the Inverse Distance Weighted (IDW) routine was the best in the validation with GPS (RMSE of 2.68 cm and RMSE of 37.10 cm, respectively). OPEN ACCESS Remote Sens. 2015, 7 8632


Introduction
Airborne Light Detection and Ranging (LiDAR), also referred to as Airborne Laser Scanning (ALS), is an active remote sensing technology which has gradually become a common tool for collecting elevation information of surface targets with high precision [1,2].Compared with the traditional photogrammetric method, the accuracies of the small footprint discrete return ALS measurements are unaffected by external light conditions [3].In addition, the high spatial resolution of ALS outperforms the use of Synthetic Aperture Radar (SAR) [4].ALS systems are capable of penetrating through vegetation and recording the terrain beneath it [5].Accordingly, ALS has been widely used for generating accurate and high spatial resolution Digital Elevation Models (DEMs) [6,7] over wide areas [8][9][10][11][12], which is essential for environmental applications.
ALS has already been adopted and accepted as a very valuable tool in forestry due to the three-dimensional nature of data [13].A wide array of vegetation structural metrics, such as tree height, biomass, crown size, leaf area index, stem volume, basal area, stand density and vertical canopy structure has been estimated (e.g., [14][15][16][17][18][19][20][21]).In this context, estimates are typically based on the height of the point cloud above a continuous gridded DEM representing the bare-Earth's surface [22].However, the raw ALS data contains a large number of points returned not only from the bare-earth surface but also from other surface objects.These non-ground/object points should be separated or classified, the so-called ALS data filtering, prior to DEM development [13].This process is the most critical step in DEM generation [23,24], which may also affect the accuracy assessment of the DEMs as some non-ground ALS points may be erroneously labelled as ground points [25].
In addition, the numerous interpolation methods developed to derive a DEM from point data vary widely in their complexity, ease of use, and computational expense, thus presenting their own advantages and disadvantages depending on the characteristics of the data sets [2,26].The fidelity with which DEMs represent the real surface has been extensively explored in the last decades [27].However, as Bilskie and Hagen [28] indicated, there is an insufficiency in the literature, as well as in available Geographical Information Systems (GIS) software to efficiently assess the vertical errors related with an interpolation method using ALS data.The selection of an appropriate interpolation algorithm and spatial resolution to generate an accurate DEM becomes an important decision, especially in uneven terrain.In fact, gridding error can comprise a very important, and often neglected, source of inaccuracy in vegetation metrics estimation [29].This may be especially relevant for canopy height model estimation in forested areas with a low ground-return sample for effective DEM surface interpolation [19,[30][31][32].In this regard, the density of ground points after filtering a point cloud varies depending on the environmental conditions.For example, in heavily wooded and vegetated areas, ground points will be particularly sparse and the DEM will typically present lower accuracy, detail, and reliability [2,22], because the laser beam penetration through the canopy can be limited.
Many local studies have explored and documented in the last decades that source data density, terrain, land cover type, interpolation method, and grid size affect DEM error [6,28,[33][34][35][36][37][38][39], although few studies have comprehensively studied the effects of all the aforementioned factors together such as Guo et al. [26] or Bater and Coops [22].Furthermore, it is relevant to study the topographic error associated with ALS data sets of low nominal point density per square meter, captured at national scale, as is the case of the new ALS data provided by the Spanish National Plan for Aerial Orthophotography (PNOA) [40].Spain has made a significant effort similar to Scandinavia and the USA to provide ALS data.In accordance with the quality levels defined by the 3D Elevation Program (3DEP) of the U.S. Geological Survey National Geospatial Program (NGP) [41], the PNOA-LiDAR project meets the Quality Level 3 (QL3), which implies a vertical accuracy of 20.0 cm RMSE and a density of 0.5 points/m 2 .This accuracy aligns with the American Society for Photogrammetry and Remote Sensing (ASPRS) 20-cm Vertical Accuracy Class.
The consideration of error distribution and error propagation in DEMs has been often neglected, perhaps because the immediacy of DEM implementation with tools commonly available in software packages that overrides any concern for accuracy and error [27].In order to gain a better understanding of the error introduced in DEM development by factors related to the territory and the ALS data acquisition, prediction uncertainty maps can be valuable tools.These maps show the spatial distribution and magnitude of potential error and can assist in the recognition of areas of low quality in the DEMs generated and derived LiDAR products (e.g., [22]).
According to Gatziolis et al. [42], nearly all evaluations of the suitability of ALS data for estimating forest structural variables have been carried out in relatively simple forest conditions with a uniform canopy, little if any understory vegetation, and gentle topography [14,15], which probably facilitate the high accuracy of ALS-derived metrics of forest canopy and bare ground extraction.However, such conditions are not common in many forested areas, including the Mediterranean forests of Spain, where little is known about the effects that their complex structure and terrain may have on the suitability of ALS-derived height estimates.Thus, the aim of this research is to assess different interpolation methods in order to generate an optimal DEM to normalize the ALS data captured by the PNOA-LiDAR mission to be used in forestry applications in the context of a study area dominated by Mediterranean pine forests and topographic variability.The research objectives of this paper are to: (1) evaluate the relative performance of six interpolation routines (natural neighbor, Triangulated Irregular Network (TIN) to raster, Inverse Distance Weighted (IDW), Australian National University DEM (ANUDEM), kriging, and point to raster) implemented in ArcGIS 10.1 (ESRI, Redlands, CA, USA) and FUSION 3.30 [43] software; (2) identify the most accurate spatial resolution for DEM creation; (3) assess the effect of terrain slope, land cover, ALS ground point density and pulse penetration on interpolation error; (4) identify the most important variables in error prediction and evaluate the error distribution applying a Classification And Regression Tree (CART) analysis; (5) and provide guidance for users of low density point clouds, as the PNOA-LiDAR, in order to select the more suitable interpolation routine.

Study Area
The study area consists of two sample sites, T1 (2 × 2 km) and T2 (4 × 2 km), located in the central Ebro valley (41°56′N, 0°56′W), sited northeastern Spain (Figure 1).The Ebro basin constitutes the northernmost semi-arid region in Europe and stretches from the Pyrenees range, in the north, to the Iberian range, in the south.This area presents a Mediterranean climate with continental features.Annual precipitation averages 350 mm and mostly occurs in autumn and spring.Moreover, the study area presents cold winters, with monthly mean temperature about 7 °C, and hot, dry summers, with temperatures about 24 °C [44].Topography is characterized by moderate to steep slopes with elevation ranging from 400 m to 750 m.a.s.l.
In the two selected sites, Aleppo pine forests (Pinus halepensis Mill.) account for 44% of total cover and pine terrace plantation only 2%.Evergreen shrub vegetation represents 25%, dominated by a mixture of Quercus coccifera L., Juniperus oxycedrus L. subsp.macrocarpa (Sibth.& Sm.) Ball and Thymnus vulgaris L., and cereal crops cover 10% of the study area.The average height of the canopy is approximately 6.5 m and the average biomass around 45 t/ha.Old stands reach heights of 12-13 m and 90 t/ha [45].During the last century, the study site has been recurrently affected by fire, with some areas being burned even twice.Two scars caused by wildfires in June 1995 and August 2008, which consumed a total of 5,300 ha of forest, are distinguishable nowadays [46].Currently, in these areas affected by fire, the vegetation is dominated by shrub species that have colonized rapidly [47].Thus, this area is characteristic of the distinctive dynamic of the Mediterranean environment, recurrently affected by wildfires [48].
In summary, the selection of T1 and T2 sites was based on the objective of our research, which is to test different interpolation methods to develop DEMs in order to normalize the PNOA-LiDAR data for forestry applications in a typical Mediterranean environment.In this regard, in order to normalize the point heights not only the evaluation of filtering procedures to classify the point cloud is relevant, but also the interpolation methods.

ALS Data Acquisition
The ALS data were provided by the Spanish National Plan for Aerial Orthophotography (PNOA) [40] and captured in several surveys conducted between 22st January and 5th February 2011, using an airborne Leica ALS60 discrete return sensor.Data were collected with up to four returns measured per pulse, and intensity values from a 1064-nm wavelength laser.The resulting ALS nominal point density was 1.5 point/m 2 with a vertical accuracy of 0.20 m.Data were delivered in three 2 km × 2 km tiles of raw data points encoded in the ASPRS laser (LAS) binary file format v. 1.1, containing x, y coordinates (UTM Zone 30 ETRS 1989) and ellipsoidal elevation z (ETRS 1989).The properties of the ALS acquisition are summarized in Table 1 [47].

Data Processing
A key step in DEM generation is the previous classification of laser returns as either on or above the ground, filtering out the aboveground returns before interpolating the ground points to generate a surface [7].In this study, the point cloud classification was performed using the algorithm designed by Evans and Hudak [49] implemented in MCC v.2.1 software.According to Montealegre et al. [50], this classification algorithm balances commission (Type II) and omission (Type I) errors and it is very suitable for forested environments.In that study, the relative performance of seven different well known filtering methods not available in proprietary software was evaluated.These methods were the progressive TIN densification algorithm (LAStools), the weighted linear least squares interpolation-based method (FUSION), the multiscale curvature classification (MCC), the interpolation-based filter (BCAL), the elevation threshold with expand window method (ETEW-ALDPAT), the progressive morphological filter (PM-ALDPAT) and the maximum local slope algorithm (MLS-ALDPAT).According to Montealegre et al. [50] results, MCC filter presented the lowest overall error (16.7%) and the more problematic cover types in filtering were sprouted scrub, stumps and woody debris, as well as terrain slopes higher than 15°.

Surface Interpolation Methods
Numerous mathematical methods for creating a raster surface from an irregular point cloud exist [13,36].In this study, we compared several commonly used interpolation methods: Natural neighbor, Triangulated Irregular Network (TIN) to raster, Inverse Distance Weighted (IDW), ANUDEM (Australian National University DEM), kriging, and point to raster.ALS-derived DEMs with spatial resolutions of 1 and 2 m were created for all interpolation routines.A brief overview of these techniques, currently available in GIS software, and their parameterization is presented below and in the Table 2.
The natural neighbor well known as the "area-stealing" or Sibson method finds the closest subset of input samples to an unknown point and applies weights to them based on proportionate areas determined by Voronoi (Thiessen) polygons to interpolate a value [22,26,51,52].
The TIN to raster is based on a set of contiguous, non-overlapping Delaunay triangles to join points in three-dimensional space.Elevation is recorded for each triangle node, while elevations between nodes can be interpolated, thus allowing the generation of a continuous surface [26].Then, the value of each output raster cell is interpolated from the TIN surface at the center of each cell [2].
The Inverse Distance Weighted (IDW) estimates the cell value by averaging the values of sample data points within its neighborhood [2,22] based on the Tobler's "first law of geography".The closer a point is to the center of the cell being estimated, the more influence, or weight, it has in the averaging process [22].The influence of known points on the interpolated values based on their distance from the output point can be controlled by defining the power.For example, a power of two is often well suited for deriving raster Canopy Height Models (CHM) considering the shape of the tree canopy and its variation in elevation.However, in order to generate DEMs, different values can be suitable.Bater and Coops [22] obtained better results using a power of three in IDW interpolation to estimate terrain elevation, than using a lower power.The characteristics of the interpolated surface can be controlled by applying a fixed or variable search radius, which limits the number of input points that can be used for calculating each interpolated cell [2,26].The topo to raster or ANUDEM uses an iterative finite difference interpolation technique specifically intended for terrain modelling that more closely represents a natural drainage surface, developed by Hutchinson [53].Although ANUDEM is capable of incorporating different additional data (e.g., contours or drainage), only ALS ground returns were used for surface development [2,36].
Kriging is an advanced geostatistical procedure that generates an estimated surface from a set of points with z values [2,54,55].It is based on the regionalized variable theory that assumes that the spatial variation in the phenomenon represented by the z values is statistically homogeneous throughout the surface.In this study ordinary kriging, one of the most commonly applied kriging approaches, was evaluated [2].It assumes that the variation in elevation values is free of any structural component (drift) [26].All parameters are determined by weighted least squares methods, which are commonly used to fit semivariogram models [26,56,57].
All the aforementioned methods were implemented in ArcGIS 10.1 software (ESRI, Redlands, CA, USA), while the following was performed with the "Gridsurfacecreate" command included in FUSION LDV 3.30 [43], given the widespread use of this software by researchers in forestry applications.This method identifies the point or points within each output raster cell and assigns an elevation value to the cell based on the averaging z value of those points.If there are no points within the output cell, this is filled using the neighboring cell heights.

DEM Accuracy Assessment
ALS ground returns were randomly divided into prediction (training) and validation (test) data sets, consisting of 7,585,872 (80%) and 1,896,468 (20%) points, respectively.The training data set percentage was selected in order to ensure the generation of DEMs with 1 and 2 m resolutions, while the test data set was used to assess vertical errors in elevation without compromising the integrity of the ALS data [28].In this sense, the ALS ground returns in the prediction data set presented a nominal point density of 1.3 points/m 2 in T1 and 0.6 points/m 2 in T2, corresponding to a point spacing of 0.86 and 1.32 m, respectively.It was not intended to test the absolute geodetic accuracy of the DEMs because the validation data were subject to the same degree of positional error as the prediction data, i.e., less than 0.30 m [22,28].
Additionally, a complementary validation was performed using a finite sample of 55 high-accuracy geodetic control points collected with the Leica VIVA GS15 CS10 GNSS real-time kinematic (RTK) global positioning system and located randomly but ensuring covering the whole variability of the study area (see location of the points in Figure 1).The 49% of the ground control points were taken on pine forest, 22% on scrub, 18% on a burned area and the remaining 10% on crops and grasslands.

Error Analysis
Once the ALS-DEMs were developed using the prediction data set, the validation ALS data set and ground truth GPS checkpoints were used to compare the biases and accuracies of the surfaces [39].The vertical error of every point in the validation data with respect to the predicted value in each DEM, generated with different methods at 1 and 2 m resolution, was calculated using the following Equation (1): where E is the error at location (x,y), Pz is the predicted value of the DEM at location (x,y), and Mz is the measured value from the validation data, both ALS and GPS data set at location (x,y).Furthermore, other global statistics to assess the overall performance of the interpolation routines, such as Mean Error (ME), Mean Absolute Error (MAE) and Root Mean Square Error (RMSE), were computed [22,28,33,37].See Equations ( 2)-(4): where n is the number of test points.
In addition, with the aim of assessing the effect on interpolation errors of terrain slope, land cover structure, return density and canopy pulse penetration, following Bater and Coops [22] approach, MAEs derived from the validation ALS data set were summarized across the range of values of each variable (Figure 1).It should be noted that this analysis was not performed with 55 GPS points because the sample is small to ensure statistical validity of the results.
In order to examine the effects of terrain slope on the error, a slope steepness model was derived from all ALS ground returns in the data set.This slope map was reclassified into three categories, specifically: 0°-5° gentle; 5°-15° moderate; >15° steep slopes (Figure 1a).
A structural analysis of the variability of land cover existing in T1 and T2 sites was performed using the land cover map from the CORINE (Coordination of Information on the Environment) program of the European Commission actualized to 2012.Five structural classes were obtained in the study area (Figure 1b): Coniferous forest, terraced reforestation, burned area, scrub, and crops and natural grassland.
With respect to the influence of ALS return density, a raster surface was generated computing the number of ground returns within each pixel of 1 m 2 .This continuous variable was then discretized into 5 ranges of ground return density: 0-0.5, 0.5-1, 1-1.5, 1.5-2, >2 points/m 2 (Figure 1c).
As for canopy pulse penetration, a canopy closure model was derived from the ALS data set using the "Cover" command included in FUSION LDV 3.30.[43].The proportion of the pulses that penetrate canopy and reach the ground was calculated using a 15 m × 15 m grid and a ground tolerance of 1 meter.Then, canopy pulse penetration was categorized into 4 classes: 0%-25%, 25%-50%, 50%-75% and 75%-100% (Figure 1d).
Finally, in order to analyze to what an extent the differences in the error obtained by each interpolation methods are statistically significant, a Kruskal Wallis analysis was performed.

Error Prediction
A CART analysis [58,59] performed with IBM SPSS Statistics 20 was used as an exploratory data technique to uncover those variables having the most influence on DEM error and to examine the spatial distribution of the prediction uncertainty.This is interesting and necessary information to control the prediction uncertainty in ALS derived products that could be used by forest managers.Following Bater and Coops [22] approach, MAE computed previously with the test data set (20%) was related to the variables described in the previous section.CART is a non-parametric modeling approach that can explain the response of a dependent variable from a set of independent continuous and/or categorical variables using binary recursive partitioning of the data [59].This leads to increasingly homogeneous subsets, based on independent variable splitting criteria using variance minimizing algorithms.The dependent data are partitioned into a series of descending left and right child nodes derived from parent nodes.Once the partitioning has concluded, the child nodes are designated as terminal nodes in which all cases have the same value for the dependent variable, i.e., they are homogeneous or "pure" nodes and do not require further splitting [58].CART is a procedure of data mining tools widely used in land use change modeling (e.g., [60]) and commonly applied to remote sensing data as a rule-based classification design (e.g., [61][62][63]).The output from a CART analysis is a series of logical if-then conditions ending in terminal nodes predicting the value of the response variable [22].These conditional rules, generated by the decision tree, were implemented in ArcGIS 10.1 software (ESRI, Redlands, CA, USA) using the raster layers presented in Figure 1, in order to produce a categorical map of prediction uncertainty.

Error Analysis
Global statistics for the DEM validation using the ALS test data set and the GPS benchmarks are presented in Tables 3 and 4, respectively.The ME shown in Table 3 was sub-centimeter using all interpolation algorithms with both spatial resolutions (1 and 2 m), with the exception of point to raster interpolation, which underestimates the prediction by more than 1 cm.Overall, the results show that interpolation biases were negligible in accordance with Bater and Coops [22] and Gallay et al. [64].Analyzing the results in more detail, with the exception of point to raster model, the DEMs derived with the rest of methodologies underestimate ground elevations at 1 m resolution, while ground elevations were over-predicted at 2 m resolution.RMSEs ranged from 2.68 to 17.67 cm, decreasing for all interpolation routines with an increase in spatial resolution from 2 to 1 m.MAE was also very consistent with respect to the method used, presenting a higher variability in both spatial resolutions.For instance, the best interpolation routine was TIN to raster with a MAE of 11.73 cm in 1 m resolution DEM, which increased to 16.94 cm in 2 m resolution DEM.Therefore, the smaller the grid sizes, the lower MAE and RMSE in accordance with the results obtained by Bilskie and Hagen [28].The highest range of error was obtained by kriging, point to raster and natural neighbor interpolation methods, while ANUDEM, IDW and TIN to raster presented the lowest ones, below 7 m.
In general, the results of Table 3 show that TIN to raster interpolation method is the optimal solution in 1 m resolution DEM generation, presenting the lowest RMSE, MAE and SD values (2.68 cm, 11.73 cm and 0.16 m, respectively).However, natural neighbor and ANUDEM obtained also good results with lower RMSE and MAE values, in comparison to IDW, kriging and point to raster interpolation methods.In fact, kriging method performed slightly better in DEMs generated at 2 m cell size, being 2.30 cm more accurate than TIN to raster method.The main drawback of kriging methodology is the processing time required which is almost three times longer than in IDW interpolation routine, as well as the flexibility of kriging, which can require a lot of decision-making [2].It should be noted that ANUDEM methodology performed relatively well in both resolutions considering that drainage enforcement was not applied.
Table 4 shows the overall vertical error statistics obtained with the 55 checkpoints captured with the high precision GPS.According to Liu [25] and considering the difficult field conditions, the GPS observations presented a high accuracy (vertical and horizontal accuracy of 2.38 and 1.32 cm, respectively).The difference between survey elevations and the ALS-derived DEM elevations was higher (around 30 cm) in terms of ME, compared to the values obtained previously and reported in the Table 3.It should be noted that the GPS sample is not affected by the horizontal and vertical errors that present the ALS training data set (see Table 1).This fact contributes significantly to the error detected with the reference GPS benchmarks validation.Although the 55 survey points represent a smaller independent sample, this approach is extremely valuable considering the time-consuming character of field surveying [25].In this regard, the RMSE values obtained in the statistical analysis confirm that 1 m resolution DEMs are reasonably better than DEMs generated at 2 m resolution using the PNOA-LiDAR data set.In this case, the IDW method achieved the best results applying both spatial resolutions (1 and 2 m), obtaining RMSEs of 37.10 cm and 40.60 cm, respectively.On the contrary, point to raster method presented the worst results, both at 1 m and 2 m resolution, with RMSEs of 50.90 and 63.00 cm, respectively.These low accuracy results match those shown in Table 3.The ME indicated a general overestimation of the elevation in all interpolation routines.This can be explained by the presence of systematic errors related with the filtering procedure.It is relatively frequent in ALS-derived DEMs of forest areas, where non-ground returns such as low vegetation and logs, are included as ground returns and subsequently in the DEM interpolation process, resulting in positive MEs [30,39].The greatest range between minimum and maximum errors was obtained with the point to raster interpolator (1.90 m), indicating an exaggeration of the elevation errors.In addition, it should be stressed that vertical accuracy of GPS points is eight times better than the ALS data set.This fact, could also explain the higher RMSE, MAE and ME values obtained with the ground control points with respect to those obtained with ALS test data set, which is much more dependent of error acquisition and filtering process.Figure 2 shows the results of the analysis of the influence of terrain slope, land cover, ground return density and canopy pulse penetration on the accuracy of the ALS-derived DEMs.In relation to the slope steepness, an increment from 0° to more than 15° increases the MAEs in all interpolation algorithms in more than 10 cm (Figure 2a).The finer spatial resolution DEMs presented higher accuracy than the coarser ones [22].The 1 m surface created using TIN to raster method shows the lowest MAE across the range of slope classes.
On the contrary, MAE decreased with an increase in the penetration rate of laser pulses reaching the ground (Figure 2b).TIN to raster method performed the best when penetration rate was 0% to 25%, with a MAE of 14.32 cm, being also the most appropriate in all penetration rates.In this regard, among the interpolation routines tested, kriging presented the worst results at 1 m resolution.It should be noted that all methods, except kriging at 1 m resolution, experiment a slightly increase of approximately 1 cm in their MAEs from the 0%-25% penetration class to the 25%-50% category.This effect is prone to be due to the error introduced by "Cover" command, inherent to the parameters chosen to map this variable for the analysis.
Land cover also affected interpolation accuracy (Figure 2c).The highest MAEs occurred in the coniferous forest class, which is characterized by the presence of Aleppo pine trees.In contrast, MAEs located in terraced reforestation were slightly lower, as a result of the gaps in the forest canopy allowing a larger proportion of ALS returns to reach the ground.In general, whether land cover structural complexity decreases, the error also decreases.As expected, the 2 m resolution DEMs present lower accuracy than 1 m resolution surfaces across all land cover classes.TIN to raster method was the best method to interpolate ALS data in presence of scrub and pine forests, as well as sprouted scrub, abandoned logs, stumps and woody debris typical of a burned area, presenting errors ranging from 8.45 to 15.60 cm.However, this method was slightly less suitable for crops and grasslands land covers than ANUDEM one, as the last presented a lower MAE of 6.50 cm, compared with the 6.69 cm of error of the TIN to raster method.
Figure 2d shows the effects of ALS ground return density on interpolation accuracy.The greatest MAE was lower than 23 cm and was produced in DEMs generated at 1 m resolution using the point to raster routine in areas with the lowest point densities, i.e., between 0 and 0.5 points/m 2 .Again, this method achieved the worst results, especially in 2 m resolution DEMs.In general, the rest of routines performed slightly better at 1 m resolution than at 2 m resolution.In addition, TIN to raster method presented the lowest MAEs with the exception of the ground return density range greater than 2 points/m 2 .In this category, ANUDEM (MAE of 9.21 cm) and IDW (MAE of 9.22 cm) had a better performance.
Finally, the results obtained after performing the Kruskal Wallis test show that the differences between the errors obtained applying different interpolations routines are statistically significant.In summary, DEMs generated at 1 m resolution present a higher accuracy that those of 2 m resolution and the TIN to raster method seems to be the most suitable one to interpolate ALS data of low point densities (<0.5 points/m 2 ) in Mediterranean forested environments characterized by a variable slope steepness and a relatively complex landscape.

Error Prediction
In order to ensure consistency in the analysis, the DEM generated at 1 m resolution using TIN to raster methodology was used in the final accuracy assessment.The error prediction map created from the CART analysis is presented in Figures 3 and 4. Similar to Bater and Coops [22], the CART analysis indicated that the most important predictor variables in interpolation error were terrain slope and point density, but also land cover, which determines the amount of returns reaching the ground.
As expected, areas with a combination of high slope (>15°) and low ground points density (≤0.30points/m 2 ) were the most prone to interpolation error.However, in terrain with slope steepness lower than 15°, prediction uncertainty was very low.This method is advantageous where the combined effects of the two predictor variables are less intuitive, for instance, high slope (>15°) and high point density (>1.07 points/m 2 ).In the end, using prediction uncertainty maps may help in the detection of potential problems with ALS-derived vegetation height estimates in those areas where the DEM surface is uncertain [22].Figure 4 shows the classification of the study area into categories of prediction uncertainty.The eastern half of T2 showed higher uncertainty, particularly compared to the southwest of this test site.Good results are shown in the flat-bottom valleys occupied by field crops and in the burned area, except in zones occupied by sprouting shrub vegetation and abandoned logs.Furthermore flight-overlaps strips can be observed in T2, which implies more point density and therefore less uncertainty in DEM surface.In general, T1 test site shows moderate error uncertainty, mainly due to the hilly relief, the high point density in the area, but also due to the presence of the pine forests and the shrub vegetation cover.Figure 4 corroborates that topographic gradient and a low point density are the main factors in DEM error.In fact, vertical errors were observed in areas where the local variability in the terrain (e.g., topographic slope) was large and when ALS point count was low (i.e., ≤0.30 points/m 2 ).

Discussion
The aim of this research was comparing performance of different interpolation techniques to derive gridded DEMs in a Mediterranean forested region in order to normalize the PNOA-LiDAR data.Our results establish the first baseline for potential users of low density point clouds, in absence of information describing the suitability of interpolation parameters in areas occupied by Aleppo pine forest mixed with evergreen shrub.This is a contribution to other researches like those of Rees [65], Lloyd and Atkinson [12] and Mitášová et al. [66].
Bater and Coops [22], and Lloyd and Atkinson [67] pointed out in their research that no interpolation method is universally superior since ground return spacing, raster pixel size, the complexity of terrain morphology, and the assumptions of a given interpolator affect the ability of interpolation routines to generate accurate DEMs.In our study, the RMSE obtained in the DEM validation using the ALS test data set, varied with different methods and resolutions.However, our research confirmed that natural neighbor, IDW, kriging, ANUDEM, and TIN to raster do not differ greatly in terms of their global RMSEs and MAEs for resolutions 1 and 2 m (see Tables 2 and 3).The ANUDEM and IDW interpolators were the most conservative routines obtaining the lowest range of error.In this sense, Bater and Coops [22], found more conservative the linear and natural neighbor methods.In our case, at 1 m resolution, TIN to raster and natural neighbor, the two best interpolators, had an RMSE of 2.68 and 2.95 cm, respectively, while point to raster, the worst interpolator, resulted in an RMSE of 6.64 cm.This indicates that the errors produced by interpolators are as significant as the measurement errors and should be considered when generating high quality DEMs from ALS data [26].On the other hand, MEs showed that the two best interpolators slightly overestimated the ground elevation.This effect is frequently encountered when working with ALS data in forested areas where is usual to find a positive bias as the point cloud is misclassified due to the presence of dense low lying vegetation under the tree canopy [5,32,39,68].This overestimation is attributed to the reduced number of ALS ground points used to interpolate each grid centroid [36] but also to the use of vegetation point as ground points to generate the DEMs.Despite this fact, overestimations were limited in our research to 0.59 cm and 0.03 cm, using TIN to raster and natural neighbor methodologies, respectively.This would likely have little impact on attributes of vegetation structure derived from these DEMs.In general the ME values around zero obtained in all tested routines to generate DEMs, suggest unbiased predictions according to Gallay et al. [64] and Bater and Coops [22].
The statistics fall below a centimeter level with the exception of point to raster interpolation, where the predictions were systematically underestimated by less than 3 cm considering the MEs values.
The supplemental validation of DEMs with the 55 ground surveyed measurements with a high precision GPS provided insight into the absolute accuracy of the bare-earth surfaces.Our results fall within the typical RMSE values reported by other empirical studies, which ranged from 0.14 to 1.50 m, depending on the operational aspects of ALS and environmental conditions [6,38].
The RMSE values showed that 1 m resolution DEMs are reasonably better than 2 m resolution DEMs.As confirmed Gonga-Saholiariliva et al. [27], it was shown that the larger the grid-cell size, the lower the accuracy in DEMs.However, the higher MEs obtained in general, ranging from 40.00 to 29.70 cm with the GPS ground control points suggest the necessity of a thorough analysis, using more control points for a higher level of confidence in the validation results.However, it was not possible to improve the number of checkpoints, as according to Liu [25] the collection with GPS of a large number of high-accuracy checkpoints was a time-consuming task, leading to an increase of the costs of the study.
The results of ALS-derived DEMs validated with the GPS benchmark tended to overestimate the reference ground elevation [6].The IDW interpolation method presented the best results, RMSE of 37.10 cm, using a power of 1 and a variable search radius with 24 minimum points.As pointed out Gallay et al. [64], the IDW is an approachable method in proprietary as well as open source GIS or statistical software and it is relatively easy to parameterize.However, Renslow [2] underlined that in some ALS data sets where the point density varies widely, this method can be a challenge since different densities of points often dictate different parameters for best results.In this regard, the amount of nearest neighbors used in IDW interpolation may seem high.However it is known that a very small number of points are prone to cause artifacts in the DEM, at least in certain cases based on the interaction with the spatial pattern of the point cloud.
Additionally, computation time should be considered when choosing the appropriate interpolation method, although the absolute computation time may change under different computation conditions, such as the computer's CPU, available memory, and software used [26].The computations for this research ran under a Windows™ server with Intel ® Core™ i5 3.10 GHz processors and 8.00 GB memory using ESRI ArcGIS 10.1 (ESRI, Redlands, CA, USA) and FUSION LDV 3.30 [43].Natural neighbor, point to raster and TIN to raster prove to be simple and fast methods (60, 120 and 198 seconds to generate each DEM, respectively), while IDW and kriging are moderate in computation time (240 and 540 seconds) compared with the ANUDEM (1140 seconds), which is the slowest one.These results are in line with those obtained by Guo et al. [26] where the simplest methods have the best processing time.
In terms of the accuracy level considering the pixel size, we have demonstrated, like Bilskie and Hagen [28] that MAE and RMSE generally increase along with larger DEM grid cells.As the DEM becomes coarser, it is unable to describe sub-scale undulations of the ground surface that are better represented by higher resolution DEMs.Nevertheless, the findings of Rees [65], Liu et al. [11] and Smith et al. [69] show that the choice of interpolation method is less influential when a surface is interpolated to coarser resolutions than the resolution of the input data.In any case, our research suggests that PNOA-LiDAR-derived DEMs with pixel size similar to point density, 1 m, achieved very good results.In this context, Behan [70] quantified the error within models produced from different interpolation algorithms and obtained the highest accuracy in surfaces created using cell sizes with spacing analogous to the original points.
Ultimately, TIN to raster interpolation is usually preferred due to the overall simplicity nature of its performance, as it has no adjustable parameters, so no user-tunable variances are introduced, and to its efficiency in processing [2].In this sense, TIN to raster interpolation is also the best option for interpolating the ground returns of PNOA-LiDAR data in a forested Mediterranean environment since it presents a consistent accuracy and relative conservative predictions.
On the other hand, analysis of the effects of terrain slope, land cover, ALS ground point density and pulse penetration on DEMs accuracy showed that all factors influence MAEs.Although it is known that ALS-derived DEMs are less sensitive to terrain slope than those derived from digital photogrammetry [6,38], topographic gradients are a significant factor in DEM error as corroborated by our research and others such as Hodgson and Bresnahan [38], Su et al. [32], Gallay et al. [64] and Bater and Coops [22].As Aguilar et al. [6] suggested, ALS planimetric error may be relatively high (up to 0.30 m for PNOA-LiDAR mission) and also may be directly translated to vertical errors on sloping surfaces.In our research, double MAE values were obtained in areas with moderate slope steepness (from 5° to 15°) in comparison with areas of low slope steepness (from 0° to 5°) in 1 m resolution DEM.Similar patterns were obtained by Hodgson and Bresnahan [38], who observed elevation error in steeper slopes (about 25°) twice of those observed on low slopes (e.g., 15°).Like Bater and Coops [22], and increase in slope steepness from 0° to more than 15° caused decimeter-level increases in MAEs in all interpolation routines, especially in point to raster method.Except for this method, the rest of interpolation methods were similar in their accuracies, although the kriging and IDW routines appeared to be more sensitive to changes in slope.
Hodgson and Bresnahan [38], Aguilar and Mills [39] and Aguilar et al. [6], also indicated the influence of land-cover in the accuracy of DEMs.In general, as the structural complexity of the land cover decreased, the MAEs obtained in the DEMs generated also diminished.Our results confirm, as Bater and Coops [22] and Hodgson and Bresnahan [38] pointed out in their studies, that higher MAEs (ranging from 15 to 40 cm) occurred in areas with tall canopy vegetation, covered by dense coniferous forest.However, Hodgson and Bresnahan [38] found the largest RMSE in areas covered by scrub.In connection with this, the presence of vegetation can limit ground detection, due to a decrease in the canopy pulse penetration.Nevertheless, this is not only a deficiency of ALS data, but also in DEMs created from stereo photogrammetry, radar, or ground surveying, where the accuracy and reliability of the surface generated is usually lower in vegetated areas than in open areas [2,25].Laser energy often fails to penetrate a dense vegetation canopy resulting in last returns that are well above the true ground surface [13].Our study confirms that the TIN to raster method was the most appropriate for the different penetration rates, as it had the lowest MAEs in the ground point density categories, ranging from 14.98 to 9.24 cm, with the exception of areas with ground penetration densities greater than 2 points/m 2 , where ANUDEM presented one centimeter more in accuracy.This support the conclusion of Hu et al. [71], who confirmed that increasing sampling density implies a decrease in interpolation error.
Finally, the CART analysis indicated that topographic variability and sampling density have significant influence on the accuracy of ALS-derived DEMs but also the land cover, which determines the amount of returns reaching the ground.As pointed out by Guo et al. [26], whether the complexity of the terrain increases, the uncertainty in the derived DEM also increases.In this regard, our results are similar to previous studies conducted by Bater and Coops [22] and Hodgson and Bresnahan [38], who indicated that the pattern of highest magnitude error was observed to occur in the areas of greater surface roughness.In addition, Aguilar et al. [72] found that morphology has the greatest influence on DEM quality, followed by the sampling density and interpolation method.Similar to Aguilar et al. [6], the error in ALS-derived DEMs is not very sensitive to change in point density in areas of low average slope as can be seen in Figure 3.According to Aguilar et al. [6], it should be noted that the total vertical error of the MDE can be disaggregated into three main components: (i) the error from ALS data capture; (ii) the error due to filtering method; and (iii) the error of interpolation method and gridding.
In order to improve the results presented in this study, we consider important for future investigation providing a replicated analysis of density reduction to suggest possible strategies for reducing ALS data sets sizes as Anderson et al. [36] performed.As data sets become more widely available across larger areas and, subsequently, data set size will often be prohibitively large, the computational requirements for handling such data will become an even greater issue [36].With a reduction in data, a more usable and operationally sized elevation data set will be possible, increasing the efficiency in terms of storage and manipulation [23].Moreover, an empirical model to quantify the relationship between DEM and the influencing factors analyzed in this paper could be developed following the trend proposed by Aguilar et al. [29] for estimating global and absolute accuracy and predicting the error budget after applying ALS data filtering and gridding processes.

Conclusions
The selection of an appropriate interpolation method and spatial resolution becomes an important decision in DEM generation.This paper focuses on the assessment of six interpolation methods to generate an optimal DEM in order to normalize the ALS data captured by the PNOA-LiDAR mission to estimate vegetation structural metrics in a Mediterranean forested landscape.The interpolation methods analyzed include natural neighbor, IDW, ordinary kriging, ANUDEM, TIN to raster, and point to raster approaches.A collection of DEMs was generated with a spatial resolution of 1 and 2 m, according to the ground point density.Then, the accuracy of the ALS-derived DEMs was assessed with a test sample of ALS points and complementary with an independent reference set of 55 ground control points collected randomly on foot with a high precision GPS.The results of the validation using the ALS test samples of points showed a higher accuracy of DEMs created with the TIN to raster interpolation method at 1 m resolution grid.On the contrary, kriging interpolation was the best at 2 m resolution DEM.Poor accuracy was achieved with point to raster routine, which is considered the simplest means of converting point data to a raster surface.The high RMSEs obtained in general with the GPS control points, showed that the IDW present the lowest RMSE both, at 1 and 2 m resolution.Overall, the results confirmed that 1 m resolution DEMs present a higher accuracy than 2 m resolution ones.Additionally, with the purpose of examining the effect of terrain slope, land cover, ground point density and pulse penetration on interpolation error, the study area was stratified by these variables.Based on the error statistics computed, we concluded that the TIN to raster interpolation was the optimal solution in any terrain slope steepness, in areas with low point densities (below 0.5 points/m 2 ) and complex land cover, such as scrub and pine forests, as well as sprouted scrub, abandoned logs, stumps and woody debris, typical of a burned area.Finally, the CART analysis allowed us to conclude that areas with a combination of high slope steepness (above 15°) and low point density (below 0.3 points/m 2 ) were the most prone to present high interpolation errors.

Figure 1 .
Figure 1.Study area with the two test sites (T1 and T2) and factors influencing DEM accuracy: (a) Terrain slope; (b) Land cover; (c) Ground return density; (d) Canopy pulse penetration.The red triangles denote the locations of the reference GPS benchmarks used in the accuracy assessment overlaying a high spatial resolution PNOA-orthophotography captured in 2009.

Figure 2 .
Figure 2. Effects of terrain slope (a); canopy pulse penetration (b); land cover (c) and (d) ground return density on MAE of interpolation.

Figure 3 .
Figure 3. Classification tree resulting from CART analysis of absolute errors for a 1 m resolution DEM created using TIN to raster interpolation.Each node (square) is labeled with average absolute error (Mean), standard deviation (S.D.) and the number (N) of points in that group.The model is read from top down until terminal nodes predicting the vertical error from selected variables appear.

Figure 4 .
Figure 4. CART-derived prediction uncertainty map over the high spatial resolution orthophotography (PNOA-2009) used as backdrop for 1 m TIN to raster DEM.For this surface, slope, ground return density and land cover were the best predictors of interpolation error.

Table 1 .
Airborne laser scanning (ALS) data specifications and acquisition properties.

Table 2 .
Interpolation routines, most important advantages and disadvantages and the parameterizations tested in this study.
(TIN) to raster It is simple and computationally efficient.If point density is lower than the output cell size, the triangle of the intermediate TIN will be transferred to the output DEM.Linear and natural neighbor methods were tested to create the raster surface from the TIN.Inverse Distance Weighted (IDW) It requires a moderate decision-making and can also be computationally intensive.Power of 0.5, 1, 1.5, 2, 2.5 and 3, and a variable search radius with 6, 12 and 24 minimum points were tested.ANUDEM It allows the incorporation of spatial restrictions in the interpolation process, such as contours, streams, etc.Its primary purpose is to create a surface suitable for hydrologic modeling.It is extremely computationally intensive.Surfaces were created with drainage enforcement both on and off.KrigingIt requires a lot of decision-making and it is very computationally intensive.The fitted model of the semivariogram was "Gaussian".Sector types of 1, 4, 4 with an offset of 45° and 8, with 2 to 5 neighbors were tested.Point to rasterIt is the simplest method and it is very computationally efficient.Mean is sensitive to extreme values/outliers, especially when the sample size is small.Not applicable.

Table 3 .
Global statistics summarizing validation errors obtained with the ALS training dataset.Only the most accurate parameterization values for all the methods applied at both resolutions (1 and 2 m) are analyzed.Mean Error (ME), Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and Standard Deviation of residuals (SD).Interpolation methods are ranked from lowest to highest RMSE.N = 1,896,468.

Table 4 .
Global statistics summarizing validation errors using the GPS benchmarks.Only the most accurate parameterization values for all the methods applied at both resolutions (1 m and 2 m) are analyzed.Mean Error (ME), Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and Standard Deviation of residuals (SD).Interpolation methods are ranked from lowest to highest RMSE.N = 55.