Optimized Stratification for Mapping Riparian Vegetation in Arid and Semiarid Environments

This paper describes a method of mapping riparian vegetation in semi-arid to arid environments using the Landsat normalized difference vegetation index (NDVI). The method successfully identified a range of riparian community types across the entire state of Nevada, USA, which spans 7 degrees of latitude and almost 4000 m of elevation. The landscape was stratified into units of similar elevation and solar exposure, and riparian areas were identified as having anomalously high NDVI within a local neighborhood. Thousands of calibration points were used in a simplex optimization to select the spatial neighborhood, the elevation and insolation strata, the minimum NDVI to be considered as potentially riparian, the number of standard deviations from the mean for an anomaly to be classified as riparian, and a limit on upslope position. Mapping of subpixel riparian corridors was improved by applying a directional high-pass filter to the NDVI data. Irrigated areas in agricultural and urban areas were removed based on land ownership maps and manual editing. The final map was tested with 400 independent test points: producer’s accuracy was 84.6% and user’s accuracy was 93.5%. This method should be broadly useful for mapping riparian features across large and complex regions with arid to semi-arid environments.


Introduction
Riparian ecosystems provide critical resources for arid to semi-arid regions [1], enhancing local plant and animal diversity [2] and supporting migratory bird populations [3]. maps derived from satellite imagery are an important tool for monitoring and managing these important ecosystems [4]. This study demonstrates a novel approach for mapping riparian habitats throughout the state of Nevada, USA. Nevada combines a number of challenging characteristics for riparian mapping, including low annual precipitation, mountainous terrain, endorheic drainage systems, a large range of latitudes and elevations, and a wide range of riparian types including wetlands, springs, gallery forests, wet meadows, and mesquite trees in washes of the mojave Desert. In Nevada's mountains, riparian vegetation may be confined to narrow, discontinuous patches within deep canyons where the amount of surface water may change dramatically over short distances due to geological and geomorphic constraints. In contrast, outflow from the mountains may be distributed in braided channels over alluvial fans and then spread widely over broad, flat playas.
Many riparian mapping studies have been performed by pragmatically generating a spatial buffer around streamlines in digital hydrographic datasets [5,6], or by more realistically applying a flooding algorithm to a digital elevation model (DEM) [7,8]. Salo et al. [9] provide a useful comparison of a number of these techniques. The goal of either approach is to constrain an analysis of complex terrain to just that part of the terrain where riparian features are expected, but both options are problematic for application across a large, dry, mountainous region. The level of detail and accuracy of available hydrographic datasets is quite variable in remote areas of Nevada, the width of riparian corridors changes dramatically over short distances, and small but important riparian patches exist that are disconnected from mapped flow lines. As such, simple distance buffers that would be reasonable for a small study area are not effective for this large region of complex terrain. For flooding algorithms, riparian corridors are often narrow relative to the resolution of the DEM. Additionally, for broad, smooth playas there would be no topographic constraint for a flooding algorithm. Finally, in the arid to semi-arid climate of Nevada, riparian vegetation may be distributed discontinuously within an arbitrary buffer zone or digitally flooded model.
The challenges to direct detection of riparian vegetation in satellite imagery include variable species composition that creates complex spectral signatures, ambiguity as to the abundance of vegetation or its relative vigor, short-term changes in vigor due to individual precipitation events, and longer term changes arising from disturbance events like wildfires. Some prior work with Landsat mapping of riparian zones has not provided encouraging results [5], though many of these studies have been based on distinguishing different plant species or communities within a riparian zone rather than presence/absence. Even when just mapping the presence of riparian vegetation with Landsat and geospatial data, low map accuracies have been reported [10]. Alaibakhsh et al. [11] used a principal components analysis of multi-date Landsat NDVI data to map riparian vegetation in the Pilbara region of Australia, and they reported 91% accuracy based on a small test area located within the overall region. No information on rates of omissions and commissions was provided. makkeasorn et al. [12] described the use of geospatial datasets, RADARSAT imagery, and Landsat-based vegetation indices for mapping riparian vegetation in Texas. Using 47 riparian ground sites and 23 non-riparian sites, they reported an accuracy of 90%. Given the strong effects of topography on radar backscatter, it is unlikely that such a method would be applicable in the complex, mountainous terrain of Nevada.
In another approach that relied on direct moisture detection, Barron et al. [13], used the normalized difference wetness index (NDWI) from Landsat as part of a strategy to map ground water dependent ecosystems. Initial examination of NDWI in Nevada proved problematic, as the absorption feature for foliar moisture is weaker than the red edge that is the basis of most vegetation indices. Because of this, background soil color often overwhelmed the moisture signal in sparsely vegetated regions. Baker et al. [14] used a combination of Landsat, elevation, and soils data to map riparian and wetland areas in the lower reaches of the Gallatin watershed in montana, USA. That study reported an accuracy of 86%, and the accuracy assessment was one of the more rigorous examples encountered. However, the situation for mapping the entire state of Nevada would be much more complex than for the broad scale hydrologic patterns in the lower reaches of a single drainage in a more mesic climate.
The mapping approach presented here takes advantage of the water-limited condition of Nevada's summer to identify riparian vegetation as locations where the normalized difference vegetation index (NDVI, [15]) for a particular elevation and level of sun exposure remains significantly more vigorous over a number of years than surrounding areas with the same elevation and exposure. multiple studies have demonstrated the ability of Landsat NDVI to discriminate relative levels of water availability within a given plant community [16][17][18][19]. Focusing on the detection of localized NDVI anomalies provided a single approach that adapted to various changes in vegetation and environmental settings across a large and complex terrain.

Study Area
The state of Nevada covers an area of more than 28.6 million hectares, spanning 35 • N to 42 • N and 114 • W to 120 • W, with eight Koppen climate zones ranging from hot desert to dry-summer subarctic. It is located predominantly within the Great Basin region, the largest contiguous area of endorheic basins in North America. The southernmost portion of the state is in the mojave Desert. Nevada's terrain is largely associated with the Basin and Range physiographic region, a vast area of horst and graben faulting that created parallel ranges of north-south oriented mountain ranges separated by arid valleys. Prevailing winds are intercepted by the Sierra Nevada mountains to the west, creating a rain shadow that makes Nevada the driest state in the country. Annual precipitation ranges from 106 mm at Las Vegas in the south to 249 mm at Elko in the northeast. Greater precipitation occurs at higher altitudes in the many mountain ranges. The majority of precipitation in the western and south-central portions of the state falls in the winter, the central and northeastern portions receive more in the spring, and some areas of the east receive the most from thunderstorms in late summer.

Measurements
This riparian map was generated using data from Landsat's Thematic mapper (TM) sensor from a number of years that had above average precipitation : 1986, 1988, 1993, 1997, 2005, and 2009. This selection of six wet years ensured that the riparian greenness signal would be strong and that disturbance events that affected certain years would not interfere with detection. A persistent drought from 2013 to 2016 followed by extensive flooding in 2017 kept Landsat 8 imagery from being used. Landsat images were acquired in late summer, between July 15 and September 30, in order to emphasize the difference in greenness between riparian areas and the surrounding dry environments, to reduce snow cover at higher elevations and possible flooding at lower elevations, and to take advantage of a period of slow phenological change. The TM Collection 1 Tier 1 data were acquired using Google Earth Engine, NDVI was calculated for each scene, and the median of NDVI values for the entire late summer season was used to represent each year. maximum value compositing has often been used to create cloud-free seasonal composites, however, in this case, the maximum value would preferentially select periods following infrequent summer precipitation events which may obscure the greenness signal of the riparian systems. In addition, maximum NDVI compositing preserves Landsat data artifacts such as sun glint anomaly and data loss, whereas the median removes these infrequent large errors [20]. Nevada has a very high proportion of cloud-free summer days, and the median of this 11 weeks period was very effective at cloud removal.
Given the mountainous terrain and limited precipitation in Nevada, many riparian areas are narrower than the spatial resolution of the Landsat imagery. While the contrast between target and background for these wet corridors is often great enough to make them visible, it does substantially diminish the greenness signal. Spatial filtering was used to boost NDVI values where the local image texture indicated a linear greenness feature. This was accomplished by subtracting a 5 × 3 directional Gaussian low-pass filter (Table 1) from the original NDVI data. The shorter dimension (3) of the filter kernel allowed the method to accommodate some meander in the linear features. In order to enhance linear features in multiple directions, the spatial filter was applied with the NDVI data rotated by 0, 45, 90, and 135 degrees, and the maximum enhanced value of any rotation was used to represent each pixel. While this enhancement increased the signal for narrow riparian corridors, it suppressed the signal of broader riparian areas, so the mapping method was applied separately to original and enhanced NDVI datasets, and the two resulting maps were merged afterwards. The method also used the 1 arc second (30 m) DEM of the USGS National Elevation Dataset (NED). Total solar radiation during the seasonal timeframe of Landsat acquisitions was calculated from the DEM as watt hours per square meter using a program in the ArcMap software package (ESRI, version 10.3) with the parameters provided in Table 2. The mapping method also used a calculation from NED that identified the difference of each pixel from a 1 km × 1 km Gaussian low-pass spatial filter. This difference from the regional mean elevation helped to identify whether pixels were in a more upslope or downslope position, and this variable will be referred to as ∆DEM1k. The accuracy of the USGS National Hydrology Dataset for Nevada was insufficient to support this analysis, however, a value-added modification of that product called NHDPlus [21] was used to identify some perennial streams. The NHDPlus product omits many smaller perennial features.
The online basemap imagery that is available through the ArcMap software package was used to help identify land cover types. ArcMap identifies the source of this digital imagery as a compilation from ESRI, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

Data Analysis
A total of 21 Landsat scenes were required to completely cover the state of Nevada. In order to manage the volume of computer memory and to allow for regional variability in the calibration of the mapping parameters, all datasets for the state were divided into eight tiles ( Figure 1) for separate processing. These processing tiles had approximately 25 km of overlap with each other in order to ensure consistency along edges in the final composited product for the state.  The accuracy of the USGS National Hydrology Dataset for Nevada was insufficient to support this analysis, however, a value-added modification of that product called NHDPlus [21] was used to identify some perennial streams. The NHDPlus product omits many smaller perennial features.
The online basemap imagery that is available through the ArcMap software package was used to help identify land cover types. ArcMap identifies the source of this digital imagery as a compilation from ESRI, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

Data Analysis
A total of 21 Landsat scenes were required to completely cover the state of Nevada. In order to manage the volume of computer memory and to allow for regional variability in the calibration of the mapping parameters, all datasets for the state were divided into eight tiles ( Figure 1) for separate processing. These processing tiles had approximately 25 km of overlap with each other in order to ensure consistency along edges in the final composited product for the state.  The mapping method was designed to identify anomalously abundant or vigorous vegetation for a given elevation and amount of sun exposure as an indicator of surface or near-surface water availability in a semi-arid or arid environment. By focusing on anomalous greenness for a specific environment, this approach avoided the complications of species-specific spectral reflectance [22], similarities in NDVI between upland forest and lowland riparian areas, confusion of surface water availability versus the reduced moisture stress of north-facing slopes, and differences in the overall amount of vegetation cover between more or less arid parts of the state. The area within each processing tile was subdivided into a grid, and anomalies were calculated separately for each grid cell. This localized analysis reduced the variability from large scale trends in rainfall, temperature, land management, and topographic rain shadows that would confound the identification of a localized anomaly. Within each grid cell, the landscape was divided into strata of elevation and solar insolation, and then deviations from the mean of each stratum were calculated. Twenty five percent overlap was added around each grid cell extent and between each stratum and its next higher and lower strata, in order to reduce discrete jumps in the calculation of the anomalies. To improve the determination of an anomaly, a map was created to mask out possible riparian or irrigated areas from calculations of the mean and standard deviation within each stratum of elevation and insolation. This mask was created by displaying NDVI and NED datasets and quickly and coarsely painting over features that appeared to possibly contain riparian or irrigated vegetation, purposefully being inclusive so that no areas would be missed. In order to compensate for disturbance events and to reduce the effect of data errors in a given year, the method was applied to multiple years of NDVI data, and pixels that were only classified as riparian in one year were removed.
The steps of the mapping method are outlined below, with the adjustable calibration parameters identified in brackets.

1.
Create the coarse mask of potential riparian and irrigated urban or agricultural areas.

2.
Refine the mask (Step 1) by unmasking any pixels that are less than [a specified NDVI value].

3.
Select one year of NDVI data for analysis.

4.
Define a [grid size] to set the scale of the neighborhood used in analyzing anomalies.
Within each grid cell, stratify the solar exposure map by a [radiation increment]. 7.
Starting at the lowest elevation stratum, and proceeding from the lowest to highest insolation strata, merge strata until a minimum number (20) of unmasked pixels are collected or no more are available. 8.
Calculate the mean and standard deviation for NDVI of the collected unmasked pixels. 9.
Calculate the NDVI anomaly of each masked pixel in the accumulated strata as the number of standard deviations from the mean of unmasked pixels. 10. Apply a threshold [number of positive standard deviations] to identify masked pixels as riparian or not. 11. Remove pixels that were classified as riparian if they exceed [a threshold value of ∆DEM1k]. 12. Repeat Steps 7-11, completing the insolation strata within an elevation stratum and then moving to the next higher stratum of elevation. 13. Repeat Steps 5-12 with remaining grid cells. 14. Repeat Steps 3-13 with remaining years of NDVI data. 15. Remove riparian pixels that only exist in one year of observations.

Parameter Optimization
A simplex optimization method was used to calibrate the following algorithm parameters:

1.
Grid size: the local neighborhood for which strata are calculated, 2.
minimum NDVI: the minimum NDVI to be considered potentially riparian,
Radiation increment: the increment determining insolation strata, and 5.
∆DEM1k threshold: maximum height above regional mean elevation that is allowed to be riparian.
Parameter optimization used the AMOEBA function in the IDL software package (Version 8.6, Harris Geospatial). The objective function for the simplex was based on the errors of omission and commission from photo-interpreted calibration points distributed throughout the state. Separate calibration points were selected for riparian areas in the original NDVI data versus the directionally filtered data. For the original data, calibration points were placed where the spatial extent of riparian vegetation was a large portion of the pixel. With the filtered data, calibration points focused on narrow, subpixel features or where a stream went through an area of generally dense vegetation. Initially, 1200 points were randomly located in order to provide broad coverage of land cover types. Additional sample locations were manually selected to focus on different types of riparian environments throughout the state and the nearby land covers that would likely cause confusions. Irrigated agricultural and urban areas were not included as calibration points, as these anomalously green areas would be removed later using a separate map of irrigated lands.
The simplex calibration was run at least four times for each processing tile using different initial values for the parameters for the original and filtered NDVI data. The highest scoring combinations of parameters for each calibration run were examined. While the highest scoring combination was usually selected, in some cases an alternative was selected that provided a large increase in either user's or producer's accuracy, accompanied by little or no penalty to the other measure. The tile was classified using the selected parameters, and additional photo-interpreted calibration points were added in locations where obvious confusions were observed. The allowed range of parameters in subsequent iterations was narrowed based on the observed convergence in parameter values. To speed up the process, the initial optimizations used three years of Landsat data: 1986, 1997, and 2009. The final optimization for each tile used all 6 years of Landsat data. By the final calibration, there was a combined total of 9317 calibration points for the two datasets.
The objective function for the simplex optimization was: where: P = producer's accuracy for riparian class; U = user's accuracy for riparian class; max = maximum of arguments in parenthesis.
This equation had the effect of prioritizing errors of commission over omission until the user's accuracy exceeded 75%, and then continuing the search for higher producer's and user's accuracies so long as the user's accuracy did not decline. This made use of the knowledge that NDVI and filtered NDVI sources were to be merged in the final product, so omissions from one source might be ameliorated by detection in another source. The value of 0.75 set a floor on errors of commission in order to limit the accumulation of error across these sources. A higher baseline value was not used because training samples were added disproportionally to areas with misclassifications in earlier iterations of the optimization. Simpler objective functions, such as overall accuracy, the mean of P and U, or the minimum of P and U, would be more sensitive to the number of riparian versus non-riparian samples and would not consider the potential for accumulated commission errors in the final product. While the optimization did provide an objective, unsupervised method for finding a good combination of parameters, the AMOEBA function would not cleanly converge to a single, optimal solution because of the discontinuous response of the error statistics to changes in parameters.

Post-Processing
The number and position of perennial stream features in the NHDPlus data product differed from those of the original USGS NHD that it was derived from, likely indicating that effort was invested in improving the accuracy of those features. NHDPlus was used to improve mapping in areas of dense vegetation, where the contrast of the riparian signal would be relatively low compared to surrounding forest. However, the NHDPlus product represents stream features at 1:100,000 scale, so the position of the stream lines deviated from the channel position observed in the Landsat and DEM data. Stream lines were typically within 3 pixels of corresponding features in the imagery. Using a moving 5 × 5 window, the highest five NDVI values that were within 3 pixels of a perennial stream line were labeled as riparian vegetation so long as they had a NDVI value greater than 0.2, a threshold that was used for riparian mapping in References [23] and [24]. The numerous stream lines that were labeled as intermittent in NHDPlus were not used, since elevated NDVI values around these features often represented unrelated changes in the amount of vegetation cover.
Results from the analysis of NDVI, the spatially filtered NDVI, and the perennial stream line buffer were merged into one map and were then masked by a map of irrigated lands in urban and agricultural areas. This agriculture/urban mask was created using existing land cover and ownership datasets, and was extensively edited based on inspection of the digital aerial photography in Arcmap and satellite imagery. In the water-limited environment of Nevada there is substantial use of riparian zones for agriculture and ranching. In cases of less intensive agricultural manipulation, these land uses of are capable of providing ecosystem services, and they also respond to the same stresses as their interconnected natural riparian systems. As such, agricultural land use in riparian zones that appeared to be based on available surface water, and for which NDVI values were comparable to neighboring unmodified riparian areas, was not masked out.
Finally, the riparian map product was superimposed on the DEM to identify anomalously green pixels that were located in unrealistic topographic positions for riparian vegetation. Pixels that were on ridgelines or short slopes without visible channel development were manually erased from the map. This editing process included inspection of the imagery and NDHPlus perennial stream lines.

Accuracy Assessment
Because of the limited spatial extent for riparian areas in Nevada, equal numbers of validation samples were selected for riparian and non-riparian areas instead of using area-weighted sampling. Additionally, since a large proportion of the state has very low amounts of vegetation cover that would generally not be confused with riparian areas, sample locations for accuracy assessment were selected from within the coarse mask of potential riparian areas that was created for the classifier. The combination of these two sampling constraints provides an accuracy assessment that is more relevant to applications in and around riparian zones, and an underestimation of the overall accuracy of the map on a statewide basis.
Within the mask of potential riparian areas, 200 pixels were randomly selected from areas that were classified as riparian and 200 that were classified as non-riparian. Given the map accuracy expectation of 85% set out in Reference [25], 400 points provides a 95% confidence interval of 81.2% to 88.2% [26]. Each validation point was examined and classified by two photo interpreters who consulted online digital photography, the DEM, Landsat imagery, and streamlines from NHDPlus. Indications of riparian vegetation included channel position, visible water, patterns of soil moisture, changes in plant species composition, and abrupt changes in plant vigor and cover density. A change in plant density alone was only used to label a point as riparian if the change was very localized in a channel, and there were no similar levels of plant density outside channels in the surrounding region. There did not have to be direct evidence of water present in the photography, and intermittent streams and springs were classified as riparian if they met the aforementioned criteria. man-made features such as impoundments on grazing allotments were labelled riparian, so long as they appeared capable of providing ecological services. Relict features, such as cut-off stream oxbows, were counted as riparian if the vegetation continued to be distinct from the surrounding environment.
Because of the mountainous terrain and differences in image scale, there was some spatial misregistration of the online digital photography, DEM, and Landsat data. In these cases, the interpreters would use matching features to reconcile the location of validation points. When the location could not be confidently reconciled, or the point was directly on a riparian boundary, the point was replaced with a new random location.

Results
The final calibrated parameters of the mapping method for each processing tile are presented in Table 3 for the original NDVI data and Table 4 for the directionally filtered data. When using three years of data, a calibration run for the larger tiles took approximately 4 h on a PC computer with a 3.4 GHz i7-4770 CPU. Typically, the optimization program could be interrupted within 2-3 h as the parameters continued circling a small range of values. Given the time to generate the ancillary datasets and four iterations of the calibration routine, approximately 7 h were required per million hectares. However, that estimate does include the efficiencies of scale associated with this very large study area.  Figure 2 presents results of the analysis for an area where streams drain from a forested mountain (center) to surrounding dry shrubland. The original NDVI data in Figure 2a portrays the challenge of separating the overlapping NDVI distributions between riparian and other types of vegetation on the topographic gradient. Figure 2b,c shows how the mapping algorithm suppresses the typical vegetation signal for a characteristic position on the landscape and isolates the anomalously green areas along drainages. Figure 2d shows the coarseness of the simple mask of potential riparian areas (white) that was used to remove most riparian pixels from the calculation of statistics for anomalies. Figure 2d also shows that pixels classified based on proximity to perennial streams in the NHDPlus dataset (black) miss a number of drainage areas where riparian vegetation is present. Using the original NDVI data, the algorithm in Figure 2e broadened the area classified using NHDPlus and identified a number of other drainages with riparian vegetation. The filtered NDVI data added some smaller drainages in Figure 2f and provided better continuity between the patches that were resolved in the original NDVI data. The broad area of anomalously high NDVI values in the upper left corner of Figure 2a  Statewide, 77.4% of pixels that were classified as riparian were identified by the original NDVI data, 20.1% were additional pixels identified by the directionally filtered data, and the remaining 2.5% were pixels were labeled with the additional use of NHDPlus. The final manual editing based on inspection of the results overlain on the DEM and Landsat imagery eliminated 0.14% of the pixels that were classified as riparian. The total area classified as riparian was 1.0% of the state of Nevada.
Results of the accuracy assessment are presented in Table 5. The overall percentage of correctly classified pixels was 88.3%, which, based on the aforementioned confidence limit, indicates that the map is confidently better than 85% accurate. User's accuracy for the riparian class was higher than producer's accuracy, indicating that the total riparian area was underpredicted somewhat, but that there were relatively few errors of misidentifying riparian areas.  Statewide, 77.4% of pixels that were classified as riparian were identified by the original NDVI data, 20.1% were additional pixels identified by the directionally filtered data, and the remaining 2.5% were pixels were labeled with the additional use of NHDPlus. The final manual editing based on inspection of the results overlain on the DEM and Landsat imagery eliminated 0.14% of the pixels that were classified as riparian. The total area classified as riparian was 1.0% of the state of Nevada.
Results of the accuracy assessment are presented in Table 5. The overall percentage of correctly classified pixels was 88.3%, which, based on the aforementioned confidence limit, indicates that the map is confidently better than 85% accurate. User's accuracy for the riparian class was higher than producer's accuracy, indicating that the total riparian area was underpredicted somewhat, but that there were relatively few errors of misidentifying riparian areas. The spatial distribution of erroneous validation points maintained a relatively uniform frequency across the state (Figure 3), and no significant clustering of error was detected by moran's I (P = 0.84) or the Getis-Ord General G test using inverse distance (P = 0.16) or inverse distance squared (P = 0.37).
The spatial distribution of erroneous validation points maintained a relatively uniform frequency across the state (Figure 3), and no significant clustering of error was detected by Moran's I (P = 0.84) or the Getis-Ord General G test using inverse distance (P = 0.16) or inverse distance squared (P = 0.37).

Discussion
While this mapping approach was based on how vegetation responds to the local energy and water balance, there are limits to the ecological conclusions that can be drawn from the calibration results. The calibration responds to a number of factors that control the distribution of riparian vegetation, and the interaction between parameters may be complex within an area as large as these processing tiles. For example, a particular elevation increment might simultaneously isolate a hydrologic constraint from a geologic formation and the lowest elevation of a non-riparian montane tree species. Acknowledging this caveat, the average elevation increment across Tables 3 and 4 was 23 m, and, including the 25% overlap with adjoining strata, that would correspond to a dry adiabatic lapse of 0.34 °C. It might be inferred that vegetation amount and vigor did not change greatly within that range of temperature change.
It is difficult to directly interpret the insolation values in Tables 3 and 4, since they relate to a seasonal total. However, if these values are compared to the total range of insolation values in a processing tile, they represent an average of 28 strata per tile. That compares to an average of 23 strata over the range of elevations in a tile. This suggests that relatively small changes in insolation were important for isolating greenness anomalies.

Discussion
While this mapping approach was based on how vegetation responds to the local energy and water balance, there are limits to the ecological conclusions that can be drawn from the calibration results. The calibration responds to a number of factors that control the distribution of riparian vegetation, and the interaction between parameters may be complex within an area as large as these processing tiles. For example, a particular elevation increment might simultaneously isolate a hydrologic constraint from a geologic formation and the lowest elevation of a non-riparian montane tree species. Acknowledging this caveat, the average elevation increment across Tables 3 and 4 was 23 m, and, including the 25% overlap with adjoining strata, that would correspond to a dry adiabatic lapse of 0.34 • C. It might be inferred that vegetation amount and vigor did not change greatly within that range of temperature change.
It is difficult to directly interpret the insolation values in Tables 3 and 4, since they relate to a seasonal total. However, if these values are compared to the total range of insolation values in a processing tile, they represent an average of 28 strata per tile. That compares to an average of 23 strata over the range of elevations in a tile. This suggests that relatively small changes in insolation were important for isolating greenness anomalies.
The average calibrated NDVI threshold from Tables 2 and 3 was 0.194, which is very close to the 0.2 minimum NDVI that has been used in other riparian studies [23,24]. The average threshold for number of standard deviations from the mean of unmasked pixels was 5.1. While this represents a large anomaly, its large magnitude is due to the use of the mask that removed potential riparian zones from the calculation of the standard deviation within a spatial grid.
In Table 3, the average ∆DEM1k was 18 m, indicating a safe height for removing anomalously green stands of upslope vegetation, such as a copse of trees in a sheltered pocket above the local tree line. In Table 4, the ∆DEM1k values for tiles 3 to 8 were much lower, largely reflecting the degree to which calibration points for the directionally filtered data focused on classifying very narrow riparian features that were located deep in canyons. For tile 3, a low threshold value for standard deviations seemed to balance an extreme value of ∆DEM1k. The larger ∆DEM1k values for tiles 1 and 2 in Table 4 reflect the more subdued relief of the mojave Desert in southern Nevada.
For the original NDVI data, tiles 7 and 8 had higher calibration accuracies than other tiles. This is likely due to broader riparian corridors that are associated with higher annual rainfall totals in the northeastern part of the state. For the directionally filtered NDVI data, the calibration producer accuracies were generally low. This occurred because the high-pass filter created a noisier product, and there was a much higher priority placed on errors of commission in the objective function. Since the final map would be a composite of multiple sources, the higher omission rate that was required to keep noise out of the final map was somewhat less of a concern.
The optimization method provided a structured way to recalculate input variables in a way that could not be replicated with a single statistical analysis of per-pixel environmental variables. The optimization method was slow, but it proceeded in an unsupervised manner and multiple tiles could be run in parallel. An alternate version of this program was tested in which the anomalies were calculated based on a neighborhood and range of elevation and radiation around each pixel, however, the computational requirements were prohibitive; computation by fixed strata was required to make the method effective. While the long compute time was problematic, the results demonstrated in Figure 2b,c were compelling.
The use of perennial streams from the NHDPlus product helped to define riparian paths through dense forest, but made a relatively small contribution at the statewide level. much of the NHDPlus component in this map already was directly detected in the NDVI data, and streams that are mapped as intermittent comprise a large fraction of Nevada's riparian vegetation.
A common concern when using vegetation indices in sparsely vegetated regions is the effect of background soil color. Despite this, NDVI was capable of producing consistent results across the state, because greenness anomalies were calculated relative to a local neighborhood where the background soil reflectance would be relatively constant. Additionally, the riparian areas that were being detected would tend to have a higher percentage of vegetation cover.
The map accuracy calculated from independent validation points was better than the user and producer accuracies for the calibration points. While this situation would typically be suspect, it occurred here because after each iteration of the parameter optimization, a disproportionate number of new calibration points were added in locations where the algorithm was having trouble. Additionally, the accuracy of the final product benefitted from the merging of maps derived from the three different datasets and some manual editing.
The LANDFIRE map [27] produced by the U.S. Geological Survey that was derived from multispectral classification is often used for broad scale vegetation mapping in the United States, including studies of riparian areas [28,29]. An accuracy assessment is available for the LANDFIRE product, organized by geographic regions called super zones [30]. In that assessment, LANDFIRE classes in the Great Basin zone that were labeled as riparian by the Ecological Systems attribute (EVT_Ecol_Sys Codes: 2012, 2154, 2155, 2159, 2160, 2164) had a combined producer's accuracy of 37.5% (16 reference points) and user's accuracy of 23.1% (26 classified points). For riparian classes that were defined as riparian using the Society for Range management (SRM) Type (EVT_SAF_SRM_Type Codes: 1235, 2203, 2418, 2422), producer's accuracy for the combined riparian classes was 72.7% (11 reference points) and user's accuracy was 30.8% (26 classified points). map accuracy using the method of stratified NDVI anomalies that is presented here is substantially better, but that is to be expected given the narrower focus of this study relative to a national scale map that makes numerous class distinctions. User and producer accuracies for the Nevada riparian map match the 85% classification accuracy standard that has traditionally been required for maps that were manually generated from aerial photography by photo-interpreters [25].
This mapping method proved quite flexible in its ability to adapt to changing conditions across more than 286,000 km 2 of rangelands, forested mountains, and desert. Among examples of riparian mapping that were found in the literature, the only study that approached the spatial extent of this effort was Alaibakhsh et al. [31] which mapped groundwater dependent ecosystems across a large region in Australia. Unfortunately, that study was not able to provide a quantitative accuracy assessment for comparison. While computationally intensive, the method presented here may prove useful in other semi-arid environments around the globe.

Conclusions
This mapping method provided a consistent approach for detecting riparian areas over a very large region with varying climate regimes, plant community types, and amounts of vegetative cover. By focusing on the relative vigor of vegetation rather than species composition, the method avoided many challenges of traditional multispectral classification. The novel optimization approach provided an objective, unsupervised method for finding the combination of neighborhood and environmental strata that would best isolate anomalous vegetative vigor. Since the method relies on contrast in water availability between riparian and adjacent environments, the method will work best in arid and semi-arid climates that have a distinct dry season.
Directional high-pass filtering was beneficial for enhancing sub-pixel scale riparian features. By optimizing the method's parameters to ensure a high user's accuracy, the filtered NDVI data was able to add 20% of the total mapped area of riparian vegetation. The method did rely on the manual creation of two maps for identifying potential riparian areas and distinguishing irrigation from natural water availability. While the mask of potential riparian areas was generated manually, there is potential that the types of flooding algorithms that have been used with digital elevation data in other riparian studies could be adapted to automate some of that process. However, there are cases where that would not work, such as when a stream courses down the ridge of an alluvial fan or spreads over just a limited portion of a broad, flat playa.
Going forward, the relatively new Sentinel-2 satellite system with its 10 m resolution in the visible and near infrared bands would be a much better candidate for this type of analysis. Sufficient temporal coverage by Sentinel-2 has now been collected so that longer term ecological states can be discriminated from shorter term changes due to drought, fire, and other disturbances.
An improved knowledge of the distribution of riparian areas will allow us to resolve a clearer history of the relative resilience of different riparian areas to drought, fire, and land management practices, as observed in more than three decades of historical Landsat imagery, thereby allowing us to better protect these resources in the future.