Seafloor Elevation and Volume Change Analyses from 2016 to 2019 Along the Florida Reef Tract, USA

Step 1: Seafloor elevation and volume change analyses were performed using methods from Yates and others (2017). Original lidar data were downloaded from NOAA’s DigitalCoast website using the site’s 'Customized Download' capability, at https://coast.noaa.gov/dataviewer/#/lidar/search/. The full extent of the 2016 data in the form of four blocks were downloaded as well as the full extent of the 2019 data which is available as a single file. The data were downloaded using the following parameters: 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. Using VDatum version 3.9, a publicly available software from NOAA (https://vdatum.noaa.gov/), the lidar tagged image file format (TIFF) files for both 2016 and 2019 datasets were transformed from the North American Horizontal Datum of 1983 (NAD83) Universal Transverse Mercator (UTM) Zone 17 North to NAD83 National Spatial Reference System (NSRS2011) horizontal datum, and the vertical datum (North American Vertical Datum of 1988 [NAVD88]) and geoid model (GEOID12B) were kept the same.

Step 2: The 2016 lidar data (blocks 01-04) were merged into a single TIFF using the “Mosaic To New Raster (Data Management)” tool in ArcMap 10.7. Blocks were merged by prioritizing the most recently collected block in areas of overlap. This was accomplished by utilizing the following parameters: blocks 01-04 were used as input rasters and were loaded in numerical order; Output Location was set to a personal workspace; Raster Dataset Name with Extension: 2016_Data_Mosaic_Last.tif; Pixel Type (optional): 32_BIT_FLOAT; Number of Bands: 1; and Mosaic Operator (optional); LAST. All other parameters were left as default.

Step 3: Footprints of both the 2016 and 2019 lidar data were then created by opening the Digital Elevation Models (DEMs) in ArcMap. The “Reclassify (Spatial Analyst)” tool was used on both DEMs by replacing all old values with ‘1’ and keeping ‘No Data’ values as ‘No Data.’ The resulting TIFF was converted into a footprint using the “Raster to Polygon (Conversion)” tool to convert the raster files to polygon shapefiles. An intersect footprint was created between the 2016 and 2019 lidar data using the “Intersect (Analysis)” tool. The output polygon (2016_2019_lidar_footprint.shp) represents the shared coverage of 2016 and 2019 DEMs. Both DEMs were then clipped to the extent of the intersect footprint using the “Clip (Data Management) tool.

Step 4: The clipped DEMs were then edited for removal of reduced quality data. Using Global Mapper 20.1, each DEM was visually inspected using the methods of Yates and others (2017) to identify areas characterized by coarse interpolation relative to surrounding areas and appearing as small groups of large triangles created during the development of each DEM. Polygons were manually drawn to encompass the areas for removal using the “Digitizer” tool and the ‘Create Area/Polygon Features’ function. By default, these polygons are generated in the ‘User Created Features’ layer within the Control Center. The ‘User Created Features’ layer was then exported as a polygon shapefile. The 2016 data has four polygon shapefiles (Block01_2016lidar_RemovedPolygons.shp, Block02_2016lidar_RemovedPolygons.shp, Block03_2016lidar_RemovedPolygons.shp, and Block04_2016lidar_RemovedPolygons.shp) and the 2019 data has one polygon shapefile (RemovedPolygons_2019lidar.shp). The “Erase (Analysis)” tool in ArcMap was used to remove the areas of reduced quality from the 2016_2019_Intersect_Footprint.shp shapefile creating the 2016_2019_edited_lidar_footprint_intersect.shp shapefile. Using the “Clip (Data Management)” tool, areas of reduced quality data were removed from the 2016 and 2019 lidar DEMs by clipping each DEM to the 2016_2019_edited_lidar_footprint_intersect.shp shapefile.

Step 5: Data limitations exist in ArcMap 10.7 for elevation change analyses performed at this geographic extent (around 1 billion square meters). Approximately 250 million data points were required to sample the 2016 and 2019 DEMs using a 2-m grid spacing. However, point shapefiles were limited to approximately 80 million data points in ArcMap 10.7. DEMs were divided into five separate blocks that were used to generate five sets of 2-m grid shapefiles containing approximately 50 million points each to reduce the time required for elevation change computations. A polygon shapefile was created using the ‘Create Features’ function in ArcMap to cover roughly one fifth of the DEM beginning at the northeastern extent and was named Block01_footprint.shp. The DEMs were then clipped to the extent of the Block01_footprint.shp shapefile using the “Clip (Data Management)” tool. The symbology panel was opened for the resulting TIFF and classification statistics were inspected to determine if the resulting coverage was between 170 and 230 million cells, the Block01_footprint.shp polygon shapefile redrawn and the classification statistics were reinspected until the target range was met. A buffer was added to the block using a clipped portion of the intersect footprint shapefile that covered the extent of the block boundary with a few hundred meters of coverage on either side of the boundary. The clipped portion was added to the Block01_footprint.shp shapefile using the “Merge (Data Management)” tool. For Block 01, this buffer was added to the southern boundary, all other blocks had their buffer added to the northern boundary. The resulting polygon was clipped to the 2016_2019_edited_lidar_footprint_intersect.shp shapefile creating the Block01_TINextensionArea.shp shapefile. The buffer was required to prevent block edge artifacts in final volume change analysis steps requiring merging of elevation change TIN models for adjacent data blocks. Moving southwest along the DEM, another arbitrary polygon was drawn beginning at the boundary of the Block01_footprint.shp polygon covering the next fifth of the DEM as well as a slight overlap with the Block01_footprint.shp shapefile. The overlap extent of the Block01_footprint.shp shapefile was removed from the newly created polygon using the “Erase (Analysis)” tool in ArcMap creating the Block02.shp shapefile. The DEM was then clipped to the extent of the Block02.shp shapefile and classification statistics were inspected for the resulting TIFF to confirm the target range cell count. Following cell count target verification, the Block02.shp shapefile was clipped to the extent of the intersect footprint shapefile using the “Clip (Analysis)” tool resulting in the Block02_footprint.shp shapefile. A second version with a buffer was created using the same method used with the Block01 shapefiles, with the buffer being added to the northern boundary, creating the Block02_TINextensionArea.shp shapefile. This process was repeated until full coverage of the DEM was achieved, and 10 shapefiles were created: five polygons without a buffer, and five polygons with a buffer. This method resulted in a polygon overlap at the intersection of the polygons with buffers and full coverage with no overlap with the non-buffered polygon version.

Step 6: 2-m grids were generated for every block’s ‘TINextensionArea.shp’ shapefile in Global Mapper Version 20.1. Each 2-m grid was generated for each block shapefile by selecting the polygon with the “Digitizer” tool, and right clicking on the polygon within the frame and selecting “GRID – Create Regular Grid of User-Specified Size/Orientation…” from the “Advanced Feature Creation Options.” In the Grid Setup window, the following parameters were used: Grid Cell Width and Grid Cell Height were set to 2; the Calculate Grid Cell Counts to Fill Rectangle Option was selected; Crop to Selected Feature(s) was selected within the Grid Bounds window; Crop Generated Grid to Selected Area Feature was enabled; Keep Area if Any Part in Crop Area was enabled; Create Points at Grid Cell Centers was enabled; and all other options were left as default. The generated 2-m grid was then manually shifted to center grid points on the DEM pixels. The 2-m grid was selected using the “Digitizer” tool in Global Mapper, and in the Move/Reshape Feature(s) window, the Shift (Offset) Selected Feature(s) tab was selected to open the Specify Offset to Apply to Point(s) menu and the following parameters were set: Units: meters; X/Longitude Offset: +0.5; and Y/Latitude Offset: +0.5. Grid points that were shifted into no-data areas were removed to avoid extraction of no-data values. The relevant ‘TINextensionArea.shp’ polygon shapefile was selected using the Digitizer tool and the 2-m grid “Options” menu was opened by right clicking the 2-m grid layer in the “Control Center,” and “Export Layers to New File(s)” was selected from the “Layer” window. Within the “Select Export Format” menu, “Data Format” was set to Shapefile. The “Export Bounds” tab was selected, and the “Crop to Selected Area Features” parameter was enabled in the “Shapefile Export Options” window to remove grid points that were in no-data areas during shapefile export. The resulting shapefile was reimported into Global Mapper for the addition of horizontal coordinates. Horizontal coordinates were added to each grid point (based on the projection of the data frame) by selecting the shapefile with the “Digitizer” tool and “Add Coordinates/Bounds Attributes to Selected Features” was selected from the “Attribute/Style Functions” window.

Step 7: Elevation values were extracted from the 2016 and 2019 DEMs at the location of the 2-m grid points generated within each block’s ‘TINextensionArea.shp’ shapefile. The 2016 DEM was opened in Global Mapper and the resampling method was changed by opening the “Options” menu of the DEM and “Resampling” was set to No Resampling (Nearest Neighbor). The block’s 2-m grid was selected using the “Digitizer” tool. Within the “Attributes/Style Functions” window, “Apply Elevations to Selected Feature(s)” was selected, with “Terrain Layers” as the only selected parameter. An attribute named ELEVATION was created containing the elevation values extracted from the 2016 DEM. Using the “Digitizer” tool, the “Edit Selected Features” window was selected to rename the ELEVATION attribute to ELEV2016. The same steps were repeated to extract elevation values from the 2019 DEM at the location of grid points in each block’s shapefile and to rename the ELEVATION attribute containing the 2019 extracted elevation values to ELEV2019. Extracted elevation values are based on the source elevation so these data are in NAD83 (NSRS2011) UTM17N and meters.

Step 8: The elevation-difference between 2016 and 2019 extracted elevation values in each block was calculated using the “Digitizer” tool. The corresponding 2-m grid was selected for each block and “Calculate/Copy Attributes for Feature Selection” was selected from the “Attribute/Style Functions” window. The elevation-difference was calculated using the following parameters: “Select Existing or Create New Attribute to Assign Calculated Values”: Diff_m; “Source Attribute”: ELEV2019; “Operation”: Subtract; and “Use Attribute Value”: ELEV2016. ELEV2019 represents the modern elevation values and ELEV2016 represents the historical elevation values. The final 2-m grid for each block was exported in two separate formats including one shapefile using the block’s buffered footprint as an export bound (for example, Block01_TINextensionArea.shp polygon), and a second shapefile using the block’s original footprint as an export bound (for example, Block01_footprint.shp polygon). The “Export Layers to New File(s)” function with “Data Format” set to Shapefile and “Crop to Selected Area Feature(s)” was selected in the Export Bounds tab. Each block’s shapefile exported using the buffer export bound was named with respect to the Block it was generated in, for example, Block01_ElevationChangePoints_TINExtension.shp. Each block’s shapefile exported using the original block footprint polygon was named in the same format, for example, Block01_ElevationChangePoints.shp. Steps 6 through 8 were repeated for each block, with one notable exception. When shifting the points for Blocks 02 through 05, the shift values were determined by matching any point generated within the buffer zone to its nearest corresponding point in the adjacent block’s buffer zone in the buffered version of the grid shapefiles. Point generation in these blocks were based on extent bounds and therefore did not use a simple half-meter shift. Latitude and longitude values were noted for two points. For example, one point was chosen within the buffer zone of the Block01_TINextensionArea.shp shapefile and the coordinates of that point were documented, the nearest point that results from the Step 6 process step for Block02 also had its coordinates documented. The difference in X and Y coordinates were used to determine the horizontal and vertical shift to be used for that dataset.

Step 9: Elevation change surface models were created in ArcGIS Pro version 2.4.0 using the calculated elevation-difference (Diff_m) points from the various ‘TINExtension.shp’ point shapefiles. The five shapefiles used for elevation surface model generation are as follows: Block01_ElevationChangePoints_TINExtension.shp, Block02_ElevationChangePoints_TINExtension.shp, Block03_ElevationChangePoints_TINExtension.shp, Block04_ElevationChangePoitns_TINExtension.shp, and Block05_ElevationChangePoints_TINExtension.shp. Due to data size limitations in ArcMap, Triangulated Irregular Network (TIN) generation is only possible using ArcGIS Pro 2.1.3 or lRWE at this time. Any attempt to generate a TIN in ArcMap version 10.7 or earlier will result in a CPU memory error. Beginning with Block 01, the Block01_ElevationChangePoints.shp shapefile was loaded into ArcGIS Pro version 2.4.0. A TIN was created from the Diff_m points using the "Create TIN (3D Analyst)" tool by specifying the Block01_ElevationChangePoints.shp shapefile as the "Input Feature Class", Diff_m as the "Height Field" and Mass_Points as the "Type", creating the Block01_Extension_TIN file. The Block01_Extension_TIN was then delineated using the “Delineate TIN Data Area (3D Analyst)” tool by specifying the Block01_Extension_TIN as the “Input TIN”, a "Maximum Edge Length" of 2.828428 (hypotenuse of a triangle with 2-m legs) and the "Method" set to ALL EDGES. It was noted that the resultant delineated TIN needed to be reloaded in a new application for changes to be applied. The tool was allowed a few minutes post completion and a new ArcGIS Pro application was used to import the delineated TIN. The delineated Block01_Extension_TIN file was then clipped to the extent of the 2016 and 2019 DEMs using the “Edit TIN (3D Analyst)” tool with the following parameters: "Input TIN": Block01_Extension_TIN file; "Input Features Class": Block01_TINextensionArea.shp; "Height Field": None; "Tag Field": None; and "Type": Hard clip. As with the delineation step, the file needed to be reloaded in a new application after tool completion for changes to be applied. The Block01_Extension_TIN file was then clipped again using clipped versions of the habitat shapefile (outlined in Step 10) under the following parameters: "Input TIN": Block01_Extension_TIN file; "Input Features Class": Block01_HabMapClip.shp; "Height Field": None; "Tag Field": None; and "Type": Hard clip. An application reset was used after this tool as well. The final TIN was exported as the Block01_ProcessedTIN file. This process was repeated for each ‘TINExtension.shp’ point shapefile, creating five TINs covering the extent of the study area. Please note the very large files and ensure you have the adequate computing compacity to access these files. These TIN files were generated and can be viewed using ArcGIS Pro version 2.1.3 or later to accommodate the large file sizes. They may not be accessible in other geoprocessing software with file size limits. Smaller TIN files can be generated for viewing in other programs using the following methods: generate a polygon shapefile covering a data extent that does not exceed available computing capabilities. Generate a point shapefile that covers the extent of polygon shapefile using the “Clip (Analysis)” tool and relevant elevation-change point shapefile (ElevationChangeShapefiles.zip). Follow process steps 7 through 9 for developing a TIN at this smaller data extent

Step 10: The original Unified Florida Reef Tract Map version 2.0 shapefile was downloaded from https://ocean.floridamarine.org/IntegratedReefMap/UnifiedReefTract.htm. Using ArcMap, the original habitat shapefile was modified using the "Clip (Analysis)" tool to clip the habitat shapefile to the extent of the 2016 and 2019 DEMs by specifying the habitat shapefile as the "Input Features" and the 2016_2019_edited_lidar_footprint_intersect.shp as the "Clip Features", creating the UR_HabMap_Clip.shp shapefile. The shapefile was then loaded into Global Mapper and exported using the 2011 realization of NAD83 as UR_Habmap_2011_GM_Export.shp. The DEMs used in this analysis cover areas that the Unified Reef map do not, so an additional polygon was integrated into the shapefile. Using the "Erase (Analysis)" tool in ArcMap, the 2016_2019_edited_lidar_footprint_intersect.shp shapefile was used as the "Input Features" and the UR_HabMap_2011_GM_Export.shp shapefile was used as the "Erase Features" creating the Unclassified_Polygons.shp shapefile. The Unclassified_Polygons.shp shapefile was edited in an edit session to add the "ClassLv2" attribute and assign its value as "Unclassified.” The edit session was closed, and the changes were saved. Using the "Merge (Data Management)" tool, the Unclassifed_Polygons.shp shapefile and the UR_HabMap_2011_GM_Export.shp shapefile were used as the "Input Datasets" creating the UR_HabMap_2011_GM_Export_UnclassifiedAdded.shp shapefile. Using the "Select by Attribute" tool, 17 individual habitat shapefiles were created from the UR_HabMap_2011_GM_Export_UnclassifiedAdded.shp shapefile by selecting one ClassLv2 habitat and exporting it as a separate shapefile. For volume analyses, this habitat file was split into five sections by using the “Clip (Analysis)” tool in ArcMap and clipping the UR_Habmap_2011_GM_Export_UnclassifiedAdded.shp file to each of the individual block footprints (Block01_footprint.shp, Block02_footprint.shp, Block03_footprint.shp, Block04_footprint.shp, and Block05_footprint.shp) creating the following files: Block01_HabMapClip.shp, Block02_HabMapClip.shp, Block03_HabMapClip.shp, Block04_HabMapClip.shp, and Block05_HabMapClip.shp.

Step 11: Five sets of elevation change statistics were determined by habitat type for each block using the elevation-difference points from the following five shapefiles: Block01_ElevationChangePoints.shp, Block02_ElevationChangePoints.shp, Block03_ElevationChangePoints.shp, Block04_ElevationChangePoints.shp, and Block05_ChangeElevationPoints.shp. In ArcGIS Pro, 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": Block01_ElevationChangePoints.shp; "Relationship": INTERSECT; "Selecting Features": Clipped Habitat shapefile; "Search Distance": left blank; and "Selection Type": NEW_SELECTION. An ArcGIS Pro model named Seafloor Elevation Change Analysis Tool (SECAT) was created to automate the process as these steps had to be repeated for 17 habitat types across each block's shapefile. Elevation change statistics were compiled by habitat type into a comma separated values (CSV) file using Microsoft Excel 2016, see Block01_Elevation.csv, Block02_Elevation.csv, Block03_Elevation.csv, Block04_Elevation.csv, Block05_Elevation, and 2016_2019_FRT_CollectiveElevationStatistics.csv included in SECAT_Analyses.zip. Excel sheet versions (XSLX) of the collective statistics are also included for convenience.

Step 12: Using ArcGIS Pro, volume change statistics per habitat type were calculated utilizing each of the TINs generated in Step 9. Surface volume changes were calculated for six 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" was set to -0.19 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.19 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.19 m threshold was determined by vertical error analysis using the uncertainties reported in the metadata of the original 2016 (0.15 m) and 2019 (0.11 m) DEMs to calculate the Root Mean Square Error (RMSE) of 0.19 m (Yates and others, 2017). An additional threshold of 0.31 m (1.65 x RMSE) was used as a conservative threshold. The "Reference Plane" was set to BELOW and the "Plane Height" was set to -0.31 m for the net erosion lower limit at this threshold (case 5). The "Reference Plane" was set to ABOVE and the "Plane Height" was set to 0.31 m for the net accretion lower limit at this threshold (case 6). Minimum net volume change was calculated by summing results from cases 1 and 3 or 1 and 5 depending on lower limit threshold. Maximum net volume change was calculated by summing results from cases 2 and 4 or 2 and 6 depending on lower limit threshold. The 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 dividing the maximum net volume for each habitat by the habitat's total area. An ArcGIS Pro model was created to automate the process, as these steps had to be repeated for 17 habitat types. Volume change statistics were generated for each individual block and the results were then compiled by habitat type in CSV format using Excel, see 2016_2019_FRT_CollectiveVolumeStatistics_019.csv and 2016_2019_FRT_CollectiveVolumeStatistics_031.csv as well as the individual block statistics at each threshold which are also included in SECAT_Analyses.zip. Excel sheet versions (XLSX) of the collective statistics are also included for convenience.