Using Random Forest Classiﬁcation and Nationally Available Geospatial Data to Screen for Wetlands over Large Geographic Regions

: Wetland impact assessments are an integral part of infrastructure projects aimed at protecting the important services wetlands provide for water resources and ecosystems. However, wetland surveys with the level of accuracy required by federal regulators can be time-consuming and costly. Streamlining this process by using already available geospatial data and classiﬁcation algorithms to target more detailed wetland mapping e ﬀ orts may support environmental planning e ﬀ orts. The objective of this study was to create and test a methodology that could be applied nationally, leveraging existing data to quickly and inexpensively screen for potential wetlands over large geographic regions. An automated workﬂow implementing the methodology for a case study region in the coastal plain of Virginia is presented. When compared to veriﬁed wetlands mapped by experts, the methodology resulted in a much lower false negative rate of 22.6% compared to the National Wetland Inventory (NWI) false negative rate of 69.3%. However, because the methodology was designed as a screening approach, it did result in a slight decrease in overall classiﬁcation accuracy compared to the NWI from 80.5% to 76.1%. Given the considerable decrease in wetland omission while maintaining comparable overall accuracy, the methodology shows potential as a wetland screening tool for targeting more detailed and costly wetland mapping e ﬀ orts.


Introduction
Wetlands are a vital natural resource providing habitat for a variety of wildlife and plants, flood and storm surge protection, water quality improvement through treatment of runoff, and recharge of aquifers [1]. However, a significant number of wetlands in the U.S. have been destroyed or repurposed for agricultural or development purposes [2]. The need to protect wetlands is widely recognized and required by federal law and regulations, specifically through Section 404 of the Clean Water Act [3]. Section 404 of the Clean Water Act sets forth a goal of maintaining the nation's remaining wetland base by avoiding adverse impacts to these ecosystems. To comply with regulations, entities including state departments of transportation (DOTs) must consider potential impacts to wetlands in their infrastructure development projects. DOTs in particular, as well as other organizations, must sufficiently prove that a selected construction plan is the Least Environmentally Damaging Practical Alternative (LEDPA) by, among other tasks, providing wetland delineations [4]. The U.S. Army Corps of Engineers (USACE) evaluates these corridors as the governing authority in wetland permitting. that addresses these shortcomings and to implement the methodology as a wetland screening tool in a widely used commercial geographic information system (GIS). As a wetland screening tool, the methodology emphasizes minimizing false negative predictions (i.e., cases of wetland omission) while also maintaining a reasonable overall wetland accuracy. We obtained verified wetland delineations created by experts for a large region (33 km 2 ) in the coastal plain of Virginia to evaluate the methodology.
We trained a classification model on a subset of the verified wetland delineation dataset and then tested the classification model using a separate subset of the verified wetland delineation dataset to evaluate its accuracy. Finally, we compared the performance of the wetland tool against the NWI given that this is the standard wetland inventory available for supporting early environmental planning efforts in the absence of more detailed, costly, and time-consuming field surveys.

Study Area
The study area is in the Hampton Roads region of Virginia and is defined by the 33 km 2 limits of Virginia Department of Transportation (VDOT) wetland delineations and the 1363 km 2 processing extent encompassing delineations ( Figure 1). The processing extent was defined by the 12-digit hydrologic unit codes (HUC 12s), unique identifiers for watersheds in the U.S. [30] ( Figure 1A), that were intersected by the VDOT delineations ( Figure 1B,C). The study area resides within the Mid-Atlantic Coastal Plain ecoregion. The Mid-Atlantic Coastal Plain is characterized by primarily flat topography and poorly drained soils, and the characteristic land cover includes forest, agriculture, and wetlands [31]. According to the VDOT delineations, wetlands are widespread throughout the area with a wetland to non-wetland ratio of 0.4.
Water 2019, 11,1158 3 of 21 a widely used commercial geographic information system (GIS). As a wetland screening tool, the methodology emphasizes minimizing false negative predictions (i.e., cases of wetland omission) while also maintaining a reasonable overall wetland accuracy. We obtained verified wetland delineations created by experts for a large region (33 km 2 ) in the coastal plain of Virginia to evaluate the methodology. We trained a classification model on a subset of the verified wetland delineation dataset and then tested the classification model using a separate subset of the verified wetland delineation dataset to evaluate its accuracy. Finally, we compared the performance of the wetland tool against the NWI given that this is the standard wetland inventory available for supporting early environmental planning efforts in the absence of more detailed, costly, and time-consuming field surveys.

Study Area
The study area is in the Hampton Roads region of Virginia and is defined by the 33 km 2 limits of Virginia Department of Transportation (VDOT) wetland delineations and the 1363 km 2 processing extent encompassing delineations ( Figure 1). The processing extent was defined by the 12-digit hydrologic unit codes (HUC 12s), unique identifiers for watersheds in the U.S. [30] ( Figure 1A), that were intersected by the VDOT delineations ( Figure 1B and 1C). The study area resides within the Mid-Atlantic Coastal Plain ecoregion. The Mid-Atlantic Coastal Plain is characterized by primarily flat topography and poorly drained soils, and the characteristic land cover includes forest, agriculture, and wetlands [31]. According to the VDOT delineations, wetlands are widespread throughout the area with a wetland to non-wetland ratio of 0.4. Figure 1. The study area defined by the wetland screening tool processing extent and by the 13 12digit hydrologic unit codes (HUC 12) watersheds (A) that encompass the Virginia Department of Transportation (VDOT) delineation area. The entire region is shown as two parts (B and C) with aerial imagery provided by ESRI basemap services. Figure 1. The study area defined by the wetland screening tool processing extent and by the 13 12-digit hydrologic unit codes (HUC 12) watersheds (A) that encompass the Virginia Department of Transportation (VDOT) delineation area. The entire region is shown as two parts (B,C) with aerial imagery provided by ESRI basemap services.

Wetland Delineations (Training and Testing Data)
Wetland delineations obtained from VDOT were the only input data used that are not available on a national scale. Similar data would be required to apply the methodology for different regions. The wetland delineations were produced through field surveys completed for a corridor project and were provided in polygon vector format. Field surveys were conducted by professional wetland scientists during the period of May-July of 2013 and 2015, and were jurisdictionally confirmed by the USACE. These data were considered to be ground truth given their creation through manual surveying by trained analysts and subsequent confirmation of the delineations by the governing authority in wetland permitting. The delineated wetlands were originally categorized by wetland type and included emergent wetlands, forested wetlands, open water, scrub/shrub wetlands, and unconsolidated bottom wetlands. However, all types were generalized into a single "wetland" class, as the wetland screening tool was configured to detect characteristics shared by all wetland types to provide a first order nomination of likely wetland areas.
The delineations were randomly split into two different datasets, a training dataset and a testing dataset, with the majority of wetlands as testing data, as is typically done for building classification models. The training dataset, which included 10% of the delineated wetland area, was used to create the classification model. The testing dataset, which included the remaining 90% of the total wetland area, was used to evaluate the accuracy of the classification model. To create the training dataset, stratified random points were generated among the wetland and non-wetland classes, spanning the entirety of the delineated area. These points were then buffered to encompass 100 m 2 each and assigned wetland and non-wetland values on a per-pixel basis, according to the underlying delineation information. These processes resulted in a randomly dispersed training dataset with class proportions that are representative of the true land cover, as suggested by Millard and Richardson [22]. Figure 2 shows a portion of the training (A) and testing (B) datasets and Table 1 provides further details describing the entirety of the data. Wetland delineations obtained from VDOT were the only input data used that are not available on a national scale. Similar data would be required to apply the methodology for different regions. The wetland delineations were produced through field surveys completed for a corridor project and were provided in polygon vector format. Field surveys were conducted by professional wetland scientists during the period of May-July of 2013 and 2015, and were jurisdictionally confirmed by the USACE. These data were considered to be ground truth given their creation through manual surveying by trained analysts and subsequent confirmation of the delineations by the governing authority in wetland permitting. The delineated wetlands were originally categorized by wetland type and included emergent wetlands, forested wetlands, open water, scrub/shrub wetlands, and unconsolidated bottom wetlands. However, all types were generalized into a single "wetland" class, as the wetland screening tool was configured to detect characteristics shared by all wetland types to provide a first order nomination of likely wetland areas.
The delineations were randomly split into two different datasets, a training dataset and a testing dataset, with the majority of wetlands as testing data, as is typically done for building classification models. The training dataset, which included 10% of the delineated wetland area, was used to create the classification model. The testing dataset, which included the remaining 90% of the total wetland area, was used to evaluate the accuracy of the classification model. To create the training dataset, stratified random points were generated among the wetland and non-wetland classes, spanning the entirety of the delineated area. These points were then buffered to encompass 100 m 2 each and assigned wetland and non-wetland values on a per-pixel basis, according to the underlying delineation information. These processes resulted in a randomly dispersed training dataset with class proportions that are representative of the true land cover, as suggested by Millard and Richardson [22]. Figure 2 shows a portion of the training (A) and testing (B) datasets and Table 1 provides further details describing the entirety of the data.   Digital Elevation Models (DEMs) were used to derive regions where there is a high likelihood for pooled water, high soil moisture, and, therefore, probable wetland areas. Elevation data were obtained from the U.S. Geological Survey (USGS) National Map Download Client [32] in raster format. The highest available resolution tiles included 1/9th (~3 m) and 1/3rd arc-second (~9 m) data [33]. The vertical accuracy of these data can vary significantly across the United States, but is estimated to be around 1.5 m root mean squared error for the entire conterminous United States [34]. The 1/3rd arc-second DEM tiles covered approximately 300 km 2 of the processing extent and were resampled using the bilinear resampling technique to match the 1/9th arc-second resolution. As the DEM data provided the highest resolution information, 3 m was used as the pixel size for all subsequent vector to raster conversion and resampling of coarser rasters.
Multispectral imagery from the Landsat 8 OLI satellite was used to detect optical and vegetative wetland characteristics. Data were acquired from the USGS Earth Explorer in raster format with a 30 m resolution [35]. Imagery was chosen that provided cloud-free, spring conditions, as done by researchers in similar Coastal Plains studies (e.g., [12,23,36]), and that was collected on dates near the creation of VDOT delineations. The Landsat scene selected was collected on 4 April 2016 with 4.0% cloud cover. The date of 4 April was selected because, being the early spring season in Virginia, it is likely to have high soil moisture conditions. Visual analyses showed there was no cloud cover over the training and testing areas. Imagery preprocessing included conversion to top of atmosphere reflectance, according to USGS guidelines [37], and pixel resampling.
Federal Emergency Management Agency (FEMA) floodplain maps were used to identify areas of water inundation for heavy storm or flood events. FEMA 100-year floodplain maps are generated at a 1:12,000 scale and were downloaded from the FEMA Flood Map Service Center in polygon vector format [38]. All flood zones designated as 1-percent annual chance flood zones [39] were merged to create the 100-year floodplain zone for the study area. These data showed that floodplains occupy approximately 160 km 2 (12%) of the processing extent and 1.6 km 2 (5%) of the VDOT delineations. There were no significant cases of categorical maps spatially conflicting during the merging process.
Soil Survey Geographic Database (SSURGO) soil data, available through the United States Department of Agriculture (USDA), provided information on the location of hydric soils, which are characteristic of wetlands. Hydric soils are defined as soil that is formed under conditions of saturation, flooding, or ponding long enough during the growing season to develop anaerobic conditions in the upper horizon [6]. SSURGO maps are generated at a 1:12,000 scale and were downloaded in polygon vector format from the Natural Resources Conservation Service's Web Soil Survey [40]. According to SSURGO data, there are approximately 410 km 2 of hydric soils within the processing extent (30% of area) and 9 km 2 within the VDOT delineations (27% of area). There were no significant cases of categorical maps spatially conflicting during the merging process.
In addition to providing information on HUC 12 boundaries, national-scale hydrography data were used to estimate riparian zones. Riparian areas are lands that occur along flow channels and waterbodies, and most riparian areas either meet the Cowardin criteria for wetlands, which is a commonly used wetland classification system [5], or share many characteristics and functions with wetlands [41]. To estimate riparian boundaries, streams and waterbodies were downloaded from the National Hydrography Dataset (NHD). Streams used in this study consisted of NHD flowlines, in polygon vector format, categorized as stream/rivers, artificial paths, and connectors. Waterbodies used in this study consisted of NHD waterbodies, in polygon vector format, categorized as swamps/marshes, lakes/ponds, and reservoirs. Both datasets were obtained using the USGS National Map Download Client at a 1:24,000 map scale [32]. According to these data, there are approximately 1030 km of streams and 67 km 2 of waterbodies within the processing extent (5% of area) and 16 km of streams and 0.5 km 2 of waterbodies within VDOT delineations (1.5% of area).
Land cover data were used to identify certain land cover types where wetlands may be more likely to occur. These data were downloaded from the 2011 National Land Cover Database (NLCD) in raster format from the USGS National Map Download Client [32]. NLCD data were resampled from 30 m resolution to match that of the DEM using the majority resampling technique. Table 2 provides the summarized distribution of NLCD classes within the processing extent (left) and VDOT delineated area (right). Note that although the NLCD data include a wetland category, these areas refer to the NWI wetland designations. Lastly, NWI wetlands were included in the wetland tool workflow under the assumption that the NWI is incomplete rather than incorrect. NWI maps are generated at a 1:24,000 scale and were downloaded from the USFWS in polygon vector format [43]. According to the NWI, there are approximately 200 km 2 of wetlands within the processing extent (15% of area) and 3 km 2 within the VDOT delineations (9% of area).

Screening Methodology Design and Implementation
The wetland screening methodology can be described as a workflow ( Figure 3) that was implemented within ArcGIS 10.4 using their Model Builder tool. The workflow is segmented into three main sections, data preparation, training, and wetland prediction, each of which is described in detail in the following subsections. Implementing the methodology as a workflow within Model Builder allows it to be executed as a succession of geoprocessing and analysis tools able to transform nationally available data to potential wetlands maps with minimal user intervention. We selected Model Builder as the platform for implementation because (i) ArcGIS is a widely used software system within VDOT and likely other state DOTs as well, (ii) the graphical user interface (GUI) provides end users with workflow transparency, and (iii) workflows in this system can be exported to the Python programming language for further enhancement and automation.

Data Preparation
Data preparation includes satellite imagery processing, DEM processing, riparian zone processing, and additional processing. The methodology for executing these processes is given below and the results of each process are shown in Figure 4.

Data Preparation
Data preparation includes satellite imagery processing, DEM processing, riparian zone processing, and additional processing. The methodology for executing these processes is given below and the results of each process are shown in Figure 4. . Input datasets to the wetland screening tool, excluding NLCD data which is described in Table 2, and the dataset source. Note, while the NWI data is not included in the classification, it is incorporated in post-classification steps.  The satellite imagery processing consists of performing a Tasseled Cap Transformation (TCT) on the Landsat 8 OLI bands 2-7. The TCT is used to reduce dimensionality from several bands to a few bands associated with physical scene characteristics; specifically, brightness, greenness, and wetness, which are calculated as where w is a scalar used for weighing bands (Table 3), derived by Baig et al. [44]. Brightness measures are related to soil and albedo, greenness describes the presence of vegetation, and wetness describes water content [44]. Panels A, B, and C of Figure 4 show the resulting brightness, greenness, and wetness rasters, respectively. The DEM processing creates a binary local depression raster. The DEM is first conditioned using rasterized NHD streams with a 3 m resolution, per the DEM raster constraints. The DEM is artificially lowered by a large depth, in this case 100 m, where NHD raster pixels exist by executing the Map Algebra Expression (4).
Con (IsNull (Stream Raster)), DEM, (DEM − 100)) (4) resulting in the reconditioned streams DEM. After Expression (4) is executed, local depressions are filled using the Planchon and Darboux method [45]. This operation is generally used to remove small imperfections or barriers in topography for hydrologic flow path analysis; however, here it is used to identify areas of local depressions that could indicate wetlands in the low relief topography. Local depressions ( Figure 4D) are identified by flagging differences in the filled reconditioned stream (FRS) DEM compared to the original DEM using the Map Algebra Expression (5).
Con (FRS DEM == DEM, 0, 1) The riparian zone processing creates a binary riparian zone raster. As previously described, riparian areas are located along waterbodies and often exhibit wetland characteristics; however, riparian extents surrounding waterbodies can vary significantly, and there is no single agreed upon buffer distance to encompass riparian zones [41]. Here, practices associated with Virginia riparian zone policies, in the context of best management practices for water quality purposes, were analyzed to estimate a standard buffer width. The Chesapeake Bay Agreement recognizes a 15 m (50 ft) wide vegetative buffer extending from theoretical centerlines as forested riparian areas, enforced in lands within the Coastal Plains portion of the Chesapeake Bay watershed [46]. Thus, a 15 m buffer width extending from each side of NHD streams and waterbodies was used to represent the riparian zone, similar to the practice used by Hancock et al. [47]. The resulting polygon vector dataset was rasterized with the 3 m pixel resolution constraint and assigned values corresponding to either within or outside of the riparian zone ( Figure 4E).
In additional processing, binary rasters were created for FEMA, SSURGO, and NWI datasets. A binary floodplain raster is created where pixels within the 100-year floodplain are set to a value of one, and all other areas are set to a value of zero ( Figure 4F). A binary soil raster is created where pixels containing hydric soils are set to a value of one and pixels containing non-hydric soils are set to a value of zero ( Figure 4G). A binary NWI raster is created where pixels within wetland areas are set to a value of one and pixels within non-wetland areas are set to a value of zero ( Figure 4H). NLCD data are incorporated without additional processing.

Training
The training portion of the workflow involves two processes: composite inputs and random forest training. In the composite inputs process, the tool combines all prepared input data, with the exception of NWI data, into a multiband raster where each band stores the attributes of a single input. In the random forest training process, the composite raster and the training dataset are used as inputs to the Train Random Trees tool to extract wetland and non-wetland signatures. Train Random Trees executes the OpenCV implementation of the random forest algorithm within the ArcGIS environment [48]. This process also produces measures of variable importance based on the mean decrease in accuracy when a variable is not used in generating a tree [20]. The Train Random Trees tool allows users to specify values for the maximum number of trees, maximum tree depth, and maximum numbers of samples per class [49]. For this preliminary demonstration of the tool, these values were held constant at the ArcGIS recommended values of 50, 30, and 1000, respectively.

Wetland Prediction
The wetland prediction section includes classification and wetland expansion. During classification, the Classify Raster tool uses training information to classify the entire composite raster. The result of this operation is a binary wetland raster. In the wetland expansion section, a final tool output of the potential wetland raster is created by merging the classified wetlands with the NWI wetlands. As previously noted, NWI wetlands should not be the sole source for wetland screening; however, here NWI data are used to expand the classified wetland class in the interest of meeting the tool's objective to screen for areas where wetlands are likely to occur. The binary NWI raster is used to supplement classification results by applying the Map Algebra Expression (6).

Accuracy Assessments
Accuracy assessments were performed using a confusion matrix summarizing only predicted areas in the testing dataset that were not used for training. In each confusion matrix, testing data are represented along columns and either the screening tool (A) or NWI (B) data are represented along rows. Wetland tool and NWI performance were defined in terms of false negative rate, false positive rate, and overall accuracy. The false negative rate represents the tendency of models to omit true wetlands, and is defined as Equation (7).
False negative rate = False nonwetland predictions Total wetlands The false positive rate represents the tendency of models to overpredict wetlands and is defined as Equation (8).
False positive rate = False wetland predictions Total nonwetlands (8) Overall accuracy, Equation (9), represents the ability of models to correctly classify the total testing area, regardless of class.
Overall accuracy = True wetland predictions + True nonwetland predictions Total predicted area (9) Note that the false positive rate is equivalent to one minus the true negative rate, and the false negative rate is equivalent to one minus the true positive rate (i.e., the producer's accuracy). In addition, the Kappa statistic [50], which quantifies model performance taking into account the estimated performance of a random classifier, was also used to compare the wetland tool and NWI output. This set of performance criteria is not exhaustive but is appropriate to evaluate the testing area given the nearly balanced true wetland and true non-wetland classes [51]. Figure 5 shows the comparison between our screening tool predicted wetlands and NWI-designated wetlands for several areas within the larger study region. The "ground-truth" wetland delineations obtained by VDOT are included for reference. The levels of agreement between both wetland maps and VDOT delineations are represented as either wetland agreement, non-wetland agreement, false positive, or false negative in Figure 6. Note that scenes in Figures 5 and 6 include both training and testing locations (i.e., the entirety of the VDOT delineated area) for clarity, although only testing locations are included in performance assessments. In Figure 5, NWI maps show a tendency to underestimate VDOT wetlands, represented by relatively large distributions of false negative predictions in Figure 6. This trend is especially prominent in scene A6, where nearly all VDOT wetlands are missed by the NWI. Figure 5 scenes also demonstrate that in instances where NWI wetlands are correct, NWI wetland boundaries align relatively precisely with VDOT wetland boundaries. This trend is translated to significantly smaller distributions of false positive predictions in Figure 6. In contrast, the potential wetland areas predicted by our method more robustly encompass extents of VDOT delineations in scenes B1-B6 of Figure 5. However, the improved coverage of true wetlands was accompanied by imprecise wetland boundaries and wetland overestimation. In scenes B1-B6 of Figure 6, the wetland screening tool results in both larger areas of wetland agreement and false positive predictions surrounding predicted wetlands, relative to NWI mapping results. Figure 5. Examples of the wetland screening tool output for selected areas of the VDOT delineation area (B1-B6) compared to NWI mapping for the same scenes (A1-A6). VDOT wetlands, used to train and assess tool output, show the true distribution of wetlands. The location of scenes within the study area is shown in the lower right, labeled 1-6. Figure 5. Examples of the wetland screening tool output for selected areas of the VDOT delineation area (B1-B6) compared to NWI mapping for the same scenes (A1-A6). VDOT wetlands, used to train and assess tool output, show the true distribution of wetlands. The location of scenes within the study area is shown in the lower right, labeled 1-6. Figure 6. Level of agreement between NWI maps (A1-A6) and tool output (B1-B6) compared to true distribution of wetlands and non-wetlands, as designated by VDOT delineations. Comparisons are summarized as wetland agreement, where predicted wetlands agree with VDOT wetlands; nonwetland agreement, where predicted non-wetlands agree with VDOT non-wetlands; false positives, where VDOT non-wetlands are predicted to be wetlands; and false negatives, where VDOT wetlands are predicted to be non-wetlands. The location of scenes within the study area is shown in the lower right, labeled 1-6.

Confusion Matrix
The confusion matrix resulting from the classification tool testing (Table 4) shows that the screening tool produced a higher kappa statistic (0.46 vs. 0.34, respectively) and predicted 76.1% of Figure 6. Level of agreement between NWI maps (A1-A6) and tool output (B1-B6) compared to true distribution of wetlands and non-wetlands, as designated by VDOT delineations. Comparisons are summarized as wetland agreement, where predicted wetlands agree with VDOT wetlands; non-wetland agreement, where predicted non-wetlands agree with VDOT non-wetlands; false positives, where VDOT non-wetlands are predicted to be wetlands; and false negatives, where VDOT wetlands are predicted to be non-wetlands. The location of scenes within the study area is shown in the lower right, labeled 1-6.

Confusion Matrix
The confusion matrix resulting from the classification tool testing (Table 4) shows that the screening tool produced a higher kappa statistic (0.46 vs. 0.34, respectively) and predicted 76.1% of the total testing area correctly, whereas the NWI correctly predicted 80.5% of the same area. Focusing more specifically on wetland prediction rates, the screening tool omitted 22.6% of wetlands from the testing dataset in its predictions, compared to 69.3% of the testing wetlands being missed by NWI maps. However, the screening tool incorrectly included 24.3% of testing nonwetland area in wetland predictions, whereas the NWI maps misidentified only 1.3% of the nonwetland area. These statistics are in line with trends observed in Figures 5 and 6, where the wetland screening tool encompassed more of the true wetland area while also over estimating the distribution and extents of VDOT wetlands. This assessment shows that, when used as a preliminary screening tool, it is expected that the method would identify 77.4% of wetlands (one minus the false negative rate) and 24.3% of non-wetland area would be unnecessarily surveyed (per the false positive rate). In contrast, using solely the NWI as a preliminary estimate of wetlands, only 30.7% of wetlands would be identified and 1.3% of the nonwetland area would be unnecessarily surveyed. Given that missing wetland areas can have significant consequences in the LEDPA environmental planning process, these statistics illustrate the potential benefit of using this wetland screening tool, along with targeted wetland identification by experts, compared to relying solely on NWI in the LEDPA process. Table 4. Confusion matrices used to assess the accuracy of screening tool predictions (A) and the NWI raster (B), where testing data classes are represented in columns and both tool and NWI classes are represented in rows. Both accuracy assessments quantify the level of agreement achieved by respective datasets within the limits of the testing dataset.

A Screening Tool Prediction Classes Testing Data Classes
Wetland (km 2 ) Non-Wetland (km 2 ) =

Ecohydrologic Insights into the Method Performance
Understanding the underlying ecohydrologic properties of the study region and how they influence the method's performance is useful for guiding future efforts to improve the method. Scene B6 of Figures 5 and 6 shows a relatively large distribution of false negative predictions. Investigation of the input data in this area shows that wetland omission may be due to the presence of non-wetland characteristics, including low TCT wetness values, scattered distributions of small local depressions and small hydric soil areas, location outside of floodplains, location outside of riparian zones, and mostly developed NLCD classification. It is possible that wetland characteristics are present on a finer scale here, but the current resolution of wetland indicators resulted in overly generalized boundaries between heterogeneous areas. In addition, VDOT wetlands in this area may be the result of data not included in the screening tool like groundwater levels, which are important in this region but difficult to quantify due to limited observational data. In addition, higher resolution topographic inputs that take into account slope, curvature, and flow accumulation have been found to successfully capture wetland hydrology (e.g., [24,[52][53][54][55][56][57][58][59]) and may improve approximation of saturated areas at a finer resolution. Similar studies have shown that finer-scale multispectral imagery would better distinguish detailed boundaries between developed areas and vegetation types (e.g., [16,17,19]), which would be especially important for identifying wetlands close to developed structures.
Although secondary to false negative rate, maintaining a low false positive rate is also important to creating a trustworthy screening tool. Scene B2 of Figures 5 and 6 demonstrates problematic false positive predictions by the screening tool. The input data here show widespread distributions of hydric soils and local depressions, both of which are likely to have contributed to wetland predictions. Disagreement between these predictions and VDOT wetlands could likely be related to overly generalized landscape due to coarse input data, as described above. However, false positive predictions that are adjacent to and surround correct wetland predictions may represent difficulty in assigning definitive extents to seasonally fluctuating wetland boundaries. While the Landsat-derived inputs generally represent the seasonal characteristics inherent to the VDOT delineations, other inputs can be considered time-averaged wetland indicators. Information from these input data likely contribute to screening tool predictions in the correct general area but with imprecise boundaries, as the true wetland boundaries are diffuse and season-dependent.

The Importance of Variable Inputs
Variable importance measures were used to investigate the value of the multisource input dataset. These measures can help target future efforts in refining the datasets most important to the wetland screening methodology. Figure 7 shows the variable importance measures for the proposed eight-input classification (classification 1), as well as measures for execution using the six most important inputs from classification 1 (classification 2) and the four most important inputs from classification 2 (classification 3).
Although secondary to false negative rate, maintaining a low false positive rate is also important to creating a trustworthy screening tool. Scene B2 of Figure 5 and Figure 6 demonstrates problematic false positive predictions by the screening tool. The input data here show widespread distributions of hydric soils and local depressions, both of which are likely to have contributed to wetland predictions. Disagreement between these predictions and VDOT wetlands could likely be related to overly generalized landscape due to coarse input data, as described above. However, false positive predictions that are adjacent to and surround correct wetland predictions may represent difficulty in assigning definitive extents to seasonally fluctuating wetland boundaries. While the Landsat-derived inputs generally represent the seasonal characteristics inherent to the VDOT delineations, other inputs can be considered time-averaged wetland indicators. Information from these input data likely contribute to screening tool predictions in the correct general area but with imprecise boundaries, as the true wetland boundaries are diffuse and season-dependent.

The Importance of Variable Inputs
Variable importance measures were used to investigate the value of the multisource input dataset. These measures can help target future efforts in refining the datasets most important to the wetland screening methodology. Figure 7 shows the variable importance measures for the proposed eight-input classification (classification 1), as well as measures for execution using the six most important inputs from classification 1 (classification 2) and the four most important inputs from classification 2 (classification 3).   Figure 7 shows that TCT brightness and greenness were among the most important input variables in each classification. These indices have also been found to be key input variables in multisource datasets used to identify wetlands in similar studies (e.g., [23,60]). Here, it is expected that TCT brightness and greenness were consistently important wetland indicators as they provide vegetative and optical characteristics that are also seasonally representative of the VDOT delineations. Since these data are derived from Landsat imagery, it is important to select Landsat images that provide ideal conditions for wetland classification. As stated in the introduction, this means spring conditions with high soil moisture and cloud cover are minimal. Without these conditions, it is possible that the method will be unable to accurately identify potential wetland locations. The importance of the hydric soils input was also consistently high, which has been observed in related studies as well (e.g., [24,60]). While it was expected that hydric soils would be successful indicators of wetlands as they are included in the wetland criteria, this trend shows that SSURGO hydric soil information was useful at our target scale despite the relatively coarse resolution of these data. The low importance of the local depressions input was unexpected since this layer was intended to indicate areas where water is likely to pool. It is likely that the elevation data used was too coarse to model the low-relief area and that additional topographic metrics would more robustly describe hydrologic drivers of wetland formation.
Additional accuracy assessments were performed for classifications 2 and 3 and are shown with the classification 1 assessment in Table 5. The accuracy assessments show that overall accuracy decreases consistently as input variables were removed from the workflow. In addition, false positive rates increase and false negative rates are variable. The overall accuracy trend shows that the screening tool is overall better able to distinguish between wetlands and non-wetlands given information from the eight originally proposed input variables. The increasing false positive rates suggest that local depressions, floodplains, riparian areas, and TCT brightness provided important information that contribute to correct non-wetland designations. Despite the false negative rate decreasing by a small margin between classifications 1 and 2, the overall higher false negative rate in classification 3 demonstrates the ability of the larger input set to capture a more robust set of wetland characteristics. It is important to note, however, that accuracy rates shown in Table 5 varied only slightly, which suggests that some of the input data provide similar spatial information, and therefore there is a relatively small cost to the accuracy when one of these layers is removed.

Potential for Additional Data and Tool Improvements
While there are several ways to advance the wetland screening method presented here, improving the quality of the input data is an important next step. In particular, using higher resolution topographic data would likely improve results. To maintain the goal of using data generally available at a national-scale, incorporating higher resolution elevation data would be the most viable next step as Light Detection and Ranging (LiDAR) elevation data is quickly growing in availability [53]. Furthermore, the ability to model small topographic changes through data products easily derived from digital elevation models (DEMs) would be particularly useful in mapping saturated areas in the low-relief terrain studied here. If high resolution multispectral imagery becomes widely available, this would also be a valuable addition to the improved wetland screening tool and may present potential for distinguishing between different classes of wetlands. Moreover, radar data may contribute additional vegetative information that is helpful for wetland mapping, as demonstrated by researchers (e.g., [23,26,61,62]). Also regarding the quality of input data, the riparian zone can be estimated in a more sophisticated way. The application of a standard riparian distance to all waterbodies could be improved by instead using variable riparian buffer distances based on the size of waterbodies and bank geomorphology, as suggested by other studies [15,24].
Additional field-verified wetland maps used to train the model and quantify its error could be obtained for other areas. Including more training data from other study areas would likely improve the robustness of wetland and non-wetland signatures detected by the random forest model. Further analysis would likely benefit the procedure used to determine the optimal train-test split of these delineated wetland areas to determine when diminishing returns of accuracy begin as more of the area is used for training. Given the unequal distribution of wetland to nonwetland area in this study region, and likely other regions as well, further work to optimize the sampling scheme between wetland and nonwetland classes of the training data, which has been shown to significantly impact wetland classifications [22], is another area for improvement.
Lastly, reconfiguring the workflow to execute using open source Python libraries would offer several advantages. Open source geoprocessing libraries such as GRASS GIS [63], GDAL [64], and TauDEM [65] would allow users without an ArcGIS license to execute necessary processes. Additionally, Scikit-Learn [66] implementation of random forest would offer more flexibility than the ArcGIS implementation in terms of random forest parameters. Following this, analyses for random forest parameter tuning should be performed to test if moving from the standard configuration for parameters improves classification accuracy.

Conclusions
This study presents a methodology designed to screen for wetlands over a large geographic region using nationally available geospatial data as input and random forest classification. The methodology was motivated by the desire to streamline environmental permitting for transportation corridor projects over large regions. By using a tool to screen for potential wetland areas, field surveying efforts could be focused to areas that are likely to contain wetlands. The methodology was implemented as an automated workflow in a commercially available geographic information system (GIS) software commonly used by DOTs. The tool was applied to identify potential wetland locations for a region in the coastal plain of Virginia, USA. The preliminary implementation of this workflow was evaluated against professionally conducted field surveys and results were compared to the commonly used NWI dataset as a benchmark for accuracy.
Results showed that, when compared to the NWI, the wetland screening methodology produced a significantly lower false negative rate (22.6% vs. 69.3%) and a higher kappa statistic (0.46 vs. 0.34). From this, we conclude that the methodology was able to capture many wetlands missed by NWI. However, this improvement in false negative predictions did result in a higher false positive rate (24.3% vs. 1.3%) and, because the study area has more non-wetland area, a slightly lower overall accuracy (76.1% vs. 80.5%). From this, we conclude that, while the method identifies significantly more true wetlands than the NWI, it comes at the cost of a slight reduction in overall prediction accuracy. This was largely by design, however, as the method purposely avoids false positives because such errors would result in missing wetlands and be costly in environmental planning. False negative errors are less costly because they could be field verified by targeted, on-the-ground surveys by experts to fine tune the wetland delineation. With additional wetland areas mapped and verified by experts, the method could be tested for other regions given that, by design, it only relies on nationally available input datasets.
While successful as a screening tool, the ultimate goal should be to achieve the highest possible overall classification accuracy at a spatial scale relevant for environmental planning purposes. Doing so would move the approach from being a wetland screening tool to a wetland mapping tool, opening up additional potential uses. For example, because the approach is largely automated now, future work could also investigate the potential for deriving time varying wetland maps using Landsat imagery and a dynamic input to capture change in wetland patterns across regions. Making this transition will likely require much higher resolution data and more data for training classification algorithms.
The wetland detection rate produced by the current screening algorithm, which again only makes use of nationally available geospatial data in order to be widely applicable, is encouraging for creating approaches able to identify wetlands at a high resolution. Future work should investigate the role of higher resolution input data, like LiDAR, and alternative parameterizations of the classification algorithm to improve wetlands predictions.