UAV and Satellite-Based Sensing to Map Ecological States at the Landscape Scale

Mapping ecological states in semi-arid rangelands is crucial for effective land management and conservation efforts because it identifies difference in the ecological conditions across a landscape. This study presents an innovative approach for mapping two ecological states, Large Shrub Grass (LSG) and Large Shrub Eroded (LSE), within the Sandy Loam Upland and Deep (SLUD) ecological sites using a combination of drone and satellite data. The methodology leverages the Largest Patch Index (LPI) as a proxy metric to estimate eroded areas and classify ecological states. The integration of unmanned aerial vehicle (UAV) data with satellite-based remote sensing provides a scalable approach that can benefit various stakeholders involved in rangeland management. The study demonstrates the potential of this methodology by generating spatial layers at the landscape scale to inform on the state of rangeland ecosystems. The workflow showcases the power of remote sensing technology to map ecological states and addresses limitations in spatial coverage by integrating UAV and satellite data. By utilizing the bare ground LPI metric, which indicates the connectedness of bare ground, the methodology enables the classification of ecological states at a regional scale. This cost-effective approach potentially offers a standardized and reproducible method applicable across different sites and regions. The accuracy of the classification process is evaluated by comparing the results to ground-based polygons, dirt roads, and water locations. While the model performs well in identifying eroded areas, misclassifications occur in regions with mixed vegetation cover or low


Introduction
There is a substantial gap between the state-of-the-art conceptual model of range management and the tools available to rangeland managers to apply that conceptual model.The core of the rangeland conceptual model is the ecological site description (ESD).As a classification system for rangeland, ESDs are grounded on recurring soil characteristics, landform variations, and geological and climatic conditions that affect the potential to produce a plant community with defined kinds, amounts, and proportions.An integral component of the ESD is the State and Transition Model (STM), a diagrammatic framework illustrating how an ecological site may transition from one stable plant community (state) to another.Thus, each ecological site is defined by the potential plant communities it can support.Transitions describe the potential trajectory from one state to another, including changes as a result of management actions [1].These concepts, collectively, offer a comprehensive understanding of land dynamics, that inform decision-making related to land management and restoration (For more details see: https://bit.ly/usdaars_esd).Bestelmeyer et al. [2] stated that STMs need to be created to describe vegetation dynamics in a particular landscape and assess the benefits of management practices in specific landscapes.STMs organize information related to management actions at specific landscape sections [3] and consist of four essential elements: states, transitions, plant communities (within states), and community pathways.This framework assists land managers with information on the state of the land [1] [4] and characterizes rangeland ecosystems and their responses to different types of disturbances [5].
Ecological sites are used by the main agencies working on rangeland management across the United States including the Bureau of Land Management, U.S. Forest Service, and Natural Resources Conservation Service) [6].The Natural Resources Conservation Service has developed National Ecological Site Handbook [7] explaining how ecological sites should be described, as well as a national database of ESDs [8].As a simple example of how a STM could be used to develop a management plan for a given area, a rangeland manager, knowing the ecological site and state, could examine the STM to identify accessible states and compare the potential costs of management actions against the benefits of changes in the flow of ecosystem services (for example forage, site stability, and wildlife habitat) for the possible alternative states.In practice, rangeland managers might have a map of ecological sites based on a soil map and with some field validation for management planning, but not a map of states.There are many complexities in making rangeland decisions, including the complexity of a landscape with many inclusions of other ecological sites, and the likelihood that within each management unit (pasture) there will be multiple ecological sites.
But without knowing states the manager cannot apply the knowledge about states and transitions in STMs except in a general way.Bestelmeyer et al. [3] (p. 364) highlight this knowledge gap: Readily available maps of ecological sites, and especially ecological states, are generally not available to assist planning.Without maps that are connected to ESDs and STMs, planners will find it difficult to use these tools.
One consequence of the lack of site-specific maps of ecological sites and states could be environmental degradation due to accelerated soil erosion.These areas would be classified as "eroded" states.Global soil resources are threatened, primarily due to anthropogenic activities [9], and rangeland areas are not the exception.Across western US rangeland, land degradation is an active line of research [10], and along with other factors, soil erosion represents one of the most common indicators for the assessment of land degradation [11] [12].Thus, monitoring and evaluating the risk of accelerated erosion by detecting "eroded" states at different scales and across rangeland ecosystems will assist policymakers, stakeholders, and ranchers in making better decisions concerning the consequences of management decisions [1] [3] [4] [5] [13].
Remote sensing has been of particular interest as a supporting tool for mapping and monitoring cover [14] [15] [16].Allred et al. [14] developed the Rangeland Analysis Program (RAP), a model based on a deep-learning approach that provides layers of perennial forbs & grass cover, annual forbs & grass cover, shrub cover, tree cover, and bare ground cover.This is a 1-dimension convolutional neural network model using 40 variables (climate, topography, and remote sensing-based data) that focuses on the temporal domain at the pixel (30 m) level for a per-pixel percent cover of each class output [14] [16].Remote sensing has also been used to map ecological sites and states [17] [18] given the ability of remote sensing to work at the landscape scale.Steele et al. [18] show an approach to map ecological states where manual methods were used to classify the states and recognize the use of remote sensing to improve these efforts.
Current practice is for rangeland conservation specialists to use vegetation patterns, soils, and climate knowledge to classify land into ecological site and state units.Soils can be sampled to determine soil series and ecological site, which is not possible with remote sensing.Range conservationists develop their skill in consultation with the narratives in ESDs, as well as quantitative information including species productivity, total annual production, canopy or foliar cover (over-and understory), vertical and horizontal structure, and ground surface cover [7].This is labor-intensive.Research in mapping ecological states has relied on manual methods [18] and acknowledged the need for more research on improving previous methods [19] including utilizing remote sensing data applied to rangeland management.As a first effort at incorporating some of the concepts to map the spatial arrangement of vegetation into remote sensing, we used UAVs and Structure from Motion (SFM) processing to capture information about the spatial arrangement of individual plants and groups of plants, which is not possible at the scale of satellite based remotely sensed imagery.Specifically, we considered 3 measures: bare ground, or the proportion of the surface not covered by vegetation; fetch, the distance from a random point to the base of the nearest plant [20] [21] and should not be considered as the distance between plants [21] (often used to assess range condition [20] [21] [22] [23] [24]); and the Largest Patch Index (LPI), or the percent of the landscape covered by the largest patch of bare ground [25] [26].
In this study, we developed a workflow to address the following objectives: 1) to evaluate landscape representation of ecological states compared to a few isolated polygons around long-term vegetation transects, 2) to assess the ability of remote sensing, enhanced with very high-resolution drone imagery, to classify states based on the season of imagery collection and spatial resolution, and 3) evaluate the accuracy across 3 bare patch metrics.We hypothesized: a) a better representation of the spatial extent and distribution of ecological states than isolated ground-based mapping; b) a better classification accuracy from higher resolution imagery; and c) that LPI metric represents better spatial distribution of bare ground.We also evaluated the ability of 3 sensing platforms to extrapolate UAV derived data to the landscape scale.

Study Area
The area of interest (AOI) for this research is in the Santa Rita Experimental Range (SRER, 31.817N,−110.85 b1W, ~1200 m) in southern Arizona on the ecological site Sandy Loam, Upland and Deep (SLUD) in Figure 1.This region typically exhibits a mean annual temperature of 18.4˚C and average annual precipitation of 358 mm [27] driven by the monsoon, with most of its rainfall occurring during the summer (Jun-Sep) and the rest during the winter (Oct-May).

Workflow
In Figure 2 we show a flow chart with the steps of the methodology presented in this work.Numbers in boxes work as a reference for the steps described in the text.

Ground Data
The baseline of this research is the field-based maps of ecological states in SRER.The extensive vegetation transect database from SRER, encompassing over 70 years of vegetation measurements carried out at intervals between one to seven years across 200 permanent transects, functioned as the primary reference point for identifying the areas that could be included in the mapping for this project.Between November 2018 and January 2019 Robinett Rangeland Resources performed fieldwork to map ecological states (ES) across the ecological SLUD sites.We refer to these as the surveyed polygons with known ecological state.A geospatial database with 92 polygons (Figure 2 box 1) was generated with a minimum polygon size of 30 m × 30 m and is available at.These polygons have irregular shapes due to the shape of ecological sites which are used as part of this process and defined using attributes of the plant community that included basal cover for both native and non-native grass species, canopy cover for the large shrubs (Prosopis velutina), median fetch values, and rangeland health indicators to determine the eroded state.Figure 1 shows a subset of the distribution of the polygons mapped by Robinet that were part of the drone survey used in this research.In the study area, the distribution of ecological states across the mapped polygons predominantly fell into three categories: large shrub non-native grass (38%), large shrub native grass (25%), and large shrub eroded (27%).The remaining polygons were classified as non-native grass (6%), native grass (2%), and large shrub (1%), comprising a smaller proportion of the total area.

UAV Data
In May and September of 2019, we collected high-resolution (1 cm) RGB drone imagery across 74 eco-state polygons within the SLUD ecological site.The objective of this data collection was to capture the seasonal variations across the study area.During the dry season in May, cover was lowest, and the green woody plants stood out for easy identification, whereas the grasses and forbs were largely senescent and challenging to identify.In contrast, in the post-monsoon season of September, the green herbaceous plants, being at their peak biomass, were more distinguishable.
For each polygon, 1 cm RGB images were produced [29].These images were used to generate a 5 cm land cover classification (Figure 2 box 3).The cover classification images from each period represent the main input used as part of the training data in the steps to build a model for scaling up the ecological states mapping.For more details regarding the protocol using a supervised Random Forest classification model implemented in Google Earth Engine (GEE).(See Gillan et al. [29] for a full list of derived products).
We performed a visual inspection of the 1 cm RGB images to assess, and we could not distinguish the invasive Lehmann lovegrass from the native species, and therefore we could not distinguish those states defined by native or nonnative grass.Instead, we used a simple grass presence to distinguish two states: large shrub eroded (LSE) and large shrub grass (LSG, combining non-native grass and large shrub native grass states).For ease of processing we decided to resample to a 5 cm cover classification, as presented in Gillan et al. [29], consisting of 4 classes, grass, shrub/tree, bare ground, and shadows, as shown in Figure 3.

Satellite Data
For this work, we used three different satellite platforms, Landsat 8, Sentinel 2, and Planet Scope.The first use of these platforms was to generate the boundaries (pixel-size grids) at which landscape metrics were calculated (Figure 2 Box 2).
We also developed predictive models, using satellite bands and satellite-derived vegetation indices as predictor variables to estimate the dependent variable (eroded state) for each sensor (Figure 2 Box 6).Table 1 shows a summary of the satellite data characteristics.We selected the best available images matching similar dates to drone surveys (late spring and late summer) in 2019.The main criteria in selecting imagery were images with the least cloud coverage in the period of interest.For late spring, we used mid-May through mid-June, and for late summer, mid-August through mid-September.We used Google Earth Engine (GEE) as the main geospatial processing platform.Landsat 8 Collection 2 and Sentinel 2 (SR) products are part of the core data catalog and are consistently ingested directly from official storage servers.For PlanetScope, we used the academic licensing offered by Planet [30] to access and download the data through their Application Programming Interface (API) planet for searching imagery and sending it out to store directly into GEE.
SR products are delivered in tiles (single images) covering large land extensions (~10,000 Km 2 ).We identified the tile codes covering our AOI, which in the case of Landsat 8 uses the World Reference System (WRS) with path 36 and row 38, and Sentinel 2 uses the Military Grid Reference System (MGRS) with a code identified as 12SWA for SLUD area.A single area can be overlapped by different tiles on different dates.Therefore, we clipped images to the extent (see Figure 1) of the SLUD in SRER, with an approximated area of 240 km 2 .To create a single image per period, we mosaicked the best images covering our area of interest.
During this step, the pixels with the highest surface reflectance value were taken in areas where two or more images overlapped.Other metadata from SR products were also used to review information such as scaling factors, data versions, and sensor and orbiting characteristics like sun and sensor view angle.
Managing these products through GEE allowed us to perform different operations on the data catalog, such as temporal searches, clipping, and metadata filtering, which was critical for assessing different images more efficiently.The scripts in GEE used to extract the images and the instructions to move data from Planet to GEE, are available in the supplementary information section.

Metrics to Support the Identification of Large Shrub Eroded State
Three metrics were calculated with the objective of identifying the best option for estimating the large shrub eroded state (referred to as eroded state for brevity), based on the 5 cm drone classification image for each period (May and September 2019) (Figure 2 box 4).Each option exhibited specific characteristics to determine the most accurate metric for classifying the eroded state.Here we used the satellite-based grids to compute each metric within each grid.The metrics calculated were the Largest Patch Index (LPI), Mean Fetch (MF) as used in [23], and a straightforward Bare ground (BG) percentage.These were chosen to represent the site classification criteria of fetch > 0.5 meter.Figure 3 illustrates an example of a grid (red) based on the Landsat resolution (30 m), with the white dots representing the centroids at each grid.The computation of these metrics was intended to act as a foundational building block for the landscape-scale classification of the eroded state.As such, this process should be understood as the initial step toward determining which pixels are classified in an eroded state.The point-based datasets (in shapefile format), representing the corresponding values of each metric for each period, were produced as a result of a set of processes developed to calculate each metric at each grid (for details see Github repository).In the subsequent sections, the application of these derived metrics to the broader landscape is described, detailing how the identified pixels were instrumental in deriving the landscape-level depiction of the eroded state.

Largest Patch Index
The Largest Patch Index (LPI) is a metric used in landscape ecology to measure a class's dominance within a spatial extent by quantifying the percentage of total landscape area comprised by the largest patch of the class [25].In reference to our eroded state, LPI provides an estimate of the proportion of the area that is covered by the largest bare ground patch.LPI is a metric of connectivity of vegetation patterns that can be related more closely to the process of characterizing ecological states [21] [22], specifically in rangeland ecosystems.
We generated a set of LPI values at each grid and used it as the variable to model in our scaling-up process that is further described.The reasoning behind the use of the LPI metric is the idea of using the connectivity of elements (i.e. 5 cm pixels) of a specific class (bare ground) to calculate a value indicating the level of pixels connected with the boundaries (satellite-based grid cell) in use. Figure 3 shows an example of a grid cell based on Landsat resolution (30 m) over the drone-derived 5 cm.classification image.LPI values were calculated in Google Earth Engine [31].The input parameters to calculate this metric are the classified image (raster), the boundary (grid), and the type of metric (LPI).LPI values were calculated at each grid cell across all the polygons surveyed and for each period, values are expressed in percentages (0 -100), see Figure 3 for an example showing the point layer with centroids (white dots).We used the class bare ground to calculate the LPI values that would support the task of discerning between the ecological states of large-shrub eroded (LSE) and large-shrub grass (LSG).The output of this process was a point-based file (shapefile) for each satellite-based grid, and it was stored in GEE.Each point has the LPI value corresponding to the grid cell used.The final sets of LPI points were clipped using an inner buffer to avoid edges where the 5 cm classification image does not fully overlap the satellite grids and could produce incomplete values of bare ground LPI.

Mean Fetch
For the Mean Fetch (MF) metric; we implemented a method presented by [23] in relation to the random point-to-plant (PTP) measures as an approach for estimating the average fetch of an area as a proxy for estimating the eroded areas.
As such, this provides an approximation of the central tendency (mean value) of the length of bare ground before encountering vegetation.It does not account for between the patchiness or size differences of bare ground amount between plants.We generated a set of randomly distributed points (n = 1000) within the delimiting area, the satellite-grid cell then took the average PTP distance value at each grid cell.This process was repeated for all the grids overlapping the surveyed polygons, and the values are expressed in meters.All the calculations to get MF values for each period were also performed in Google Earth Engine, and a dataset was exported out in a table format (for details see Github repository).
The resulting dataset has the same structure as that presented for the LPI; in this way, we were able to just add the MF output as another attribute to the existing LPI collection.

Bare Ground
The Bare ground (BG) percentage area is a common measure available in remote sensing products [10].Like MF, this metric does not account for spatial patchiness of bare ground patch sizes.In this research, the BG metric was obtained by getting the percentage of the area covered by the class bare ground at each grid cell and across the polygon surveyed and for the two periods of time where the classified image was available.This metric was computed in GEE, and values are expressed in percentages (0 -100).The BG dataset was incorporated into the structure used for the previous two metrics.This combined data is stored in a metric data table (MDT) format, the structure of which is detailed in Appendix A of the supplementary information.
The main input to calculate the metrics listed in the MDT is the 5 cm classification image derived from drone surveys along with the grids for each sensor.
The metrics from this dataset were considered as the "ground-data" it was used as a training dataset to build the Random Forest model, which is described later.

Assessment of Metrics
We evaluated the three different metrics (LPI, MF, BG) across each satellite platform to determine which metric best classified the surveyed polygons with known ecological site.We assessed these metrics using a classification assessment for both datasets, the ground-based data (drone-based data) and the modeled data.We utilized the information from the MDT to determine the threshold values at which each grid cell metric can be classified as eroded.
The first step in using a proxy metric to determine the state of the land (eroded areas) was to determine at which value (threshold) a piece of land (i.e.pixel size) should be classified as eroded.Since the distribution of our metrics was primarily unimodal, we decided to simplify our approach by using the mean value of the metric distribution as the threshold value for discerning between the LSE and LSG states.The average value is commonly used as the initial approach to test in several thresholding methods [32], especially when some of the thresholding techniques have very specific distribution requirements.
To do this, we obtained the distribution of each metric, for each resolution and time period, and from there, the mean value was calculated and used as the threshold value.This threshold was used to split grid cell values from the dataset into Large-Shrub-Eroded [LSE] (equal to or greater than the threshold) and Large-Shrub-Grass [LSG] (lower than the threshold).
Then, from the entire set of polygons surveyed, we filtered out the small polygons, those having an area smaller than half a hectare.The reasoning behind this filter was our consideration that polygons with less than half a hectare because they would few pixels and they are too small to warrant management attention.Given the spatial resolutions of Sentinel-2 (10 m) and Landsat (30 m), only a minimal number of pixels fall within these smaller polygons, in some cases just one or two pixels.This sparse pixel distribution introduces an imbalance, making the comparison across different sensors uneven and less reliable.Within each of the remaining polygons, we were able to calculate the percentage of grid cell values within each one labeled as LSE and LSG even though the polygons are irregular and grid cells were based on three different satellite resolutions.
In the case of the modeled data (RandomForest), we followed the same assessment (Figure 2 box 9) as in the drone-based data described above.

Random Forest Modeling Landscape-Scale Estimates
We built a model to scale up an erosion metric as a tool to identify the prevalent ecological states (Large Shrub Eroded and Land Shrub Grass) in the ecological site SLUD.We used the non-parametric decision tree-based model Random Forest (RF) [33] in regression mode (Figure 2 box 7) since we needed to model a quantitative metric.The RF modeling was performed using the implementation in GEE [31] based on Li [34].To model an erosion metric using a RF model we used the satellite data from the surface reflectance bands, red, green, blue, and near-infrared; and vegetation indexes Enhanced Vegetation Index (EVI), Enhanced Vegetation Index version 2 (EVI2), and Normalized Difference Vegetation Index (NDVI) as predictors.All the remote sensing features are available on the GEE platform.
By using an RF approach, we were able to model any of the three erosion metrics presented.The predictors were the same in all the cases, using information from three different satellite platforms, Landsat 8, Sentinel-2, and PlanetScope.
The training and testing datasets for each one of the metrics were generated using the dataset shown in the section above describing the metrics.We built the RF-based LPI predictor to run across the SLUD to model values of LPI at a pixel level allowing us to generate LPI maps for each period for the entire SLUD extent, as shown in appendix A. We performed a random selection of 75% of the drone-derived LPI (ground data) for use in training the model in regression mode, and the remainder was used to assess the accuracy of LPI values.One of the metrics to assess the accuracy available in GEE is the Out-of-bag error (OOB).Appendix A also shows the resulting variable importance of the modeling for each satellite and a table with the model OOB errors based on the training process.

Predicting Large Shrub Eroded Areas
The RF model allowed us to run predictions across SLUD and maps of LPI-bare ground were generated based on each satellite platform resolution (i.e. 30, 10, 3 m) (Figure 2 box 8). Figure 4 shows the example of the results based on modeling using Landsat 8 resolution (30 m).
The resulting prediction model is a GEE object that can be executed over any specific area of interest or re-trained with other datasets.This will simplify follow-up research intended to incorporate remote sensing datasets from previous years.The capacity of GEE allowed us to run these predictions across this extent (SLUD) in a few seconds, and results are stored as assets.

Connectivity of Eroded Areas
Using the threshold value, we masked out pixels falling below that threshold and then vectorize the remaining pixels (the process of converting raster data into polygons).The resulting vectorized layer consisted of a set of polygons (Figure 5(a)) in which we calculated the area and then filtered out areas smaller than 1  hectare (to represent the smallest management unit) and kept the remaining areas as part of the significant extent in Large Shrub Eroded state (Figure 2 box 10). Figure 5 shows an example of a smaller region with the RF model output (i.e.LPI) map (a) and the filtered version, showing polygons that build up (connected pixels using queen's case contiguity [35]) to an area of at least 1 ha or larger (b).The green polygon shown is the only LSE polygon of an area greater than equal to 1 ha.
We converted the resulting polygons back into raster images for better under-standing and faster visualization and refer to this layer as the LSE state layer (Figure 2 box 11).The GEE platform played a critical role in this process, mainly with Planet Scope and Sentinel 2 spatial resolutions, which required more computing resources and datasets always remain an asset inside this platform.The vectorization process produces small and unconnected polygons, as depicted in Figure 5(a); therefore, applying the 1-hectare or larger area threshold allows us to generate a cleaner and more homogeneous layer, as well as identify areas that are large enough to justify management attention.

Comparing the Estimates
We compared the classification of ecological states across the different remote sensing platforms to assess the influence of spatial resolution of the sensors.To do this, we first scaled up the estimation of the eroded state by deploying the model across the entire Sandy Loam Upland and Sandy Loam, Deep area in SRER.Then, the eroded areas of at least 1 ha or larger based on the model outputs from each period (May and September) were collapsed into a single layer per sensor by choosing only those classified as eroded regions in both periods.
We focused our comparison on the results using the LPI metric for the collapsed coverage showing eroded coverage for both the May and September images.This method was based on the LPI being the most accurate among the three metrics, and the collapsed coverage is the most parsimonious depiction of the eroded state, with agreement across seasons.
This strategy of collapsing the two seasons helped us to diminish some errors in the classification led by some ground cover characteristics, such as dry grass during early spring, where spectral signatures and vegetation indexes used in the model can mislead the classification as an eroded state.

Assessment and Landscape-Scale Validation Methods
The accuracy assessment involved ground-based known polygons that were identified by a rangeland specialist (Figure 2 Boxes 6 and 9).The specialist surveyed the study area and identified polygons that represented different ecological states, thus we refer to these as know polygons.Each polygon was labeled with its corresponding ecological state.These polygons were then used to evaluate the classification model's accuracy in identifying areas with eroded ecological state.The known polygons were overlaid on the classification model output to determine the percentage area of each polygon covered by eroded state.A minimum threshold of 50% coverage was used to determine whether a polygon was classified correctly by the model.
The final assessment compared the landscape-scale precision of the modeled distribution of the eroded ecological among the three satellite platforms.This used a visual comparison of the precision of the mapped depiction of the distribution of the eroded state.
The first landscape-scale validation used the main dirt roads across SRER as the basis for evaluation.This approach was chosen because the large dirt roads have no vegetation cover and should be classified as eroded state in the model-based representation.The location and extent of the roads was accurately digitized by manually drawing polygon shapes.The roads selected for this study were Santa Rita Road, Helvetia, and Highway 62.The modeled eroded areas from each sensor and time period were extracted using the road polygons to assess the accuracy of the results.The modeled eroded areas layers were then clipped using the road polygons to determine the percentage of eroded area.
The second landscape-scale validation used the location of livestock drinking water developments because livestock congregate in these areas which creates bare ground (http://bit.ly/srer_spatial).To validate our model outputs, we filtered our database for water sources located within modeled area and defined a box with an area of 2-ha around each water location.We then overlaid these boxes with the outputs from the model and extracted the percentage of area covered by eroded state model output within each box.Specifically, we considered a match when at least 50% of the area within the square showed eroded state.

Assessment of the Classification Using Drone-Based Metrics
Here we present the results of assessing the model in classifying ground-based polygons using three different metrics: the largest patch index (LPI), mean fetch (MF), and percentage bare ground (BG).First, we assess how well the metrics calculated using drone data performed in classifying the polygons as eroded (Large Shrub-Eroded -LSE) or non-eroded (Large Shrub-Grass LSG).We evaluated the three metrics to estimate the eroded areas by performing a classification assessment.From this assessment, for each polygon surveyed, we calculated the percentage of pixels in the eroded state based on each of the metrics, for each period and platform resolution.Figure 6 shows the results for the ground-based dataset (drone-5 cm classification).
In Figure 6, each bar represents a polygon surveyed with its corresponding known ground-based classification, LSG (Large Shrub Grass) in green and LSE (Large Shrub Eroded) in red.The blue dashed line is the 50% eroded-pixels value and is used as the threshold to classify the two ecological states and perform the accuracy assessment.We considered that 50% of eroded pixels or greater within a polygon is a significant portion of land to classify a polygon as an LSE.
Therefore, to consider a correct ecological state classification, LSE polygons are expected to have values above the dashed line and LSG below.
Polygons without a bar indicate that no metric values within a polygon crossed its threshold (mean of the entire metric distribution).In some cases, during the dry period (May) the polygon is classified correctly as "eroded" (% pixels eroded above the blue dashed line) but not the case during the wet period (September), this happens because of an increase of vegetative cover following the summer rains.

Drone-Based Metrics Classification Accuracy
We assessed the classification performance of each erosion metric using the statistical classification values for accuracy, sensitivity, and specificity.Table 2 shows the numerical values for the different statistics.
To compare the erosion metric across the statistical values presented in Table 2, we selected a value of at least 10% (0.1) as a threshold to identify patterns of difference across metrics, periods, and resolutions.In addition, in Table 2, the period labeled as "collapsed" represents the assessed values of only those pixels (grid-cells) that match the state (LSE and LSG) in both periods.The rationale behind this criterion was to address instances where senescent vegetation cover may be mistakenly interpreted as eroded.
Landsat in May for LPI shows a higher accuracy (0.91) compared to Landsat in May for MF (0.94) and BG (0.85).The Landsat in September for LPI (0.94) is also higher than Landsat in September for MF (0.91) and BG (0.79).This suggests that the LPI metric performs better than MF and BG for Landsat in both May and September.We can also see that Landsat May LPI value has a higher sensitivity (0.91) than Landsat in May for MF (0.69) and BG (0.95), indicating that LPI is better at identifying true positives compared to MF and BG for Landsat Resolutions: Landsat 30 m, Sentinel-2 (Sent 2) 10 m, PlanetScope (PScope) 3 m.Accuracy is a measure of the percentage of correctly classified as eroded by the classification method.Sensitivity tells us the percentage of polygons that actually are eroded and were correctly identified by the process as eroded.Specificity is a measure of how good the process is at ruling out polygons that are not eroded as the percentage of polygons that are not eroded and were correctly identified by the process as not being eroded.
in May.Similarly, Landsat in September for LPI (0.9) has a higher sensitivity than Landsat in September for MF (0.62) and BG (0.95).The specificity values for Landsat in May and September are generally high across all metrics, indicating that the models are performing well in identifying true negatives.
In this dataset, the use of a 10% threshold allowed us to identify patterns of differences across the different dimensions (sensor, period, metric) where polygons classification was performed.The dataset compiled involved different resolutions, periods, and metrics of bare ground amount and spatial pattern (LPI, MF, and BG) and the statistical assessment shows that the LPI metric has a substantially better performance in classifying the two ecological states.This is supported by the fact that for a metric like BG, there are several cases where statistical accuracy values are significantly lower than MF and LPI.In the case of MF, when compared with LPI, the cases where accuracy and sensitivity have the lowest values, are always higher than LPI.LPI never came out as the lowest value in any statistic across all the dimensions.This suggests that the relative size of bare patches is a better representation of the difference between the ecological states, than metrics that ignore spatial pattern of bare ground (MF and BG).

Assessment of the Classification of Polygons Using Model-Based Metrics
Here we present the results of assessing the model in classifying ground-based polygons based on Random Forest model data using the three different metrics: the largest patch index (LPI), mean fetch (MF), and percentage bare ground (BG).Figure 7 shows the results based on the classification using model data (Random Forest). Figure 7 presents the percentage of eroded pixels for each polygon based on the modeled data, following the same format as Figure 6.Specifically, each bar represents a polygon, and the height of the bar indicates the percentage of pixels classified as eroded.To classify a polygon as eroded, the threshold of at least 50% of the area covered by pixels in the eroded state must be met.

Model-Based Metrics Classification Accuracy
Across all sensor types, the LPI metric performs the best, consistently achieving high accuracy, sensitivity, and specificity values.For Landsat and Sentinel 2, the LPI metric achieves an accuracy of 0.97 in the May and collapsed periods and an accuracy of 0.91 in the September period.The accuracy for PlanetScope using LPI is slightly lower, with a maximum accuracy of 0.94 in the May and collapsed periods.
The MF metric shows a mixed performance across the different sensor types and periods, with some achieving high accuracy and sensitivity but lower specificity.The BG metric consistently shows the lowest performance across all metrics, with the highest accuracy value of only 0.91 for Sentinel 2 in the May period.
Overall, Table 3 highlights the importance of metric selection for accurately identifying eroded areas, with LPI consistently performing the best.It also shows the variability in performance across different sensors and time periods, suggesting the need for careful assessment in the processes involving the spatial scale up to areas beyond the SLUD area.Again, like the drone-based model of ecological states, this suggests that the relative size of bare patches is a better representation of the difference between the ecological states, than metrics that ignore spatial pattern of bare ground (MF and BG).

Assessment of Classification Precision among Satellite Platforms
The map in Figure 8 shows the results of the modeled coverage of the eroded state (LSE) across SLUD using the LPI metric for each satellite platform.In general, the eroded states are concentrated at the lower elevations (elevation increases from left to right), and are more discrete and isolated at higher elevations, and along roads (see Supplemental Information).Comparisons across three maps in Figure 8 of the SLUD region revealed varying eroded area coverage, influenced significantly by the resolution of the satellite platforms used.The Landsat-based map indicated more eroded area (5156 total eroded hectares), more continuous eroded coverage, and less precise representation of roads, which is consistent with the relative course 30-m resolution.
Whereas the depiction using the more precise 3-m resolution from Planet-Scope has less area in eroded state (3578 hectares), more discrete coverage rather than continuous, and more precise representation of roads.

Roads as Large Bare Ground Patches
Table 4 presents the results from the eroded areas calculated based on the roads polygons.
As expected, we found that the eroded area percentages for the roads were high (Table 4).However, the width of the road and season influenced the eroded area predicted by the model.The widest road, Santa Rita Road, showing the highest percentages of eroded area across all sensors and periods.However, the Helvetia and Highway 62 roads showed lower accuracy percentages, the lowest during September when green vegetation was abundant following the summer rains.The satellite resolution is most obvious in May, when Landsat was the least accurate.This can be attributed to the fact that these roads are narrower and can fit fewer pixels, resulting in less connectivity between eroded pixels.As a result, smaller, less continuous polygons of eroded area were predicted, leading to a lower percentage of eroded area compared to Santa Rita Road.

Livestock Waters as Bare Ground Patches
As expected, we found that the eroded area percentages around the livestock waters were high in 9 of the 11 cases (Table 5). 1) The water point locations (Water_Devs database) were used to identify eroded areas nearby, and then polygons were manually delineated around the affected regions using satellite imagery from Google Maps.A two hectares square was generated around the centroid of the polygon created.2) Only polygons covering an area larger than a Landsat pixel (900 m 2 ) were used.3) Those polygons with at least 50% eroded area were considered as eroded in the accuracy assessment.
This level of validation was most consistent among the satellite platforms in May when green vegetation was less common before the summer rains.In general, it was very likely (>70% accuracy) that the livestock waters would be identified as the LSE (eroded state) for all platforms and seasons.

Discussion
In this study, we present a workflow for mapping the ecological state of Large Shrub Eroded (LSE) areas within a specific ecological site, Sandy Loam Upland and Deep (SLUD), using a combination of drone, satellite data, and proxy metrics of bare ground extent and spatial distribution.Our results demonstrate the potential of this approach for scaling up the depiction of ecological states to a landscape scale, rather than the typical < 1 ha ground-based assessment along 20 -50 transects.This landscape scale depiction will be valuable for stakeholders in rangeland management, including conservation agencies, ranchers, and scientists because it identifies extent and location of important ecological states [3] [4] [13], like the eroded state that are easily unobserved in typical field assessment that are confined to small portions of the landscape referred to as key areas.Thus, we have satisfied our first objective to evaluate the landscape repre-sentation of ecological states compared to isolated polygons along long-term transects, and we support our hypothesis that the former is superior.By providing detailed and accurate information at the landscape-scale, our approach has the potential to improve management decisions and ultimately contribute to the long-term health and sustainability of these important landscapes, although a lot of work remains to develop a cost-effective mapping approach across large landscapes.

Using UAV and Satellite-Based Data
The utility of satellite-based remote sensing for mapping ecological states, ecosystem site potential, and plant communities has been well-demonstrated in several studies [18] [36] [37].Our research goes a step further by combining this approach with unmanned aerial vehicle (UAV) derived land cover data to map ecological states in semi-arid rangelands.Thus, we have satisfied our second objective to assess the ability of remote-sensing enhanced by very high-resolution drone imagery to classify states at the landscape-level.As we hypothesized, the drone-based imagery was critical to developing the classification model using each of the three land cover metrics (LPI, MF, and BG).However, the use of UAV systems for mapping ecological states is limited by the spatial extent they can cover mainly due to energy requirements.Addressing this limitation by integrating UAV data with publicly available satellite imagery can have a significant impact on managing rangeland resources, making this research an important contribution to the field.

Using a Proxy to Identify Eroded Ecological States
Our analysis showed that the LPI (Largest Patch Index) metric was superior to the MF (Mean Fetch) and BG (Bare Ground percentage) metrics in accurately predicting the eroded ecological state.In the context of mapping eroded areas, the LPI offers an advantage over MF and BG by not only quantifying the extent of bare ground but also its spatial configuration [4] [38].This unique characteristic allows for the detection of large, connected patches indicative of the significance of a class [26], which could be overlooked when using MF and BG, as they primarily capture the average exposure and overall proportion of bare ground respectively, without considering the spatial pattern.
This suggests that the relative size of bare patches is a better representation of the difference between the eroded and non-eroded ecological states, than metrics that ignore spatial pattern of bare ground (MF and BG).Thus, we satisfied our final objective to evaluate the accuracy of ecological state classification across three metric of vegetation cover and bare ground, and supported our hypothesis that the LPI would be the superior metric of the three.
By spatially predicting LPI values, setting a threshold, and defining 1 hectare as the minimal spatial unit to map, we generate rasters that can be used in rangeland management practices.Our goal was not to map ecological states at satellite resolutions (e.g., 3, 10, 30 m), but to develop a protocol based on a quantita-tive measurement, the LPI (proxy), that allows us to model areas with significant distances between plants, indicating the presence of bare ground.This approach provides a useful tool for classifying ecological states and can aid in rangeland management efforts.It is important to note that while the use of a proxy to estimate eroded areas may not provide information as detailed as other methods, such as field measurements or very high-resolution satellite imagery [18], it allows for a practical and cost-effective way to classify ecological states over large areas.This approach can help to identify areas that are particularly vulnerable to soil degradation and prioritize management efforts accordingly.Moreover, the use of a proxy like the LPI metric allows for a standardized and reproducible method that can be used across different sites and regions, enhancing the comparability and transferability of the results.

The Use of Remote Sensing Imagery for Scaling up
Mapping ecological states is a complex and ongoing challenge in rangeland management [4] [18].Our study focused on a specific ecological site, Sandy Loam Upland and Deep, and used a proxy to estimate eroded areas at a regional scale.However, mapping at the Ecological State level requires accounting for a wide range of factors that vary across different sites and landscapes, including soil type, topography, vegetation, climate, and disturbance history [39].These factors can interact in complex ways and result in unique ecological states and transitions that are difficult to map accurately.While remote sensing technology offers the potential to scale up mapping efforts, as reported in Ludwig et al. [17], it is important to recognize that there are limitations and challenges in applying these methods to different ecological sites.Thus, more research is needed to develop robust and transferable methods for mapping ecological states that can be applied across a range of landscapes and regions.
The most important feature that remote sensing technology offers is maximizing resources and generating scalable information that can meet the needs of different entities.However, the complexity of mapping at the Ecological State level requires the development of new methods, such as the use of proxy metric [18], as we have done in this study.It is not just enough to rely on the bands or vegetation indices of satellite or drone data; instead, developing new metrics that capture the important ecological features can lead to improved mapping efforts.

On the Accuracy of the Classification Process
Evaluating the accuracy of the classification process to ensure that the resulting ecological state maps are reliable and useful is a critical first step.In our study, we evaluated the accuracy of the classification process by comparing our results to ground-based polygons surveyed by an expert, as well as to known locations of dirt roads and water features.We found that our model performed well in identifying eroded areas, which is critical in identifying ecological states in this ecological site and more broadly across other rangelands.However, we also identified some misclassifications, particularly in areas with mixed vegetation cover or vegetation cover with low biomass and also a dependency on the phenological states, where models misclassify areas with senescent grasses.These limitations highlight the challenges in accurately mapping ecological states and the need for continued research and development of new methods to improve the accuracy of these maps, and especially the need to represent the full range of conditions to be mapped in training datasets.However, we also recognize that accuracy of the classification is inherently limited by the accuracy and precision of the base information.For example, the ecological site map has a spatial resolution that can have inclusions of nonconforming ecological sites (different soils or slopes) as large as 15 ha.Thus, we should expect some level of inaccuracy given any differences in the spatial accuracy and precision of the base data.

Validating Outputs
The validation process of our methodology indicated satisfactory results overall, with high accuracy in the maps produced compared to ground-truth polygons, dirt roads, and water locations.However, the accuracy varied by season of acquisition of the drone and satellite imagery, suggesting how important it is to identify the most parsimonious data sets that can best transcend seasonal variations.
However, it is important to acknowledge that the validation process has limitations in our workflow due to the limited number of ground-truth polygons surveyed, which could have contributed to misclassifications in certain regions of SLUD.In addition, the validation was based on a single date of field observations, which may not fully capture the temporal dynamics of the rangeland ecosystems.Therefore, further validation using more ground-truth polygons and incorporating multiple satellite observations from different dates into the model building process could improve the accuracy assessment of the mapping outputs.
While the validation process can be seen as the weakest step in our workflow, we are confident that our methodology provides a valuable contribution to the mapping of ecological states in rangeland ecosystems.

Implications for Interpreting Ecological States
Figure 8 and Appendix B highlight the fact that much of the eroded (LSE) state is located at lower elevations, which are also where temperatures are the highest, precipitation is the lowest, creating a high vapor pressure deficit that stresses plants, and limits plant abundance compared to higher and wetter elevations.
First, the fact that our model recognizes this gradient of plant growth potential is a measure of its validation.Second, and more important, recognizing this gradient is critical to the interpretation of the extent and location of the eroded state across the landscape.For example, attention should be focused to detect the expansion of the eroded state from lower to higher elevations, and from small patches to large continuous patches.

Further Work
The methodology presented in this study offers a promising approach for scaling up ecological states mapping in semi-arid rangelands with similar ecological sites, such as MLRA 41-3.By leveraging the power of machine learning algorithms, proxies for the abundance and distribution of bare ground (LPI, MF and BG), and remote sensing data, our workflow can serve as a blueprint for other research efforts and help advance the field of ecological states mapping.
In addition to refining the single season or year accuracy, the next direction for research will be the time series of ecological state extent and spatial distribution using the archive of more than 3 years of Landsat imagery, and shorter time series for Sentinal 2 and PlanetScope imagery.This retrospective analytical approach can be assessed using known land use, fire and climate information for accuracy assessments.Having such a temporal record of change across a landscape will identify the most vulnerable and most resistant parts of the landscape to ecological state deterioration.Thus, directing managers towards those more sensitive areas rather than spreading the few resources widely across the landscape.

Conclusions
This research represents an incremental step forward in mapping the ecological states of the Sandy Loam Upland and Deep (SLUD) site.By integrating drone and satellite data with proxy metrics for bare ground extent and spatial distribution, we expanded the traditional scope of ground-based assessments.This method enables us to study ecological states on a landscape scale, an approach which provides more comprehensive insights for rangeland management stakeholders.
We were successful in overcoming the spatial extent limitations of UAV systems by combining UAV-derived land cover data with satellite imagery.Furthermore, our study found the Largest Patch Index (LPI) to be a more suitable metric in predicting eroded ecological states than either bare ground or mean fetch.This highlights the significance of considering both the size and spatial configuration of bare ground when recognizing area with potential and expression of eroded conditions.
While the validation process of our approach yielded satisfactory results, it also underscored its limitations, thus emphasizing the need for continued refinement of our methods.The path ahead includes refining classification accuracy, leveraging the power of satellite imagery archives, and deepening our understanding of landscape-scale vulnerability to erosion.
for his valuable collaboration.Additionally, we would like to express our gratitude towards Alessandra Gorlier, Research Scientist at the University of Arizona, for providing guidance in the use of SRER data.Their contributions have significantly enriched this manuscript and are sincerely appreciated.The results of the percent of eroded area across precipitation ranges show a expected correlation between precipitation levels and the eroded areas.The precipitation ranges were divided into thirds, with ranges of less than or equal to 362 mm, greater than 362 and less than 390 mm, and greater than 390 mm (see Figure S2).Results from the three different satellites show an agreement between the platforms across these ranges, the drier the region, the larger the eroded areas.
The highest eroded areas were observed in the lowest precipitation range.In contrast, the highest precipitation range had the lowest percentages of eroded area based on the three platforms (Figure S1).In the case of annual mean temperature, the results also showed a clear relationship with eroded areas, as expected.We also partitioned the temperature values into thirds ranges to assess the gradients in relation to the percent of eroded areas.The ranges of temperature were ≤18.5, >18.5 -19, and > 19 (˚C) as in Figure S4.Areas with the highest temperature range, greater than 19˚C, had the most significant amount of eroded areas.
The areas with temperatures between 18.5˚C and 19˚C had a lower percentage of eroded area, between ~ 40% -65% of the total area.The lowest percentage of eroded area was observed in temperature less than or equal to 18.5˚C, with a very low presence of eroded areas (Figure S3).The results of modeled eroded areas and elevation gradients showed that, as expected, areas in the lower elevation zones of SLUD had the largest percentage of eroded area.Figure S6 confirms that across all sensors, eroded areas in the low elevation zones exceeded 70%, while those in higher elevation zones were significantly lower.

Figure 1 .
Figure 1.The study area is in the Santa Rita Experimental Range (black), the Sandy Loam Upland and Sandy Loam Deep polygons (blue), and the ground-based ecological states map polygons (green) surveyed by Robinett Rangeland Resources in 2019 (Source: https://bit.ly/srer_es_polys).

Figure 2 .
Figure 2. The workflow for assessing the importance of satellite resolution, season, and vegetation metrics for modeling ecological states.Numbered boxes provide a reference for the steps described in the text.(BG: Percent Bare Ground, LPI: Largest Patch Index, MF: Mean Fetch).

Figure 3 .
Figure 3.An example of a drone-derived 5 cm cover classification for one of the polygons with an approximate area of 4 hectares.The four classes identified are grass, shrub/tree, bare ground, and shadows.A Landsat image-based grid (red dashed line) with its centroids (white) where each grid cell represents a 30 m pixel was used as the boundary to calculate erosion metrics.

Figure 4 .
Figure 4.An example of a full extent modeled LPI per pixel for the entire SLUD based on Landsat 8 resolution (30 m) for the September period.Higher values of LPI are shown in warmer colors.

Figure 5 .
Figure 5. Example of the eroded state map produced with the RF model (a) and the vectorized LPI format of the filtered (≥1 ha) polygons (b).

Figure 6 .
Figure 6.The percentage of eroded pixels within each polygon (≥0.5 ha) for the three different metrics, (a) LPI, (b) Mean Fetch, and (c) Bare ground for the periods May, Sep and Collapsed (both).The x-axis shows the polygon size (ha.) as the label for each polygon.LSE is represented by the red bars and in green LSG.The dashed blue line is at 50% eroded pixels percentage.

Figure 7 .
Figure 7. Model-based data on the percentage of eroded pixels within each polygon (≥0.5 ha) for the three different modeled metrics, (a) LPI, (b) Mean Fetch, (c) Bare ground for the periods May, Sep and Collapsed.The x-axis shows the polygon size (ha.) as the label for each polygon.LSE is represented by the red bars and in green LSG.The dashed blue line shows the 50% eroded pixels percentage.

Figure 8 .
Figure 8. Predicted eroded state across SLUD based on Landsat 8 ((a), blue), Sentinel 2 ((b), red), and PlanetScope ((c), green) based on the collapse of May and September results.The color scheme has been deliberately utilized to clearly distinguish between the predictions yielded from each individual satellite source.

Figure S1 .
Figure S1.Precipitation gradients across SLUD.The figure displays the partitioning of the study area of SLUD into three precipitation ranges, which include low precipitation in red (below 362 mm), intermediate precipitation green (362 mm -390 mm), and high precipitation in blue (above 390 mm).

Figure S2 .
Figure S2.Percent of eroded area using the model outputs for each platform (Landsat, Sentinel-2, and PlanetScope) across each precipitation ranges.

Figure S3 .
Figure S3.Mean temperature gradients across SLUD.The figure displays the partitioning of the study area of SLUD into three mean temperature ranges, which include higher temperature in red (>19 (˚C)), intermediate temperature in green (>18.5 -19 (˚C)), and lower temperature in blue (≤18.5 (˚C)).

Figure S4 .
Figure S4.Percent of eroded area using the model outputs for each platform (Landsat, Sentinel-2, and PlanetScope) across each temperature ranges.

Figure S5 .
Figure S5.Elevation gradients across SLUD showing the partitioning of the study area of SLUD into three elevation ranges, which include low elevation in red (below 1061 m), intermediate elevation in green (1061 m -1150 m), and high elevation in blue (above 1150 m).

Figure S6 .
Figure S6.Percent of eroded area across elevation ranges for results from each satellite at SLUD.

Table 1 .
Characteristics of the satellite platforms.Surface reflectances from all the platforms utilized in this work.

Table 2 .
Summary statistics for drone-based ecological states mapping.

Table 3 .
Summary statistics for model-based ecological states mapping.

Table 4 .
Eroded area at each road in hectares and the percentage of the total area.

Table 5 .
Percent of eroded area nearby water locations based on the resulting layers of eroded state for each sensor.