Get 20M+ Full-Text Papers For Less Than $1.50/day. Subscribe now for You or Your Team.

Learn More →

Comparative Quality and Trend of Remotely Sensed Phenology and Productivity Metrics across the Western United States

Comparative Quality and Trend of Remotely Sensed Phenology and Productivity Metrics across the... remote sensing Article Comparative Quality and Trend of Remotely Sensed Phenology and Productivity Metrics across the Western United States 1 , 1 1 2 Ethan E. Berman * , Tabitha A. Graves , Nate L. Mikle , Jerod A. Merkle , 3 3 Aaron N. Johnston and Geneva W. Chong U.S. Geological Survey, Northern Rocky Mountain Science Center, Glacier Field Station, 38 Mather Drive, West Glacier, MT 59936, USA; tgraves@usgs.gov (T.A.G.); natemikle@gmail.com (N.L.M.) Department of Zoology and Physiology, University of Wyoming, Department 3166, 1000 E University Ave, Laramie, WY 82071, USA; jmerkle@uwyo.edu U.S. Geological Survey, Northern Rocky Mountain Science Center, 2327 University Way, Suite 2, Bozeman, MT 59715, USA; ajohnston@usgs.gov (A.N.J.); Geneva_chong@usgs.gov (G.W.C.) * Correspondence: bermane@gmail.com Received: 23 June 2020; Accepted: 30 July 2020; Published: 7 August 2020 Abstract: Vegetation phenology and productivity play a crucial role in surface energy balance, plant and animal distribution, and animal movement and habitat use and can be measured with remote sensing metrics including start of season (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), and integrated vegetation indices. However, for most metrics, we do not yet understand the agreement of remotely sensed data products with near-surface observations. We also need summaries of changes over time, spatial distribution, variability, and consistency in remote sensing dataset metrics for vegetation timing and quality. We compare metrics from 10 leading remote sensing datasets against a network of PhenoCam near-surface cameras throughout the western United States from 2002 to 2014. Most phenology metrics representing a date (SOS, PIRGd, POS, and EOS), rather than a duration (length of spring, length of growing season), better agreed with near-surface metrics but results varied by dataset, metric, and land cover, with absolute value of mean bias ranging from 0.38 (PIRGd) to 37.92 days (EOS). Datasets had higher agreement with PhenoCam metrics in shrublands, grasslands, and deciduous forests than in evergreen forests. Phenology metrics had higher agreement than productivity metrics, aside from a few datasets in deciduous forests. Using two datasets covering the period 1982–2016 that best agreed with PhenoCam metrics, we analyzed changes over time to growing seasons. Both datasets exhibited substantial spatial heterogeneity in the direction of phenology trends. Variability of metrics increased over time in some areas, particularly in the Southwest. Approximately 60% of pixels had consistent trend direction between datasets for SOS, POS, and EOS, with the direction varying by location. In all ecoregions except Mediterranean California, EOS has become later. This study comprehensively compares remote sensing datasets across multiple growing season metrics and discusses considerations for applied users to inform their data choices. Keywords: phenology; productivity; forage; green up; growing season; phenocam; timing; land surface phenology; elk; deer 1. Introduction Vegetation phenology and productivity are important drivers of ecosystem function by influencing processes as varied as surface energy balance [1] and plant species distribution [2]. Shifts in the timing of plant seasonality occur largely due to annual weather patterns and climate change [3], Remote Sens. 2020, 12, 2538; doi:10.3390/rs12162538 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 2 of 27 and these dynamics have consequences throughout terrestrial ecosystems. For example, phenology and productivity patterns strongly influence animal behavior, survival, and population dynamics [4]. Shifts in the timing and amount of vegetation can lead to trophic mismatch [5] and put stress on migratory species resilience and adaptive capacity [6]. Thus, to understand changes to the scale, rate, spatial configuration, and variability of ecological processes, we need to accurately measure vegetation phenology and productivity at broad spatial and temporal scales [7]. In recent decades, land surface phenology (LSP), the study of vegetation phenology and productivity from remote sensing, has revolutionized our understanding of ecological responses to phenological change [8]. LSP observations broaden the spatial scale of data to previously unattainable extents, enabling analyses of regional and continental vegetation patterns over time, with certain datasets extending back more than 30 years. Studies of LSP metrics in recent decades in North America have shown trends toward later senescence in the fall but inconclusive evidence for earlier spring green-up [9]. Green-up trends within ecological communities are complex, with plant species in the same community often showing opposite responses in both sign and magnitude to general warming patterns [10]. These trend analyses, along with annual LSP metrics, provide users (researchers, biologists, ecologists, natural resource managers, etc.) from a variety of fields with important insights into ecosystem processes and changing landscape dynamics driven by weather, climate, and human and natural disturbances. However, few studies have provided regional summaries of change that are accessible to users who are managing resources [7,11] or provided information on trends in variability of phenology over time. Insights into changing variability is important for wildlife researchers, as environmental predictability may shift habitat use, population density, and movement patterns [12]. Wildlife biologists began using LSP metrics in the early 2000s to better understand the influence of vegetation on the dynamics and distribution of animal populations including birds [13,14] and ungulates such as elk and deer [4,8,15–19]. For example, the timing of spring has implications for the fitness and body condition of ungulates [16,20–22] and peak instantaneous rate of green-up date (PIRGd), the date half way between start of season (SOS) and peak of season (POS) and representative of optimal forage quality, explains migratory patterns of ungulates surfing the green wave [23–25]. In addition, spatial heterogeneity of plant phenology, which may be declining due to warming temperatures, relates to the reproduction rates of caribou [26]. Autumn phenology, which has received considerably less attention than that of spring, can influence body mass and overwinter survival of mule deer [4] and migratory patterns of elk and red deer [18,27,28]. In addition to phenology, vegetation productivity is closely correlated to greenness indices such as Normalized Di erence Vegetation Index (NDVI) [8] and explains ungulate habitat use [24,29], health characteristics [30], and demographic parameters [16]. Changes in the seasonal timing of LSP metrics and the development of new datasets from various satellite platforms have received substantial focus from the remote sensing field [9,31–34]. However, despite the importance of LSP metrics to ecological applications, few studies have tested the quality of competing datasets in measuring LSP metrics against ground or near-surface observations or examined their relative agreement across various land cover types (but see [32,35]). Because LSP metrics synthesize information from millions of individual plants and multiple species, large biases across diverse vegetation communities can occur, and metrics may not directly represent the biological processes of the vegetation of interest [36]. Di erences in the processing algorithm of LSP data, such as the use of logistic curve fitting techniques versus splines, can also greatly impact performance [37]. Although users choose datasets based on their desired temporal and spatial scale, often several datasets exist with similar resolution yet di erent processing methodologies. A comparison of the quality of freely accessed and commonly derived LSP datasets against near-surface observations would assist users in selecting the best dataset for their application. Here, we provide a comparison of commonly used phenology and productivity metrics derived from 10 freely available LSP datasets. We examine correlation and bias in the western United States from 2002 to 2014 with near-surface observations from PhenoCam, a network of cameras spread Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 33 Remote Sens. 2020, 12, 2538 3 of 27 phenology and productivity metrics derived from remote sensing an indication of the strengths and weaknesses of different datasets, especially as related to land cover. We also compare long-term throughout North America providing data multiple times per day [38]. We provide users of phenology (1982–2016) trends in phenology from two leading LSP data products to identify spatial patterns, and productivity metrics derived from remote sensing an indication of the strengths and weaknesses assess the agreement about changing vegetation dynamics over time, and describe changing of di erent datasets, especially as related to land cover. We also compare long-term (1982–2016) trends variability in phenology. in phenology from two leading LSP data products to identify spatial patterns, assess the agreement about changing vegetation dynamics over time, and describe changing variability in phenology. Figure 1. The western United States (11 states) overlaid with Level 1 Ecoregions [39] and the locations and vegetation types of 29 PhenoCam sites. Note that some PhenoCam locations are slightly jittered to Figure 1. The western United States (11 states) overlaid with Level 1 Ecoregions [39] and the locations prevent overlap and therefore may not represent exact location (but are always in same ecoregion). and vegetation types of 29 PhenoCam sites. Note that some PhenoCam locations are slightly jittered to prevent overlap and therefore may not represent exact location (but are always in same ecoregion). 2. Materials and Methods We evaluated six phenology and three productivity metrics across the western United States, 2. Materials and Methods which represents a diverse range of ecosystems and land cover types (Figure 1). To evaluate agreement We evaluated six phenology and three productivity metrics across the western United States, with PhenoCam data, we focused on 10 datasets from 2002 to 2014 (Table 1), comparing overall which represents a diverse range of ecosystems and land cover types (Figure 1). To evaluate agreement and agreement within land cover types. We then assessed long-term trends using a subset agreement with PhenoCam data, we focused on 10 datasets from 2002 to 2014 (Table 1), comparing of more strongly associated datasets from 1982 to 2016. overall agreement and agreement within land cover types. We then assessed long-term trends using Phenology and productivity metrics can be derived from reflectance values recorded in optical a subset of more strongly associated datasets from 1982 to 2016. satellite imagery. Optical satellite imagery captures di erences in the reflectance of vegetation through Phenology and productivity metrics can be derived from reflectance values recorded in optical phenology cycles. Most commonly users measure vegetation indices as a ratio between reflectance in satellite imagery. Optical satellite imagery captures differences in the reflectance of vegetation the visual and near-infrared parts of the electromagnetic spectrum, with slightly di erent formulas through phenology cycles. Most commonly users measure vegetation indices as a ratio between for reflectance NDVI and in Enhanced the visual and near Vegetation Index -infrared parts (EVI; [40of ]). the el Fromectrom curvesag describing netic spectrum, wi the cycleth slightl of change, y users different for extract and mulas for ND evaluate key VI LSP and phenology Enhanced Ve and getatio productivity n Index (E metrics VI; [40] including ). From cuthose rves d evaluated escribing ther he e cycle of change, users extract and evaluate key LSP phenology and productivity metrics including (Figure 2): SOS, PIRGd, POS, end of season (EOS), length of spring (LOSp), length of growing season those evaluated here (Figure 2): SOS, PIRGd, POS, end of season (EOS), length of spring (LOSp), (LGS), integrated vegetation index (IVI), peak vegetation index (PVI), and amplitude of vegetation length of growing season (LGS), integrated vegetation index (IVI), peak vegetation index (PVI), and index (AVI). amplitude of vegetation index (AVI). Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 4 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 33 Figure 2. Six phenology and three productivity metrics shown through an example growing season Figure 2. Six phenology and three productivity metrics shown through an example growing season by by day of year (DOY). Note that these metrics are described through values of Normalized Difference day of year (DOY). Note that these metrics are described through values of Normalized Di erence Vegetation Index (NDVI)/Enhanced Vegetation Index (EVI) but certain datasets may produce metrics Vegetation Index (NDVI)/Enhanced Vegetation Index (EVI) but certain datasets may produce metrics using other methods/algorithms. using other methods/algorithms. Most datasets we evaluated use NDVI or EVI calculated from the surface reflectance of the Most datasets we evaluated use NDVI or EVI calculated from the surface reflectance of the Moderate Moderate Resolution Imaging Spectroradiometer (MODIS) or Advanced Very High Resolution Resolution Imaging Spectroradiometer (MODIS) or Advanced Very High Resolution Radiometer Radiometer (AVHRR) sensors (Table 1). We also evaluated one dataset, the National Phenology (AVHRR) sensors (Table 1). We also evaluated one dataset, the National Phenology Network first leaf Network first leaf spring index, hereinafter referred to as NPN, that models spatially explicit spring index, hereinafter referred to as NPN, that models spatially explicit temperature measurements temperature measurements parametrized via an extensive network of in situ phenological parametrized via an extensive network of in situ phenological observations [41]. This is the only observations [41]. This is the only dataset evaluated that does not incorporate optical satellite dataset evaluated that does not incorporate optical satellite imagery, but we included it because it is imagery, but we included it because it is readily available and provides annual SOS estimates across readily available and provides annual SOS estimates across the contiguous United States (CONUS). We the contiguous United States (CONUS). We also assessed two net primary productivity (NPP) also assessed two net primary productivity (NPP) datasets (CONUS Landsat NPP and CONUS MODIS datasets (CONUS Landsat NPP and CONUS MODIS NPP) that combine remote sensing with NPP) that combine remote sensing with meteorological data and plant physiological parameters [42]. meteorological data and plant physiological parameters [42]. Because they represent the total Because they represent the total productivity across the season, we include NPP in the IVI comparisons. productivity across the season, we include NPP in the IVI comparisons. AVI and PVI both measure AVI and PVI both measure the height of the curve at the peak, but AVI only includes the height from the height of the curve at the peak, but AVI only includes the height from the baseline to the peak. the baseline to the peak. Spatial resolutions spanned 30 m to approximately 5 km (0.05 degrees) and all Spatial resolutions spanned 30 m to approximately 5 km (0.05 degrees) and all data were freely available. Most datasets required minimal code and memory because the relevant phenology and data were freely available. Most datasets required minimal code and memory because the relevant productivity metrics were pre-processed and could be downloaded directly (e.g., SOS, POS, and EOS; phenology and productivity metrics were pre-processed and could be downloaded directly (e.g., SOS, Table 1). When necessary, we derived metrics (e.g., PIRGd and LOSp) from available metrics (see POS, and EOS; Table 1). When necessary, we derived metrics (e.g., PIRGd and LOSp) from available Appendix A Table A1 and Equations (A1)–(A7) for a complete description of metric availability, metrics (see Appendix A Table A1 and Equations (A1)–(A7) for a complete description of metric source, and derivation formula). We calculated all metrics for one dataset based on MODIS availability, source, and derivation formula). We calculated all metrics for one dataset based on MODIS MOD09Q1 NDVI (referred to as DLC for “double logistic curve”), fitting a double logistic curve and MOD09Q1 NDVI (referred to as DLC for “double logistic curve”), fitting a double logistic curve and accounting for snow cover following methods of [24], because of its use in many wildlife applications. accounting for snow cover following methods of [24], because of its use in many wildlife applications. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 5 of 27 Table 1. Abridged data table summarizing the 10 freely available datasets evaluated against PhenoCam. See Table A1 and Equations (A1)–(A7) for more information including methods used to derive individual metrics and data acquisition information. Significant Post-processing Dataset Base Input Spatial Resolution Temporal Coverage Accessible Online? Under Production? Spatial Coverage Reference Time Requirement? AVHRRP NDVI 1 km 1989–2014 yes no no CONUS [43] eMODIS NDVI 250 m 2001–2018 yes no yes CONUS [44] MCD12Q2 EVI 500 m 2001–2017 yes no yes Global [45] CMGLSP EVI 0.05 1982–2016 no no yes Global [46] VIPPHEN-EVI2 EVI 0.05 1981–2016 yes no no Global [47] VIPPHEN-NDVI NDVI 0.05 1981–2016 yes no no Global [47] NPN (First Leaf Spring Index) Temp, in situ observations 0.0417 1981–2019 yes no yes CONUS [41] Landsat NPP NDVI, meteorological inputs 30 m 1986–2018 yes no yes CONUS [42] MODIS NPP NDVI, meteorological inputs 250 m 2001–2018 yes no yes CONUS [42] DLC NDVI 250 m 2000–present yes yes yes Global [23] PhenoCam Green chromatic coordinate Point 2000–present yes no yes U.S. and Canada [48] 1 2 Collection 6 of eMODIS starts in 2003, but Collection 5 data is available from 2001. This analysis used Version 5 for 2002; Phenocam also has a limited number of sites in Panama, Hawaii, and Europe Remote Sens. 2020, 12, 2538 6 of 27 To compare datasets, we evaluated the agreement of LSP metrics from each dataset with metrics from PhenoCam data [48] at 29 sites (106 site-years) across the western United States (Figure 1) from 2002 to 2014. Using the R package phenocamr [49], we downloaded 3 day window PhenoCam data and calculated the 90th percentile green chromatic coordinate (GCC; Equation (1)). The 90th percentile value e ectively accounts for day-to-day variation due to weather and illumination patterns [50]. We defined rising and falling phenophases with a 10% amplitude threshold (to derive SOS and EOS) as described in PhenoCam documentation [48,50]. We chose the 10% threshold because it more closely resembles SOS than a 25% or 50% threshold and minimizes the bias between PhenoCam transition dates and MODIS transition dates [35]. GCC is a vegetation index derived from RGB camera imagery and is defined as: DNG GCC = (1) (DNR + DNG + DNB) where DNG, DNR, and DNB represent the digital number (i.e., strength) of the green, red, and blue channels, respectively. To assess the agreement of LSP metrics with PhenoCam, we compared correlation using the coecient of determination (R ) and the magnitude of di erence using mean bias. These common measures of statistical agreement have been previously used to compare LSP and near-surface phenology [36,51,52]. Mean bias is defined as: Mean bias = (est obs ) (2) i i i=1 where est and obs are the ith estimate (from LSP dataset) and observation (from PhenoCam) i i respectively. For productivity metrics, scales di er across datasets and thus we focused on correlation. We compared both overall agreement and agreement by landcover type, using the following vegetation types defined for PhenoCam sites (PhenoCam metadata: grasslands (GR, 45 site-years), shrublands (SH, 7 site-years), deciduous/broadleaf (DB, 27 site-years), evergreen (EN, 25 site-years), and wetlands (2 site-years). We excluded two sites in the Mediterranean California ecoregion that displayed a non-northern-temperate seasonal signal (SOS > 225) because green-up begins in late fall with rainfall and ends in early spring with drought. We excluded wetlands from the analysis by land cover classification due to the limited number of sites. We analyzed long-term trends (1982–2016) across the western U.S. using CMGLSP and VIPPHEN_EVI2, which agreed best with PhenoCam of the datasets that extended back to the 1980s. Of the six phenology metrics, we reported trends for SOS, POS, and EOS, because the other three are derived directly from these metrics and PIRGd spatial patterns are similar to SOS and POS. Using the Theil–Sen slope, the median of slopes between all pairs of observations within a pixel [53], we assessed each metric by pixel, including only pixels with data for at least 18 of 35 years. As a non-parametric test the Theil–Sen is more robust to outliers and provides higher statistical power when non-normality exists [54]. Negative slopes indicate change to earlier dates and positive slopes indicate change to later dates. We did not screen pixels for disturbances as the goal of this work is to identify patterns and magnitudes of change, rather than make inferences about the causes of phenological change. We evaluated overall variation within each pixel based on the standard deviation for each metric and whether variability has increased across years by applying the Theil–Sen slope to the absolute values of the residuals of the trend estimate against time. In addition, we measured consistency between the two datasets, defined as the by-pixel agreement that phenology dates are trending earlier or later. Lastly, to assess agreement in regional patterns of change, we report the proportion of area where both datasets agreed on direction of change by ecoregions using the United States Environmental Protection Agency Level 1 Ecoregions of North America dataset [39]. We used the statistical computing environment R for all analyses [55] and R package rkt [56] to calculate the Theil–Sen slope estimator. Remote Sens. 2020, 12, 2538 7 of 27 3. Results 3.1. Agreement of LSP metrics with PhenoCam We compared LSP metrics to 106 site-years of PhenoCam data from 29 sites. We found substantial variation in agreement between remotely sensed datasets and near-surface observations (Figure 3; Figure 4; Appendix B Tables A2–A6 for full results). Across datasets and six phenology metrics, R ranged from 0.03 to 0.37 (SOS), 0.20 to 0.55 (PIRGd), 0.25 to 0.54 (POS), 0.16 to 0.45 (EOS), 0.01 to 0.11 (LOSp), and 0 to 0.10 (LGS). Absolute value of mean bias (in days) ranged from 4.39 to 25.49 (SOS), 0.38 to 17.66 (PIRGd), 7.05 to 23.82 (POS), 0.52 to 37.92 (EOS), 12.45 to 44.66 (LOSp), 3.78 to 58.47 (LGS). Mean bias was lowest for SOS and PIRGd. Four datasets matched SOS and PIRGd from PhenoCam within 7 days (MCD12Q2, VIPPHEN_EVI2, DLC, and eMODIS for SOS; CMGLSP, MCD12Q2, VIPPHEN_EVI2, and VIPPHEN_NDVI for PIRGd). Generally, R values were higher for PIRGd and POS than SOS and EOS. Datasets were least correlated with PhenoCam transition dates for LOSp and LGS. In pairing a high R value with a low mean bias, MCD12Q2 and CMGLSP had best overall agreement with PhenoCam. Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 33 Figure 3. R and 2 mean bias (in days) for six phenology metrics classified by land cover type identified Figure 3. R and mean bias (in days) for six phenology metrics classified by land cover type identified in PhenoCam metadata and overall for all data points. The phenology metrics are start of season (SOS), in PhenoCam metadata and overall for all data points. The phenology metrics are start of season peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), length of (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), spring length of (LOSp) spring (LOSp) and length of growing seas and length of growing season (LGS). The on (LG classifications S). The classif areic grasslands ations are g (GR), rasslands shrublands (GR), (SH), shru deciduous blands (S /br H), oadleaf decidu (DB), ous/broadleaf (D evergreen (EN) B), e and verg overall reen (EN) and (OV). In over general, all (OV) this illustrates . In general, thi how the s overall illustrate evaluat s how ion the overal falls between l eval the uation fall landcover s be -specific tween thvalues e landcfor over-spec each dataset ific valu and es for metric eachcombination. dataset and metric combination. Phenology metrics had large di erences in correlation and bias by land cover type (Figure 3). Most datasets agreed moderately well with PhenoCam in grasslands, shrublands and deciduous/broadleaf forests for the four key phenology dates (SOS, PIRGd, POS, and EOS). Duration Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam data and overall for all data points. The productivity metrics are integrated vegetation index (IVI), peak vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 33 Remote Sens. 2020, 12, 2538 8 of 27 metrics (LOSp and LGS) correlated poorly with PhenoCam across all datasets. At grassland sites (n = 45), dates for most datasets correlated well with PhenoCam dates. At shrubland sites (n = 7), eMODIS was best for SOS, PIRGd, and POS, whereas CMGLSP was best for EOS. At deciduous/broadleaf sites (n = 27), CMGLSP and NPN agreed best for SOS, DLC for POS, and MCD12Q2 and eMODIS for EOS. Most datasets agreed well for PIRGd in deciduous/broadleaf forests. Datasets generally showed poor agreement with PhenoCam at evergreen forest sites (n = 25), with the exception of MCD12Q2. Figure 3. R and mean bias (in days) for six phenology metrics classified by land cover type identified Overall, datasets better represented phenology than productivity metrics (Figure 4). in PhenoCam metadata and overall for all data points. The phenology metrics are start of season Across datasets, R for the three productivity metrics ranged from 0.00 to 0.16 (IVI), 0.00 to 0.15 (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), (PVI) and 0.00 to 0.34 (AVI), with low values in overall R masking high correlations in a few land length of spring (LOSp) and length of growing season (LGS). The classifications are grasslands (GR), cover types. For the IVI metric, MODIS_NPP had an R value of 0.9 at deciduous/broadleaf locations shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). In general, this and VIPPHEN_EVI2 and VIPPHEN_NDVI both agreed well. MCD12Q2 was highly correlated with illustrates how the overall evaluation falls between the landcover-specific values for each dataset and PhenoCam at evergreen locations for PVI and AVI. metric combination. Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam data data and overall for all data points. The productivity metrics are integrated vegetation index (IVI), and overall for all data points. The productivity metrics are integrated vegetation index (IVI), peak peak vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). 3.2. Historical Trend Analysis for Western United States Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Of the 10 datasets evaluated, five extend back to 1982 (CMGLSP, VIPPHEN_EVI2, VIPPHEN_NDVI, AVHRRP and NPP) and therefore provide an appropriate historical record for relatively long-term trend analysis. Based on the comparison with PhenoCam, we chose to analyze CMGLSP and VIPPHEN_EVI2, plus NPN for SOS (NPN results in Appendix C Figures A2–A4). CMGLSP most often had the ideal combination of high R and low mean bias across land cover types and VIPPHEN_EVI2 also agreed well with PhenoCam, especially in grasslands and shrublands. Long-term trends indicate large changes in key phenology dates in some places by more than 70 days (2 days per year for 35 years) between 1982 and 2016 (Figure 5; see Figures A1 and A5–A7 for mean and standard deviation maps and histograms of changes). Trend patterns had high spatial heterogeneity, with substantial interspersion of pixels with delayed and advanced phenology in both datasets and all metrics. The strength of the trends varied spatially, but in both datasets the areas of greatest change in date occurred in New Mexico and Arizona (later SOS), low-lying areas of California (earlier EOS), and the Great Basin (earlier SOS, mixed strong trends for EOS). EOS in the Southwest was also generally later but was more variable in California. Higher elevation areas across Montana, Washington, Idaho and Colorado showed trends toward later SOS but mixed results for EOS. Over time, variation in phenology dates slightly increased across most regions, with very large increases in variability corresponding to areas with larger changes in date, namely the Southwest for Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 33 3.2. Historical trend analysis for western United States Of the 10 datasets evaluated, five extend back to 1982 (CMGLSP, VIPPHEN_EVI2, VIPPHEN_NDVI, AVHRRP and NPP) and therefore provide an appropriate historical record for Remote Sens. 2020, 12, 2538 9 of 27 relatively long-term trend analysis. Based on the comparison with PhenoCam, we chose to analyze CMGLSP and VIPPHEN_EVI2, plus NPN for SOS (NPN results in Appendix C Figures A2–A4). SOS and POS, and scattered areas around the edges of the Great Basin as well as California for EOS. CMGLSP most often had the ideal combination of high R and low mean bias across land cover types The spatial patterns of agreement on the direction of trend across CMGLSP and VIPPHEN-EVI2 were and VIPPHEN_EVI2 also agreed well with PhenoCam, especially in grasslands and shrublands. also relatively heterogenous, with the least consistency occurring in Wyoming for SOS and in the Great Long-term trends indicate large changes in key phenology dates in some places by more than 70 Basin for EOS. days (2 days per year for 35 years) between 1982 and 2016 (Figure 5; see Figure A1 and Figure A5– The consistency results over the entire study region showed ~60% agreement in trend direction for A7 for mean and standard deviation maps and histograms of changes). Trend patterns had high all three metrics (Figure 6). For SOS and POS ~30% of pixels agreed on earlier and later dates in trends, spatial heterogeneity, with substantial interspersion of pixels with delayed and advanced phenology whereas for EOS, 44.4% of pixels in both datasets agreed on a later EOS (16.2% agreed on an earlier in both datasets and all metrics. The strength of the trends varied spatially, but in both datasets the EOS). In assessing consistency by ecoregion, we found agreement between the datasets in the range of areas of greatest change in date occurred in New Mexico and Arizona (later SOS), low-lying areas of 50–70% of pixels for most of the regions and metrics as well as agreement on the predominant direction California (earlier EOS), and the Great Basin (earlier SOS, mixed strong trends for EOS). EOS in the of change within ecoregions (Figure 6). For SOS, the two datasets agreed that more area trended Southwest was also generally later but was more variable in California. Higher elevation areas across towards earlier dates in Marine West Coast Forests, Mediterranean California, North American Deserts, Montana, Washington, Idaho and Colorado showed trends toward later SOS but mixed results for and Northwestern Forested Mountains, and later dates in Great Plains, Southern Semiarid Highlands, EOS. Over time, variation in phenology dates slightly increased across most regions, with very large and Temperate Sierras. For POS, the datasets agreed on more area trending toward earlier dates in increases in variability corresponding to areas with larger changes in date, namely the Southwest for Marine West Coast Forests, Mediterranean California, and North American Deserts, and later dates in SOS and POS, and scattered areas around the edges of the Great Basin as well as California for EOS. Great Plains, Northwestern Forested Mountains, Southern Semiarid Highlands, and Temperate Sierras. The spatial patterns of agreement on the direction of trend across CMGLSP and VIPPHEN-EVI2 were We found the most dataset agreement across ecoregions for EOS, in which all regions agreed on more also relatively heterogenous, with the least consistency occurring in Wyoming for SOS and in the area trending towards later dates except Mediterranean California. Great Basin for EOS. Figure 5. Change in date, consistency, and change in variation of CMGLSP and VIPPHEN_EVI2 Figure 5. Change in date, consistency, and change in variation of CMGLSP and VIPPHEN_EVI2 across three phenology metrics from 1982 to 2016. Negative values correspond to earlier dates (or fewer days) across three phenology metrics from 1982 to 2016. Negative values correspond to earlier dates (or and fewer day positive s) and pos values toitiv later e va dates lues to (or lat mor er date e days). s (orConsistency more days). Consisten shows agreement cy shows agreement in th in the direction ofe change between datasets. Change in date is the regressed trend slope multiplied by the number of direction of change between datasets. Change in date is the regressed trend slope multiplied by the years. Change in variation is the regressed trend slope, multiplied by the number of years, of the number of years. Change in variation is the regressed trend slope, multiplied by the number of years, absolute of the ab value solute of value of the residuals the residuals of the of the trend estimate trend est against imate a time. gainst time. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 33 The consistency results over the entire study region showed ~60% agreement in trend direction for all three metrics (Figure 6). For SOS and POS ~30% of pixels agreed on earlier and later dates in trends, whereas for EOS, 44.4% of pixels in both datasets agreed on a later EOS (16.2% agreed on an earlier EOS). In assessing consistency by ecoregion, we found agreement between the datasets in the range of 50–70% of pixels for most of the regions and metrics as well as agreement on the predominant direction of change within ecoregions (Figure 6). For SOS, the two datasets agreed that more area trended towards earlier dates in Marine West Coast Forests, Mediterranean California, North American Deserts, and Northwestern Forested Mountains, and later dates in Great Plains, Southern Semiarid Highlands, and Temperate Sierras. For POS, the datasets agreed on more area trending toward earlier dates in Marine West Coast Forests, Mediterranean California, and North American Deserts, and later dates in Great Plains, Northwestern Forested Mountains, Southern Semiarid Highlands, and Temperate Sierras. We found the most dataset agreement across ecoregions Remote Sens. for EO 2020 S, ,in w 12, 2538 hich all regions agreed on more area trending towards later dates except Mediterranea 10 n of 27 California. Figure 6. Consistency of CMGLSP and VIPPHEN_EVI2 long-term (1982–2016) trend direction over the Figure 6. Consistency of CMGLSP and VIPPHEN_EVI2 long-term (1982–2016) trend direction over entire study area and classified by ecoregion. the entire study area and classified by ecoregion. 4. Discussion 4. Discussion In our comparison of 10 leading LSP datasets against a network of PhenoCam near-surface cameras, In our comparison of 10 leading LSP datasets against a network of PhenoCam near-surface we found promising results in terms of R and mean bias in some landcover types, but substantial cameras, we found promising results in terms of R and mean bias in some landcover types, but substantial variation between datasets and metrics (similar to [32] for SOS estimates). This highlights variation between datasets and metrics (similar to [32] for SOS estimates). This highlights a need for a need for improved communication and more extensive dialogue on the strengths, weaknesses and improved communication and more extensive dialogue on the strengths, weaknesses and development development goals of LSP datasets, especially with the proliferation of applied data users from a goals of LSP datasets, especially with the proliferation of applied data users from a variety of research variety of research and management practices. and management practices. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing When choosing an appropriate dataset for analysis, applied users should consider study location, years needed, spatial resolution, processing capacity, and quality in di erent land cover classes. For phenology metrics, our results indicate MCD12Q2 and CMGLSP best matched PhenoCam observations derived using a 10% amplitude threshold and the 90th GCC percentile (Figure 3). For those interested in phenology prior to the 2000s, CMGLSP has global coverage and extends back to 1982, but has a coarse 0.05 degree (~5 km) spatial resolution and can only be acquired by request. VIPPHEN_EVI2 has the same temporal coverage and spatial resolution, can be downloaded directly, and agreed almost as well. For users working with high resolution data (such as from GPS collars in the wildlife field) and conducting analyses after 2001, MODIS-based datasets at 250 to 500 m spatial resolution provide finer-scale observations. MCD12Q2, with 500 m spatial resolution and global coverage, agreed best with PhenoCam and was the only dataset that had high correlations in evergreen forests. eMODIS also agreed well and has 250 m spatial resolution but is only available in the CONUS. The DLC method can be applied globally (we used MOD09Q1 as the base input) but requires substantial processing which may be prohibitive depending on study area size. This method also provides daily values of NDVI and instantaneous rate of green-up (IRG), which are useful for certain movement ecology and habitat questions such as those relating to the green-wave hypothesis [21,23–25,57]. Remote Sens. 2020, 12, 2538 11 of 27 Readily available datasets from MODIS (e.g., MOD13A1, MOD13Q1, and eMODIS NDVI) only provide cleaned and pre-processed weekly to bi-weekly NDVI/EVI values and therefore could provide a coarser resolution IRG time series. Based on the di erences in dataset correlation and bias we found between land cover types, users should consider the predominant land cover of their study area in choosing a phenology or productivity dataset (Figure 3; Figure 4). For instance, eMODIS agreed best for PIRGd in shrublands, yet MCD12Q2 was better in grasslands and evergreen forests. Large di erences in correlation were observed between productivity datasets, which correlated poorly overall but well in certain landcover types. For IVI, DLC agreed well in grasslands, but MODIS_NPP agreed best in deciduous/broadleaf forests. The lowest correlations between LSP estimates and PhenoCam observations were in evergreen forests, where identifying phenological metrics is challenging because of the small annual change in greenness amplitude and over-saturation [35,58]. Using the red chromatic coordinate to process PhenoCam transition dates, as opposed to GCC, has shown potential to improve predictions of SOS and EOS in these environments [59]. Of all the LSP datasets, MCD12Q2 agreed best with PhenoCam at evergreen sites, possibly because it is EVI based [40]. Ecologists studying wildlife in areas with multiple land cover types should be aware that phenology metrics in some land cover types are more reliable than others, which could influence results for analyses comparing wildlife GPS locations across di erent land cover classes [17,60,61]. Quality of LSP metrics may di er geographically, and therefore our results are not necessarily indicative of dataset quality in other parts of North America or around the world [62]. Quantifying phenology and productivity trends through time is crucial to improve our understanding of ecosystem processes, such as how changing forage in critical habitat or management units a ects wildlife and how to manage for future climate and phenology scenarios [8,26,63,64]. We found a 44.4% agreement between CMGLSP and VIPPHEN_EVI2 trends toward a later EOS across the Western U.S and only a 16.2% agreement toward an earlier EOS, which is consistent with other studies signaling a general pattern toward a later EOS [9,65,66]. Trends are influenced by a variety of factors, including temperature and precipitation, species succession, human and natural disturbance, and land management practices (as observed by [37] from 1982 to 2006). The complex interactions between these factors make it dicult to interpret the ecological significance of large and spatially heterogeneous trends, especially across diverse ecosystems. The large phenology trends and high variability we observed over a 35 year period likely reflect real changes in remote sensing metrics, but highlight the need to understand what exactly is changing throughout time, the role of data processing choices, and whether the changes are biologically significant. We know that patterns of temperature and precipitation in some areas are complex, leading to variability in the timing, seasonality, and spatial heterogeneity of snow coverage [67], soil moisture, drought, and storms that may all influence the timing, variability, and heterogeneity of vegetation phenology and productivity. Variability can be an important component to ecosystem processes and to understanding the degree and mechanism of trends. For instance, vegetation strongly a ected by climate conditions, such as invasive grasses responding to rainfall, can cause high year-to-year variability and therefore introduce uncertainty into the direction of overall trend, as well as have large e ects on the ecosystem [68,69]. In our results, some regions in which LSP metrics changed by more than 70 days often showed high and increasing variability. Though these trends may be showing real changes in phenology throughout time within a pixel, it is challenging to interpret whether these changes are driven by variability caused by weather, invasive species, or other factors, or actually represent the changing dynamics of vegetation of interest, such as important forage species. This overall assessment is a first step in understanding the overall patterns and degree of change across the CONUS. Several challenges exist in producing accurate LSP estimates and assessing their quality. Large di erences in LSP metrics and trend are related to sensor-specific characteristics, such as revisit time and spectral resolution, as well as processing algorithms used to extract metrics, which may include Remote Sens. 2020, 12, 2538 12 of 27 cloud and atmospheric correction, gap filling and curve fitting. Greenness (NDVI/EVI/GCC) values di er between satellite platforms, as well as with near-surface cameras, and are therefore not always comparable without applying translation equations [70]. Future work may scale PhenoCam and satellite derived productivity metrics to enable a comparison of mean bias and other statistics. Di erences in consistency (by-pixel agreement on trend direction) between CMGLSP and VIPPHEN_EVI2 result solely from di erent processing algorithms, as both CMGLSP and VIPPHEN_EVI2 are based on EVI2 calculated from AVHRR and MODIS at a 0.05 degree spatial resolution (for CMGLSP, see [46]; for VIPPHEN_EVI2, see [47]). These two datasets yielded conflicting yet large trends in phenology around the Great Basin and Eastern Washington and Oregon, specifically for EOS. These areas experienced large fires, landcover changes, and other disturbances. It is possible that the di erent approaches to fitting NDVI/EVI curves and extracting metrics for these datasets are a ected di erently by these landcover type changes and thus lead to these conflicting signals. In general, LSP datasets use di erent processing methodologies and di erent threshold values [36,37,71], and dataset development goals rarely correspond to specific biological events important to data users. For instance, biological events recorded by a botanist to mark the start of growing season may include first leaf or flower dates, whereas an LSP dataset such as MCD12Q2 captures SOS as the date when EVI crosses the 15% threshold of AVI [45]. Ref. [32] used in situ first leaf observations to validate SOS from various LSP datasets and found high root mean squared error (~20 days), indicating a need to better understand the link between biological events and image reflectance values. Our use of near-surface cameras removes variability that may occur when comparing LSP metrics with the timing of biological processes from in situ observations. Although some research assesses LSP quality through comparison with other kinds of ground-based datasets [72,73], this may confound di erences resulting from discrepancies between biological events and their reflectance values with di erences in image quality, which can vary inter-annually and regionally [62]. Datasets also vary in how they deal with growing seasons that span calendar years or have multiple annual cycles. Most datasets calculate phenology metrics based on a continuous time series yet report metrics in discrete annual data layers. The one exception in our analysis is the DLC dataset, which fits a curve to a single year of data. When phenology cycles span calendar years, such as when SOS occurs in the fall and/or when EOS occurs in late winter/early spring, users should ensure they understand the approach used to store the metrics within an annual data layer. Datasets handle patterns of multiple or missing annual wetness cycles in di erent ways, with some datasets reporting up to three annual cycles while others only report phenology dates from a single cycle that is not indexed to a general time of year. We most commonly see multiple annual cycles of green-up in arid and monsoonal environments, such as in the Southwest [74] and Mediterranean California, but this varies by year.. When the number of cycles varies, such as when drought years lead to missed wet periods [75], or seasonal weather patterns dictate two to three annual cycles, analyses of trends by year, as we have conducted, could suggest large changes in dates that might not reflect the complexity of the ecological impacts well. We briefly evaluated these potential issues, and did find strong positive trends in SOS and POS in areas prone to multiple annual cycles, but overall found that both varying numbers of cycles and small amplitude cycles had more limited extent than could easily explain the large changes in trend we found. Inherently, the scale of near-surface PhenoCam observations does not match those of remote sensing pixels, which include heterogenous vegetation and even land cover class [36,38,76]. The mismatch can be substantial even for our comparisons of PhenoCam with MODIS-based 250 or 500 m spatial resolution products versus CMGLSP at ~5 km resolution, in which the synthesis of large pixels could include billions of plants. LSP observations e ectively provide only a broad picture of landscape phenology, and the scale of imagery has important implications for users. For example, CMGLSP at ~5 km spatial resolution may not provide the spatial variability needed to assess the behavior of a white-tailed deer or other animals with home ranges within one or two pixels, whereas it may be useful for understanding movement patterns of long distance migrating animals such as eagles or Remote Sens. 2020, 12, 2538 13 of 27 mule deer. Given the heterogeneity of vegetation within large pixels, we were surprised at the high agreement of CMGLSP with PhenoCam, at a spatial resolution of 0.05 degrees. The dynamics of LSP metrics in large pixels are complex and do not necessarily represent the average of smaller pixels [77]. For SOS, [78] found that the heterogeneity of landcover and SOS, as well as vegetation growth speed during spring, all influenced the scale e ect. Variation in spatial footprints of PhenoCam data and resolution of LSP datasets could also add noise to the comparison of LSP and near-surface datasets. New fine-resolution datasets, such as daily 30 m EVI products developed using fusion between MODIS and Landsat [79,80] or combinations of Landsat and Sentinal-2 imagery to increase temporal resolution and enable 30 m phenology retrievals [81], could better match the scale of reference data and user analysis objectives. Applied data users will benefit from future developments in processing capacity and dataset construction, a better understanding of phenology drivers, and greater understanding of how algorithms intersect with phenology predictions. The growing popularity of platforms able to process large data quickly (such as Google Earth Engine) may render daily, global datasets derived from the DLC method readily available in the future. In terms of improving LSP estimates, the use of mechanistic models to predict key metrics [82–84] may help to address data quality issues and discrepancies across land cover types and ecoregions. These models couple remote sensing data with local observations or other models of elevation, temperature, precipitation, and plant phenology to improve phenology and productivity metrics. While still limited, in the future these models may be developed into a single framework to derive more relevant phenology and productivity estimates across diverse regions using localized data. Similarly, an increased understanding of the drivers of phenological changes and variability on a local or regional scale may help users to better interpret and plan for long-term patterns of change. For example, [85] showed that temperature, growing-season precipitation, and snowpack had larger e ects than most management and habitat treatments on annual phenology metrics in sagebrush communities in Southwestern Wyoming, but identified that treatments had some impact on IVI and late season phenology metrics. Finally, deeper examination of how processing and curve-fitting techniques influence the resulting NDVI/EVI time series and phenology metrics will improve understanding of how to interpret large di erences in phenology that are predicted from LSP data. 5. Conclusions Land surface phenology (LSP) data provide the means to assess large spatial and temporal patterns in environmental change and understand a variety of ecological questions, such as those related to surface energy balance and animal fitness, movement, and habitat use. The quality of these data depend on factors including cloud cover, atmospheric e ects, processing algorithm, land cover type, spatial resolution and regionality. We found highly variable agreement between 10 leading LSP datasets and PhenoCam and outlined dataset choices based on spatial and temporal coverage, processing capability, and agreement across ecosystems, providing a basis for informed and justifiable decisions about data choices. Most datasets we examined and many of the highest-quality datasets are readily available and do not require users to extract metrics, as is often performed by individual users using raw NDVI or EVI values. Given the popularity and importance of phenology and productivity data in a wide range of research and management fields [61,86], a better understanding and justification for the use of certain datasets can only help to bolster the e ectiveness and validity of results. Author Contributions: Conceptualization, T.A.G., N.L.M., J.A.M., A.N.J. and G.W.C.; methodology, E.E.B., T.A.G. and N.L.M.; software, E.E.B. and N.L.M.; validation, E.E.B.; formal analysis, E.E.B. and N.L.M.; resources, E.E.B. and T.A.G.; writing—original draft preparation, E.E.B. and N.L.M.; writing—review and editing, T.A.G., N.L.M., J.A.M., A.N.J. and G.W.C.; visualization, E.E.B.; supervision, T.A.G.; project administration, T.A.G.; funding acquisition, T.A.G.; E.E.B., T.A.G., and N.L.M. contributed equally. All authors have read and agreed to the published version of the manuscript. Remote Sens. 2020, 12, 2538 14 of 27 Funding: This work was supported by the U.S. Geological Survey Northern Rocky Mountain Science Center and the North Central Climate Adaptation Science Center. The Wyoming Landscape Conservation Initiative provides support to Chong and Johnston. The work of Ethan Berman was done under contract to the U.S. Geological Survey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Acknowledgments: We would like to thank the Wyoming Game and Fish Department, Bureau of Land Management Kemmerer, Montana Department of Fish, Wildlife and Parks, and the National Park Service for their support in the development of this research. In addition, we appreciate the help of Jill Randall, Troy Fieseler, Brent Jamison, Arvid Aase, Kelly Prott, Justin Gude, as well as David Wood and anonymous reviewers for their feedback and comments. Conflicts of Interest: The authors declare no conflict of interest. The external funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Appendix A Detailed table of phenology and productivity metrics along with formulas for metric calculations where not readily available. Remote Sens. 2020, 12, 2538 15 of 27 Table A1. Detailed table of phenology and productivity metrics. NPN First Leaf Spring Dataset AVHRRP eMODIS MCD12Q2 CMGLSP VIPPHEN-EVI2 VIPPHEN-NDVI Landsat NPP MODIS NPP DLC PhenoCam Index Temperature, weather MODIS NDVI MODIS NDVI Calibrated events, in situ composites, composites, MOD09Q1 Calibrated radiance data Time series of AVHRR and AVHRR and MODIS AVHRR and MODIS RGB photo time Input data radiance data observations of one lilac meteorological meteorological surface reflection (level 1-B). Aqua only NBAR-EVI2 MODIS EVI2 EVI2 NDVI series (level 1-B). and two honeysuckle inputs, biome inputs, biome 8 day composites species specific contraints specific contraints Data are smoothed Data are smoothed 2 band EVI Cleaned, NDVI with 2 step filter and 3 with 2 step filter and 3 Cleaned, NDVI (EVI2) calculated, and Eliminating year moving window year moving window Calculated at PRISM Daily GPP and Daily GPP and calculated, and determination combined into 7 outliers and approach, AVHRR approach, AVHRR minimum and maximum maintenance maintenance Cleaned, NDVI Raw data processing details combined into 7 day from raw data. GCC day composites filling dormant and MODIS continuity and MODIS continuity temperature values reach respiration values respiration values calculated composites using max Smoothed with using max NDVI period values is based on linear is based on linear a “stable” state. calculated calculated NDVI and quality Savitsky-Golay and quality regressions between regressions between filter overlapping data. overlapping data. SOS NDVI value NPP extracted from NPP extracted from PELT method, Cubic spline fit, phenological Average of the first day 2 piece logistic compared with average Phenological changes Phenological changes daily GPP, daily GPP, AIC-selected SOS is 15% of changes are each year that the model. SOS/EOS Did not receive of previous 36 days and are identified using are identified using maintenance maintenance smoothing spline, Metric extraction details AVI, EOS is 15% identified using properties of the three are local information SOS is when a trend half-max method (at half-max method (at respiration, and respiration, and SOS at 10% of AVI, of AVI hybrid piecewise individual species models min/max from shift is identified. only 35%) only 35%) annual growth annual growth EOS at 10% of green-down logistic model are met (SOS only) 2nd derivatives Opposite for EOS respiration respiration green-down curve Spatial resolution 1 km 250 m 500 m 0.05 0.05 0.05 0.0417 30 m 250 m 250 m Point Temporal coverage 1989–2014 2001–2018 2001–2017 1982–2016 1981–2016 1981–2016 1981–2019 1986–2018 2001–2018 2000–present 2000–present Accessible online? yes yes yes no yes yes yes yes yes yes yes Significant post-processing no no no no no no no no no yes no time requirement? Under production? no yes yes yes no no yes yes yes yes yes Spatial coverage CONUS CONUS Global Global Global Global CONUS CONUS CONUS Global U.S. and Canada Years acquired 2002–2014 2002–2014 2002–2014 1982–2016 1982–2016 1982–2016 1982–2016 2002–2014 2002–2014 2002–2014 2002–2014 Reference [43] [44] [45] [46] [47] [47] [41] [42] [42] [23] [48] http: https: https://www.ntsg. https://www.usgs.gov/ xiaoyang. //files.ntsg.umt.edu/ //www.usgs.gov/ https://earthdata. https: https: https://www.usanpn.org/ umt.edu/project/ https://earthdata. Retrieval location land-resources/eros/ zhang@sdstate. data/NTSG_ phenocamr land-resources/ nasa.gov/ //earthdata.nasa.gov/ //earthdata.nasa.gov/ data/spring_indices landsat/landsat- nasa.gov/ phenology edu Products/MOD17/ eros/phenology productivity.php MODIS_250/ Metrics available or calculated by authors? (A is available, C is calculated, N is not available/used) SOS A A A A A A A N N C C PIRGd C C C C C C N N N C C POS A A A A A A N N N C C EOS A A A A A A N N N C C LOSp C C C C C C N N N C C LGS A A C C A A N N N C C IVI A A A A A A N A A C C PVI A A C A A A N N N C C AVI A A A C A A N N N C C 1 2 Collection 6 of eMODIS starts in 2003, but Collection 5 is available from 2001. This analysis used Collection 5 for 2002; Some sites in Panama, Hawaii, and Europe. Remote Sens. 2020, 12, 2538 16 of 27 Metrics reflecting di erent components of forage phenology and productivity were derived (or readily provided) from the 10 datasets used in this analysis. In total we assessed ten di erent metrics, six related to phenology and three related to productivity. Many of the datasets provide the desired metrics, and when available these values were used. In other cases, the metrics were derived using the following equations (with the exception of the Bischof dataset, for which methods can be found in [23]: POS SOS PIRGd = SOS + (assuming logistic growth) (A1) EOS POS = MAX Greenness (A2) i=SOS LOSp = POS SOS (A3) LGS = EOS SOS (A4) EOS IVI = Greenness (A5) i=SOS PVI = AVI + SMV (A6) AVI = PVI SMV (A7) where PIRGd is peak instantaneous rate of spring greenup date, SOS is start of spring date, POS is peak of season, MAX is the maximum value for the given period, EOS is end of season, greenness is the daily value of greenness (whether EVI, NDVI, or GCC), LOSp is length of spring, LGS is length of growing season, IVI is integrated vegetation index, PVI is peak vegetation index, AVI is amplitude of vegetation index and SMV is spring minimum value. All dates were represented as Julian date for calculation Purp. Appendix B R and mean bias results of 10 datasets tested against PhenoCam observation. The number of observations is displayed in parenthesis. Table A2. Overall results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.07 (98) 0.20 (98) 0.33 (106) 0.22 (79) 0.02 (98) 0.00 (72) 0.00 (104) 0.13 (106) 0.29 (106) CMGLSP 0.37 (102) 0.52 (101) 0.39 (105) 0.20 (104) 0.03 (101) 0.06 (100) 0.03 (100) 0.00 (105) 0.00 (104) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.10 (88) NA (NA) NA (NA) MCD12Q2 0.35 (76) 0.55 (76) 0.54 (76) 0.35 (76) 0.02 (76) 0.02 (76) 0.01 (76) 0.01 (76) 0.02 (76) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (78) NA (NA) NA (NA) NPN 0.03 (104) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.26 (106) 0.42 (106) 0.32 (106) 0.16 (106) 0.01 (106) 0.04 (106) 0.15 (106) 0.13 (106) 0.34 (103) VIPPHEN_NDVI 0.36 (106) 0.44 (106) 0.37 (106) 0.23 (106) 0.09 (106) 0.10 (106) 0.16 (106) 0.15 (106) 0.26 (97) DLC 0.35 (106) 0.50 (106) 0.53 (106) 0.45 (106) 0.11 (106) 0.09 (106) 0.00 (106) 0.05 (106) 0.05 (106) eMODIS 0.26 (99) 0.33 (99) 0.25 (103) 0.22 (101) 0.05 (99) 0.05 (97) 0.02 (103) 0.11 (106) 0.19 (103) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 9.88 (98) 15.92 (98) 20.54 (106) 37.92 (79) 14.09 (98) 43.44 (72) CMGLSP 8.41 (102) 0.38 (101) 12.10 (105) 0.52 (104) 17.75 (101) 3.78 (100) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 4.39 (76) 0.83 (76) 7.05 (76) 19.92 (76) 12.45 (76) 24.32 (76) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 19.10 (104) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 4.39 (106) 2.69 (106) 10.76 (106) 23.58 (106) 16.15 (106) 27.97 (106) VIPPHEN_NDVI 25.49 (106) 4.16 (106) 18.17 (106) 32.98 (106) 44.66 (106) 58.47 (106) DLC 5.63 (106) 17.66 (106) 21.32 (106) 4.13 (106) 27.95 (106) 9.76 (106) eMODIS 4.61 (99) 15.74 (99) 23.82 (103) 23.26 (101) 24.26 (99) 25.15 (97) Remote Sens. 2020, 12, 2538 17 of 27 Table A3. Grasslands (GR) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.10 (41) 0.13 (41) 0.32 (45) 0.52 (38) 0.03 (41) 0.15 (34) 0.30 (45) 0.05 (45) 0.32 (45) CMGLSP 0.38 (44) 0.52 (44) 0.55 (45) 0.56 (44) 0.00 (44) 0.01 (43) 0.00 (43) 0.00 (45) 0.05 (44) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.15 (37) NA (NA) NA (NA) MCD12Q2 0.45 (37) 0.69 (37) 0.78 (37) 0.51 (37) 0.00 (37) 0.00 (37) 0.00 (37) 0.05 (37) 0.40 (37) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (34) NA (NA) NA (NA) NPN 0.07 (45) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.41 (45) 0.64 (45) 0.69 (45) 0.50 (45) 0.04 (45) 0.24 (45) 0.03 (45) 0.12 (45) 0.48 (45) VIPPHEN_NDVI 0.35 (45) 0.50 (45) 0.53 (45) 0.68 (45) 0.00 (45) 0.18 (45) 0.00 (45) 0.05 (45) 0.43 (45) DLC 0.57 (45) 0.69 (45) 0.58 (45) 0.72 (45) 0.07 (45) 0.28 (45) 0.43 (45) 0.04 (45) 0.02 (45) eMODIS 0.48 (42) 0.60 (42) 0.56 (44) 0.60 (44) 0.08 (42) 0.18 (42) 0.07 (44) 0.14 (45) 0.13 (44) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 7.17 (41) 4.32 (41) 15.56 (45) 62.47 (38) 24.98 (41) 74.35 (34) CMGLSP 9.68 (44) 3.02 (44) 3.22 (45) 18.93 (44) 11.32 (44) 11.40 (43) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 8.57 (37) 1.11 (37) 7.35 (37) 40.32 (37) 16.92 (37) 48.89 (37) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 20.18 (45) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 15.07 (45) 0.31 (45) 16.69 (45) 46.44 (45) 32.76 (45) 61.51 (45) VIPPHEN_NDVI 33.00 (45) 8.67 (45) 16.67 (45) 54.44 (45) 50.67 (45) 87.44 (45) DLC 16.04 (45) 24.19 (45) 19.71 (45) 25.04 (45) 36.76 (45) 41.09 (45) eMODIS 5.43 (42) 8.80 (42) 20.64 (44) 39.02 (44) 30.45 (42) 49.74 (42) Table A4. Shrublands (SH) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.01 (7) 0.07 (7) 0.00 (7) 0.20 (4) 0.03 (7) 0.40 (4) 0.02 (7) 0.16 (7) 0.08 (7) CMGLSP 0.22 (7) 0.28 (7) 0.28 (7) 0.62 (7) 0.12 (7) 0.20 (7) 0.01 (7) 0.19 (7) 0.10 (7) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (7) NA (NA) NA (NA) MCD12Q2 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.01 (7) NA (NA) NA (NA) NPN 0.13 (7) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.35 (7) 0.28 (7) 0.50 (7) 0.39 (7) 0.59 (7) 0.07 (7) 0.09 (7) 0.46 (7) 0.14 (7) VIPPHEN_NDVI 0.31 (7) 0.22 (7) 0.00 (7) 0.15 (7) 0.17 (7) 0.35 (7) 0.00 (7) 0.19 (7) 0.12 (7) DLC 0.55 (7) 0.53 (7) 0.43 (7) 0.54 (7) 0.44 (7) 0.26 (7) 0.00 (7) 0.14 (7) 0.07 (7) eMODIS 0.68 (7) 0.89 (7) 0.85 (7) 0.03 (6) 0.52 (7) 0.22 (6) 0.05 (7) 0.22 (7) 0.11 (7) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 38.43 (7) 30.64 (7) 23.86 (7) 24.50 (4) 13.57 (7) 34.00 (4) CMGLSP 25.43 (7) 6.79 (7) 10.86 (7) 8.71 (7) 35.29 (7) 16.71 (7) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 11.00 (1) 11.00 (1) 10.00 (1) 2.00 (1) 2.00 (1) 9.00 (1) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 108.57 (7) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 25.71 (7) 20.14 (7) 15.57 (7) 7.86 (7) 9.14 (7) 17.86 (7) VIPPHEN_NDVI 22.00 (7) 11.79 (7) 2.57 (7) 13.71 (7) 18.43 (7) 35.71 (7) DLC 7.29 (7) 10.79 (7) 6.14 (7) 24.00 (7) 0.14 (7) 31.29 (7) eMODIS 23.00 (7) 21.43 (7) 20.86 (7) 15.00 (6) 1.14 (7) 6.67 (6) Remote Sens. 2020, 12, 2538 18 of 27 Table A5. Deciduous/Broadleaf (DB) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.37 (23) 0.50 (23) 0.21 (27) 0.10 (26) 0.02 (23) 0.05 (23) 0.26 (25) 0.10 (27) 0.34 (27) CMGLSP 0.65 (24) 0.37 (23) 0.00 (26) 0.01 (26) 0.00 (23) 0.00 (23) 0.00 (23) 0.23 (26) 0.26 (26) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.24 (20) NA (NA) NA (NA) MCD12Q2 0.28 (26) 0.31 (26) 0.19 (26) 0.52 (26) 0.07 (26) 0.29 (26) 0.00 (26) 0.13 (26) 0.25 (26) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.90 (10) NA (NA) NA (NA) NPN 0.90 (27) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.04 (27) 0.45 (27) 0.13 (27) 0.21 (27) 0.17 (27) 0.01 (27) 0.54 (27) 0.00 (27) 0.33 (27) VIPPHEN_NDVI 0.59 (27) 0.33 (27) 0.07 (27) 0.21 (27) 0.07 (27) 0.13 (27) 0.59 (27) 0.16 (27) 0.25 (23) DLC 0.28 (27) 0.51 (27) 0.77 (27) 0.33 (27) 0.45 (27) 0.07 (27) 0.22 (27) 0.20 (27) 0.32 (27) eMODIS 0.02 (25) 0.36 (25) 0.26 (26) 0.59 (26) 0.02 (25) 0.03 (25) 0.24 (26) 0.05 (27) 0.59 (26) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 13.61 (23) 5.28 (23) 0.30 (27) 17.15 (26) 18.65 (23) 39.83 (23) CMGLSP 19.79 (24) 18.33 (23) 28.92 (26) 37.04 (26) 2.13 (23) 6.30 (23) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 15.58 (26) 7.46 (26) 1.65 (26) 4.08 (26) 18.23 (26) 19.65 (26) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 18.74 (27) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 17.56 (27) 8.39 (27) 1.78 (27) 1.15 (27) 20.33 (27) 16.41 (27) VIPPHEN_NDVI 49.63 (27) 27.22 (27) 3.81 (27) 10.63 (27) 46.81 (27) 60.26 (27) DLC 20.19 (27) 27.54 (27) 18.52 (27) 11.85 (27) 39.70 (27) 8.33 (27) eMODIS 8.08 (25) 0.08 (25) 8.92 (26) 11.15 (26) 18.00 (25) 21.48 (25) Table A6. Evergreen (EN) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.03 (25) 0.03 (25) 0.04 (25) 0.45 (9) 0.02 (25) 0.49 (9) 0.01 (25) 0.01 (25) 0.04 (25) CMGLSP 0.33 (25) 0.40 (25) 0.28 (25) 0.11 (25) 0.11 (25) 0.08 (25) 0.01 (25) 0.13 (25) 0.01 (25) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.04 (22) NA (NA) NA (NA) MCD12Q2 0.44 (11) 0.54 (11) 0.64 (11) 0.31 (11) 0.08 (11) 0.07 (11) 0.15 (11) 0.91 (11) 0.75 (11) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.02 (25) NA (NA) NA (NA) NPN 0.00 (25) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.16 (25) 0.06 (25) 0.01 (25) 0.00 (25) 0.04 (25) 0.06 (25) 0.06 (25) 0.01 (25) 0.09 (22) VIPPHEN_NDVI 0.24 (25) 0.19 (25) 0.04 (25) 0.17 (25) 0.12 (25) 0.09 (25) 0.06 (25) 0.04 (25) 0.52 (20) DLC 0.10 (25) 0.07 (25) 0.00 (25) 0.09 (25) 0.00 (25) 0.02 (25) 0.04 (25) 0.00 (25) 0.16 (25) eMODIS 0.09 (23) 0.04 (23) 0.00 (24) 0.09 (23) 0.00 (23) 0.00 (22) 0.07 (24) 0.00 (25) 0.11 (24) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 55.56 (25) 52.88 (25) 51.20 (25) 14.33 (9) 3.36 (25) 23.56 (9) CMGLSP 29.44 (25) 9.52 (25) 9.40 (25) 5.04 (25) 37.84 (25) 24.40 (25) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 32.64 (11) 22.68 (11) 13.73 (11) 11.36 (11) 17.91 (11) 44.00 (11) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 7.52 (25) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 24.52 (25) 15.70 (25) 7.88 (25) 18.84 (25) 15.64 (25) 5.68 (25) VIPPHEN_NDVI 3.96 (25) 26.14 (25) 49.32 (25) 39.36 (25) 46.36 (25) 35.40 (25) DLC 25.64 (25) 2.06 (25) 28.68 (25) 5.60 (25) 4.04 (25) 31.24 (25) eMODIS 32.35 (23) 44.93 (23) 47.17 (24) 10.17 (23) 27.17 (23) 8.32 (22) Remote Sens. 2020, 12, 2538 19 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 25 of 33 Appendix C Appendix C Additional figures relating to the long-term trend analysis of CMGLSP, VIP_EVI, and NPN. Additional figures relating to the long-term trend analysis of CMGLSP, VIP_EVI, and NPN. Figure A1. Mean and standard deviation of historical (1982–2016) data from CMGLSP and Figure A1. Mean and standard deviation of historical (1982–2016) data from CMGLSP and VIPPHEN_EVI2. VIPPHEN_EVI2. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 20 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 33 Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 33 Figure A2 Figure A2 . . Change in SOS da Change in SOS da te and te and variatio variatio n of NPN n of NPN from 1982 to from 1982 to 2016 2016 alongside alongside consist consist ee ncy with ncy with Figure A2. Change in SOS date and variation of NPN from 1982 to 2016 alongside consistency with CMG CMGLLSP and SP and VIPPH VIPPHEN_EVI EN_EVI2. N 2. Neeggaative tive va valu lues c es coorrespond to rrespond to earlier earlier da dates (or tes (or fewer fewer day dayss) and ) and CMGLSP and VIPPHEN_EVI2. Negative values correspond to earlier dates (or fewer days) and positive positiv positiv ee values to later date values to later date s ( s ( oor more days). Co r more days). Consistency shows dataset a nsistency shows dataset a gg reement with regards to reement with regards to values to later dates (or more days). Consistency shows dataset agreement with regards to the direction the direct the direct ion of ion of change. Change in date change. Change in date is is the r the r ee gre gre ssed trend ssed trend slope m slope m u u ltipl ltipl ie ie d d by the number of years. by the number of years. of change. Change in date is the regressed trend slope multiplied by the number of years. Change in Change in variation is the reg Change in variation is the reg rr essed essed trend sl trend sl ope of ope of the abso the abso lu lu te valu te valu e of th e of th e residu e residu al al s mu s mu ltiplie ltiplie d by d by variation is the regressed trend slope of the absolute value of the residuals multiplied by the number the number of years. the number of years. of years. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Rem Rem oo te Sens te Sens . .2020 2020 , , 12 12 , x; , x; doi: FOR doi: FOR PEER PEER RE RE VIEW VIEW www. www. mdpi.com/ mdpi.com/ journal/remotese journal/remotese nn si si ng ng Remote Sens. 2020, 12, 2538 21 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 27 of 33 Figure A4. Change in date and variation of CMGLSP, VIPPHEN_EVI2, and NPN across three phenology Figure A4. Change in date and variation of CMGLSP, VIPPHEN_EVI2, and NPN across three metrics (only SOS for NPN) from 1982 to 2016. Shown here with a more robust color scheme to better phenology metrics (only SOS for NPN) from 1982 to 2016. Shown here with a more robust color indicate change. Negative values correspond to earlier dates (or fewer days) and positive values to scheme to better indicate change. Negative values correspond to earlier dates (or fewer days) and later dates (or more days). Change in date is the regressed trend slope multiplied by the number of positive values to later dates (or more days). Change in date is the regressed trend slope multiplied years. Change in variation is the regressed trend slope of the absolute value of the residuals multiplied by the number of years. Change in variation is the regressed trend slope of the absolute value of the by the number of years. residuals multiplied by the number of years. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 22 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 28 of 33 Remote Sens. 2020, 12, x FOR PEER REVIEW 28 of 33 Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 23 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 29 of 33 Figure A7. Histogram of CMGLSP and VIPPHEN_EVI2 EOS change in date from 1982 to 2016. Figure A7. Histogram of CMGLSP and VIPPHEN_EVI2 EOS change in date from 1982 to 2016. References References 1. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, 1. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. 2013, 169, 156–173. [CrossRef] 2013, 169, 156–173. 2. Chuine, I. Why does phenology drive species distribution? Philos. Trans. R. Soc. B Biol. Sci. 2010, 365, 2. Chuine, I. Why does phenology drive species distribution? Philos. Trans. R. Soc. B Biol. Sci. 2010, 365, 3149– 3149–3160. [CrossRef] [PubMed] 3. Schwartz, M.D. Green-wave phenology. Nature 1998, 394, 839. [CrossRef] 3. Schwartz, M.D. Green-wave phenology. Nature 1998, 394, 839. 4. Hurley, M.A.; Hebblewhite, M.; Gaillard, J.M.; Dray, S.; Taylor, K.A.; Smith, W.K.; Zager, P.; Bonenfant, C. 4. Hurley, M.A.; Hebblewhite, M.; Gaillard, J.M.; Dray, S.; Taylor, K.A.; Smith, W.K.; Zager, P.; Bonenfant, C. Functional analysis of normalized di erence vegetation index curves reveals overwinter mule deer survival Functional analysis of normalized difference vegetation index curves reveals overwinter mule deer is driven by both spring and autumn phenology. Philos. Trans. R. Soc. B Biol. Sci. 2014, 369. [CrossRef] survival is driven by both spring and autumn phenology. Philos. Trans. R. Soc. B Biol. Sci. 2014, 369, 5. Post, E.; Forchhammer, M.C. Climate change reduces reproductive success of an Arctic herbivore through doi:10.1098/rstb.2013.0196. trophic mismatch. Philos. Trans. R. Soc. B Biol. Sci. 2008, 363, 2369–2375. [CrossRef] 5. Post, E.; Forchhammer, M.C. Climate change reduces reproductive success of an Arctic herbivore through 6. Rickbeil, G.J.M.; Merkle, J.A.; Anderson, G.; Atwood, M.P.; Beckmann, J.P.; Cole, E.K.; Courtemanch, A.B.; trophic mismatch. Philos. Trans. R. Soc. B Biol. Sci. 2008, 363, 2369–2375. Dewey, S.; Gustine, D.D.; Kau man, M.J.; et al. Plasticity in elk migration timing is a response to changing 6. Rickbeil, G.J.M.; Merkle, J.A.; Anderson, G.; Atwood, M.P.; Beckmann, J.P.; Cole, E.K.; Courtemanch, A.B.; environmental conditions. Glob. Chang. Biol. 2019, 25, 2368–2381. [CrossRef] Dewey, S.; Gustine, D.D.; Kauffman, M.J.; et al. Plasticity in elk migration timing is a response to changing 7. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response to environmental conditions. Glob. Chang. Biol. 2019, 25, 2368–2381. global change. Trends Ecol. Evol. 2007, 22, 357–365. [CrossRef] 7. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response 8. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.M.; Tucker, C.J.; Stenseth, N.C. Using the satellite-derived to global change. Trends Ecol. Evol. 2007, 22, 357–365. NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 2005, 20, 503–510. [CrossRef] 8. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.M.; Tucker, C.J.; Stenseth, N.C. Using the satellite-derived 9. Zhu, W.; Tian, H.; Xu, X.; Pan, Y.; Chen, G.; Lin, W. Extension of the growing season due to delayed autumn NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 2005, 20, 503–510. over mid and high latitudes in North America during 1982–2006. Glob. Ecol. Biogeogr. 2012, 21, 260–271. 9. Zhu, W.; Tian, H.; Xu, X.; Pan, Y.; Chen, G.; Lin, W. Extension of the growing season due to delayed autumn [CrossRef] over mid and high latitudes in North America during 1982–2006. Glob. Ecol. Biogeogr. 2012, 21, 260–271. 10. Cook, B.I.; Wolkovich, E.M.; Parmesan, C. Divergent responses to spring and winter warming drive 10. Cook, B.I.; Wolkovich, E.M.; Parmesan, C. Divergent responses to spring and winter warming drive community level flowering trends. Proc. Natl. Acad. Sci. USA 2012, 109, 9000–9005. [CrossRef] community level flowering trends. Proc. Natl. Acad. Sci. USA 2012, 109, 9000–9005. 11. Brown, J.F.; Ji, L.; Gallant, A.; Kau man, M. Exploring relationships of spring green-up to moisture and 11. Brown, J.F.; Ji, L.; Gallant, A.; Kauffman, M. Exploring relationships of spring green-up to moisture and temperature across Wyoming, USA. Int. J. Remote Sens. 2019, 40, 956–984. [CrossRef] temperature across Wyoming, USA. Int. J. Remote Sens. 2019, 40, 956–984. 12. Riotte-Lambert, L.; Matthiopoulos, J. Environmental Predictability as a Cause and Consequence of Animal 12. Riotte-Lambert, L.; Matthiopoulos, J. Environmental Predictability as a Cause and Consequence of Animal Movement. Trends Ecol. Evol. 2020, 35, 163–174. [CrossRef] [PubMed] Movement. Trends Ecol. Evol. 2020, 35, 163–174. 13. St-Louis, V.; Pidgeon, A.M.; Clayton, M.K.; Locke, B.A.; Bash, D.; Radeloff, V.C. Satellite image texture and a vegetation index predict avian biodiversity in the Chihuahuan Desert of New Mexico. Ecography 2009, 32, 468–480. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 24 of 27 13. St-Louis, V.; Pidgeon, A.M.; Clayton, M.K.; Locke, B.A.; Bash, D.; Radelo , V.C. Satellite image texture and a vegetation index predict avian biodiversity in the Chihuahuan Desert of New Mexico. Ecography 2009, 32, 468–480. [CrossRef] 14. Hurlbert, A.H. Species-energy relationships and habitat complexity in bird communities. Ecol. Lett. 2004, 7, 714–720. [CrossRef] 15. Searle, K.R.; Rice, M.B.; Anderson, C.R.; Bishop, C.; Hobbs, N.T. Asynchronous vegetation phenology enhances winter body condition of a large mobile herbivore. Oecologia 2015, 179, 377–391. [CrossRef] [PubMed] 16. Pettorelli, N.; Pelletier, F.; Von Hardenberg, A.; Festa-Bianchet, M.; Côté, S.D. Early onset of vegetation growth vs. rapid green-up: Impacts on juvenile mountain ungulates. Ecology 2007, 88, 381–390. [CrossRef] 17. Hebblewhite, M.; Merrill, E.; McDermid, G. A Multi-Scale Test Of The Forage Maturation Hypothesis In A Partially Migratory Ungulate Population. Ecol. Monogr. 2008, 78, 141–166. [CrossRef] 18. Middleton, A.D.; Kau man, M.J.; Mcwhirter, D.E.; Cook, J.G.; Cook, R.C.; Nelson, A.A.; Jimenez, M.D.; Klaver, R.W. Animal migration amid shifting patterns of phenology and predation: Lessons from a Yellowstone elk herd. Ecology 2013, 94, 1245–1256. [CrossRef] 19. Stoner, D.C.; Sexton, J.O.; Choate, D.M.; Nagol, J.; Bernales, H.H.; Sims, S.A.; Ironside, K.E.; Longshore, K.M.; Edwards, T.C. Climatically driven changes in primary production propagate through trophic levels. Glob. Chang. Biol. 2018, 24, 4453–4463. [CrossRef] 20. Monteith, K.L.; Klaver, R.W.; Hersey, K.R.; Holland, A.A.; Thomas, T.P.; Kau man, M.J. E ects of climate and plant phenology on recruitment of moose at the southern extent of their range. Oecologia 2015, 178, 1137–1148. [CrossRef] 21. Middleton, A.D.; Merkle, J.A.; McWhirter, D.E.; Cook, J.G.; Cook, R.C.; White, P.J.; Kau man, M.J. Green-wave surfing increases fat gain in a migratory ungulate. Oikos 2018, 127, 1060–1068. [CrossRef] 22. Hamel, S.; Garel, M.; Festa-Bianchet, M.; Gaillard, J.M.; Côté, S.D. Spring Normalized Di erence Vegetation Index (NDVI) predicts annual variation in timing of peak faecal crude protein in mountain ungulates. J. Appl. Ecol. 2009, 46, 582–589. [CrossRef] 23. Bischof, R.; Loe, L.E.; Meisingset, E.L.; Zimmermann, B.; Van Moorter, B.; Mysterud, A. A Migratory Northern Ungulate in the Pursuit of Spring: Jumping or Surfing the Green Wave? Am. Nat. 2012, 180, 407–424. [CrossRef] [PubMed] 24. Merkle, J.A.; Monteith, K.L.; Aikens, E.O.; Hayes, M.M.; Hersey, K.R.; Middleton, A.D.; Oates, B.A.; Sawyer, H.; Scurlock, B.M.; Kau man, M.J. Large herbivores surf waves of green-up during spring. Proc. R. Soc. B Biol. Sci. 2016, 283, 20160456. [CrossRef] 25. Aikens, E.O.; Kau man, M.J.; Merkle, J.A.; Dwinnell, S.P.H.; Fralick, G.L.; Monteith, K.L. The greenscape shapes surfing of resource waves in a large migratory herbivore. Ecol. Lett. 2017, 20, 741–750. [CrossRef] 26. Post, E.; Pedersen, C.; Wilmers, C.C.; Forchhammer, M.C. Warming, plant phenology and the spatial dimension of trophic mismatch for large herbivores. Proc. R. Soc. B Biol. Sci. 2008, 275, 2005–2013. [CrossRef] 27. Mikle, N.L.; Graves, T.A.; Olexa, E.M. To forage or flee: Lessons from an elk migration near a protected area. Ecosphere 2019, 10, e02693. [CrossRef] 28. Rivrud, I.M.; Bischof, R.; Meisingset, E.L.; Zimmermann, B.; Loe, L.E.; Mysterud, A. Leave before it’s too late: Anthropogenic and environmental triggers of autumn migration in a hunted ungulate population. Ecology 2016, 97, 1058–1068. [CrossRef] 29. Mueller, T.; Olson, K.A.; Fuller, T.K.; Schaller, G.B.; Murray, M.G.; Leimgruber, P. In search of forage: Predicting dynamic habitats of Mongolian gazelles using satellite-based estimates of vegetation productivity. J. Appl. Ecol. 2008, 45, 649–658. [CrossRef] 30. Pettorelli, N.; Gaillard, J.M.; Mysterud, A.; Duncan, P.; Stenseth, N.C.; Delorme, D.; Van Laere, G.; Toïgo, C.; Klein, F. Using a proxy of plant productivity (NDVI) to find key periods for animal performance: The case of roe deer. Oikos 2006, 112, 565–572. [CrossRef] 31. Zhang, X.; Tarpley, D.; Sullivan, J.T. Diverse responses of vegetation phenology to a warming climate. Geophys. Res. Lett. 2007, 34, 1–5. [CrossRef] 32. Peng, D.; Zhang, X.; Wu, C.; Huang, W.; Gonsamo, A.; Huete, A.R.; Didan, K.; Tan, B.; Liu, X.; Zhang, B. Intercomparison and evaluation of spring phenology products using National Phenology Network and AmeriFlux observations in the contiguous United States. Agric. For. Meteorol. 2017, 242, 33–46. [CrossRef] Remote Sens. 2020, 12, 2538 25 of 27 33. Keenan, T.F.; Richardson, A.D. The timing of autumn senescence is a ected by the timing of spring phenology: Implications for predictive models. Glob. Chang. Biol. 2015, 21, 2634–2641. [CrossRef] 34. Schwartz, M.D.; Reed, B.C.; White, M.A. Assesing satellite-derived start-of-season measures in the conterminous USA. Int. J. Climatol. 2002, 22, 1793–1805. [CrossRef] 35. Richardson, A.D.; Hufkens, K.; Milliman, T.; Frolking, S. Intercomparison of phenological transition dates derived from the PhenoCam Dataset V1.0 and MODIS satellite remote sensing. Sci. Rep. 2018, 8, 1–12. 36. Klosterman, S.T.; Hufkens, K.; Gray, J.M.; Melaas, E.; Sonnentag, O.; Lavine, I.; Mitchell, L.; Norman, R. Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery. Biogeosciences 2014, 11, 4305–4320. [CrossRef] 37. White, M.A.; de Beurs, K.M.; Didan, K.; Inouye, D.W.; Richardson, A.D.; Jensen, O.P.; O’Keefe, J.; Zhang, G.; Nemani, R.R.; van Leeuwen, W.J.D.; et al. Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982–2006. Glob. Chang. Biol. 2009, 15, 2335–2359. [CrossRef] 38. Richardson, A.D.; Hufkens, K.; Milliman, T.; Aubrecht, D.M.; Chen, M.; Gray, J.M.; Johnston, M.R.; Keenan, T.F.; Klosterman, S.T.; Kosmala, M.; et al. PhenoCam Dataset v1.0: Vegetation Phenology from Digital Camera Imagery, 2000–2015. 2017. [CrossRef] 39. Omernik, J.M.; Grith, G.E. Ecoregions of the Conterminous United States: Evolution of a Hierarchical Spatial Framework. Environ. Manag. 2014, 54, 1249–1266. [CrossRef] 40. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.; Gao, X.; Ferreira, L. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [CrossRef] 41. Crimmins, T.M.; Marsh, R.L.; Switzer, J.R.; Crimmins, M.A.; Gerst, K.L.; Rosemartin, A.H.; Weltzin, J.F.; Jewell, S. USA National Phenology Network Gridded Products Documentation; United States Geological Survey: Reston, VA, USA, 2017; p. 34. 42. Robinson, N.P.; Allred, B.W.; Smith, W.K.; Jones, M.O.; Moreno, A.; Erickson, T.A.; Naugle, D.E.; Running, S.W. Terrestrial primary production for the conterminous United States derived from Landsat 30 m and MODIS 250 m. Remote Sens. Ecol. Conserv. 2018, 4, 264–280. [CrossRef] 43. Meier, G.A.; Brown, J.F. Remote Sensing of Land Surface Phenology; United States Geological Survey: Reston, VA, USA, 2014. 44. Jenkerson, C.B.; Maiersperger, T.; Schmidt, G. eMODIS: A User-Friendly Data Source; United States Geological Survey: Reston, VA, USA, 2010. 45. Friedl, M.; Gray, J.; Sulla-Menashe, D. MCD12Q2 MODIS/Terra+Aqua Land Cover Dynamics Yearly L3 Global 500m SIN Grid V006 [Data Set]; NASA: Washington, DC, USA, 2019. 46. Zhang, X. Reconstruction of a complete global time series of daily vegetation index trajectory from long-term AVHRR data. Remote Sens. Environ. 2015, 156, 457–472. [CrossRef] 47. Didan, K.; Munoz, A.B.; Miura, T.; Tsend-Ayush, J.; Zhang, X.; Friedl, M.; Gray, J.; Van Leeuwen, W.; Czapla-Myers, J.; Jenkerson, C.; et al. Multi-Sensor Vegetation Index and Phenology Earth Science Data Records Algorithm Theoretical Basis Document and User Guide; Vegetation Index and Phenology Lab, University of Arizona: Tucson, AZ, USA, 2015. 48. Richardson, A.D.; Hufkens, K.; Milliman, T.; Aubrecht, D.M.; Chen, M.; Gray, J.M.; Johnston, M.R.; Keenan, T.F.; Klosterman, S.T.; Kosmala, M.; et al. Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery. Sci. Data 2018, 5, 1–24. [CrossRef] [PubMed] 49. Hufkens, K.; Basler, D.; Milliman, T.; Melaas, E.K.; Richardson, A.D. An integrated phenology modelling framework in r. Methods Ecol. Evol. 2018, 9, 1276–1285. [CrossRef] 50. Sonnentag, O.; Hufkens, K.; Teshera-Sterne, C.; Young, A.M.; Friedl, M.; Braswell, B.H.; Milliman, T.; O’Keefe, J.; Richardson, A.D. Digital repeat photography for phenological research in forest ecosystems. Agric. For. Meteorol. 2012, 152, 159–177. [CrossRef] 51. Melaas, E.K.; Sulla-Menashe, D.; Gray, J.M.; Black, T.A.; Morin, T.H.; Richardson, A.D.; Friedl, M.A. Multisite analysis of land surface phenology in North American temperate and boreal deciduous forests from Landsat. Remote Sens. Environ. 2016, 186, 452–464. [CrossRef] 52. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O’Connor, B.; Paganini, M. Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island. Remote Sens. Environ. 2018, 215, 517–529. [CrossRef] Remote Sens. 2020, 12, 2538 26 of 27 53. Sen, P.K. Estimates of the Regression Coecient Based on Kendall’s Tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [CrossRef] 54. De Beurs, K.M.; Henebry, G.M. Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in Kazakhstan. Remote Sens. Environ. 2004, 89, 497–509. [CrossRef] 55. R Core Team, R. A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019. 56. Marchetto, A. rkt: Mann-Kendall Test, Seasonal and Regional Kendall Tests; R Package Version 1.5; Verbania Pallanza, Italy, 2017. 57. Geremia, C.; Merkle, J.A.; Eacker, D.R.; Wallen, R.L.; White, P.J.; Hebblewhite, M.; Kau man, M.J. Migrating bison engineer the green wave. Proc. Natl. Acad. Sci. USA 2019, 116, 25707–25713. [CrossRef] 58. Filippa, G.; Cremonese, E.; Migliavacca, M.; Galvagno, M.; Sonnentag, O.; Humphreys, E.; Hufkens, K.; Ryu, Y.; Verfaillie, J.; di Cella, U.M.; et al. NDVI derived from near-infrared-enabled digital cameras: Applicability across di erent plant functional types. Agric. For. Meteorol. 2018, 249, 275–285. [CrossRef] 59. Liu, Y.; Wu, C.; Sonnentag, O.; Desai, A.R.; Wang, J. Using the red chromatic coordinate to characterize the phenology of forest canopy photosynthesis. Agric. For. Meteorol. 2020, 285–286, 107910. [CrossRef] 60. Borowik, T.; Pettorelli, N.; Sönnichsen, L.; Jedr ˛ zejewska, B. Normalized di erence vegetation index (NDVI) as a predictor of forage availability for ungulates in forest and field habitats. Eur. J. Wildl. Res. 2013, 59, 675–682. [CrossRef] 61. Pettorelli, N.; Ryan, S.; Mueller, T.; Bunnefeld, N.; Jedrzejewska, B.; Lima, M.; Kausrud, K. The Normalized Di erence Vegetation Index (NDVI): Unforeseen successes in animal ecology. Clim. Res. 2011, 46, 15–27. [CrossRef] 62. Zhang, X.; Liu, L.; Yan, D. Comparisons of global land surface seasonality and phenology derived from AVHRR, MODIS, and VIIRS data. J. Geophys. Res. Biogeosci. 2017, 122, 1506–1525. [CrossRef] 63. Lendrum, P.E.; Anderson, C.R.; Monteith, K.L.; Jenks, J.A.; Bowyer, R.T. Relating the movement of a rapidly migrating ungulate to spatiotemporal patterns of forage quality. Mamm. Biol. 2014, 79, 369–375. [CrossRef] 64. Garroutte, E.L.; Hansen, A.J.; Lawrence, R.L. Using NDVI and EVI to map spatiotemporal variation in the biomass and quality of forage for migratory elk in the Greater Yellowstone Ecosystem. Remote Sens. 2016, 8, 404. [CrossRef] 65. Myneni, R.B.; Keeling, C.D.; Tucker, C.J.; Asrar, G.; Nemani, R.R. Increased plant growth in the northern high latitudes from 1981 to 1991. Nature 1997, 386, 698–702. [CrossRef] 66. Zhou, L.; Tucker, C.J.; Kaufmann, R.K.; Slayback, D.; Shabanov, N.V.; Myneni, R.B. Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981–1999. J. Geophys. Res. 2001, 106, 20069–20083. [CrossRef] 67. Berman, E.E.; Bolton, D.K.; Coops, N.C.; Mityok, Z.K.; Stenhouse, G.B.; Moore, R.D. Daily estimates of Landsat fractional snow cover driven by MODIS and dynamic time-warping. Remote Sens. Environ. 2018, 216, 635–646. [CrossRef] 68. Bradley, B.A.; Mustard, J.F. Identifying land cover variability distinct from land cover change: Cheatgrass in the Great Basin. Remote Sens. Environ. 2005, 94, 204–213. [CrossRef] 69. Bradley, B.A.; Mustard, J.F. Comparison of phenology trends by land cover class: A case study in the Great Basin, USA. Glob. Chang. Biol. 2008, 14, 334–346. [CrossRef] 70. Van Leeuwen, W.J.D.; Orr, B.J.; Marsh, S.E.; Herrmann, S.M. Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications. Remote Sens. Environ. 2006, 100, 67–81. [CrossRef] 71. Luo, Y.; El-Madany, T.S.; Filippa, G.; Ma, X.; Ahrens, B.; Carrara, A.; Gonzalez-Cascon, R.; Cremonese, E.; Galvagno, M.; Hammer, T.W.; et al. Using near-infrared-enabled digital repeat photography to track structural and physiological phenology in Mediterranean tree-grass ecosystems. Remote Sens. 2018, 10, 1293. [CrossRef] 72. Morisette, J.T.; Richardson, A.D.; Knapp, A.K.; Fisher, J.I.; Graham, E.A.; Abatzoglou, J.; Wilson, B.E.; Breshears, D.D.; Henebry, G.M.; Hanes, J.M.; et al. Tracking the rhythm of the seasons in the face of global change: Phenological research in the 21st century. Front. Ecol. Environ. 2009, 7, 253–260. [CrossRef] 73. Browning, D.M.; Rango, A.; Karl, J.W.; Laney, C.M.; Vivoni, E.R.; Tweedie, C.E. Emerging technological and cultural shifts advancing drylands research and management. Front. Ecol. Environ. 2015, 13, 52–60. [CrossRef] Remote Sens. 2020, 12, 2538 27 of 27 74. Crimmins, T.M.; Crimmins, M.A.; Bertelsen, D.; Balmat, J. Relationships between alpha diversity of plant species in bloom and climatic variables across an elevation gradient. Int. J. Biometeorol. 2008, 52, 353–366. [CrossRef] 75. Zhang, X.; Goldberg, M.; Tarpley, D.; Friedl, M.A.; Morisette, J.; Kogan, F.; Yu, Y. Drought-induced vegetation stress in southwestern North America. Environ. Res. Lett. 2010, 5, 024008. [CrossRef] 76. Hufkens, K.; Friedl, M.; Sonnentag, O.; Braswell, B.H.; Milliman, T.; Richardson, A.D. Linking near-surface and satellite remote sensing measurements of deciduous broadleaf forest phenology. Remote Sens. Environ. 2012, 117, 307–321. [CrossRef] 77. Zhang, X.; Wang, J.; Gao, F.; Liu, Y.; Schaaf, C.; Friedl, M.; Yu, Y.; Jayavelu, S.; Gray, J.; Liu, L.; et al. Exploration of scaling e ects on coarse resolution land surface phenology. Remote Sens. Environ. 2017, 190, 318–330. [CrossRef] 78. Liu, L.; Cao, R.; Shen, M.; Chen, J.; Wang, J.; Zhang, X. How does scale e ect influence spring vegetation phenology estimated from satellite-derived vegetation indexes? Remote Sens. 2019, 11, 2137. [CrossRef] 79. Baumann, M.; Ozdogan, M.; Richardson, A.D.; Radelo , V.C. Phenology from Landsat when data is scarce: Using MODIS and Dynamic Time-Warping to combine multi-year Landsat imagery to derive annual phenology curves. Int. J. Appl. Earth Obs. Geoinf. 2017, 54, 72–83. [CrossRef] 80. McClelland, C.J.R.; Coops, N.C.; Berman, E.E.; Kearney, S.P.; Nielsen, S.E.; Burton, A.C.; Stenhouse, G.B. Detecting changes in understorey and canopy vegetation cycles in West Central Alberta using a fusion of Landsat and MODIS. Appl. Veg. Sci. 2019, 23, 223–238. [CrossRef] 81. Bolton, D.K.; Gray, J.M.; Melaas, E.K.; Moon, M.; Eklundh, L.; Friedl, M.A. Continental-scale land surface phenology from harmonized Landsat 8 and Sentinel-2 imagery. Remote Sens. Environ. 2020, 240, 111685. [CrossRef] 82. Rodriguez-Galiano, V.F.; Dash, J.; Atkinson, P.M. Intercomparison of satellite sensor land surface phenology and ground phenology in Europe. Geophys. Res. Lett. 2015, 42, 2253–2260. [CrossRef] 83. Czernecki, B.; Nowosad, J.; Jabłonska, ´ K. Machine learning modeling of plant phenology based on coupling satellite and gridded meteorological dataset. Int. J. Biometeorol. 2018, 62, 1297–1309. [CrossRef] 84. Tang, J.; Körner, C.; Muraoka, H.; Piao, S.; Shen, M.; Thackeray, S.J.; Yang, X. Emerging opportunities and challenges in phenology: A review. Ecosphere 2016, 7, 1–17. [CrossRef] 85. Johnston, A.N.; Beever, E.A.; Merkle, J.A.; Chong, G. Vegetation responses to sagebrush-reduction treatments measured by satellites. Ecol. Indic. 2018, 87, 66–76. [CrossRef] 86. Schwartz, M.D. (Ed.); Phenology: An. Integrative Environmental Science; Springer: Dordrecht, The Netherlands, © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Remote Sensing Unpaywall

Comparative Quality and Trend of Remotely Sensed Phenology and Productivity Metrics across the Western United States

Loading next page...
 
/lp/unpaywall/comparative-quality-and-trend-of-remotely-sensed-phenology-and-7SUqVPPdm0

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
2072-4292
DOI
10.3390/rs12162538
Publisher site
See Article on Publisher Site

Abstract

remote sensing Article Comparative Quality and Trend of Remotely Sensed Phenology and Productivity Metrics across the Western United States 1 , 1 1 2 Ethan E. Berman * , Tabitha A. Graves , Nate L. Mikle , Jerod A. Merkle , 3 3 Aaron N. Johnston and Geneva W. Chong U.S. Geological Survey, Northern Rocky Mountain Science Center, Glacier Field Station, 38 Mather Drive, West Glacier, MT 59936, USA; tgraves@usgs.gov (T.A.G.); natemikle@gmail.com (N.L.M.) Department of Zoology and Physiology, University of Wyoming, Department 3166, 1000 E University Ave, Laramie, WY 82071, USA; jmerkle@uwyo.edu U.S. Geological Survey, Northern Rocky Mountain Science Center, 2327 University Way, Suite 2, Bozeman, MT 59715, USA; ajohnston@usgs.gov (A.N.J.); Geneva_chong@usgs.gov (G.W.C.) * Correspondence: bermane@gmail.com Received: 23 June 2020; Accepted: 30 July 2020; Published: 7 August 2020 Abstract: Vegetation phenology and productivity play a crucial role in surface energy balance, plant and animal distribution, and animal movement and habitat use and can be measured with remote sensing metrics including start of season (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), and integrated vegetation indices. However, for most metrics, we do not yet understand the agreement of remotely sensed data products with near-surface observations. We also need summaries of changes over time, spatial distribution, variability, and consistency in remote sensing dataset metrics for vegetation timing and quality. We compare metrics from 10 leading remote sensing datasets against a network of PhenoCam near-surface cameras throughout the western United States from 2002 to 2014. Most phenology metrics representing a date (SOS, PIRGd, POS, and EOS), rather than a duration (length of spring, length of growing season), better agreed with near-surface metrics but results varied by dataset, metric, and land cover, with absolute value of mean bias ranging from 0.38 (PIRGd) to 37.92 days (EOS). Datasets had higher agreement with PhenoCam metrics in shrublands, grasslands, and deciduous forests than in evergreen forests. Phenology metrics had higher agreement than productivity metrics, aside from a few datasets in deciduous forests. Using two datasets covering the period 1982–2016 that best agreed with PhenoCam metrics, we analyzed changes over time to growing seasons. Both datasets exhibited substantial spatial heterogeneity in the direction of phenology trends. Variability of metrics increased over time in some areas, particularly in the Southwest. Approximately 60% of pixels had consistent trend direction between datasets for SOS, POS, and EOS, with the direction varying by location. In all ecoregions except Mediterranean California, EOS has become later. This study comprehensively compares remote sensing datasets across multiple growing season metrics and discusses considerations for applied users to inform their data choices. Keywords: phenology; productivity; forage; green up; growing season; phenocam; timing; land surface phenology; elk; deer 1. Introduction Vegetation phenology and productivity are important drivers of ecosystem function by influencing processes as varied as surface energy balance [1] and plant species distribution [2]. Shifts in the timing of plant seasonality occur largely due to annual weather patterns and climate change [3], Remote Sens. 2020, 12, 2538; doi:10.3390/rs12162538 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 2 of 27 and these dynamics have consequences throughout terrestrial ecosystems. For example, phenology and productivity patterns strongly influence animal behavior, survival, and population dynamics [4]. Shifts in the timing and amount of vegetation can lead to trophic mismatch [5] and put stress on migratory species resilience and adaptive capacity [6]. Thus, to understand changes to the scale, rate, spatial configuration, and variability of ecological processes, we need to accurately measure vegetation phenology and productivity at broad spatial and temporal scales [7]. In recent decades, land surface phenology (LSP), the study of vegetation phenology and productivity from remote sensing, has revolutionized our understanding of ecological responses to phenological change [8]. LSP observations broaden the spatial scale of data to previously unattainable extents, enabling analyses of regional and continental vegetation patterns over time, with certain datasets extending back more than 30 years. Studies of LSP metrics in recent decades in North America have shown trends toward later senescence in the fall but inconclusive evidence for earlier spring green-up [9]. Green-up trends within ecological communities are complex, with plant species in the same community often showing opposite responses in both sign and magnitude to general warming patterns [10]. These trend analyses, along with annual LSP metrics, provide users (researchers, biologists, ecologists, natural resource managers, etc.) from a variety of fields with important insights into ecosystem processes and changing landscape dynamics driven by weather, climate, and human and natural disturbances. However, few studies have provided regional summaries of change that are accessible to users who are managing resources [7,11] or provided information on trends in variability of phenology over time. Insights into changing variability is important for wildlife researchers, as environmental predictability may shift habitat use, population density, and movement patterns [12]. Wildlife biologists began using LSP metrics in the early 2000s to better understand the influence of vegetation on the dynamics and distribution of animal populations including birds [13,14] and ungulates such as elk and deer [4,8,15–19]. For example, the timing of spring has implications for the fitness and body condition of ungulates [16,20–22] and peak instantaneous rate of green-up date (PIRGd), the date half way between start of season (SOS) and peak of season (POS) and representative of optimal forage quality, explains migratory patterns of ungulates surfing the green wave [23–25]. In addition, spatial heterogeneity of plant phenology, which may be declining due to warming temperatures, relates to the reproduction rates of caribou [26]. Autumn phenology, which has received considerably less attention than that of spring, can influence body mass and overwinter survival of mule deer [4] and migratory patterns of elk and red deer [18,27,28]. In addition to phenology, vegetation productivity is closely correlated to greenness indices such as Normalized Di erence Vegetation Index (NDVI) [8] and explains ungulate habitat use [24,29], health characteristics [30], and demographic parameters [16]. Changes in the seasonal timing of LSP metrics and the development of new datasets from various satellite platforms have received substantial focus from the remote sensing field [9,31–34]. However, despite the importance of LSP metrics to ecological applications, few studies have tested the quality of competing datasets in measuring LSP metrics against ground or near-surface observations or examined their relative agreement across various land cover types (but see [32,35]). Because LSP metrics synthesize information from millions of individual plants and multiple species, large biases across diverse vegetation communities can occur, and metrics may not directly represent the biological processes of the vegetation of interest [36]. Di erences in the processing algorithm of LSP data, such as the use of logistic curve fitting techniques versus splines, can also greatly impact performance [37]. Although users choose datasets based on their desired temporal and spatial scale, often several datasets exist with similar resolution yet di erent processing methodologies. A comparison of the quality of freely accessed and commonly derived LSP datasets against near-surface observations would assist users in selecting the best dataset for their application. Here, we provide a comparison of commonly used phenology and productivity metrics derived from 10 freely available LSP datasets. We examine correlation and bias in the western United States from 2002 to 2014 with near-surface observations from PhenoCam, a network of cameras spread Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 33 Remote Sens. 2020, 12, 2538 3 of 27 phenology and productivity metrics derived from remote sensing an indication of the strengths and weaknesses of different datasets, especially as related to land cover. We also compare long-term throughout North America providing data multiple times per day [38]. We provide users of phenology (1982–2016) trends in phenology from two leading LSP data products to identify spatial patterns, and productivity metrics derived from remote sensing an indication of the strengths and weaknesses assess the agreement about changing vegetation dynamics over time, and describe changing of di erent datasets, especially as related to land cover. We also compare long-term (1982–2016) trends variability in phenology. in phenology from two leading LSP data products to identify spatial patterns, assess the agreement about changing vegetation dynamics over time, and describe changing variability in phenology. Figure 1. The western United States (11 states) overlaid with Level 1 Ecoregions [39] and the locations and vegetation types of 29 PhenoCam sites. Note that some PhenoCam locations are slightly jittered to Figure 1. The western United States (11 states) overlaid with Level 1 Ecoregions [39] and the locations prevent overlap and therefore may not represent exact location (but are always in same ecoregion). and vegetation types of 29 PhenoCam sites. Note that some PhenoCam locations are slightly jittered to prevent overlap and therefore may not represent exact location (but are always in same ecoregion). 2. Materials and Methods We evaluated six phenology and three productivity metrics across the western United States, 2. Materials and Methods which represents a diverse range of ecosystems and land cover types (Figure 1). To evaluate agreement We evaluated six phenology and three productivity metrics across the western United States, with PhenoCam data, we focused on 10 datasets from 2002 to 2014 (Table 1), comparing overall which represents a diverse range of ecosystems and land cover types (Figure 1). To evaluate agreement and agreement within land cover types. We then assessed long-term trends using a subset agreement with PhenoCam data, we focused on 10 datasets from 2002 to 2014 (Table 1), comparing of more strongly associated datasets from 1982 to 2016. overall agreement and agreement within land cover types. We then assessed long-term trends using Phenology and productivity metrics can be derived from reflectance values recorded in optical a subset of more strongly associated datasets from 1982 to 2016. satellite imagery. Optical satellite imagery captures di erences in the reflectance of vegetation through Phenology and productivity metrics can be derived from reflectance values recorded in optical phenology cycles. Most commonly users measure vegetation indices as a ratio between reflectance in satellite imagery. Optical satellite imagery captures differences in the reflectance of vegetation the visual and near-infrared parts of the electromagnetic spectrum, with slightly di erent formulas through phenology cycles. Most commonly users measure vegetation indices as a ratio between for reflectance NDVI and in Enhanced the visual and near Vegetation Index -infrared parts (EVI; [40of ]). the el Fromectrom curvesag describing netic spectrum, wi the cycleth slightl of change, y users different for extract and mulas for ND evaluate key VI LSP and phenology Enhanced Ve and getatio productivity n Index (E metrics VI; [40] including ). From cuthose rves d evaluated escribing ther he e cycle of change, users extract and evaluate key LSP phenology and productivity metrics including (Figure 2): SOS, PIRGd, POS, end of season (EOS), length of spring (LOSp), length of growing season those evaluated here (Figure 2): SOS, PIRGd, POS, end of season (EOS), length of spring (LOSp), (LGS), integrated vegetation index (IVI), peak vegetation index (PVI), and amplitude of vegetation length of growing season (LGS), integrated vegetation index (IVI), peak vegetation index (PVI), and index (AVI). amplitude of vegetation index (AVI). Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 4 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 33 Figure 2. Six phenology and three productivity metrics shown through an example growing season Figure 2. Six phenology and three productivity metrics shown through an example growing season by by day of year (DOY). Note that these metrics are described through values of Normalized Difference day of year (DOY). Note that these metrics are described through values of Normalized Di erence Vegetation Index (NDVI)/Enhanced Vegetation Index (EVI) but certain datasets may produce metrics Vegetation Index (NDVI)/Enhanced Vegetation Index (EVI) but certain datasets may produce metrics using other methods/algorithms. using other methods/algorithms. Most datasets we evaluated use NDVI or EVI calculated from the surface reflectance of the Most datasets we evaluated use NDVI or EVI calculated from the surface reflectance of the Moderate Moderate Resolution Imaging Spectroradiometer (MODIS) or Advanced Very High Resolution Resolution Imaging Spectroradiometer (MODIS) or Advanced Very High Resolution Radiometer Radiometer (AVHRR) sensors (Table 1). We also evaluated one dataset, the National Phenology (AVHRR) sensors (Table 1). We also evaluated one dataset, the National Phenology Network first leaf Network first leaf spring index, hereinafter referred to as NPN, that models spatially explicit spring index, hereinafter referred to as NPN, that models spatially explicit temperature measurements temperature measurements parametrized via an extensive network of in situ phenological parametrized via an extensive network of in situ phenological observations [41]. This is the only observations [41]. This is the only dataset evaluated that does not incorporate optical satellite dataset evaluated that does not incorporate optical satellite imagery, but we included it because it is imagery, but we included it because it is readily available and provides annual SOS estimates across readily available and provides annual SOS estimates across the contiguous United States (CONUS). We the contiguous United States (CONUS). We also assessed two net primary productivity (NPP) also assessed two net primary productivity (NPP) datasets (CONUS Landsat NPP and CONUS MODIS datasets (CONUS Landsat NPP and CONUS MODIS NPP) that combine remote sensing with NPP) that combine remote sensing with meteorological data and plant physiological parameters [42]. meteorological data and plant physiological parameters [42]. Because they represent the total Because they represent the total productivity across the season, we include NPP in the IVI comparisons. productivity across the season, we include NPP in the IVI comparisons. AVI and PVI both measure AVI and PVI both measure the height of the curve at the peak, but AVI only includes the height from the height of the curve at the peak, but AVI only includes the height from the baseline to the peak. the baseline to the peak. Spatial resolutions spanned 30 m to approximately 5 km (0.05 degrees) and all Spatial resolutions spanned 30 m to approximately 5 km (0.05 degrees) and all data were freely available. Most datasets required minimal code and memory because the relevant phenology and data were freely available. Most datasets required minimal code and memory because the relevant productivity metrics were pre-processed and could be downloaded directly (e.g., SOS, POS, and EOS; phenology and productivity metrics were pre-processed and could be downloaded directly (e.g., SOS, Table 1). When necessary, we derived metrics (e.g., PIRGd and LOSp) from available metrics (see POS, and EOS; Table 1). When necessary, we derived metrics (e.g., PIRGd and LOSp) from available Appendix A Table A1 and Equations (A1)–(A7) for a complete description of metric availability, metrics (see Appendix A Table A1 and Equations (A1)–(A7) for a complete description of metric source, and derivation formula). We calculated all metrics for one dataset based on MODIS availability, source, and derivation formula). We calculated all metrics for one dataset based on MODIS MOD09Q1 NDVI (referred to as DLC for “double logistic curve”), fitting a double logistic curve and MOD09Q1 NDVI (referred to as DLC for “double logistic curve”), fitting a double logistic curve and accounting for snow cover following methods of [24], because of its use in many wildlife applications. accounting for snow cover following methods of [24], because of its use in many wildlife applications. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 5 of 27 Table 1. Abridged data table summarizing the 10 freely available datasets evaluated against PhenoCam. See Table A1 and Equations (A1)–(A7) for more information including methods used to derive individual metrics and data acquisition information. Significant Post-processing Dataset Base Input Spatial Resolution Temporal Coverage Accessible Online? Under Production? Spatial Coverage Reference Time Requirement? AVHRRP NDVI 1 km 1989–2014 yes no no CONUS [43] eMODIS NDVI 250 m 2001–2018 yes no yes CONUS [44] MCD12Q2 EVI 500 m 2001–2017 yes no yes Global [45] CMGLSP EVI 0.05 1982–2016 no no yes Global [46] VIPPHEN-EVI2 EVI 0.05 1981–2016 yes no no Global [47] VIPPHEN-NDVI NDVI 0.05 1981–2016 yes no no Global [47] NPN (First Leaf Spring Index) Temp, in situ observations 0.0417 1981–2019 yes no yes CONUS [41] Landsat NPP NDVI, meteorological inputs 30 m 1986–2018 yes no yes CONUS [42] MODIS NPP NDVI, meteorological inputs 250 m 2001–2018 yes no yes CONUS [42] DLC NDVI 250 m 2000–present yes yes yes Global [23] PhenoCam Green chromatic coordinate Point 2000–present yes no yes U.S. and Canada [48] 1 2 Collection 6 of eMODIS starts in 2003, but Collection 5 data is available from 2001. This analysis used Version 5 for 2002; Phenocam also has a limited number of sites in Panama, Hawaii, and Europe Remote Sens. 2020, 12, 2538 6 of 27 To compare datasets, we evaluated the agreement of LSP metrics from each dataset with metrics from PhenoCam data [48] at 29 sites (106 site-years) across the western United States (Figure 1) from 2002 to 2014. Using the R package phenocamr [49], we downloaded 3 day window PhenoCam data and calculated the 90th percentile green chromatic coordinate (GCC; Equation (1)). The 90th percentile value e ectively accounts for day-to-day variation due to weather and illumination patterns [50]. We defined rising and falling phenophases with a 10% amplitude threshold (to derive SOS and EOS) as described in PhenoCam documentation [48,50]. We chose the 10% threshold because it more closely resembles SOS than a 25% or 50% threshold and minimizes the bias between PhenoCam transition dates and MODIS transition dates [35]. GCC is a vegetation index derived from RGB camera imagery and is defined as: DNG GCC = (1) (DNR + DNG + DNB) where DNG, DNR, and DNB represent the digital number (i.e., strength) of the green, red, and blue channels, respectively. To assess the agreement of LSP metrics with PhenoCam, we compared correlation using the coecient of determination (R ) and the magnitude of di erence using mean bias. These common measures of statistical agreement have been previously used to compare LSP and near-surface phenology [36,51,52]. Mean bias is defined as: Mean bias = (est obs ) (2) i i i=1 where est and obs are the ith estimate (from LSP dataset) and observation (from PhenoCam) i i respectively. For productivity metrics, scales di er across datasets and thus we focused on correlation. We compared both overall agreement and agreement by landcover type, using the following vegetation types defined for PhenoCam sites (PhenoCam metadata: grasslands (GR, 45 site-years), shrublands (SH, 7 site-years), deciduous/broadleaf (DB, 27 site-years), evergreen (EN, 25 site-years), and wetlands (2 site-years). We excluded two sites in the Mediterranean California ecoregion that displayed a non-northern-temperate seasonal signal (SOS > 225) because green-up begins in late fall with rainfall and ends in early spring with drought. We excluded wetlands from the analysis by land cover classification due to the limited number of sites. We analyzed long-term trends (1982–2016) across the western U.S. using CMGLSP and VIPPHEN_EVI2, which agreed best with PhenoCam of the datasets that extended back to the 1980s. Of the six phenology metrics, we reported trends for SOS, POS, and EOS, because the other three are derived directly from these metrics and PIRGd spatial patterns are similar to SOS and POS. Using the Theil–Sen slope, the median of slopes between all pairs of observations within a pixel [53], we assessed each metric by pixel, including only pixels with data for at least 18 of 35 years. As a non-parametric test the Theil–Sen is more robust to outliers and provides higher statistical power when non-normality exists [54]. Negative slopes indicate change to earlier dates and positive slopes indicate change to later dates. We did not screen pixels for disturbances as the goal of this work is to identify patterns and magnitudes of change, rather than make inferences about the causes of phenological change. We evaluated overall variation within each pixel based on the standard deviation for each metric and whether variability has increased across years by applying the Theil–Sen slope to the absolute values of the residuals of the trend estimate against time. In addition, we measured consistency between the two datasets, defined as the by-pixel agreement that phenology dates are trending earlier or later. Lastly, to assess agreement in regional patterns of change, we report the proportion of area where both datasets agreed on direction of change by ecoregions using the United States Environmental Protection Agency Level 1 Ecoregions of North America dataset [39]. We used the statistical computing environment R for all analyses [55] and R package rkt [56] to calculate the Theil–Sen slope estimator. Remote Sens. 2020, 12, 2538 7 of 27 3. Results 3.1. Agreement of LSP metrics with PhenoCam We compared LSP metrics to 106 site-years of PhenoCam data from 29 sites. We found substantial variation in agreement between remotely sensed datasets and near-surface observations (Figure 3; Figure 4; Appendix B Tables A2–A6 for full results). Across datasets and six phenology metrics, R ranged from 0.03 to 0.37 (SOS), 0.20 to 0.55 (PIRGd), 0.25 to 0.54 (POS), 0.16 to 0.45 (EOS), 0.01 to 0.11 (LOSp), and 0 to 0.10 (LGS). Absolute value of mean bias (in days) ranged from 4.39 to 25.49 (SOS), 0.38 to 17.66 (PIRGd), 7.05 to 23.82 (POS), 0.52 to 37.92 (EOS), 12.45 to 44.66 (LOSp), 3.78 to 58.47 (LGS). Mean bias was lowest for SOS and PIRGd. Four datasets matched SOS and PIRGd from PhenoCam within 7 days (MCD12Q2, VIPPHEN_EVI2, DLC, and eMODIS for SOS; CMGLSP, MCD12Q2, VIPPHEN_EVI2, and VIPPHEN_NDVI for PIRGd). Generally, R values were higher for PIRGd and POS than SOS and EOS. Datasets were least correlated with PhenoCam transition dates for LOSp and LGS. In pairing a high R value with a low mean bias, MCD12Q2 and CMGLSP had best overall agreement with PhenoCam. Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 33 Figure 3. R and 2 mean bias (in days) for six phenology metrics classified by land cover type identified Figure 3. R and mean bias (in days) for six phenology metrics classified by land cover type identified in PhenoCam metadata and overall for all data points. The phenology metrics are start of season (SOS), in PhenoCam metadata and overall for all data points. The phenology metrics are start of season peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), length of (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), spring length of (LOSp) spring (LOSp) and length of growing seas and length of growing season (LGS). The on (LG classifications S). The classif areic grasslands ations are g (GR), rasslands shrublands (GR), (SH), shru deciduous blands (S /br H), oadleaf decidu (DB), ous/broadleaf (D evergreen (EN) B), e and verg overall reen (EN) and (OV). In over general, all (OV) this illustrates . In general, thi how the s overall illustrate evaluat s how ion the overal falls between l eval the uation fall landcover s be -specific tween thvalues e landcfor over-spec each dataset ific valu and es for metric eachcombination. dataset and metric combination. Phenology metrics had large di erences in correlation and bias by land cover type (Figure 3). Most datasets agreed moderately well with PhenoCam in grasslands, shrublands and deciduous/broadleaf forests for the four key phenology dates (SOS, PIRGd, POS, and EOS). Duration Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam data and overall for all data points. The productivity metrics are integrated vegetation index (IVI), peak vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 33 Remote Sens. 2020, 12, 2538 8 of 27 metrics (LOSp and LGS) correlated poorly with PhenoCam across all datasets. At grassland sites (n = 45), dates for most datasets correlated well with PhenoCam dates. At shrubland sites (n = 7), eMODIS was best for SOS, PIRGd, and POS, whereas CMGLSP was best for EOS. At deciduous/broadleaf sites (n = 27), CMGLSP and NPN agreed best for SOS, DLC for POS, and MCD12Q2 and eMODIS for EOS. Most datasets agreed well for PIRGd in deciduous/broadleaf forests. Datasets generally showed poor agreement with PhenoCam at evergreen forest sites (n = 25), with the exception of MCD12Q2. Figure 3. R and mean bias (in days) for six phenology metrics classified by land cover type identified Overall, datasets better represented phenology than productivity metrics (Figure 4). in PhenoCam metadata and overall for all data points. The phenology metrics are start of season Across datasets, R for the three productivity metrics ranged from 0.00 to 0.16 (IVI), 0.00 to 0.15 (SOS), peak instantaneous rate of green-up date (PIRGd), peak of season (POS), end of season (EOS), (PVI) and 0.00 to 0.34 (AVI), with low values in overall R masking high correlations in a few land length of spring (LOSp) and length of growing season (LGS). The classifications are grasslands (GR), cover types. For the IVI metric, MODIS_NPP had an R value of 0.9 at deciduous/broadleaf locations shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). In general, this and VIPPHEN_EVI2 and VIPPHEN_NDVI both agreed well. MCD12Q2 was highly correlated with illustrates how the overall evaluation falls between the landcover-specific values for each dataset and PhenoCam at evergreen locations for PVI and AVI. metric combination. Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam Figure 4. R for three productivity metrics classified by land cover type identified in PhenoCam data data and overall for all data points. The productivity metrics are integrated vegetation index (IVI), and overall for all data points. The productivity metrics are integrated vegetation index (IVI), peak peak vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands vegetation index (PVI) and amplitude of vegetation index (AVI). The classifications are grasslands (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). (GR), shrublands (SH), deciduous/broadleaf (DB), evergreen (EN) and overall (OV). 3.2. Historical Trend Analysis for Western United States Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Of the 10 datasets evaluated, five extend back to 1982 (CMGLSP, VIPPHEN_EVI2, VIPPHEN_NDVI, AVHRRP and NPP) and therefore provide an appropriate historical record for relatively long-term trend analysis. Based on the comparison with PhenoCam, we chose to analyze CMGLSP and VIPPHEN_EVI2, plus NPN for SOS (NPN results in Appendix C Figures A2–A4). CMGLSP most often had the ideal combination of high R and low mean bias across land cover types and VIPPHEN_EVI2 also agreed well with PhenoCam, especially in grasslands and shrublands. Long-term trends indicate large changes in key phenology dates in some places by more than 70 days (2 days per year for 35 years) between 1982 and 2016 (Figure 5; see Figures A1 and A5–A7 for mean and standard deviation maps and histograms of changes). Trend patterns had high spatial heterogeneity, with substantial interspersion of pixels with delayed and advanced phenology in both datasets and all metrics. The strength of the trends varied spatially, but in both datasets the areas of greatest change in date occurred in New Mexico and Arizona (later SOS), low-lying areas of California (earlier EOS), and the Great Basin (earlier SOS, mixed strong trends for EOS). EOS in the Southwest was also generally later but was more variable in California. Higher elevation areas across Montana, Washington, Idaho and Colorado showed trends toward later SOS but mixed results for EOS. Over time, variation in phenology dates slightly increased across most regions, with very large increases in variability corresponding to areas with larger changes in date, namely the Southwest for Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 33 3.2. Historical trend analysis for western United States Of the 10 datasets evaluated, five extend back to 1982 (CMGLSP, VIPPHEN_EVI2, VIPPHEN_NDVI, AVHRRP and NPP) and therefore provide an appropriate historical record for Remote Sens. 2020, 12, 2538 9 of 27 relatively long-term trend analysis. Based on the comparison with PhenoCam, we chose to analyze CMGLSP and VIPPHEN_EVI2, plus NPN for SOS (NPN results in Appendix C Figures A2–A4). SOS and POS, and scattered areas around the edges of the Great Basin as well as California for EOS. CMGLSP most often had the ideal combination of high R and low mean bias across land cover types The spatial patterns of agreement on the direction of trend across CMGLSP and VIPPHEN-EVI2 were and VIPPHEN_EVI2 also agreed well with PhenoCam, especially in grasslands and shrublands. also relatively heterogenous, with the least consistency occurring in Wyoming for SOS and in the Great Long-term trends indicate large changes in key phenology dates in some places by more than 70 Basin for EOS. days (2 days per year for 35 years) between 1982 and 2016 (Figure 5; see Figure A1 and Figure A5– The consistency results over the entire study region showed ~60% agreement in trend direction for A7 for mean and standard deviation maps and histograms of changes). Trend patterns had high all three metrics (Figure 6). For SOS and POS ~30% of pixels agreed on earlier and later dates in trends, spatial heterogeneity, with substantial interspersion of pixels with delayed and advanced phenology whereas for EOS, 44.4% of pixels in both datasets agreed on a later EOS (16.2% agreed on an earlier in both datasets and all metrics. The strength of the trends varied spatially, but in both datasets the EOS). In assessing consistency by ecoregion, we found agreement between the datasets in the range of areas of greatest change in date occurred in New Mexico and Arizona (later SOS), low-lying areas of 50–70% of pixels for most of the regions and metrics as well as agreement on the predominant direction California (earlier EOS), and the Great Basin (earlier SOS, mixed strong trends for EOS). EOS in the of change within ecoregions (Figure 6). For SOS, the two datasets agreed that more area trended Southwest was also generally later but was more variable in California. Higher elevation areas across towards earlier dates in Marine West Coast Forests, Mediterranean California, North American Deserts, Montana, Washington, Idaho and Colorado showed trends toward later SOS but mixed results for and Northwestern Forested Mountains, and later dates in Great Plains, Southern Semiarid Highlands, EOS. Over time, variation in phenology dates slightly increased across most regions, with very large and Temperate Sierras. For POS, the datasets agreed on more area trending toward earlier dates in increases in variability corresponding to areas with larger changes in date, namely the Southwest for Marine West Coast Forests, Mediterranean California, and North American Deserts, and later dates in SOS and POS, and scattered areas around the edges of the Great Basin as well as California for EOS. Great Plains, Northwestern Forested Mountains, Southern Semiarid Highlands, and Temperate Sierras. The spatial patterns of agreement on the direction of trend across CMGLSP and VIPPHEN-EVI2 were We found the most dataset agreement across ecoregions for EOS, in which all regions agreed on more also relatively heterogenous, with the least consistency occurring in Wyoming for SOS and in the area trending towards later dates except Mediterranean California. Great Basin for EOS. Figure 5. Change in date, consistency, and change in variation of CMGLSP and VIPPHEN_EVI2 Figure 5. Change in date, consistency, and change in variation of CMGLSP and VIPPHEN_EVI2 across three phenology metrics from 1982 to 2016. Negative values correspond to earlier dates (or fewer days) across three phenology metrics from 1982 to 2016. Negative values correspond to earlier dates (or and fewer day positive s) and pos values toitiv later e va dates lues to (or lat mor er date e days). s (orConsistency more days). Consisten shows agreement cy shows agreement in th in the direction ofe change between datasets. Change in date is the regressed trend slope multiplied by the number of direction of change between datasets. Change in date is the regressed trend slope multiplied by the years. Change in variation is the regressed trend slope, multiplied by the number of years, of the number of years. Change in variation is the regressed trend slope, multiplied by the number of years, absolute of the ab value solute of value of the residuals the residuals of the of the trend estimate trend est against imate a time. gainst time. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 33 The consistency results over the entire study region showed ~60% agreement in trend direction for all three metrics (Figure 6). For SOS and POS ~30% of pixels agreed on earlier and later dates in trends, whereas for EOS, 44.4% of pixels in both datasets agreed on a later EOS (16.2% agreed on an earlier EOS). In assessing consistency by ecoregion, we found agreement between the datasets in the range of 50–70% of pixels for most of the regions and metrics as well as agreement on the predominant direction of change within ecoregions (Figure 6). For SOS, the two datasets agreed that more area trended towards earlier dates in Marine West Coast Forests, Mediterranean California, North American Deserts, and Northwestern Forested Mountains, and later dates in Great Plains, Southern Semiarid Highlands, and Temperate Sierras. For POS, the datasets agreed on more area trending toward earlier dates in Marine West Coast Forests, Mediterranean California, and North American Deserts, and later dates in Great Plains, Northwestern Forested Mountains, Southern Semiarid Highlands, and Temperate Sierras. We found the most dataset agreement across ecoregions Remote Sens. for EO 2020 S, ,in w 12, 2538 hich all regions agreed on more area trending towards later dates except Mediterranea 10 n of 27 California. Figure 6. Consistency of CMGLSP and VIPPHEN_EVI2 long-term (1982–2016) trend direction over the Figure 6. Consistency of CMGLSP and VIPPHEN_EVI2 long-term (1982–2016) trend direction over entire study area and classified by ecoregion. the entire study area and classified by ecoregion. 4. Discussion 4. Discussion In our comparison of 10 leading LSP datasets against a network of PhenoCam near-surface cameras, In our comparison of 10 leading LSP datasets against a network of PhenoCam near-surface we found promising results in terms of R and mean bias in some landcover types, but substantial cameras, we found promising results in terms of R and mean bias in some landcover types, but substantial variation between datasets and metrics (similar to [32] for SOS estimates). This highlights variation between datasets and metrics (similar to [32] for SOS estimates). This highlights a need for a need for improved communication and more extensive dialogue on the strengths, weaknesses and improved communication and more extensive dialogue on the strengths, weaknesses and development development goals of LSP datasets, especially with the proliferation of applied data users from a goals of LSP datasets, especially with the proliferation of applied data users from a variety of research variety of research and management practices. and management practices. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing When choosing an appropriate dataset for analysis, applied users should consider study location, years needed, spatial resolution, processing capacity, and quality in di erent land cover classes. For phenology metrics, our results indicate MCD12Q2 and CMGLSP best matched PhenoCam observations derived using a 10% amplitude threshold and the 90th GCC percentile (Figure 3). For those interested in phenology prior to the 2000s, CMGLSP has global coverage and extends back to 1982, but has a coarse 0.05 degree (~5 km) spatial resolution and can only be acquired by request. VIPPHEN_EVI2 has the same temporal coverage and spatial resolution, can be downloaded directly, and agreed almost as well. For users working with high resolution data (such as from GPS collars in the wildlife field) and conducting analyses after 2001, MODIS-based datasets at 250 to 500 m spatial resolution provide finer-scale observations. MCD12Q2, with 500 m spatial resolution and global coverage, agreed best with PhenoCam and was the only dataset that had high correlations in evergreen forests. eMODIS also agreed well and has 250 m spatial resolution but is only available in the CONUS. The DLC method can be applied globally (we used MOD09Q1 as the base input) but requires substantial processing which may be prohibitive depending on study area size. This method also provides daily values of NDVI and instantaneous rate of green-up (IRG), which are useful for certain movement ecology and habitat questions such as those relating to the green-wave hypothesis [21,23–25,57]. Remote Sens. 2020, 12, 2538 11 of 27 Readily available datasets from MODIS (e.g., MOD13A1, MOD13Q1, and eMODIS NDVI) only provide cleaned and pre-processed weekly to bi-weekly NDVI/EVI values and therefore could provide a coarser resolution IRG time series. Based on the di erences in dataset correlation and bias we found between land cover types, users should consider the predominant land cover of their study area in choosing a phenology or productivity dataset (Figure 3; Figure 4). For instance, eMODIS agreed best for PIRGd in shrublands, yet MCD12Q2 was better in grasslands and evergreen forests. Large di erences in correlation were observed between productivity datasets, which correlated poorly overall but well in certain landcover types. For IVI, DLC agreed well in grasslands, but MODIS_NPP agreed best in deciduous/broadleaf forests. The lowest correlations between LSP estimates and PhenoCam observations were in evergreen forests, where identifying phenological metrics is challenging because of the small annual change in greenness amplitude and over-saturation [35,58]. Using the red chromatic coordinate to process PhenoCam transition dates, as opposed to GCC, has shown potential to improve predictions of SOS and EOS in these environments [59]. Of all the LSP datasets, MCD12Q2 agreed best with PhenoCam at evergreen sites, possibly because it is EVI based [40]. Ecologists studying wildlife in areas with multiple land cover types should be aware that phenology metrics in some land cover types are more reliable than others, which could influence results for analyses comparing wildlife GPS locations across di erent land cover classes [17,60,61]. Quality of LSP metrics may di er geographically, and therefore our results are not necessarily indicative of dataset quality in other parts of North America or around the world [62]. Quantifying phenology and productivity trends through time is crucial to improve our understanding of ecosystem processes, such as how changing forage in critical habitat or management units a ects wildlife and how to manage for future climate and phenology scenarios [8,26,63,64]. We found a 44.4% agreement between CMGLSP and VIPPHEN_EVI2 trends toward a later EOS across the Western U.S and only a 16.2% agreement toward an earlier EOS, which is consistent with other studies signaling a general pattern toward a later EOS [9,65,66]. Trends are influenced by a variety of factors, including temperature and precipitation, species succession, human and natural disturbance, and land management practices (as observed by [37] from 1982 to 2006). The complex interactions between these factors make it dicult to interpret the ecological significance of large and spatially heterogeneous trends, especially across diverse ecosystems. The large phenology trends and high variability we observed over a 35 year period likely reflect real changes in remote sensing metrics, but highlight the need to understand what exactly is changing throughout time, the role of data processing choices, and whether the changes are biologically significant. We know that patterns of temperature and precipitation in some areas are complex, leading to variability in the timing, seasonality, and spatial heterogeneity of snow coverage [67], soil moisture, drought, and storms that may all influence the timing, variability, and heterogeneity of vegetation phenology and productivity. Variability can be an important component to ecosystem processes and to understanding the degree and mechanism of trends. For instance, vegetation strongly a ected by climate conditions, such as invasive grasses responding to rainfall, can cause high year-to-year variability and therefore introduce uncertainty into the direction of overall trend, as well as have large e ects on the ecosystem [68,69]. In our results, some regions in which LSP metrics changed by more than 70 days often showed high and increasing variability. Though these trends may be showing real changes in phenology throughout time within a pixel, it is challenging to interpret whether these changes are driven by variability caused by weather, invasive species, or other factors, or actually represent the changing dynamics of vegetation of interest, such as important forage species. This overall assessment is a first step in understanding the overall patterns and degree of change across the CONUS. Several challenges exist in producing accurate LSP estimates and assessing their quality. Large di erences in LSP metrics and trend are related to sensor-specific characteristics, such as revisit time and spectral resolution, as well as processing algorithms used to extract metrics, which may include Remote Sens. 2020, 12, 2538 12 of 27 cloud and atmospheric correction, gap filling and curve fitting. Greenness (NDVI/EVI/GCC) values di er between satellite platforms, as well as with near-surface cameras, and are therefore not always comparable without applying translation equations [70]. Future work may scale PhenoCam and satellite derived productivity metrics to enable a comparison of mean bias and other statistics. Di erences in consistency (by-pixel agreement on trend direction) between CMGLSP and VIPPHEN_EVI2 result solely from di erent processing algorithms, as both CMGLSP and VIPPHEN_EVI2 are based on EVI2 calculated from AVHRR and MODIS at a 0.05 degree spatial resolution (for CMGLSP, see [46]; for VIPPHEN_EVI2, see [47]). These two datasets yielded conflicting yet large trends in phenology around the Great Basin and Eastern Washington and Oregon, specifically for EOS. These areas experienced large fires, landcover changes, and other disturbances. It is possible that the di erent approaches to fitting NDVI/EVI curves and extracting metrics for these datasets are a ected di erently by these landcover type changes and thus lead to these conflicting signals. In general, LSP datasets use di erent processing methodologies and di erent threshold values [36,37,71], and dataset development goals rarely correspond to specific biological events important to data users. For instance, biological events recorded by a botanist to mark the start of growing season may include first leaf or flower dates, whereas an LSP dataset such as MCD12Q2 captures SOS as the date when EVI crosses the 15% threshold of AVI [45]. Ref. [32] used in situ first leaf observations to validate SOS from various LSP datasets and found high root mean squared error (~20 days), indicating a need to better understand the link between biological events and image reflectance values. Our use of near-surface cameras removes variability that may occur when comparing LSP metrics with the timing of biological processes from in situ observations. Although some research assesses LSP quality through comparison with other kinds of ground-based datasets [72,73], this may confound di erences resulting from discrepancies between biological events and their reflectance values with di erences in image quality, which can vary inter-annually and regionally [62]. Datasets also vary in how they deal with growing seasons that span calendar years or have multiple annual cycles. Most datasets calculate phenology metrics based on a continuous time series yet report metrics in discrete annual data layers. The one exception in our analysis is the DLC dataset, which fits a curve to a single year of data. When phenology cycles span calendar years, such as when SOS occurs in the fall and/or when EOS occurs in late winter/early spring, users should ensure they understand the approach used to store the metrics within an annual data layer. Datasets handle patterns of multiple or missing annual wetness cycles in di erent ways, with some datasets reporting up to three annual cycles while others only report phenology dates from a single cycle that is not indexed to a general time of year. We most commonly see multiple annual cycles of green-up in arid and monsoonal environments, such as in the Southwest [74] and Mediterranean California, but this varies by year.. When the number of cycles varies, such as when drought years lead to missed wet periods [75], or seasonal weather patterns dictate two to three annual cycles, analyses of trends by year, as we have conducted, could suggest large changes in dates that might not reflect the complexity of the ecological impacts well. We briefly evaluated these potential issues, and did find strong positive trends in SOS and POS in areas prone to multiple annual cycles, but overall found that both varying numbers of cycles and small amplitude cycles had more limited extent than could easily explain the large changes in trend we found. Inherently, the scale of near-surface PhenoCam observations does not match those of remote sensing pixels, which include heterogenous vegetation and even land cover class [36,38,76]. The mismatch can be substantial even for our comparisons of PhenoCam with MODIS-based 250 or 500 m spatial resolution products versus CMGLSP at ~5 km resolution, in which the synthesis of large pixels could include billions of plants. LSP observations e ectively provide only a broad picture of landscape phenology, and the scale of imagery has important implications for users. For example, CMGLSP at ~5 km spatial resolution may not provide the spatial variability needed to assess the behavior of a white-tailed deer or other animals with home ranges within one or two pixels, whereas it may be useful for understanding movement patterns of long distance migrating animals such as eagles or Remote Sens. 2020, 12, 2538 13 of 27 mule deer. Given the heterogeneity of vegetation within large pixels, we were surprised at the high agreement of CMGLSP with PhenoCam, at a spatial resolution of 0.05 degrees. The dynamics of LSP metrics in large pixels are complex and do not necessarily represent the average of smaller pixels [77]. For SOS, [78] found that the heterogeneity of landcover and SOS, as well as vegetation growth speed during spring, all influenced the scale e ect. Variation in spatial footprints of PhenoCam data and resolution of LSP datasets could also add noise to the comparison of LSP and near-surface datasets. New fine-resolution datasets, such as daily 30 m EVI products developed using fusion between MODIS and Landsat [79,80] or combinations of Landsat and Sentinal-2 imagery to increase temporal resolution and enable 30 m phenology retrievals [81], could better match the scale of reference data and user analysis objectives. Applied data users will benefit from future developments in processing capacity and dataset construction, a better understanding of phenology drivers, and greater understanding of how algorithms intersect with phenology predictions. The growing popularity of platforms able to process large data quickly (such as Google Earth Engine) may render daily, global datasets derived from the DLC method readily available in the future. In terms of improving LSP estimates, the use of mechanistic models to predict key metrics [82–84] may help to address data quality issues and discrepancies across land cover types and ecoregions. These models couple remote sensing data with local observations or other models of elevation, temperature, precipitation, and plant phenology to improve phenology and productivity metrics. While still limited, in the future these models may be developed into a single framework to derive more relevant phenology and productivity estimates across diverse regions using localized data. Similarly, an increased understanding of the drivers of phenological changes and variability on a local or regional scale may help users to better interpret and plan for long-term patterns of change. For example, [85] showed that temperature, growing-season precipitation, and snowpack had larger e ects than most management and habitat treatments on annual phenology metrics in sagebrush communities in Southwestern Wyoming, but identified that treatments had some impact on IVI and late season phenology metrics. Finally, deeper examination of how processing and curve-fitting techniques influence the resulting NDVI/EVI time series and phenology metrics will improve understanding of how to interpret large di erences in phenology that are predicted from LSP data. 5. Conclusions Land surface phenology (LSP) data provide the means to assess large spatial and temporal patterns in environmental change and understand a variety of ecological questions, such as those related to surface energy balance and animal fitness, movement, and habitat use. The quality of these data depend on factors including cloud cover, atmospheric e ects, processing algorithm, land cover type, spatial resolution and regionality. We found highly variable agreement between 10 leading LSP datasets and PhenoCam and outlined dataset choices based on spatial and temporal coverage, processing capability, and agreement across ecosystems, providing a basis for informed and justifiable decisions about data choices. Most datasets we examined and many of the highest-quality datasets are readily available and do not require users to extract metrics, as is often performed by individual users using raw NDVI or EVI values. Given the popularity and importance of phenology and productivity data in a wide range of research and management fields [61,86], a better understanding and justification for the use of certain datasets can only help to bolster the e ectiveness and validity of results. Author Contributions: Conceptualization, T.A.G., N.L.M., J.A.M., A.N.J. and G.W.C.; methodology, E.E.B., T.A.G. and N.L.M.; software, E.E.B. and N.L.M.; validation, E.E.B.; formal analysis, E.E.B. and N.L.M.; resources, E.E.B. and T.A.G.; writing—original draft preparation, E.E.B. and N.L.M.; writing—review and editing, T.A.G., N.L.M., J.A.M., A.N.J. and G.W.C.; visualization, E.E.B.; supervision, T.A.G.; project administration, T.A.G.; funding acquisition, T.A.G.; E.E.B., T.A.G., and N.L.M. contributed equally. All authors have read and agreed to the published version of the manuscript. Remote Sens. 2020, 12, 2538 14 of 27 Funding: This work was supported by the U.S. Geological Survey Northern Rocky Mountain Science Center and the North Central Climate Adaptation Science Center. The Wyoming Landscape Conservation Initiative provides support to Chong and Johnston. The work of Ethan Berman was done under contract to the U.S. Geological Survey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Acknowledgments: We would like to thank the Wyoming Game and Fish Department, Bureau of Land Management Kemmerer, Montana Department of Fish, Wildlife and Parks, and the National Park Service for their support in the development of this research. In addition, we appreciate the help of Jill Randall, Troy Fieseler, Brent Jamison, Arvid Aase, Kelly Prott, Justin Gude, as well as David Wood and anonymous reviewers for their feedback and comments. Conflicts of Interest: The authors declare no conflict of interest. The external funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Appendix A Detailed table of phenology and productivity metrics along with formulas for metric calculations where not readily available. Remote Sens. 2020, 12, 2538 15 of 27 Table A1. Detailed table of phenology and productivity metrics. NPN First Leaf Spring Dataset AVHRRP eMODIS MCD12Q2 CMGLSP VIPPHEN-EVI2 VIPPHEN-NDVI Landsat NPP MODIS NPP DLC PhenoCam Index Temperature, weather MODIS NDVI MODIS NDVI Calibrated events, in situ composites, composites, MOD09Q1 Calibrated radiance data Time series of AVHRR and AVHRR and MODIS AVHRR and MODIS RGB photo time Input data radiance data observations of one lilac meteorological meteorological surface reflection (level 1-B). Aqua only NBAR-EVI2 MODIS EVI2 EVI2 NDVI series (level 1-B). and two honeysuckle inputs, biome inputs, biome 8 day composites species specific contraints specific contraints Data are smoothed Data are smoothed 2 band EVI Cleaned, NDVI with 2 step filter and 3 with 2 step filter and 3 Cleaned, NDVI (EVI2) calculated, and Eliminating year moving window year moving window Calculated at PRISM Daily GPP and Daily GPP and calculated, and determination combined into 7 outliers and approach, AVHRR approach, AVHRR minimum and maximum maintenance maintenance Cleaned, NDVI Raw data processing details combined into 7 day from raw data. GCC day composites filling dormant and MODIS continuity and MODIS continuity temperature values reach respiration values respiration values calculated composites using max Smoothed with using max NDVI period values is based on linear is based on linear a “stable” state. calculated calculated NDVI and quality Savitsky-Golay and quality regressions between regressions between filter overlapping data. overlapping data. SOS NDVI value NPP extracted from NPP extracted from PELT method, Cubic spline fit, phenological Average of the first day 2 piece logistic compared with average Phenological changes Phenological changes daily GPP, daily GPP, AIC-selected SOS is 15% of changes are each year that the model. SOS/EOS Did not receive of previous 36 days and are identified using are identified using maintenance maintenance smoothing spline, Metric extraction details AVI, EOS is 15% identified using properties of the three are local information SOS is when a trend half-max method (at half-max method (at respiration, and respiration, and SOS at 10% of AVI, of AVI hybrid piecewise individual species models min/max from shift is identified. only 35%) only 35%) annual growth annual growth EOS at 10% of green-down logistic model are met (SOS only) 2nd derivatives Opposite for EOS respiration respiration green-down curve Spatial resolution 1 km 250 m 500 m 0.05 0.05 0.05 0.0417 30 m 250 m 250 m Point Temporal coverage 1989–2014 2001–2018 2001–2017 1982–2016 1981–2016 1981–2016 1981–2019 1986–2018 2001–2018 2000–present 2000–present Accessible online? yes yes yes no yes yes yes yes yes yes yes Significant post-processing no no no no no no no no no yes no time requirement? Under production? no yes yes yes no no yes yes yes yes yes Spatial coverage CONUS CONUS Global Global Global Global CONUS CONUS CONUS Global U.S. and Canada Years acquired 2002–2014 2002–2014 2002–2014 1982–2016 1982–2016 1982–2016 1982–2016 2002–2014 2002–2014 2002–2014 2002–2014 Reference [43] [44] [45] [46] [47] [47] [41] [42] [42] [23] [48] http: https: https://www.ntsg. https://www.usgs.gov/ xiaoyang. //files.ntsg.umt.edu/ //www.usgs.gov/ https://earthdata. https: https: https://www.usanpn.org/ umt.edu/project/ https://earthdata. Retrieval location land-resources/eros/ zhang@sdstate. data/NTSG_ phenocamr land-resources/ nasa.gov/ //earthdata.nasa.gov/ //earthdata.nasa.gov/ data/spring_indices landsat/landsat- nasa.gov/ phenology edu Products/MOD17/ eros/phenology productivity.php MODIS_250/ Metrics available or calculated by authors? (A is available, C is calculated, N is not available/used) SOS A A A A A A A N N C C PIRGd C C C C C C N N N C C POS A A A A A A N N N C C EOS A A A A A A N N N C C LOSp C C C C C C N N N C C LGS A A C C A A N N N C C IVI A A A A A A N A A C C PVI A A C A A A N N N C C AVI A A A C A A N N N C C 1 2 Collection 6 of eMODIS starts in 2003, but Collection 5 is available from 2001. This analysis used Collection 5 for 2002; Some sites in Panama, Hawaii, and Europe. Remote Sens. 2020, 12, 2538 16 of 27 Metrics reflecting di erent components of forage phenology and productivity were derived (or readily provided) from the 10 datasets used in this analysis. In total we assessed ten di erent metrics, six related to phenology and three related to productivity. Many of the datasets provide the desired metrics, and when available these values were used. In other cases, the metrics were derived using the following equations (with the exception of the Bischof dataset, for which methods can be found in [23]: POS SOS PIRGd = SOS + (assuming logistic growth) (A1) EOS POS = MAX Greenness (A2) i=SOS LOSp = POS SOS (A3) LGS = EOS SOS (A4) EOS IVI = Greenness (A5) i=SOS PVI = AVI + SMV (A6) AVI = PVI SMV (A7) where PIRGd is peak instantaneous rate of spring greenup date, SOS is start of spring date, POS is peak of season, MAX is the maximum value for the given period, EOS is end of season, greenness is the daily value of greenness (whether EVI, NDVI, or GCC), LOSp is length of spring, LGS is length of growing season, IVI is integrated vegetation index, PVI is peak vegetation index, AVI is amplitude of vegetation index and SMV is spring minimum value. All dates were represented as Julian date for calculation Purp. Appendix B R and mean bias results of 10 datasets tested against PhenoCam observation. The number of observations is displayed in parenthesis. Table A2. Overall results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.07 (98) 0.20 (98) 0.33 (106) 0.22 (79) 0.02 (98) 0.00 (72) 0.00 (104) 0.13 (106) 0.29 (106) CMGLSP 0.37 (102) 0.52 (101) 0.39 (105) 0.20 (104) 0.03 (101) 0.06 (100) 0.03 (100) 0.00 (105) 0.00 (104) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.10 (88) NA (NA) NA (NA) MCD12Q2 0.35 (76) 0.55 (76) 0.54 (76) 0.35 (76) 0.02 (76) 0.02 (76) 0.01 (76) 0.01 (76) 0.02 (76) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (78) NA (NA) NA (NA) NPN 0.03 (104) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.26 (106) 0.42 (106) 0.32 (106) 0.16 (106) 0.01 (106) 0.04 (106) 0.15 (106) 0.13 (106) 0.34 (103) VIPPHEN_NDVI 0.36 (106) 0.44 (106) 0.37 (106) 0.23 (106) 0.09 (106) 0.10 (106) 0.16 (106) 0.15 (106) 0.26 (97) DLC 0.35 (106) 0.50 (106) 0.53 (106) 0.45 (106) 0.11 (106) 0.09 (106) 0.00 (106) 0.05 (106) 0.05 (106) eMODIS 0.26 (99) 0.33 (99) 0.25 (103) 0.22 (101) 0.05 (99) 0.05 (97) 0.02 (103) 0.11 (106) 0.19 (103) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 9.88 (98) 15.92 (98) 20.54 (106) 37.92 (79) 14.09 (98) 43.44 (72) CMGLSP 8.41 (102) 0.38 (101) 12.10 (105) 0.52 (104) 17.75 (101) 3.78 (100) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 4.39 (76) 0.83 (76) 7.05 (76) 19.92 (76) 12.45 (76) 24.32 (76) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 19.10 (104) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 4.39 (106) 2.69 (106) 10.76 (106) 23.58 (106) 16.15 (106) 27.97 (106) VIPPHEN_NDVI 25.49 (106) 4.16 (106) 18.17 (106) 32.98 (106) 44.66 (106) 58.47 (106) DLC 5.63 (106) 17.66 (106) 21.32 (106) 4.13 (106) 27.95 (106) 9.76 (106) eMODIS 4.61 (99) 15.74 (99) 23.82 (103) 23.26 (101) 24.26 (99) 25.15 (97) Remote Sens. 2020, 12, 2538 17 of 27 Table A3. Grasslands (GR) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.10 (41) 0.13 (41) 0.32 (45) 0.52 (38) 0.03 (41) 0.15 (34) 0.30 (45) 0.05 (45) 0.32 (45) CMGLSP 0.38 (44) 0.52 (44) 0.55 (45) 0.56 (44) 0.00 (44) 0.01 (43) 0.00 (43) 0.00 (45) 0.05 (44) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.15 (37) NA (NA) NA (NA) MCD12Q2 0.45 (37) 0.69 (37) 0.78 (37) 0.51 (37) 0.00 (37) 0.00 (37) 0.00 (37) 0.05 (37) 0.40 (37) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (34) NA (NA) NA (NA) NPN 0.07 (45) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.41 (45) 0.64 (45) 0.69 (45) 0.50 (45) 0.04 (45) 0.24 (45) 0.03 (45) 0.12 (45) 0.48 (45) VIPPHEN_NDVI 0.35 (45) 0.50 (45) 0.53 (45) 0.68 (45) 0.00 (45) 0.18 (45) 0.00 (45) 0.05 (45) 0.43 (45) DLC 0.57 (45) 0.69 (45) 0.58 (45) 0.72 (45) 0.07 (45) 0.28 (45) 0.43 (45) 0.04 (45) 0.02 (45) eMODIS 0.48 (42) 0.60 (42) 0.56 (44) 0.60 (44) 0.08 (42) 0.18 (42) 0.07 (44) 0.14 (45) 0.13 (44) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 7.17 (41) 4.32 (41) 15.56 (45) 62.47 (38) 24.98 (41) 74.35 (34) CMGLSP 9.68 (44) 3.02 (44) 3.22 (45) 18.93 (44) 11.32 (44) 11.40 (43) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 8.57 (37) 1.11 (37) 7.35 (37) 40.32 (37) 16.92 (37) 48.89 (37) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 20.18 (45) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 15.07 (45) 0.31 (45) 16.69 (45) 46.44 (45) 32.76 (45) 61.51 (45) VIPPHEN_NDVI 33.00 (45) 8.67 (45) 16.67 (45) 54.44 (45) 50.67 (45) 87.44 (45) DLC 16.04 (45) 24.19 (45) 19.71 (45) 25.04 (45) 36.76 (45) 41.09 (45) eMODIS 5.43 (42) 8.80 (42) 20.64 (44) 39.02 (44) 30.45 (42) 49.74 (42) Table A4. Shrublands (SH) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.01 (7) 0.07 (7) 0.00 (7) 0.20 (4) 0.03 (7) 0.40 (4) 0.02 (7) 0.16 (7) 0.08 (7) CMGLSP 0.22 (7) 0.28 (7) 0.28 (7) 0.62 (7) 0.12 (7) 0.20 (7) 0.01 (7) 0.19 (7) 0.10 (7) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.00 (7) NA (NA) NA (NA) MCD12Q2 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) 0.00 (1) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.01 (7) NA (NA) NA (NA) NPN 0.13 (7) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.35 (7) 0.28 (7) 0.50 (7) 0.39 (7) 0.59 (7) 0.07 (7) 0.09 (7) 0.46 (7) 0.14 (7) VIPPHEN_NDVI 0.31 (7) 0.22 (7) 0.00 (7) 0.15 (7) 0.17 (7) 0.35 (7) 0.00 (7) 0.19 (7) 0.12 (7) DLC 0.55 (7) 0.53 (7) 0.43 (7) 0.54 (7) 0.44 (7) 0.26 (7) 0.00 (7) 0.14 (7) 0.07 (7) eMODIS 0.68 (7) 0.89 (7) 0.85 (7) 0.03 (6) 0.52 (7) 0.22 (6) 0.05 (7) 0.22 (7) 0.11 (7) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 38.43 (7) 30.64 (7) 23.86 (7) 24.50 (4) 13.57 (7) 34.00 (4) CMGLSP 25.43 (7) 6.79 (7) 10.86 (7) 8.71 (7) 35.29 (7) 16.71 (7) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 11.00 (1) 11.00 (1) 10.00 (1) 2.00 (1) 2.00 (1) 9.00 (1) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 108.57 (7) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 25.71 (7) 20.14 (7) 15.57 (7) 7.86 (7) 9.14 (7) 17.86 (7) VIPPHEN_NDVI 22.00 (7) 11.79 (7) 2.57 (7) 13.71 (7) 18.43 (7) 35.71 (7) DLC 7.29 (7) 10.79 (7) 6.14 (7) 24.00 (7) 0.14 (7) 31.29 (7) eMODIS 23.00 (7) 21.43 (7) 20.86 (7) 15.00 (6) 1.14 (7) 6.67 (6) Remote Sens. 2020, 12, 2538 18 of 27 Table A5. Deciduous/Broadleaf (DB) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.37 (23) 0.50 (23) 0.21 (27) 0.10 (26) 0.02 (23) 0.05 (23) 0.26 (25) 0.10 (27) 0.34 (27) CMGLSP 0.65 (24) 0.37 (23) 0.00 (26) 0.01 (26) 0.00 (23) 0.00 (23) 0.00 (23) 0.23 (26) 0.26 (26) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.24 (20) NA (NA) NA (NA) MCD12Q2 0.28 (26) 0.31 (26) 0.19 (26) 0.52 (26) 0.07 (26) 0.29 (26) 0.00 (26) 0.13 (26) 0.25 (26) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.90 (10) NA (NA) NA (NA) NPN 0.90 (27) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.04 (27) 0.45 (27) 0.13 (27) 0.21 (27) 0.17 (27) 0.01 (27) 0.54 (27) 0.00 (27) 0.33 (27) VIPPHEN_NDVI 0.59 (27) 0.33 (27) 0.07 (27) 0.21 (27) 0.07 (27) 0.13 (27) 0.59 (27) 0.16 (27) 0.25 (23) DLC 0.28 (27) 0.51 (27) 0.77 (27) 0.33 (27) 0.45 (27) 0.07 (27) 0.22 (27) 0.20 (27) 0.32 (27) eMODIS 0.02 (25) 0.36 (25) 0.26 (26) 0.59 (26) 0.02 (25) 0.03 (25) 0.24 (26) 0.05 (27) 0.59 (26) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 13.61 (23) 5.28 (23) 0.30 (27) 17.15 (26) 18.65 (23) 39.83 (23) CMGLSP 19.79 (24) 18.33 (23) 28.92 (26) 37.04 (26) 2.13 (23) 6.30 (23) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 15.58 (26) 7.46 (26) 1.65 (26) 4.08 (26) 18.23 (26) 19.65 (26) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 18.74 (27) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 17.56 (27) 8.39 (27) 1.78 (27) 1.15 (27) 20.33 (27) 16.41 (27) VIPPHEN_NDVI 49.63 (27) 27.22 (27) 3.81 (27) 10.63 (27) 46.81 (27) 60.26 (27) DLC 20.19 (27) 27.54 (27) 18.52 (27) 11.85 (27) 39.70 (27) 8.33 (27) eMODIS 8.08 (25) 0.08 (25) 8.92 (26) 11.15 (26) 18.00 (25) 21.48 (25) Table A6. Evergreen (EN) results. R Squared (N) Phenology Metrics Productivity Metrics Dataset SOS PIRGd POS EOS LOSp LGS IVI PVI AVI AVHRRP 0.03 (25) 0.03 (25) 0.04 (25) 0.45 (9) 0.02 (25) 0.49 (9) 0.01 (25) 0.01 (25) 0.04 (25) CMGLSP 0.33 (25) 0.40 (25) 0.28 (25) 0.11 (25) 0.11 (25) 0.08 (25) 0.01 (25) 0.13 (25) 0.01 (25) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.04 (22) NA (NA) NA (NA) MCD12Q2 0.44 (11) 0.54 (11) 0.64 (11) 0.31 (11) 0.08 (11) 0.07 (11) 0.15 (11) 0.91 (11) 0.75 (11) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) 0.02 (25) NA (NA) NA (NA) NPN 0.00 (25) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 0.16 (25) 0.06 (25) 0.01 (25) 0.00 (25) 0.04 (25) 0.06 (25) 0.06 (25) 0.01 (25) 0.09 (22) VIPPHEN_NDVI 0.24 (25) 0.19 (25) 0.04 (25) 0.17 (25) 0.12 (25) 0.09 (25) 0.06 (25) 0.04 (25) 0.52 (20) DLC 0.10 (25) 0.07 (25) 0.00 (25) 0.09 (25) 0.00 (25) 0.02 (25) 0.04 (25) 0.00 (25) 0.16 (25) eMODIS 0.09 (23) 0.04 (23) 0.00 (24) 0.09 (23) 0.00 (23) 0.00 (22) 0.07 (24) 0.00 (25) 0.11 (24) Mean Bias (N) Phenology Metrics Dataset SOS PIRGd POS EOS LOSp LGS AVHRRP 55.56 (25) 52.88 (25) 51.20 (25) 14.33 (9) 3.36 (25) 23.56 (9) CMGLSP 29.44 (25) 9.52 (25) 9.40 (25) 5.04 (25) 37.84 (25) 24.40 (25) Landsat_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) MCD12Q2 32.64 (11) 22.68 (11) 13.73 (11) 11.36 (11) 17.91 (11) 44.00 (11) MODIS_NPP NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) NPN 7.52 (25) NA (NA) NA (NA) NA (NA) NA (NA) NA (NA) VIPPHEN_EVI2 24.52 (25) 15.70 (25) 7.88 (25) 18.84 (25) 15.64 (25) 5.68 (25) VIPPHEN_NDVI 3.96 (25) 26.14 (25) 49.32 (25) 39.36 (25) 46.36 (25) 35.40 (25) DLC 25.64 (25) 2.06 (25) 28.68 (25) 5.60 (25) 4.04 (25) 31.24 (25) eMODIS 32.35 (23) 44.93 (23) 47.17 (24) 10.17 (23) 27.17 (23) 8.32 (22) Remote Sens. 2020, 12, 2538 19 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 25 of 33 Appendix C Appendix C Additional figures relating to the long-term trend analysis of CMGLSP, VIP_EVI, and NPN. Additional figures relating to the long-term trend analysis of CMGLSP, VIP_EVI, and NPN. Figure A1. Mean and standard deviation of historical (1982–2016) data from CMGLSP and Figure A1. Mean and standard deviation of historical (1982–2016) data from CMGLSP and VIPPHEN_EVI2. VIPPHEN_EVI2. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 20 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 33 Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 33 Figure A2 Figure A2 . . Change in SOS da Change in SOS da te and te and variatio variatio n of NPN n of NPN from 1982 to from 1982 to 2016 2016 alongside alongside consist consist ee ncy with ncy with Figure A2. Change in SOS date and variation of NPN from 1982 to 2016 alongside consistency with CMG CMGLLSP and SP and VIPPH VIPPHEN_EVI EN_EVI2. N 2. Neeggaative tive va valu lues c es coorrespond to rrespond to earlier earlier da dates (or tes (or fewer fewer day dayss) and ) and CMGLSP and VIPPHEN_EVI2. Negative values correspond to earlier dates (or fewer days) and positive positiv positiv ee values to later date values to later date s ( s ( oor more days). Co r more days). Consistency shows dataset a nsistency shows dataset a gg reement with regards to reement with regards to values to later dates (or more days). Consistency shows dataset agreement with regards to the direction the direct the direct ion of ion of change. Change in date change. Change in date is is the r the r ee gre gre ssed trend ssed trend slope m slope m u u ltipl ltipl ie ie d d by the number of years. by the number of years. of change. Change in date is the regressed trend slope multiplied by the number of years. Change in Change in variation is the reg Change in variation is the reg rr essed essed trend sl trend sl ope of ope of the abso the abso lu lu te valu te valu e of th e of th e residu e residu al al s mu s mu ltiplie ltiplie d by d by variation is the regressed trend slope of the absolute value of the residuals multiplied by the number the number of years. the number of years. of years. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Figure A3. SOS Mean and standard deviation of historical (1982–2016) data from NPN. Rem Rem oo te Sens te Sens . .2020 2020 , , 12 12 , x; , x; doi: FOR doi: FOR PEER PEER RE RE VIEW VIEW www. www. mdpi.com/ mdpi.com/ journal/remotese journal/remotese nn si si ng ng Remote Sens. 2020, 12, 2538 21 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 27 of 33 Figure A4. Change in date and variation of CMGLSP, VIPPHEN_EVI2, and NPN across three phenology Figure A4. Change in date and variation of CMGLSP, VIPPHEN_EVI2, and NPN across three metrics (only SOS for NPN) from 1982 to 2016. Shown here with a more robust color scheme to better phenology metrics (only SOS for NPN) from 1982 to 2016. Shown here with a more robust color indicate change. Negative values correspond to earlier dates (or fewer days) and positive values to scheme to better indicate change. Negative values correspond to earlier dates (or fewer days) and later dates (or more days). Change in date is the regressed trend slope multiplied by the number of positive values to later dates (or more days). Change in date is the regressed trend slope multiplied years. Change in variation is the regressed trend slope of the absolute value of the residuals multiplied by the number of years. Change in variation is the regressed trend slope of the absolute value of the by the number of years. residuals multiplied by the number of years. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 22 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 28 of 33 Remote Sens. 2020, 12, x FOR PEER REVIEW 28 of 33 Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A5. Histogram of CMGLSP and VIPPHEN_EVI2 SOS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Figure A6. Histogram of CMGLSP and VIPPHEN_EVI2 POS change in date from 1982 to 2016. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 23 of 27 Remote Sens. 2020, 12, x FOR PEER REVIEW 29 of 33 Figure A7. Histogram of CMGLSP and VIPPHEN_EVI2 EOS change in date from 1982 to 2016. Figure A7. Histogram of CMGLSP and VIPPHEN_EVI2 EOS change in date from 1982 to 2016. References References 1. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, 1. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. 2013, 169, 156–173. [CrossRef] 2013, 169, 156–173. 2. Chuine, I. Why does phenology drive species distribution? Philos. Trans. R. Soc. B Biol. Sci. 2010, 365, 2. Chuine, I. Why does phenology drive species distribution? Philos. Trans. R. Soc. B Biol. Sci. 2010, 365, 3149– 3149–3160. [CrossRef] [PubMed] 3. Schwartz, M.D. Green-wave phenology. Nature 1998, 394, 839. [CrossRef] 3. Schwartz, M.D. Green-wave phenology. Nature 1998, 394, 839. 4. Hurley, M.A.; Hebblewhite, M.; Gaillard, J.M.; Dray, S.; Taylor, K.A.; Smith, W.K.; Zager, P.; Bonenfant, C. 4. Hurley, M.A.; Hebblewhite, M.; Gaillard, J.M.; Dray, S.; Taylor, K.A.; Smith, W.K.; Zager, P.; Bonenfant, C. Functional analysis of normalized di erence vegetation index curves reveals overwinter mule deer survival Functional analysis of normalized difference vegetation index curves reveals overwinter mule deer is driven by both spring and autumn phenology. Philos. Trans. R. Soc. B Biol. Sci. 2014, 369. [CrossRef] survival is driven by both spring and autumn phenology. Philos. Trans. R. Soc. B Biol. Sci. 2014, 369, 5. Post, E.; Forchhammer, M.C. Climate change reduces reproductive success of an Arctic herbivore through doi:10.1098/rstb.2013.0196. trophic mismatch. Philos. Trans. R. Soc. B Biol. Sci. 2008, 363, 2369–2375. [CrossRef] 5. Post, E.; Forchhammer, M.C. Climate change reduces reproductive success of an Arctic herbivore through 6. Rickbeil, G.J.M.; Merkle, J.A.; Anderson, G.; Atwood, M.P.; Beckmann, J.P.; Cole, E.K.; Courtemanch, A.B.; trophic mismatch. Philos. Trans. R. Soc. B Biol. Sci. 2008, 363, 2369–2375. Dewey, S.; Gustine, D.D.; Kau man, M.J.; et al. Plasticity in elk migration timing is a response to changing 6. Rickbeil, G.J.M.; Merkle, J.A.; Anderson, G.; Atwood, M.P.; Beckmann, J.P.; Cole, E.K.; Courtemanch, A.B.; environmental conditions. Glob. Chang. Biol. 2019, 25, 2368–2381. [CrossRef] Dewey, S.; Gustine, D.D.; Kauffman, M.J.; et al. Plasticity in elk migration timing is a response to changing 7. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response to environmental conditions. Glob. Chang. Biol. 2019, 25, 2368–2381. global change. Trends Ecol. Evol. 2007, 22, 357–365. [CrossRef] 7. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response 8. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.M.; Tucker, C.J.; Stenseth, N.C. Using the satellite-derived to global change. Trends Ecol. Evol. 2007, 22, 357–365. NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 2005, 20, 503–510. [CrossRef] 8. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.M.; Tucker, C.J.; Stenseth, N.C. Using the satellite-derived 9. Zhu, W.; Tian, H.; Xu, X.; Pan, Y.; Chen, G.; Lin, W. Extension of the growing season due to delayed autumn NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 2005, 20, 503–510. over mid and high latitudes in North America during 1982–2006. Glob. Ecol. Biogeogr. 2012, 21, 260–271. 9. Zhu, W.; Tian, H.; Xu, X.; Pan, Y.; Chen, G.; Lin, W. Extension of the growing season due to delayed autumn [CrossRef] over mid and high latitudes in North America during 1982–2006. Glob. Ecol. Biogeogr. 2012, 21, 260–271. 10. Cook, B.I.; Wolkovich, E.M.; Parmesan, C. Divergent responses to spring and winter warming drive 10. Cook, B.I.; Wolkovich, E.M.; Parmesan, C. Divergent responses to spring and winter warming drive community level flowering trends. Proc. Natl. Acad. Sci. USA 2012, 109, 9000–9005. [CrossRef] community level flowering trends. Proc. Natl. Acad. Sci. USA 2012, 109, 9000–9005. 11. Brown, J.F.; Ji, L.; Gallant, A.; Kau man, M. Exploring relationships of spring green-up to moisture and 11. Brown, J.F.; Ji, L.; Gallant, A.; Kauffman, M. Exploring relationships of spring green-up to moisture and temperature across Wyoming, USA. Int. J. Remote Sens. 2019, 40, 956–984. [CrossRef] temperature across Wyoming, USA. Int. J. Remote Sens. 2019, 40, 956–984. 12. Riotte-Lambert, L.; Matthiopoulos, J. Environmental Predictability as a Cause and Consequence of Animal 12. Riotte-Lambert, L.; Matthiopoulos, J. Environmental Predictability as a Cause and Consequence of Animal Movement. Trends Ecol. Evol. 2020, 35, 163–174. [CrossRef] [PubMed] Movement. Trends Ecol. Evol. 2020, 35, 163–174. 13. St-Louis, V.; Pidgeon, A.M.; Clayton, M.K.; Locke, B.A.; Bash, D.; Radeloff, V.C. Satellite image texture and a vegetation index predict avian biodiversity in the Chihuahuan Desert of New Mexico. Ecography 2009, 32, 468–480. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2538 24 of 27 13. St-Louis, V.; Pidgeon, A.M.; Clayton, M.K.; Locke, B.A.; Bash, D.; Radelo , V.C. Satellite image texture and a vegetation index predict avian biodiversity in the Chihuahuan Desert of New Mexico. Ecography 2009, 32, 468–480. [CrossRef] 14. Hurlbert, A.H. Species-energy relationships and habitat complexity in bird communities. Ecol. Lett. 2004, 7, 714–720. [CrossRef] 15. Searle, K.R.; Rice, M.B.; Anderson, C.R.; Bishop, C.; Hobbs, N.T. Asynchronous vegetation phenology enhances winter body condition of a large mobile herbivore. Oecologia 2015, 179, 377–391. [CrossRef] [PubMed] 16. Pettorelli, N.; Pelletier, F.; Von Hardenberg, A.; Festa-Bianchet, M.; Côté, S.D. Early onset of vegetation growth vs. rapid green-up: Impacts on juvenile mountain ungulates. Ecology 2007, 88, 381–390. [CrossRef] 17. Hebblewhite, M.; Merrill, E.; McDermid, G. A Multi-Scale Test Of The Forage Maturation Hypothesis In A Partially Migratory Ungulate Population. Ecol. Monogr. 2008, 78, 141–166. [CrossRef] 18. Middleton, A.D.; Kau man, M.J.; Mcwhirter, D.E.; Cook, J.G.; Cook, R.C.; Nelson, A.A.; Jimenez, M.D.; Klaver, R.W. Animal migration amid shifting patterns of phenology and predation: Lessons from a Yellowstone elk herd. Ecology 2013, 94, 1245–1256. [CrossRef] 19. Stoner, D.C.; Sexton, J.O.; Choate, D.M.; Nagol, J.; Bernales, H.H.; Sims, S.A.; Ironside, K.E.; Longshore, K.M.; Edwards, T.C. Climatically driven changes in primary production propagate through trophic levels. Glob. Chang. Biol. 2018, 24, 4453–4463. [CrossRef] 20. Monteith, K.L.; Klaver, R.W.; Hersey, K.R.; Holland, A.A.; Thomas, T.P.; Kau man, M.J. E ects of climate and plant phenology on recruitment of moose at the southern extent of their range. Oecologia 2015, 178, 1137–1148. [CrossRef] 21. Middleton, A.D.; Merkle, J.A.; McWhirter, D.E.; Cook, J.G.; Cook, R.C.; White, P.J.; Kau man, M.J. Green-wave surfing increases fat gain in a migratory ungulate. Oikos 2018, 127, 1060–1068. [CrossRef] 22. Hamel, S.; Garel, M.; Festa-Bianchet, M.; Gaillard, J.M.; Côté, S.D. Spring Normalized Di erence Vegetation Index (NDVI) predicts annual variation in timing of peak faecal crude protein in mountain ungulates. J. Appl. Ecol. 2009, 46, 582–589. [CrossRef] 23. Bischof, R.; Loe, L.E.; Meisingset, E.L.; Zimmermann, B.; Van Moorter, B.; Mysterud, A. A Migratory Northern Ungulate in the Pursuit of Spring: Jumping or Surfing the Green Wave? Am. Nat. 2012, 180, 407–424. [CrossRef] [PubMed] 24. Merkle, J.A.; Monteith, K.L.; Aikens, E.O.; Hayes, M.M.; Hersey, K.R.; Middleton, A.D.; Oates, B.A.; Sawyer, H.; Scurlock, B.M.; Kau man, M.J. Large herbivores surf waves of green-up during spring. Proc. R. Soc. B Biol. Sci. 2016, 283, 20160456. [CrossRef] 25. Aikens, E.O.; Kau man, M.J.; Merkle, J.A.; Dwinnell, S.P.H.; Fralick, G.L.; Monteith, K.L. The greenscape shapes surfing of resource waves in a large migratory herbivore. Ecol. Lett. 2017, 20, 741–750. [CrossRef] 26. Post, E.; Pedersen, C.; Wilmers, C.C.; Forchhammer, M.C. Warming, plant phenology and the spatial dimension of trophic mismatch for large herbivores. Proc. R. Soc. B Biol. Sci. 2008, 275, 2005–2013. [CrossRef] 27. Mikle, N.L.; Graves, T.A.; Olexa, E.M. To forage or flee: Lessons from an elk migration near a protected area. Ecosphere 2019, 10, e02693. [CrossRef] 28. Rivrud, I.M.; Bischof, R.; Meisingset, E.L.; Zimmermann, B.; Loe, L.E.; Mysterud, A. Leave before it’s too late: Anthropogenic and environmental triggers of autumn migration in a hunted ungulate population. Ecology 2016, 97, 1058–1068. [CrossRef] 29. Mueller, T.; Olson, K.A.; Fuller, T.K.; Schaller, G.B.; Murray, M.G.; Leimgruber, P. In search of forage: Predicting dynamic habitats of Mongolian gazelles using satellite-based estimates of vegetation productivity. J. Appl. Ecol. 2008, 45, 649–658. [CrossRef] 30. Pettorelli, N.; Gaillard, J.M.; Mysterud, A.; Duncan, P.; Stenseth, N.C.; Delorme, D.; Van Laere, G.; Toïgo, C.; Klein, F. Using a proxy of plant productivity (NDVI) to find key periods for animal performance: The case of roe deer. Oikos 2006, 112, 565–572. [CrossRef] 31. Zhang, X.; Tarpley, D.; Sullivan, J.T. Diverse responses of vegetation phenology to a warming climate. Geophys. Res. Lett. 2007, 34, 1–5. [CrossRef] 32. Peng, D.; Zhang, X.; Wu, C.; Huang, W.; Gonsamo, A.; Huete, A.R.; Didan, K.; Tan, B.; Liu, X.; Zhang, B. Intercomparison and evaluation of spring phenology products using National Phenology Network and AmeriFlux observations in the contiguous United States. Agric. For. Meteorol. 2017, 242, 33–46. [CrossRef] Remote Sens. 2020, 12, 2538 25 of 27 33. Keenan, T.F.; Richardson, A.D. The timing of autumn senescence is a ected by the timing of spring phenology: Implications for predictive models. Glob. Chang. Biol. 2015, 21, 2634–2641. [CrossRef] 34. Schwartz, M.D.; Reed, B.C.; White, M.A. Assesing satellite-derived start-of-season measures in the conterminous USA. Int. J. Climatol. 2002, 22, 1793–1805. [CrossRef] 35. Richardson, A.D.; Hufkens, K.; Milliman, T.; Frolking, S. Intercomparison of phenological transition dates derived from the PhenoCam Dataset V1.0 and MODIS satellite remote sensing. Sci. Rep. 2018, 8, 1–12. 36. Klosterman, S.T.; Hufkens, K.; Gray, J.M.; Melaas, E.; Sonnentag, O.; Lavine, I.; Mitchell, L.; Norman, R. Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery. Biogeosciences 2014, 11, 4305–4320. [CrossRef] 37. White, M.A.; de Beurs, K.M.; Didan, K.; Inouye, D.W.; Richardson, A.D.; Jensen, O.P.; O’Keefe, J.; Zhang, G.; Nemani, R.R.; van Leeuwen, W.J.D.; et al. Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982–2006. Glob. Chang. Biol. 2009, 15, 2335–2359. [CrossRef] 38. Richardson, A.D.; Hufkens, K.; Milliman, T.; Aubrecht, D.M.; Chen, M.; Gray, J.M.; Johnston, M.R.; Keenan, T.F.; Klosterman, S.T.; Kosmala, M.; et al. PhenoCam Dataset v1.0: Vegetation Phenology from Digital Camera Imagery, 2000–2015. 2017. [CrossRef] 39. Omernik, J.M.; Grith, G.E. Ecoregions of the Conterminous United States: Evolution of a Hierarchical Spatial Framework. Environ. Manag. 2014, 54, 1249–1266. [CrossRef] 40. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.; Gao, X.; Ferreira, L. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [CrossRef] 41. Crimmins, T.M.; Marsh, R.L.; Switzer, J.R.; Crimmins, M.A.; Gerst, K.L.; Rosemartin, A.H.; Weltzin, J.F.; Jewell, S. USA National Phenology Network Gridded Products Documentation; United States Geological Survey: Reston, VA, USA, 2017; p. 34. 42. Robinson, N.P.; Allred, B.W.; Smith, W.K.; Jones, M.O.; Moreno, A.; Erickson, T.A.; Naugle, D.E.; Running, S.W. Terrestrial primary production for the conterminous United States derived from Landsat 30 m and MODIS 250 m. Remote Sens. Ecol. Conserv. 2018, 4, 264–280. [CrossRef] 43. Meier, G.A.; Brown, J.F. Remote Sensing of Land Surface Phenology; United States Geological Survey: Reston, VA, USA, 2014. 44. Jenkerson, C.B.; Maiersperger, T.; Schmidt, G. eMODIS: A User-Friendly Data Source; United States Geological Survey: Reston, VA, USA, 2010. 45. Friedl, M.; Gray, J.; Sulla-Menashe, D. MCD12Q2 MODIS/Terra+Aqua Land Cover Dynamics Yearly L3 Global 500m SIN Grid V006 [Data Set]; NASA: Washington, DC, USA, 2019. 46. Zhang, X. Reconstruction of a complete global time series of daily vegetation index trajectory from long-term AVHRR data. Remote Sens. Environ. 2015, 156, 457–472. [CrossRef] 47. Didan, K.; Munoz, A.B.; Miura, T.; Tsend-Ayush, J.; Zhang, X.; Friedl, M.; Gray, J.; Van Leeuwen, W.; Czapla-Myers, J.; Jenkerson, C.; et al. Multi-Sensor Vegetation Index and Phenology Earth Science Data Records Algorithm Theoretical Basis Document and User Guide; Vegetation Index and Phenology Lab, University of Arizona: Tucson, AZ, USA, 2015. 48. Richardson, A.D.; Hufkens, K.; Milliman, T.; Aubrecht, D.M.; Chen, M.; Gray, J.M.; Johnston, M.R.; Keenan, T.F.; Klosterman, S.T.; Kosmala, M.; et al. Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery. Sci. Data 2018, 5, 1–24. [CrossRef] [PubMed] 49. Hufkens, K.; Basler, D.; Milliman, T.; Melaas, E.K.; Richardson, A.D. An integrated phenology modelling framework in r. Methods Ecol. Evol. 2018, 9, 1276–1285. [CrossRef] 50. Sonnentag, O.; Hufkens, K.; Teshera-Sterne, C.; Young, A.M.; Friedl, M.; Braswell, B.H.; Milliman, T.; O’Keefe, J.; Richardson, A.D. Digital repeat photography for phenological research in forest ecosystems. Agric. For. Meteorol. 2012, 152, 159–177. [CrossRef] 51. Melaas, E.K.; Sulla-Menashe, D.; Gray, J.M.; Black, T.A.; Morin, T.H.; Richardson, A.D.; Friedl, M.A. Multisite analysis of land surface phenology in North American temperate and boreal deciduous forests from Landsat. Remote Sens. Environ. 2016, 186, 452–464. [CrossRef] 52. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O’Connor, B.; Paganini, M. Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island. Remote Sens. Environ. 2018, 215, 517–529. [CrossRef] Remote Sens. 2020, 12, 2538 26 of 27 53. Sen, P.K. Estimates of the Regression Coecient Based on Kendall’s Tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [CrossRef] 54. De Beurs, K.M.; Henebry, G.M. Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in Kazakhstan. Remote Sens. Environ. 2004, 89, 497–509. [CrossRef] 55. R Core Team, R. A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019. 56. Marchetto, A. rkt: Mann-Kendall Test, Seasonal and Regional Kendall Tests; R Package Version 1.5; Verbania Pallanza, Italy, 2017. 57. Geremia, C.; Merkle, J.A.; Eacker, D.R.; Wallen, R.L.; White, P.J.; Hebblewhite, M.; Kau man, M.J. Migrating bison engineer the green wave. Proc. Natl. Acad. Sci. USA 2019, 116, 25707–25713. [CrossRef] 58. Filippa, G.; Cremonese, E.; Migliavacca, M.; Galvagno, M.; Sonnentag, O.; Humphreys, E.; Hufkens, K.; Ryu, Y.; Verfaillie, J.; di Cella, U.M.; et al. NDVI derived from near-infrared-enabled digital cameras: Applicability across di erent plant functional types. Agric. For. Meteorol. 2018, 249, 275–285. [CrossRef] 59. Liu, Y.; Wu, C.; Sonnentag, O.; Desai, A.R.; Wang, J. Using the red chromatic coordinate to characterize the phenology of forest canopy photosynthesis. Agric. For. Meteorol. 2020, 285–286, 107910. [CrossRef] 60. Borowik, T.; Pettorelli, N.; Sönnichsen, L.; Jedr ˛ zejewska, B. Normalized di erence vegetation index (NDVI) as a predictor of forage availability for ungulates in forest and field habitats. Eur. J. Wildl. Res. 2013, 59, 675–682. [CrossRef] 61. Pettorelli, N.; Ryan, S.; Mueller, T.; Bunnefeld, N.; Jedrzejewska, B.; Lima, M.; Kausrud, K. The Normalized Di erence Vegetation Index (NDVI): Unforeseen successes in animal ecology. Clim. Res. 2011, 46, 15–27. [CrossRef] 62. Zhang, X.; Liu, L.; Yan, D. Comparisons of global land surface seasonality and phenology derived from AVHRR, MODIS, and VIIRS data. J. Geophys. Res. Biogeosci. 2017, 122, 1506–1525. [CrossRef] 63. Lendrum, P.E.; Anderson, C.R.; Monteith, K.L.; Jenks, J.A.; Bowyer, R.T. Relating the movement of a rapidly migrating ungulate to spatiotemporal patterns of forage quality. Mamm. Biol. 2014, 79, 369–375. [CrossRef] 64. Garroutte, E.L.; Hansen, A.J.; Lawrence, R.L. Using NDVI and EVI to map spatiotemporal variation in the biomass and quality of forage for migratory elk in the Greater Yellowstone Ecosystem. Remote Sens. 2016, 8, 404. [CrossRef] 65. Myneni, R.B.; Keeling, C.D.; Tucker, C.J.; Asrar, G.; Nemani, R.R. Increased plant growth in the northern high latitudes from 1981 to 1991. Nature 1997, 386, 698–702. [CrossRef] 66. Zhou, L.; Tucker, C.J.; Kaufmann, R.K.; Slayback, D.; Shabanov, N.V.; Myneni, R.B. Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981–1999. J. Geophys. Res. 2001, 106, 20069–20083. [CrossRef] 67. Berman, E.E.; Bolton, D.K.; Coops, N.C.; Mityok, Z.K.; Stenhouse, G.B.; Moore, R.D. Daily estimates of Landsat fractional snow cover driven by MODIS and dynamic time-warping. Remote Sens. Environ. 2018, 216, 635–646. [CrossRef] 68. Bradley, B.A.; Mustard, J.F. Identifying land cover variability distinct from land cover change: Cheatgrass in the Great Basin. Remote Sens. Environ. 2005, 94, 204–213. [CrossRef] 69. Bradley, B.A.; Mustard, J.F. Comparison of phenology trends by land cover class: A case study in the Great Basin, USA. Glob. Chang. Biol. 2008, 14, 334–346. [CrossRef] 70. Van Leeuwen, W.J.D.; Orr, B.J.; Marsh, S.E.; Herrmann, S.M. Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications. Remote Sens. Environ. 2006, 100, 67–81. [CrossRef] 71. Luo, Y.; El-Madany, T.S.; Filippa, G.; Ma, X.; Ahrens, B.; Carrara, A.; Gonzalez-Cascon, R.; Cremonese, E.; Galvagno, M.; Hammer, T.W.; et al. Using near-infrared-enabled digital repeat photography to track structural and physiological phenology in Mediterranean tree-grass ecosystems. Remote Sens. 2018, 10, 1293. [CrossRef] 72. Morisette, J.T.; Richardson, A.D.; Knapp, A.K.; Fisher, J.I.; Graham, E.A.; Abatzoglou, J.; Wilson, B.E.; Breshears, D.D.; Henebry, G.M.; Hanes, J.M.; et al. Tracking the rhythm of the seasons in the face of global change: Phenological research in the 21st century. Front. Ecol. Environ. 2009, 7, 253–260. [CrossRef] 73. Browning, D.M.; Rango, A.; Karl, J.W.; Laney, C.M.; Vivoni, E.R.; Tweedie, C.E. Emerging technological and cultural shifts advancing drylands research and management. Front. Ecol. Environ. 2015, 13, 52–60. [CrossRef] Remote Sens. 2020, 12, 2538 27 of 27 74. Crimmins, T.M.; Crimmins, M.A.; Bertelsen, D.; Balmat, J. Relationships between alpha diversity of plant species in bloom and climatic variables across an elevation gradient. Int. J. Biometeorol. 2008, 52, 353–366. [CrossRef] 75. Zhang, X.; Goldberg, M.; Tarpley, D.; Friedl, M.A.; Morisette, J.; Kogan, F.; Yu, Y. Drought-induced vegetation stress in southwestern North America. Environ. Res. Lett. 2010, 5, 024008. [CrossRef] 76. Hufkens, K.; Friedl, M.; Sonnentag, O.; Braswell, B.H.; Milliman, T.; Richardson, A.D. Linking near-surface and satellite remote sensing measurements of deciduous broadleaf forest phenology. Remote Sens. Environ. 2012, 117, 307–321. [CrossRef] 77. Zhang, X.; Wang, J.; Gao, F.; Liu, Y.; Schaaf, C.; Friedl, M.; Yu, Y.; Jayavelu, S.; Gray, J.; Liu, L.; et al. Exploration of scaling e ects on coarse resolution land surface phenology. Remote Sens. Environ. 2017, 190, 318–330. [CrossRef] 78. Liu, L.; Cao, R.; Shen, M.; Chen, J.; Wang, J.; Zhang, X. How does scale e ect influence spring vegetation phenology estimated from satellite-derived vegetation indexes? Remote Sens. 2019, 11, 2137. [CrossRef] 79. Baumann, M.; Ozdogan, M.; Richardson, A.D.; Radelo , V.C. Phenology from Landsat when data is scarce: Using MODIS and Dynamic Time-Warping to combine multi-year Landsat imagery to derive annual phenology curves. Int. J. Appl. Earth Obs. Geoinf. 2017, 54, 72–83. [CrossRef] 80. McClelland, C.J.R.; Coops, N.C.; Berman, E.E.; Kearney, S.P.; Nielsen, S.E.; Burton, A.C.; Stenhouse, G.B. Detecting changes in understorey and canopy vegetation cycles in West Central Alberta using a fusion of Landsat and MODIS. Appl. Veg. Sci. 2019, 23, 223–238. [CrossRef] 81. Bolton, D.K.; Gray, J.M.; Melaas, E.K.; Moon, M.; Eklundh, L.; Friedl, M.A. Continental-scale land surface phenology from harmonized Landsat 8 and Sentinel-2 imagery. Remote Sens. Environ. 2020, 240, 111685. [CrossRef] 82. Rodriguez-Galiano, V.F.; Dash, J.; Atkinson, P.M. Intercomparison of satellite sensor land surface phenology and ground phenology in Europe. Geophys. Res. Lett. 2015, 42, 2253–2260. [CrossRef] 83. Czernecki, B.; Nowosad, J.; Jabłonska, ´ K. Machine learning modeling of plant phenology based on coupling satellite and gridded meteorological dataset. Int. J. Biometeorol. 2018, 62, 1297–1309. [CrossRef] 84. Tang, J.; Körner, C.; Muraoka, H.; Piao, S.; Shen, M.; Thackeray, S.J.; Yang, X. Emerging opportunities and challenges in phenology: A review. Ecosphere 2016, 7, 1–17. [CrossRef] 85. Johnston, A.N.; Beever, E.A.; Merkle, J.A.; Chong, G. Vegetation responses to sagebrush-reduction treatments measured by satellites. Ecol. Indic. 2018, 87, 66–76. [CrossRef] 86. Schwartz, M.D. (Ed.); Phenology: An. Integrative Environmental Science; Springer: Dordrecht, The Netherlands, © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Journal

Remote SensingUnpaywall

Published: Aug 7, 2020

There are no references for this article.