An Automated Model to Classify Barrier Island Geomorphology Using Lidar Data †

Limited research has studied the use of Lidar in mapping coastal geomorphology. The purpose of this project was to build on existing research and develop an automated modeling approach to classify the coastal geomorphology of barrier islands and test this at four sites in North Carolina. Barrier islands are shaped by natural coastal processes, such as storms and longshore sediment transport, as well as by human influences, such as beach nourishment and urban development. An automated geomorphic classification model was developed to classify Lidar data into 10 geomorphic types over four time-steps from 1998 to 2014. Tropical storms and hurricanes had the most influence on change and movement. On the developed islands, there was less influence of storms, owing to the inability of features to move because of coastal infrastructure. Beach nourishment was the dominant influence on developed beaches, because this activity ameliorated the natural tendency of an island to erode. Understanding how natural and anthropogenic processes influence barrier island geomorphology is critical to predicting an island’s future response to changing environmental factors such as sea-level rise. The development of an automated model equips policy makers and coastal managers with information to make development and conservation decisions, and the model can be implemented at other barrier islands.


Introduction
Beautiful beaches and expensive properties are found on barrier islands, which are linear features that parallel the coastline and protect the mainland from waves and storms.Their location and sandy composition make barrier islands both economically valuable and physically vulnerable.These systems are dynamic and studies have shown that over five years, a barrier island can migrate over 100 m and experience a 50% change in volume [1,2].Understanding the evolution of barrier island geomorphology can assist policy makers and coastal managers to make decisions regarding future land use development.
The geomorphology of a barrier island depends on several factors such as underlying geology, regional geologic framework, tide regime, and wave energy [3,4].A typical barrier island system is composed of a gently sloping continental shelf, a sandy island, a back-barrier marsh that extends into an estuary, and individual barrier islands are separated by tidal inlets [5].Features on an island include: intertidal and supratidal beach; dunes which are usually the highest elevations on the island and are vegetated with grasses; unvegetated or sparsely vegetated swales and channels between dunes; hummocks (or sometimes referred to as relic dunes) usually vegetated with grasses and scrub/shrub; upland areas which are bare sand or vegetated areas of grasses, scrub/shrub and occasionally maritime forest; overwash fans on the back barrier created when sand is transported via washover events from a large storm; and intertidal and supratidal marsh and scrub/shrub.
In North Carolina, the coastline is fronted by a chain of barrier islands.In particular, the barrier islands in the southern part of the state are distinctively different from those in the north, which is largely due to differences in subsurface geology, coastline orientation, wave height/energy and tidal range [6].In the north, the area is underlain by a thick layer of Quaternary sediments and four large rivers that drain the large coastal plain that is slowly subsiding [6].The mainland is separated from the barrier islands by very large and relatively shallow sounds.This geography and geology, as well as wave energy having more influence than tides, has led to the development of long barrier islands that are divided by only four inlets.In the southern coast, the area is underlain by Cretaceous and Pliocene rocks with relatively little Quaternary sediment, very little subsidence, and a steeper coastal plain in comparison to the northern area.This setting is a mixed energy coast where the islands are influenced by both tides and waves, and has therefore led to the development of barrier islands that are located much closer to the mainland, many inlets and islands that are substantially shorter than those in the north [4].
A variety of processes impact changes in barrier island geomorphology.Sea level and sediment availability influence longshore sediment transport and the width of barrier islands.For example, in regions with abundant sediment supply and slow relative sea-level rise, barrier islands prograde/accrete seaward and can increase in width.Alternatively, low sediment supply and moderate or fast sea−level rise can lead to landward migration and island narrowing [7].In North Carolina, the size, frequency and track of storms has resulted in coastal erosion (shoreline and back barrier) as well as washover that leads to an accumulation of overwash fans [2].
On developed, or urbanized, islands, buildings and infrastructure can trap sand in one location and prevent transport to other locations.Therefore, on developed islands, natural sediment transport is reduced which can result in accretion or erosion depending on the location of infrastructure.Additionally, developed islands prevent natural storm-induced overwash from reaching the back barrier, causing increased erosion and island narrowing [8].Beach nourishment temporarily alters beach topography, and initially protects the island from storm impacts, however it hinders overwash, preventing natural sediment transport, and the effects are obliterated if the area is impacted by several successive storms [2].Nourished segments of beaches have also been shown to experience the largest amount of storm induced erosion [9].
Researchers have mapped the vegetation of barrier islands using a variety of imagery and remote sensing techniques [7,[10][11][12][13][14], but habitat maps do not provide information about the geomorphology which is important with regards to understanding how the islands can change through time.Coastal morphology has been studied using Lidar and other remotely sensed data [1,2,8,[15][16][17][18][19][20]; however, these studies have mapped individual features such as shorelines, dunes, and marshes.To our knowledge, there has been little research addressing a comprehensive mapping method for classifying all geomorphic features on a barrier island, from the ocean to the sound.Therefore, the purpose for this study was to develop and test an automated model for classifying barrier island geomorphology.The first objective was to develop a model using only one data source (Lidar) in order to test whether or not these data can accurately map all barrier island features.If successful, the second objective was to process several time-steps of Lidar data to quantify change over time and correlate results with human and environmental processes that impact barrier islands.
To account for the different geologic settings and anthropogenic influences, this study mapped two islands in the north (one developed and one undeveloped) and two in the south (also one developed and one undeveloped) (Figure 1).These study sites enable a comparison between geologic settings and wave dominated versus mixed energy (wave and tide) environments.It was hypothesized that by testing the model in multiple locations, and is successful, then the model can be implemented in other coastal settings.To account for the different geologic settings and anthropogenic influences, this study mapped two islands in the north (one developed and one undeveloped) and two in the south (also one developed and one undeveloped) (Figure 1).These study sites enable a comparison between geologic settings and wave dominated versus mixed energy (wave and tide) environments.It was hypothesized that by testing the model in multiple locations, and is successful, then the model can be implemented in other coastal settings.

Methods
To address the research objectives, the following steps were undertaken: (1) compile a list of geomorphic features and their characteristics, (2) design and conduct fieldwork, (3) create and test a geospatial model for deriving maps of barrier island features, (4) assess the accuracy of the model, (5) quantify changes through time, and (6) compare changes with coastal processes that may impact the islands.
The following types of geomorphic features were studied: 1. Intertidal: Region that is inundated daily due to tides.Land above Mean Sea level (MSL) and below Mean Higher High Water (MHHW).MSL is defined as the arithmetic mean of hourly water heights of each tidal day over a 19-year period, known as the National Tidal Datum Epoch (NTDE).MHHW is defined as the average of all the higher-high water heights of each tidal day over the NTDE (see NOAA Tides and Currents, https://tidesandcurrents.noaa.gov/datum_options.html).2. Supratidal: Region that is inundated occasionally due to astronomically high tides or severe weather events.Land above MHHW and below Highest Astronomical Tide (HAT) (see NOAA Tides and Currents, https://tidesandcurrents.noaa.gov/datum_options.html).

Methods
To address the research objectives, the following steps were undertaken: (1) compile a list of geomorphic features and their characteristics, (2) design and conduct fieldwork, (3) create and test a geospatial model for deriving maps of barrier island features, (4) assess the accuracy of the model, (5) quantify changes through time, and (6) compare changes with coastal processes that may impact the islands.
The following types of geomorphic features were studied: 1.
Intertidal: Region that is inundated daily due to tides.Land above Mean Sea level (MSL) and below Mean Higher High Water (MHHW).MSL is defined as the arithmetic mean of hourly water heights of each tidal day over a 19-year period, known as the National Tidal Datum Epoch (NTDE).MHHW is defined as the average of all the higher-high water heights of each tidal day over the NTDE (see NOAA Tides and Currents, https://tidesandcurrents.noaa.gov/datum_options.html).

3.
Dune: Linear feature that is parallel to the shoreface and has the highest elevation on the island.

4.
Hummock: Relic dune, usually located behind the primary dune, and is lower elevation than dunes, but higher elevation than other surrounding features (usually upland).They are usually a round shape from erosion due to wind and rain.

5.
Overwash: Slightly elevated and flat areas located in the back barrier and created from sediment transport from the oceanfront.Also known as overwash fans.
Channel: Low depression, cut by water, located adjacent to the supratidal region.8.
Upland: Flat portions of the barrier island, behind the primary dune; all land that is not classified as one of the other features listed above.

Fieldwork
Fieldwork was conducted in May 2016 and October through December, 2016 to collect ground reference points (GRPs) to build and test model classification accuracy.The purpose for the field effort was to randomly sample each study area so that these data would be statistically useful for both building the classification model as well as testing the accuracy of the results.Therefore, each study area was segmented by transects, cast 100 m apart, perpendicular to the island centerline and then 25 transects were randomly selected for the field survey (an example is given in Figure 2).We chose 25 transects because this number resulted in a large sample size that was geographically spaced across each study area and it was estimated that we could map 3-5 transects per day resulting in approximately 30 days of fieldwork (not including travel time).
GRPs were collected using a Trimble R8 Real Time Kinematic (RTK) GPS.Along each transect, the center of each geomorphic feature was recorded with position and elevation (X, Y, Z), feature type, and a GoPro Hero2 was used to collect video.The Trimble R8 has 8 mm horizontal and 15 mm vertical accuracy [21].After fieldwork, the GRPs were post-processed and differentially corrected in Trimble Office, exported as CSV files and imported into ArcGIS.The North Carolina Geodetic Survey maintains a comprehensive network of base stations which enables users to connect to the network and improve spatial positioning of GPS coordinates using differential correction [22].The total number of GRPs at each study area were 657 at Corolla, 583 at Currituck, 540 at Wrightsville, and 300 at Masonboro.

Lidar Data
Lidar data were acquired from the National Oceanic Atmospheric Administration's (NOAA) Digital Coast using the Data Access Viewer tool (www.coast.noaa.gov/dataviewer/#/lidar/search/).The US Army Corps of Engineers, NOAA, and the North Carolina Department of Emergency Management have been mapping the North Carolina coastline for several years and numerous Lidar data are available from 1996 to 2014 with varying spatial extents, accuracies and point densities.Data were downloaded in laser file (LAS) format, with all returns, in North Carolina State Plane North American Datum of 1983 (NAD83), North American Vertical Datum of 1988 (NAVD88), and units in meters.Each dataset was examined for point spacing and accuracy.The highest quality data

Lidar Data
Lidar data were acquired from the National Oceanic Atmospheric Administration's (NOAA) Digital Coast using the Data Access Viewer tool (www.coast.noaa.gov/dataviewer/#/lidar/search/).The US Army Corps of Engineers, NOAA, and the North Carolina Department of Emergency Management have been mapping the North Carolina coastline for several years and numerous Lidar data are available from 1996 to 2014 with varying spatial extents, accuracies and point densities.Data were downloaded in laser file (LAS) format, with all returns, in North Carolina State Plane North American Datum of 1983 (NAD83), North American Vertical Datum of 1988 (NAVD88), and units in meters.Each dataset was examined for point spacing and accuracy.The highest quality data sets were selected for analysis.Different Lidar data were used for the two northern sites and the two in the south because there were no single data sets that covered all four areas.For Currituck and Corolla (northern The LAS files were imported into ArcGIS and LAS datasets were created.Research has tested the spatial resolution for examining volume change in coastal features determined that 1-2 m is optimal [23] and the inverse distance weighted (IDW) interpolation technique was best for producing raster surfaces [9,24].A spatial sensitivity test was conducted using the Lidar ground returns (bare earth) where we tested the type of interpolation (spline and IDW) and multiple neighbors (search radius), maximum distance, and derived cell size.As previous research found, the IDW interpolation technique was the most accurate along with a 10-point search radius, 10 m maximum distance and 1 × 1 m cell size.The derived DEM was tested for accuracy by comparing interpolated elevation values to field collected GRP data.An average elevation difference of less than 10 cm was considered acceptable based on RTK accuracy and the time span from Lidar data collection to fieldwork (2 years).For developed areas (Wrightsville and Corolla), anthropogenic features were extracted from the ground return DEM prior to classifying the geomorphic features.

Geomorphic Classification
The classification model consists of a series of steps that identifies each type of geomorphic feature.The model requires several inputs: (1) a DEM (in this case it was derived from Lidar), (2) a study area boundary polygon, (3) an ocean front line in order to determine marine water from estuarine water, (4) MHHW and HAT tide height measurements (in meters) for the study area and corrected to NAVD88, and (5) several parameters unique to the study area and described in more detail below.The result is a polygon dataset (feature class in an ArcGIS geodatabase) with attributes for each type of geomorphic feature in the study area.
First, the model classifies Intertidal and Supratidal features based on the Mean Sea Level (elevation = 0), the elevation of Mean Higher High Water (MHHW), and the elevation of Highest Astronomical Tide (HAT).Importantly, most Lidar data are not collected at low tide and therefore the full extent of the intertidal zone is not calculated.In this model, the intertidal class extends to mean sea level, not mean lower low water, which would be the full extent of the intertidal habitat.
Second, topographic characteristics can be used to identify features based on the relative change in elevation.The Topographic Position Index (TPI) was developed to identify peaks and valleys in a terrain model [25,26] and in this study, the TPI was used to classify the hummocks and dunes as topographic peaks, swales and channels as topographic valleys, and the overwash and upland areas are relatively flat.
The equation to calculate TPI is Equation (1): For each cell (row i and column j), the focal mean elevation is computed and then compared to the elevation of the cell.A cell that has an elevation that is higher than the focal mean has a positive TPI value, while a cell that has a lower elevation than the focal mean is negative value.Therefore, the focal mean is the average elevation of all cells within a specified neighborhood, or distance, around the cell.Small features have a small neighborhood distance and larger features have larger neighborhoods [26].In this study, the grid cell size was 1 × 1 m and distance sensitivity tests were conducted to identify neighborhood distances for each type of feature where: 12 m = 6 neighbors, 40 m = 20 neighbors, and 200 m = 100 neighbors.The neighborhood distances for each type of feature are given in Table 1.These distances were derived from many tests during model development and therefore this is a parameter that is dependent on the study area.Therefore, the distance to calculate the focal mean is a parameter that the user can change and test for any particular study area.The TPI index values were standardized using the mean and standard deviation of the TPI values for the dataset.This scaling enables comparison across study areas so that the classification model uses the same TPI values for deriving each type of feature.The equation to calculate the scaled TPI is Equation ( 2 Third, proximity and shape metrics were used to refine the feature classification.In this study, the swales and channels have the same focal mean and similar TPI values, so in order to distinguish these features; the proximity to supratidal areas was calculated.Channels are valleys where water cuts through the barrier island, usually perpendicular to the beach, and these features intersect the supratidal region.Swales are valleys between dunes, usually parallel to the beach, so they are adjacent to dunes and do not intersect the supratidal areas. Overwash fans are located close to the back barrier.Therefore, overwash fans were classified by calculating the distance to the oceanfront shoreline and those that were greater than 0.5× the standard deviation of the distance to the shoreline were classified as overwash and other areas closer to the ocean were reclassified as either hummock or dune.Dunes and hummocks have similar TPI values, so a shape index was used to separate these classes.Dunes are generally oval shaped while hummocks are circular.The equation to calculate the shape index for each polygon (p) is Equation (3): Knowledge of the study area is critical to determining the appropriate values for neighborhood distance, TPI, and shape index.Therefore, the model enables users to specify these parameters so they can be appropriate to any study area.Shape index values ≤0.6 were classified as dunes and values >0.6 were classified as hummocks.Not all dunes are long and linear, so the distance between dunes and hummocks was calculated and hummocks that were located closer than 0.5× the standard deviation of the distance to dunes, were reclassified as dunes.The result was a more accurate depiction of the "dune field".Table 1 summarizes the classification parameters for each geomorphic feature type in this study and Figure 3 illustrates the model workflow.and transition from one type of feature to another.Using these matrices, the expected amount of change was computed using the percentage of the study for each type of feature and then compared with the observed amount of change to identify statistically significant transitions from one feature to another [7,27].
Several methods were used to capture feature movement.Oceanfront shoreline dynamics were compared using the line shared by supratidal and intertidal and then analyzed using AMBUR [28].AMBUR is similar to the Digital Shoreline Analysis System (DSAS) [29] where the difference in two shorelines are compared at equally spaced transects and the rate of change (erosion or accretion) is calculated [28].Movement/migration of dune features was calculated by using the "Detect Feature Change" and "Near" tools in ArcMap [30] in order to measure the distance between features through time.The difference in spatial position (distance) was classified according to the amount of change: (1) no change (<3 m), ( 2) movement (3 m ≤ distance ≤ 25 m), ( 3) deletion (feature completely eroded), and (4) new dune (>25 m).

Results
The geomorphic classification model was developed in Esri's ModelBuilder [30].The model was run 16 times (four study areas and four dates), resulting in classified maps for each study area (Figure 4).  to create cross-tabulation, or classification change, matrices which summarize the amount of gain, loss, and transition from one type of feature to another.Using these matrices, the expected amount of change was computed using the percentage of the study for each type of feature and then compared with the observed amount of change to identify statistically significant transitions from one feature to another [7,27].
Several methods were used to capture feature movement.Oceanfront shoreline dynamics were compared using the line shared by supratidal and intertidal and then analyzed using AMBUR [28].AMBUR is similar to the Digital Shoreline Analysis System (DSAS) [29] where the difference in two shorelines are compared at equally spaced transects and the rate of change (erosion or accretion) is calculated [28].Movement/migration of dune features was calculated by using the "Detect Feature Change" and "Near" tools in ArcMap [30] in order to measure the distance between features through time.The difference in spatial position (distance) was classified according to the amount of change:

Results
The geomorphic classification model was developed in Esri's ModelBuilder [30].The model was run 16 times (four study areas and four dates), resulting in classified maps for each study area (Figure 4).

Accuracy Assessment
First, we tested the accuracy of the elevations from the most recent Lidar data (2014) by comparing them with the field GRPs collected in 2016 and there was very good correspondence.The mean difference was 0.02 m at Currituck, 0.06 m at Masonboro, 0.07 m at Wrightsville, and 0.12 m at Corolla.These results were substantially smaller than the 20 cm vertical accuracy of the Trimble RTK.A small percentage of GRPs had elevation differences greater than 1 m: 8/300 (2.67%) on Masonboro Island, 12/540 (2.22%) on Wrightsville Beach, 20/583 (3.43%) on Currituck, and 24/657 (3.65%) on Corolla.The greatest/largest positive elevation difference (fieldwork GRP minus Lidar) was on dunes and largest negative elevation was on intertidal features.Based on these results, the fieldwork GRPs were comparable to the Lidar data.
Next, the model classification accuracy was computed using the feature types recorded at the GRPs in comparison to the model output.There was an average classification accuracy of 83.14%; however, class accuracies varied where some classes were classified more accurately.For example, the most accurate features were buildings, roads, intertidal, dunes, and swales versus, hummocks, overwash, channels, and upland were less accurate (Table 2).These results were expected since defining the boundary of these features can be difficult.The overall percentage of each type of feature within the study areas varies greatly and therefore the number of GRPs also varied.In particular, other than Masonboro, there were very few hummocks along the transects at the other study areas so there were far fewer GRPs for this feature class (28) and more than half were classified as a dune.Overwash was similar with 62 GRPs and 17 were misclassified as upland.Two classes with plenty of GRPs had less than ideal classification accuracies: supratidal 28% omission error and channel had 24% omission error.For supratidal, most GRP error was channel and for the channel class, the most error was misclassifying these as swale and upland.Given the overall good model accuracy (83%), which exceeded our goal of 75% accuracy, the earlier Lidar data sets were then processed through the model and classified.

Accuracy Assessment
First, we tested the accuracy of the elevations from the most recent Lidar data (2014) by comparing them with the field GRPs collected in 2016 and there was very good correspondence.The mean difference was 0.02 m at Currituck, 0.06 m at Masonboro, 0.07 m at Wrightsville, and 0.12 m at Corolla.These results were substantially smaller than the 20 cm vertical accuracy of the Trimble RTK.A small percentage of GRPs had elevation differences greater than 1 m: 8/300 (2.67%) on Masonboro Island, 12/540 (2.22%) on Wrightsville Beach, 20/583 (3.43%) on Currituck, and 24/657 (3.65%) on Corolla.The greatest/largest positive elevation difference (fieldwork GRP minus Lidar) was on dunes and largest negative elevation was on intertidal features.Based on these results, the fieldwork GRPs were comparable to the Lidar data.
Next, the model classification accuracy was computed using the feature types recorded at the GRPs in comparison to the model output.There was an average classification accuracy of 83.14%; however, class accuracies varied where some classes were classified more accurately.For example, the most accurate features were buildings, roads, intertidal, dunes, and swales versus, hummocks, overwash, channels, and upland were less accurate (Table 2).These results were expected since defining the boundary of these features can be difficult.The overall percentage of each type of feature within the study areas varies greatly and therefore the number of GRPs also varied.In particular, other than Masonboro, there were very few hummocks along the transects at the other study areas so there were far fewer GRPs for this feature class (28) and more than half were classified as a dune.Overwash was similar with 62 GRPs and 17 were misclassified as upland.Two classes with plenty of GRPs had less than ideal classification accuracies: supratidal 28% omission error and channel had 24% omission error.For supratidal, most GRP error was channel and for the channel class, the most error was misclassifying these as swale and upland.Given the overall good model accuracy (83%), which exceeded our goal of 75% accuracy, the earlier Lidar data sets were then processed through the model and classified.

Geomorphology Change through Time
Changes through time were measured by analyzing: (1) area and elevation, (2) statistically significant changes using cross-correlation matrices and (3) shoreline and dune movement.

Change in Area and Elevation
Upland was the largest feature on Wrightsville and Currituck at 30% of the area and on Corolla, upland was 40%.Alternatively, on Masonboro, intertidal and supratidal were the largest area at ~50% of the island.As expected, at all four study sites, intertidal and supratidal features had the lowest average elevation and dunes had the highest average elevation.Dunes, overwash, and upland areas experienced the most change in the percentage of each study area and change in elevation (Figure 5).Linear regression analysis revealed that the change in area (percentage) had negative slopes (decreasing area) for dunes and upland, and positive slopes for overwash (increasing area) (Table 3).
The elevation of dunes, overwash and upland was higher in the northern areas (Corolla and Currituck) than the southern areas (Masonboro and Wrightsville).The northern areas also experienced more change in elevation through time whereas the southern areas did not change as much.For example, at Masonboro, elevation changed by only +/− 0.1 m for all features except dunes, and they only changed by an average of 0.2 m.The most change in elevation was during the 2005-2009/2010 period, when most upland features (dunes, hummocks, overwash, channels, swales, and upland) had a consistent drop in elevation, especially in the northern areas.The average decrease in elevation for upland features was 1.54 m at Corolla and 1.1 m at Currituck.

Geomorphology Change through Time
Changes through time were measured by analyzing: (1) area and elevation, (2) statistically significant changes using cross-correlation matrices and (3) shoreline and dune movement.

Change in Area and Elevation
Upland was the largest feature on Wrightsville and Currituck at 30% of the area and on Corolla, upland was 40%.Alternatively, on Masonboro, intertidal and supratidal were the largest area at ~50% of the island.As expected, at all four study sites, intertidal and supratidal features had the lowest average elevation and dunes had the highest average elevation.Dunes, overwash, and upland areas experienced the most change in the percentage of each study area and change in elevation (Figure 5).Linear regression analysis revealed that the change in area (percentage) had negative slopes (decreasing area) for dunes and upland, and positive slopes for overwash (increasing area) (Table 3).
The elevation of dunes, overwash and upland was higher in the northern areas (Corolla and Currituck) than the southern areas (Masonboro and Wrightsville).The northern areas also experienced more change in elevation through time whereas the southern areas did not change as much.For example, at Masonboro, elevation changed by only +/− 0.1 m for all features except dunes, and they only changed by an average of 0.2 m.The most change in elevation was during the 2005-2009/2010 period, when most upland features (dunes, hummocks, overwash, channels, swales, and upland) had a consistent drop in elevation, especially in the northern areas.The average decrease in elevation for upland features was 1.54 m at Corolla and 1.1 m at Currituck.Comparing the amount of gain and loss at the four study areas provides more insight into how the areas changed than a simple percent change in feature.For example, Masonboro had substantially more gains and losses than the other three areas especially during from 2005 to 2010 when there was a 28% gain in intertidal area and 40% loss in supratidal areas (Figure 6).In the same time period (2005-2010) at Wrightsville, upland had the largest amount of gain (12%) and loss (−14%).This is an example of the creation of new upland area (gain) while simultaneously losing this class in other places in the study area.Therefore, the gains and losses illustrate how much the features are changing within the study areas.From 1998 to 2005, at Masonboro, there was a large gain in supratidal (25%) while the intertidal had substantial loss (−20%).These results compliment the change in elevation where Masonboro had the least change in elevation while it had the most change in gains and losses.On this undeveloped island, there is relative stability, or equilibrium, in the elevation, but substantial change in the oceanfront and back barrier intertidal and supratidal habitats.The other study areas had substantially more gains and losses in the other features in comparison with Masonboro.For example, Corolla and Currituck had the largest amount of gains and losses in dunes and upland from 2005 to 2009.Wrightsville had the least amount of gains and losses, which corroborates the lack of change on highly developed islands.

Geomorphology Change with Respect to Gain and Loss
Comparing the amount of gain and loss at the four study areas provides more insight into how the areas changed than a simple percent change in feature.For example, Masonboro had substantially more gains and losses than the other three areas especially during from 2005 to 2010 when there was a 28% gain in intertidal area and 40% loss in supratidal areas (Figure 6).In the same time period (2005-2010) at Wrightsville, upland had the largest amount of gain (12%) and loss (−14%).This is an example of the creation of new upland area (gain) while simultaneously losing this class in other places in the study area.Therefore, the gains and losses illustrate how much the features are changing within the study areas.From 1998 to 2005, at Masonboro, there was a large gain in supratidal (25%) while the intertidal had substantial loss (−20%).These results compliment the change in elevation where Masonboro had the least change in elevation while it had the most change in gains and losses.On this undeveloped island, there is relative stability, or equilibrium, in the elevation, but substantial change in the oceanfront and back barrier intertidal and supratidal habitats.The other study areas had substantially more gains and losses in the other features in comparison with Masonboro.For example, Corolla and Currituck had the largest amount of gains and losses in dunes and upland from 2005 to 2009.Wrightsville had the least amount of gains and losses, which corroborates the lack of change on highly developed islands.

More Change Than Expected
Polygon overlay and cross-tabulation change matrices were generated for each study area and time period.Although it is important to identify the largest changes (in percent) where geomorphic features transition from one type to another, it is also important that we identify transitions that have occurred more than the expected amount.We do this by first calculating the expected amount of change using the percentage of each feature in time 1, and then calculate the deviation as observed minus expected [27].Lastly, the rate, or a measure of significance, is calculated as deviation divided by expected.Using this technique enables the transitions/changes to be identified that are more than expected regardless of the size, or proportion of the study area, that is dominated by the feature [27].In all four study areas and time periods, the intertidal and supratidal swapped more than expected.This indicates how consistently variable these features can be.

More Change Than Expected
Polygon overlay and cross-tabulation change matrices were generated for each study area and time period.Although it is important to identify the largest changes (in percent) where geomorphic features transition from one type to another, it is also important that we identify transitions that have occurred more than the expected amount.We do this by first calculating the expected amount of change using the percentage of each feature in time 1, and then calculate the deviation as observed minus expected [27].Lastly, the rate, or a measure of significance, is calculated as deviation divided by expected.Using this technique enables the transitions/changes to be identified that are more than expected regardless of the size, or proportion of the study area, that is dominated by the feature [27].In all four study areas and time periods, the intertidal and supratidal swapped more than expected.This indicates how consistently variable these features can be.At Corolla and Currituck, almost all transitions were between intertidal, supratidal and water for the 2001-2005 and 2005-2010 time periods, but then from 2009-2014 there was much more change in the upland features (e.g., hummocks to dunes, upland and overwash).At Masonboro, upland features dominated the transitions: upland to swale, hummock and dune; dune to swale; swale to channel; and hummocks to dune and channel (Figure 7).At all four study areas, the number of significant transitions increased substantially in the last time period.The change matrices for all study areas and time periods, including expected and deviations, are given in the supplementary materials.
The average rate of deviation to expected change was calculated for each study area and time period (Figure 8).Interestingly, the deviation/expected values in the northern region (Currituck and Corolla) were much higher than the two southern areas (Masonboro and Wrightsville) indicating these areas experienced much more change.These results compliment the change in elevation results, which were greatest in the northern areas.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 21 hummocks to dunes, upland and overwash).At Masonboro, upland features dominated the transitions: upland to swale, hummock and dune; dune to swale; swale to channel; and hummocks to dune and channel (Figure 7).At all four study areas, the number of significant transitions increased substantially in the last time period.The change matrices for all study areas and time periods, including expected and deviations, are given in the supplementary materials.
The average rate of deviation to expected change was calculated for each study area and time period (Figure 8).Interestingly, the deviation/expected values in the northern region (Currituck and Corolla) were much higher than the two southern areas (Masonboro and Wrightsville) indicating these areas experienced much more change.These results compliment the change in elevation results, which were greatest in the northern areas.

Shoreline Change and Dune Movement
Oceanfront shoreline positions were compared and the northern region was accreting from 2001 to 2005 when Corolla had 78% of the shoreline accreting and Currituck had 66%.From 2005 to 2009, even more of the shoreline was accreting at Corolla and Currituck.Conversely, from 2009 to 2014, almost all (more than 90%) of the shoreline experienced a large amount of erosion.In the southern areas, from 1998 to 2005, 72% of the shoreline eroding at Masonboro and 52% of the shoreline was eroding at Wrightsville.Shoreline accretion rates increased and erosion rates decreased from 2005 to 2010.From 2010 to 2014, the shoreline was again dominated by erosion where Masonboro had 72% of the shoreline eroding and Wrightsville had 83% of the shoreline eroding (Figure 9).

Shoreline Change and Dune Movement
Oceanfront shoreline positions were compared and the northern region was accreting from 2001 to 2005 when Corolla had 78% of the shoreline accreting and Currituck had 66%.From 2005 to 2009, even more of the shoreline was accreting at Corolla and Currituck.Conversely, from 2009 to 2014, almost all (more than 90%) of the shoreline experienced a large amount of erosion.In the southern areas, from 1998 to 2005, 72% of the shoreline eroding at Masonboro and 52% of the shoreline was eroding at Wrightsville.Shoreline accretion rates increased and erosion rates decreased from 2005 to 2010.From 2010 to 2014, the shoreline was again dominated by erosion where Masonboro had 72% of the shoreline eroding and Wrightsville had 83% of the shoreline eroding (Figure 9).Movement/migration of dune features was calculated by measuring the difference in spatial position through time.At all study areas, very few dunes had no change and approximately 25% of the dunes were new or deleted.Therefore, the majority (approximately 75%) of dune change was movement between 3 and 25 m.Unlike the previous results which highlighted the differences between the northern and southern study areas, the differences in dune movement were much greater between developed (Corolla and Wrightsville) and undeveloped (Currituck and Masonboro) (Table 4).The developed areas had much higher dune movement, especially during the first two time periods when the shorelines were accreting, and the most recent period had the least amount of movement at all study areas, when the shorelines were eroding.

Discussion
Previous studies of barrier island geomorphology have delineated shorelines and dunes by digitizing aerial photographs [8,31,32] or subtracting DEMs [2,9,17,19,23,33].Others have developed automated processes to map shoreline positions from DEMs [34] and used Lidar to compute beach profiles [35], but these studies did not map all barrier island features.In this project, we have Movement/migration of dune features was calculated by measuring the difference in spatial position through time.At all study areas, very few dunes had no change and approximately 25% of the dunes were new or deleted.Therefore, the majority (approximately 75%) of dune change was movement between 3 and 25 m.Unlike the previous results which highlighted the differences between the northern and southern study areas, the differences in dune movement were much greater between developed (Corolla and Wrightsville) and undeveloped (Currituck and Masonboro) (Table 4).The developed areas had much higher dune movement, especially during the first two time periods when the shorelines were accreting, and the most recent period had the least amount of movement at all study areas, when the shorelines were eroding.

Discussion
Previous studies of barrier island geomorphology have delineated shorelines and dunes by digitizing aerial photographs [8,31,32] or subtracting DEMs [2,9,17,19,23,33].Others have developed automated processes to map shoreline positions from DEMs [34] and used Lidar to compute beach profiles [35], but these studies did not map all barrier island features.In this project, we have developed and tested an automated Lidar classification model for mapping all barrier island features from the ocean to the sound.The overall map accuracy for the model was over 83%, but that also means that the model does not do a good job with 17% of the GRPs.The model uses spatial metrics, such as proximity and shape, but future work could further develop these and other metrics to more accurately map the geomorphology.There was an assumption that randomly selected field transects would provide ample ground referencing data, but this resulted in some types of features (e.g., hummock, overwash, and channel) with substantially fewer reference points.Now that the study areas have been mapped, future research could employ a different sampling technique, such as spatially balanced points, which would provide more samples regardless of the size or prevalence of each feature [36].
Lidar is not the only type of elevation data that can be used in the model.Future work could test the classification of barrier island geomorphology using unmanned aerial vehicle (UAV) imagery; however, a digital terrain model (DTM) is not equivalent to a DEM and therefore fieldwork should be conducted to correctly calibrate UAV imagery [37][38][39].The areas mapped in this study have plenty of Lidar data; however, this is not always the case and so UAV imagery may be a useful alternative.Conversely, the cost of a drone, sensors, and software that is capable of producing mapping-grade products is rather expensive so even though technology is quickly advancing, there will always be a concern for cost.
Results from the shoreline change analysis documented a prevalence of oceanfront erosion, but there are also places where the shoreline is accreting.The dune movement analysis confirmed this with a large amount of lateral movement in the southwest direction.Change matrices showed significant change (more than expected) in intertidal, supratidal, dunes, overwash, and channels, but less significant changes in hummocks, upland, swales, and water (estuarine and ocean).These measurements provided evidence for where islands are narrowing, widening, and migrating landward.
Regional differences between the north (Currituck and Corolla) and the south (Masonboro and Wrightsville) were much larger than differences between developed and undeveloped barrier islands.There are two primary reasons why the north is different from the south: geologic setting [6,[40][41][42][43][44] and beach nourishment [45].Developed islands had less overall change in geomorphology because development prevents natural processes such as washover and rollover (island migration landward).This results in the developed islands geomorphic features being "locked" in place which creates increased shoreline erosion and island narrowing [6].Beach nourishment temporarily increases the amount of sediment and overall elevation of the oceanfront shoreline [2].However, nourishment has been shown to result in the largest amount of storm induced erosion [9].Beach nourishment only occurred in the southern areas (Masonboro and Wrightsville in 2006 and 2010), just prior to the Lidar data collection.From 2005 to 2010 the average shoreline change rate was +1.2 m/year at Masonboro and +0.5 m/year at Wrightsville.This accretion was immediately followed by a period of high erosion from 2010 to 2014 at −1.7 m/year at Masonboro and −3.5 m/year at Wrightsville.The location of accretion corresponded to the areas of highest erosion.
Storms were a dominant process influencing the four study areas (Figure 10).Post-storm analysis computed the average shoreline change of −3 m/year (Masonboro), −0.45 m (Wrightsville), −3.4 m/year (Currituck), and −2.7 m/year (Corolla).Feature change analysis documented dune erosion and the transition of dunes to supratidal, intertidal and channels.On the undeveloped islands, there was an increase in overwash.In the southern study areas (Masonboro and Wrightsville), the stormiest period was 1998 to 2005 when the area was impacted by eight major storms, four of which were hurricanes.Only a few small storms impacted the northern study areas between 2001 and 2009, and then in 2011 Hurricane Irene, a category 3, passed directly through the region, likely responsible for the changes observed on Currituck and Corolla.The model was developed and tested at four islands in North Carolina (two developed and two undeveloped); however, one of the primary goals for this project was to create a model that could be applied to other barrier islands.As such, the model parameters can be adapted to reflect the characteristics of other study areas.In order to reflect local characteristics, the model requires tide information and expert knowledge for the determination of: (1) size of geomorphic features to determine the neighborhood distance for calculating the focal mean used in TPI calculation, (2) sensitivity testing to most accurately determine the breaks in TPI index values for each type of feature, (3) distance/proximity to supratidal features to separate channels from swales, (4) distance to oceanfront shoreline to separate overwash from dunes and hummocks, (5) shape index values to separate dunes from hummocks and (6) distance between dunes and hummocks in order to reclassify hummocks as dunes of they are very close to each other.Although there are many parameters to consider, this is the first comprehensive barrier island geomorphology model and future research can build on this approach by testing more fully automated methods such as data mining and pattern recognition.

Conclusions
One of the benefits of using Lidar data for studying coastal systems is that it provides elevation; therefore, in sandy coastal environments, where many different features have similar reflectance properties, features can be distinguished on the basis of topography [19].This study developed an automated model to classify barrier island geomorphic features.On average, the model accuracy was 83%, which is acceptable given that it can be difficult to precisely measure the boundaries of different types of features that gently vary over the terrain.Historical Lidar can be used to analyze change through time, and geospatial techniques can measure the distance and direction that features have moved.The model was tested at undeveloped and developed islands and islands with different geologic settings.When storms occur, they are the dominant force influencing change.The model was developed and tested at four islands in North Carolina (two developed and two undeveloped); however, one of the primary goals for this project was to create a model that could be applied to other barrier islands.As such, the model parameters can be adapted to reflect the characteristics of other study areas.In order to reflect local characteristics, the model requires tide information and expert knowledge for the determination of: (1) size of geomorphic features to determine the neighborhood distance for calculating the focal mean used in TPI calculation, (2) sensitivity testing to most accurately determine the breaks in TPI index values for each type of feature, (3) distance/proximity to supratidal features to separate channels from swales, (4) distance to oceanfront shoreline to separate overwash from dunes and hummocks, (5) shape index values to separate dunes from hummocks and (6) distance between dunes and hummocks in order to reclassify hummocks as dunes of they are very close to each other.Although there are many parameters to consider, this is the first comprehensive barrier island geomorphology model and future research can build on this approach by testing more fully automated methods such as data mining and pattern recognition.

Conclusions
One of the benefits of using Lidar data for studying coastal systems is that it provides elevation; therefore, in sandy coastal environments, where many different features have similar reflectance properties, features can be distinguished on the basis of topography [19].This study developed an automated model to classify barrier island geomorphic features.On average, the model accuracy was 83%, which is acceptable given that it can be difficult to precisely measure the boundaries of different types of features that gently vary over the terrain.Historical Lidar can be used to analyze change through time, and geospatial techniques can measure the distance and direction that features have moved.The model was tested at undeveloped and developed islands and islands with different geologic settings.When storms occur, they are the dominant force influencing change.Between stormy periods, other activities, such as beach nourishment, can temporarily increase the oceanfront elevation, which later leads to the greatest rates of erosion.Lastly, urban development reduces the amount of change across the entire island because natural processes (e.g., washover) are prohibited from moving sediment to adjacent back-barrier areas.However, shorelines and dune movement changed the most on developed islands.Geologically, the study areas in the north had the greatest decrease in elevation in comparison to the southern areas which were generally more stable through time.

Figure 1 .
Figure 1.Study areas: Wrightsville Beach and Masonboro Island are located in New Hanover County in southern North Carolina (NC) and Currituck and Corolla are in Currituck County in the north.

Figure 1 .
Figure 1.Study areas: Wrightsville Beach and Masonboro Island are located in New Hanover County in southern North Carolina (NC) and Currituck and Corolla are in Currituck County in the north.

Figure 2 .
Figure 2. Example fieldwork transects (spaced 100 m apart) at Wrightsville.Twenty-five randomly selected transects, shown in cyan, were surveyed and ground reference points (GRP) were collected for each type of geomorphic feature.

Figure 2 .
Figure 2. Example fieldwork transects (spaced 100 m apart) at Wrightsville.Twenty-five randomly selected transects, shown in cyan, were surveyed and ground reference points (GRP) were collected for each type of geomorphic feature.

Figure 3 .
Figure 3. Model flowchart to derive barrier island geomorphology.Input and output are shown in blue, processing steps are yellow, intermediate outputs are light green and final classification outputs are dark green.

2. 4 .
Temporal Change Analysis Classified maps of each study area were analyzed for change over time.In the northern areas, the time steps were 2001-2005, 2005-2009, and 2009-2014, and in the southern areas the time steps were 1998-2005, 2005-2010, and 2010-2014.Polygon overlay and change statistics were calculated to create cross-tabulation, or classification change, matrices which summarize the amount of gain, loss,

Figure 3 .
Figure 3. Model flowchart to derive barrier island geomorphology.Input and output are shown in blue, processing steps are yellow, intermediate outputs are light green and final classification outputs are dark green.

Figure 5 .
Figure 5.The average elevation (in m) and percentage of each study area for dunes (A); overwash (B) and upland (C).The area (measured in percent) is shown in the bars with the left axis and the elevation (measured in m) is shown in lines with the right axis.

Figure 5 .
Figure 5.The average elevation (in m) and percentage of each study area for dunes (A); overwash (B) and upland (C).The area (measured in percent) is shown in the bars with the left axis and the elevation (measured in m) is shown in lines with the right axis.
At Corolla and Currituck, almost all transitions were between intertidal, supratidal and water for the 2001-2005 and 2005-2010 time periods, but then from 2009-2014 there was much more change in the upland features (e.g.,

Figure 8 .
Figure 8.Average ratio of deviation/expected at the four study areas and three time periods.

Figure 8 .
Figure 8.Average ratio of deviation/expected at the four study areas and three time periods.

21 Figure 8 .
Figure 8.Average ratio of deviation/expected at the four study areas and three time periods.3.2.4.Shoreline Change and Dune MovementOceanfront shoreline positions were compared and the northern region was accreting from 2001 to 2005 when Corolla had 78% of the shoreline accreting and Currituck had 66%.From 2005 to 2009, even more of the shoreline was accreting at Corolla and Currituck.Conversely, from 2009 to 2014, almost all (more than 90%) of the shoreline experienced a large amount of erosion.In the southern areas, from 1998 to 2005, 72% of the shoreline eroding at Masonboro and 52% of the shoreline was eroding at Wrightsville.Shoreline accretion rates increased and erosion rates decreased from 2005 to 2010.From 2010 to 2014, the shoreline was again dominated by erosion where Masonboro had 72% of the shoreline eroding and Wrightsville had 83% of the shoreline eroding (Figure9).

Figure 10 .
Figure 10.Named tropical storms from 1998 through 2014.Also shown are the dates for Lidar data and when beach nourishment took place at Masonboro and Wrightsville.

Figure 10 .
Figure 10.Named tropical storms from 1998 through 2014.Also shown are the dates for Lidar data and when beach nourishment took place at Masonboro and Wrightsville.

Table 2 .
Confusion matrix for each type of geomorphic feature.CE is confusion error and OE is omission error.

Table 3 .
Regression analysis for area (percent) dune, overwash and upland at all four study areas.Values with higher R 2 are highlighted in gray.

Table 3 .
Regression analysis for area (percent) dune, overwash and upland at all four study areas.Values with higher R 2 are highlighted in gray.

Table 4 .
Average shoreline change (m/year) and dune movement (m) by study area and time period.

Table 4 .
Average shoreline change (m/year) and dune movement (m) by study area and time period.