The Use of LiDAR Terrain Data in Characterizing Surface Roughness and Microtopography

The availability of light detection and ranging data (LiDAR) has resulted in a new era of landscape analysis. For example, improvements in LiDAR data resolution may make it possible to accurately model microtopography over a large geographic area; however, data resolution and processing costs versus resulting accuracy may be too costly. We examined two LiDAR datasets of differing resolutions, a low point density (0.714 points/m spacing) 1m DEM available statewide in Pennsylvania and a high point density (10.28 points/m spacing) 1m DEM research-grade DEM, and compared the calculated roughness between both resulting DEMs using standard deviation of slope, standard deviation of curvature, a pit fill index, and the difference between a smoothed splined surface and the original DEM.These results were then compared to field-surveyed plots and transects ofmicroterrain. Using both datasets, patterns of roughnesswere identified, whichwere associatedwith different landforms derived fromhydrogeomorphic features such as stream channels, gullies, and depressions. Lowland areas tended to have the highest roughness values for all methods, with other areas showing distinctive patterns of roughness values across metrics. However, our results suggest that the high-resolution research-grade LiDAR did not improve roughness modeling in comparison to the coarser statewide LiDAR. We conclude that resolution and initial point density may not be as important as the algorithm and methodology used to generate a LiDAR-derived DEM for roughness modeling purposes.


Introduction
Over the past several decades, geomorphologists, soil scientists, ecologists, foresters, and hydrologists have increasingly utilized terrain data for landscape classification [1][2][3][4], predicting forest communities [5], predicting soil properties [6][7][8][9], and understanding riparian zones and their stream networks [10].Due to improvements in data acquisition, computing power, and storage capacity, terrain data has become increasingly available at finer and finer resolutions and at broader scales, from National Elevation Dataset (NED) and Shuttle Radar Topography Mission (SRTM) to LiDAR.
Although LiDAR-derived DEMs have been shown to be extremely accurate when compared to non-LiDAR generated DEMs [11], the accuracy of LiDAR-derived DEMs for measuring landscape microtopography is debated [12].This can be due to data interpretation difficulties arising from abiotic (such as slope complexity) and biotic terrain factors (such as evergreen vegetation and coarse woody debris) [13,14].LiDAR processing methods may also result in terrain interpretation difficulties.For example, some researchers have found LiDAR-derived DEMs to be oversmoothed [12], which can minimize surface roughness and result in less topographic complexity.In contrast, others have found LiDARderived DEMs effective at identifying features such as landslides, which can have complex roughness patterns [15].
Microtopography is an important variable to measure for modeling water movement [16], geomorphology [17], vegetation dynamics [18,19], riparian communities [20,21], and surface roughness [22][23][24] or geologic features such as landslides and alluvial fan deposits [25,26].Modeling microtopography and its associated topographic variables may result in improved estimates of water storage and infiltration, a greater capability to identify wildlife habitat and facilitate high-resolution digital soil mapping.For example, pit-andmound topography is a type of microtopography very common in natural-forested ecosystems that is commonly driven by historical incidents of wind throw [27,28].Pit and mound topography is known to occur more frequently in areas prone to shallow rooting such as landscapes with a shallow soil depth to a water table or other restricting layer [29].Pit and mound topography may also result from differences in tree species, ages of stand, and other vegetation-based variables.The application of LiDAR data in identification of such landscapes may result in improved understanding of cooccurring soil and vegetation at site-specific to regional scales.
Various methods have been proposed to measure surface roughness including analysis of fractal dimensions, [30] identifying eigenvectors parallel to microsurfaces [15], and analyzing variograms associated with surfaces at multiple scales [31,32].Other methods include measuring the standard deviation or range of elevation over a particular scale, measuring variation of slope or curvature over multiple scales [33], and calculating variability over a small scale while removing the effect of the broader scale topography [25,26].Unfortunately, a consistent, preferred method of delineating surface roughness has not emerged given the diverse array of user needs.Scale is also an important consideration, as variation of slope or curvature calculated using a 10 m resolution DEM is measuring something very different than variation of slope or curvature calculated using a 1 m resolution DEM [34].
LiDAR data is becoming increasingly available for broad scale coverage of terrains.Often, this data is collected using a low initial point density (0.714 points/m 2 spacing) [35] when compared to research-grade datasets (10.28 points/m 2 spacing) [36].A goal of this research was to first evaluate two LiDAR datasets of differing initial point densities for relative accuracy using a field survey as a control.A second goal of this study was to evaluate roughness metrics derived using multiple methodologies to assess the effectiveness of both LiDAR datasets for characterizing surface roughness and microtopography.In addition, microtopographic signatures of several landforms within our research site are modeled and described.Surface roughness is modeled using a 1 m DEM; surface roughness and microtopography are used interchangeably.Although one of the questions investigated in this research is whether different soil units may exhibit different patterns of roughness, because this research was undertaken at a larger scale than available soil survey data, boundaries from the soil survey were not utilized to aggregate roughness data.

Study Area.
This study took place in a 119 ha watershed located in the Ridge and Valley Province of Pennsylvania, USA (Figure 1) known as Leading Ridge Watershed One.The elevation of the watershed ranges from 260 m at the mouth of the watershed to 512 m at the top of its northwestern border.
Since this watershed is in the Ridge and Valley Province, its hydrologic network tends to form a trellis pattern instead of the more common dendritic pattern typical of sandstone and shale bedrock [37].The geologic formation underlying the watersheds consists of deeply dipping strata ranging from resistant Tuscarora quartzite and sandstone at the top of the watershed to less resistant Rose Hill shale that comprises the valley area [37] (Figure 2).This terrain is consistent with that which makes up most side slopes and ridges of the Ridge and Valley Province, which is generally characterized by canoeshaped valleys and long, linear, parallel ridges formed by differential erosion [37].Ridges tend to be extremely steep and rocky, and the topography is generally well drained [37].This landscape has also been largely influenced by periglacial processes of the late Pleisticene [38].The watersheds contain a mature oak/hickory forest approximately 100 years of age.
This study site was chosen due to its history as an experimental watershed and the vast dataset that currently exists for the site.The Leading Ridge Watershed research units were established in Penn State's Stone Valley Experimental Forest of central Pennsylvania in 1959 as paired watersheds to study the hydrologic response of different forest practices [40].This watershed is also within the larger Shaver's Creek watershed which is being used as part of the Susquehanna/Shale Hills Critical Zone Observatory (CZO), an interdisciplinary observatory toward quantitatively predicting creation, evolution, and structure of regolith as a function of the geochemical, hydrologic, biologic, and geomorphologic processes operating in a temperate, forested landscape [41,42].
Digital U.S. Department of Agriculture Soil Survey Geographic (SSURGO) data [43] and field verification with soil pits and auger holes show that soil series in the upper half of the watershed are derived from sandstone colluvium and consist of the Hazelton and Dekalb series (Typic Dystrudepts).Sandstone and shale colluvium are the parent material for the Laidig series (Typic Fragiudult) found on the lower part of the watershed, the Buchanan series (Aquic Fragiudult) found across the valley bottom on the southwest portion of the watershed, and the Andover series (Typic Fragiaquult) across the valley bottom on the southeast portion of the watershed.Shale residuum is the parent material for the Berks (Typic Dystrudept) and Weikert series (Lithic Dystrudept) found on shale hills in the southern end of the watershed [43] (Figure 3).

DEM Data Sources.
Two sources of LiDAR-derived DEM data were utilized for this project.The first data set was collected in 2007 during leaf-off conditions as part of the Pennsylvania base map (PAMAP) LiDAR program.After spacing for the LiDAR returns used to generate the 1 m DEM was 1.4 meters, with a target vertical RMSE of 18.5 cm in open areas and 37 cm in vegetated or forested areas.Points were first classified as either ground or nonground points, with ground points being thinned to create a TIN that fit the final specifications by an independent vendor, BAE Systems.Using company-specific proprietary methods, an approximately 1 m resolution DEM was produced using the TIN.All finished products were checked for quality and accuracy and were  shown to meet or exceed target Pennsylvania state-defined RMSE requirements [35].
The second data set, CZO LiDAR, was collected in the winter of 2010-2011 by the National Center for Airborne Laser Mapping (NCALM).Initial LiDAR point density was approximately 10 points/m 2 , with bare earth point spacing of approximately 4 points/m 2 .Bare earth points were isolated using TerraScan and then converted to a DEM using Golden Software's Surfer 8 Kriging algorithm using a linear variogram model with a nugget variance of 0.15 m and a search radius of 25 m or 40 m.Complete specifications can be found in the 2010 NCALM project report [36].The final DEM resolution of this dataset was also 1 m.Final field-based RMSE measurements were not calculated for this DEM, but errors on initial point values were between 5 and 25 mm horizontally and 15-55 mm vertically.
Although the LiDAR-derived DEMs from the PAMAP program and the CZO LiDAR both have a 1 m resolution, many differences exist between the two datasets based on their initial point density and subsequent processing techniques.The PAMAP LiDAR was converted to a DEM using a process based on creating a triangular irregular network (TIN) from points classified as ground points using a proprietary algorithm.This TIN was then converted to a DEM.Conversely, the CZO LiDAR was converted to a DEM using a kriging technique.Figure 4 presents shaded relief maps created from each data set.The difference between the kriged and TIN-based DEMs is expressed in the faceted appearance of the shaded relief map of the PAMAP LiDAR dataset when compared to the CZO-generated dataset.

Modeling.
Surface roughness was assessed using four methods: standard deviation of slope, standard deviation of curvature, standard deviation of residual topography, and a pit fill metric.These were chosen based on their appropriateness for use at this scale and for geomorphometric analysis [33].ArcGIS (ESRI, Redlands, CA, USA) was utilized for all roughness modeling.Standard deviation of slope was calculated using focal statistics over a 5 m by 5 m moving window of a slope layer measured in percent slope.A 5 m by 5 m window was chosen in order to emphasize the fine scale of microtopography that meets the mapping objective for this study.This scale was chosen with consideration to the resolution of the DEM and the types of features  identified in the field that reflect microtopography of the site.Similarly, a second roughness metric was calculated using the standard deviation of curvature using a 5 m moving window.Curvature was calculated for the LiDAR-derived 1 m DEMs using the curvature tool that measures a combination of both plan and profile curvatures using the method of Zevenbergen and Thorne [34].Calculating curvature at this scale is representing the microscale curvature of the surface, not dominant surface features that would be traditionally identified as curvature.
A third roughness dataset was created by characterizing the difference between local elevation and residual topography [25,26,33].This was done by generating a new, smoothed surface using a thin-plated spline on a resampled 10 m DEM.First, the LiDAR-derived 1 DEM was thinned to 10 m resolution.A regularized spline with a weight of 0 was fitted to the thinned data, and the 10 m DEM was interpolated back to a 1 m DEM.The difference between the LiDAR-derived 1 m DEM and the resampled/splined DEM was calculated to show localized differences from the broader scale topography.
The fourth metric of microtopography was the pit fill metric, which measured the difference between a hydrologically corrected DEM (pits filled DEM) and the original DEM.This method was proposed because it may be able to identify pit and mound topography, vernal pools, and other features of high ecological significance.In order to calculate this layer, a filled DEM was created using the fill tool [44] with the LiDAR-derived 1 m DEMs.The original DEM was then subtracted from the filled DEM, and the values were summed over a 10 m area using the block statistics tool to improve visualization of data.

Field Verification.
Elevation values from the LiDARderived DEMs were validated in the field with a total station used to survey transects and microplots throughout the watershed (Figure 5).Four transects were located perpendicular to landform breaks throughout the watershed.Each transect was approximately 100 m in length, with one end being in one SSSURGO [43] soil map unit and landform, and the other end in a different soil map unit and landform.Across each transect, elevations were measured at every slope break (approximately 80 points per 100 m transect).Two microplots (10 m by 10 m squares) were chosen on each transect in which a topographic survey was conducted across all areas of slope change (Figure 5).

Statistical Analysis.
In order to test and compare the accuracy of the LiDAR datasets to actual elevations, root mean square difference (RMSD) between the DEM-modeled elevations and the surveyed elevations was calculated for the CZO LiDAR dataset and the PAMAP LiDAR dataset using the field survey points as control points.There were over 700   points that were surveyed from benchmarks.These points are the same as those used in field verification.Each point was converted to a 0.3 m raster cell using the mean value of points to assign a raster value.RMSD of the LiDAR-derived DEMs was calculated using the formula RMSD = sqrt((1/) * sum(( −   ) ∧ 2)), where  equals the number of cells/points,  is the surveyed value, and   is the LiDAR DEM value.Roughness metrics generated from the PAMAP LiDAR were analyzed by landform position to improve the understanding of the surface expressions of particular soil types.Landform position was identified by analyzing slope, slope position, soil characteristics, and microtopography patterns and was classified as being either Top of Ridge, Top Slope, Lower Slope, Valley Bottom, or Shale Hill.The different soils found in the watershed could be expected to have different surface characteristics that could be represented as microtopography.Landform position was chosen because it approximately represents both geology and soil, and soil boundaries are not available at a fine enough scale to differentiate the watershed.

Field-Based Modeling.
Root mean square differences from both datasets tended to mimic one another although there were some differences.The CZO LiDAR RMSD was slightly larger than the RMSD for the PAMAP LiDAR with values of 0.417 m and 0.410 m, respectively (Table 1).Differences between DEM elevations and surveyed elections ranged within about 1.5 meters for both datasets (Table 2).RMSD values for both datasets were greater on hillslopes (Figure 6) with slopes approaching 100%.RMSD values for both datasets were also greater directly along the stream channel.The mean difference between surveyed and LiDARderived DEM elevations was −0.3 meters (surveyed elevations were on average about a third of a meter lower than the DEM elevations).

Roughness
Modeling.Due to differences in processing algorithms of the two data sets and their respective DEMs, the algorithms used to calculate the roughness metrics generated roughness maps for each LiDAR product that appeared very different at a fine scale, although the broad patterns of high and low roughness were more consistent.In Figure 7, the pit fill metric is shown for the PAMAP and CZO DEMs.Both datasets presented similar patterns, with depressions occurring on both tops of ridges and along valley bottoms.Many of these depressions tend to be located along stream areas in both datasets, and there is a clear line across the middle of the slope where the prevalence of these depressions increases.There is also an inverse relationship between the rate of closed depressions and slope; areas on the watershed with a higher slope also have a lower pit fill metric value.
The second roughness metric analyzed was standard deviation of curvature over a 5 meter moving window (Figure 8).This metric resulted in a very clear striping artifact pattern approximately aligned to the dominant slope in the CZO DEM.When the actual curvature values and shaded relief maps were analyzed at a fine scale, it became apparent that this striping was caused by the initial DEM processing and its use of a kriging algorithm to generate the DEM from the LiDAR points.Despite this striping, broad scale patterns in roughness values were visible from this layer.The PAMAP LiDAR did not display such striping artifacts and showed similar broader patterns in roughness.Areas of high standard deviation of curvature are found along steeper rocky slopes and along linear features such as roads, stream channels, and slope breaks.Features perpendicular to the slope are prominent, particularly when compared to the standard deviation of slope layer (Figure 9).
Patterns in the standard deviation of slope (Figure 9) are similar to the patterns found in the standard deviation of curvature (Figure 8), with both methods producing very high values along linear features such as streams and roads.The last roughness metric (spline) was based on the degree that the local topography differed from the regional topography, and was calculated by thinning the DEM and creating a splined surface, then subtracting the original from the Applied and Environmental Soil Science splined surface and measuring the absolute value to remove negative numbers (Figure 10).Differences between local and regional topography are in part due to artifacts from each data sets' respective interpolation algorithm, particularly in the PAMAP LiDAR this time; however, in both data sets, there is similarity in landscape patterns with high and low roughness values.For example, high roughness values tend to correspond to areas of high slope, and an area of low roughness was highlighted along the top of Leading Ridge and in the mid-slope area.
Visual analysis of all roughness metrics together allowed for the identification of contrasting roughness patterns (Table 3); descriptive statistics calculated for each roughness metric per delineated landform and means are shown in Table 4 (these values were calculated for the PAMAP DEM).For example, the Top of Ridge position was characterized by having high pit fill metric values, low values of standard deviation of slope, low values of standard deviation of curvature, and low values for the difference between splined surface and regular surface.The Top Slope position exhibited an opposite pattern.The Lower Slope was characterized by large amount of variability in all of the roughness metrics, along with a high value for the pit fill metric.The Valley Bottom contained the highest values for all roughness metrics, while the Shale Hill contained relatively low roughness values for all of the metrics.Landforms in general have boundaries similar to Table 3: Qualitative interpretive table used for visual interpretation.Landforms were delineated using these measurements.For example, Top of Ridge was delineated based on a low value for SD of curvature, a high value for the pit fill metric, a low value for SD of slope, and a low value for the spline.

Discussion
We found that surveyed elevations were approximately 0.3 meters lower than the elevations measured by both LiDAR datasets, with the largest differences occurring in areas of high slope.This agrees with Spaete et al. [14] who found high LiDAR-derived DEM error rates in areas with high slopes.
The presence of vegetation in the form of tree roots, coarse woody debris, and evergreen vegetation could be the reason the surveyed elevations were slightly lower than the modeled elevations, as some LiDAR signals may have been reflecting off of tree roots or coarse woody debris (CWD) instead of the ground.Tenenbaum et al. [45] suggests that tree roots may affect DEM results by creating noise in the original point cloud data.This suggests that in forested settings, particularly densely forested settings such as eastern deciduous forests, there may be error in LiDAR-derived DEMs caused by vegetation in the understory, coarse woody debris and roots of trees, and leaf litter, even in leaf off conditions.In our study, RMSD results were similar across the watershed, but this could vary in different vegetation or landform settings.More research should be conducted to explore whether this may be due to topography or vegetation.Both datasets represented roughness differently, which may indicate that there is no clear advantage to researchgrade LiDAR for calculating roughness metrics unless accompanying resolution is also improved, such as moving to a 0.5 m resolution DEM.We suggest that for DEM generation and roughness calculations, initial point density is less important than the algorithm type used to generate the DEM, at least at the 1 m resolution.A previous study analyzing thinning of LiDAR ground points [46] found that initial LiDAR point data can be thinned by at least 50% with minimal effect to DEM accuracy, but that study was basing their results on density of ground points, while this study examined density of all points.Mitasova et al. [47] questioned the appropriateness of using kriging methods to generate DEMs from LiDAR point clouds, in part due to the high initial density of LiDAR point clouds, which is also reflected in the results of this study.The kriging algorithm definitely generated some periodic errors in the ground surface, as expressed by the standard deviation of curvature metric.This highlights that more research could be done to quantify the impact of algorithm on roughness metric generation and more care used in generating DEMs for different purposes.
By analyzing patterns of the different roughness metrics (Table 3), we can delineate physiographic features based on roughness.Figure 11 shows landforms with a graph of the various roughness values for each landform.For improved visualization, roughness metrics were converted to relative roughness indices by dividing each mean roughness by the highest mean roughness value for that metric.In Figure 11, Valley Bottom features have the highest roughness values for all four roughness metrics.Ridge Top features have relatively low roughness values except for the pit fill metric, which is high.The pit fill metric and the spline both indicate rougher surfaces in the Valley Bottom than in other sites, while the SD of curvature and the SD of slope are more similar across formations.The Top Slope area has higher relative roughness values for SD of slope and SD of curvature, while the Bottom Slope and Shale Hills are fairly similar.The Bottom Slope still has a relatively high value for pit fill metric, and the Shale Hill has almost no filled pits.
The Ridge Top feature has very high values in the pit fill roughness metric which measures the difference between the filled and original DEM.This may be due to terrain rockiness or an interaction between the LiDAR beam and the high density of ericaceous shrubs (e.g., blueberry, huckleberry, and laurel) which may have created error in the DEM.The ridge top soils are known to contain a high rock fragment content [48], which may be reflected in these results.Even though LiDAR data was collected during leaf-off deciduous conditions, dense branching patterns of deciduous shrubs could have resulted in erroneous pits being represented on the landscape due to false ground points that have a higher elevation.Some watershed areas, particularly rocky terrain and mid-slope areas of the watershed which have visually complex terrain patterns, did not have a high surface roughness as might be expected.In addition, talus areas in the watershed, with rocks approximately 0.5 m on a side, did not have high roughness values using any of the metrics.This suggests that other studies that have successfully modeled roughness [25,26] may have been detecting features on the ground such as landslides and alluvial fan deposits that are seen on the surface as larger features than the talus slopes in this study area.This may also indicate that the scale of the LiDAR, although extremely fine, is still not sufficient for predicting talus areas under a canopy cover.Another factor that may complicate roughness measurements at a finer scale is the effect of coarse woody debris on the LiDAR signal.The noise generated by forested terrain may be enough to obscure the signal of rocks and other small surface features.Features that have high roughness values using these roughness metrics are larger features such as major slope breaks, stream channels, roads and trails, and other man-made features.Textural features such as rockiness occur on too fine of a scale to be detected and measured using a 1 m DEM.Much of the modeled roughness on Leading Ridge seems to be from intermittent and ephemeral flow channels which are very dense in the lower portion of the watershed.We found differences in the boundaries of soil series and landforms associated with soil series boundaries when using LiDAR-derived DEMs, which reflects mapping scale differences between soil polygons as delineated from USDA-NRCS SSURGO data and our roughness metrics.However, by incorporating roughness metrics into our analysis, we may be able to refine our soil mapping polygons and thus improve the soil survey.For example, the presence or absence of a fragipan could be expressed on the landscape as an increase in local features created from increased pit and mound topography caused by increased windthrow due to more shallow rooting depth in these sites (due to a fragipan) [29].

Conclusions
When compared to a high-resolution ground survey, both the CZO DEM and the PAMAP DEM had a RMSD of approximately 0.4 m, with surveyed elevations being on average 0.3 m lower than DEM-modeled elevations.Therefore, research suggests that research-grade LiDAR was not any more accurate than the statewide LiDAR dataset.Additionally, we suggest that the high-resolution research-grade LiDAR did not improve roughness modeling in comparison to the coarser statewide LiDAR.No single roughness method stood out as the most effective at delineating physiographic landforms, but when viewed simultaneously, roughness patterns relating to the presence or absence of hydrogeomorphic features were visible from the data and associated with landforms.A question requiring further research is the effect that the algorithm used to generate a DEM can have on resulting patterns of surface roughness.

Figure 1 :
Figure 1: Site of the Leading Ridge Watersheds in Pennsylvania, USA.
T u s c a r o r a q u a r t z i t e C a s t a n e a s a n d s t o n e R o s e H i l l s h a l e K e e f e r S a n d s t o n e R o c h e s t e r s h a l

Figure 2 :
Figure 2: Map of the watershed showing approximate locations of contacts and the approximate arrangement of geology in the Leading Ridge Watershed 1 (adapted from shields 1966 [39]).Rock formations are almost vertical in the watershed.
loam, from 8 to 25 percent slopes Hazleton-Dekalb association, moderately steep Hazleton-Dekalb association, steep Hazleton-Dekalb extremely stony sandy loams, from 0 to 8 percent slopes Laidig extremely stony loam, from 8 to 30 percent slopes Soil map units in Leading Ridge Watershed 1 Microplot locations Transect locations Andover extremely stony loam, from 0 to 8 percent slopes Berks shaly silt loam, from 8 to 15 percent slopes Berks-Weikert association, steep

Figure 3 :
Figure 3: USDA-NRCS SSURGO mapping units for the Leading Ridge Watershed with transect and plot locations.

Figure 4 : 9 Figure 5 :
Figure 4: Shaded relief maps generated from PAMAP LiDAR (a) and CZO LiDAR (b).The TIN-based algorithm was used for the PAMAP LiDAR and is visible in the shaded relief map.Small trails and streams are clearly more visible in the CZO LiDAR image.

Figure 6 :
Figure 6: Close-up view of two plots and transect.Values represent the difference between surveyed elevations and LiDAR-derived elevations.Images are overlain over shaded relief maps of corresponding DEMs.Values in the key are in feet to match the original units of the DEMs.

Figure 8 :Figure 9 :
Figure 8: Standard deviation of curvature values for CZO LiDAR (a) and PAMAP LiDAR (b).White areas correspond to high roughness values.

Figure 10 :
Figure 10: Difference between splined surface and original DEM for CZO LiDAR (a) and PAMAP LiDAR (b).White areas correspond to high roughness values.

Figure 11 :
Figure 11: Roughness landforms with bar graphs showing relative value of roughness indices for each landform.

Table 1 :
RMSD for PAMAP LiDAR and CZO LiDAR in the Leading Ridge Watershed (in meters).

Table 2 :
Difference between surveyed points and the CZO and PAMAP LiDAR datasets (in meters).

Table 4 :
Mean values of roughness metrics by landform delineated by roughness metrics.