Palaseanu-Lovejoy, M.; Danielson, J.; Thatcher, C.; Foxgrover, A.; Barnard, P.; Brock, J., and Young, A., 2016. Automatic delineation of seacliff limits using lidar-derived high-resolution DEMs in Southern California. In: Brock, J.C.; Gesch, D.B.; Parrish, C.E.; Rogers, J.N., and Wright, C.W. (eds.), Advances in Topobathymetric Mapping, Models, and Applications. Journal of Coastal Research, Special Issue, No. 76, pp. 162–173. Coconut Creek (Florida), ISSN 0749-0208.
Seacliff erosion is a serious hazard with implications for coastal management and is often estimated using successive hand-digitized cliff tops or bases (toe) to assess cliff retreat. Even if efforts are made to standardize manual digitizing and eliminate subjectivity, the delineation of cliffs is time-consuming and depends on the analyst's interpretation. An automatic procedure is proposed to extract cliff edges from high-resolution lidar-derived bare-earth digital elevation models, generalized coastal shoreline vectors, and approximate measurements of distance between the shoreline and the cliff top. The method generates orthogonal transects and profiles with a minimum spacing equal to the digital elevation model resolution. The method also extracts the xyz coordinates for each profile for the cliff top and toe, as well as second major inflections along the profile. Over 75% of the automated cliff top points and 78% of the toe automated points are within 95% confidence interval of the hand-digitized top and toe lines, and over 79% of the digitized top points and 84% of the digitized toe points are within the 95% confidence interval of the automated top and toe lines along a stretch of coast in Del Mar, California. Outlier errors were caused by either the failure to remove all vegetation from the bare-earth digital elevation model or errors of interpretation. The automatic method was further applied between Point Conception and Los Angeles Harbor, California. This automatic method is repeatable, takes advantage of detailed topographic information within high-resolution digital elevation models, and is more efficient than hand-digitizing.
Seacliff erosion is a coastal hazard with impacts on coastal management, infrastructure, safety, and the local economy. In California, many coastal communities are located on top of uplifted marine terraces that are vulnerable to cliff-line retreat, which occurs episodically during large storm and wave events (Benumof et al., 2000; Earlie et al., 2015; Griggs and Johnson, 1979; Kuhn and Shepard, 1984; Sallenger et al., 2002). Modern erosion rates of coastal California cliffs over a 70-year period, based on aerial photos and more recent lidar data, averaged ~0.3 m/yr for southern California, ~0.5 m/yr for central California, and ~0.7 m/yr for northern California (Hapke and Reid, 2007). While the volumetric change and failure dynamics of the seacliff face can be captured directly by successive lidar scans (Collins and Sitar, 2008; Rosser et al., 2013; Sallenger et al., 2002; Young and Ashford, 2006), traditionally the recession of the cliff top or cliff base is obtained from hand-digitizing aerial photographs, topographic maps, or in situ surveys (e.g., Hapke and Reid, 2007; Moore and Griggs, 2002; Young et al., 2009). However cliff erosion is categorized, the positions of the cliff top and cliff base are key indicators of cliff behavior and are useful to coastal managers.
Hapke and Reid (2007) proposed a standardized procedure to hand-digitize cliff edges based on visually interpreting the break in slope from high-resolution hillshade renderings of lidar-based digital elevation models (DEMs). For very complex cliff morphology, such as where roads or terraces cut into an existing cliff slope or features associated with differential erosion of cliff-forming strata are present, interpreting the location of the cliff edge can be difficult (Hapke and Reid, 2007). Although efforts were made to standardize and eliminate as much as possible any digitizing subjectivity, the delineation of cliffs and other shoreline features is time-consuming and depends on the analyst's interpretation (Ruggiero et al., 2013). Rutzinger et al. (2012) demonstrated that breaklines manually digitized with the same digitizing procedure as described in Hapke and Reid (2007) vary in coverage and spatial accuracy, reflecting the analyst's skill, perception, and intention.
Automatic procedures are fully reproducible and comparable, resulting in comprehensive features derived independently of human skill. Since cliff edges are linear features detectable in both lidar data and lidar-derived DEMs, automated methods to extract breaklines can potentially be adapted to extract cliff edges. The automated methods of breakline extraction can be grouped into three major categories: deriving lines from intersecting planes, extracting lines through a neighborhood analysis of DEM elevation values, or applying edge detection filters and segmentation methods developed for image processing. Briese (2004) defined the breaklines as an intersection of a series of two planes, each on either side of the breakline. To obtain the final breakline, a spline function was fitted through the original intersection results. Brzank et al. (2005) improved Briese's (2004) method by using hyperbolic tangent functions recognizing that not all features have a planar representation, such as in the case of cliff faces. Brzank et al. (2008) developed a methodology for extracting structure lines using a combination of edge extraction methods applied to a lidar-derived DEM followed by fitting a pre-defined surface model described by a 2D version of the hyperbolic tangent curve. A derivation of the intersecting plane methodology is presented by Choung et al. (2013) in which angular values defined by two normal vectors on Delaunay triangles on both sides of a straight cliff edge are examined against a threshold chosen according to the local cliffs' characteristics. To obtain the final cliff edge all identified segments are merged.
Automated methods to extract linear features from lidar data based on curvature are investigated by Rutzinger et al. (2012) by looking at the influence of scale on the results. The automated results are compared with cliff tops digitized by nine different analysts. Manually digitized lines were found to have a low positional accuracy and tend to favor longer structures rather than shorter ones, in part due to a more limited information content of the shaded relief raster used as a base for digitization (Rutzinger et al., 2012). However, the automatic procedure was able to generate reproducible results, although by no means exhaustive, meaning that features that were not well defined in the DEM could be misinterpreted by the automatic procedure. While the curvature-based methodology uses the second order derivative of the surface elevations to identify significant linear features (Rutzinger et al., 2012), the Least Cost Path (LCP) analysis is based on cost functions (impedance) derived from DEMs in order to isolate features of interest. Although LCP may not be suitable in finding the top of the cliff since it is not always the highest elevation in the DEM and thus an impedance raster may not be able to isolate the top, it may be appropriate in defining the toe of the cliff if a modified version of the LCP is used that considers a more general distance from a stretch surface connecting the shoreline to the ridgeline (Hardin et al., 2012; Mitasova et al., 2011).
Processing the lidar-derived DEM as imagery and employing edge detection filters and segmentation procedures to extract breaklines was proven successful by several researchers. Sui (2002) employed a series of image processing methods to a grayscale representation of lidar-elevation data and applied edge detection techniques to extract structural lines. Richter et al. (2013) used a combination of Sobel and Laplace filters on a lidar-derived DEM, taking advantage of the broad edge detection of the Sobel filter and the exact edge positions of the Laplace raster derivative in order to determine the actual cliff position. Di et al. (2003) used the mean-shift-segmentation method for shoreline extraction from high-resolution IKONOS satellite imagery, while Lee et al. (2009) extended the same methodology to the integration of lidar data with color aerial orthoimages.
One method to automatically generate cliff top and toe edges that cannot be classified in any of the three categories of structure line detection mentioned above is the method developed by Liu et al. (2009). This method is based on elevation profile extraction across the cliffs and the observation that generally the variation of the slope along the elevation profile is commonly greater at the top and the toe of the cliff than anywhere else along the profile. However, this observation may not be substantiated in the case of complex cliffs with roads or terraces cut through the cliff gradient, cliffs with different erosional profiles or slope gradients, or cliffs formed at the base of hills.
The automatic methods described above can in principle be used to define cliff edges, but none of them is fully suitable in dealing with complex cliffs in very different geomorphic settings, and they are not always location independent. In order to accommodate both simple and complex cliffs, this study presents an automatic procedure based on profile extraction from lidar-derived bare-earth DEMs that does not involve variation in slopes between the profile points, is robust in different geomorphic settings, can determine secondary points of inflection along the profile besides the top and the toe of the cliff, and is location independent. This new procedure takes advantage of the full resolution of the lidar-derived bare-earth DEM, and all parameters involved are linked to the DEM resolution and accuracy.
The automatic methodology was developed on a 2.5-km stretch of coast in the Del Mar area of southern California and refined further by being applied on a more morphologically complex coast between Point Conception and Los Angeles Harbor, California (Figure 1).
The cliffs in the Del Mar area have less morphological complexity with a mean height of about 15–20 m and an average slope of 40–45 degrees, with a maximum slope of slightly over 70 degrees. The cliffs between Point Conception and Los Angeles Harbor are 10 to over 50 m in height, with few over 100 m, and are locally either terraced or intersected by roads and can be almost vertical in some parts (Figure 2).
A comparison of summary statistics of cliff heights and slope from the two sites is presented as a series of violin plots in Figure 3. A violin plot is essentially a box plot that incorporates information about the underlying distribution of the data, adding a rotated kernel density plot on each side (Hintze and Nelson, 1998).
The cliffs in these areas are mostly composed of Miocene and Pliocene marine and nonmarine sedimentary rocks (sandstone, shale, siltstone, and conglomerate, mostly moderately consolidated) and Pleistocene-Holocene surficial terrace deposits that are unconformable on the Miocene-Pliocene formations, although some gradational contact was noticed in places (Griggs and Savoy, 1985; Gutierrez et al., 2010; Kennedy, 1975). The terrace deposits are mostly unconsolidated detrital material (alluvium, lacustrine, and playa deposits) and the gradational Pliocene-Pleistocene deposits are mostly loosely consolidated sandstone, shale, and gravel. The high coastal bluffs and vertical seacliffs are cut on the sides of high hills and are divided by canyons and valleys, sometimes backing small pocket beaches (permanent or seasonal) separated by rocky headlands, or can linearly stretch for kilometers.
Wave erosion, inundation, and rainstorms are major contributors to cliff retreat, but human activity such as construction of breakwaters and other engineering structures can significantly impact wave patterns and alongshore sediment flux that can in turn exacerbate cliff and bluff erosion. Episodic collapse of sea caves can produce nearly instantaneous retreat of cliffs up to tens of meters (Griggs and Savoy, 1985).
The challenge was to define an automatic procedure that performed adequately on any type of cliff morphology. The automatic procedure of delineating the cliff top and toe was first developed on the less morphologically complex area of Del Mar. Carefully hand-digitized limits on the same dataset used to develop this automatic procedure were available for comparison (Young et al., 2009). In contrast, the area between Point Conception to Los Angeles Harbor was much larger (more than 100 times greater in length than the Del Mar site), with a complex geomorphology of pocket beaches and rocky promontories, interspaced with stretches of linear cliffs.
The data for the Del Mar area were collected using an Optech Inc. Airborne Laser Terrain Mapper (ALTM) 1225. The survey's metadata reported a vertical accuracy of 19.6 cm (root mean square error [RMSE] of 10 cm) and a horizontal accuracy of 1 m (SCRIPPS, 2015).
The data used for the Point Conception to Los Angeles Harbor cliff analysis were collected in 2010 by TerraPoint using an ALTM lidar scanner. The reported vertical RMSE for this survey was 9.25 cm (Dewberry, 2011) and the estimated horizontal RMSE was 0.5 m (TerraPoint, 2004).
The cliff edge automated extraction procedure is based on a bare-earth DEM obtained from high-resolution lidar data and a generalized shoreline vector that was used as a reference feature. The results are dependent on the quality of the bare-earth DEM and the lidar data classification, including how well vegetation was removed from the bare earth surface. The shoreline should be free of tight bends and as much as possible parallel with the general direction of the cliffs. Furthermore, the shoreline was split in segments up to 2,000 m using the split polyline tool from XTools Pro 11 in ArcGIS 10.x (Data East Soft, 2015), and for each segment, an approximate distance in meters that was slightly larger than the distance between the shoreline and the approximate top of the cliff was recorded in the shapefile attribute table. Generally, this distance can fluctuate a few tens of meters without changing the delineation of the position of the top and the toe of the cliff, as is demonstrated later in the paper, and it is the single un-automated step. This distance was used to generate buffers around the shoreline segments, and these polygons were used as masks to clip the bare-earth DEM for each shoreline segment. An ArcPy script (Esri, Redlands, CA) was utilized to automate the generation of separate shapefiles for each shoreline segment, buffer polygons, and the clipped DEMs, and to perform basic file management. The buffer should be large enough to include the cliff and small enough so that the cliff is the only major geomorphic feature included in the DEM.
All these data (shoreline segments, buffers, and clipped DEMs) were loaded into a script written in R, an open source statistical programing environment (R Development Core Team, 2014). Custom functions were written to generate pseudo-perpendicular lines that intersect the cliff. The buffer was generated with rounded ends (the default option in ArcGIS), and the lines that connect these ends were called the buffer margins (Figure 4a, lines in magenta and cyan). Points placed at a user-defined distance apart (usually equal with the DEM resolution) were generated along the longest buffer margin, and the same number of equally spaced points was established on the shorter margin. The next step was to identify for each point on the longer buffer margin the corresponding point from the shortest buffer margin that was the minimum distance away. This could result in a many-to-one relationship for some of the points, and generate transects too far apart along the main direction of the cliff (Figure 4b, transects in thin magenta line). In order to correct for that, for all the points not selected on the short buffer margin, minimum distance points on the longer margin that were between two already consecutive transects were identified (Figure 4c, transects in thin cyan line). This was to ensure that transects would never intersect across the cliff itself and they would be pseudo-perpendicular to the cliff direction. If the many-to-one situation generated transects that were too dense across the cliff direction, then the final step required was to select only one in 2 or 3 existing transects as final. The resulting cliff top (cyan line) and toe (magenta line) are shown in Figure 4d.
Furthermore, xyz coordinates were extracted along these pseudo-perpendicular transects to generate profiles spaced apart at a distance approximately equal to the resolution of the DEM. Each profile was analyzed to resolve the top and toe of the cliff and to determine the second major inflection along the cliff face. For any given profile, a straight line connects the first point on the profile with the last point. Each profile was reoriented in such a way that the first point on the left of the profile was always toward the ocean and the last point on the profile to the right was always toward land. Next, the distance from each point on the profile to this straight line was calculated following the convention that all points that were above this line have a positive distance and all points below this line have a negative distance. The top of the cliff always has the largest positive distance to the straight line that connects the first and last point of the profile, while the toe of the cliff has the smallest negative distance to that same line (Figure 5).
The method is robust enough that the position of the top and toe of the cliff does not change with the width of the buffer/length of the profile as long as the cliff was the most prominent geomorphic feature present (Figure 6).
The same procedure was applied again between the identified top and the toe of the cliff. The point of inflection with the largest positive distance to the straight line between the top and the toe of the cliff represents the second major positive inflection. This inflection point is important when the cliff profiles are more complex, cut by narrow roads, terraced for stability, or have structural terraces due to differential rates of erosion.
All the xyz coordinates of the profiles were saved as ASCII files for future analysis, and the top, toe, and secondary inflections were saved as points and line shapefiles. Since all the information was recorded for each profile across the cliffs, further analysis of cliff face geometry (concavity/convexity) and generalized cliff slopes could be computed. An overall procedure diagram is presented in Figure 7.
The concavity/convexity of the cliff face was determined against the straight line between the top and the toe of the cliff. By convention, if the cliff profile was above this straight line (positive distance) then that portion of the cliff face was considered convex, and if it was below the straight line (negative distance) it was concave. The point on the convex profile that is farthest away from the straight line corresponds to the previously identified second major positive inflection point. A cliff profile can be entirely convex or concave, or both if the cliff profile intersects the straight line between the top and the toe. Furthermore, for each profile, summary statistics (minimum, first quartile, mean, median, third quartile, maximum, and standard deviation) of these distances were calculated. These values can provide relative information on how pronounced the concavity or convexity is, or if the cliff profile is more convex than concave overall, or the reverse. In order to visualize more easily these data and more quickly identify areas of interest (e.g., profiles that are particularly concave or convex), a heat map type of raster was generated in which the individual values contained in the matrix were represented as colors (Kinney, 1991). In this case, the raster x-axis signifies the cliff profile's identification number that can be associated with the respective profile location, and on the y-axis the summary statistics values were symbolized by two different contrasting color scales, one for concave values (yellow-green to dark blue) and one for convex values (light yellow to dark brown, Figure 8). The darker the color on both scales, the more pronounced the concavity or convexity. Profiles with darker contrasting colors representing the summary statistics along the profile indicate locations of potential instability, while light colors reflect a state of equilibrium in which the cliff face profile closely approximates a straight line.
RESULTS AND DISCUSSION
The automated procedure saves point and line shapefiles for the top and the toe of the cliffs, point shapefiles for secondary inflections along the profiles together with a classifier that identifies if the secondary inflection is above a certain threshold and thus significant, and ASCII tables of each profile xyz coordinates. The cliff top and toe data were further processed by applying an outlier filter to eliminate some of the inconsistencies due to DEM errors, remnant vegetation, or sudden transitions between adjacent shorter and higher cliffs. This type of filter does not smooth the original result but rather eliminates outliers using the box-and-whisker method (Tukey, 1997) in an overlapping window (Figure 9a, b).
Comparison of Automatic Results with Hand-digitized Results from the Same DEM Dataset
The cliffs in the Del Mar area were previously studied by Young et al., (2009), who very carefully hand-digitized the cliff top and toe using the same 2006 lidar-derived DEM used in this study. Also, cliff tops were delineated for this area as part of the USGS National Assessment of Shoreline Change project based on lidar-derived DEMs surveys from 1998–2002 (Hapke and Reid, 2007). Due to the regional to national scale of the National Assessment data and the fact that cliff toe data are absent from this dataset, only a visual comparison is presented (Figure 10).
The comparison between the hand-digitized and automatic data is done by either calculating the shortest distance between the digitized points for both top and toe to the respective automated top and toe lines, or calculating the shortest distance between the automated top and toe points to the respective digitized lines. Both analyses are necessary because the hand-digitized top has 599 points and the toe 716 points, while the automated top has 2,403 points and toe 2,461 points. The automated procedure generated an equal number of points (2,486 points) for top and toe, but the digitized top and toe lines are shorter, and the top lines are not continuous, so all the automated points that do not correspond to a digitized line were eliminated. The results are presented in Figure 11 as violin plots and the errors and statistics of minimum distance between the automated top and toe cliff points and their respective digitized line in Table 1.
Error budget and statistics of minimum distance between automated top and toe cliff points and respective digitized lines; MAE = Mean Absolute Error (m), Median = Median (m) RMSE = Root Mean Square Error (m), SEM = Standard Error of the Mean (m), St.Dev = Standard Deviation (m), Min = Minimum value (m), Max = Maximum value (m), % <= 1 m = Percent of points within horizontal accuracy of 1 m.
Considering the Del Mar DEM horizontal accuracy of 1 m (Scripps, 2008), distances between hand-digitized and automatic points within this interval cannot be considered statistically different. Consequently, 75.41% of automated top points and 78.18% of automated toe points are within 1 m of the respective top and toe digitized lines. Looking at the outliers defined by the box-plots (Figure 11), 7.0% of the cliff top points and 3.6% of the toe points are in this category. The cliff top outliers range between 2 to 14 m and the distribution has a major break at 5.3 m, whereas the toe outliers' range is between 1.9 and 7.6 m with a major break in the outliers' distribution at 3.2 m. All the top outliers above the 5.3 m threshold (1.1% of all top points) were due to tall vegetation that was not completely removed from the bare earth DEM. One exception is a location where the cliff profile is relatively complex and the digitized line identifies a secondary inflection on the profile instead of the cliff top (Figure 12a, b).
The presence of the extreme outliers translated in a root mean square error (RMSE) and standard deviation greater than 1 m (1.38 m and 1.17 m, respectively, Table 1); however, the mean absolute error (MAE) and the standard error of the mean (SEM) were less than 1 m (0.53 m and 0.03 m, respectively, Table 1). If the extreme outliers greater than 5.3 m were eliminated, then the RMSE would drop from 1.38 m to 1.04 m. The rest of the top outliers (5.9%) are due generally to more complex profiles where the digitized top identifies a secondary inflection on the profile instead of the top, and to some extent to digitizing errors or incomplete removal of vegetation from the DEM (Figure 13a, b, c).
There were only 7 toe outliers (0.3%) with distances greater than 3.2 m. The most extreme outlier, with a distance to the digitized toe line of 7.6 m, was due to the presence of landslide talus at the base of the cliff. The digitized line is in front of the talus while the automated toe point is located between the talus and the cliff, indicating two different viewpoints, one that includes the talus in the active cliff profile and one that does not. The automated procedure identifies the end of the talus as a secondary toe inflection, coinciding with the digitized line (Figure 14).
The rest of the extreme outliers are due to more complex profiles where either vegetation, landslide talus, or digitizing errors are present (Figure 15).
Except in the case of obvious digitizing errors, the digitized line is identified by the automated procedure as a secondary toe inflection. The rest of the toe outliers (3.3%) are due to potentially small landslide debris that are sometimes identified by the automated procedure as secondary toe inflections, depending on their size, and potential errors in the digitized toe line (Figure 16a, b). Even with the extreme outliers, all the automated toe error statistics were below 1 m, with an RMSE of 0.77 m, MAE of 0.34 m, and SEM of 0.01 m (Table 1).
The minimum distance between the digitized top and toe points to the automated top or toe lines was also compared and the error statistics are presented in Table 2. The number of digitized points within 1 m of the automated top or toe lines increased by 4 to 5% compared to the previous calculation, because only the digitized vertices were considered, not the digitized line itself. The cliff top and toe were digitized in less detail with 3 to 4 times fewer vertices than the automated method, and the line connecting the digitized vertices did not always follow the sinuous shape of the cliff as closely as the automated method. For these reasons, the digitized top and toe points were closer to the automated lines than the automated vertices to the digitized line. For both digitized top and toe points, the error statistics were below 1 m (Table 2).
Error budget and statistics of minimum distance between digitized top and toe cliff points and respective automated lines; MAE = Mean Absolute Error (m), Median = Median (m) RMSE = Root Mean Square Error (m), SEM = Standard Error of the Mean (m), St.Dev = Standard Deviation (m), Min = Minimum value (m), Max = Maximum value (m), % <= 1 m = Percent of points within horizontal accuracy of 1 m.
This paper presents an automated procedure to extract cliff edges (top and toe) from a high-resolution lidar-derived bare-earth DEM. This automated procedure is flexible enough to accommodate the highly complex geomorphology of the southern California coast, adjusting from almost vertical cliffs with sharply defined top and toe inflection points to complex cliff profiles. Within the southern California study area, cliff complexity is related to roads and terraces cut into the cliff gradient, cliffs with different slope gradients or differential erosional profiles, or cliffs formed at the base of hillslopes. This procedure extracts transverse profiles across the cliffs spaced at any distance interval specified by the user (up to the DEM resolution). Besides the top and the toe of the cliff, the procedure defines other major inflections along the profile that may correspond to points of major change in general slope, structural terraces due to differential erosion, or landslides edges. The results are repeatable with high precision, comprehensive, and robust regarding analyst decisions on parameters. However, the results of this automated procedure depend on the quality of the DEM used, and especially how well the vegetation was removed. The threshold parameters used in the automated procedure are location independent, but the threshold parameters are dependent on the accuracy of the DEM. Even if, generally speaking, the manually digitized cliff edges have a low location precision, and are not reproducible, over 75% of automated generated cliff edge points are within 95% accuracy interval of the manually digitized lines. Outliers are due to tall vegetation not removed from the DEM, complex cliff profiles in which the cliff edges were confused with secondary inflection points during manual digitizing, landslide debris, or simply errors of interpretation, and are less than 10% of the total points. Incorporated in the automated procedure are routines that will assess the convexity/concavity of the cliff face, which can indicate instability hotspots along the cliff shoreline.
Future work will focus on rigorously validating the automated cliff edge results using field data collected with real-time kinematic (RTK) global positioning system (GPS) along the cliff top and toe, and a Riegl VZ-1000 Terrestrial Laser Scanner. Scans will be conducted from multiple locations on the cliff tops and beach in order to provide complete coverage of cliffs tops and toes, cliff faces, and the beach itself. The terrestrial lidar data will focus on morphologically complex cliffs with overhangs and sea caves that are not resolved by aerial lidar data.
We appreciate reviewers' constructive comments. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.