Seabeach Amaranth Presence-Absence Data, Assateague Island National Seashore, 2014

This and subsequent process steps were performed in ArcGIS Pro’s Python window. The original code was written by Emily J. Sturdivant, and the code was translated into process steps by Julia L. Heslin.

1) DEM Processing: A few steps were required to process the DEM from raw data (DEM_raw). It is important to note that all source data were projected to UTM zone 18N NAD83 prior to processing. The DEM tile(s) were copied from their native GeoTIFF format to a file geodatabase (Data Management Tools >> Raster >> Raster Dataset >> Copy Raster. The study area may be comprised of multiple DEM tiles which were mosaicked into a single mosaic dataset (Data Management Tools >> Raster >> Raster Dataset >> Mosaic) (DEM_preprocessed) . Any NoData values that were present in the DEM were set to null (Spatial Analyst Tools >> Conditional >> Set Null). Unless otherwise noted, the following steps were performed by one person, Julia L. Heslin.

2) Morphology processing: A dataset with morphological characteristics, was available for the extraction of multiple variables, specifically dune toe, dune crest, and shoreline point location information for Assateague Island. For 2014, these points were already converted and available from this data release (https://www.sciencebase.gov/catalog/item/5daa377fe4b09fd3b0c9ce84).

3) Barrier island boundary processing: A polygon delineating the barrier island boundary was created using the DEM (DEM_preprocessed) and morphology source data (SL_pts) (see process steps 1 and 2).
Several tools were needed to convert the previously created boundary data to the barrier island boundary used in this project. Using the Con tool (Geoprocessing Tools >> Spatial Analyst >> Conditional >> Con), the DEM was reclassified where anything above the local MHW value was assigned as 1. The raster was converted to a polygon feature class (Geoprocessing Tools >> Conversion >> From Raster >> Raster to Polygon) and aggregated (Geoprocessing Tools >> Cartography >> Generalization >> Aggregate Polygons) to merge any features at least 300 square meters in size within a 10 meter distance of each other. These same steps were performed to create polygon features at the mean tide level (MTL), which is MHW (0.34) plus the local mean low water (MLW) (-0.13) divided by two. The purpose of this is to capture the extent of the bay side shoreline of the island.
The MHW and MTL polygons were then combined into one bayside polygon. Before this, inlet lines were manually drawn to capture the edges of the barrier island using the basemap imagery from Esri. These inlet lines were saved as a feature class (inlet_lines). Then, the Symmetrical Difference tool (Geoprocessing Tools >> Analysis >> Overlay >> Symmetrical Difference) was used to find where the MHW and MTL polygons intersected and keep features where they do not overlap (symdiff). Using the Feature to Polygon tool (Geoprocessing Tools >> Data Management >> Features >> Feature to Polygon), any extraneous shoreline polygon features were split using the symdiff and inlet_lines feature classes (split). This split feature class was then spatially joined with the SL_pts features (extracted from process step 2). Any of the split polygon features that completely contained shoreline points were kept (split_join). Then, any features from the symdiff and split_join feature classes that overlapped were removed using the Erase tool (Geoprocessing Tools >> Analysis >> Overlay >> Erase).
The bayside polygon was combined with the MHW polygon using the Union tool (Geoprocessing Tools >> Analysis >> Overlay >> Union). All features from this output were dissolved to make a single part polygon (Geoprocessing Tools >> Data Management >> Generalization >> Dissolve). Any extraneous features that did not reflect the barrier island boundary (barrier_boundary) were manually deleted.
The last step in creating the barrier island boundary polygon was to snap the polygon to the shoreline points within a 25 meter distance where they don’t already match. The Densify tool (Geoprocessing Tools >> Editing >> Densify) was used to add new vertices to the barrier island boundary at 25 meter increments. These vertices were then snapped to the shoreline points as long as they were within 25 meters of each other. After this step, any holes that existed in the barrier island boundary polygon were deleted using the Eliminate Polygon Part tool (Geoprocessing Tools >> Data Management >> Generalization >> Eliminate Polygon Part). Any part less than 99% of the total polygon area that was contained in the barrier island boundary was deleted. This is to ensure that all holes in the polygon were eliminated. Finally, the Dissolve tool was used to aggregate the individual polygons into a single multi-part feature.

4) DEM raster extraction: Using the barrier island boundary as the extent (see step 3), the DEM was then resampled to a five meter pixel size using a bilinear resampling technique (Data Management Tools >> Raster >> Raster Processing >> Resample). It was then adjusted to the local mean high water (MHW) for the study area (0.34m) by subtracting 0.34 from the DEM (3D Analyst Tools >> Raster Math >> Minus)(local_MHW_DEM).
Additional rasters were calculated using the DEM, including aspect (3D Analyst Tools >> Raster Surface >> Aspect) and slope (with the output measurement set as percent rise) (3D Analyst Tools >> Raster Surface >> Slope).

5) Morphology raster extraction: Using the barrier island boundary as the extent (see step 3), distance rasters using the Euclidean Distance tool (Geoprocessing Tools >> Spatial Analyst >> Distance >> Euclidean Distance) were created using the preprocessed morphology points.
Two rasters were created related to the shoreline morphology: distance to the full island shoreline; and distance to the oceanside shoreline. In order to retrieve the distance to the full island shoreline, the barrier_boundary was converted to a line using the Polygon to Line tool (Geoprocessing Tools >> Data Management >> Features >> Polygon to Line). This line feature was used to calculate the distance to the full boundary shoreline. To calculate the distance to the oceanside shoreline, first the barrier_boundary and inlet_lines were converted using the Feature to Line tool (Geoprocessing Tools >> Data Management >> Features >> Feature to Line) indicating where the inlets split the shoreline. Then, the SL_pts file was spatially joined to the split shoreline feature class. Any line features that completely contained the SL_pts were kept. The shoreline was then aggregated to one oceanside shoreline feature using the Dissolve tool (Geoprocessing Tools >> Data Management >> Generalization >> Dissolve) and used to generate the oceanside shoreline distance raster (distance_to_oceanside_shoreline).
Additional distance rasters using the dune toe, dune crest, and inlet features were created. The DL_pts and DH_pts were used to derive the distance_to_DL and distance_to_DH rasters (with 200 meters being the maximum threshold value). The inlet at Ocean City (which was manually selected and saved as a separate feature class) was used to create the distance_to_OceanCity_inlet raster.

6) Field data processing: Field data documenting locations of seabeach amaranth plants (field_data_raw) were used in the creation of this dataset. The field dataset has records of seabeach amaranth locations, grazed status, and other attributes from 2001-2018.
The seabeach amaranth field data, originally formatted as an Excel spreadsheet, were imported as a Python Pandas data frame (amaranth_data_frame). Since the field data contains data from multiple years, rows from the year in question were selected.
The data frame was converted to a feature class by first keeping all columns from the data frame and converting them to a NumPypy array. The NumPyArrayToFeatureClass function (Python >> Data Access Module >> Functions >> NumPyArrayToFeatureClass) was used to convert the array to a feature class.
A training dataset with random points and plant observations (training_pts) was created first by generating random points using the Create Random Points tool (Geoprocessing Tools >> Data Management >> Sampling >> Create Random Points). The random points distributed within the barrier island boundary polygon (created in step 3) no less than 5 meters apart. The random points were merged with the plant observations into one feature class. A flag field called plant_presence was created to distinguish plants from random points (a value of 1 indicates plant presence and 0 indicates a random point).
A feature class of plants from the previous year (lastyr_pts) was created by selecting the rows in the field data data frame from the previous year, which for this dataset was 2011 The selected rows in the data frame were converted to a NumPy array then to feature class using the NumPyArrayToFeatureClass function (Python >> Data Access Module >> Functions >> NumPyArrayToFeatureClass). A tally field (called cnt) was added to the lastyear_pts feature class.

7)Field data raster extraction:
A raster dataset indicating the number of last year’s (2013) points within a 30 meter radius of a grid cell was created. The NbrCircle tool was used to specify this radius value (Python >> Spatial Analyst module >> Classes >> Neighborhood classes >> NbrCircle). The Point Statistics tool (Geoprocessing Tools >> Spatial Analyst >> Neighborhood >> Point Statistics) was used to create a raster of the number of last year’s points in the 30 meter radius.
Distance and direction rasters to the nearest plant from the previous year were also created. The Euclidean Distance tool was used to calculate: a distance raster with 200 meters being the maximum threshold distance the raster values cannot exceed (distance_to_plant); and a direction raster (direction_to_plant). Another direction raster relative to the adjacent mean high water shoreline was generated as well (direction_to_shore).
The difference between direction_to_plant and direction_to_shore was calculated with the Minus tool (3D Analyst Tools >> Raster Math >> Minus) (direction_to_plant_alongshore). The values from this raster indicate updrift (0 to 180) and downdrift (less than 0 or greater than 180). The dir2plant_alongshore raster was reclassified into these two categories (1 for updrift and –1 for downdrift) (Geoprocessing Tools >> Spatial Analyst >> Reclass >> Reclassify). The reclassified raster was multiplied by the distance_to_plant raster to get the final distance raster (distance_to_plant_final), where positive values indicate updrift and negative indicate downdrift.

8)Vegetation processing: The original vegetation data (vegetation_raw) were in a shapefile format and were first imported into a file geodatabase for processing. The vegetation feature class was then converted to a raster (Conversion Tools >> To Raster >> Polygon to Raster) (vegetation_raster).

9)Raster to spreadsheet: From the raster datasets created from the process steps above(vegetaion, elevation, and distance rasters), the pixel values were extracted to the corresponding training points feature class (training_pts) (Spatial Analyst Tools >> Extraction >> Extract Multi Values to Points). Extracted values can be found in the field attributes of the final csv.
The training points feature class (training_pts) was then converted to a data frame by first converting the points to a NumPy array (Python >> Data Access Module >> Functions >> FeatureClassToNumPyArray), then to a Python dictionary, then a data frame using the Pandas library.
The plants data frame (amaranth_data_frame) described in process step 6 was then joined with this data frame described in the previous paragraph (data frame from training_pts). Fields for lat and lon were added as attributes to this data frame and calculated using the PyProj library (https://pyproj4.github.io/pyproj/v2.2.0rel/api/transformer.html#pyproj-transform). The data frame was then written to an Excel spreadsheet using the Python Pandas library (https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.to_excel.html).

10) Conversion of csv to shapefile. The .csv file was converted to a shapefile by right clicking the dataset in ArcMap, and selecting "Display XY Data" where x_utm and y_utm were chosen for the x and y fields respectively, and the coordinate system NAD_1983_UTM_Zone_18N was selected for the dataset.

11) Coordinate system updated from projected coordinates in "NAD_1983_UTM_Zone_18N" to geographic "GCS_North_American_1983" for publication using ArcToolbox > Data Management Tools > Projections and Transformations > Project.