Multitemporal sets of lidar data provide a unique opportunity to analyze and quantify changes in topography in rapidly evolving landscapes. Methodology for geospatial analyses of lidar data time series was developed to investigate patterns of coastal terrain evolution, including the beach and dune systems. The diverse lidar-point data density, noise, and systematic errors were first quantified, and the results were used to compute a consistent series of high-resolution digital elevation models using spline-based approximation with optimized parameters. Raster-based statistical analysis was applied to the elevation-model time series to derive maps representing multiyear trends in spatial patterns of elevation change, to quantify dynamics at each cell using standard deviation maps, and to extract the core surface below which the elevation has never decreased. The methodology was applied to a North Carolina barrier island that was mapped by a sequence of 13 lidar surveys during the past decade, using several different lidar systems. Assessment of vertical differences between the lidar data sets using stable structures such as a road, was shown to be essential for correct quantification of coastal terrain change and its pattern. The analysis revealed the highly dynamic nature of foredunes, the trend toward inland sand transport, and the impact of anthropogenic sand disposal on that trend.
INTRODUCTION
After more than 10 years of lidar mapping, the first decadal time series of elevation data are becoming available for sections of the U.S. Atlantic coast. They provide a unique opportunity to analyze spatial patterns of coastal terrain dynamics and to detect locations where the trends in terrain evolution may lead to breakdown of the landform system—an increasingly important task given global climate change and projected sea-level rise. Previous studies have demonstrated advantages of lidar surveys for evaluating beach changes and assessment of shoreline erosion (Sallenger et al., 2003; Stockdon et al., 2002). Lidar-based digital elevation models (DEMs) are also widely used for quantification of beach and dunes volume change (e.g., Mitasova et al., 2004; Overton et al., 2006; White and Wang, 2003; Woolard and Colby, 2002), especially to study impacts of major storms or hurricanes (Sallenger et al., 2006). The high density of lidar data points and almost-annual frequency of mapping offer the possibility of extracting additional information about spatial patterns of coastal dynamics using raster-based techniques. However, such analysis, when applied to decadal lidar time series, needs to address the issues of accurate data integration because of the rapid evolution of lidar technology during the past decade produced data sets with different accuracies, scanning patterns, and point densities. To address these issues, this article presents a geospatial-analysis method for a time series of lidar data that integrates data from different lidar sensors and derives new types of maps for assessing coastal change.
METHODS
The proposed methodology for analysis of coastal topographic change from diverse, multitemporal lidar data was developed using open-source geospatial tools implemented within GRASS GIS (Neteler and Mitasova, 2008). The methodology employs the following workflow:
the point-cloud data acquired from various sources are integrated within a single coordinate system;
the per-cell statistical analysis of point data is performed at a hierarchical set of resolutions, and the results are used to select common DEM resolution;
the spatial extent of each survey and a mask for the study area are derived from preliminary, lower-resolution DEMs computed using the mean elevation value for each cell;
more detailed, smoothed DEMs are computed for the masked study area using spatial approximation;
the DEMs are compared with high accuracy ground-based data to remove potential systematic errors and verify each DEM accuracy; and
multitemporal per-cell statistics are applied to the corrected time series of DEMs to derive maps that characterize various spatial aspects of terrain change.
The methodology was developed using a multitemporal set of lidar data that covers a narrow, approximately 700-m-wide and 15-km-long stretch of coast in the U.S. Fish and Wildlife Service Pea Island National Wildlife Refuge, along the Outer Banks of North Carolina (Figure 1). This section of a barrier island was repeatedly mapped by lidar over the past decade (Table 1). It was chosen because of the lack of civil infrastructure seaward of the oceanfront, two-lane road, North Carolina State Highway 12 (NC12), and the observed changes are limited to sandy surface.
Table 1.
Characteristics of the lidar surveys based on the available metadata.*
Lidar and Real-Time, Kinematic Global Positioning System Surveys
The lidar data used to demonstrate the methodology were acquired by several agencies for a variety of mission objectives: the U.S. Army Corps of Engineers Topo/Bathy mapping project in years 2004 and 2005; the National Oceanic and Atmospheric Administration (NOAA) lidar project in the summer of 2003; the National Aeronautics and Space Administration (NASA)/U.S. Geological Survey (USGS) Experimental Advanced Airborne Research Lidar (EAARL) Hurricane Isabel survey in the fall of 2003; the North Carolina Flood-plain Mapping Program (NCFMP) in year 2001; and the NOAA/NASA/USGS Airborne Lidar Assessment of Coastal Erosion project surveys carried out between 1996 and 1999. The point data were downloaded from two online distribution sites (NCFMP, 2006; NOAA Coastal Services Center, 2006) in the North Carolina State Plane North American Datum 83 coordinate system with North American Vertical Datum 88 and units in U.S. feet. Further details pertaining to each of the lidar surveys are found in Table 1, and metadata records are available for download from the distribution Web sites.
To verify the published accuracy and consistency between the lidar data sets, a real-time, kinematic (RTK) global positioning system (GPS) survey of a 15-km section of NC12 (Figure 1) was performed in July 2006. The survey employed a Trimble 5700 RTK GPS, mounted on the side of the survey vehicle. A detailed site-calibration was performed using a network of National Geodetic Survey, North Carolina Geodetic Survey (NCGS), and North Carolina Department of Transportation (NCDOT) benchmarks. Elevation of the highway centerline was surveyed in three passes; each road curb was measured in a single pass. Accuracy of the RTK GPS data was evaluated, using a set of 11 NCGS and NCDOT benchmarks located along the NC12 centerline, by computing the differences between the elevations of the RTK GPS points and the benchmarks. The RTK GPS points located within a 2-m distance from benchmarks were used for the evaluation. The root mean square error was 0.023 m, with a mean difference of −0.027 m, indicating a small systematic error. The accuracy of the RTK GPS survey along the stable road was sufficiently higher than the published accuracy of the lidar surveys (Table 1), making the RTK GPS data suitable for assessment of the actual lidar data accuracy.
Computation of Multitemporal Set of DEMs
To evaluate the properties of the point data sets obtained by differing lidar technologies, per-cell point-data statistical analysis (Bowman, 2007; Neteler and Mitasova, 2008) was performed at a sequence of resolutions: 0.30, 0.91, 1.83, 4.57, and 9.14 m (1, 3, 6, 15, 30 ft). Maps representing the number of points within each grid cell and maps representing within-cell point-elevation ranges were then used to select the resolution for the multitemporal set of DEMs. The resolution was based on the condition that the mean within-cell elevation range was smaller than the published vertical accuracy of the lidar data to ensure that important topographic features, such as the ridges of foredunes were preserved. This analysis was also helpful for selection of an approximation method that produced consistent elevation surfaces by evaluating the need for smoothing and interpolation. A series of preliminary DEMs for rapid visual analysis of available data was created using the mean elevation value for each grid cell at 4.57 m (15 ft) resolution. Those DEMs were also used to verify that the selected study area had complete coverage by all data sets and to create a mask for further processing and analyses.
A more detailed, smoothed set of DEMs and topographic parameters were computed using the regularized spline with tension (RST) method (Mitasova, Mitas, and Harmon, 2005). The tension and smoothing parameters for each data set were optimized by minimizing the deviation from the RTK GPS data while reducing the impact of noise and ensuring a comparable level of detail in each DEM (see Mitasova, Mitas, and Harmon, 2005a; Mitasova et al., 2004; or Neteler and Mitasova, 2008 for more details on optimizing the RST parameters for lidar data).
The interpolated DEMs were compared with high-accuracy RTK GPS data to identify and remove potential systematic errors that could distort the analysis of terrain change, to verify the published data accuracy, and to ensure that interpolation had not introduced unacceptable errors. If systematic error was detected along the road centerline, a correction was applied to the entire DEM, assuming that the systematic error was uniformly distributed. This assumption was based on the study area having relatively uniform land cover, composed of sand, pavement, and prevailing low slopes.
Terrain Evolution Analysis
To characterize coastal terrain evolution, a set of standard and novel coastal state indicators was defined. The standard indicators included shoreline and sand-volume change, with summary values and trends estimated for beach and foredune sections defined by cross-shore profiles (Overton et al., 2006).
This article focuses on the analysis of spatial patterns of terrain-surface evolution based on multitemporal, per-cell raster statistics (Wegmann and Clements, 2004; Neteler and Mitasova, 2008). The approach generates a set of maps, where each output cell value is a function of the values in the corresponding cells in the entire time series. Depending on the application, various functions can be used, from simple statistics, such as mean, median, mode, minimum, maximum, standard deviation, and diversity, to more complex relationships, such as linear regression slope, offset, and coefficient of determination, computed for each cell. The temporal aspect of terrain evolution can be analyzed using a map that represents the time when the given cell was at its lowest elevation and a map representing the time when each cell was at its highest elevation.
Terrain evolution can also be represented as a volume, with time used as third dimension, and elevation as a grid-cell value. Shoreline evolution is then derived as an isosurface with the elevation value equal to the mean high-water line. Isosurfaces extracted at gradually increasing elevations facilitate identification of elevation levels that are highly dynamic and those that are stable over time.
RESULTS
Multitemporal Data Integration and Systematic Errors
Per-cell statistical analysis of lidar point data confirmed that the mean point densities ranged between the low of 1 point per 3-m grid cell for the 1996 and 2001 surveys to a high of 1 point per 0.3-m grid cell for the 2004 survey (Table 1). The data with the highest point density significantly oversampled most terrain features but provided excellent representation of sharp edges and breaklines. The within-cell mean range of elevation data was more than 1 m at 9.14-m resolution, 0.2–0.6 m at 4.57-m resolution, and 0.3–0.02 m at 1.83-m resolution. At 1.83-m resolution, the within-cell mean range was less than the published data accuracy for all but the 2004 data set; therefore, the 1.83-m resolution was selected for the multitemporal set of DEMs. At this resolution, some detail captured by the 2004 data was lost, and interpolation was needed to fill in gaps in the rest of the DEMs.
Evaluation of data accuracy, especially identification of possible systematic errors, was essential. Comparison of the interpolated DEMs with the on-ground RTK GPS survey revealed vertical shifts in the lidar DEMs (Figure 2). The largest shifts (mean difference) of −0.23 m and −0.25 m were found for the February 2001 and September 16, 2003, surveys, whereas the September 18, 1999, and June 22, 2003, surveys were very accurate and had less than 0.01-m mean difference. The shifts were verified by computing elevation differences among the closest pairs of RTK GPS points and the original lidar points to eliminate potential contribution of interpolation to the shift and by displaying the DEM profiles along the RTK GPS road centerline data (Figure 2). Although the vertical shifts were mostly within the published lidar data accuracy, they had the potential to significantly influence the quantification of terrain change, especially when the change was computed using the data shifted in the opposite direction (e.g., one data set floating above the elevation surface, whereas the other one pushed below). The shifts were treated as a uniform systematic error, and the median value of the shift (Figure 2b) was applied to the relevant DEMs. Although the difference between the mean and median was relatively small (less than 0.04 m for all data sets, except for September 7, 1998, when it was 0.10 m), median was considered more appropriate because of its lower sensitivity to outliers present in some of the data sets. The improvement gained by the correction was verified by computing the differences between the DEMs (original and corrected) and the NCDOT benchmarks (Figure 3).
Figure 2.
Evaluation of systematic error using the NC12 highway centerline: (a) RTK GPS and DEM profiles for the September 18, 1999, and February 2001 surveys, with 2001 DEM floating below the RTK GPS profile; arrows show locations of the first and last NCDOT benchmarks used to verify the RTK GPS survey accuracy and (b) mean, median, and interquartile range of elevation differences between the DEMs and the RTK GPS survey.
Spatial Pattern of Terrain Dynamics
Multitemporal, per-cell statistical analysis resulted in a set of novel maps that highlight various aspects of terrain dynamics in the studied area bordered by the NC12 highway and the 2005 mean high-water line. Minimum and maximum elevation surfaces provided information about the three-dimensional envelope, within which, the elevation surface has evolved over the studied period (Figure 4a). The maximum elevation ranged between the low of 8.90 m (29.19 ft) and the high of 12.05 m (39.55 ft), with a mean elevation of 2.15 m (7.05 ft) for the entire period and for the entire studied section. The minimum elevation surface represents the core of the island's studied section present in all data sets (Figure 4a); its spatial extent was 349,100 m2, only 67% of the 520,000-m2 spatial extent covered by the maximum elevation surface. The volume of the core was 46% of the volume of the maximum elevation surface, 67% of the latest 2005 surface, and 78% of the September 9, 1999, volume (after hurricane Dennis, the smallest volume observed during the past 10 y). The standard deviation map with the mean value of 0.40 m (1.30 ft) identified the most dynamic areas in terms of elevation change, with the highest values on foredunes, in the area of the 2003 overwash (Overton et al., 2006) and around the narrow central sections (Figure 4b). The dynamics of the topography was reflected in 92% of the area having a standard deviation greater than 0.15 m (above the published accuracy of the lidar data).
Figure 4.
Raster time-series analysis: (a) cross-sections through the elevation-surface envelope defined by the core surface (red) and outer envelope (green) and (b) standard deviations map draped on the 2005 elevation surface with the 2005 shoreline.
Temporal aspects of terrain evolution were captured by a map representing the time when each grid cell was at its lowest elevation (Figure 5a) and a map representing the time when each grid cell was at its highest elevation (Figure 5b). Comparison of these two maps indicates that most of the beach and road side adjacent to the dunes had the highest elevation in 2005, whereas the dune ridges were highest at the beginning of study period in 1997 and 1998. The over-wash area in the southern section of the study area was lowest in 2003, after hurricane Isabel. The maps provide spatial information that complements the total volume estimates (Overton et al., 2006).
Figure 5.
Maps representing (a) time of minimum elevation, (b) time of maximum elevation, and (c) linear-regression slope map with elevation-increase trends in blue areas and elevation-decrease trends in red areas. Inset adds transparency as a function of r2 to the regression-slope map (areas with high r2 are transparent and show solid colors, whereas areas with low r2 are washed out, leading to white color in areas with r2 < 0.3).
The per-cell linear-regression slope and coefficient of determination r2 maps were used for detailed identification of areas with different trends in elevation change (Figure 5c). The resulting maps show that the ocean side of the foredunes has a trend of losing elevation, whereas the road side of the dunes is gaining elevation, especially in the northern section of the study area, indicating the potential impact of wind transport. The coefficient of determination r2 is highest on the dunes, ranging between 0.50–0.98. The beach does not exhibit a clear trend, as reflected by the low values of r2 < 0.3 because of the impact of sand disposal in 2003 that compensated for the loss in the late 1990s. The central section of the study area shows a trend toward losing elevation on both the ocean and road sides of the dunes, indicating high vulnerability.
The elevation isosurfaces derived from volume representation of terrain evolution (time used as a third variable and elevation as the grid-cell value), provided interesting insight into the dynamics of different elevation contours during the studied time period (Figure 6). Isosurfaces representing evolution of contours close to the shoreline (z < 0.5 m) had very complex geometry, indicating high dynamics because of shoreline erosion and sand disposal. Isosurfaces with elevations typical for foredunes (z > 5 m) were discontinuous because of variations in elevation over time caused, most likely, by wind transport and overwash. The isosurfaces for elevation values of 1.5 m < z < 4 m (typical for the upper part of the backshore beach) had simple, almost planar, geometry, indicating stable elevation levels.
DISCUSSION
Application of the proposed methodology to the section of a barrier island within the Pea Island National Wildlife Refuge revealed spatial patterns of change that identified the foredunes as the most dynamic feature in the area in terms of variation in elevation. The spatial pattern of elevation gain/loss indicates an important, but often overlooked, role of wind in sand redistribution. The dunes in the northern section of the study area are constantly supplied with windblown sand from the large, flat groin fillet, just to the north, but outside of the study area. Further south, the change in the dune topography is driven by storm overwash. Both mechanisms are important to terrain evolution and should be incorporated into modeling efforts.
In terms of the proposed methodology, assessment of the vertical shifts in lidar data using stable structures, such as roads, was shown to be essential for correct quantification of trends in coastal-terrain change; however, it has rarely been included in coastal-change studies using time series of lidar data. For example, in a section of our study area, the shoreline erosion estimated from the 1999 and 2001 uncorrected DEMs was 12 m, whereas after correction it was only 4 m.
The linear trend for elevation, estimated from uncorrected DEMs between 1997–2003, would have represented sinking of the NC12 highway; however, NCDOT has resurveyed the road benchmarks and confirmed that the road is stable, as indicated by the more recent 2004 and 2005 lidar surveys. Thus, lowering of elevation was indeed due to relatively small systematic error in lidar data that had gradually shifted from the positive to the negative for most (but not all) surveys between 1996 and 2003 (Figure 2). Overall, the actual accuracy of the lidar data on the road was within the published accuracy or better.
CONCLUSIONS
A novel approach to analyzing time series of lidar data with focus on spatial pattern of coastal evolution was proposed and applied to a small section of the North Carolina barrier-island system. The approach is simple to implement, and the resulting maps provide unique information about the location and spatial extent of the most dynamic terrain features and possible stable areas. It also introduces the concept of the core surface, below which, the elevation has never decreased during the study period. The observed spatial patterns of terrain evolution indicated that wind is a major contributor to sand transport in the studied area and identified vulnerable locations with negative elevation trends. The analysis provided detailed geospatial information that can improve the effectiveness of coastal management, such as dune maintenance, targeted beach nourishment, and sand-disposal activities.
Acknowledgments
The support by the U.S. Army Research Office is gratefully acknowledged. We would like to thank Glynn Clements for the development of the new GRASS GIS module r.series and its specific enhancements for this article and Hamish Bowman for the new module r.in.xyz. The lidar data were provided by USGS/NOAA/NASA Coastal Change Program, NCFMP, and Joint Airborne Lidar Bathymetry Center of Expertise.
