Shrub Encroachment Early Detection System (SEEDS): a rangeland conservation tool

Landsat 8 Vegetation Indices and Woody Cover Modeling

Woody vs non-woody classification accuracy and error matrix for the 2015 NAIP imagery is provided in Table 1. A total of 140 sample points across 7 cover classes (20 points per class) were used to investigate the relationship between NAIP-derived woody cover (true color imagery example = Fig. 4A; Fig. 4B = derived woody cover) and Landsat 8 NDVI and MSAVI (Fig. 5). These indices were chosen due to their documented performance and application within arid/semi-arid ecosystems including the SRER (Huete & Jackson, 1987; Hueng et al., 2018; Holifield-Collins et al., 2015). The best-fit lines for both NDVI (RMSE = 7.89%, R² = 0.8, p ≤ 0.001) and MSAVI (RMSE = 10.8%, R² = 0.72, p ≤ 0.001) were quadradic, wherein shrub cover became asymptotic at ~70%. Areas with high cover were concentrated near washes and lower montane elevation zone. Landsat 8 NDVI was subsequently used for further modeling given its lower RMSE compared to MSAVI (Fig. 6). Furthermore, high-resolution NAIP imagery (Fig. 4A) availability is limited in time and space with a relatively short collection history compared to the moderate-resolution Landsat, which has been collected every 8 or 16 days since 1982. This flexibility in the Landsat dataset allowed for broader geographic applications of the final SEEDS model beyond the SRER and over longer time-periods. With a 7.9% woody cover error, the Landsat-based woody cover estimates are elevated slightly compared to the cover estimate obtained from classified imagery (Fig. 4B). However, high versus low woody cover areas were accurately modeled using only the Landsat 8 NDVI values.   

Figure 4. Natural color 1-m 2015 NAIP imagery (A) and resulting 30-m shrub cover (%) (B) raster-based shrub/non-shrub classifications. The shrub cover raster was used in conjunction with Landsat 8 vegetation indices (e.g., NDVI and MSAVI) to develop algorithms to quantify patterns of shrub cover on historical Landsat imagery.

Table 1. Error matrix accuracy assessment of the classification performed on July 2, 2015 NAIP images of the Santa Rita Experimental Range. The reference data for the accuracy assessment was the 2017 NEON Airborne Observation Platform imagery (10-cm spatial resolution).

Figure 5. Classified USDA NAIP (National Agricultural Imagery Program) and Landsat 8 NDVI (Normalized Difference Vegetation Index) (A) and MSAVI (Modified Soil Adjusted Vegetation Index) quadratic regression (B). Points (n = 140) were equally stratified along 7 classes in 10% shrub cover increments for model development. The NDVI relationship was used for modeling woody cover across the Santa Rita Experimental Range (see Fig. 6).

Figure 6. Classified woody cover in 2015 (9 cover classes in 10% increments) modeled using NAIP (National Agricultural Imagery Program) woody cover and Landsat 8 OLI NDVI (Normalized Difference Vegetation Index). The maximum and minimum percent woody cover values were 88% and 0%, respectively. The densely vegetated areas (e.g., washes and foothills) appear orange to dark red due to their high shrub cover. Blue areas are bare soil patches and/or Pleistocene aged surfaces with high clay content and low shrub cover (see also Fig. 1A).

Modeling Woody Cover Based on Biophysical Variable Importance

Landsat 8 NDVI values show patterns of shrub cover but do not provide insight into the biophysical/ecological variables that determine cover site potential. A total of 31 independent variables pertaining to seasonal and annual mean precipitation (PPT) and temperature (T), topography, soil characteristics, and ecological site were used in a stepwise regression, support vector machine (SVM), and random forest (RF) modeling procedures for determining 1) the best/most accurate variables governing woody cover potential; 2) the relative importance of the 31 biophysical variables potentially influencing maximum woody cover; and 3) the SEEDS/ShrubRisk model best suited to represent woody plant encroachment risk. A graphical comparison of the three models is shown in Fig. 7.

Figure 7. Stepwise regression (A), support vector machine (B), and random forest (C) woody cover model predictions versus 2015 Landsat NDVI-based woody cover for the Santa Rita Experimental Range. Per-pixel cross-validation (n = 45,839) was used for assessing each model’s prediction performance. Blue lines represent a 1:1 fit; red line is the linear regression of the actual fit.

The stepwise regression modeling (RMSE = 5.6%, R² = 0.457) revealed 12 variables to be the most important: slope aspect; depth to restrictive layer/bedrock; soil clay content over 1-m depth; depth to maximum clay content; maximum 1980-2010 T; distance to the nearest wash; mean winter PPT; minimum spring and autumn T and maximum spring and summer T. Although the stepwise error was relatively low (RMSE = 5.6%), the prediction utility was poor (Fig. 7A), and the temperature variable/model variability was unconstrained as reflected in a low R² (R² = 0.46). The SVM model improved the overall model prediction errors (Mean Absolute Error [MAE] = 4.2%; RMSE = 5.5%) and accounted for more of variation (R² = 0.63) compared to the stepwise regression model. The best modeling technique was found to be the 2015 random forest with the lowest modeling errors (MAE = 3.22%; RMSE = 4.2%), best fit (R² = 0.70), and strongest prediction utility (Fig. 7C). The most important biophysical variables/predictors for determining woody cover potential, as determined via Gini Index of Impurity, were elevation, soil clay content, slope aspect, autumn PPT, winter PPT, slope inclination, and distance to the nearest wash (Fig. 8).

Figure 8. Gini Index of Impurity ranking of the importance of variables predicting woody cover on the Santa Rita Experimental Range.

Model Testing Back-in-Time

Using the same biogeophysical predictor variables as the 2015 random forest model (e.g., elevation, soil clay content, slope aspect, autumn PPT, winter PPT, slope inclination, and distance to the nearest wash), the 1995 RF regression model (Fig. 9A) had strong predictive power with a low error (RMSE = 3.93%) and a model variability constraint R² of 0.74. The variable importance of the 1995 RF model vs the 2015 random forest model ranked elevation higher than autumn PPT and the importance of soil clay content was reduced compared to topographic variables (e.g., slope aspect and slope inclination). The differences in predictor variable importance are possibly due to PPT patterns (e.g., drought vs wet years) during the autumn and winter months across the SRER. Drought conditions can underestimate woody cover estimations based on Landsat VIs that are dependent on greenness to detect a vegetation signal. Alternatively, wet years with increased herbaceous and woody plant growth can make discerning the two cover types more difficult. The 2005 random forest model (Fig. 9B) also ranked topographic variables higher than the 2015 model, but its performance was still quite robust (RMSE = 2.64%, R² = 0.76). The consistency in low error rates and greater amount of variance explained in the multiple random forest compared to that of SVMs or stepwise regression techniques models argues for their use in modelling woody cover.

Figure 9. Random forest models of woody cover potential on the Santa Rita Experimental Range in (A) 1995 (max = 82.6%, min = 0%), (B) 2005 (max = 80%, min = 0%), and (C) 2015 (max = 77%, min = 0%). The densely vegetated areas (e.g., washes and foothills) appear orange to dark brown, reflecting their high shrub cover. Yellow areas are dominated by bare soil patches and/or Pleistocene-aged soils with high clay content and low shrub cover.

The random forest model predictions for maximum potential woody cover were used in finalizing the proposed SEEDS/ShrubRisk models (Fig. 10A-C) for determining a site’s woody plant encroachment (WPE) risk. Risk was calculated as the ratio of a site’s shrub cover on a given date (estimated by the Enhanced Vegetation Index for a given year: 1995, 2005, 2015) to its potential shrub cover. The wildfire that occurred in 1994 in the center of the SRER was detectable in the 1995 potential woody cover model (Fig. 9A, light yellow) and in the resulting WPE risk map (Fig. 10A, dark brown). Much of the burned region was classified as moderately high to high WPE risk in the 1995 due to the potential for velvet mesquite, a vigorous sprouter tolerant of top-kill, to re-establish post-fire. The moderately high to high risk areas within the area burned in 1995 then decreased in their risk level slightly by 2005 (Fig.10B), whereas many of the high risk burned areas transitioned to moderate/moderately high. By 2015 (Fig. 10C) woody cover was approaching the site’s woody cover potential, and its risk level declined compared to 2005. Regions dominated by bare soil patches and/or Pleistocene-aged soils with high clay content and low shrub cover potential remained at low to moderately low risk of shrub encroachment over the entire 20-year period (Fig.10).

Figure 10. Classified woody plant encroachment risk calculated as the ratio of a site’s shrub cover on a given date (estimated by the Enhanced Vegetation Index for a given year: 1995, 2005, 2015) to its potential shrub cover. The five risk classes range from Low to High and were defined by Natural Breaks (Jenks, 1967). The areas at most risk for shrub encroachment appear orange to dark brown; yellow to blue areas are low to moderate risk areas reflective of low shrub cover potential and/or areas that have reached their potential woody cover.

The areal extent of the various WPE risk levels for 1995, 2005, and 2015 is summarized in Table 2. Areas at the most risk have been decreasing in their areal extent over time ostensibly reflecting the progression toward maximum potentials on the various topoedaphic settings. Accordingly, the low to moderately low levels have been increasing in area. It should be noted that the areal extent of high risk areas in 1995 (5.6% of the SRER land area) is likely inflated due to the 1994 wildfire, while moderately high risk areas remained rather stable between 1995 and 2005 (24.6% vs 24.9%, respectively). The greatest change in WPE risk occurred between 1995 and 2015, wherein moderately low risk areas increased by 10.0% or ~5,200 acres (Table 3) reflective of the reduction in moderate to high risk sites

Table 2. Percentage of the ~21,000-hectare Santa Rita Experimental Range classified in various woody plant encroachment risk level in 1995, 2005, 2015.

Table 3. Absolute percent change in WPE risk level area over three time steps (e.g., 2015 vs 2005, 2015 vs 1995, and 2005 vs 1995). Calculated from percentages provided in Table 1.

Percent area calculations (Table 2) and acreage estimates (Table 4) for each WPE risk level can be used by stakeholders to estimate costs to implement brush management treatments on a landscape, watershed, pasture or allotment. Costs will vary depending on treatment method chosen (e.g., individual plant, chemical, mechanical, etc.), but an accurate estimate of acreage for potential treatment will aid in budget calculations and future planning. Furthermore, assessment of WPE risk level change through time across a landscape allows for understanding the trend of WPE in a given area of interest. For example, a decrease in high WPE risk acreage overtime (1995-2015 in Table 4) indicates more of the landscape has approached it’s carrying capacity for shrubs. While this may seem like a positive to decrease in high WPE risk-level acreage, these lands will require more intensive management treatments to reduce shrub cover and maintain a grassland state as these sites are close to or at their maximum shrub cover. In other words, it will cost more to treat more shrubs, especially if using a mechanical or chemical technique.

Table 4. Areal extent (acres) of landscapes at various levels of risk for woody plant encroachment on the Santa Rita Experimental Range in 1995, 2005 and 2015.

The random forest modeled woody cover for 2015 was validated using independent field data (Fig. 11) obtained from long-term transects (1953 – 2021; n = 131). Samples of predicted woody cover obtained by averaging the pixels that intersected transect locations using the ‘Sample’ tool in ArcGIS Pro (v2.6.1). We then calculated MAE, RMSE, and R² for model vs. transect cover estimates. Matching ground-based woody cover transect data were available only for model year 2015.

Figure 11. Random forest modeled woody cover across the Santa Rita Experimental Range in 2015. Points denote locations of ground-based woody cover measurements from 2015 on long-term transects (n = 131). White areas on maps are sandy wash (or drainage) and rock outcrop areas with missing maximum clay data.

While random forests emerged as the best machine learning technique, the ground- vs model-based shrub cover validation comparison suggested improvements could be made (Fig. 12). A model estimated vs field-observed woody cover RMSE of 13.9% is comparable however to similar continuous woody cover mapping studies. Xian et al. (2013) estimated shrub cover with RMSEs of up to 15.8% in Arizona using regression trees, where Jones et al. (2018) used random forests to model shrub cover across the entire Western US with RMSE values ranging from 7.1% to 11%. The slight RMSE increase of ~2.9% in our random forest model as compared to Jones et al. 2018 (e.g., 11% with SageSTEP ground-based observations) may be attributable to variation in transect measurement protocols and the spatial distribution of transects across the SRER (Fig. 11), as well as the transect orientation relative to estimated woody cover raster cells. Numerous types of error are associated with field-based vegetation cover sampling which is also highly influenced by the observer-based error. Taken together, SEEDS/ShrubRisk model estimates are not necessarily better or worse than field-based estimates and on-the-ground monitoring may be warranted to attain additional site-level vegetation cover information that’s beyond the scope of SEEDS/ShrubRisk. Regional and local vegetation mapping with random forests has proven beneficial in this application in development of a WPE risk framework, and supports the utility of machine learning, model-based approaches for application in conserving Sonoran Desert grasslands.

Figure 12. Final random forest modeled woody cover (%) relative to ground-based estimated woody cover (%) in 2015. Red line is the linear fit between modeled and observed values. Inset: random forest estimated woody cover evaluation metrics (number of samples, N; mean absolute error, MAE; root mean square error, RMSE; and coefficient of determination, R²) calculated using the Santa Rita Experimental Range long-term ground-based transects for 2015.

WPE Risk and Disturbance

The WPE risk/SEEDS framework sensitivity to fire was assessed on an area experiencing wildfire, where wildfire would ostensibly reduce shrub cover potential and thus increase the risk for WPE. Following years after a wildfire, a site would move closer to its potential cover and WPE risk would be reduced. The 2017 Sawmill fire (Fig. 13) burned 3,340 ac in the spring with only ~5% of the fire experiencing moderate burn severity and remaining areas at low severity or left unburned (Fig. 13 inset). Three years post-fire (2020), there remained 302 ac of moderately high and 213 ac of high WPE risk relative to 2015. Many areas of moderately low WPE risk were replaced by moderate risk from 2015 to 2020, but many of the moderately high to high risk areas in or near sandy washes changed little through time (Fig. 9). The high clay (maximum = ~47.5%), Pleistocene-aged alluvial fan in the southern portion of the Sawmill Fire remained at low to moderately low WPE risk over ~25 years.

Figure 13. Woody plant encroachment (WPE) risk maps for 1995, 2005, 2015, and 2020 within the 1,352-ha area burned by the Sawmill Fire (April 23, 2017) that occurred on the Santa Rita Experimental Range. Inset: Burn severity from the Monitoring Trends in Burn Severity remote sensing-based dNBR (differenced Normalized Burn Ratio) products for the 2017 Sawmill Fire. White areas in the burn area indicate missing data. White areas on WPE risk maps are sandy wash and rock outcrop areas with missing maximum clay data.

The lack of a moderate-to-high risk conversion with the Sawmill Fire suggests that SEEDS/WPE risk’s response to fire disturbance is relatively short-lived (e.g., ≤ 3-y post-fire). The differences in response to fire within the burn boundary could be due to contrasting soil clay content, where the loamy upland areas (southwest portion of the burn area; Fig. 13) within the Sawmill Fire experienced little increased WPE risk as compared to the sandy loam upland areas. The loamy, southwestern portion of the Sawmill Fire was also subjected to a lower severity burn. A higher severity burn may have facilitated a conversion of low to high WPE risk classes. Follow-up comparative work between ecological site shrub cover responses to prescribed fire and fire severity may inform more precise timing for implementing controlled burns as effective brush management strategies.

Conclusion

With the advancements in cloud-computing and the robust accuracy and reliability of machine learning modeling, a new decision support tool for proactive, as opposed to traditionally reactive, brush management was developed. Our WPE risk/SEEDS framework offers a cost-effective, and preventative alternative to reactive brush management on western rangelands. Furthermore, assessment of a site’s risk of encroachment helps to increase effectiveness and reduce the cost (money and time) of rangeland monitoring effort by pinpointing specific sites most at-risk, while simultaneously facilitating planning efforts for the location, type (e.g., prescribed fire), and timing of brush management. In an effort to benefit our stakeholders, we provide the most important biogeophysical determinants of shrub cover as indicated in our random forest model across Major Land Resource Areas 41 (Southeastern AZ Basin and Range) and 40.1 (Upper Sonoran Desert) including the 2015 WPE risk map for the SRER in an interactive online tool, called ShrubRisk (ShrubRisk.app).

Participation Summary

This project does not contain a stakeholder academic interview component for use in analysis per se, but instead included general participation to guide research applicability and development of the SEEDS/ShrubRisk application. Ranchers/producers and agriculture professionals participated in local meetings and field days/tours, where this project was highlighted and discussed, including the 2020-2022 Altar Valley Conservation Alliance Community Meetings and the 2019-2022 Las Cienegas National Conservation Area Bio Planning Field Days (see Educational and Outreach Activities section). Results on the project have been shared with the Rangelands Partnership at their Annual Meetings and with regional stakeholders through technology panel discussions and field days. The latter includes the Altar Valley Conservation Alliance and other stakeholders from the Arizona State Land Department; US Fish and Wildlife Service; Bureau of Land Management; Pima County, AZ Natural Resources, Parks, and Recreation; and The Nature Conservancy. Collaboration with CCAST’s Grassland Restoration Community of Practice has allowed for completion of a new Case Study and ArcGIS StoryMap including a PDF handout for dissemination. The web-based Brush Management topic pages are finalized, and the Rangelands Gateway website updated with new text and online resource links/URLs. The prototypic geospatial application (ShrubRisk) is completed, and the official version 1.0 was released June 2022. Stakeholder feedback will continually be solicited through the Rangelands Partnership newsletter and various University of Arizona email lists, with feedback incorporated by the CCT group.

Appendix A