Seafloor elevation change from the 1930s to 2016 along the Florida Reef Tract, USA

Step 1: Acquisition and preparation of digital elevation model (DEM) data and elevation- and volume-change analyses was performed using Yates and others (2017) methods.

Step 2: The original 2016 DEMs were acquired from NOAA’s Digital Coast website (https://coast.noaa.gov/digitalcoast/) using the Data Access Viewer (DAV) tool. The elevation search option was used to download four topobathymetric datasets: 2016 NOAA NGS Topobathy Lidar: Florida Keys Outer Reef Blocks 01, 02, 03, and 2017 NOAA NGS Topobathy Lidar: Florida Keys Outer Reef Block 04. The DAV was used to generate a custom DEM from the point cloud data. The four DEMs were created and downloaded in September 2017 and correspond to the North American Datum of 1983 (NAD83) National Spatial Reference System (NSRS2007) horizontal datum, and North American Vertical Datum of 1988 (NAVD88) GEOID12B vertical datum. The data were downloaded from the DAV with the following specifications; Projection: UTM, Zone: Zone 17 Range 084W-078W, Horizontal Datum: NAD83, Horizontal Units: Meters, Vertical Datum: NAVD88, Vertical Units: Meters, File Format: TIFF 32-bit Float, Bin Method: TIN, Bin Size: 1.0, Bin Units: Meters, Data Classification: Bathymetric Lidar Points, Data Returns: Any Points, Ancillary Data: No Ancillary Data and Geoid Name: GEOID12B.

Step 3: The original DEMs were visually inspected and suspect data characterized by coarse interpolation relative to surrounding areas were removed using the methods of Yates and others (2017). Polygons encompassing removed data for each block are provided in Yates and others (2021), see 2016_Lidar_RemovedAreas.zip. The DEMs were opened in ascending order starting with the Block01 DEM TIFF and ending with the Block04 DEM TIFF to prioritize the most recent DEM where there was overlapping data. Using Esri ArcGIS Desktop Advanced version 10.6 (ArcMap), the four edited DEMs were merged together. The merged DEM was created using the "Mosaic to New Raster" tool with the following parameters: Pixel Type: 32_BIT_FLOAT; Number of Bands: 1; and Mosaic Operator: LAST, creating the FRT_Merged_LidarDEM TIFF.

Step 4: A footprint of the merged DEM was created using the "Reclassify (Spatial Analyst)" tool to convert all old values to 1 and leaving 'No Data' values as 'No Data'. Then, the output raster was converted into a polygon SHP file using the "Raster to Polygon (Conversion)" tool, creating the FRT_Merged_LidarDEM_Footprint SHP file.

Step 5: The original digitized H-sheet XYZ files were downloaded from the NOAA Bathymetric Data Viewer (https://www.ngdc.noaa.gov/maps/bathymetry/#). Using VDatum version 3.9, a publicly available software from NOAA (https://vdatum.noaa.gov/), the H-sheet’s vertical datum was transformed from mean lower low water (MLLW) to NAVD88 (GEOID12B), and the horizontal datum was transformed from geographic NAD83 (1986) to UTM Zone 17N NAD83 (NSRS2007). The H-sheet data were also converted from soundings in meters (m) to height (m). For more information on downloading, prepping, and using the H-sheet data, see Yates and others (2017).

Step 6: The output XYZ files from VDatum were converted to comma-separated values (CSV) files and displayed as XY data. Each H-sheet CSV was opened in ArcMap and Add Data > Add XY Data was selected from the File menu. Then, in the Add XY Data window, an H-sheet XYZ file was selected, and LON, LAT, and DEPTH were used for the X, Y, and Z fields, respectively, creating a point SHP file. This process was repeated for each H-sheet, creating a total of 12 individual point SHP files. Then, the individual point SHP files were merged together using the "Merge (Data Management)" tool, creating the FRT_HistoricalPoints SHP file. For more information on prepping and using the H-sheet data, see Yates and others (2017).

Step 7: The FRT_HistoricalPoints SHP file depth values were corrected for sea level rise by applying relative sea level rise (RSLR) estimates collected by NOAA at the Key West and Vaca Key, Florida tide stations. The RSLR correction was calculated by adding a field to the attribute table of the FRT_HistoricalPoints SHP file using the Field Calculator and the expression RSLR_2016 = DEPTH - ((2016 - [year]) * 0.00302), where 0.00302 is the mean annual rate of sea level rise computed from the tide stations. For additional information on applying the sea level rise correction, see Yates and others (2017).

Step 8: The historical points SHP file was clipped to the extent of the merged DEM using the "Clip (Analysis)" tool by specifying the FRT_HistoricalPoints SHP file as the 'Input Features' and the FRT_Merged_LidarDEM_Footprint SHP file as the 'Clip Features', creating the FRT_HistoricalPoints_Clip SHP file.

Step 9: Values from the merged DEM were extracted at the location of the historical points using the "Extract Values to Points (Spatial Analyst)" tool by specifying the FRT_HistoricalPoints_Clip SHP file as the 'Input Point Features' and the FRT_Merged_LidarDEM TIFF as the 'Input Raster', creating the FRT_IntersectPoints SHP file. The elevation difference (Diff_m) between the H-sheets and modern lidar data were calculated by adding a field to the attribute table of the FRT_IntersectPoints SHP file using the Field Calculator and the expression Diff_m = RASTERVALU (lidar) – RSLR_2016 (H-sheets).

Step 10: The original Unified Florida Reef Tract Map version 2.0 SHP file was downloaded from http://ocean.floridamarine.org/IntegratedReefMap/UnifiedReefTract.htm. Using ArcMap, the original habitat SHP file was modified using the "Clip (Analysis)" tool to clip the habitat SHP file to the extent of the merged DEM by specifying the habitat SHP file as the 'Input Features' and the FRT_Merged_LidarDEM_Footprint SHP file as the 'Clip Features', creating the FRT_HabitatClip SHP file. Using the "Select by Attribute" tool, 14 individual habitat SHP files were created from the FRT_HabitatClip SHP file by selecting one ClassLv2 habitat and exporting it as a separate SHP file.

Step 11: Elevation change statistics were determined by habitat type using the XYZ points from the FRT_IntersectPoints SHP file. The "Select Layer by Location (Data Management)" tool was used to extract points within or on the boundary of a specific habitat type by using the following parameters: Input Feature Layer: FRT_IntersectPoints; Relationship: INTERSECT; Selecting Features: Habitat SHP file; Search Distance: left blank; and Selection type: NEW_SELECTION. An ArcMap model was created to automate the process, since these steps had to be repeated for 14 habitat types. Elevation change statistics were compiled by habitat type into a CSV file using Microsoft Excel 2016, see FRT_Elevation_Statistics.csv.

Step 12: An elevation change surface model was created using the "Create TIN (3D Analyst)" tool by specifying the FRT_IntersectPoints SHP file as the 'Input Feature Class', Diff_m as the 'Height Field' and Mass_Points as the 'Type', creating the Intersect_TIN file. Then, the Intersect_TIN file was delineated using the "Delineate TIN Data Areas (3D Analyst)" tool by specifying the Intersect_TIN as the 'Input TIN', a 'Maximum Edge Length' of 400 m (Yates and others, 2017) and the 'Method' set to ALL. A polygon footprint of the delineated TIN was created using the "TIN Domain (3D Analyst)" tool by specifying the Intersect_TIN file as the 'Input TIN' and the 'Output Feature Class Type' as POLYGON, creating the Intersect_TINDomain SHP file. A polygon representing the intersection of the Intersect_TIN file and merged DEM was generated using the "Clip (Analysis)" tool by specifying the Intersect_TINDomain SHP file as the 'Input Features' and the FRT_Merged_LidarDEM_Footprint SHP file as the 'Clip Features', creating the FRT_TINDomainIntersect SHP file. Then, the Intersect_TIN file was clipped to the extent of the intersect footprint SHP file using the "Edit TIN (3D Analyst)" tool with the following parameters: Input TIN: Intersect_TIN; Input Features Class: FRT_TINDomainIntersect SHP file; Height Field: None; and Type: Hard clip, creating the final TIN file.

Step 13: Volume-change statistics per habitat type were calculated using the final TIN file. Surface volume changes were calculated for four cases using the "Surface Volume (3D Analyst)" tool. To calculate the net erosion lower limit (case 1) the 'Reference Plane' was set to BELOW and the 'Plane Height' set to -0.5 m. For the net erosion upper limit (case 2) the 'Reference Plane' was set to BELOW and the 'Plane Height' set to 0 m. For the net accretion lower limit (case 3) the 'Reference Plane' was set to ABOVE and the 'Plane Height' was set to 0.5 m. For the net accretion upper limit (case 4) the 'Reference Plane' was set to ABOVE and the 'Plane Height' was set to 0 m. A 0.5 m threshold was determined by vertical error analysis using the uncertainties reported in the metadata of the original lidar (0.15 m) and the H-sheet (0.23 m) products to calculate the Root Mean Square Error (RMSE) of 0.29 m. The RMSE of 0.29 was multiplied by 1.56 to compute a RMSE of 0.47 m to include 90% of the variance. This value was rounded to 0.5 to produce the final RMSE used as a threshold in this study (Yates and others, 2017). Minimum net volume was calculated by summing results from cases 1 and 3. Maximum net volume was calculated by summing results from cases 2 and 4. Area normalized volume change lower limit was calculated by dividing the minimum net volume change for each habitat by the habitat's total area. The area normalized volume change upper limit was calculated by diving the maximum net volume for each habitat by the habitat's total area. An ArcMap model was created to automate the process, since these steps had to be repeated for 14 habitat types. Volume change statistics were compiled by habitat type in CSV format using Excel, see FRT_Volume_Statistics.csv.