Multi-Temporal Trend Analysis of Coastal Vegetation Using Metrics Derived from Hyperspectral and LiDAR Data

Abstract

Monitoring and modeling of coastal vegetation and wetland systems are considered major challenges, especially when considering environmental response to hazards, disturbances, and management activities. Remote sensing applications can provide alternatives and complementary approaches to the often costly and laborious field-based collection methods traditionally used for coastal ecosystem monitoring. New and improved sensors and data analysis techniques have become available, making remote sensing applications attractive for evaluation and potential use in monitoring coastal vegetation properties and ecosystem conditions and change. This study involves the extraction of vegetation metrics from airborne LiDAR (Light Detection and Ranging) and hyperspectral imagery (HSI) to quantify coastal dune vegetation characteristics and assesses landscape-level trends from those derived metrics. HSI- and LiDAR-derived elevation (digital elevation model) and vegetation metrics (canopy height model, leaf area index, and normalized difference vegetation index) were used in conjunction with per-pixel linear regression and hot spot analyses to evaluate hurricane-induced spatial and temporal changes in elevation and vegetation properties. These assessments showed areas with greatest decreases in vegetation metric values were associated with direct tropical storm energies and processes (i.e., overwashing events eroding beach and dune features), while those with the greatest increases in vegetation metric values were in areas where overwashed sediments were distributed. This study narrows existing gaps in dune vegetation data by advancing new methodologies to classify, quantify, and estimate critical coastal vegetation metrics. The tools and methods developed in this study will ultimately improve future estimates and predictions of nearshore dynamics and impacts from disturbance events.

1. Introduction

Coastal landscape features (e.g., beaches, dunes, wetlands) provide benefits that range from regulating services (floods, drought, and land degradation), supporting services (soil formation and nutrient cycling), and provisioning services (food and fresh water), to maintaining high biological productivity and serving as critical habitat for fish and wildlife. Coastal dunes also yield substantial economic benefits by reducing damage to residential and commercial infrastructure during significant storm events. However, these dynamic environments are inherently vulnerable and geomorphologically unstable by nature of their position at the interface of land and sea [1]. While coastal beach, dune, and wetland features frequently undergo nominal changes due to natural fair-weather pressures (i.e., wind, waves, and tides), the effects of climate change, especially large storms and irregular weather patterns, can induce significant shifts in the state, evolution, and dynamics of these features, thereby increasing the susceptibility of these systems to destructive erosional forces [2,3]. Therefore, many restoration activities and coastal reconstruction projects focus on dune resiliency to help reduce the risk to coastal communities [4].

Vegetation is a critical biotic component of coastal ecosystems since they have direct and indirect impacts on system stability and resilience [1]. Plants can contribute to erosion control through their above- and belowground components. For example, high leaf-area and stem densities in aboveground plant mass can reduce rain and wave splash impacts, mitigate storm surge, and diminish near-bed shear stresses, as well as facilitate vertical and horizontal accretion through organic matter accumulation and sediment trapping [5]. Belowground root mass can also increase the sheer strength of the substrate and thereby prevent erosion by resisting structural failure and creating shore bathymetric profiles that assist in dissipating wave energy [6]. Therefore, within the coastal ecosystem, dune vegetation productivity and processes play a complex and vital role in shoreline stability and wetland resilience [6]. Estimating vegetation characteristics (i.e., biomass) can be essential for monitoring and predicting wetland productivity, resilience, and ecosystem health [7,8]. Thus, rapidly quantifying vegetation dynamics across the dune profile (from the dune toe to beyond the more elevated and heavily vegetated primary dune) is an important aspect for understanding dune morphology and associated impacts from coastal processes [9].

Recent studies have focused on evaluating physical and biological processes and the long-term evolution of coastal systems as functions of natural processes and anthropogenic activities [7,10,11]. Although some emphasis has been on coastal vegetation monitoring and influence on erosion potential, fewer studies have directly evaluated correlations to vegetation structural properties (i.e., vegetative densities), presumably in part due to limited and insufficient data to support these efforts. Similarly, there are data and knowledge gaps regarding required vegetation inputs to coastal numerical models, particularly for large-scale or regional coastal projects, resulting in many models incorporating synthetic vegetation data or constant values as inputs [12]. Despite ongoing improvements, monitoring and modeling are still considered major challenges in anticipating environmental response to hazards, disturbances, and management activities.

Some major challenges associated with monitoring of coastal wetland landscapes are the cost and laborious nature of traditional ground-based data collections. Remote sensing systems and applications have been shown to provide cost-effective and more efficient alternatives and complementary approaches to ground-based collections and assessments [13]. However, monitoring wetland landscapes and measuring vegetation biophysical and biochemical properties using traditional broadband sensors (usually up to six spectral bands) [14] can be limited (due to spectral resolution) and ineffective (due to the underlying soil types, canopy and leaf properties, and atmospheric conditions), leading to asymptotic saturation for some vegetation metrics [8,15]. Other challenges associated with the monitoring of coastal landscapes are related to the quantification of changes in vegetation structure or function at scales (i.e., spatial, spectral, and temporal) necessary to adequately evaluate system stability. However, studies have shown that narrowband data (i.e., hyperspectral imagery [HSI]) are more sensitive to differences in energy wavelengths and have higher potential to detect subtle variations in vegetation characteristics, as compared to broadband multispectral imagery [15,16,17,18]. Higher spectral information is useful for estimating important vegetation photosynthetic activities and canopy structural variations [19]. Similarly, Light Detection and Ranging (LiDAR) data can provide efficient and effective alternatives to ground-based structural measurements. Some LiDAR sensors can penetrate vegetation cover to obtain estimates of relative heights above ground and detailed vegetation structure parameters (depending on point spacing) [20,21,22]. Using high-resolution HSI imagery in conjunction with LiDAR data can assist in discriminating plant species, determining biomass distribution, and extracting key coastal dune vegetation metrics [23,24,25,26], all of which are useful in monitoring and modeling changes in expansive wetland ecosystems [27].

The goal of this study was to develop simple but effective methods for mapping and monitoring coastal vegetation properties within the dune complex by capitalizing on recent sensor improvements, analytical techniques, and the temporal resolution of data collected by the U.S. Army Corps of Engineers’ (USACE) National Coastal Mapping Program (NCMP). The USACE NCMP provides regional, repetitive, high-resolution, high-accuracy elevation and imagery data that are useful for assessing changes and trends in geomorphological and environmental properties. Therefore, the objectives of this study were to (1) develop tools to semi-automate metric extraction from NCMP LiDAR and HSI data to quantify dune vegetation biological characteristics; and (2) use a pilot study to demonstrate how derived metrics can be utilized to evaluate and compare spatial and temporal changes in vegetation properties, whether due to episodic events (i.e., hurricane impacts) or natural processes. This study advances existing dune and wetland system knowledge and remote sensing techniques by utilizing spatially and temporally coincident HSI and LiDAR data for developing new tools to more efficiently classify, quantify, and estimate critical dune vegetation metrics, thereby improving future estimates and predictions of storm processes and impacts.

2. Materials and Methods

2.1. Study Area

The study site (Figure 1), which is located on Pea Island in the Outer Banks of North Carolina, United States (approximately 170 km east of Greenville), was chosen to evaluate methods for extracting dune vegetation (DuneVeg) metrics, develop an automated workflow for remote extraction techniques, and perform a trend analysis using the derived vegetation metrics. The site is part of the US Coastal Research Program’s DUring Nearshore Event eXperiment (DUNEX), which is a multi-agency, academia, and stakeholder collaboration to study nearshore processes during coastal storm events [19,21,22]. The Pea Island site is part of the Pea Island National Wildlife Refuge and is approximately 19 km in length and 2.4 km wide.

2.2. Vegetation Metric Extraction

2.2.1. Airborne Hyperspectral Imagery and LiDAR Data

This study utilized coastal HSI and LiDAR data routinely collected for the USACE NCMP to evaluate extraction methods of coastal vegetation characteristics. These data are collected by the Joint Airborne LiDAR Bathymetry Technical Center of Expertise [24] to capture existing beach and near-shore conditions along the sandy coastlines of the U.S. The NCMP uses the Coastal Zone Mapping and Imaging Lidar (CZMIL) system, which integrates a Teledyne Optech LiDAR sensor with topographic (70 kHz pulse rate) and bathymetric (10 kHz pulse rate) capabilities, a digital camera, and the Itres Compact Airborne Spectrographic Imager (CASI)-1500 on a single remote sensing platform mounted on a fixed-wing aircraft. The CZMIL HSI and LiDAR data are acquired spatially and temporally coincident and at a consistent scale and resolution in the coastal zone [27]. Thus, they represent a unique opportunity to use open-access, freely available imagery (in this case HSI) with high-fidelity LiDAR data to extract regional-scale vegetation metrics at project-relevant spatial resolutions, which is not currently possible with standard cameras or other available imagery sources.

The HSI data were collected at an altitude of 400 m across a 300 m ground swath. Although the CASI-1500 imager is a programmable sensor capable of hundreds of narrow spectral bands (288 channels at a spectral resolution of 3.5 nanometer [nm] Full Width Half Max [FWHM]), the data were collected in a 48-band configuration with a spectral range of 362 to 1045 nm at a spectral resolution of 14.4 nm (FWHM). The spectral resolution of each band was 14.4 nm, assigned at the center of the FWHM. The CZMIL system used Teledyne’s HydroFusion software suite to perform radiometric calibration processing, which includes conversion of Digital Number (DN) to at-sensor radiance, glint/ripple correction to smooth water surface effects, and conduct atmospheric corrections which convert surface radiance to surface reflectance, prior to mosaicking into 1 m spatial resolution imagery [25,26,28]. The true-color three band (RGB) collection returned higher spatial resolution imagery (approximately 5 cm ground resolution for 400 m altitude collection), which were used as ancillary data for visual accuracy and feature verification [25].

The integrated system, which was developed specifically for use in coastal mapping and charting activities, was designed to meet or exceed U.S. Geological Survey Quality Level 2 (QL2) standards for topographic LiDAR data collections, and produces a vertical accuracy ≤ 10 cm RMSEZ on land surfaces and shallow water [24]. Vertical positions from the LiDAR data were referenced to the NAD83 ellipsoid (GRS 1980) and were converted to NAVD88 orthometric heights using geoid model 12B. Seamless topobathymetric first return and bare-earth digital elevation model (DEM) grids (1 m spot spacing and approximately 19.5 cm elevation accuracy at a 95% confidence level) were produced through CZMIL processing [26,28,29].

Table 1. Summary of CZMIL collection type, date, location and purpose in the North Carolina Outer Banks from 2016 to 2019.

Data Type	Collection Date	Location	Purpose
Hyperspectral (1 m), LiDAR (1 m), RGB (5 cm)	22 November 2016	Pea Island, NC, USA	Post-hurricane Matthew
Hyperspectral (1 m), LiDAR (1 m), RGB (5 cm)	15 October 2018	Pea Island, NC, USA	Post-hurricane Florence
Hyperspectral (1 m), LiDAR (1 m), RGB (5 cm)	1 October 2019	DUNEX Extent	DUNEX

provides a summary of all qualifying NCMP data from the late-growing season (2016 NCMP data were collected several weeks past the typical growing season for North Carolina) [30] across the study period (2016 to 2019). These data were used to develop DuneVeg tools [31] to semi-automate vegetation metric extraction and to demonstrate how derived metrics can be used to evaluate and compare spatial and temporal changes in key vegetation properties.

2.2.2. Dune Vegetation Metric Extraction

A suite of custom geospatial processing tools (DuneVeg Toolbox, hereafter) was developed to automate the extraction of metrics from the NCMP data [31]. Using standard NCMP imagery and LiDAR products enable workflow automation to extract prioritized dune vegetation metrics in an efficient and repeatable way. The software used to develop the methodological framework for process automation was ENVI version 5.5 (L3Harris Geospatial Solutions, Inc., Broomfield, CO, USA). The ENVI Application Programming Interface (API) in IDL was used to customize ENVI functionality for automating metric extraction processes and the ENVI Py for ArcGIS module allowed for integrated access using ArcGIS Pro (version 2.3) software (ESRI, Redlands, CA, USA). The software integration allows extended capabilities that make use of spatial analytics and spectral band math algorithms to derive the vegetation metrics.

Table 2. Summary of dune vegetation metrics generated from the DuneVeg Toolbox used for conducting multi-temporal trend analysis. Included are National Coastal Mapping Program (NCMP) collection dates (year) used for each metric.

Metric	NCMP Data Product	NCMP Collection Year
Normalized Difference Vegetation Index (NDVI)	HSI	2016, 2018, 2019
Vegetation Cover (VC)	HSI	2016, 2018, 2019
Leaf Area Index (LAI)	HSI	2016, 2018, 2019
Digital Elevation Model (DEM)	LiDAR	2016, 2018, 2019
Digital Surface Model (DSM)	LiDAR	2016, 2018, 2019
Canopy Height Model (CHM)	LiDAR	2016, 2018, 2019

provides a list of NCMP data used in conjunction with the DuneVeg Toolbox to generate elevation and vegetation metrics.

2.2.3. Vegetation Metric: Normalized Difference Vegetation Index (NDVI)

The normalized difference vegetation index (NDVI) provides an efficient means of monitoring changes in vegetation biomass (relative biomass of live, green plants) and health. Previous studies have used NDVI data to evaluate vegetation impacts and recovery from anthropogenic, natural, and invasive disturbances [8,32,33]. Chlorophyll, the pigment in plant leaves that is used in photosynthesis, strongly absorbs visible light (from 400 to 700 nm). Conversely, the cell structure of the leaves strongly reflects near-infrared (NIR) light (from 700 to 1100 nm). The more leaves a plant has, the more these wavelengths of light are affected [34]. Previous studies have demonstrated the use of high spectral resolution HSI-derived NDVI for vegetation monitoring [35,36,37] and have shown optimal contrast between chlorophyll pigment absorption (red band) and the high reflectivity of plant material (NIR band). Moreover, when leaves become stressed or otherwise diseased, they turn yellow and as such reflect significantly less in the NIR range. Because of these photosynthetic influences, the red and NIR bands are standard for calculating NDVI [38], as shown in Equation (1):

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}$$

(1)

where NDVI is the normalized difference vegetation index, NIR is HSI band 35 (856 nm), the near infrared spectral band center value, and Red is HSI band 21 (656 nm), the red spectral band center value. Band assignments were based on the wavelength closest to the band center within each spectral range according to ENVI software functionality and spectral index requirements (L3Harris Geospatial Solutions, Inc.). NDVI values range from −1 to 1, where lower values (closer to −1) generally represent non-vegetated features and higher values (closer to 1) tend to represent healthy, green vegetation [39]. Therefore, the higher the NDVI value, the higher the biomass, productivity, and vigor inferred from these calculations.

2.2.4. Vegetation Metric: Density Estimation (Vegetation Cover)

Delineation of temporal dune vegetation presence and extraction of vegetation density estimates from the atmospherically corrected NCMP HSI, were performed using the ENVI Spectral Processing Exploitation and Analysis Resource (SPEAR) vegetation delineation tool. The SPEAR tool performs a density slice method based on NDVI-derived images to generate density categories. This method was used since NDVI has a strong correlation to aboveground biomass estimation [39,40]. Therefore, inferences can be made about the level of vigor and vegetation density using NDVI value ranges. The density categories estimated from this method were (1) dense vegetation (NDVI = 0.70–1.00), (2) moderate vegetation (NDVI = 0.50–0.69), (3) sparse vegetation (NDVI = 0.25–0.49), and (4) non-vegetation (NDVI = −1.0–0.24).

The output format for this analysis was ESRI polygon features (shapefiles) indicative of the NDVI threshold estimates from each vegetation density category. Once the output densities were derived, a spatial merge (geoprocessing task) was used to join the three vegetation densities (dense, moderate, and sparse) into one category representing all vegetation (polygon feature layer). Subsequently, the output of the merge operation represented vegetation presence for the image scene. Likewise, the non-vegetation layer represented vegetation absence. Wetland, submerged, or inundated areas appearing to have brown or dying (non-green) vegetation were included as non-vegetation for purposes of deriving vegetation density estimates. The resultant products of this analysis were used as mask layers for subsequent vegetation metric extraction.

2.2.5. Vegetation Metric: Leaf Area Index (LAI)

A commonly used variable for quantifying vegetation canopy is the leaf area index (LAI), which is a measure of green leaf area (m²) per unit of ground surface area (m²). In basic terms, it is the estimated (image-derived) or measured (field-derived) amount of foliage in a plant canopy. Leaf surfaces are the primary driver of energy and mass exchange. Therefore, important processes such as canopy interception, evapotranspiration, and photosynthesis are directly proportional to LAI [41]. Moreover, the amount of leaf area directly influences net primary productivity and biomass potential, and in many systems aboveground biomass can be used to estimate belowground biomass [42].

For this study, LAI was estimated from optical remote methods using the NCMP HSI and select spectral bands. The calculations were made using ENVI band algebra and a spectral indices tool specifically designed for LAI. The equation used for calculating this vegetation index is shown in Equation (2) [43]:

$$\mathrm{L}\mathrm{A}\mathrm{I}=3.618\times \mathrm{E}\mathrm{V}\mathrm{I}-0.118$$

(2)

where EVI is the Enhanced Vegetation Index and is an important variable for estimating the LAI using remote sensing methods, especially in areas with high relative biomass. EVI is a broadband greenness index, which was calculated using Equation (3) [44]:

$$\mathrm{E}\mathrm{V}\mathrm{I}=\mathrm{G}\times \frac{\left(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}\right)}{(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{C}1\times \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{C}2\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+\mathrm{L})}$$

(3)

where G is the gain factor (2.5), NIR is HSI band 35, Red is HSI band 21, Blue is HSI band 8 (471 nm, the blue spectral band center value; all spectral bands were atmospherically corrected surface reflectance values), and the numerical coefficients (C1, C2) and canopy background adjustment (L) that ENVI uses as part of the EVI equation. The EVI is useful for optimizing the vegetation signal in areas of high LAI (very dense vegetation), and where NDVI may saturate, which can be a limitation of using NDVI alone [45]. This algorithm makes use of the blue band to correct for soil background signals and to reduce atmospheric influences such as aerosol scattering (L3Harris Geospatial Solutions, Inc., Boulder, CO, USA).

2.2.6. Vegetation Metric: Canopy Height Model

The LiDAR data were processed to generate gridded models of the topography and surface features. A 1 m DEM was created from points classified as ground returns from the classified LAS files. A DEM is a regularly spaced raster grid of elevation values of a surface terrain with non-ground points removed so the output is a surface representation of bare ground void of vegetation and other aboveground features. A Triangulated Irregular Network (TIN) interpolation method was used with rural area filtering to produce the DEM. Similarly, a 1 m Digital Surface Model (DSM) was created from points classified as non-ground using the first return points from the classified LAS file. A DSM is a regularly spaced raster grid of surface elevation values that includes vegetation and other aboveground features. The LAS files were processed using an interpolation method in ENVI that is optimized for vegetation analysis and filters points specifically for extracting vegetation heights. A vegetation canopy height model (CHM) was created by performing a subtraction (DSM—DEM) routine using band math in ENVI to derive a height model with 1 m spatial resolution. The non-vegetation polygon feature (from the density estimate calculation) was used as a mask (geoprocessing task) to exclude non-vegetation features. The CHM contains vegetation canopy height for all areas identified as vegetation within the LiDAR extent footprint.

2.2.7. Vegetation Metric: Statistical Analysis

Vertical accuracies of LiDAR-based elevation models were assessed using ground control points (GCPs) surveyed as part of an associated DUNEX study [46]. The field-based elevation measurements on Pea Island were collected using a real-time kinematic (RTK) Global Positioning System at designated island feature (foredune, dune toe, back toe, highest dune, extra dune, and landward limit) locations along transects with 100 m spacings. The mean error and root mean square error (RMSE) were calculated using the LiDAR-derived DEM elevations relative to the RTK ground survey elevations [47]. The RMSE has been previously used to quantify the accuracy of LiDAR-derived DEMs for open and sparsely vegetated coastal environments [48,49] and the mean error is a good indicator of vertical offsets in wetland landscapes [50,51].

In the absence of in situ measurements of vegetation density and biomass, ancillary imagery and visual assessments were used to evaluate accuracies in HSI-derived NDVI (as well as vegetation density and LAI, by inference). The criterion standards for validating density and biomass were constructed through vegetation presence and greenness ratings based on independent visual examination of project HSI and RGB imagery, as well as ancillary aerial photographs (e.g., National Agriculture Imagery Program imagery collected in close temporal proximity to project HSI). To assess the accuracies of the extracted vegetation metrics, the ArcGIS Desktop v10.7.1 Confusion Matrix tool was used to compute user and producer accuracies, derive Cohen’s Kappa statistic for measuring agreement, and generate an overall accuracy score. Producer’s accuracy was used to assess the probability that real on-the-ground features (based on independent visual examination) were correctly estimated using HSI-based extraction methods. User’s accuracy was used to assess the probability that features identified using HSI-based extractions were present on the ground (per individual visual estimates of density and greenness). Cohen’s Kappa was used as a general measure of inter-rater agreement or reliability [52].

Moreover, to attain comparability among vegetation metrics, statistical analyses were performed using Statistical Analysis System software version 9.2. The PROC General Linear Models (GLM) procedure was used to do a one-way analysis of variance (ANOVA) and a means separation test (Tukey’s, α = 0.01) to evaluate the significance of differences between years for each assessment unit.

2.3. Pilot Study: Multi-Temporal Trend Analysis

Storm Impacts

Deriving vegetation metrics that bracket major disturbance events (i.e., tropical storms) provides useful data for assessing episodic impacts and ecosystem resilience and recovery. A pilot study was conducted to demonstrate the use of the derived vegetation metrics for identifying and evaluating differences and trends across a range of storm-induced impacts. For trend analyses, select areas of interest (AOIs, also referred to as assessment units; Figure 1) from within the Pea Island study site were used to evaluate recent changes in elevation and vegetation properties. The AOIs were selected based on coastal and geomorphological changes related to recent storms, i.e., overwash. An overwash event occurs when storm surge overtops the foredune ridge or when the foredune ridge is eroded by prolonged wave energies [53].

Table 3. Hurricanes that made landfall from 2016 to 2019 on or near the Pea Island study site.

Hurricane	Date (Landfall)	Wind Speed (Landfall)	Maximum Surge
Hermine	3 September 2016	20.6 ms−1	1.06 m
Matthew	9 October 2016	21.9 ms−1	1.09 m
Florence	14 September 2018	11.6 ms−1	0.90 m
Michael	12 October 2018	23.2 ms−1	1.65 m
Dorian	6 September 2019	20.1 ms−1	1.33 m

provides a list of major storms that impacted the Pea Island study site between 2016 and 2019. These storms, which ranged in wind speed from 11.6 ms⁻¹ to 23.2 ms⁻¹ at Pea Island, and had a maximum storm surge from 0.9 m to 1.65 m, include Hurricanes Hermine, Matthew, Florence, Michael, and Dorian [23]. Unlike major storms, which typically provide more direct and measurable impacts, these moderate and higher frequency storms can cause more subtle yet extended impacts that can be more difficult to monitor.

Very high spatial resolution airborne true-color imagery and high-resolution LiDAR-based elevation data (described in the Airborne Hyperspectral Imagery and LiDAR Data section) were used to identify hurricane-related overwash areas within the Pea Island site. The AOIs selected for this study consisted of 2 km long bins along island dunes—with widths extending from the dune toe to the furthest landward dune [46]—that experienced moderate or low overwash impacts related to tropical storms between 2016 and 2019 (Figure 1).

Vegetation-based analyses were performed for each AOI (excluding non-vegetated areas) using summary statistics and per-pixel linear regression and hot spot analyses. To evaluate and compare spatial and temporal changes in vegetation properties, the Spatial Analyst Zonal Statistics as Table tool was applied in ArcGIS to generate summary statistics (i.e., mean, min, max, range, standard deviation) for each vegetation metric, by AOI and date combination. Trends analyses were performed using the Curve Fit v10.1 [54] ArcGIS extension to conduct linear regressions across the period of record (2016 to 2019). The Curve Fit extension, which utilized the range of vegetation metric values at each pixel location (regressed over time), was used to generate raster surfaces of regression model parameter estimates, standard errors, and goodness-of-fit for each DuneVeg metric [54]. The per-pixel linear regression data, which can provide spatial representation of vegetation change over time, were used with a Getis–Ord Gi* hot spot analysis to identify areas that have been significantly impacted by disturbance events and/or may be susceptible to future pressures.

3. Results

3.1. DuneVeg Toolbox

3.1.1. Dune Vegetation Metric Products

The DuneVeg Toolbox was developed as a means of automating methodologies and standardizing target vegetation metric and elevation-based products as part of the National Coastal Mapping Program. Figure 2 shows example NCMP data and vegetation metric products from the Pea Island South assessment unit (moderate overwash). Figure 2A,B show the NCMP HSI and LiDAR (represented as digital elevation model) data, and panels C-F show examples of the canopy height model (CHM), normalized difference vegetation index (NDVI), vegetation density (VD), and leaf area index (LAI), within a portion of the Pea Island South AOI. The HSI in Figure 2A shows, from west-to-east, the Pea Island back barrier wetlands, North Carolina Highway 12, back-island dunes, fore-island dunes, and beach, respectively. Figure 2B illustrates an example of Pea Island ground elevation, represented as a DEM (last return), derived from NCMP LiDAR data. Minimum and maximum elevations (NAVD88) occurred along the back barrier wetlands and the fore-island dunes, respectively. The elevations between these extremes were occupied by the back-island dunes and beach features, which are located at moderately high and moderately low elevations, respectively.

Example CHM data, which were generated by subtracting project DEM (last return/ground) from DSM (first return) data, are shown in Figure 2C. The minimum CHM elevations correspond to unvegetated island features (i.e., beach and developed; represented by beige in Figure 2C) and the maximum elevations correspond to shrub species located between the back-island dunes and wetland zones. The locations of these shrub features, represented by the yellow and blue colors in Figure 2C, are corroborated by estuarine scrub/shrub wetlands identified in NOAA’s 10 m Coastal Change Analysis Program (C-CAP) Regional Land Cover data [55].

Figure 2D–F represent the other vegetation metrics that were generated using the DuneVeg Toolbox. Figure 2D shows NDVI, which has a well-established correlation to photosynthetic activity and plant health. Orange and red areas indicate vegetation with low plant vigor, while green and blue areas represent vegetation with higher biomass and vigor (Figure 2D). The lower NDVI areas consist of low biomass beachgrasses, while the higher NDVI areas consist of high biomass shrubs and emergent wetland plants. Figure 2E shows an example of vegetation density data that were generated using density thresholds and the NDVI data. Blue represents high density vegetation, yellow represents moderate density vegetation, and red represents sparse vegetation (Figure 2E). These vegetation density zones correspond to regions of greenness and density observed in the HSI (Figure 2A). LAI, which is an estimate of leaf productivity, is shown in Figure 2F. Red areas represent vegetation with low LAI (low leaf productivity) and greens and blues represent vegetation with high LAI (high leaf productivity). LAI values in Figure 2F correspond well to observed regions of high and low vegetation productivity (greenness and density, Figure 2A).

3.1.2. Elevation and Vegetation Metric Accuracy Assessments

The overall average vertical accuracies of the LiDAR-derived DEMs for all dates and assessment units—expressed as the mean error and RMSE—were −0.006 m (n = 1228) and 0.214 m, respectively. The mean errors for the 2016, 2018, and 2019 DEMs were −0.03 m (n = 421), −0.07 m (n = 400), and 0.08 m (n = 407), respectively. The RMSE values for the 2016, 2018, and 2019 DEMs were 0.146 m, 0.225 m, and 0.257 m, respectively. Across all dates, the average vertical accuracies of DEMs ranged from a low RMSE of 0.15 m (secondary dune features) to a high RMSE of 0.28 m (dune toe). All other features had RMSE ≤ 0.25 m. These derived errors fall within the expected range for LiDAR-based DEMs, and those reported as part of other barrier island studies [56,57].

To assess the accuracy of the HSI-derived NDVI data, a total of 180 assessment points (stratified random sample) were used across all dates and assessment units to construct a confusion matrix (

Table 4. Confusion matrix of vegetation classes derived from hyperspectral imagery.

Observed
HIS Derived	Class	Sparse	Moderate	Dense	Total	User’s acc.
	Sparse	58	4	1	63	0.92
	Moderate	2	51	4	57	0.89
	Dense	0	5	55	60	0.92
	Total	60	60	60	180	0.00
	Producer’s acc.	0.97	0.85	0.92	0.00	0.87

). Accuracy assessments of NDVI returned minimum and maximum accuracies of 0.85 and 0.97 for the moderate and low NDVI ranges, respectively. The confusion matrix returned an overall accuracy and Cohen’s Kappa value of 0.91 and 0.87, respectively. Most misclassifications occurred in the moderate vegetation class. This is due to difficulties in differentiating moderate from sparse or dense vegetation when greenness and biomass approach class thresholds.

3.1.3. Vegetation Metric Summary Statistics

The Pea Island North and South AOIs had Low and Moderate overwash events, respectively.

Table 5. Raster-derived mean area (Total Area), vegetation cover (Veg Cover), elevation, canopy height (CHM), normalized difference vegetation index (NDVI), and leaf area index (LAI) within the Pea Island assessment units.

Location	Overwash	Date	Total Area	Veg Cover	Elevation (m)		CHM (m)		NDVI		LAI (m2/m2)
			Km2	Percent	Mean	Std	Mean	Std	Mean	Std	Mean	Std
North	Low	2016	0.47	10.52	3.91	2.09	0.77	1.21	0.52	0.18	5.07	3.61
North	Low	2018	0.47	57.36	3.85	2.09	0.23	0.61	0.37	0.16	2.94	1.49
North	Low	2019	0.47	22.03	3.84	2.14	0.42	1.03	0.45	0.23	3.08	1.72
South	Moderate	2016	0.27	10.61	4.68	2.16	0.28	0.45	0.45	0.15	4.15	2.32
South	Moderate	2018	0.27	48.29	4.66	2.19	0.17	0.29	0.39	0.14	2.94	1.50
South	Moderate	2019	0.27	23.29	4.69	2.26	0.12	0.36	0.38	0.18	2.53	1.25

provides the overwash designation, along with vegetation metric summary statistics (i.e., means and standard deviations) for Total Area, Vegetation Cover (VC), CHM, NDVI, and LAI, all derived from the DuneVeg Toolbox for the North and South assessment units. Similarly, Figure 3 provides the means and trends of DEM (A), CHM (B), NDVI (C), and LAI (D) values for the Pea Island units. These represent all qualifying data, from late-growing season, across the period of record (2016 to 2019).

The total area for the Pea Island North and South AOIs were 0.47 km² and 0.27 km², respectively (Table 5). The percentage of each AOI that was vegetated at any time ranged from a minimum of 10.52% (North AOI in 2016) to a maximum of 57.36% (North AOI in 2018). NOAA’s C-CAP land cover data show the vegetation in 2016 consisted of 84% emergent and 16% shrub in the North AOI, and 95% emergent and 5% shrub in the South AOI [55]. The summary statistics show all AOIs experienced similar vegetation cover trends, increasing significantly from 2016 to 2018, then experiencing smaller decreases between 2018 and 2019.

The elevation ranged from a low of 3.84 m to a high of 4.69 m (NAVD88), for the North AOI (2016) and South AOI (2019), respectively (Table 5 and Figure 3A). While the mean elevation in the South AOI remained relatively unchanged from 2016 through 2019, the North AOI experienced slight decreases in mean elevations, across the period of analysis. The total mean elevation for the Pea Island site (across all AOIs) remained relatively unchanged from year-to-year, decreasing from 4.29 m in 2016 to 4.25 m 2018, then increasing slightly to 4.26 m by 2019 (Figure 3A). Given the known overwash events within the AOIs, these mean elevations are indicative of sediment that was either reworked within the system or moved out of the system and replaced by sediment introduced via wind and tidal forces.

The range of CHM values for Pea Island was 0.1 m to 20 m. The mean CHM values ranged from a low of 0.12 m for the South AOI in 2019 to a high of 0.77 m for the North AOI in 2016 (Table 5 and Figure 3B). The North and South AOIs both experienced decreasing trends in CHM values across the entire period of analysis (2016 to 2019). The initial loss was driven largely by a reduction in canopy height in the North AOI, which experienced a −0.54 m reduction in mean CHM from 2016 to 2018. The North AOI did experience a minimal increase in mean CHM values from 2018 to 2019, while the South AOI experienced decreases across both time periods (2016–2018 and 2018–2019) (Figure 3B). All AOI-specific mean CHM values were significantly different between years (Tukey’s HSD [honestly significant difference]: p < 0.01).

Individual NDVI values (pixel-wise) ranged from 0 to 1 across all vegetated areas within the Pea Island study site, where low and high values were related to lower and higher levels of vegetation greenness (for vegetated areas), respectively. The mean NDVI values ranged from a low of 0.37 for the North AOI in 2018 to a high of 0.52 for the same AOI in 2016 (Table 5 and Figure 3C). The Pea Island sites experienced significant decreases in mean NDVI values between 2016 and 2018, decreasing from 0.52 to 0.37 and 0.45 to 0.39 for the North and South AOIs, respectfully (Figure 3C). While the North AOI did experience a moderate increase (recovery) in NDVI by 2019, the South AOI underwent a nominal decrease. By comparison, the decade (2006 to 2015) preceding this period of analysis (2016 to 2019) encountered larger but less frequent storms (Hurricane Irene, 2011; 25.9 ms⁻¹ sustained winds and 2.35 m surge) and experienced moderately increasing NDVI trends (based on Landsat-derived NDVI data for the Pea Island study site). This is corroborated by previous studies which have shown vegetation typically recovers by the next full growing season after a major disturbance event [58].

The range of pixel-wise LAI values was 0 to ~10 (m²/m²) for the dune vegetation across all AOIs, where higher values were related to higher leaf area. Mean LAI values ranged from a low of 2.53 for the South AOI in 2019 to a high of 5.07 for the North AOI in 2016. The mean total LAI values for the Pea Island study site decreased from 4.38 in 2016, to 3.44 in 2018, to 3.09 in 2019 (Figure 3D). LAI experienced similar overall trends to those of CHM and NDVI, however, there was much smaller increase (recovery) in LAI in the North AOI in 2019, as compared to CHM and NDVI.

3.2. Multi-Temporal Trend Analysis

3.2.1. Elevation and Vegetation Metrics

Although mean elevation and vegetation attribute values can provide general measures of wetland disturbance and recovery, they can also lack the detail required to evaluate the spatial variability and patterns associated with short- and long-term stressors. Data derived from the DuneVeg Toolbox were used in conjunction with the Curve Fit tool to evaluate changes in elevation and vegetation over space (across areas of low and moderate hurricane impacts) and time (across the period of analysis). Figure 4 provides examples of the rate of change and spatial variability of elevation (Panels A, B) and vegetation productivity (Panels C, D), on a pixel-by-pixel basis. For the Northern AOI (Figure 4A), most elevation changes occurred oceanside along the beach and dune toe, with linear erosion along the northern beach (orange and red colors in Figure 4A) and sediment deposition in the central and southern portions (blue color in Figure 4A) of the AOI. The erosion in the Southern AOI was non-linear and indicative of sediment erosion and distribution induced by an overwash event. This is evidenced by the fragmented areas of erosion (orange and red in Figure 4B) throughout the northern and southern dunes in the Southern AOI. The central region of the South AOI was not overwashed and therefore experienced increasing rates of elevation (greens and blue in Figure 4B) as a result of sediment transport and deposition. The overwash events within the Southern AOI resulted in higher rates of accretion along the inter- and back-dunes (green color in Figure 4B) due to sediment transport beyond the primary dune.

Figure 4 also illustrates the vegetation biomass (NDVI) rates of change and spatial patterns on a per-pixel basis within the Pea Island AOIs. The NDVI rates of change are color-ramped (Figure 4, Panels C and D), where the orange-to-red colors represent decreasing rates of NDVI, and the yellow-to-green colors represent increasing rates. The Pea Island North AOI (Panel C) and South AOI (Panel D) encountered similar patterns of change, with the majority of each unit experiencing slight decreasing rates of NDVI (represented by a light orange color), while smaller regions experienced higher decreasing (linear features to the north, represented by darker shades of red) and higher increasing (disaggregated regions to the south, represented by darker shades of green) rates of NDVI. The areas of distinct decreasing NDVI were primarily along the beach and foredune features in each AOI while the areas of increasing NDVI occurred primarily at leeward positions within inter- or along secondary-dune features (Figure 4C,D). These change rates are consistent with the NDVI summary statistics for Pea Island, which showed an overall decreasing trend in NDVI within the period of analysis (Table 5 and Figure 3). The vegetation-based regression data (Figure 4C,D) correspond to the elevation data (Figure 4A,B), where areas of distinct reductions in NDVI (along dune features) coincide with distinct reductions in elevation. Likewise, areas with moderate increasing NDVI (interdune features) coincide with areas that experienced moderate elevation increases. The areas of distinct reductions are indicative of direct hurricane impacts (caused by wind and wave energies) in barrier island systems. The areas of increasing NDVI are indicative of indirect hurricane impacts in wetland systems, where overwashed sediments can result in elevation benefits, including flood stress reduction, and increases in productivity and aboveground biomass [59]. Similar trend analysis results were observed with LAI and CHM data, and are therefore not presented here, for brevity.

3.2.2. Hot Spot Analysis: Vegetation Metrics

The Getis–Ord Gi* hot spot analysis was performed for each vegetation metric (CHM, LAI, NDVI). This process identifies pixels with significant values and assigns them to associated Gi* bins based on p-value confidence levels (

Table 6. Hot spot analysis of vegetation metrics using the Getis–Ord Gi* and corresponding z-scores and p-values.

Cluster Type	CHM		LAI		NDVI		Multi-Metric		Gi* Bin	Z-Score	p-Value
	Number	Area (ha)	Number	Area (ha)	Number	Area (ha)	Number	Area (ha)	(Single, Combined)
Hot Spot	251	2.51	149	1.49	274	2.74	13	0.13	(−3, −7–−9)	−4.97	p < 0.001
Hot Spot	128	1.28	74	0.74	128	1.28	141	1.41	(−2,− 4–−6)	−2.27	p < 0.01
Hot Spot	97	0.97	62	0.62	119	1.19	778	7.78	(−1, −1–−3)	−1.81	p < 0.05
Not Significant	7167	71.67	7554	75.54	6629	66.29	6152	61.52	(0, 0, 0)	0.02
Cold Spot	35	0.35	17	0.17	157	1.57	704	7.04	(1, 1–3)	1.81	p < 0.05
Cold Spot	53	0.53	39	0.39	191	1.91	140	1.40	(2, 4–6)	2.26	p < 0.01
Cold Spot	205	2.05	41	0.41	438	4.38	8	0.08	(3, 7–9)	5.73	p < 0.001

). Bins for single metric analyses are determined at the p < 0.05, 0.01, and 0.001 levels, and assigned values of 1, 2, 3, and −1, −2, −3, for the Cold and Hot Spots, respectively. To perform a more comprehensive evaluation of vegetation change, the output layers from each individual hot spot analysis (CHM, LAI, and NDVI) were overlayed (summed) to provide a more robust vegetation hot spot assessment. This combined approach resulted in Gi* bins that ranged from 1 to 9 and −1 to −9, for cold and hot spots, respectively. Table 6 shows the results of the Getis–Ord Gi* hot spot analyses. Single metric assessments show only small portions of each AOI were classified as hot or cold spots, with a majority of those occurring at the p < 0.001 designations (3 and −3). For the multi-metric assessment, most of the hot and cold spot regions occurred at the p < 0.05 designation. This indicates most pixels that experienced significant vegetation-based rates of change either experienced those at the p < 0.001 level for only one vegetation metric, or at the p < 0.01 and/or 0.05 levels for multiple metrics.

Figure 5 shows the overall vegetation hot spot analysis, where the hot spot bins are represented by light orange (−1, −2, −3), dark orange (−4, −5, −6), and dark red (−7, −8, −9), respectively; while the cold spot bins are represented by light green (1, 2, 3), medium green (4, 5, 6), and dark green (7, 8, 9), respectively. Most of the hot spot areas (orange and red pixels, Figure 5) are concentrated within linear features at the northern reach of each Pea Island AOI or scattered throughout the southern portions of the AOIs. These hot spot areas represent regions within Pea Island where vegetation biomass and health experienced significant reductions across the entire period of analysis (2016 to 2019). The cold spot areas (Figure 5, blue pixels), which represent regions where vegetation biomass and health experienced significant increases across the period of analysis, are found primarily leeward of the primary dune along southern portions of the AOIs, or intermittently scattered throughout central portions of the AOIs.

4. Discussion

Dune vegetation plays a crucial role in the processes that control erosion, dune stability, and overall coastal resilience. Therefore, it is important for coastal managers to rapidly and accurately evaluate coastal vegetation properties and dynamics, especially in the face of current and future climate change trends. However, collecting or generating vegetation data that are useful for ecosystem management, or as inputs to coastal numerical models, can be challenging in coastal and island systems. Even when traditional data collections or monitoring techniques can be applied, they often lack the detail (i.e., spatial and temporal) or coverage required to evaluate finer scale changes, especially in complex and dynamic landscapes. This research effort demonstrated the use of DuneVeg Tool generated elevation and vegetation metric data to assess hurricane-induced trends and changes in coastal morphology and habitat. Summary statistics were shown to provide general elevation and vegetation trends across the period of analysis, but per-pixel regression and hot spot analyses provide the spatial detail necessary for assessing areas or features of concern. Although the period of analysis used to demonstrate the utility of the developed tools and metrics is relatively short (evaluating short-term and direct impacts from hurricanes), the application should be equally useful in assessing indirect impacts over longer time periods.

The DuneVeg tools and metrics may likewise prove useful in coastal numerical modeling. Models rely greatly on estimates of vegetation properties [60]. For instance, wave attenuation through vegetation is a highly variable function of not only hydrodynamics, but also of general vegetation characteristics. Recent studies have updated models that account for wave interactions using spatially explicit and variable vegetation characteristics as inputs [12]. However, model input data are limited, particularly for large scale or regional coastal projects, and therefore many models currently incorporate synthetic vegetation data or constant values as inputs. The vegetation data gaps that exist for many coastal numerical models could be narrowed by way of the vegetation metrics generated through the DuneVeg Tool. Future studies should evaluate the use of these vegetation metrics to refine numerical models and quantify the performance and accuracy of those models relative to existing data model inputs.

5. Conclusions

This study advances existing dune and coastal wetland system knowledge involving remote sensing techniques for classifying, quantifying, and estimating elevation and critical dune vegetation metrics. It also describes the use of summary statistics and regression analyses to evaluate spatial and temporal changes in vegetation properties. Future studies should evaluate the ability of the derived methods and metrics to quantify and monitor long-term indirect changes in coastal wetland vegetation. Future assessments should also investigate the power of fusing DuneVeg Toolbox data (combining geomorphology and environmental characteristics) for a holistic analysis paradigm (as opposed to individual assessments). The ability to combine and analyze spectral and LiDAR data together opens possibilities to explore landscapes as mosaics [61,62]. Still, the methodology and products derived as part of this study can aid researchers and managers in improving predictions of future disturbance, which would help in the advancement of coastal wetland monitoring, restoration, and adaptive management.