Long-term trajectories of fractional component change in the Northern Great Basin, USA

. The need to monitor change in sagebrush steppe is urgent due to the increasing impacts of climate change, shifting ﬁ re regimes, and management practices on ecosystem health. Remote sensing provides a cost-effective and reliable method for monitoring change through time and attributing changes to drivers. We report an automated method of mapping rangeland fractional component cover over a large portion of the Northern Great Basin, USA, from 1986 to 2016 using a dense Landsat imagery time series. Our method improved upon the traditional change vector method by considering the legacy of change at each pixel. We evaluate cover trends strati ﬁ ed by climate bin and assess spatial and temporal relationships with climate variables. Finally, we statistically evaluate the minimum time density needed to accurately characterize temporal patterns and relationships with climate drivers. Over the 30-yr period, shrub cover declined and bare ground increased. While few pixels had > 10% cover change, a large majority had at least some change. All fractional components had signi ﬁ cant spatial relationships with water year precipitation (WYPRCP), maximum temperature (WYTMAX), and minimum temperature (WYTMIN) in all years. Shrub and sagebrush cover in particular respond positively to warming WYTMIN, resulting from the largest increases in WYTMIN being in the coolest and wettest areas, and respond negatively to warming WYTMAX because the largest increases in WYTMAX are in the warmest and driest areas. The trade-off of lowering temporal density against removing cloud-contaminated years is justi ﬁ ed as temporal density appears to have only a modest impact on trends and climate relationships until n ≤ 6, but multi-year gaps are proportionally more in ﬂ uential. Gradual change analysis is likely to be less sensitive to n than abrupt change. These data can be used to answer critical questions regarding the in ﬂ uence of climate change and the suitability of management practices.


INTRODUCTION
High evaporation rates and variable precipitation make spatially extensive arid and semiarid regions sensitive to climate variability and other perturbations (Hufkens et al. 2008, Mitchell et al. 2015. In the sagebrush steppe of the U.S. Great Basin, sagebrush (Artemisia tridentata spp.) historically occurred in uniform high-density stands (Blaisdell 1953) or in mosaic with grasslands (Bukowski and Baker 2013). A complex interaction of biotic and abiotic change agents including climate change, invasive annual grasses, livestock grazing, shrub control vegetation treatments (physical and chemical), industrial development, and shifting fire regimes has reduced sagebrush cover extent (Wisdom et al. 2005, Miller et al. 2011) and increased fragmentation (Bukowski and Baker 2013). Bromus tectorum, in particular, has dramatically altered ecosystem function (Chambers et al. 2015), causing a positive feedback loop with fire (Boyte and Wylie 2016). Restoration and creation of landscapes resistant to cheatgrass spread and resilient after fire are key management concerns in the sagebrush ecosystem (Chambers et al. 2015, Integrated Rangeland Fire Management Strategy Actionable Science Plan Team 2016).
The need to monitor change in sagebrush steppe is urgent to understand the impacts of change agents, improve carbon and water flux models, and help predict future changes (Anderson andInouye 2001, Olsoy et al. 2016). However, field-based methods are insufficient due to their limited spatial and temporal footprint (Sant et al. 2014) and cost. Remote sensing provides a cost-effective and reliable method of monitoring change, a prerequisite to interpretation of change agents (Vogelmann et al. 2012). Opening the Landsat archive with 40 yr of freely accessible data has driven research in change analysis (Roy et al. 2014). The archive contains more than 950 images in some regions of the western USA (Wulder et al. 2016). Numerous methods have been developed to detect change among a series of multispectral images assuming land surface change results in altered radiance values across one or more spectral bands, including Vegetation Change Tracker (VCT; Huang et al. 2009), Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr; Kennedy et al. 2010), Breaks for Additive Seasonal and Trend (BFAST; Verbeselt et al. 2010), and Continuous Change Detection and Classification (CCDC; Zhu and Woodcock 2014). The increasing number of change detection methods obfuscates the selection of the best for each particular application and the trade-offs among the methods (Tewkesbury et al. 2015). The challenge with all image change detection methods is to detect desired image change while eliminating noise related to phenological change, and image issues of atmospheric condition, geometry, and solar angle. This challenge is compounded when attempting to understand that a change trend is not just a difference between two points in time, but a continual process (Kennedy et al. 2010). Equally challenging is identifying the change agent (wildfire, vegetation treatment, disease, etc.) responsible for altering surface conditions and determining the magnitude and direction of the vegetation change.
While many remote sensing change detection techniques have been developed, very few have been designed for use in rangelands where a relatively weak vegetation signal, complex scattering, and a high amount of bare ground, soil reflectance, and fraction of woody material obfuscate change detection (Hufkens et al. 2008, Sant et al. 2014, Olsoy et al. 2016. Phenological variation in semiarid regions is still a problem in the CCDC approach, as Julian date is not synonymous with phenophase, and thus, phenological variation may be flagged as change. Some have elected to forgo change detection altogether, and simply produce yearly maps of fractional component cover (Jones et al. 2018), an approach that can result in erratic trends due to the same factors that limit change detection.
Gradual change is by far the most pervasive, yet least studied, type of vegetation change (Vogelmann et al. 2012, Shi et al. 2018. While gradual change processes are slow, their cumulative impacts can be substantial (Vogelmann et al. 2012) and provide clues to the impact of various change agents. Homer et al. (2015) demonstrated a method, later refined by Shi et al. (2018), that uses Landsat imagery in a change vector (CV) framework to model gradual changes through time in a suite of fractional rangeland attributes (shrub, sagebrush, herbaceous, bare ground, and litter) measured as 0-100% cover, hereafter fractional components.
Here, we focus on capturing gradual change measured as fractional components using a CV approach. Fractional component maps provide a direct way to quantify surface cover in units both managers and researchers can understand (Jones et al. 2018). Hence, quantifying change using fractional components provides units of change magnitude that are easily understood. Although fractional components are not a direct measure of rangeland health, they function as a useful surrogate to understand longer-term trends and the effectiveness of past management practices and can inform future management (Sant et al. 2014). For example, the Integrated Rangeland Fire Management Strategy Actionable Science Plan Team (2016) identified 37 priority science needs in the sagebrush steppe, many of which can be addressed by a long-term record of fractional vegetation cover. Long-term records of rangeland cover change can provide important insights into ecosystem condition through space and time. Although CV methods have some limitations, such as potentially providing ambiguous results on the direction and magnitude of change (Tewkesbury et al. 2015) and only evaluating change between one pair of images at a time (Verbeselt et al. 2010), we present a modified CV procedure that mitigates these issues by considering the history of change at each pixel across time and reduces the influence of binary thresholds.
Our three primary objectives for this research were to, first, implement a more comprehensive operational approach for quantifying fractional component change from 1986 to 2016 over a much larger area than previously tested (Shi et al. 2018). Our second objective was to evaluate fractional component temporal trends with corresponding climate variables over a large area to assess significant spatial and temporal relationships of interest. Specifically, we evaluated the correlation between fractional components temporal slopes and long-term mean precipitation and temperature. Our hypothesis is that change would be most pronounced in the warmest and driest portions of the study area. Third, while many recent efforts have included many/all available images, few have explicitly evaluated the necessity and impact of doing so. We explore the sensitivity of the temporal sample size (number of years included in analysis) and evaluate its impacts on affecting significant climate and fractional component trends. Our hypothesis is that trends and climate relationships will be largely insensitive to the number of years included but will destabilize at an n < 10. Secondly, we hypothesize that the climate relationships of the herbaceous and annual herbaceous component cover will be more sensitive to n than shrub or sagebrush.

Study site
We analyzed nine Worldwide Reference System (WRS) path-rows in the Northern Great Basin, including parts of California, Nevada, Utah, Oregon, and Idaho (Fig. 1). The total study area comprises 208,093 km 2 , of which 187,013 km 2 was rangeland. Based on Daymet data (Thornton et al. 2014) from 1981 to 2015, the average water year (October to September) precipitation ranges from 40 to 1570 mm, with a mean of 306 mm. Mean water year maximum temperatures ranges from 4°to 21°C, with a mean of 14.9°C, while average minimum temperatures range from À5°to 4°C, with an average of 0.4°C. Elevation is significantly correlated with precipitation (R 2 = 0.41, P < 0.05) and temperature (R 2 = 0.69) in the study area. Elevation averages 1646 m, with a range of 696 to 3293 m. Based on Monitoring Trends in Burn Severity (MTBS), 25.3% of the study area burned at least once from 1984 to 2016.
There are six general vegetation zones within the study area, and we list them in order of decreasing aridity. Shadscale (Atriplex confertifolia)/salt scrub occurs in saline and arid portions of the study area with little perennial grass cover. Desert sagebrush: low-density sagebrush stands with low perennial grass cover. Sagebrush steppe: low-to moderate-density sagebrush with moderate cover of perennial grass, including Sandberg bluegrass (Poa secunda), bluebunch wheatgrass (Pseudoroegneria spicata), thickspike wheatgrass (Elymus lanceolatus), needle-andthread (Hesperostipa comate), bluegrass (Poa spp.), and sedges (Carex spp.; Olsoy et al. 2016). Sagebrush steppe is the most common vegetation zone in the study area. Cheatgrass is most abundant in this zone, often forming extensive monocultures and filling the interspaces between shrubs. Pinyon-juniper woodlands: commonly found in association with mountain mahogany (Cercocarpus ledifolius), bitterbrush (Purshia tridentate; Mitchell et al. 2015), scattered sagebrush, and generally low herbaceous cover. Coniferous forest: consisting of ponderosa pine (Pinus ponderosa), Douglas fir (Pseudotsuga menziesii), with limited amounts of sagebrush and more abundant manzanita (Arctostaphylos spp.) and mountain mahogany. Subalpine to alpine grasslands: predominately meadows but includes stands of mountain and low sagebrush.

Method overview
While our method is fundamentally similar to that of Shi et al. (2018, see flowchart in Fig. 2 of ❖ www.esajournals.org that manuscript), we have made several key improvements to better detect change, address issues with cloud and shadow cover, and developed code to automate the process. All processing was done within individual Landsat images, which were subsequently mosaicked to complete the study area.

Data selection and preparation
Ancillary data.-Our base maps represent the fractional component cover of shrub, sagebrush, bare, herbaceous, annual herbaceous (mostly cheatgrass), and litter in circa 2014 , available at www.mrlc.gov. These base maps were used as the training data source for all other mapped years (target years). Non-rangeland pixels were identified in circa 2014 corresponding with base mapping and were excluded from analysis. Portions of barren land cover from National Land Cover Database (NLCD) 2011 corresponding to large playas were identified by hand and added to the non-rangeland mask. Next, elevation and aspect were used to exclude areas where annual herbaceous cover did not occur or occurred in insignificant amounts, typically above~2100 m, including 7.3% of the mapped area, where annual herbaceous was assigned a cover of 0%. Both the land cover and annual herbaceous mask were static through time so years before or after 2013 may have slight omission or commission of masking (e.g., a parcel of farmland in 2013 was rangeland in 1988).
Topographic variables of aspect (treated as categorical), topographic position index (using a 7 9 7 matrix; Weiss 2001), and slope were derived from a 30-m Digital Elevation Model. We  projection. We selected six bands for analysis: blue (0.435-0.451 lm), green (0.533-0.590 lm), red (0.636-0.673 lm), near-infrared (0.851-0.879 lm), short-wave infrared (SWIR) 1 (1.566-1.651 lm), and SWIR 2 (2.107-2.294 lm). Landsat data were calculated as top of atmosphere reflectance. Target year images were selected on the following criteria: <30% cloud and snow cover and <30 d of difference in date of year between target and base year seasonal image. Our initial analysis showed that matching the date of year between the target and base was important (Singh 1989, Zhu andWoodcock 2014), with more spurious change evident if the dates were more than 2 weeks different due to differences in solar angle, soil moisture, and vegetation phenology (Chen et al. 2003). This finding was balanced against the need for low cloud cover in the context of limited image availability, especially for 2012. Years in which no reasonable image could be found were excluded from analysis on a path-row by path-row basis. Base year images that matched those used in the production of the base maps were preferred, but we deviated from those images in some cases to facilitate better date matching in target years.
We elected to use only two seasons (summer and fall) of imagery for each year in an effort to include more target years. Spring was determined to be least critical as it tends to be the cloudiest and has the highest snow cover in our study area, which would result in the exclusion of more target years. Further, Olsoy et al. (2016) and Mitchell et al. (2015) suggest that the influence of herbaceous cover on the shrub cover models is lower in late summer and early fall, following senescence of herbaceous plants; importantly, tall grasses were found to interfere with shrub cover predictions.
Image processing.-Cloud and shadows pose significant challenges to accurate change detection and mapping, but automated masking products (Zhu et al. 2015) often have considerable omission error. To remedy this issue, we chose to augment automated masking products with manual identification and exclusion of clouds, cloud shadows, and snow cover from all images. These areas were treated as missing data and ignored in further analysis. To further enhance vegetation discrimination beyond individual spectral bands, we calculated three indices for each target year image: Soil Adjusted Vegetation Index (SAVI; Huete 1988), Normalized Difference Built-up Index (NDBI;Xu 2008), and Normalized Difference Water Index (NDWI; Gao 1996). Calculated as where NIR = near-infrared (0.851-0.879 lm) and SWIR band 1 = short-wave infrared (1.566-1.651 lm).

Change detection
Change vector analysis.-We conducted a spectral band-by-band regression comparison between each target image and the corresponding seasonal base Landsat 8 image in areas without cloud cover or recent burns and applied these regressions to normalize target images. Normalization serves to minimize differences related to phenological variability, sun angle, and atmospheric contamination. It is important to note that normalized imagery was used only for change detection.
Next, we compared each normalized target and base image pair (e.g., summer 1988 and summer 2014) using a CV approach ), producing the difference in spectral values by band. These are summed across bands to produce a change vector magnitude value for each image pair. Magnitude values for each pixel are compared relative to the target year image standard deviation within each NLCD 2011 land cover (Table 1), set to be 25% more sensitive to change than previous efforts (Shi et al. 2018). If a pixel has a change vector magnitude beyond this standard deviation threshold in a given target year image, it is flagged as changed. Previously, Shi et al. (2018) required 2-out of 3-season agreement in change. In the current effort, we chose to use 2-out of 2-season agreement for a pixel to be flagged as changed between target and base year. This more conservative approach as compared to 1-out of 2-season change agreement increases the amount of available training data, reduces cloud contamination, provides less commission in change, and alleviates some concern related to date matching between target and base images.
Change fraction.-The CV approach has been reported to be quite robust (Xian et al. 2009); however, it has several limitations. (1) It does not consider the history of change at each pixel, (2) it is binary, producing strongly juxtaposed maps of changed and unchanged areas (Tewkesbury et al. 2015), and (3) it lacks flexibility in its application, often based on subjective binary thresholds (Singh 1989, Chen et al. 2003. To overcome these challenges, we use an approach that considers the range of CV values through time at each pixel. Then, we determine how each yearly CV magnitude compares to the range, allowing us to define the degree of confidence in that change. We define this metric as the change fraction (CF).
The first step in calculating CF is to clean the seasonal CV magnitude stacks by excluding clouds and outliers. CV magnitude values are capped at a specified threshold (CV = 2000), deemed to remove spurious data. From the cleaned CV magnitude stacks, we determine the 10th and 90th percentile values across each season through the time series, respectively. Finally, we take the 90th percentile value minus the 10th percentile value to give the seasonal CV magnitude range. The CV magnitude range represents the extent of variability through time at each pixel: Those pixels with stable conditions will have a low value, while those with disturbance will have a high value.
The CF considers three weighted attributes in change detection. First, the CF considers the seasonal CV magnitude score for a given image in relation to the CV magnitude range per pixel, which constituted 50% of the total CF 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 CF considers the CV magnitude range per pixel, relative to the maximum in each path-row, which constituted 25% of the total CF 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 CF considers the CV magnitude per pixel relative to the maximum CV magnitude in each path-row, each year, which constituted 25% of the total CF score. This gives a high value to pixels with a high amount of spectral change in a given image relative to the maximum change in that given image. The weighting proportion of each of the three attributes was determined based on evaluation of multiple locations known to have changed and others with no change evident in the imagery. We tested various iterations of weighting proportions and selected the scenario that best captured the known change and minimized the detection of false change. The total score for each season is calculated, and the maximum of the two seasons is given as the final CF. Values of CF range from 0 (high confidence unchanged) to 100 (high confidence changed).
Final change area. -The final change area in each target year considered both the CV and CF using a convergence of evidence approach, wherein two scenarios can indicate change. Different scenarios of CV and CF thresholds were tested on multiple areas that either are known to have changed or not changed. We selected the combination of CV and CF scenarios that best captured known change and minimized false change. Specifically, there were two ways in which a pixel could be flagged as changed: (1) CV indicates change and CF > 34 or (2) CV does not indicate change and CF > 80.
If either scenario is satisfied, then the pixel is flagged as change in that target year. This procedure produces a time series of changed/unchanged areas that allow changed areas to be labeled with the appropriate fractional component cover and unchanged areas to be used to train target year predictions.

Training data
Target year fractional component training is obtained from the base year fractional components. First, we identify the potential training area for each target year, in which the base year fractional cover values are appropriate to train the target year. This is defined as pixels that are unchanged from the base in a given target year; those that do not meet either criteria of change in Eq. 1. Those pixels that are unchanged relative to the base year in a given target year are defined as the unique training area. Pixels that are unchanged relative to the base year through all target years are defined as the common training area. Second, we randomly select pixels from within the unique and common potential training area for each target year. The pixels that are sampled from these potential areas are defined as the unique and common training, respectively. Third, for shrub and sagebrush only, we add a training pool for recently burned areas. These training pools are combined in each target year to serve as training for fractional component models.
Unique training.-Unique training constitutes the majority of the data used to train our target year fractional cover components, consisting of 210,000 randomly selected pixels within areas that are spectrally unchanged between a given target year relative to the base year. The unique training area is based on pixels that do not meet either change criteria in Eq. 1. From this area, we also removed pixels that burned in each target year and the preceding 4 yr and clouds and shadows in the base and target years. Finally, a clump and sieve approach was used to exclude clusters of <5 pixels that met all other criteria in order to favor larger patches of unchanged area in the training data pool for each target year. In larger patches (≥5 pixels) of unchanged areas, we can be more confident that the pixels are indeed unchanged, and therefore, the base year predictions are appropriate to use as training.
Common training.-To make target year fractional cover components more spatially and temporally consistent, we also added a stricter training region from which the same 90,000 training points are used for all predictions throughout the time series. We called this region the common area, which was produced by intersecting the yearly unchanged area, and further refined by excluding all 1984-2016 burns and clouds in any target or base year. This defines the area in which we are confident that minimal component cover change has occurred during the study period. We did allow pixels with one year of change to be included in the common area (assuming this change was spurious), provided that it was not in the last target year, in burned area, or under cloud cover.
Burned area training.-We obtained a third pool of training data to ensure accurate predictions of shrub and sagebrush cover in burned areas. Specifically, we selected 30,000 random points within the combined extent of moderate and severe burns of each target year and the preceding 4 yr and assigned a training data value of 0% shrub and sagebrush cover to each. These training points were appended to the unique and common area points for each target year. Research in the Great Basin has shown that 0% shrub and sagebrush cover is a reasonable assumption in these recent burns (Hanna and Fulgham 2015). Other fractional components are based on only the unique and common training pool, since considerable variability among fire outcomes exists in those fractional components.

Target year fractional component estimates
Training data points from the unique and common area were pooled for each fractional component. Next, we imported training data into a statistical procedure that randomly thins the most common values based on histogram bins (Wylie et al. 2018). Following these procedures, approximately 90,000 points were available for each fractional component. These data were used as the dependent variable in Cubist software (RuleQuest Research 2008), which is based on regression tree algorithms, to create regression models (Xian et al. 2013) with slope, aspect, topographic position index, non-normalized target year Landsat images, spectral indices SAVI, NDWI, and NDBI, and summer thermal band used as independent variables. A second estimate was run for burned areas in shrub and sagebrush, where we appended the burned area training to the training pool and ran the Cubist estimate as above. Next, for shrub and sagebrush, we inserted the burned area estimate into pixels where a fire occurred that year or the preceding 4 yr and kept the primary estimate elsewhere. Finally, we filled in cloud and shadow areas in either season of the target year with the estimate from the following target year, and if ❖ www.esajournals.org that year was also cloudy, then we used the preceding target year.

Change labeling
Following completion of the target year estimates and change detection metrics, we labeled change in change pixels considering both the CV and CF and kept the base year fractional component cover in unchanged pixels. The CV was run through a 3 9 3 focal window to filter changed pixels, keeping the change only if >4 of the adjacent 9 pixels were flagged as change. Final change indication was based on the formula in Eq. 1. If a pixel was flagged as change, it was labeled with the target year estimated percent cover, if not, then the base year percent cover was kept. Final labeled change was calculated as

Post-processing
Following the initial labeling of change, we developed several post-processing steps to remove potential noise and ensure consistent relationships among predictions. We want our products to represent a combination of empirical data and logic-based rules, and therefore, we implemented several cleanup procedures. (1) Based on expert knowledge, we set limits of fractional component change based on the number of years between target and base (Table 2). If component change in any given target year was beyond these limits, change was lowered to the limit, otherwise we kept the change.
(2) Ideally, we want the sum of the primary fractional components to equal 100% in every target year, but we tolerated a range between 90% and 110%. We identified pixels with a summation beyond these tolerances for each target year and adjusted the bare ground value to force the summation to be 90-110%. (3) Next, we identified the local concentration of AE1% change, likely to be noise, using a 5 9 5 focal moving window. In pixels with a high amount of AE1% change (>80% of pixels in the focal window) or with a low amount (<20% of pixels), the base fractional component cover was inserted, while the target year estimate was maintained in the remaining area. (4) We determined the total number of fractional components with change in each target year. If only 1 out of 6 fractional components had change, the base year fractional component cover was maintained for all fractional components. We presumed that if one fractional component changes, then other component(s) should change in response. (5) Finally, we insured that fractional component relationships are intact (shrub > sagebrush, herbaceous > annual herbaceous).
Next, 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, which resets to zero with every new burn and accrues with subsequent years without a burn. Pixels are treated as unburned until the time of first fire. We created second-order polynomial regressions of shrub and sagebrush cover limits by time since fire based on the literature and expert knowledge Notes: n is the number of years between the target and base year; y is the number of years since most recent fire.
( Table 2), which allow more change than in the unburned area. If the target year estimates for shrub or sagebrush were above the limit at a given time, they were lowered to the limit. Target year estimates below the limit were kept. Finally, we insured that successive years of estimates in burned areas had values greater than or equal to the previous year.

Mosaicking and trends analysis
To produce seamless fractional component mosaics across path-rows in our study area, we first selected years to be included in the mosaics. Data missing from path-rows in each year were filled with the closest available subsequent date. To minimize the amount of temporal discrepancy, the decision to include/exclude years in the data stack was a balance between temporal density and data purity. We elected to use only years in which greater than 75% of path-rows had data available.
Mosaics of fractional component data for each selected year (n = 19) were produced using ERDAS Imagine (https://www.hexagongeospatia l.com/products/power-portfolio/erdas-imagine) mosaic express. Overlap conflicts in fractional component values were resolved by calculating an automatic most nadir seamline and feathering the overlapping data from that seamline for a distance of 990 m. Feathering performs a distance-weighted averaging function that considers distance from the generated seamline in its calculation. Of the entire study area, 3992 km 2 or 1.9% was potentially averaged in each year. Cubic convolution resampling was utilized in the mosaicking process. Filling data in missing path-rows can result in estimates that are temporally out of sequence with burn history. To correct this, we ran the fire post-processing model on shrub and sagebrush once more to ensure reasonable post-fire succession. Finally, each fractional component mosaic and selected climate datasets were run through a linear trends analysis model to produce the temporal slope, t-score, standard error, and P-value.

Data analysis
Validation.-Fractional component data were validated using three methods: independent validation, Bureau of Land Management Assessment, Inventory, and Monitoring (AIM) data, and pseudo-invariant pixels. Independent field validation data were collected concurrent with our base mapping in circa 2014 . At the independent points widely dispersed across our study area (n = 156), we calculated the coefficient of determination and root mean square error (RMSE) between field-measured data and modeled fractional components in the base year. AIM data collected in 2011, 2013, 2014, and 2016, corresponding with our mapping years, were pooled (n = 1794). We calculated the pooled (spatiotemporal) coefficient of determination and RMSE between AIM data and modeled fractional components. Using Google Earth, we identified large-diameter (>60 m) pseudo-invariant surfaces consisting of talus slopes, lava flows, and shear rock faces, which should bẽ 100% bare ground. Playas were excluded as water-level and moisture fluctuations can drive spurious predictions. At these locations (n = 80), we selected the Landsat pixel nearest the centroid of the surface to which we extracted the time-series fractional cover value and temporal slope. We then calculated the mean fractional component cover value and tested whether the mean temporal slope differed significantly from zero.
Spatial trends.-From a random sample of 88,035 observations, including 65,770 points in unburned areas, we calculated the average yearly value for each fractional component and weather variable and calculated the slope of the temporal trend. Next, we converted the fractional component slope into total net change by multiplying by the total number of years in the study (n = 30), binned the change into 11 increments, and calculated the proportional area of each fractional component within each bin. Finally, in each observation, we took the maximum value of each fractional component through time minus the minimum value to give the range and we calculated the histograms of range. For example, a pixel with a maximum shrub cover value of 50% and minimum of 40% through the time series would have a range of 10%.
Climate influence on fractional component slope.-To evaluate the degree to which climate influenced fractional component slopes, we divided the 1986-2016 mean climate data into bins and then averaged (total n = 65,770) the temporal slope of all fractional components and WYPRCP, WYTMAX, and WYTMIN in each bin in unburned areas. WYPRCP was divided into 50mm increments, WYMAX into 1°C increments, and WYTMIN into 0.5°C increments. Bins that occupied <0.1% of the analyzed area were excluded. The correlation (r) between fractional component slope and climate bin values was then evaluated.
Sensitivity to number of years (n) and temporal gaps.-We calculated the temporal slope and coefficient of variation (CV) for each fractional component including all years (n = 19), then randomly withheld an increasing number of years until 4 yr remained. This procedure was repeated over 10 iterations, and the average slope and CV among the iterations were calculated and plotted against the n of the regression. Our hypothesis is that slope values will destabilize and CV will increase as more years are withheld. The point at which slope and CV increase in variance represents the minimum n required over our study period to capture salient trends, which may vary among fractional components. A similar exercise was undertaken to evaluate the impact of temporal gaps of 2-4 yr on fractional component slope and standard deviation.
Using similar logic, we evaluated the relationship between each fractional component and selected climate driver (WYPRCP, WYTMIN, and WYTMAX) with randomly withheld years, decreasing until 3 yr remained, repeated over 10 iterations. The average fractional component to climate correlation (r) and CV were calculated among the 10 iterations.

Validation
Our base year-independent validation reveals that modeled fractional component values used as training are moderately related to field-measured data, with an R 2 ranging from 0.19 to 0.43 (Table 3). AIM data validation showed similar, but somewhat weaker correlations with modeled fractional component values (Table 3). This weaker relationship may be due to differences in observation methods and spatial discrepancies. Similar to Shi et al. (2018), we found that no fractional components in pseudo-invariant pixels had a mean temporal trend significantly different from zero, and that bare ground was substantial, averaging 81.5% (Table 3). Moreover, Shi et al. (2018) found that fractional components responded as expected to burns, with greater cover change in severe burn intensity as opposed to low and moderate; bare ground and litter increased and all other fractional components decreased in the year of burn relative to the prior year.

Spatial trends
Both abrupt and gradual change was evident in the fractional component predictions through time (Fig. 2). Abrupt and gradual changes were largely confined to the burned and unburned areas, respectively. Most change was associated Notes: Independent validation results demonstrate the coefficient of determination and RMSE of field-measured data in the circa 2014 base year against the base year prediction. Assessment, Inventory, and Monitoring (AIM) data demonstrate the coefficient of determination and RMSE of pooled data from 2011, 2013, 2014, and 2016 against the corresponding target year modeled component cover. Pseudo-invariant results show the mean fractional cover value. The P-value is for a t-test of significant difference from zero temporal trend, fractional component. with disturbance and recovery from burns before or during the study period; however, a significant amount of change occurred in areas with no known disturbance. Over the entire study area shrub (À0.04% per decade), sagebrush (À0.02% per decade), litter (À0.02% per decade), herbaceous (À0.01% per decade), and annual herbaceous cover (À0.006% per decade) decreased and bare ground (0.01% per decade) increased (Fig. 3). Precipitation was fluctuated substantially throughout the study period, with a slight increase in water year precipitation evident (Fig. 3). The net trend of TMAX was~0, while TMIN markedly increased (+0.47°C per decade).
Shrub and sagebrush trends in unburned areas diverged from the study area average with increases of 0.03% per decade and 0.01% per decade, respectively. The increase in bare ground (0.03% per decade) in unburned areas was less dramatic, and the decrease in herbaceous cover more dramatic (0.03% per decade) than the entire study area. Not only did the trends diverge between unburned and all areas, but the fractional component values also differed. Shrub, sagebrush, and bare ground were higher in unburned areas throughout the study period, while herbaceous, annual herbaceous, and litter cover were lower in unburned areas. Fire occurred in 25.3% of the study area through time, yet 50% of the significant (P < 0.05) declines in shrub cover, 36% of the significant increase in bare ground, and 43% of the significant increase in herbaceous cover occurs in this area. Many fractional component trends are inter-related, the most ecologically significant relationships are (1) the negative association between bare ground and herbaceous cover (R 2 = 0.46, P < 0.05), (2) negative association between bare ground and shrub cover (R 2 = 0.38, P < 0.05), and (3) positive association between shrub and herbaceous cover (R 2 = 0.04, P > 0.05). In unburned areas, the bare ground to herbaceous cover relationship is even stronger (R 2 = 0.75, P < 0.05), but the bare ground to shrub cover association is not significant, and the shrub and herbaceous cover association is reversed to a negative direction (R 2 = 0.16, P < 0.05).
Overall, the percentage of area with a net trend of less than AE1% cover change varied from 31% in bare ground to 52% in annual herbaceous cover; while in unburned areas, 39% of the pixels had a net trend of <AE1% in bare ground and 61% in annual herbaceous (Fig. 4). In areas that burned at least once from 1984 to 2016, only 10% of bare ground and 25% of annual herbaceous pixels had a net trend of <AE1%, respectively. The histograms of cover change in shrub and sagebrush cover change are skewed to the right (i.e., more decrease than increase), especially in burned areas. Put into context, 87% and 91% of shrub and sagebrush pixels, respectively, had a net change of within AE5% over the study period. Herbaceous cover over the study area has slightly more area of decrease than increase, but in burned areas, the opposite was the case.
The net change over the study period obscures some of the interannual variability in fractional component cover. When considering this variability, the proportion of pixels with 0% cover change through time decreases (Table 4). Bare ground, for example, has only 9.9% pixels with a cover range through time equal to 0 (Table 4). This value area is considerably smaller than the overall 31.5% of bare ground pixels with a net trend of <AE1%. Average range varies from 4.9% for annual herbaceous cover to 11.2% for bare ground. Maximum range with a fractional component is lowest in sagebrush at 48% and highest with bare ground at 88%.  (Fig. 5A brown line). Interestingly, herbaceous cover slope decreases significantly with WYPRCP, while annual herbaceous slope significantly increases (Fig. 5A light green and dark green lines, respectively). Sagebrush cover slope has a non-linear relationship with WYPRCP, peaking at about 550 mm. Shrub cover slope significantly decreases with WYTMAX and WYTMIN ( Fig. 5C and E blue line), conversely bare ground slope significantly increases ( Fig. 5C and E, brown line). Herbaceous cover slope significantly increases with WYTMAX, but not WYTMIN ( Fig. 5C and E, light green line, respectively). Annual herbaceous cover slope decreases non-significantly with WYTMAX and WYTMIN ( Fig. 5C and E green line). Finally, sagebrush slope decreases significantly with both WYTMAX and WYTMIN ( Fig. 5C and E black line). Shrub and herbaceous cover differ in spatial response to WYPRCP, WYTMAX, and WYTMIN (Fig. 5A, C, and E blue and light green lines, respectively).

Slope response to temporal gaps in data
The random temporal holdout analysis confirmed our hypothesis that slope values and CV destabilize with more temporal periods withheld (Fig. 6). With one-year periods withheld (Fig. 6A), average slope values remained relatively consistent until more than 15 yr were removed, after which slopes became more erratic. Slope CV among iterations, however, continually increased with greater date withholding then abruptly increased at around 15 yr removed. Comparable results were found in withholding two-, three-, and four-year gaps (Fig. 6B-D), with average slope values tending to destabilize and CV increase with less observations. From a perspective of total years withheld, larger blocks appear to be more problematic. The slope CV is generally larger for the three-and four-year blocks for a similar total year withholding than in one-or two-year blocks. Moreover, the stability of slope values with equivalent year withholdings tends to be less with larger blocks.

Climate relationship response to temporal gaps in data
Similar to temporal slopes, fractional component relationships with climate variables are impacted by the temporal density of observations (Fig. 7). Fractional component relationships with WYPRCP and WYTMAX tend to destabilize with more observations withheld; this pattern is most evident in relationships with WYTMIN, particularly in sagebrush and shrub cover. Climate relationship CV increased non-linearly with years withheld, especially those with WYTMAX and WYTMIN. Herbaceous and annual herbaceous cover tended to have the greatest CV in correlation with WYPRCP among iterations at any given level of year withholding, bare ground had the greatest CV with WYTMIN, and shrub and sagebrush cover had the greatest CV in relationship with WYTMAX.

Temporal trends
Across the study area, shrub cover significantly declined at a rate of À0.04% per decade and bare ground increased at 0.01% per decade (Fig. 3). This finding is consistent with, but less dramatic, than the trends reported by Shi et al. (2018) in a subset of the current study area who found shrub declined at À0.14% per decade and bare ground increased 0.05% per decade. While Shi et al. (2018) reported little trend in litter cover, we report a significant decrease of À0.02% per decade. Many areas are still susceptible to initial cheatgrass invasion and densification (Miller et al. 2011); however, we show a slight downward trend in annual herbaceous cover, possibly due to the fact that 50-60% of the historic sagebrush steppe range either possessed an understory of exotic annual grass (i.e., cheatgrass) or has already been converted to an annual grassland (West 2000). Indeed,  reported that most significant trends in cheatgrass performance from 2000 to 2010 in the Northern Great Basin were negative. Moreover, when looking at only burned areas, annual herbaceous cover increased at a rate of 0.03% per decade, as expected in burned areas (Hanna and Fulgham 2015). Unburned area differences are consistent with ecological expectations, in that shrub and sagebrush cover is increasing, driven largely by secondary succession (Hanna and Fulgham 2015), while decreases in the study area overall are mostly driven by burns , Shi et al. 2018, due to the slow recovery of sagebrush (Baker 2011, Bukowski andBaker 2013). Rigge et al. (2019) demonstrated that postfire recovery of shrubs in the Great Basin required 10-15 yr, using a space-for-time approach with the same base year data and sharing the same geography as the current study. Concurrently, the increase in bare ground is less in unburned areas as compared to burned areas. The greater decline of herbaceous cover in the unburned areas compared to burned areas likely results from the replacement of shrubs with herbaceous cover in burns (Blaisdell 1953, Shi et al. 2018. Additionally, the significantly greater decline in WYPRCP in unburned areas (À8.5 mm per decade) as compared to burned areas (À4.4 mm per decade) creates less favorable bioclimatic conditions for herbaceous cover.
The negative temporal association of bare ground with herbaceous and shrub cover reveals a strong negative relationship wherein vegetative cover responds positively to WYPRCP (bare ground R 2 = 0.34, P < 0.05, herbaceous cover R 2 = 0.26, P < 0.05), and negatively to WYTMAX and WYTMIN, thereby modulating bare ground.  The vacant niches in burned areas result in a weak positive association between shrub and herbaceous cover, suggesting only a limited amount of competition. Conversely, in unburned areas competition increases, resulting in a negative relationship between shrub and herbaceous cover (R 2 = 0.16, P = 0.09).
Similar to previous findings (Vogelmann et al. 2012, Shi et al. 2018, we find limited change when averaged across the study area, due to diverse landscapes and a large area . Vegetation fractional components and bare ground have a net decrease and increase, respectively (Fig. 4). All fractional components had >80% of pixels within AE5% net change over the study period. However, when looking at temporal range in cover (Table 4), considerably more change is evident and only 26-62% of pixels have AE5% range. The implications of this difference are twofold: It implies there may be some pixels (varying from 29% to 53% by fractional component) in which a linear model does not fit the temporal trend and that may be better assessed with a non-linear or break-point method (Vogelmann et al. 2012, Forkel et al. 2013. Secondly, it suggests that real or anomalous temporal outliers exist in the data that result in a higher temporal range of cover values than simply extrapolating the slope across the time period (e.g., a pixel that had a maximum bare ground cover of 90%, and minimum of 60%, thus a range of 30%, but a slope of only 5.0% per decade).

Fractional component slopes by climate
The influence of climate is complex and sitespecific, and with differing impacts by species and life stage, seedlings often being most sensitive (Schlaepfer et al. 2014). Indeed, the ecohydrological niche of sagebrush ecosystems was most strongly determined by weather, with edaphic factors of lesser importance (Schlaepfer et al. 2012a). We show that fractional component cover slopes were significantly related to 1986-2016 climate in linear and non-linear relationships (Fig. 5). While sagebrush tolerates a range of conditions (Still and Richardson 2015), its cover tended to decrease the most in the warmest (Fig. 5B-C), mostly southern, portions of the study area, where ecohydrological and climatebased models predicted sagebrush distribution to decrease (Schlaepfer et al. 2012b). Sagebrush and shrub cover slope increased with WYPRCP at values below 550 mm (Fig. 5A), with greater increases in cover in wetter sites that also tend to be cooler. Above 550 mm of WYPRCP, sagebrush and shrub cover slope response to WYPRCP diverges, with shrub slope following the expected increasing pattern, and sagebrush slope unexpectedly decreasing. This difference in patterns is ostensibly related to the pattern of WYPRCP slopes, specifically that WYPRCP is decreasing at the slowest rate around 550 mm, above which the rate of decline sharpens (Fig. 5B). Moreover, the rate of WYTMIN increase is strongly related to annual average WYPRCP, where wetter areas are most rapidly warming (Fig. 5B). This suggests that sagebrush responds differently to warming differently in wet locations than does shrub cover, as indicated by its more gradual increase in sites with an average WYTMAX of less than~14°C (Fig. 5C). Another possible explanation is model confusion between tree and shrub cover, but not sagebrush cover, in the wettest 5% of the study area.
WYTMAX and WYPRCP slope is positively related to average WYTMIN and WYTMIN slope is negatively related to average WYTMIN (Fig. 5D). Fig. 5C shows that shrub and sagebrush cover increase in sites with the cooler average WYMIN below~17°C and decrease in warmer sites. Cool WYTMIN sites have the greatest increase in WYTMIN, smallest increase in WYTMAX, and greatest decrease in WYPRCP (Fig. 5D). This pattern suggests that increased WYTMIN is more influential than decreased WYTMAX in cool sites for shrub and sagebrush. Secondly, it demonstrates that warming in cool sites creates a positive cover response, and vice versa in warm sites, similar to Kleinhesselink and Adler (2018) who reported that colder portions of the sagebrush range have a positive response in cover in warmer years, while warmer areas have a negative response. The positive influence of warmth at cold sites is more significant than its negative influence at hot sites (Fig. 5C), again concurring with Kleinhesselink and Adler (2018). In summary, shrub and sagebrush cover respond positively to warming WYTMIN in our study area, resulting from the largest increases in WYTMIN being in the coolest and wettest areas, and negatively to warming WYTMAX since the largest increases in WYT-MAX are in the warmest and driest areas. Sagebrush in sites with cool and wet climates within the Great Basin are likely to be most persistent with future climate changes (Creutzberg et al. 2015) and disturbance (Wisdom et al. 2005).
Patterns of bare ground slope by climate bin are inverted in most cases from those of shrub and sagebrush (Fig. 5A, C, D). While bare ground slope is negatively related to WYPRCP bins, WYPRCP itself does not appear to be the driver. Rather, similar but opposite to shrub cover, bare ground responds negatively to warming WYTMIN, resulting from the largest increases in WYTMIN being in the coolest and wettest areas, and positively to warming WYT-MAX since the largest increases in WYTMAX are in the warmest and driest areas. Herbaceous cover slope patterns are also opposite to shrub cover, but to a lesser degree than in bare ground, with WYPRCP slope being more critical to herbaceous cover than other fractional components. Herbaceous cover slope is likely a displacement response to the direction of shrub cover slope at a site. This analysis partly confirms our hypothesis but underscores the site-specific nature of climate relationships (Schlaepfer et al. 2014). Moreover, it highlights the importance of considering the direction and magnitude of change for multiple climate variables.

Temporal density
The time period assessed in analysis can affect the strength and direction of temporal trends, as can the number of observations per year (Wulder et al. 2016) or in the current study, the number of years observed. With 36 yr available in the Landsat (5, 7, 8) archive, some researchers (Jones et al. 2018) elect to produce maps and corresponding trends over all available years (Vogelmann et al. 2012). Others have used a sparse temporal density, especially prior to opening the Landsat archive for free download, such as Walston et al. (2009) who used only five Landsat images to evaluate trends in the distribution of big sagebrush in Wyoming over a 22-yr period. Our analysis of temporal density influence of fractional component slopes shows that net results vary little until n ≤ 6 ( Fig. 6A), confirming our hypothesis that trends and climate relationships will be largely insensitive to the number of years included but will destabilize at n < 10. On the other hand, temporal slope CV increases continually with decreasing n (Fig. 6B). Multi-year temporal gaps tend to disrupt results and increase CV more per total years withheld than with single years (Fig 6B-D).
Reduction of n for increased processing speed and/or cost reduction then appears to be reasonable for analysis of gradual change, especially those with a consistent linear change through time. The analysis of withholding observations from climate relationships (Fig. 7) underscores this point, as the impact of missing dates on climate relationships is positively correlated with the temporal variability of the climate metric. Temporal CV of WYTMIN is highest, followed by WYTMAX and WYPRCP. The climate relationship sensitivity to n withheld (i.e., CV) is highest with WYTMIN across all fractional components; moreover, those with the greatest temporal variation are generally more sensitive overall. In other words, it appears that fractional component relationship with climate variables sensitivity to n is a function of both the temporal variability of the climate variable and fractional component in question.
Gaps in data of >2 yr should be avoided. In our case, missing years result from the lack of phenologically consistent or cloud-free imagery between the target and base year. The trade-off would then be temporal completeness versus data quality and consistency; we deemed the latter to be more critical. The results show this prioritization was justified, and moreover that n could be further reduced with no meaningful change in results. Further, our results indicate we can be confident that our temporal slopes (Fig. 3) and climate relationships (Fig. 5) contain a sufficient density of observations to report accurate trends. For this analysis, we randomly withheld years; however, certain years are likely to be more influential to overall trends and relationships, particularly those near the beginning and end of the time series. Analysis of abrupt change will require greater temporal density or non-linear regression (Vogelmann et al. 2012).

Methodology
Time-series analysis of change, in particular gradual change, is subject to issues such as atmospheric influence, sensor consistency, change detection performance, and lack of validation data (Vogelmann et al. 2016). Our fractional component predictions were successful in providing a consistent, objective, and broad view ) of the vegetation conditions in the Northern Great Basin through space and time. Shi et al. (2018) demonstrated that our cover predictions were consistent in pseudoinvariant pixels. Moreover, Shi et al. (2018) found that products responded as expected to burns with greater cover change in severe burn intensity as opposed to low and moderate; bare ground and litter increased and all other fractional components decreased in the year of burn relative to the prior year.
Our method was designed to be more aggressive in finding and labeling change as compared to previous efforts; we found 92.8% of pixels had change in one or more fractional component as compared to Shi et al. (2018) who found 36.1%. Increased sensitivity to change allows the identification of a more complete trend and facilitates analysis of subtler change. These benefits come at the risk of introducing noise. This risk is moderated by the CF approach, which considers the pixel response through time. This is important since land cover change tends to be persistent through time (Zhu and Woodcock 2014). With thematic mapping, more conservative change detection thresholds such as the requirement of three consecutive dates of change are acceptable, but with our fractional mapping representing the maximum vegetative cover within a year, doing so would negate any fractional component-climate relationship. Moreover, we would rather falsely map a non-changed pixel as changed (i.e., commission of change) than to falsely map a changed pixel as non-changed (i.e., omission of change), similar to Zhu and Woodcock (2014).
Our goal is to use this approach for multi-ecosystem trend analysis; hence, we plan on expanding production of historical fractional component mapping over a substantial portion of the West. Further methodological improvements are underway including improving the temporal domain accuracy of our products, developing better seamless composite image products, exploring non-linear trend analysis, and linking historical cover predictions with ecosystem site potential (Rigge et al. 2019) to make maps of condition trends. We also plan on developing more robust validation methods using repeat observations of high-resolution satellite imagery.

CONCLUSIONS
Temporal patterns of shrubland fractional component cover from 1986 to 2016 over the Northern Great Basin, USA, were assessed using a remote sensing-based automated modeling procedure. These data can be used to evaluate the influence of climate change, the role of weather on vegetation cover dynamics, and the suitability of management practices. Increasing the overall aggressiveness of change detection results in more pixels being labeled as change, while the addition of the change fraction method has improved the accuracy of change detection. The combined result was an improved time series of fractional component cover and tighter climate relationships.
Over the 30-yr study period, mean shrub cover declined and mean bare ground increased across the entire study region. While relatively few pixels had >10% cover change, a large majority had at least some change, suggesting that gradual change was more common than abrupt change. A minority of pixels had temporal patterns that were not well-fit by linear models. All fractional components had significant spatial relationships with WYPRCP, WYTMAX, and WYTMIN in all years. In particular, shrub and sagebrush cover respond positively to warming WYTMIN, resulting from the largest increases in WYTMIN being in the coolest and wettest areas, and negatively to warming WYTMAX because the largest increases in WYTMAX are in the warmest and driest areas.
Removal of years with phenological inconsistencies and pervasive cloud cover will result in a more accurate temporal trend than when included in analysis. The trade-off of a lower n is justified as temporal density appears to have only a modest impact on trends and climate relationships until n ≤ 6, but multi-year gaps are proportionally more influential. The minimum n needed will vary among regions, analysis objectives, what is being mapped, and fractional versus thematic mapping units. Gradual change analysis is likely to be less sensitive to n than abrupt change, and similarly fractional change less sensitive than thematic change. Future work will expand this mapping over a large region of the western United States.