30 meter Esri binary grids of coastal response type probabilities with respect to projected sea levels for the Northeastern U.S. from Maine to Virginia for the 2020s, 2030s, 2050s and 2080s (Albers, NAD 83)

Interpolate surfaces for sea-level projections specific to: 1) time step and 2) percentile using Esri ArcGIS (v. 10.2.1). Sea-level data were exported from MATLAB (v. R2013b) in XYZ ASCII format as a .CSV file. XY data were provided in geographic coordinates, and data were imported to ArcGIS > add xyz data in a geographic projection (World Geodetic System 1984 or WGS84). Point data were converted to a point feature class within a geodatabase in ArcCatalog by right-mouse clicking the point data file > Create Feature Class > From XY table. Point data were then interpolated to generate a surface of sea-level estimates using natural neighbor interpolation: ArcToolbox > Toolboxes > 3D Analyst Tools.tbx > Natural Neighbor (option: automatic densification). The resulting raster surface provided cells at 0.66 degrees (or approximately 74,000 m) resolution, and was imported into a geodatabase (with sea-level and VLM rasters). Sea-level grid values were used in calculations of coastal response type probabilities. Extraction of sea-level values for use in modeling is described in metadata sections to follow. This process step and all subsequent process steps performed by Sawyer Stippa.

Vertically converted elevation datasets from NAVD 88 to mean high water (MHW) were included to take into account tidal variation. The vertical datum was adjusted to MHW along the coastline using a one-step conversion grid generated through the National Oceanic and Atmospheric Administration (NOAA) Vertical Datum Transformation (VDatum, v. 3.1) tool which was added to the NED data in ArcGIS (v. 10.2.1) > Spatial Analyst Tools > Plus. Vertical conversion was not deemed necessary for the CRM data (see vertical accuracy). The conversion grid is provided at a resolution of 0.0017 degrees, which is the full resolution as supplied with VDatum. Vertical datum shifts to MHW were propagated landward via Euclidean allocation, which smoothes vertical ledges at edges of VDatum zones. Elevation grid values were used in calculations of coastal response type probabilities. Extraction of elevation values for use in modeling is described in metadata sections to follow.

Vertical land motion estimates in the form of point data (rates) were used to interpolate a surface to provide continuous VLM rate estimates across the region using ArcGIS v. 10.2.1. Point data (XYZ format) were downloaded from Sella and others, 2007 and Zervas and others, 2013 and imported to Esri ArcGIS (add XYZ data) and a feature class of the points was created. The points (geographic coordinates, projected in WGS84) were used to interpolate a raster surface using the multiquadratic radial basis function: ArcToolbox > Toolboxes > Geostatistical Analyst Tools.tbx > Radial Basis Functions (selected tool options: Multiquadratic Function). The resulting surface was gridded at a per cell resolution of 0.23 degrees (or approximately 24,000 m) and imported to the geodatabase containing sea-level projections and elevation. VLM grid values were used in calculations of coastal response type probabilities. Extraction of VLM values for use in modeling is described in metadata sections to follow.

Land cover data was assigned to generalized categories of similar ecology, geomorphology and/or anticipated physical response to sea-level rise. Every land cover class, as determined from ESM Plus name and/or description, was assigned a code corresponding to one of the following six categories: 1) Subaqueous: bays, lakes, rivers, marine and estuarine subtidal, and deepwater. 2) Marsh: Salt and freshwater marshes, bogs, swamps, fens, wetland forests, intertidal aquatic beds and reefs. 3) Rock: rocky outcrops and shores, marine and estuarine intertidal rock bottom. 4) Beach: dune and swale/sandy beach (including bluffs), marine and estuarine intertidal unconsolidated shore. 5) Forest: forests, woodlands, grasslands, agricultural, shrublands. 6) Developed: all NLCD developed classes (open space, low, medium, and high density), roads, active and abandoned railroad tracks. Assignments were made by joining an Excel table of the classified land cover categories with the attribute table of the raster land cover data using Esri ArcGIS (v. 10.2.1): ArcToolbox > Toolboxes > Data Management Tools.tbx > Add Join.

Generate a database of points at land cover locations from which model inputs might be culled. An ESM Plus database (EPD) of points was created using the extent of the converted 1/3 arc-second DEM at -15 to 10 m. The 1/3 data were selected because they provided continuous regional coverage. The 1/3 DEM data was reclassified by converting all data within the stated ranges (-15 to 10) to 1 and all other to null: ArcToolbox > Toolboxes > Spatial Analyst Tools.tbx > Reclassify (options: min to -15 no data. 10 to max no data). The resulting raster was converted to a polygon to be used in clipping the ESM Plus raster: ArcToolbox > Toolboxes > Conversion Tools.tbx > Raster To Polygon (option: no simplify). The ESM Plus raster clip was converted to points using Esri ArcGIS software (v. 10.2.1) to form the EPD: ArcToolbox > Toolboxes > Conversion Tools.tbx > Raster to Point. Note the resulting EPD retained the land cover values after conversion.

Use elevation and land cover information to delineate study area extent. Elevation and land cover information was used to define the coastal extent as follows: a) at or below 10 m and at or above -10 m elevation for the BEACH land cover category and b) at or below 5 m and at or above -10 m for all other general land cover categories. The delineation of this area required the creation of several polygons generated from these raster inputs that were then merged together to define the full study area extent. Footprints of the NED and CRM datasets were created using polygon clips of the EPD. The NED footprint was created by reclassifying the 1/3 arc-second DEM within 0 to 5 m: ArcToolbox > Toolboxes > Spatial Analyst Tools.tbx > Reclassify (options: min to 0 no data. 5 to max no data). The resulting raster was converted to a polygon: ArcToolbox > Toolboxes > Conversion Tools.tbx > Raster To Polygon (option: no simplify). A polygon of the ESM Plus BEACH land cover category was created by reclassifying the ESM Plus to only show the beach classification: ArcToolbox > Toolboxes > Spatial Analyst Tools.tbx > Reclassify (options: 1 to 2 no data. 4 to 6 no data). The resulting raster was converted to a polygon: ArcToolbox > Toolboxes > Conversion Tools.tbx > Raster To Polygon (option: no simplify). The beach polygon was clipped further using a reclassified 1/3 arc-second DEM polygon at 5 to 10 m. The beach polygon was merged with the reclassified 1/3 arc-second DEM polygon at 0 to 5m: ArcToolbox > Toolboxes > Data Management Tools.tbx > Merge. A third polygon was merged with the NED footprint that included a reclassified segment of the NAVD88 1/3 arc-second DEM (prior to MHW conversion) at 0 to 3.5 m: ArcToolbox > Toolboxes > Data Management Tools.tbx > Merge. This segment was used to provide continuous coverage starting slightly below MHW. The CRM footprint was created by reclassifying the CRM dataset within -10 to 1 m: ArcToolbox > Toolboxes > Spatial Analyst Tools.tbx > Reclassify (options: min to -10 no data. 1 to max no data). The resulting raster was converted to polygon: ArcToolbox > Toolboxes > Conversion Tools.tbx > Raster To Polygon (option: no simplify).CRM and NED footprints were used separately to allow for faster processing times, though future performance increases may allow for the combination of each footprint.

Cull points from the EPD within the study area extent and extract grid values at coincident locations from VLM, elevation, sea-level projection, and land cover rasters. A NED point database was created by clipping the EPD with the NED footprint. The resulting dataset was consequently divided into 10 subregions to make processing more manageable: ArcToolbox > Toolboxes > Data Management Tools > Create Fishnet (options: 10 rows 1 column): ArcToolbox > Toolboxes > Analysis Tools.tbx > Clip. The CRM point database was created by clipping the EPD with the CRM footprint: ArcToolbox > Toolboxes > Data Management Tools.tbx > Clip. CRM database points that overlapped NED database points were removed: ArcToolbox > Toolboxes > Data Management Tools.tbx > Select Layer By Location (options: intersect new selection): ArcToolbox > Toolboxes > Data Management Tools.tbx > Delete Rows. The resulting CRM point database was divided into 59 regions based on increments of 500,000 points using a Python, v. 2.7 script. For each NED and CRM subregion, raster data from each sea-level projection year and percentile range, vertical land movement, elevation and land cover layers were extracted at each point location: ArcToolbox > Toolboxes > Spatial Analyst Tools.tbx > Extract Multi Values To Points. A script (Python, v. 2.7) was developed to populate an elevation column with 1/3 arc-second data that were overwritten with 1/9 arc-second values if present for all NED subregions. An error column was also populated to distinguish 1/3 arc-second and 1/9 arc-second values for each NED subregion. The NED subregions were further truncated to the 0 to 10 m range for BEACH land cover category and 0 to 5 m range for all other land cover categories, this time based on the extent defined by the 1/9 arc second data where present. CRM points outside of the -10 to 1 m range were removed. Geographic coordinate information (latitude and longitude, WGS84) for each point was also added to the NED and CRM subregions before being exported to csv tables.

The Bayesian network (BN) model was used to generate region-wide geospatial values for CR outputs using elevation data, sea-level projections for selected time intervals, vertical land movement rates, and land cover information. A BN offers the capability to generate robust, spatially-explicit predictions and a corresponding confidence interval with every predicted outcome. Parameters that influence coastal change may all be included in the network design to make a prediction (an updated conditional probability) given a set of observations. The BN results in two outputs: 1) the likelihood of each adjusted elevation range (PAE) being reached given a range of sea-level rise scenarios and 2) the probability of dynamic vs. static coastal response type (CR). Full details of the BN, including conceptual diagrams and probabilistic equations used to describe the model may be found elsewhere (Lentz and others, 2015). Point data from each subregion were exported as tables and imported (sans headers, object ID and shape geometry fields) into MathWorks MATLAB software (v. R2013B) for processing in MATLAB and Netica. Input datasets consisted of 21 primary columns: SL1020, SL1030, SL1050, SL1080, SL2520, SL2530, SL2550, SL2580, SL7520, SL7530, SL7550, SL7580, SL9020, SL9030, SL9050, SL9080 (sea-level predictions where the first two digits indicate percentile range, the second two the prediction year in 20XX); VLM (vertical land motion); DEM (1/3 or 1/9 arc-second DEM values or CRM DEM values); LC (land class category 1-6); latitude; and longitude.

Equations which previously appeared in this process step (describing Bayes theorem) have been modified to provide more explanatory content specific to our Bayesian network function. The updated equations appear in the Open-File Report cross-reference Lentz and others, 2015, version 2. The modifications to the equations do not change the results, but just provide a bit more detail and specifics on the probabilistic approach being used. Because of special characters and other considerations, please reference the Open-File Report for these equations.

CR prediction matrices from completed model runs were combined with lat/long of input matrices and additional Matlab code was used to convert outputs to shapefile format. The resulting shapefiles were brought back into ArcGIS as rasters using a python script (v. 2.7). Point to raster conversion did not require interpolation because the horizontal spacing between input points retained the 30 m horizontal resolution of the original land cover grid. All corresponding rasters were then merged ArcToolbox > Toolboxes > Data Management Tools.tbx > Mosaic to New Raster to generate a single raster of region-wide coverage.

The metadata was modfied to edit the 8th process step which contained equations describing the Bayes theorem specific to the Bayesian network function used in this work. The updated equations appear in the cross-reference Lentz and others, 2015, version 2.

USGS Thesaurus keywords added to the keyword section. Additionally, a DOI number was assigned and the metadata updated to reflect this. Additionally, the page containing the data was moved to a new location, with the data file downloads pointing to the new location.

Added keywords section with USGS persistent identifier as theme keyword (20200908). Fixed a USGS Thesaurus term (20240712).