Seafloor elevation change from 2002 to 2016 in the Upper Florida Keys

Step 1: The original 2002 Experimental Advanced Airborne Research Lidar (EAARL) Submarine Topography: Biscayne National Park and EAARL Submarine Topography: Northern Florida Keys Reef Tract digital elevation models (DEMs) were downloaded from two U.S. Geological Survey (USGS) Open-File Reports (Brock and others, 2006 and Brock and others, 2007). Horizontal coordinates were provided in the North American Horizontal Datum of 1983 (NAD83) Universal Transverse Mercator (UTM) Zone 17 North, GRS80 ellipsoid; and vertical positions were referenced to the North American Vertical Datum of 1988 (NAVD88), in meters (m). Global Mapper version 19.1 was used to export the two DEMs as a single, merged lidar DEM in the geographic tagged image file format (GEOTIFF). Using VDatum version 3.9, a publicly available software from the National Oceanic and Atmospheric Administration (NOAA) (https://vdatum.noaa.gov/, accessed on September 7, 2018), the merged 2002 lidar DEM was transformed from NAD83 to NAD83 National Spatial Reference System of 2007 (NSRS2007) National Readjustment horizontal datum, and from GEOID03 to the GEOID12B geoid model.

Step 2: The original 2016 DEMs were acquired from NOAA’s Data Access Viewer (https://coast.noaa.gov/dataviewer/#/lidar/search/) using the “Custom Download” capability of NOAA’s Digital Coast website. The elevation search option was used to download three topobathymetric datasets: 2016 NOAA NGS Topobathy Lidar: Florida Keys Outer Reef Blocks 02, 03, and 2017 NOAA NGS Topobathy Lidar: Florida Keys Outer Reef Block 04 (National Oceanic and Atmospheric Administration, 2017a, b; and National Oceanic and Atmospheric Administration, 2018). The data were downloaded with the following parameters: “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. Using Global Mapper, the three DEMs were merged into a single DEM and transformed into the NAD83 (NSRS2007) National Readjustment horizontal datum.

Step 3: The 2016 DEM contained a small region of sub-aerial elevations that were manually removed in Global Mapper. Polygons encompassing the sub-aerial regions were created using the “Digitizer” tool and the “Create New Feature Area” function. The polygons were selected, and the “Cut Out Currently Selected Polygon(s)” parameter was enabled within the “Cropping” tab of the 2016 DEM “Elevation Options” window to remove the section of the 2016 DEM that intersects with the polygons encompassing the sub-aerial elevations.

Step 4: Using Esri ArcGIS Desktop Advanced version 10.6 (ArcMap), footprints of the original 2002 and 2016 DEMs were created with the “Reclassify (Spatial Analyst)” tool. To create each raster file, all old data values were replaced with the number 1 and “No Data” values were left as “No Data”. The “Raster to Polygon (Conversion)” tool was then used to create a footprint of each original DEM by converting the raster files to polygon shapefiles, creating the 2002_UFK_Footprint and 2016_UFK_Footprint shapefiles.

Step 5: A polygon shapefile of the geometric intersection between the two DEMs was created with the “Intersect (Analysis)” tool by adding the 2002_UFK_Footprint and 2016_UFK_Footprint shapefiles as “Input features”, creating the Intersect_Footprint shapefile. Then, the 2002 and 2016 DEMs were clipped to the extent of the Intersect_Footprint shapefile using the “Clip (Data Management)” tool by specifying one of the DEMs as the “Input Raster”, the Intersect_Footprint shapefile as the “Output Extent”, and box checked for “Use Input Features for Clipping Geometry”, creating the 2002_UFK_LidarClip and 2016_UFK_LidarClip GEOTIFFs.

Step 6: A 2-m grid was generated in Global Mapper using the “Digitizer” tool. To define the correct geographic location and extent of the 2-m grid, the Intersect_Footprint shapefile was opened. When the Intersect_Footprint shapefile was opened, it automatically displayed in the NAD83 horizontal datum, despite being set to the NAD83 (NSRS2007) National Readjustment horizontal datum. Using the “Configure” tool, the correct horizontal datum was specified. The “Digitizer” tool was used to select the Intersect_Footprint shapefile, and within the “Advanced Feature Creation Options” window, “Create Regular Grid of User Specified Size/Orientation” was selected to create a 2-m spaced grid with the following parameters: “Grid Cell Width”: 2; “Grid Cell Height”: 2; “Calculate Grid Cell Counts to Fill Rectangle”: enabled; “Crop Generated Grid to Selected Area Feature”: enabled; “Keep Area if Any Part in Crop Area”: enabled; “Create Points at Grid Cell Centers”: enabled; and “Crop to Selected Area Feature(s)”: enabled. The remaining default parameters were left the same, and the tool was run, creating the 2-m grid at the extent of the Intersect_Footprint shapefile.

Step 7: The 2-m grid was shifted to center the 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.

Step 8: Grid points that were shifted into no-data areas were removed to avoid the extraction of no-data values. The Intersect_Footprint shapefile was selected using the “Digitizer” tool, and the 2-m grid was left unselected. 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. By exporting the 2-m grid with this parameter enabled, grid points that were in no-data areas were removed, creating the 2m_grid shapefile. Before the next steps, the 2m_grid shapefile was reimported into Global Mapper.

Step 9: Horizontal coordinates (x,y) were added by selecting the 2m_grid shapefile with the “Digitizer” tool and “Add Coordinates/Bounds Attributes to Selected Features” was selected from the “Attribute/Style Functions” window. Horizontal coordinates were added to each grid point based on the projection of the data frame.

Step 10: Elevation values were extracted from the 2002 and 2016 DEMs at the location of the 2m_grid shapefile. The 2002 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 DEM and 2m_grid shapefile were 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 2002 DEM. Using the “Digitizer” tool, the “Edit Selected Features” window was selected to rename the ELEVATION attribute to ELEV2002. The same steps were repeated to extract elevation values from the 2016 DEM at the location of the 2m_grid shapefile and to rename the ELEVATION attribute containing the 2016 extracted elevation values to ELEV2016.

Step 11: The elevation-difference between the 2002 and 2016 extracted elevation values was calculated using the “Digitizer” tool. The 2m_grid shapefile was selected 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”: ELEV2016; “Operation”: Subtract; and “Use Attribute Value”: ELEV2002. ELEV2016 represents the modern elevation values and ELEV2002 represents the historical elevation values. The final 2-m grid was exported using the “Export Layers to New File(s)” function with “Data Format” set to Shapefile, creating the 2002_2016_UFK_ElevationChange_Points shapefile.

Step 12: The original Unified Florida Reef Tract Map version 2.0 shapefile was downloaded from http://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 2002 and 2016 DEMs by specifying the habitat shapefile as the “Input Features” and the Intersect_Footprint shapefile as the “Clip Features”, creating the 2002_2016_UFK_HabitatClip shapefile. Using the “Select by Attribute” tool, 13 individual habitat shapefiles were created from 2002_2016_UFK_HabitatClip shapefile by selecting one ClassLv2 habitat and exporting as a separate shapefile.

Steps 13: Elevation change statistics were determined by habitat type using the elevation-difference points from the 2002_2016_UFK_ElevationChange_Points shapefile in ArcGIS Pro version 2.1.3. 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”: 2002_2016_UFK_ElevationChange_Points; “Relationship”: INTERSECT; “Selecting Features”: Habitat shapefile; “Search Distance”: left blank; and “Selection Type”: NEW_SELECTION. An ArcGIS Pro model was created to automate the process, since these steps had to be repeated for 13 habitat types. Elevation change statistics were compiled by habitat type into a comma separated values (CSV) file using Microsoft Excel 2016, see 2002_2016_UFK_ElevationStatistics.csv.

Step 14: An elevation change surface model was created in ArcMap using the calculated elevation-difference (Diff_m) points from the 2002_2016_UFK_ElevationChange_Points shapefile. Due to the memory limitation ArcMap encounters when attempting to edit large Triangulated Irregular Networks (TINs), a Python script was developed to create and edit a TIN following the Yates and others (2017) TIN creation method. Using the Python script, a TIN was created from the Diff_m points using the “Create Tin (3D Analyst)” tool by specifying the 2002_2016_UFK_ElevationChange_Points shapefile as the “Input Feature Class”, Diff_m as the “Height Field” and Mass_Points as the “Type”, creating the 2002_2016_UFK_ElevationChange_TIN file. The 2002_2016_UFK_ElevationChange_TIN file was then delineated using the “Delineate TIN Data Area (3D Analyst)” tool by specifying the 2002_2016_UFK_ElevationChange_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. The delineated 2002_2016_UFK_ElevationChange_TIN file was clipped to the extent of the 2002 and 2016 DEMs using the “Edit TIN (3D Analyst)” tool with the following parameters: “Input TIN”: 2002_2016_UFK_ElevationChange_TIN file; “Input Features Class”: Intersect_Footprint shapefile; “Height Field”: None; “Tag Field”: None; and “Type”: Hard clip. Due to the large file size of the final 2002_2016_UFK_ElevationChange_TIN file, ArcMap couldn’t render or run geoprocessing tools on the TIN, and therefore, ArcGIS Pro was used for subsequent processing steps.

Step 15: Volume change statistics per habitat type were calculated using the final 2002_2016_UFK_ElevationChange_TIN file in ArcGIS Pro. 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.24 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.24 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.24 m threshold was determined by vertical error analysis using the uncertainties reported in the metadata of the original 2002 (0.20 m) and 2016 (0.15 m) DEMs to calculate the Root Mean Square Error (RMSE) of 0.24 m. Minimum net volume change was calculated by summing results from cases 1 and 3. Maximum net volume change 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 dividing the maximum net volume for each habitat by the habitat's total area. An ArcGIS Pro model was created to automate the process, since these steps had to be repeated for 13 habitat types. Volume change statistics were compiled by habitat type in CSV format using Excel, see 2002_2016_UFK_VolumeStatistics.csv.