Seafloor Elevation Change from 2017 to 2019 at Looe Key, Florida Keys—Impacts from Hurricane Irma

Step 1: Seafloor elevation and volume change analyses were performed using methods from Yates and others (2017). The original 2017–2018 multibeam data were downloaded from Fredericks and others (2019). Using VDatum v.3.9, the multibeam data were transformed from their native World Geodetic System of 1984 (WGS84) horizontal datum and ellipsoid heights to the North American Datum of 1983 (NAD83 [2011]) horizontal datum and North American Vertical Datum 1988 (NAVD88) vertical datum, applying the GEOID12B model. The transformed XYZ points were loaded into Global Mapper version 18.2 and gridded using the "Create Elevation Grid from 3D Vector/Lidar Data" tool in the Analysis menu. Grid spacing was manually set to 1 m for the X- and Y-axes and the 'Elevation Grid No Data Distance Criteria' was set to 3.0. The resultant 2017_LooeKey_Multibeam_Clip DEM was exported as a Tagged Image File Format (TIFF, .tif) with the following parameters: 32-bit floating point samples, Sample Spacing of 1 m for both X- and Y-axes, Always Generate Square Pixels, LZW Compression, Generate TFW (World) File, and Generate PRJ File. The 2019 lidar data were downloaded from National Geodetic Survey (2020). The data were downloaded with 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

Step 2: Using ArcGIS Pro version 3.4.3, footprints of both the 2017–2018 multibeam data and 2019 lidar data were made using the "Reclassify (Spatial Analyst)" tool to replace all valid values with 1 and setting the 'No Data' value to NoData to create a raster mask of data coverage for each source dataset. The "Raster to Polygon (Conversion)" tool was then used to convert those raster files to polygon shapefiles (SHP) representing the footprints of each dataset. The intersection of these footprints was generated using the "Intersect (Analysis)" tool and using the multibeam and lidar shapefile footprints as inputs to create the intersection footprint (Intersect_footprint SHP file). Both digital elevation models (DEMs) were then clipped to the extent of the intersect footprint using the "Clip (Data Management)" tool, creating the 2017_LKR_multibeam_19lidarExtent.tif and 2019_LKR_lidar_17multibeamExtent.tif files.

Step 3: A 2-m grid was then created using the "Create Fishnet (Data Management)" tool, using the following parameters: Template extent: Intersect_footprint SHP file (created in Step 2); Cell size width: 2; Cell size height: 2; Number of Rows: left blank, Number of Columns: left blank; Geometry type: POLYLINE and box checked for 'Create Label Points'. The 2-m grid label points SHP file was clipped to the extent of the Intersect_footprint SHP file using the "Clip (Analysis)" tool by specifying the 2-m grid label points SHP file as the 'Input features' and the Intersect_footprint SHP file as the 'Clip features.' XY coordinates were added to the 2-m grid SHP file using the "Add XY Coordinates (Data Management)" tool, to create the 2m_grid point SHP file with coordinate fields POINT_X and POINT_Y in Universal Transverse Mercator (UTM) units.

Step 4: Elevation values were extracted from the 2017 and 2019 DEMs at the location of the 2-m grid points in Global Mapper. The clipped 2017 DEM and the 2-m grid shapefile were loaded into Global Mapper, and the resampling method was changed to No Resampling (Nearest Neighbor) within the "Options" menu of the DEM. The 2-m grid was then selected using the "Digitizer" tool and 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 2017 DEM values. Using the "Digitizer" tool, the "Edit Selected Features" window was selected to rename the ELEVATION attribute to ELEV2017. The same steps were repeated for the clipped 2019 DEM, creating the ELEV2019 field.

Step 5: Elevation-differences between the ELEV2017 and ELEV2019 values were calculated using the "Digitizer" tool in Global Mapper. "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”: ELEV2017. ELEV2019 represents the modern elevation values and ELEV2017 represents the historical elevation values in this analysis. The final 2-m grid was then exported as a shapefile, 1719_LKR_IntersectPoints.shp.

Step 6: Elevation change surface models were created in ArcGIS Pro using the calculated elevation-difference (Diff_m) field within the 1719_LKR_IntersectPoints.shp point shapefile. A triangular irregular networks (TIN) file was created from the Diff_m points using the "Create TIN (3D Analyst)" tool by specifying the 1719_LKR_IntersectPoints.shp shapefile as the "Input Feature Class", Diff_m as the "Height Field" and Mass_Points as the "Type", creating the 1719_LK_Diffm_TIN file. The TIN was then delineated using the "Delineate TIN Data Area (3D Analyst)" tool by specifying the 1719_LK_Diffm_TIN file 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.

Step 7: The Unified Florida Reef Tract Map Version 2.0 shapefile was downloaded from Florida Fish and Wildlife Conservation Commission (2015). Using ArcGIS Pro, the habitat shapefile was modified using the "Clip (Analysis)" tool to clip the habitat shapefile to the extent of the 2017 and 2019 DEMs by specifying the habitat shapefile as the "Input Features" and the Intersect_footprint SHP file from step 2 as the "Clip Features", creating the LK_HabMap_Clip.shp file. The shapefile was then loaded into Global Mapper and exported using the 2011 realization of NAD83 as LK_Habmap_2011_GM_Export.shp.

Step 8: The DEMs used in this analysis cover additional areas outside of the Unified Florida Reef Tract Map Version 2.0, so an additional polygon needed to be integrated into the shapefile. Using the "Erase (Analysis)" tool in ArcGIS Pro, the Intersect_footprint SHP file was used as the "Input Features" and the LK_HabMap_2011_GM_Export.shp file was used as the "Erase Features" creating the LK_HabMap_2011_GM_Export_UnclassifiedAdded file. 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.

Step 9: Using the "Merge (Data Management)" tool, the Unclassifed_Polygons.shp file and the LK_HabMap_2011_GM_Export.shp shapefile were used as the "Input Datasets" creating the 1719_LKR_habitat_clipped.shp file. Using the "Select by Attribute" tool, 11 individual habitat shapefiles were created from the LK_HabMap_2011_GM_Export_UnclassifiedAdded.shp file by selecting one ClassLv2 habitat and exporting it as a separate shapefile.

Step 10: Elevation change statistics were determined by habitat type from the LooeKey_1719_IntersectPoints.shp file. In ArcGIS Pro, the "Select Layer by Locations (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": LooeKey_1719_IntersectPoints.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; Zieg and Zawada, 2021) was created to automate the process, as these steps were repeated for 11 habitat types across the extent of the study area. Elevation change statistics were compiled by habitat type into a single comma separated values (CSV) file, see 1719_LKR_Elevation_and_Volume_Statistics.zip for this file.

Step 11: Volume change statistics per habitat type were calculated utilizing the TIN generated in Step 6. Surface volume changes were calculated using the "Surface Volume (3D Analyst)" tool in ArcGIS Pro. To calculate the net erosion lower limit, the "Reference Plane" was set to BELOW and the "Plane Height" was set to -0.19 m. For the net erosion upper limit, the "Reference Plane" was set to Below, and the "Plane Height" to 0 m. Net accretion lower and upper limits were calculated by setting the "Reference Plane" to ABOVE and the "Plane Height" to 0.19 m and 0 m for each case, respectively. The 0.19 m threshold was determined via vertical error analysis of the reported vertical uncertainties for the original 2017 multibeam (0.15 m) and 2019 lidar (0.11 m) to calculate a root mean square error (RMSE) of 0.19 m (Yates and others, 2017). 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. The SECAT model (Zieg and Zawada, 2021) was used to automate the process, as these steps were repeated for 11 habitat types. The volume statistics were compiled by habitat type into a single CSV file, see 1719_LKR_Elevation_and_Volume_Statistics.zip for this file.