Conceptual marsh units of salt marshes on the Eastern Shore of Virginia

This process step and all subsequent process steps were performed by the same person, Kate Ackerman, in ArcGIS Pro (ver. 2.8.8), unless otherwise stated. For simple operations the name of the geoprocessing tool used is given in capital letters; for complex operations, the name of geoprocessing tool used is given in capital letters and any critical parameters used are given in parentheses, separated by a semicolon, immediately after the tool name. The input and output file names are provided in [square brackets] when necessary. Units for length and area calculations are meters (m) and square meters (m^2) unless otherwise stated.
a) Set the data frame coordinate system and projection to NAD 1983 UTM Zone 18N.
b) Prepare elevation dataset. Combine all downloaded tiles into one raster, using MOSAIC TO NEW RASTER, and PROJECT RASTER to NAD 1983 UTM Zone 18N.
c) Create a polygon that covers the area of interest [extent.shp]. CLIP the mosaicked elevation dataset created in step (b) by [extent.shp] to create [elev_mosaic.tif].

This process step creates two mask polygons to define analysis boundaries over the area of interest. In order to account for the possible influence of the surrounding terrain on the analysis, a mask area larger than the salt marsh extent is created [UVVR_mask.shp]. This larger mask area is clipped to the boundaries of the salt marsh area to create [MU_mask.shp]. This process step explains how [UVVR_mask.shp] and [MU_mask.shp] are created and edited to determine which interior polygons are removed.
a) Create a shapefile of the salt marshes as defined by NWI. EXPORT features (use the same coordinate system as the data frame) from NWI dataset of all wetlands for the area (referred to here as [NWI_entire.shp]) after SELECT (ATTRIBUTE LIKE 'E2%EM%' OR ATTRIBUTE LIKE 'E2AB3%' OR ATTRIBUTE LIKE 'E2EM%' OR ATTRIBUTE LIKE 'E2SS%' OR ATTRIBUTE LIKE 'E2US4%') to select from estuarine intertidal areas of 1) emergent wetland, 2) scrub-shrub area, 3) rooted vascular aquatic bed, 4) organic unconsolidated shore classes to create [NWI_exclusive_select.shp].
b) Create the preliminary masks that will be edited in the subsequent steps to make the final masks. BUFFER salt marsh classes (Input features=[NWI_exclusive_select.shp]; Linear unit=5 meters; Dissolve type=All) to obtain the mask polygon [buff5m.shp] and ELIMINATE POLYGON PART (Input features= [buff5m.shp]; Condition=Percentage; Percentage=99; Eliminate contained part only) to remove small interior polygons to obtain [UVVR_mask01.shp]. Apply BUFFER (Input features=[UVVR_mask01.shp]; Linear unit= -5 meters; Dissolve type= All) to obtain [MU_mask01.shp].
c) EXPORT features from NWI dataset after SELECT (ATTRIBUTE NOT LIKE 'E%' OR ATTRIBUTE LIKE '%FO') to obtain features that are not estuarine or are forested estuarine [NWI_inclusive_select.shp]. MULTIPART TO SINGLEPART to separate individual polygons [NWI_inclusive_single.shp].
d) To exclude the voids in the NWI map during the analysis, these interior polygons need to be removed from the mask. ERASE [NWI_exclusive_select.shp] from [extent.shp] and MULTIPART TO SINGLEPART to get [NWI_exc_voids_single.shp].
To create the final mask shapefiles ([MU_mask.shp] and [UVVR_mask.shp]), areas that should not be included are removed (e.g., forest areas), and areas that should be included are incorporated (estuarine and marine wetland systems that are completely surrounded by salt marsh (e.g., interior ponds)):
e)SELECT BY LOCATION from [NWI_inclusive_single.shp] that are completely within the [MU_mask01.shp] and SELECT BY LOCATION from [NWI_exc_voids_single.shp] that are completely within the [MU_mask01.shp]. UNION all feature sets to get polygons to be excluded [erase_inside.shp] from the final masks.
f) Spatially join the [NWI_entire.shp] to [erase_inside.shp] in order to assign wetland attributes to the [erase_inside.shp] polygons. SPATIAL JOIN (target features= [erase_inside.shp]; join features= [NWI_entire.shp]; output feature class= [erase_inside_NWI_info.shp]; join operation= one to one; keep all target features; match option= intersect; search radius= -1 m).
g) Select the polygons that are marine or non-forested estuarine wetlands from [erase_inside_NWI_info.shp]: SELECT (ATTRIBUTE LIKE 'E1%' OR ATTRIBUTE LIKE 'E2AB%' OR ATTRIBUTE LIKE 'E2RF%' OR ATTRIBUTE LIKE 'E2SB%' OR ATTRIBUTE LIKE 'E2RS%' OR ATTRIBUTE LIKE 'E2US%' OR ATTRIBUTE LIKE 'E2EM%' OR ATTRIBUTE LIKE ‘E2SS%’ OR ATTRIBUTE LIKE ‘M1%’ OR ATTRIBUTE LIKE ‘M2%) to select all marine systems, estuarine subtidal systems, estuarine intertidal areas of aquatic bed, reef, streambed, rocky shore, unconsolidated shore, emergent or scrub-shrub [erase_inside_keepIntPonds.shp].
h) Remove the polygons from [erase_inside_NWI_info.shp] that should be retained (these are the interior ponds): ERASE (input feature= [erase_inside_NWI_info.shp]; erase feature= [erase_inside_keepIntPonds.shp]; output feature= [erase_inside_final.shp]).
i) If any additional polygons need to be removed from the mask, include them in a new feature layer [erase_mask_manual.shp] and MERGE with [erase_inside_final.shp].
j) Remove any unwanted polygons from the original [MU_mask01.shp] to make the final [MU_Mask.shp]: ERASE (input feature= [MU_mask01.shp]; erase feature= [erase_inside_final.shp]; output feature= [MU_mask.shp]).
k) BUFFER [MU_mask.shp] by 5m (side type= full; end type= round; method: geodesic; dissolve type: single feature) to obtain [UVVR_mask.shp].

Define preliminary boundaries for marsh units by basin analysis:
a) EXTRACT BY MASK from elevation dataset [elev_mosaic.tif] using [UVVR_mask.shp] polygon.
b) FILL (no Z limit) sinks in extracted elevation raster.
c) Calculate FLOW DIRECTION (do not force edge cells to flow outward) raster.
d) Calculate BASIN based on flow direction raster [basin_buff5m.tif].
e) Apply MAJORITY FILTER (Number of neighbors=8) to clean the raster [basin_buff5m.tif] and convert RASTER TO POLYGON [basin_b5m_poly.shp]. Create a 1-meter buffered marsh units mask [MU_mask_buff1m.shp] and use this to CLIP the polygon to create [basin_b5m_polyclip.shp]. This additional buffer is necessary to prevent any polygon smoothing algorithm from changing the exterior boundaries and will be removed by clipping to the [MU_mask.shp] later. Using [basin_b5m_polyclip.shp] as the input, MULTIPART TO SINGLEPART and ADD GEOMETRY ATTRIBUTES(Geometry properties=Area_geodesic; Area unit=Square_meters) to get [basin_single.shp].

This step merges the "orphan marsh units," units that are smaller than 5000 m^2, to the nearest to "parent marsh units," units that are larger than 5000 m^2, to obtain the preliminary marsh units. A Python script that iteratively calls a list of ArcGIS tools to aggregate orphan units on parent units is used.
a) Create orphan marsh unit features by exporting marsh units from [basin_single.shp] with an area less than 5000 m^2 [inputfeatures.shp].
b) Create parent marsh unit features by exporting marsh units with an area greater than or equal to 5000 m^2 [nearfeatures.shp].
c) Run Python script (hydUnitloop.py). At each iteration step, the script uses NEAR tool to find orphan marsh units within 1 meter of a parent marsh unit and merges them to the parent using UNION and DISSOLVE tools. The orphan units merged with a parent unit are removed from the orphan units dataset. The script will iterate until there is no change in the number of parent or orphan marsh units. The remaining orphan marsh units are those that are more than 1 meter away from any parent. DISSOLVE them to create larger orphan units by grouping the orphans attached together into larger units. MULTIPART TO SINGLEPART to have single part polygons with unique identifiers. UNION(with gaps) single part features with parent units to obtain preliminary marsh units [pMUopt1.shp]. SMOOTH POLYGON(Smoothing algorithm=PAEK; Smoothing tolerance=15 meters) to get the smoothed marsh units polygon [pMUopt1_smooth.shp].
d) CLIP the extra buffer around the smoothed marsh units using the marsh units mask [MU_mask.shp]. If this results in any additional marsh units with an area less than 5000 m^2, dissolve into the nearest parent units repeating the same steps above to obtain [pMUopt1_clip_final.shp].
e) Compute the final marsh units by removing any artifacts from [pMuOpt1_clip_final.shp] using ELIMINATE and REPAIR GEOMETRY tools. ELIMINATE(Expression="Area_geo" less than 900; Eliminating by border) polygon slivers smaller than 900 m^2 by merging them to the main polygon. MULTIPART TO SINGLEPART to get conceptual marsh units [CMU.shp]. This only applies to polygons that are bordering other polygons. Polygons that are not bordering other polygons are not affected as those are more likely to be real features as opposed to geoprocessing artifacts. To [CMU.shp], ADD FIELD “FID_CMU” and CALCULATE FIELD as FID_CMU= FID x 1. Delete fields from [CMU.shp] so that only area_geo and FID_CMU remain. REPAIR GEOMETRY of [CMU.shp].
f) DISSOLVE the conceptual marsh units to get the final outlines for the marsh complex [CMU_mask_final.shp].

Determine the vegetated and unvegetated areas by image processing.
a) MOSAIC TO NEW RASTER the NAIP tiles.
b) EXTRACT BY MASK from the elevation raster [elev_mosaic.tif] and the NAIP imagery raster using [UVVR_mask.shp] polygon and append “UVVR_mask” to the file name. Rescale the elevation raster values to the same range with 8-bit NAIP imagery (0 to 255) using min-max scaling in RASTER CALCULATOR (map algebra expression: (("elev_mosaic_UVVR_mask.tif"-(min value of the input raster))*(255-0)/(max value of the input raster)-( min value of the input raster)))+0). Perform ISO CLUSTER UNSUPERVISED CLASSIFICATION with 5 input bands: the Near Infrared, Blue, Green, and Red from the NAIP imagery and the rescaled elevation raster; with 32 classes and a minimum class size of 5000 cells. RECLASSIFY the classified raster by visually comparing the NAIP imagery to obtain the unvegetated-vegetated raster [UVVc.tif]. In some regions, a couple of the 32 classes are classified as vegetated in one area but unvegetated in another area. For these classes, the area can be split into separate sections to allow one section to be classified as vegetated and the other as unvegetated.
c) Dissolve unvegetated regions smaller than a threshold value to the surrounding vegetated regions and vice versa to clean the UVV raster. For this purpose, use the REGION GROUP(Number of neighbors=4; Zone grouping method=Within; Add Link field to output) to get the [UVV_region.tif], and TEST("Count" < 9) to set a threshold value of 9 raster cells [lt9.tif]. Use RASTER CALCULATOR to toggle the value of the Link field in the region raster [Map algebra expression=Con("%lt9.tif%", ~Lookup("%UVVc_region.tif%","LINK"),Lookup("%UVVc_region.tif%","LINK"))], where the test raster indicates regions with areas smaller than the threshold.
d) RASTER TO POLYGON(Simplify polygons=False) and CLIP with the final marsh complex outline [CMU_mask_final.shp] to obtain [UVVc_filt_clip.shp].
e) REPAIR GEOMETRY, ADD GEOMETRY ATTRIBUTES(Geometry properties=Area_geodesic; Area unit=Square_meters), CALCULATE FIELD(Python_9.3 expression; Field name=gridcode, Expression=filt(!gridcode!, !AREA_GEO!), Code Block= def filt(gridcode, AREA_GEO):/if gridcode== 0 and AREA_GEO < 9:/return 1/else:/return gridcode) and DISSOLVE(Create multipart features=False) to get filtered, clipped and dissolved unvegetated-vegetated polygons [UVVc_filt_clip_diss.shp].
f) INTERSECT the final conceptual marsh units [CMU.shp] with unvegetated-vegetated polygons [UVVc_filt_clip_diss.shp], REPAIR GEOMETRY, and calculate the area of unvegetated and vegetated polygons in each marsh unit. To do this ADD GEOMETRY ATTRIBUTES (Geometry properties=Area_geodesic; Area unit=Square_meters) to get [UVV_poly.shp], and DISSOLVE(Input features=[UVV_poly.shp]; Dissolve fields=FID_CMU, gridcode; Statistics field=gridcode.MEAN) and REPAIR GEOMETRY to get [UVV_poly_diss.shp]. Add fields "TYP", "APGN_M2", and "ATOT_M2" for type of polygon indicating vegetated or unvegetated, surface area of the polygon, and total surface area of the marsh unit, respectively, and calculate "APGN_M2" and "ATOT_M2" areas.
g) Visually inspect and manually remove marsh units that are fully forested due to polygon artifacts that are a result of buffering. Remove marsh units that are completely unvegetated, unattached to other marsh units (i.e., islands) and less than 25,000 m2, as these are artifacts of the NWI dataset. Remove marsh units that are linear features (edges parallel to each other) and unattached to other marsh units (i.e., islands), as these are artifacts of the NWI dataset.

Finalize the conceptual marsh units.
a) Add field "COLORID" to cycle through integers from 0 to 9 to be used for coloring the marsh units. This field is used for display purposes only to distinguish between marsh units.
b) Add field "FLG" to flag anomalous marsh units based on absence of vegetation, elevation, surface area and tidal range extrapolation. The flags are calculated in their respective shapefile (i.e., elevation flags are calculated in the elevation shapefile). Set FLG to (-1) no vegetated area in the marsh unit; (-10) marsh unit elevation higher than the 99.8 percentile; (-100) marsh unit elevation less than the 0.2 percentile; (-200) greater than 25 percent of the marsh unit elevation is hydro-flattened; (-1000) marsh unit surface area less than 900 m^2; (-10000) mean tidal range of marsh unit is extrapolated; and (0) no flag. Combination of negative values indicates a combination of flags. For example, -1001 indicates no vegetated area and area less than 900 m^2.
c) Rearrange field names and change the projection for better performance of web services with online base maps. PROJECT(Output coordinate system=WGS 1984 Web Mercator Auxiliary Sphere) the feature dataset to obtain the final conceptual marsh units [CMU_ESV].