Rangeland Fractional Components Across the Western United States from 1985 to 2018

Abstract

Monitoring temporal dynamics of rangelands to detect and understand change in vegetation cover and composition provides a wealth of information to improve management and sustainability. Remote sensing allows the evaluation of both abrupt and gradual rangeland change at unprecedented spatial and temporal extents. Here, we describe the production of the National Land Cover Database (NLCD) Back in Time (BIT) dataset which quantified the percent cover of rangeland components (bare ground, herbaceous, annual herbaceous, litter, shrub, and sagebrush (Artemisia spp. Nutt.) across the western United States using Landsat imagery from 1985 to 2018. We evaluate the relationships of component trends with climate drivers at an ecoregion scale, describe the nature of landscape change, and demonstrate several case studies related to changes in grazing management, prescribed burns, and vegetation treatments. Our results showed the net cover of shrub, sagebrush, and litter significantly (p < 0.01) decreased, bare ground and herbaceous cover had no significant change, and annual herbaceous cover significantly (p < 0.05) increased. Change was ubiquitous, with a mean of 92% of pixels with some change and 38% of pixels with significant change (p < 0.10). However, most change was gradual, well over half of pixels have a range of less than 10%, and most change occurred outside of known disturbances. The BIT data facilitate a comprehensive assessment of rangeland condition, evaluation of past management actions, understanding of system variability, and opportunities for future planning.

1. Introduction

Rangelands are a spatially extensive land cover type worldwide and in the Western United States. Many ecosystem goods and services are provided by rangelands including wildlife habitat, livestock forage, erosion control, recreational, and cultural resources. At the same time rangelands are subject to increasing pressure from climate change, alterations to historical fire dynamics, energy development, and recreation. Moreover, many western United States rangelands have large interannual variations in climate conditions (e.g., [1,2]), gradual change related to succession and long-term climate change [3] and are affected by perhaps the second worst drought in 1200 years [4]. The high evaporation rates coupled with the low, and highly variable, precipitation make arid and semiarid regions acutely sensitive to climate variability and various other perturbations.

Remote sensing offers the means of evaluating abrupt and gradual rangeland change at unprecedented spatial and temporal extents [5]. Dense time-series data over regional extents provide opportunities to comprehensively analyze the condition of rangeland resources, evaluate the success of past management actions, understand system variability, and plan for the future. Many of these efforts have relied on Landsat imagery (e.g., [6,7] whose 30-m resolution typically reflects the mean landscape structure [8] to quantitatively monitor change at ecosystem scales. This approach increases product accuracy, increases the detection of subtle changes related to ecosystem health and condition, and facilitates the determination of change drivers [9]. Data at these scales allows researchers to better match the spatial and temporal scale of monitoring data to that of the phenomena being studied [3,8].

Ecosystems are in a constant state of flux, some of this change may be considered noise in one application, but may be critical signal in another [9]. For most rangeland applications, fractional cover maps are needed to successfully parameterize state-and-transition models (e.g., calculate transition probabilities among states). Fractional cover maps can capture abrupt (e.g., fire, and land clearing), gradual (e.g., succession, and drought), and seasonal (e.g., interannual variability due to weather) changes [10]. Though not a direct measure of rangeland health, fractional cover maps can be used to determine the effectiveness of past management practices and to improve future management, when evaluated over long time scales.

Several recent studies have created time-series fractional component datasets across the western United States. Jones et al. [6] used National Resources Inventory (NRI) and Assessment, Inventory, and Monitoring (AIM) field data, Random Forest models, and Google Earth Engine to predict the fractional cover of components across the western United States over the Landsat 5–8 archive. The National Land Cover Database (NLCD) Back in Time (BIT) dataset quantifies the percent cover of rangeland components (bare ground, herbaceous, annual herbaceous, litter, shrub, and sagebrush (Artemisia spp. Nutt.) across the western United States using Landsat imagery at an annual time-step from 1985 to 2018 [3]. The BIT dataset was produced using extensive field observations, change detection, and Regression Tree (RT) models. Previous generations of BIT predictions produced only for the Great Basin were analyzed by Rigge et al. [3] and Shi et al. [11]. The current generation of BIT components have been robustly validated through comparison to high resolution imagery [12] and long-term field monitoring plots [13].

Relative to the previous generation of BIT products described by Rigge et al. [3], we have made several key improvements in the current iteration. First, we have pooled training data through space and time to make more comprehensive and robust models for each component. Next, we implemented more aggressive change detection and simultaneously used algorithms to identify and remove noise in the time-series predictions. Our main objectives are to fully describe the methodology used to produce the BIT data and investigate component trends over the time-series. The overall BIT approach emphasized (1) a temporally harmonious product and (2) detection of gradual change. Secondary goals are to evaluate correspondence in ecoregion-scale climate and component cover trends, assess component relationships with disturbance, and describe the nature of landscape change. We demonstrate case studies of BIT data application in several resource management concerns.

2. Materials and Methods

2.1. Study Area

Our study area encompasses nearly 3 million km² across much of the western United States, of which 71% was identified as rangeland [14] (Figure 1). We include most of the base mapping area described by Rigge et al. [14] but exclude several Channel Islands and small sections of the California coast. The region spans 33 U.S. Environmental Protection Agency (EPA) level III (i.e., moderate-scale) ecoregions [15], 16 states, and a large range of biophysical conditions and vegetation communities. Our mapping region includes nearly all the sagebrush biome, as defined by Jeffries and Finn [16]. Much of the region is federally managed by the Bureau of Land Management (BLM) or U.S. Forest Service (USFS). Fire is the most widespread disturbance, with 14.2% burned between 1985 and 2018 according to Monitoring Trends in Burn Severity (MTBS) data. Other frequent drivers of change are interannual climate variation, energy development, mining, infrastructure (pipelines, roads, etc.) development, vegetation treatments, and grazing.

2.2. Method Overview

The BIT suite quantified the percent cover of bare ground, herbaceous, annual herbaceous, litter, shrub, and sagebrush (Artemisia spp. Nutt) components across the western United States using 30-m Landsat imagery from 1985 to 2018. We use an automated method, extending from Rigge et al. [3] to identify change in spectral conditions between each year in the Landsat archive and the circa 2016 base map [14] (Figure 2). Next, we use the unchanged portion of the circa 2016 base map to train RT models predicting component cover in each year. We pooled training data in areas identified as having no spectral change relative to the base across all years in the time-series and used a series of procedures to remove the most spatially and temporally common values. Yearly fractional component cover outputs were inserted in the changed area while the base map was maintained in the unchanged area. We used a series of post-processing models to ensure accurate and unbiased post-burn trajectories and eliminate noise and illogical change in the predictions. Processing occurred on individual Landsat path and rows which were subsequently mosaicked to a regional level.

2.3. Base Map-Training Data Development

The training data for the BIT time-series were obtained from the circa 2016 base fractional component cover maps developed across the western United States [14]. Briefly, these data were produced by training a series of 2-m high-resolution satellite images (n = 331) using observations collected in the field. High-resolution predictions for each component were completed using Cubist RT models [17]. For a complete description of components, see Supplement 1 of [14]. Because component cover can change throughout the course of a growing season [5], our predictions were targeted to capture the highest vegetation cover within the year. Using the degraded high-resolution training, field plots collected directly at a Landsat scale (n = 5382), and BLM AIM (AIM) plots (n = 3226) as training, we predicted component cover at a regional scale using three seasons of Landsat imagery. Non-rangeland pixels (forest, urban, etc.) were masked out from the final products [14]. The circa 2016 base fractional cover products used to train the BIT were validated using a series of independent field data collected across the western United States (n = 1860) and with cross-validation relative to training data. Component accuracy at independent validation sites across all components had a mean coefficient of determination (R²) of 0.46 and a root mean squared error (RMSE) of 10.37%, and cross-validation accuracy had a mean R² value of 0.72 and an RMSE of 5.09% For more detail on the methods used to prepare the base fractional components, see [14].

2.4. BIT Process

2.4.1. Image Processing

We selected a summer and fall Landsat 5, 7, or 8 image for each year in the 1985–2018 study period. We used two seasons of imagery in the BIT process compared to three in the production of the 2016 base as it was determined that doing so sufficiently captured change. The year 2012 was excluded because only scan line corrector-off Landsat 7 imagery were available. Images were selected based on a balance of several competing criteria. First, the date (but more importantly the phenology) of the summer and fall images of a given target year correspond to that of the base year summer and fall images. Second, cloud, shadow, snow cover, and active burns within the analysis extent of the image are minimized. A total of 172 Worldwide Reference System (WRS) path-rows were used for the mapping region covering the western United States. Selected imagery were Landsat Collection 1 Level-2, top of atmosphere products at 30-m resolution in the Albers Equal Area map projection. The Quality Assessment (QA) band for each selected image was used to exclude clouds, shadows, haze, and snow cover.

2.4.2. Ancillary Data

In addition to Landsat imagery, we included several ancillary variables in the regression modelling of fractional components. We obtained topographic data from a 30-m digital elevation model from the National Elevation Dataset, including slope, aspect, and topographic position index. For each selected Landsat image, three spectral indices were calculated, the Soil Adjusted Vegetation Index (SAVI), Normalized Difference Built-up Index (NDBI), and Normalized Difference Water Index (NDWI). Monitoring Trends in Burn Severity (MTBS) and GeoMac burn data were used to produce a suite of burned area products. We obtained 800-m gridded climate data from PRISM (Parameter-elevation Relationships on Independent Slopes Model) [18] to conduct climate analysis and accuracy assessment. The monthly PRISM data were used to calculate the growing season (GS) (April–September) and non-growing season (NGS) (October–March) mean minimum and maximum temperature (TMIN and TMAX, respectively), and the GS and NGS precipitation (PRCP) totals for each year from 1985 to 2018. This process yielded six annual climate variables GSPRCP, NGSPRCP, GSTMAX, NGSTMAX, GSTMIN, and NGSTMIN.

2.4.3. Land Cover Masking

Products were masked using the base land cover mask [14] which included pixels identified as non-rangeland in circa 2016. Some pixels identified as rangeland in the base year were previously (or later) non-rangeland land covers. We developed a series of ERDAS Imagine (Hexagon Geospatial) models to detect and mask out pixels that were inundated with water or were actively farmed (i.e., agriculture) at any point in the time-series. Specifically, to exclude water, we calculated the maximum NDWI through the time-series, adding pixels with 0% topographic slope and a high NDWI value to the land cover mask. Next, to exclude potential agriculture, we highlighted pixels with both a high maximum SAVI value and a high amount of mean variance between summer and fall. We used elevation, land ownership, and ecoregion data to remove commission from pixels identified as agriculture. The water and agriculture masks were added to the base mask and applied throughout the process for each path-row and to the final components.

2.4.4. Normalization

An approach was developed to normalize the individual scenes in the same season to the base year image [11,19]. This approach improves the harmonization of images from different Landsat sensors and reduces differences among scenes within a time period. Linear regression models were fit for the spatial relationship between band x of a given summer/fall target year image and band x of the summer/fall image from the base year. Pixels with cloud, shadow, or snow cover and those that were recently burned were removed from the normalization regressions. This procedure is completed for all six spectral bands and for both seasons. The regression parameters of slope, y-intercept, and correlation were used to adjust the spectral band values of the target year image, so that spatial mean matched that of the base year. If the R² value of a given seasonal pair regression was less than 0.30, the scene was flagged, and either additional cloud masking was applied, or a new image was selected. Normalization serves to minimize spectral differences related to difference among Landsat sensors, phenological variability, sun angle, and atmospheric contamination. The normalized imagery was only used for identifying change pixels, with non-normalized imagery reserved for component predictions. One consequence of normalization is that it can subdue the map depiction of real change occurring on the landscape; however, in most cases the change removed by normalization is overwhelmingly noise.

2.4.5. Change Detection

Each normalized target and base image pair (i.e., summer target and summer base, fall target and fall base) is compared using a Change Vector (CV) approach, finding differences in spectral values by band [3,11,19]. These are summed across bands to produce a change vector magnitude value for each image pair. If a pixel has a CV magnitude beyond a standard deviation threshold in a given target image, for a given NLCD 2016 land cover class, it is flagged as changed [3,11,19]. Two out of two seasons (summer and fall) agreement is needed for a pixel to be flagged as changed between a target and base year. The requirement for two season agreements is critical to remove cloud, shadow, and haze influence missed by the QA bands.

In addition to the CV value approach, we include a secondary change detection method that considers the range of CV values through time at each pixel. We determine how each yearly CV magnitude compares to the range, allowing us to better contextualize the nature of change. We define this metric as Change Fraction (CF), which considers three weighted attributes. First, the seasonal CV magnitude score for a given image in relation to the CV magnitude range per pixel constituted 50% of the total score. This gives a high value to pixels with a high amount of spectral change relative to the total range of variation in the time-series. Second, the CV magnitude range per pixel, relative to the maximum in each path-row, was 25% of the total score. This gives a high value to pixels with a high total variation in the time-series relative to those pixels with lower variation through time. Third, the CV magnitude per pixel relative to the maximum CV magnitude in each path-row, each year, was 25% of the total score. This gives a high value to pixels with a high amount of spectral change relative to the maximum change in that given image.

The total CF score for each season is calculated in addition to the difference in score between summer and fall. Because pixels with a high amount of difference in CF score between seasons in a year are often indicative of haze/atmospheric contamination (in addition to real change), we calculate two versions of the CF, one for defining changed areas (CF_detection), and a second for labeling change (CF_label). Specifically, in this advancement in methodology relative to Rigge et al. [3], we use a more aggressive approach (i.e., detect more change) to define changed areas, CF_detection, given as the maximum CF across seasons. To define areas in which component change should be labeled, CF_label, we again take the maximum of the seasonal CF scores, unless the difference between seasons is large, in which case we take the minimum of the two seasons. In both CF approaches, areas of cloud cover are ignored in the calculation. The benefits of creating two versions of CF are that doing so reduces the impacts of haze/clouds/shadow/other atmospheric issues in the training data (i.e., it creates more pure training data), and it gives us the ability to retain the base predictions in these same areas where the evidence of change is low/suspect. In other words, it creates spectral change “buffer” between pixels used as training and those labeled as change.

We used a convergence of evidence approach to produce the final change area in each target year, considering both the CV and CF_detection. The combination of CV and CF_detection scenarios that best captured known change and minimized false change were selected based on testing multiple combinations in areas known to have changed and not changed. Specifically, change was captured if the CV indicates change and CF_detection > 30 or if the CV does not indicate change and CF_detection > 80. If either scenario is satisfied, then the pixel is flagged as changed in that target year. Determination of change in pixels with cloud cover in one season is based on the CV and CF_detection values from the non-cloudy season, if both seasons are cloudy, no change can be detected in that target year. This procedure produces a time-series of changed/unchanged areas which allow changed areas to be labeled with the appropriate fractional component cover and unchanged areas to be used to train target year predictions.

2.4.6. Training

Within the potential training area, we randomly sampled 200,000 pixels for every target year. To these points, we extract the dependent and independent variables and the SAVI temporal mean per pixel. Independent variables included slope; aspect; topographic position index; non-normalized target year Landsat images; SAVI; NDWI; NDBI; and summer thermal band of each target year. Next, we append all yearly sampled data files per components into a combined spatiotemporal database. We apply two procedures to the training database to improve prediction spatiotemporal dynamic range. First, to improve dynamic range through space, we randomly remove half the observations within one standard deviation of the mean for each component [20]. Next, to improve temporal dynamic range, we compare the ratio of the yearly summer SAVI conditions relative to the long-term mean. Wet years will have values over 100, dry years lower than 100. We again remove half of the observations within one standard deviation of the mean. The overall objective being to stretch out component histograms through space and time. For shrub and sagebrush only, we placed 30,000 points in recently burned areas (those burned in the target year and preceding four years) to which we assigned a dependent variable value of 0% cover. We then appended those observations to train burned areas in the shrub and sagebrush components. The presence of non-zero observations within the main data pool ensures that we do often predict >0% cover of shrub and sagebrush in burned areas, especially in unburned islands.

The training data selection procedure produces a unique changed/unchanged footprint for each year, with the mean component cover value potentially varying over through time within the unchanged (i.e., training data) area. Because training data are pooled through time, the differences in training data available in each year, a potentially confounding effect, are eliminated. Instead, the differences in component cover mapped through time are related only to differences in the spectral conditions in the Landsat images. Pooling training data through time also has the benefit of exposing the RT models to a greater dynamic range of values; additionally sites that were changed in a given year still influence the prediction in that year if they were unchanged in other years. There are two important and unavoidable limitations to our training data approach. First, all training data are derived from the well-trained and validated (see Section 2.3) 2016 base map; however, they are linked to unique spectral conditions in each target year. Second, there is some spectral change below our thresholds in the “unchanged” training area in each target year. This change below the thresholds is small but is critical to calibrate training data to the specific yearly conditions and results in predictions that are more robust to the range of conditions occurring through time.

2.4.7. Component Prediction

We developed Cubist RT models for each component using the spatiotemporal training pool and the full set of independent variables. We developed the Cubist models for the spatiotemporal training dataset but applied the model to data from each target year to produce component cover predictions. A second set of predictions was run for shrub and sagebrush, using training data that included the burned area observations. This burned area prediction for shrub and sagebrush was inserted in pixels that burned in the target year or preceding 4 years, while the standard prediction was kept elsewhere. In pixels with cloud or shadow cover in either season of the target year, we filled the target year prediction with that of the prediction from the preceding target year.

2.4.8. Change Labeling

In pixels meeting the criteria of change, described in Section 2.4.5, except using CF_label data, we apply the target year prediction. If the pixels are identified as non-changed, the base year prediction was maintained. Next, we limited component cover change in three situations. First, if component change between a given target and base year was greater than thresholds described in Table 2 of Rigge et al. [3], we lowered it to the maximum of the limit range. Second, we summed the primary fractional components (bare ground, litter, herbaceous, and shrub) in each pixel, which should sum to near 100%. We identified pixels with a primary component summation over 110% or less than 90% and adjusted the bare ground value so the sum of primary components was 110 and 90%, respectively. For example, if the summation of primary components was 115%, and bare ground was 55%, we lowered the bare ground to 50% so the summation would then be 110%. In cases where the sum was ≥90 and ≤110, no correction to bare ground occurred. Third, we determined the total number of components with change in each target year. If only 1 out of 6 components in a pixel had change, the base year fractional component cover was maintained for all fractional components. The logic is that if one component changes, then other component(s) should change in response.

2.5. Post Processing

We created a series of models to ensure accurate post-burn temporal trajectories in shrub and sagebrush. First, we created a time since most recent fire layer where the value (years) accrues with subsequent years without a burn, resets to zero with every new burn, and accrues again thereafter. We created second-order polynomial regressions of shrub and sagebrush cover limits [3] by time since fire based on expert knowledge which allow more change than in the unburned area. If a shrub or sagebrush cover prediction is greater than this regression line, it is lowered to the limit; in practice this occurred in less than 10% of cases. It is important to note that these limits are somewhat subjective and based on mean recovery conditions across the study area. However, this procedure is very conservative in lowering the rate of recovery in shrub and sagebrush following burns, influencing ~1% of burned area predictions. Shrub and sagebrush cover would recover too quickly in some cases in post-fire conditions based on spectral data alone (i.e., without this procedure), as the burn signature vanishes more quickly than the ecological impact. Such post-fire procedures are not possible for other components as they tend to not follow a recovery pattern as predictable as shrub and sagebrush.

To correct for noise in the time-series we developed a model to detect and remove outliers, defined under two conditions. First, if pixel component cover in target year x is unchanged relative to the base, but the previous and following target years are changed, we replaced the component value in target year x with the mean of the previous and following years. Second, we compared each target year prediction to the prediction of the previous year, to determine if change is beyond broad limits set by expert knowledge. The limits for increase between years was 100% for bare ground, 5% for shrub and sagebrush, 100% for litter, and 15% for herbaceous and annual herbaceous. The limits of decrease between years was −20% for bare ground, −100% for shrub and sagebrush, −10% for litter, and −100% for herbaceous and annual herbaceous. If change was beyond these limits, we replaced the predicted value with the limit. These noise removal procedures were implemented only in pixels that never burned during the time-series.

2.6. Mosaicking

We used the mosaicking functions in ERDAS Imagine to mosaic data from all 172 Landsat path-rows together to produce seamless products. First, automated seamlines were generated in overlap areas with the lowest component value difference. From these seams, feathering was applied to 990-m on either side of the cutline. After the ERDAS processing we did a visual inspection and addressed any seam lines or data voids within the mapping region. The pixel counts were reconciled for all components, for the mask value of 101 and the area outside the mapping region with a value of 102. Finally, we ensured that fractional component relationships were intact, (shrub > sagebrush and herbaceous > annual herbaceous). If not, we forced the relationships to be intact by lowering sagebrush and/or annual herbaceous cover to be equal to shrub and herbaceous cover, respectively.

2.7. Validation

Two robust validation approaches have been previously applied by Rigge et al. [12] and Shi et al. [13] to evaluate the accuracy of the BIT product generation described in the current manuscript. Rigge et al. [12] compared BIT data to predictions of component cover on high-resolution satellite images over multiple locations and years (n = 77) in Wyoming, Nevada, and Montana. Results showed strong temporal correlations with a mean R² of 0.63 and RMSE of 5.47% and robust spatio-temporal correlations with a mean R² of 0.52 and RMSE of 7.89% across components.

The second major validation approach compared BIT results to data from two long term monitoring sites in southwest Wyoming [13]. The field data included 10 years of consistent observation at 126 plots over the period of 2006–2018. The spatial-temporal correlation across all 126 plots showed robust correlations for all components, with an R² (and RMSE) of 0.69 (10.8) for bare ground, 0.40 (7.5) for shrub cover, and a mean of 0.46 (8.1) across components. Moreover, the BIT and field data generally showed similar responses to interannual variation in weather. Changes observed in both the field and BIT data at these sites were typically gradual, within-state, changes that are most difficult to resolve, but were successfully captured. In both BIT validation studies [12,13], the authors found no evidence that accuracy varies as a function of time. In the current study, we additionally validate the BIT data through (1) comparison to weather data and (2) comparison to ancillary data in several case studies (methods described in Section 2.9).

2.8. Trends Analysis

We calculated temporal trends in each component across the time-series using linear models applied to the mosaiced products. For each pixel, we calculated the slope, t-statistic, and standard error using ordinary least squares linear regression. The slope represents the mean percent cover change per year over the times-series, given in units of (cover change%·year⁻¹).

The standard error measures the degree of departure from the linear slope model across all yearly observations in each pixel. Higher values indicate more variance about the mean. The t-statistic reflects the confidence of change in each pixel relative to the null-hypothesis of no-change and is coded negative and positive according to the direction of change. The t-statistic can be used in a Student’s t-test to determine if a pixel has statistically significant change in component cover over the time-series.

2.9. Statistics and Case Study Methods

We evaluated yearly component cover change and annual climate in a sample of 150,000 randomly located points across the study area. From this set of points we calculated the proportion of pixels that had changed at any point in the time-series, the proportion with significant (p < 0.10) increase or decrease in component cover, and the mean component cover and climate variable trends for each level III ecoregion present in the study area. Next, we correlated ecoregion mean climate trends against mean component cover trends.

We demonstrated the local scale accuracy and applicability of our fractional product time-series in three case studies. For the first case study we leveraged a study evaluating experimental prescribed (Rx) burn and control treatments conducted in Hart Mountain National Antelope Refuge, Oregon [21]. In this study, four pastures were burned in fall 1997 and four adjacent pastures were unburned throughout our time-series. We digitized the pasture boundaries using Google Earth and calculated the mean component cover among treated and untreated pastures in each year of our study period. Next, we evaluated component change in riparian areas related to the removal of domestic livestock and wild horses. We obtained a repeat-photo (1990 and 2013) and plot location from Batchelor et al. [22] and compared our component cover values at these same time periods.

In our second case study we evaluate component cover changes related to woody encroachment by honey mesquite (Prosopis glandulosa Torr.), vegetation treatments in the Land Treatment Digital Library [23], and energy development in southeast New Mexico. We demonstrate the departure in component cover related treated compared to untreated sites.

In our third case study, we investigate the impacts of wildfire on component cover values and the trajectories of recovery following burns (n = 4) in southwest Idaho and North Central Nevada. We calculate the mean shrub, herbaceous, and bare ground cover in each fire for all years in the time-series, independently. Next, for clearer visualization, we calculated the departure from the time-series mean component cover in each burned area separately, where positive values indicate higher component cover than the long-term mean in that year.

3. Results

3.1. Overall Trends

Component cover varied substantially though space and time (Figure 3 and Figure 4). Our results showed the net cover of shrub, sagebrush, and litter significantly (p < 0.01) decreased, bare ground and herbaceous cover had no significant change, and annual herbaceous cover significantly (p < 0.05) increased (Figure 4). Temporal trends were somewhat different in pixels that did not burn during the study period, specifically shrub and annual herbaceous cover had no significant change and sagebrush cover increased significantly (p < 0.10). Most pixels changed at least once over the study period (Supplement 1), ranging from 58.8% in sagebrush to 98.9% in bare ground (

Table 1. Summary statistics for linear trends in fractional component cover observed over the 1985–2018 study period. Total percentage of pixels with any change in component cover, significant (p < 0.10) decrease, significant increase, with significant change in either direction, and the ratio of significant decrease pixels to significant increase (where positive values indicate decrease is more common).

Linear Trend	Shrub	Sage.	Herb.	Annual Herb.	Bare Ground	Litter
Changed (%)	98.1	58.8	97.2	98.5	98.9	97.5
Decrease (p < 0.10) (%)	26.7	14.8	22.2	12.6	20.9	24.5
Increase (p < 0.10) (%)	19.2	11.3	18.6	18.6	21.1	16.7
Total (p < 0.10) change (%)	45.9	26.1	40.9	31.2	42.0	41.3
Ratio of decrease to increase	1.39	1.31	1.19	0.68	0.99	1.47
Change (p < 0.10) in burns (%) 1	22.3	23.4	18.6	33.8	17.3	18.5

¹ Proportion of significant change occurring in areas burned during the 1985–2018 study period.

). Significant temporal trends (p < 0.10) occurred in a substantial minority of pixels, ranging from 26.1% in sagebrush to 45.9% in shrub. The ratio of significant change direction varied among components—decrease was predominant in shrub, sagebrush, herbaceous, and litter, and increase more common in annual herbaceous and bare ground (Table 1). Across the study period a slight net increase in GSPRCP (0.75 mm) and slight decrease in NGSPRCP (−5.4 mm) was observed on average across the study area. Temperature variables were more temporally dynamic with a large net increase in NGSTMAX (+1.12 °C), GSTMIN (+0.94 °C), and especially NGSTMIN (+1.54 °C).

Spatial patterns of component trends were highly inter-related; for example, bare ground slope was negatively related to herbaceous slope (r = −0.46, p < 0.05), litter slope (r = −0.57, p < 0.05), and shrub slope (r = −0.41, p < 0.05). Clear patterns in the temporal slope of component cover are evident in burned areas; for example, shrub cover has dramatic increases within the 1988 burned areas in Yellowstone National Park in northwest Wyoming (Figure 3). Conversely, in the more recent large burns in California chaparral (e.g., Mendocino Complex and Camp Fire in 2018), shrub cover decreased, and herbaceous cover and bare ground increased. While change in burned areas is common, it represents a minority of the overall significant change in the study area, with only 17.3% and 22.3% of significant change in the bare ground and shrub components related to burns, respectively (Table 1). Regional patterns in component cover trends are also evident, with increased bare ground in the desert Southwest and Southern Great Plains and decreased bare ground in the Northern Great Plains (Northwestern Great Plains and Northwestern Glaciated Plains) (Figure 3).

3.2. Level III Ecoregion Trends

Component cover (Figure 5) and climate (Figure 6) trends strongly differed by level III ecoregions. Shrub cover decreased considerably in ecoregions with extensive burn histories, namely the Canadian Rockies, Idaho Batholith, Sierra Nevada, and Southern California Mountains. Shrub and sagebrush cover increased in the Northern Great Plains and Wyoming, and Middle Rockies with more limited burn history and increasing GSPRCP and NGSPRCP. Herbaceous cover tended to increase, and bare ground tended to decrease, in many of the same areas. Large increases in bare ground were observed in Arizona, New Mexico, and Texas, especially in the Chihuahuan Desert. This increase in bare ground appears to be related to decreased herbaceous cover in the same regions, itself related to decreased GSPRCP and NGSPRCP (Figure 6). Annual herbaceous cover increased in most ecoregions, except in some portions of California (e.g., Southern California/Northern Baja Coast, and Central California Foothills and Coastal Mountains). In many ecoregions (13/33) annual herbaceous cover increased despite a decrease in herbaceous cover overall.

A strong dichotomy exists in precipitation trends where NGSPRCP decreased in the south and increased in the north (Figure 6). GSPRCP trends showed increases in the northeast and southwest and decrease elsewhere. Temperature variables increased nearly everywhere across the study period, with the largest increases observed in the Southern Rockies, Middle Rockies, Idaho Batholith, and Sierra Nevada ecoregions (Figure 6). Conversely, GSTMIN and GSTMAX decreased in the Northern Great Plains. Ecoregion average climate trends were related to component trends, with the strongest drivers across components being GSTMAX, NGSTMAX, and GSPRCP. Bare ground was by far the most sensitive to climate, especially to GSTMAX and NGSTMAX, which are critical drivers of plant vapor pressure deficit. At an ecoregion scale, we found GSPRCP trends were negatively associated with bare ground trends (r = −0.34), and positively related with shrub cover (r = 0.30). The impact of NGSPRCP tends to be somewhat weaker than GSPRCP, with a weak association with shrub cover (r = 0.06), but stronger influence on herbaceous trends (r = 0.23) than GSPRCP (r = 0.03). Maximum temperature variable trends were more strongly related to component cover trends (r = 0.28) than minimum temperature trends (r = 0.16) across all components.

3.3. Case Studies

We evaluated component change related to fire, vegetation treatments, and grazing management in three case studies to demonstrate the local relevance and applicability of our data in disparate environments and change drivers. The first case study evaluates change in component cover related to a prescribed burn and control treatment study implemented in Hart Mountain National Antelope Refuge in southeast Oregon [21] (Figure 7). Prescribed fire was implemented in 1997 in four units, with four unburned as a control. No other wildfire occurred within these treatments. Within the burn treatments, the severity of burn varied substantially (Figure 7B inset, [21]). Our data show a downward trend in shrub (Figure 7D) and upward trend in herbaceous cover in burned units. Bare ground increased in some burned units, and decreased in others, but changed more than in control units (Figure 7C). Averaged across the control and burn treatments our data show a clear decline in shrub and increase in herbaceous cover immediately following the 1997 fires, and little net change in the control units (Figure 7E). Cover differences among treatments persisted for the remainder of the time-series, with a gradual recovery of shrub (and decline in herbaceous) cover starting ~2005. Within the same management unit, livestock grazing was eliminated by 1991, and wild horse gathers occurred [22]. Our data show decreasing bare ground cover and increasing shrub cover in the riparian areas (linear red feature in upper left to center of Figure 7D) within the case study area. Two inset photos (Figure 7C,D) indicate the change in riparian vegetation between 1990 (grazed by livestock) and 2013 (no grazing). At the corresponding pixel, our data show a 17% reduction in bare ground, and a 22% and 4% increase in herbaceous and shrub cover, respectively.

Our second case study in southeast New Mexico demonstrates component cover change related to the invasion/densification of honey mesquite, the vegetation treatments intended to control its spread, and energy development (Figure 8). The vegetation treatments in the case study area consisted of herbicide application, targeted to eliminate shrubs in general. Our data demonstrate mixed success in these treatments (implemented from 1981 to 2010), with overall decreasing shrub cover and bare ground in treated areas (Figure 8A,B). Patchy areas of increasing shrub cover remain within the treated areas, possibly a result of missed areas in treatment or the re-sprouting of shrubs including Havard/Shinnery oak (Quercus havardii Rydb.). In a plot treated in 2006, bare ground and shrub cover decreased by 1% and 11%, respectively, and herbaceous increased by 6% from 1985 to 2018. In a plot untreated throughout the study period, bare ground, and herbaceous cover both decreased by 3% and shrub cover increased by 11% from 1985 to 2018. The impacts of energy development are also clear in the component cover trends, where oil pad sites and roads have a dramatically declining shrub and increasing bare ground cover.

Our third case study demonstrates the impact of large (2300 to 20,500 ha) burns (n = 4) on component cover, and subsequent recovery in southwest Idaho and north central Nevada (Figure 9). In all fires we document an increase in bare ground cover of 5 to 8% in the 5 to 10 years following fire. Shrub cover decreased by 6 to 10% immediately after fire, to near 0% cover. Recovery of shrubs occurred over the following 5 to 10 years. Herbaceous cover increased following only the Wilson Complex fire in 2005, in other fires it appeared that vegetation cover was lower overall, replaced with higher amounts of bare ground. Our analysis encompassed the entire burned perimeters, including many areas identified in the MTBS burn severity products as unburned areas or lightly burned which may impact the strength of changes related to fire.

3.4. Change Fraction Sensitivity

This multi-level CF approach is possible due to the usage of two seasons of Landsat imagery per target year and is instrumental in reducing the influence of haze and other contamination (Supplement 2). Throughout our data we found no evidence of cloud/shadow/haze contamination. The lower values of CF_label in cloud/haze/shadow cover relative to the CF_detection results in a lower likelihood of false change being labeled in the final products. However, if there is meaningful change in the non-cloudy season, it will be captured. Moreover, the CF is instrumental to our focus of capturing gradual change. Figure 10 shows the sensitivity of the CF to modeled component cover difference in WRS path/row 37/31. Component cover absolute difference between base and target years increase robustly and nearly linearly with CF_label values. These data indicate that on average, a 4.5% and 2.1% cover change is the minimum detection limit for bare ground and shrub cover, respectively. On the other extreme, we expect to capture bare ground and shrub cover change of >6.1% and >2.9%, respectively, in most cases. These data are based on means, our minimum detection limit could vary in practice through space and time. Annual herbaceous cover is limited in this path/row, resulting in a weak relationship.

4. Discussion

Several biotic and abiotic mechanisms (e.g., climate changes, cheatgrass invasion, livestock fencing, industrial/energy development, land use change, and fire regimes) have driven declines in the quality and extent of sagebrush habitat (and rangelands overall) [24]. In many cases component cover change may be associated with change in community state, making restoration of these lands difficult, especially in areas of low ecosystem resilience [24]. Modern management practices and land use decisions are often applied at the scale of individual pastures or ranches, which limits the opportunities related to landscape scale processes like fire, animal migrations, and metapopulation dynamics [25]. Our shrubland component cover estimates provide a wide lensed view of ecosystems, necessary for assessing landscape and ecosystem scale patterns [24], that is also relevant to local decision making (e.g., gradients within a pasture). While both error and bias exist in our predicted component cover values, the latter is remediated when looking across relative changes through time, thus still providing valuable information.

Net changes over the study area resulted in decreased shrub, sagebrush, herbaceous, and litter cover by 14,828, 2238, 3567, and 6575 km², respectively and increased annual herbaceous and bare ground of 3707 and 5036 km², respectively (Figure 3 and Figure 4). Averaging the component change across the study area does cancel out many regional specific patterns (Table 1). Moreover, temporal variation in vegetation tends to be lower than that of spatial variability which widely varies by orders of magnitude across scales of meters to kilometers [8]. The ratio and intensity of change between increase and decrease drives the overall mean, with litter and shrub cover skewed most strongly to decreases and annual herbaceous skewed most strongly to increases (Table 1). Indeed, when looking at pixels with significant (p < 0.10) change in either direction, between 661,313 km² in annual herbaceous and 972,893 km² in shrub experienced a change in component cover. Change is ubiquitous (e.g., Supplement 1); however, most of the change is gradual and/or small. Indeed, given the 33 years in the time-series a gradual change could be statistically significant. Looking at the range in component cover over the time series, well over half of pixels have a range of less than 10%, except for bare ground which has higher mean cover values (

Table 2. Proportion (%) of study area by component cover range throughout the 1985–2018 time-series. Range was calculated as the minimum cover value subtracted from the maximum cover value.

Range (% Cover)	Shrub	Sage.	Herb.	Annual Herb.	Bare Ground	Litter
0	1.46	0.67	2.34	0.53	0.66	1.99
1 to 5	27.87	36.58	25.57	24.58	7.57	44.27
6–10	35.84	39.34	27.43	44.55	21.74	37.91
11–15	16.76	15.08	19.88	17.20	22.08	12.38
16–20	7.45	5.93	12.17	8.23	17.21	2.81
21–25	3.89	1.91	6.56	4.11	12.08	0.51
26–30	2.29	0.43	3.19	2.21	7.81	0.09
31–35	1.77	0.05	1.51	1.10	4.74	0.02
>35	2.66	0.01	1.35	1.22	6.11	0.00
Mean	10.54	7.93	11.18	10.35	17.14	6.59

). Given the prevalence of relatively small cover ranges, paucity of pixels with no change, and that only 14.2% of the study area burned at any point between 1985 and 2018, clearly much of the change observed is occurring outside known disturbances (Table 1). Thus, it is likely that most component cover trends are related to gradual change such as range management practices, change in species composition, long-term climate trends, interannual climate variation, or invasion of annual grasses [10,19].

4.1. Ecoregion Level Trends

Component cover trends varied dramatically among (Figure 5) and within ecoregions (Figure 3). Summarizing trends at an ecoregion level facilitates a gross understanding of ecosystem scale change. One of the major ecoregion-scale patterns was the increase in shrub and sagebrush cover in the Northern Great Plains, Wyoming, and Middle Rockies (Figure 5). In a study combining long-term field observations and a BIT data in southwest Wyoming, Shi et al. [13] found that shrub and sagebrush cover were increasing, ostensibly related to climatic changes in the region, confirming in part our trends in the region. Bromley et al. [26] investigated climate trends over the Northern Great Plains from 1970–2015, finding increased temperatures in winter months, but cooling trends in the spring, countering the prevailing global patterns (as seen in Figure 6 GSTMAX). Moreover, vapor pressure deficit mostly decreased, and precipitation increased (Figure 6) indicating better growing conditions for plants, especially in grasslands where productivity is strongly related to precipitation [27,28]. The anomalous patterns in northern Great Plains, relative to global patterns, may be due to the reduction of summer fallow and expansion of agricultural land cover [26]. The climatic changes in the northern Great Plains are paralleled by significant long-term trends of increasing productivity and water and Nitrogen use efficiency, especially in the driest and warmest portions of this region [28,29]. Reeves et al. [28] also suggest that some increases in productivity in the northern Great Plains may be driven by woody encroachment, which is suggested by the increasing shrub cover observed in the BIT over the same ecoregions (Figure 5).

Severe droughts have occurred during the study period across the West [4], notably in the Chihuahuan Desert [2] and surrounding ecoregions. Droughts in the Chihuahuan Desert have reached such intensity that livestock grazing has been temporarily halted in some locations, or reduced as a hedge against drought, and perennial grass cover declined by 50% during drought years [2]. Perennial grass production decreased 43% in the Chihuahuan Desert of New Mexico from 1967–1992 to 1993–2018, with lower and more variable precipitation and increased temperature in the latter time period [30]. At the same time, Asner et al. [31] found up to a 500% increase in woody cover in untreated areas from 1937 to 1999 in northern Texas, and −75% to +280% woody cover change in areas that had been managed with herbicides, mechanical treatment, or fire. Our data echo these findings in the increase in bare ground cover and decrease in herbaceous cover in the Chihuahuan Desert, and adjacent ecoregions such as the Madrean Archipelago and Sonoran Basin and Range (Figure 5). The ecoregion scale climate trends show decreases in both GSPRCP and NGSPRCP and increases in all seasonal temperature variables in these ecoregions (Figure 6). Our findings support predictions that increased temperature and drought, and lower and more variable precipitation will decrease perennial grass cover in the southwest United States [30].

Some ecoregion-scale patterns were predominately driven by fire history. Ecoregions with extensive fire history earlier in the time-series tend to demonstrate increasing shrub and sagebrush cover, and vice versa for those with fires mostly in the second half of the time-series (e.g., Idaho Batholith, Canadian Rockies, Sierra Nevada, Arizona/New Mexico Mountains, and Central California Foothills and Coastal Mountains) (Figure 5). In these ecoregions, net trends are related to the direct impact of fires and succession thereafter (e.g., [11,32]. Indeed, the proportion of each ecoregion burned over the study period is strongly negatively related to the mean shrub cover slope (r = −0.81, p < 0.05). Similarly, annual herbaceous cover is increasing overall due to increases in its cover in burned areas. The overall declining herbaceous trend countered with the overall increase in annual herbaceous shows that total herbaceous cover is becoming increasingly annual (Figure 4).

4.2. Going Forward

We plan to update the source of imagery to Landsat Analysis Ready Data (ARD) composited at a regional or larger scale for the production of data from 2019, 2020, and onward. We are investigating other methodological improvements to the process, such as incorporating high-resolution training data into the BIT process and use of other machine learning approaches. Our annual herbaceous cover data from 2016 to 2018 was integrated with data from Rangelands Analysis Platform [6] and U.S. Geological Survey Harmonized Landsat Sentinel exotic annual grass maps [33] to produce weighted mean estimates of annual herbaceous cover [34]. This map is a foundational component of the Western Governors Association toolkit for invasive grass management and shows the promise of combining multiple data streams. The current BIT data are available for download [35] and a web application to visualize data and download custom extents is available at www.mrlc.gov/rangeland-viewer/.

5. Conclusions

Our results show that the net study area cover of shrub, sagebrush, herbaceous, and litter decreased, and bare ground and annual herbaceous cover increased over the study period. Most pixels changed at least once over the study period, and significant temporal trends (p < 0.10) occurred in a substantial minority of pixels. While change was ubiquitous, most of the change is gradual and/or small, with well over half of pixels having a range of less than 10% over the time-series. Thus, it is likely that most component cover trends are related to generally slow/gradual processes such as range management practices, change in species composition, long-term climate trends, interannual climate variation, or invasion of annual grasses. Xian et al. [19] reported that such gradual change constituted greater than 60% of total change in southwest Wyoming.

The BIT fractional component time-series provides a wide lensed view of ecosystems, necessary for assessing landscape and ecosystem scale patterns, which is also relevant to local decision making. As with all remote sensing and monitoring datasets error does exist in our products, which needs to be carefully considered. BIT components have been robustly validated through comparison to high resolution imagery [12] and to long-term field observations [13], lending confidence to our data. Moreover, the changes observed in relation to climate variation and management/disturbance (e.g., the case studies presented here) show the BIT data reasonably track these patterns. Our case studies demonstrate but a few of the potential applications of the BIT to track long-term changes. Data facilitate a comprehensive assessment of rangeland condition, evaluation of past management actions, understanding of system variability, and opportunities for long-term planning.