Mapping Seafloor Relative Reflectance and Assessing Coral Reef Morphology with EAARL-B Topobathymetric Lidar Waveforms

Topobathymetric lidar is becoming an increasingly valuable tool for benthic habitat mapping, enabling safe, efficient data acquisition over coral reefs and other fragile ecosystems. In 2014, a novel topobathymetric lidar system, the Experimental Advanced Airborne Research Lidar-B (EAARL-B), was used to acquire data in priority habitat areas in the U.S. Virgin Islands (USVI), spanning the 0–44-m depth range. In this study, new algorithms and procedures were developed for generating seafloor relative reflectance, along with a suite of shape-based waveform features from EAARL-B. Waveform features were then correlated with percent cover of coral morphologies, domed and branched, and total cover of hard and soft corals. Results show that the EAARL-B can be used to produce useful seafloor relative reflectance mosaics and also that the additional waveform shape-based features contain additional information that may benefit habitat classification—specifically, to aid in distinguishing among hard corals and their coral morphologies, domed and branched. Knowing the spatial extent of changes in coral communities is important to the understanding of resiliency of coral reefs under stress from human impacts.


Introduction
Coral reef ecosystems provide essential ecosystem services to millions of people around the world (Hughes et al. 2017). The habitat and food provided by these reefs are essential for many coastal communities, including the US Caribbean. In the U.S.
Virgin Islands (USVI), these ecosystems are estimated to provide economic benefit over $187 million annually to the local economy, by supporting tourism, providing coastal protection from storms, and providing habitat for commercially important fisheries (Van Beukering et al. 2011). Over the last several decades, coral reefs have undergone an unprecedented rate of Communicated by Brian B. Barnes decline in the USVI and worldwide (Gardner et al. 2003;Pandolfi et al. 2003;Bellwood et al. 2004). It is likely this decline will endure as these ecosystems continue to be impacted by anthropogenic stressors from ocean warming and acidification (Hughes et al. 2017). In the face of these threats, the ability to map and characterize coral reef ecosystems is a critical tool for coastal managers to monitor these rapid changes on a broad scale (Mumby and Harborne 1999;Brown et al. 2011;Monaco et al. 2012). Such maps would facilitate policies designed to improve the resiliency of these important ecosystems and to sustain the value services these ecosystems provide for coastal communities.
Benthic habitat maps that depict the spatial extent and distribution of coral reefs and other seafloor habitats are valuable to coastal management and policy makers in managing coastal ecosystems and assessing change over time. Mapping of these habitats using divers is infeasible, due to inability of divers to access dangerous or challenging locations and to the time it would take to create maps of sufficient spatial extent. While acoustic techniques are most effective in temperate ecosystems or in deeper waters, airborne bathymetric lidar has increasingly gained recognition as a viable technology for benthic habitat mapping and characterization (Collin et al. 2008;Narayanan et al. 2009;Costa et al. 2009). Past habitat mapping efforts have used lidar for consistent classification of broad functional groups (seagrass, coral, etc.). However, the development of topo-bathymetric lidar systems that record waveform metrics presents an opportunity to explore their use for finer classification of coral reef communities. Linking lidar waveforms metrics to biological characteristics of coral reef habitats and the seafloor may provide new or unique information that will help to capture fundamental changes in these habitats. This may provide another tool to better determine optimal sites for species restoration projects or to focus limited resources on areas that may be of national or conservation value.
This study included developing and testing procedures for generating relative reflectance mosaics and additional waveform features, including area under the curve, skewness, and standard deviation, from the Experimental Advanced Airborne Research Lidar-B (EAARL-B) to benefit mapping of coral reef habitats. The EAARL-B is well suited to acquire spatially dense data in the depth ranges of interest for benthic habitat mapping . However, the EAARL-B system and its processing software, ALPS (Airborne Lidar Processing System) (Nagle and Wright 2016), did not previously (prior to this study) provide functionality for generating seafloor reflectance products or bottom return waveform shape-based features. Reflectance mosaics may greatly enhance the value of the system to benefit benthic habitat mapping, while the addition of other shape-based waveform features may further facilitate assessment of coral communities and benthic composition (Collin et al. 2011;Rogers et al. 2015). The procedures were implemented and investigated using EAARL-B data collected over two priority locations in the USVI in 2014 and assessed using in situ seafloor reflectance spectra collected from a small boat. Correlations between waveform metrics and coral reef communities were performed using 100-m 2 photo mosaics of sites surrounding Flat Cay.

Study Locations
Two study locations in the USVI (Fig. 1) were selected for this research, as they had a diverse range of seafloor habitat types, bottom complexity, and bathymetric relief. The first location was relatively small (25 km 2 ) and included the Buck Island Reef National Monument (BIRNM), north of St. Croix. The second location, situated due south of St. Thomas, was larger (120 km 2 ) and included Flat Cay island and surrounding waters. These locations are of high interest to marine conservation managers because they include several federal and territorial marine protected areas, including Virgin Islands National Park (VINP), St. Thomas East End Reserve (STEER), and Cas Cay-Mangrove Lagoon Marine Reserve and Wildlife Sanctuary (CCMLMR). These areas provide a home to several species protected under the federal Endangered Species Act, including Hawksbill Turtles (Eretmochelys imbricata), Elkhorn Coral (Acropora palmata), and Staghorn Coral (Acropora cervicornis) (Pittman et al. 2008, NOAA OPR 2019. For these reasons, these locations have been studied intensively in the past, and as a result, they contain a variety of remotely sensed and in situ data sets that could help support the research described here. An additional benefit of the two locations is that the data for Buck Island were useful for testing (specifically, for comparisons of the output relative reflectance mosaics against in situ reflectance spectra), while Flat Cay was used to test if waveform metrics could be used to distinguish or characterize coral reef communities.

Data Collection
Bathymetric lidar data were collected by the USGS using the EAARL-B. The data were collected on 11 separate days between March 7 and March 21, 2014. The airborne data acquisition parameters are listed in Table 1. Processing of lidar point clouds and digital elevation models of the region was conducted by United States Geological Survey (USGS), as described in Fredericks et al. (2015).
To process the lidar data, we adopted the definitions of lidar radiometric processing levels given in Kashani et al. (2015), wherein level 0 = raw intensity; level 1 = intensity correction (i.e., correction for range, angle of incidence); level 2 = intensity normalization (i.e., histogram normalization to match adjacent flight strips or data collected across different days, sites, following the level 1 processing); and level 3 = full, rigorous radiometric correction and calibration to obtain "true" surface reflectance (generally unattainable, due to lack of manufacturer-proprietary system information and full environmental characterization). With reference to these processing levels, seafloor relative reflectance, as defined in this study, is a level 2 product, whereas true reflectance corresponds to level 3.
The reference data for testing the seafloor relative reflectance mosaics consisted of underwater spectral reflectance measurements acquired from a small boat in July, 2012, at Buck Island (Fig. 2). These reference underwater spectral reflectance measurements (see Pe'eri et al. 2013, for details on data collection) were acquired for assessing the lidar-derived relative reflectance. While this fieldwork was performed specifically for this experiment, delays in fielding the EAARL-B led to the~20-month gap between the field and airborne data collections. A consequence of the temporal offset between the field and airborne data collection is that some change in benthic habitat type occurred. Specifically, some of the seagrass bed boundaries were observed to have changed. These seagrass bed distribution changes were relatively easy to identify using our own aerial imagery collected with the EAARL-B lidar data, Google Earth imagery, and Esri World Imagery. Reference spectra collected in areas of seagrass bed migration were removed from the analysis, as well as a few additional spectra that were collected outside of the EAARL-B coverage extents, which were not precisely known at the time of the fieldwork. The number of remaining points (14) was fewer than desired but spanned a range of seafloor reflectance values and habitat types (including seagrass, coral, and sand) in the 1-10-m depth range (Fig. 3).
The ground truth for the reef cover characterization consisted of generating 100-m 2 photo mosaics of coral reef communities around Flat Cay. Underwater video footage of 9 sites, ranging in depth from~2 to~17 m, was collected by swimming in a lawnmower pattern along transects placed on the seafloor between September 4 and 9, 2016. Overlapping still frames were extracted from the video and stitched together into a single composite image using texture based video mosaic Gu and Rzhanov 2006). To create a species map for each dive site, each mosaic was georeferenced and viewed on a high-resolution computer screen. Corals and macroalgae were identified and manually segmented down to the lowest possible taxonomic level, typically genus or species. Adobe Photoshop's Magic Wand tool was used to isolate each coral head and patch of macroalgae and mask them with a species-specific color. These masks were used to calculate percent cover. To account for difficulties in identification to species level (especially among octocoral genera) and the effects of less-abundant species, we created five functional groups. Hard corals were divided into "domed" and "branched" groups based on colony growth form in order to represent the varied habitat types they provide and expected differences in response to lidar signals. The domed coral group includes boulder, brain, hill, pillar, and star corals, while the branched coral group includes finger, fire, and staghorn corals. The remaining three groups comprise the octocorals, sponges, and macroalgae turf (mostly genus Dictyota).

Signal Processing
The bathymetric lidar equation, which is given in various forms in the published literature (e.g., Wang and Philpot 2002;Tuell and Park 2004;Collin et al. 2008;Narayanan et al. 2009;Tuell and Carr 2013), relates the received optical power for a laser return pulse to parameters related to the lidar system, the acquisition geometry, and the environment. While various formulations differ slightly, a general form is as follows: where P R is the received optical power; P T is the transmitted Fig. 2 Acquisition of underwater reflectance spectra, which served as reference data for assessing the lidar-derived seafloor relative reflectance mosaics. Left: camera frame and spectrometer probe on deck of boat. Middle: deploying camera frame over the side (note: fiber optic cable connects spectrometer probe to the instrument, which remains on vessel).
Right: underwater image showing one of the seagrass bed sites, with the white reference panel rotated out of the spectrometer's field of view. The quadrat at the base of the camera frame is 0.3 m × 0.3 m, and the graduations are 2 cm power; η is the system optical efficiency factor; ρ is the reflectance of bottom; F p is the loss due to insufficient FOV; A r is the effective area of the receiver optics; θ is the off-nadir transmit angle; n w is the refractive index of the water; H is the altitude of the lidar above the water; D is the water depth; n(s, ω 0 , θ) is the pulse stretching factor; s is the scattering coefficient; ω 0 is the single scattering albedo; K is the diffuse attenuation coefficient of the water; and ϕ is the off-nadir angle of the lidar beam after refraction at the air-water interface. Note that the term D sec ϕ is the slant range of the laser pulse from the water surface to seafloor. For all wavelength-dependent parameters, such as ρ, K, and ω 0 , it is understood that the wavelength, λ, at which the parameter is evaluated is that of the lidar system, which is 532 nm for nearly all current bathymetric lidar systems (Guenther 2007). Typically, it is not possible to directly solve Eq. 1 for ρ, the reflectivity of the seafloor at the laser wavelength (532 nm), due to unknown lidar system parameters (e.g., P T , η, F P ) and environmental parameters (e.g., s, ω 0 , θ, K, n). It is for this reason that most studies of topo-bathymetric lidar reflectance mapping, including this one, aim to produce relative reflectance, rather than "true" or "absolute" seafloor reflectance. In this work, a data-driven approach was taken to derive corrections to lidar bottom intensity (i.e., the data itself was used to drive the determination of correction coefficients) to obtain seafloor relative reflectance. The full, end-to-end workflow for generating the relative reflectance mosaics is depicted in Fig. 4. It is important to note that this workflow was designed to be applied to very large data sets, covering up to hundreds of square kilometers of the seafloor and encompassing tens of millions of lidar points collected over a period of several days to weeks. Therefore, key considerations in developing the workflow included the following: (1) reducing human operator time, (2) reducing computer processing time, and (3) reducing seamline artifacts at the junctions of flight lines and acquisition dates.
The input to the relative reflectance mapping process depicted in Fig. 4 consisted of georeferenced EAARL-B lidar point clouds created with the USGS ALPS software. (Readers interested in the details of this step and the algorithms implemented in ALPS are referred to Nagle and Wright et al. (2016).) Importantly for this work, each bottom return lidar Fig. 3 Locations of underwater reflectance spectra (green triangles) that served as the reference data set for evaluating the EAARL-B relative reflectance mosaic. Note that the EAARL-B coverage extends to the northeast, as shown in Fig. 1, but the seafloor reflectance spectra were limited to this region of the site. Fortunately, however, these sites include a range of different habitat types and depths. The inset (locator map) shows the location of the project site in the Caribbean and in relation to the U.S. East Coast. The red line indicates the southern border of the project site. point in each point cloud had an intensity value, I, which was taken to be the peak amplitude of the detected bottom return and, in turn, proportional to the received optical power. The first step in our procedure was to perform pre-processing or cleaning of the lidar data set, which entailed removing areas of land, as well as obvious noise points. The next step was to apply intensity corrections, corresponding to level 1 processing, as defined in Kashani et al. (2015). Corrections were applied for the following: (1) depth (or, perhaps more precisely stated, for the range-dependent attenuation of radiant flux in the water column) and (2) angle of incidence. The depth correction was obtained by first considering a simplified form of the bathymetric lidar equation (Eq. 1), adapted from Guenther (1985): The form of the lidar equation given in Eq. 2 is based on the simplifying assumptions that (a) the above-water flying height is relatively large compared to the depth, (b) pulse stretching can be neglected, and (c) other unknown parameters in Eq. 1 can be assumed constant and combined into a constant system-loss term, W. From Eq. 2 (and with the further assumptions that the transmit power and system losses are constant throughout the data collection), it can be seen that the natural log of the bottom reflectance is linearly related to depth (or to slant range through the water column), with linear parameters that are functions of the diffuse attenuation coefficient and seabed reflectance. As the goal of the depth correction was to remove the depth dependence to obtain better estimates of ρ, we performed a linear regression of uncorrected intensities on depth for areas of constant bottom type and water clarity. The linear transformation parameters were then used in the correction to remove the depth dependence. Figure 5 shows a heatmap created using points from entire day of data collection with the resulting linear best-fit line. The color scale can be interpreted as "hotter" regions being those of high point concentration. This type of heatmap was found very beneficial in this study for visualizing trends in large volumes of data.
Using the parameters of the linear fit depicted in Fig. 5, the depth correction is given by the following: Fig. 4 Processing workflow used in this study. The color scheme is as follows: green denotes an input to a process; yellow denotes an output of a process (i.e., a final, processed data product); gray denotes an intermediate step; and orange denotes the use of these products to meet science objectives where I′ is the corrected intensity, I is the input (uncorrected) intensity, a and b are the coefficients of the linear fit described above, D is water depth, and ϕ is the off-nadir angle of the laser beam in the water. As before, the product D sec ϕ is the slant range of the laser pulse from the water surface to seafloor. The next step in the process was the incidence angle correction. This correction is extremely important, since, unlike other bathymetric lidar systems, such as the Optech CZMIL (Feygels et al. 2013) and SHOALS (Collin et al. 2008), the EAARL-B does not attempt to maintain a constant off-nadir transmit pulse angle but instead scans back and forth across the field of view, passing nearly through nadir. This created a pronounced reduction in intensity towards the outer edges of the swath. Fig. 6 Heatmap scatterplot of depth-corrected intensity, I′, and incidence angle, along with fitted curve (blue line) Fig. 7 Seafloor relative reflectance mosaic for the 120-km 2 St. Thomas site The incidence angle correction was computed in a data-driven approach, similar to the depth correction. The form of the incidence angle correction, based on the Phong reflectance model (Phong 1975;Jutzi and Gross 2009;Hasegawa 2006) is: In a manner similar to the depth correction, the parameters α and β (which relate to the specular or non-Lambertian nature of the seafloor) were determined empirically through a curve-fitting procedure (Fig. 6).
An underlying assumption in the correction procedures described above is that the points used to derive the correction parameters had the same "true" reflectance; hence, it was important that the subsets of points used as input were collected from a homogeneous bottom type. When and where possible, homogeneous regions were identified with the aid of imagery and/or existing habitat maps. For the EAARL-B deep receiver channel, an initial approximation of correction parameters was made using all of the points calculated for given day. Then, the resulting point cloud was used to assist in delineation of uniform bottom type, typically sand. The points of uniform bottom were then used to determine final correction parameters.
Continuing with the workflow depicted in Fig. 4, a normalization step was next performed, corresponding to level 2 Fig. 8 Before-and-after images, showing the results of the processing performed to generate seafloor relative reflectance from the raw lidar intensity data. Top: raw intensity, which is the peak amplitude of the detected bottom return. Bottom: seafloor relative reflectance generated through our processing procedures. Salient artifacts have been greatly reduced or even eliminated. The dark feature visible near the middle of the area, which is visible in both the raw intensity data and the seafloor reflectance mosaic, was identified as a wreck in the largest-scale NOAA Nautical chart of the area: Chart 25649-1, St. Thomas Harbor processing, as defined in Kashani et al. (2015). This step consisted of first matching points from overlapping point clouds within 1 m of each other. The distributions of the corrected intensities of the matched points were analyzed, and a linear transformation (shifting and scaling of the intensities) was performed on the second point cloud, such that its mean and standard deviation were made to equal that of the first point cloud.
Next, the level 2 intensities were interpolated to a regular grid. Although any of a number of interpolation algorithms could have been used in this step, based on experimentation, we used inverse distance weighting (IDW), which was found to reduce seamline artifacts between adjacent flight lines and to generally create a more uniform representation of regions in which the angle of incidence correction has been either overor under-applied, while also keeping processing times within practical limits.
Next, a second histogram normalization was performed, such that the level 2 intensity rasters could be combined for multiple flight lines and days, while minimizing seamlines. This was achieved using custom software developed as part of this research. This software applies semi-automated histogram scaling and shifting to adjust the contrast and brightness across adjacent rasters using a graphical user interface (GUI). Overlapping gridded data were adjusted by the user until overlapping regions visually matched. The output relative reflectance mosaic was generated in Esri ASCII raster format, with values linearly scaled to the 0-255 range (i.e., 8-bit rasters), for compatibility with other coastal GIS data layers. The resulting relative reflectance mosaics were then assessed visually and through quantitative comparison with the reference spectra acquired in the Buck Island site.
The final signal processing step was to generate additional features (metrics) that relate to the shape of the identified bottom return in the EAARL-B lidar waveform. The specific waveform features generated in this study and further described in Parrish et al. (2014) were as follows: (1) standard deviation (a measure of the width or "spread" of the bottom return pulse), (2) area under the curve (a measure of the total energy returned from the bottom), and (3) skewness (a measure of the asymmetry of the bottom return waveform). The key idea underlying the use of these features is that coral, seagrass, or other cover types are theorized to modify the shape of the bottom return waveform; hence, these features that relate to the bottom return waveform shape could be useful predictors of cover type. Further description of and justification for the use of these types of shape-based waveform features for seabed habitat analysis is provided by Collin et al. (2008). The three features are computed as described in Parrish et al. (2014): Distance along cross-secƟon (m) Distance along cross-secƟon (m) Fig. 9 Example cross section taken perpendicular to flight paths, colored by collection date. The top cross section shows raw intensities, while the bottom cross section shows the same points after depth and angle of incidence corrections and histogram normalization where N is the length, in samples, of the subset of the full waveform identified as the bottom return by the ALPS software, y[n] is the digitized bottom return waveform; n = 0, 1, 2, …, N − 1 is the waveform sample number, A is the area under the curve, n is the bottom return waveform mean, σ is the standard deviation (width metric for the bottom return), and γ is the skewness (asymmetry metric for the bottom return).
Using the underwater image mosaics generated from the in situ (diver) data, we performed linear regressions to assess the associates between the percent cover of each of these groups (along with a combined hard coral group) and the lidar waveform metrics. Three of the nine sites were below a depth of 10 meters and were excluded from this analysis, for two key reasons. First, the manner by which the waveform metrics were interpolated to into GIS-compatible formats meant that the increasing footprint size of the lidar waveform with increasing depth was not represented, so the waveform metrics include higher amounts of error at greater depths. Second, the deeper sites were dominated by soft corals and experienced greater amounts of surge-the constant motion of the soft corals within the video meant that the percent cover metrics calculated from those photo mosaics were less accurate, as the same moving soft corals may appear in multiple locations or be excluded entirely. In situ reflectance (%) RelaƟve Reflectance (DN) In situ reflectance vs. relaƟve reflectance Fig. 10 Investigation of the strengths of the linear relationship between lidar-derived intensity or relative reflectance and the in situ reflectance data. The plots show the coefficient of determination from linear regression of (top) in situ reflectance on raw intensity and (bottom) in situ reflectance on relative reflectance

Results
The seafloor relative reflectance mosaic produced for the1 20-km 2 St. Thomas site is shown in Fig. 7. Data from 218 flight lines collected on eight separate project days were used in generating the St. Thomas mosaic. A similar mosaic was created for the~25-km 2 Buck Island site, encompassing data created from 30 EAARL-B flight lines on three separate days. Visual comparisons of the before-and-after data (i.e., raw intensity vs. relative reflectance), such as shown in Fig. 8, confirmed that the procedure greatly reduced or removed salient artifacts (e.g., seamlines, falloff at swath edges, rapid falloff with depth) and improved the ability to detect real seafloor features.
Depth and incidence angle corrections applied to the raw lidar intensity data had the desired effects reducing the impacts of these variables. An example is shown in Fig. 9, in which it can be seen that the corrections had the effect of flattening the curves (raw intensity or relative reflectance vs. along-track distance, color coded by flight line), indicating the expected trend of greater agreement between adjacent swaths and between different acquisition dates after applying the corrections. All such plots examined were consistent in showing substantially greater agreement after applying the depth and incidence angle corrections and normalizations. The strength of the linear relationship between lidar relative reflectance and in situ reflectance (at the laser wavelength of 532 nm) was also higher after applying the corrections, with the R 2 value improving from 0.46 to 0.73 (Fig. 10).
A pairwise heatmap plot, shown in Fig. 11, depicts the correlation between each of the variables for both the full depth range (bottom left portion of matrix) and just the 0-10 m subset (top right portion of matrix). A high correlation between the additional waveform metrics and relative reflectance or depth would indicate that the additional metrics did not add much useful information. For example, area under the curve is highly correlated (R 2 = 0.94) with relative reflectance, indicating that it may be of limited value if included, along with relative reflectance, in a benthic habitat classification procedure. Skewness and standard deviation, however, are not strongly correlated with depth or relative reflectance, suggesting that these waveform metrics provide additional, uncorrelated data, useful for benthic habitat classification.
Finally, the reef cover ground truth data for Flay Cay, obtained from the diver data, were analyzed along with the lidar waveform features to identify correlations between the lidar waveform features and cover type (Fig. 12). For sites shallower than 10 m depth, the standard deviation of Fig. 11 Pairwise heatmaps depict the correlation between water depth (D), relative reflectance (ρ rel ), area under the curve (A), skewness (γ), and standard deviation (σ). The main diagonal contains a histogram of the distribution of each variable for the full data set. The off-diagonal elements represent a heatmap comparing the two variables depicted by the labels for that row and column. The upper right off-diagonal elements contain the regressions for just the 0-10 m depth range, which is the primary depth range considered in this work, while the bottom left plots are for all depths. The solid gray lines are regression lines waveform skewness positively correlates with octocoral cover and negatively correlates with hard coral cover, explaining 62.6% and 36.5% of the variation in cover, respectively. Within hard coral cover, the contributions of domed and branched corals can be separated by other metrics. Mean waveform dispersion positively correlates with branched coral cover and explains 39.4% of variation, while standard deviation of area under curve positively correlates with domed coral cover and explains 41.3% of variation. No metric could meaningfully explain macroalgae cover or sponge cover. Dictyota turf was only highly abundant at one site, and negligible at most others. Sponges were not abundant at any site, but more so at sites > 10 m that were excluded from the analysis.

Discussion
Visual analysis and in situ seafloor reflectance comparisons showed that the lidar-derived relative reflectance data are substantially free of artifacts, such as seamlines, falloff at swath edges, and falloff with depth, which are salient in the raw intensity data. Profiles from adjacent swaths and different acquisition dates exhibit much greater agreement after our correction and normalization procedures and show substantially higher agreement with ground truth seafloor reflectance, with R 2 values improving from 0.46 to 0.73. Importantly, the procedures were designed to be efficient, such that they could be applied over large spatial extents, and, in the future, potentially even larger areas. The achieved processing performance for generating the seafloor relative reflectance mosaics was approximately one day of processing time per day of lidar data acquisition and could be further reduced by optimizing the data read/write functions within the software.
One result from our analysis of the relative reflectance data was deemed unusual and merits additional discussion. It was found that sun glint (specular solar reflection from the water surface) can impact the lidar bottom return intensities and also the relative reflectance values, if it is not accounted for in the processing. This effect was discovered in investigating intensities for one collection date, which exhibited across-flight line artifacts not consistent with the other two different collection dates over approximately the same region, even with parallel flight paths. Examination of aerial photographs taken at the time of collection for the anomalous flight line revealed substantial sun glint. Because this sun glint artifact is dependent on collection geometry (i.e., the elevation angle and azimuth of the sun at the time of acquisition), the fitted curve used to correct for angle of incidence artifacts had a less pronounced trend to model (Fig. 13). This led to some across-flight line artifacts persisting in the final results for this day. Previous research in bathymetric lidar has indicated that specular reflection from the water surface can change the gain in the photomultiplier tube(s) in other lidar systems, including the NASA Airborne Oceanographic Lidar (AOL) (Hoge et al. 1986;Hoge et al. 1988). We suspect that this is the cause of the change in the EAARL-B relative reflectance when sun glint is apparent in the concurrently collected Fig. 12 Relationships of percent cover of functional groups and highly correlated select lidar waveform metrics at the Flat Cay site RGB imagery, although this merits further investigation. Interestingly, this same effect was observed independently by Quantum Spatial Inc., using our procedures adapted for an entirely different topo-bathymetric lidar system: the Riegl VQ-880-G (E. Silvia, personal correspondence).

Conclusions
This study developed and tested new algorithms and techniques for producing relative reflectance mosaics and waveform metrics from the novel EAARL-B topobathymetric lidar. These new data sets where then investigated as an aid in distinguishing among hard corals and their morphologies. Results of the analysis performed using the waveform metrics suggest that these metrics may help to detect changes in the morphological composition of coral reef communities. Morphological types of branched and domed corals were positively correlated with standard deviation of the area under the curve and mean dispersion, respectively. This may be highly beneficial for groups who are engaged in creating benthic habitat maps for this region. Previous studies have shown that radiometric and geometric artifacts in imagery (whether optical or acoustic) can degrade the quality and accuracy of the resulting benthic habitat maps (Mumby et al. 1998, Costa and Battista 2013, Kumpumäki et al. 2015, Lecours et al. 2017). Thus, having artifact free, normalized seafloor reflectance will be critical for developing a quality habitat map for use by marine managers of the Buck Island Reef National Monument, Virgin Islands National Park, St. Thomas East End Reserve, and the Cas Cay-Mangrove Lagoon Marine Reserve and Wildlife Sanctuary. Furthermore, these methods could be extended to a wide range of coastal areas with clear, shallow (≤ 40m) water, to support a variety of management and regulatory agencies (state and federal), as well as other coastal stakeholder groups. Regarding the extensibility to other areas, it is important to note that the methods are in no way specific to coral reefs but can be applied to oyster reefs, seagrass beds, and any number of substrates and bedforms of ecological or management interest.
This study has also led to the identification of avenues for ongoing and future work. First, while the current procedures are reasonably efficient, it is anticipated that further efficiency gains-particularly, reduction of human time in the processcan be obtained by automating the determination of the coefficients in the depth and incidence angle corrections. The correction parameters (a, b) from the depth correction are functions of water clarity, while the parameters (α and β) from the incidence angle correction are functions of the seafloor reflectance properties. Through analysis of a large number of sites, Fig. 13 Example of artifacts in the relative reflectance mosaics caused by sun glint. The top image (a) shows sun glint visible in the imagery for one particular flight line. The bottom image (b) shows the corresponding cross section containing artifacts (discontinuities, or large differences in the curves) in relative reflectance products collected at a time when sun glint is present (green), but not in overlapping data sets (blue and red) collected with the sun at a lower elevation angle it may be possible to both obtain reasonable starting estimates of the parameters for new sites and also to predict when, where, and how frequently new parameters need to be determined. Second, the encouraging results from analyzing the lidar waveform metrics suggest that, combined with the relative reflectance mosaics, these waveform features may further benefit ecological assessments performed using lidar. Prior studies suggest that reefs are degrading and shifting from hard to soft coral or macroalgae (Fabricius et al. 2011;Inoue et al. 2013). Methods described in this study provide several metrics from which it may be possible to observe this degradation over 100s to 1000s of meters. Waveform metrics of standard deviation of the skewness and standard deviation of the area under the curve can provide qualitative assessment of changes in reef communities that can help to prioritize in situ, detailed monitoring of selected sites. Lastly, the interesting effect of specular solar reflections on lidar intensity and derived reflectance values merits further investigation. By incorporating an image processing step which identifies areas of sun glint into the procedure, it will likely be possible to better anticipate and correct for this effect.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.