On the Use of Standardized Multi-Temporal Indices for Monitoring Disturbance and Ecosystem Moisture Stress across Multiple Earth Observation Systems in the Google Earth Engine

: In this work we explore three methods for quantifying ecosystem vegetation responses spatially and temporally using Google’s Earth Engine, implementing an Ecosystem Moisture Stress Index (EMSI) to monitor vegetation health in agricultural, pastoral, and natural landscapes across the entire era of spaceborne remote sensing. EMSI is the multitemporal standard (z) score of the Normalized Difference Vegetation Index (NDVI) given as I, for a pixel (x,y) at the observational period t. The EMSI is calculated as: z xyt = (I xyt − µ xyT )/ σ xyT , where the index value of the observational date (I xyt ) is subtracted from the mean ( µ xyT ) of the same date or range of days in a reference time series of length T (in years), divided by the standard deviation ( σ xyT ), during the same day or range of dates in the reference time series. EMSI exhibits high signiﬁcance (z > |2.0 ± 1.98 σ |) across all geographic locations and time periods examined. Our results provide an expanded basis for detection and monitoring: (i) ecosystem phenology and health; (ii) wildﬁre potential or burn severity; (iii) herbivory; (iv) changes in ecosystem resilience; and (v) change and intensity of land use practices. We provide the code and analysis tools as a research object, part of the ﬁndable, accessible, interoperable, reusable (FAIR) data principles.

Before the inception of Google's Earth Engine (GEE) [21], the ability to interrogate decades of EOS data spanning multi-petabyte sized archives stored at numerous separate data centers (with little or no publicly available computation capacity) was effectively impossible. Today, GEE users can explore anywhere on the Earth's surface in near realtime, across the entire history of space-based remote sensing. Researchers can compare observations to trends and changes seen in instrumental data, weather interpolations, and climate projections.
We have three integrated objectives in this paper: (1) introduce a novel application of multi-temporal standard (z) scores as both an "Ecosystem Moisture Stress Index" (EMSI) for measuring and monitoring ecosystem health and for detecting rapid changes from disturbance, over the time scale of EOS; (2) provide empirical applications for the EMSI with data from legacy and continuance EOS missions over the past 36 years; and (3) take an open-science approach whereby the data and analysis code follow the findable, accessible, interoperable, reusable (FAIR) data principles [22] for replication of results and new analyses by others.
This paper is organized as follows: (1) a perspective on the use of standardized indices and a brief history of past work that led to the formulation of the EMSI; (2) a description of the calculation of EMSI at the pixel scale and raster scale, as well as in GEE; (3) examples of EMSI in different ecosystems across the world; and (4) a discussion about the limitations and a conclusion with potential future applications of the EMSI.

Justification
Current challenges to measuring the 'health' of an ecosystem from EOS remote sensing relate to the establishment of what is the 'normal' or 'average' condition for a given location at any point in time in the present or the past. The justification for a multi-temporal standardized score is that it characterizes ecosystem health in terms of standardized greenness, rather than relative greenness. Areas that are always green, e.g., rainforests, or only briefly reach their peak greenness, e.g., desert grasslands, change dynamically, requiring a standardized value for their health within a season or period of time.
The availability of legacy remote sensing data for Earth systems study is both an asset and a liability. Decades of spaceborne remote sensing supply a rich source of data on systems dynamics; but spectral, spatial, and temporal resolutions are often mixed. The application of multitemporal standard scores serves to address such challenges and limitations in remote sensing data. In the present study, we apply a multitemporal standard score to ecosystem monitoring, adopting the term "Ecosystem Moisture Stress Index". The EMSI represents a subset of potential applications of multitemporal standard scores, which we believe can apply more broadly to land use changes produced, for example, by urbanization and wildland fire disturbances. Our chief aim here is monitoring effects of climate forcing (i.e., temperature changes and precipitation deficits) through time. Variations in regional climate affect ecosystems differently at annual, multi-year, and decadal time scales [23,24]. Interannual variations in temperature and precipitation mediate ecosystem function and can emerge over time as climatic events. 'Hot droughts' are climatic episodes both drier and warmer than the long-term average. Hot drought increases stress on plants through elevated vapor pressure deficit [25][26][27][28], in addition to decreased overall moisture availability and photosynthetic efficiency [28]. Pluvial, in contrast, are climatic events both wetter and cooler than long term averages, leading to more robust vegetation and greater growth rates over decadal or longer time scales [29,30]. While evidence of global warming exacerbating ecosystem drought stress is emerging from remote sensing studies at the interannual time scale [31][32][33], new tools are needed to attribute changes in ecosystem health to new climate extremes, episodic disturbance, or part of the natural range of variability in any given locale.
Establishment of whether a location or an area are within their historical range of variability requires a long-term EOS based time series of that area, normalized across multiple platforms with varying spectral, spatial, and temporal resolutions.
Short term trajectories of plant physiological conditions are generally regulated by weather and radiative forcing, within the characteristic seasonal climatic variability in their local range [30,34]. Over semi-to decadal time frames these seasonal variations may also include climatic extremes such as pluvial or drought, as well as disturbances and their legacies [17,29,32,33]. Prediction of the prevailing health of vegetation in an agricultural field, pasture, or natural landscape requires understanding present and past ecophysiological conditions, as well as their climatic variability. Direct stressors of plant condition include chronic and episodic heat, changes in radiative forcing, increased atmospheric evaporative demand, and soil moisture deficit from drought [25,26,35]. Wildfires, disease, herbivory, and insect outbreaks are additional sources of plant stress [19,32,33,36].

Fuel Moisture Stress Index & Standard Vegetation Index
The use of multi-temporal standardized indices in remote sensing of vegetation health goes back at least to the early 2000s. The fuel moisture stress index (FMSI) [37] was conceived originally as part of a NASA Regional Earth Sciences Applications Center cohosted by The University of Arizona and the University of California, Berkeley. The FMSI is a spectral proxy of vegetation moisture stress, derived from remotely sensed data and thought to represent effects of climate forcing on live vegetation (i.e., fuel) moistures. It was tested using Advanced Very High-Resolution Radiometer (AVHRR) [38] data to map regional-scale fuel moisture and enhance wildfire burn scars [37]. The FMSI has been applied as a means to assess live fuel moisture deficits that contribute to the Energy Release Component (ERC) [39,40] and potential flammability of vegetation during a fire season [41]. Similarly, Peters et al. [42] derived a Standard Vegetation Index (SVI) for drought monitoring using the same temporal z-score framework as Yool [37] with AVHRR data. The Yool [37]  In the present work, we expand these original motivations for the FMSI in fire management, and the SVI for drought monitoring, to reflect more broadly the open science applications of standard z scores for all EOS legacy data. We improve upon the foundational work cited above by operationalizing the community's ability to calculate multi-temporal z-scores for a range of spectral indices, including the Enhanced Vegetation Index (EVI) and Normalized Difference Water Index (NDWI), in Google's Earth Engine at any resolution across a range of dates for the entire EOS record until the present day.

An Open Science Approach
Data-intensive scientific discovery [44] requires a reproducible research approach which becomes a component of the "research object" (RO) [45,46]. For digital reproducibility, data should be made available; protocols should be documented in computational notebooks, version control systems should be used to track software and code changes over time, and versioned software containers should be hosted. With publicly available cloud computing services such as CyVerse [47][48][49], and a version control system (GitHub), other scientists can reproduce our results using our example code and R notebooks; or otherwise calculate multi-temporal z-scores anywhere else using the data and scripts provided for GEE. We follow an open-science approach whereby the data and analyses attempt to follow the FAIR data principles [22] in the effort to support a complete RO.

Earth System Observations
The term ecosystem phenology refers to mean seasonal variability and by proxy the viability of ecosystem health [50,51]. Variations in ecosystem phenology over multi-annual observation periods (e.g., five to ten years) may under-sample for long-term climatic variation. For example, semi-decadal and multi-decadal oscillations in synoptic scale meteorology [23,36,51] require longer time series to capture entire cycles. To determine whether a credible assessment of ecosystem phenology from climate forcing can be accomplished, we compiled up to forty years of legacy EOS data from 1977 until the present, for our computation of EMSI at different locations. The National Aeronautics and Space Administration (NASA) constellation of EOS satellites offer nearly five decades of remote sensing data ( Figure 1, Table 1). GEE (re)hosts both NASA and European Space Agency (ESA) data via its application programming interface (API) and web-hosted graphic user interface (GUI) ( Table 1).  EOS data that stretch back to the early 1970s are available at a variety of spatial, temporal, and spectral resolutions ( Figure 1, Table 1). While data are available directly from NASA data services for Landsat [52], MODIS [53], AVHRR [38,54], and ESA data service for Copernicus (Sentinels) [55,56], Google Cloud provides these data as image collections for its GEE cyberGIS. Researchers can also upload their own data to GEE via the Google Cloud Platform and create private or public collections. For example, commercial PlanetScope data are available from Planet Labs [57] at higher spatial and temporal resolution than other public EOS data, and can be utilized in GEE as private collections.
Gridded climate data hosted on GEE include: Daymet [58], PRISM [59], World-Clim2 [60], and NASA NEX-DCP30 [61]. These data help to establish the magnitude of past climate forcing at the level of an individual pixel in GEE. In our code repository Appendix A, we provide examples of how to extract GridMet data [62] and compare a wildfire danger index and energy release component (ERC) to the EMSI.
Our analyses used the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) Surface Reflectance (SR) product [63][64][65] for Landsat 4, 5, 7, and Landsat Surface Reflectance Code (LaSRC) SR product for Landsat 8 [66]. For Sentinel-2 we use the Level-1C data which support cloud masking. MODIS (Terra and Aqua) Vegetation Indices 16-Day Global 250 m, VIIRS daily surface reflectance (VNP09GA) Version 1 product [67], AVHRR Normalized Difference Vegetation Index (NDVI) (NOAA CDR v4) [38], also have bit flags for cloudy scenes. Surface reflectance products provide an atmospherically corrected estimate of vegetation spectral response over time, unlike Top of Atmosphere (TOA) or raw imagery. Surface Reflectance Landsat data in GEE also supports the removal of pixels with cloud, cloud shadow, water, and snow via quality assessment (QA) bit-masks [66]. The Landsat LEDAPS and LaSRC products in GEE have bit flags for cloud, cloud shadow, snow, and water. We used those flags to remove pixels with altered NDVI values. These QA masks were applied to the collections to remove outlier pixels before sampling NDVI.

Normalized Difference Vegetation Index (NDVI)
The Normalized Difference Vegetation Index (NDVI) [68] is a well-established linear combination of the red and infrared radiometric bands (Table 1). Over the last 40 years, NDVI has become a key metric for monitoring vegetation health and moisture stress [69,70], land-cover change [5], wildfire location and burn severity [12], agricultural health and production [71], and to track climate change across ecosystems [72][73][74][75]. The NDVI is a derived measure of 'greenness' within a given pixel, equation: where ρNIR is reflectance in the near infrared band (typically between 725 and 1100 nm) and ρRED is the reflectance in the red band (typically between 580 and 684 nm) ( Table 1). Within any vegetation type, NDVI varies across time depending upon composition and complexity, phenology, and cover percentage relative to the pixel size. Peak values for NDVI are generally observed following wet seasons and the period of greatest net primary production (NPP), which corresponds generally to the seasonal precession of the Earth. Minimum NDVI values are reported in periods of senescence, e.g., winters or dry seasons. NDVI is also used to detect interannual variability in greenness and NPP.

Calculating EMSI from NDVI Time Series
The EMSI is calculated for individual pixels during the same Julian day or time period for every year in a time series (Figure 2). The EMSI takes the equation: where, x xyt is NDVI at the observed date or time period t of interest of pixel xy, µ xyT is the mean of the NDVI time series T of the same date or period nearest the observation time step t, and σ xyT is the standard deviation of the NDVI time series or period of the pixel xy. The normalization coefficient ρ is positive; thus, a high EMSI value corresponds to robust greenness beyond the pixel-scale average ( Figure 2). EMSI is essentially the standard score, or z-score, of an observed NDVI xy at time period t equation: where z is the standard score and µ and σ are the NDVI mean and standard deviation of the population, respectively. In traditional inferential statistics, standard scores facilitate analysis of variations within an x,y matrix of values, in which the mean and standard deviation are derived from all cells in the spatial plane of the matrix. In EMSI, we use standard scores of the temporal plane: the values used to represent the mean and standard deviation are derived from values in a single cell of this matrix through time, where the matrix is composed of pixels in a remotely sensed image. An EMSI score of ±1 corresponds to one standard deviation from that pixel's time series average; thus, a score of ±1.96 is a 95% confidence interval (CI), and a score of ±2.58 is 99% CI. We propose the use of benchmark values in establishing the condition of an ecosystem over the reference time series as for example: stable (−1.0 < z < 1.0) characterized as having no significant variation from the average, degrading (−2.0 < z < −1.0), collapsing (z < −2.0), improving (1.0 < z < 2.0), and robust (e.g., 2.0 < z) ( Table 2). Pixel scale calculations of the EMSI with smoothing were done in the RStudio [76,77], which was run using Docker [78,79], with Linux-based images released by the Rocker Group [80,81]. Analyses in R were run on remote workstations and on the CyVerse cloud service [82].
For the pixel scale analyses, we used locally weighted regression (loess method; [101], to calculate a moving average (mean) and standard deviation of the NDVI for each EOS time series. Reference periods for estimating the natural variability of the system are dependent upon the temporal and spatial properties of the system of interest.

Computations in the GEE CyberGIS versus Local Computing
We pixel-and area-sampled EOS platform time series data for their surface reflectance adjusted NDVI using the GEE CyberGIS framework. This protocol enabled us to explore any location on the Earth's surface over the duration of EOS operation. GEE empowers scientific analyses through its extensive collection of algorithms and filtering routines. We wrote an interactive 'Explorer' window ( Table 3) that enables the user to select a single pixel, a sub-pixel scaled reduction, or a polygon 'Area of Interest' (AOI) from any of the EOS collections in the GEE web browser interface ( Figure 3). In depth details about our GEE analyses and tutorial workflows, including code, are provided in our Appendix A on Github.  We examined three methods of calculation for the multi-temporal standardized z index: (1) at a pixel scale, using time series from GEE pixel extractions of surface reflectance normalized NDVI at a single locale; (2) at area scale directly in the GEE without extracting data; and (3) at area scale with rasters extracted from GEE and calculated in R using GDAL. The pixel based NDVI extraction uses a daily interpolation of NDVI between periods (dates) by using a smoothing spline. The raster based NDVI uses mosaicking of scenes at weekly to monthly periods. Surface reflectance NDVI from the Landsat 5, 7, and 8, and Sentinel-2 were exported from GEE as CSV and GeoTiff ( Figure 3, Table 3, Appendix B). EMSI is also calculated directly in GEE using. .js (Appendix A). For our open-science approach, we hosted our analysis on a version control system, GitHub (https://github.com/tyson-swetnam/emsi, accessed 11 March 2021). Our analysis workflows ( Figure 3) are written in a combination of JavaScript (.js), Python (.py), and R Markdown (.rmd) [99]. The Earth Engine (EE) Code Editor (http://code.earthengine. google.com, accessed on 11 March 2021) is a web-based integrated development environment (IDE) for the Earth Engine JavaScript API. We use .js files as extractors for pixel sampling the EOS legacy data ( Table 2). We use .py scripts to extract select location data from GEE and moved these data onto CyVerse [48] Data Store. Virtual machines running instances of RStudio-Server [76,82] on CyVerse, can be used to re-analyze these data using our provided R Markdown notebooks. We provide instructions in our Appendices A and B for launching extraction and on-the-fly analyses of the EMSI in GEE, as well as instructions for re-running our statistical analyses using an NSF supported cloud native data science workbench, CyVerse, or other platforms such as National Data Service Labs Work Bench or the Google Colaboratory.

Study Areas
We selected four model ecosystems (Figure 4, Table 4) across the planet to demonstrate how annual variations in NDVI differ for the resulting EMSI. Within each ecosystem we selected two sites with different seasonal phenology, vegetation, or land use practice. In Botswana, Africa, we selected sites near the Okavango River Delta: an upland thornscrub area and a lowland wetland area in the river delta. In North America, we selected a semi-arid grassland separated by an international border wall between Arizona, United States, and Sonora, Mexico with different cattle grazing practices. In Acre, Brazil, we selected an area of rural agriculture and nearby primary tropical rainforest. In the Sahka Republic, Russian Federation, we selected a taiga forest composed of deciduous (larch) and evergreen (pine) conifers and an 'alaas' thermokarst meadow consisting of grasses, lichens, and sedges [102].

EMSI Spatial Variablity over Time
With up to a maximum of 36 years for designating the reference period across our four example ecosystems, the resulting EMSI generally ranges between z ± 2.0 (Equation (3)) with some extreme z ± 3.0 values (Figures 5 and A1). Extreme values both positive and negative tended to be associated with major type conversions, e.g., rainforest to agriculture, extreme drought, or recent wildfire events ( Figure A2, Appendix B). For the pixel scale extractions, time series are interpolated between dates using a smoothing function in R, allowing a daily estimate of EMSI across the year ( Figure 5). In regions with winter snow cover, the EMSI was not calculated. This is due to the pixel flag removal for 'white' flagged pixels; this problem was most dominant only in the Yakutia site in the present study but would be evident in any snow-covered surface globally.

Seasonal Flooding of an Inland Delta Lowland versus Localized Drought in a Thornscrub
Upland. Botswana, Africa NDVI in the Okavango ecosystems exhibited a sinusoidal curve related to seasonal rainfall ( Figure 5), the curve has greater variance in the upland site than the lowland site (which remained greener year-round; Figure A3). EMSI in the Okavango exhibits greater variation during austral summer, corresponding to the summer rainy season ( Figure 5). Over the 34-year reference period ( Figures 5 and A3), there was ± 0.15 NDVI from the long-term average for the Landsat series ( Figure 5, black line). The resulting EMSI ranges between z ± 1.8 (Equation (3)) with extreme values outside z ± 1.8 (Figures 5b and A4).
EMSI captures the variability of vegetation heatlh in areas with low NDVI, e.g., the southerly portions seen in Figure A3 are consistently much lower in March than the areas north in the swampy delta lowlands. In the EMSI however, this variation is revealed as a much wider range (Figures 6 and A4). An ancient linear dune formation with thornscrub and savannah uplands to the west of the Okavango delta exhibits a different EMSI pattern over the time series than the lowland areas of the delta ( Figure 6). Periods of drought and wet surplus were also different for the lowland versus the upland sites. The lowland (swamp) has a greater amplitude over the reference period than the upland (shrubland) (Figure 7a

Grazing Pressure in a Semi-Arid Grassland, Arizona, USA and Sonora, Mexico
NDVI in the US-Mexico border region exhibits a modal bump; the variation was nearly identical on the US side and Mexican side ( Figure A1), reflecting the identical vegetation and climate of the two points which are spatially close together but separated by a border wall. The resulting EMSI exhibited greater variation during summer ( Figure A1) corresponding to North American monsoon summer rainfall. EMSI ranged between −1.8 < z < 2.4 (Equation (3)) with extreme values approaching z < 2.5 in very wet summer monsoon conditions ( Figure A2, Appendix B). Over the 36-year reference period (1984-2020), years with severe drought (i.e., 1994, 1995, 2005) and very wet periods (1990,1999,2006,(2015)(2016)) epitomized the wide range of variation in this ecosystem (Figures 8 and 9). On the Mexican side of the international border, there appears to be higher grazing pressure (lower maximum NDVI), as further evidenced by the horizontal demarcation line across the center of the scene (Figure 8, lower right panel) representing the border fence (now, "wall"). Individual fenced pastures on the US side of the border can also be differentiated for grazing pressure by year ( Figure 8). In some years, pastures appear to be ungrazed, while in other years they are robustly high in their EMSI, reflecting a strong monsoon. In other years, specific pastures appear degraded, indicating their use for cattle grazing during weak or failed monsoonal rainfall.  Areawide EMSI histograms which take into account all of the pixels in the Landsat scene for that period indicate a normal distribution with a moving average above and below the z = 0 (Figure 9). Years with above normal to extreme z scores correspond to years in which the seasonal (July-August) monsoon rainfall was below or above average. Deeply stressed EMSI years, e.g., 1994, 1995, or 2005, represent years in which the monsoon was late, below average, or both. Exceptional positive EMSI years are years with above average monsoonal moisture.

Drought in Rainforests and Rapid Agricultural Expansion in a Tropical Rainforest, Acre, Brazil
The NDVI of rainforest in the state of Acre, Brazil region exhibited continuous high values > 0.6, with the exception of some deep drought years ( Figures A1 and A2). The EMSI exhibited the highest variation during the dry season, which is also one of the only cloud free times in this region. The resulting EMSI ranged between −2.2 < z < 1.5 (Equation (3)) with extreme values approaching −3 < z. The reference period was set between 1986-2020 ( Figure A6). In 1993, and again in 2010, Brazil underwent severe droughts [104,105] all of these are captured in the EMSI (Figures 10 and A7).
In our EMSI series, the impact of the droughts on the surrounding rainforest and agricultural pastures is most evident in the spring (Figure 10). The agricultural areas began to be cleared in the mid 2000s, with rapid expansion up through the current year (2020) ( Figure A6). When we establish an EMSI reference period before the clearing for agriculture, these areas appear as 'disturbances' which have much lower NDVI, and therefore extremely low EMSI. Evaluating the agricultural pixels alone, reveals that they have greater sensitivity to changing conditions than the surrounding rainforest ( Figure A2).

Winter and Repeated Wildfires in a Taiga Forest, Yakutia, Russian Federation
NDVI in the Yakutian region exhibits a rapid green up following winter snow melt around Julian day 130, with values jumping from near 0.0 to 0.85, and then rapidly declining again in the fall prior to the return of winter ( Figure A1). The EMSI of the example areas exhibited the greatest variation during the summer growing season between −2.2 < z < 1.0 (Equation (3)). Yakutia had limited years of Landsat data availability in GEE collections (~1999-2020) ( Figure A8). This area was particularly also affected by the lack of cloudfree and snow-free dates for imagery. Those dates that typically fell between Julian days 110-310, which corresponded to the boreal summer and the absence of snow cover on the surface ( Figure A2). Despite the shortened availability of imagery, significant seasonal responses are observed. Multiple wildfire scars are visible in the center and the east and south of the scene. This area experiences multiple fires thoroughout the observed time series (Figure 11a). The alaas (herbaceous sedges and reindeer moss lichen dominated open meadows) were more sensitive in their EMSI response to drought years than the surrounding coniferous forest ( Figure A2, Appendix B). Changes in river depth of the Vilyuy, visible in the northwestern edge of the scene, reveal the migration of a sand-bar oxbow following seasonal ice breakup (Figure 11b). The EMSI becomes more apparent with the green up of algae in the river. The numerous wildfires that burned through the Scots pine forests are visible in the center of the scene (Figure 11), many of which exhibit classic Huygens ellipses, indicative of a wind-driven crown fire.

Discussion
The EMSI captures and describes variation in NDVI independent of: (1) long-term changes in vegetation condition related to succession; (2) spatial complexity of vegetation in complex topographic relief; (3) exogenous disturbances, e.g., wildfire or land clearing; and (4) ecosystem phenology relative to climate forcing. Our worked examples demonstrate that these responses may in many cases depend in part upon long term changes in climatology (e.g., periodic wet or dry periods or other periodic variation in the ocean-atmosphere system), land use practice (e.g., agriculture or grazing versus native ecosystem), and disturbances (wildfires and insect outbreaks).

Establishing the Appropriate Reference Period
For the highlighted ecosystems and land use cases presented, an ordinal ranking criterion to characterize the trajectory of each ecosystem is used to contextualize their ecosystem phenology (Table 2). This criterion is dependent upon the establishment of the reference period.
The establishment of each reference period is critical to assessing natural variation in ecosystem health [106]. Ideally, the baseline period for EMSI should be representative of a given pixel's variability but be independent of disturbances and extreme departures. For example, setting the reference period during a mesoscale climate anomaly, e.g., a pluvial, a very strong El Niño/La Niña, or immediately after a major disturbance, will alter the z-score at any given time of interest before or afterwards. As climate and ecological events are asynchronous across the global surface, analyses may benefit from site-specific reference periods in lieu of a single global reference period covering the entire epoch of EOS.
Choosing the appropriate reference period can be specific to a particular research question. Researchers may be interested in comparing vegetation health from the 1970s or 1980s to the 2000s and 2010s, or vice versa. A reference period set in the 1980s for the southwestern US (one of the wettest decades on record [107]) would result in z-scores that are exceptionally low during the 2000s and 2010s (a historically hot and dry period [107,108]). Similarly, a reference period set in the southwestern US in the 1990s would show less significant change to the 2010s but show the 1980s as being scored as exceptionally wet. Given the entire record of EOS, the choice of reference period must also take into account known disturbance histories, such that these locations do not hide the intensity of historical disturbance events. If the goal is to establish variations in disturbance impacts over the same areas, then the best practice may be to use the maximum time period available.

Ecosystem Phenology and Disturbance Detection
Seasonal variation in the NDVI from intermittent flooding in the Okavango River Delta is representative of regional influences on precipitation and temperature in the neighboring Angolan and Namibian portions of the river basin, relative to the upland thornscrub areas in Botswana above the delta, which is dependent upon local precipitation. In the Okavango river basin, the year 1987 was an average water year for the upper watershed in Angola [109], but drier than average conditions for the upland thornscrub areas surrounding the delta in Botswana (Figure 5a). The resulting EMSI for March of 1987 shows significant declines in vegetation health on the dry uplands, relative to the slightly above average conditions in the river fed delta. The 1994-1996 period was drier than average for the Okavango River Basin (http://okavangodata.ub.bw/ori/monitoring/ water/, accessed on 6 April 2021). During these dry periods there are minor fluctuations in NDVI in the delta. In that example, EMSI was more responsive to these drought signatures (Figure 6a,c) than is the NDVI. EMSI in the delta responded positively to the following period of wetness in water years 2000-2001 (Figures 5b and 7b). We found the thornscrub upland is less sensitive to these variations in its EMSI, which may be attributed to local rainfall patterns (Figure 7a). Periodic wildfires burning throughout the Okavango scene are also apparent throughout the time series as significantly negative EMSI values in particular the years 2002, 2003 ( Figure A4).
In the North American case study, a recently expanded and relatively impassable fence divides the US-Mexico international border creating an experiment. The climate, geology, plant species, and ecological conditions are highly similar, but local ranchers graze cattle on either side of the border differently. Consequently, variations in NDVI and the EMSI should relate to local weather and topography. Specifically, differences in EMSI that are apparent in the figure are due to variation in the land use. In this portion of North America, late August is associated with peak greenness related to the North American Monsoon. EMSI is less sensitive to topographic-induced variations than NDVI; for example, shrubs and agriculture along drainage bottoms vs. grass covered slopes and ridgelines (Figure 8). The EMSI reveals which pastures (or paddocks) cattle are grazing in during this period of high greenness, as evidenced by NDVI (Figure 8c).
In the Amazon basin, the NDVI of rainforest is generally very high (>0.8) for the entire year ( Figure A1). The EMSI is however better able to characterize minute variations in the rainforest's NDVI changes by calculating its variability during the dry season scenes shown in detail in Figure 10. Because of the heavy cloud cover in Acre state, the winter dry season is the only time of year we have to view cloud-free images over these sites. The influence of the 2010 heatwave is apparent in both the rainforest and the agricultural areas in Figures 10, A6 and A7.
The influence of heatwaves and drought are present in the Yakutian thermokarst Alaas meadows, but less so in the Taiga forest areas ( Figure A2). Periodic wildfires have burned across the reference period, altering the estimated EMSI for the burned areas to be positive prior to the fire, and negative afterward (Figure 11a).

Utility of a Multitemporal Standardized Index
While NDVI provides valuable information about spatial and temporal variation in greenness, its interpretation can be problematic because productivity varies widely in different ecosystem types and climate conditions [110][111][112]. For example, arid and semi-arid ecosystems are invariably lower in biomass than high-productivity tropical or temperate forests, so the proportional variation in NDVI is not evident. This is the motivation for a standardized index, versus a relative greenness index, which characterizes the variation in NDVI.
The relative greenness (RG) index [113,114] and departure from average (DA) [37,42] are descriptive statistics originally developed for drought and wildland fire applications [115,116]. RG and DA estimate inter-pixel moisture in live vegetation by relating current NDVI observations to average NDVI across a time series. RG compares NDVI values (maximums and minimums) for a given period, across the time series, to the NDVI value for that period. DA divides the per pixel NDVI for a given period by the per pixel mean NDVI for that period, across the time series. In addition to indicating high or low values for a given pixel, relative values over the same periods for multiple years distinguish changes in NDVI not due to annual phenological patterns. RG and DA indices have been useful in estimating live vegetation moisture stress in some ecosystems. However, it has proven difficult to set thresholds for these indices, leading to inaccurate estimates [117]. A multi-temporal standardized index, by contrast, enables analysts to establish such thresholds using standardized scores.
The EMSI rescales NDVI values into standardized greenness: unlike RG and DA, the EMSI for each pixel can be used to build a temporally determined distribution of standard (i.e., z) scores ( Figure 2). The two most significant advances contributed by the EMSI are as follows. First, the EMSI standardizes spectral differences that arise from variations in the position and width of NIR and red spectral bands among different sensor systems ( Figure 1, Table 1). This enables analysts to mix sensor system data in a time series analysis. Second, the multi-temporal standard scores enable comparative studies of stress between different ecosystems. Rescaling to standardized indices establishes thresholds that enable comparisons of stress among different ecosystems [43]. Although we base the EMSI on the well-established NDVI, other spectral indices, e.g., the normalized difference water index (NDWI) [70,118], which is sensitive to the water content in vegetation, or EVI, could be used to formulate analogous standardized stress indices (see Appendices A and B).

Limitations
For scene-level data extracts the number of pixels (e.g., one to twenty-five million pixels per raster) with a daily scale interpolation would be prohibitive in GEE over the entire global surface. These can be run in other cyberGIS, but these will require significant data storage hosting and compute resources.
The EMSI is not, and will not, be available as a pre-calculated image collection within GEE. The number of independent factors listed above, including the establishment of a sitespecific reference period which is unbiased by disturbance or mesoscale climate anomaly and the interpolation of NDVI between annual repeat dates which do not align exactly, variously constrain the possibility of generating a static and global EMSI product. EMSI must therefore be parameterized and calculated dynamically once the researcher sets the prior parameter values. Future inclusion of dynamically calculated EMSI and derivative multitemporal z-score based measures of alternate spectral indices in a cyberGIS portal, such as the Climate Engine [119], or GEE would be highly beneficial to the community.
For Landsat and Sentinel data, QA bit flags were useful in removing a large quantity of noise, although this process was imperfect. Some pixels still have shadowing or other artifacts present within them, resulting in variations in the NDVI per pixel. The Landsat MSS 1, 2, 3, and Landsat 4 had unresolved issues related to clouds and shadows which we were not able to correct in our code. Future work that cleans up those earliest data to derive a useful SR product that can be compared to the rest of the Landsat series will strengthen the computation and use of EMSI over longer time series.
The use of non-parametric statistical derivations for establishing EMSI may yield more sensitive results than ours in the present study. For example, Anees & Aryal [120] developed a non-parametric method using z-scores and MODIS data to delineate beetle infestation in pine forests. Rather than apply the mean (Equation (3)), as we have described, the application of the median may be necessary in some systems and conditions. Theirs and others' application of non-parametric statistics for evaluating EOS data, such as from MODIS [120][121][122] and Landsat [123], could provide more precise and earlier detection of variation from historical averages. The application of these methods warrants future investigation in the GEE.

Conclusions
The accelerating progression of ecosystem degradation in response to changes in the global climate system requires the ability to monitor and track ecosystem condition at a simultaneously high spatial and temporal resolution. The EMSI introduced in this paper is a statistically straightforward technique for capturing, within the GEE dataspace, ecosystem dynamics at fine temporal scale and for characterizing health trajectories over time. The GEE-hosted EMSI can for example be integrated with gridded climate [58,59,61] and wildfire danger index data [62] to monitor moisture stress in vegetative fuels (Appendix A). We believe going forward that such integrated technology can be generalized for a variety of earth and life science applications within the GEE cyberGIS.
In this paper we featured an open science approach using a version control system (GitHub), unifying data from legacy and contemporary remote sensing systems in the GEE cyberGIS. We demonstrated in four model ecosystems how the GEE tools and popular open-source software packages for statistical analyses can be used to derive a seamless spectral trajectory of ecosystem performances.
Among several future applications of the EMSI, ecosystem recovery following disturbances such as wildland fire is exemplary because fire is such an important component of ecosystem health; fire seasons are starting sooner, ending later, and burning with greater severity than in at least the past century [24,[124][125][126][127]. The growing length of fire seasons has been linked to antecedent weather conditions including: (a) quantity and phase (snow vs. rain) of winter precipitation [128,129]; (b) increasing spring and summer temperatures [130]; (c) increasing vapor pressure deficits [131,132]; and (d) hot drought [17,18,30]. Each of these changes represents a potential use case for multi-temporal standardized indices such as the EMSI.
Forests are degrading globally, due chiefly to the combined effects of severe moisture stress related to hotter droughts [18,133,134] and increasing landscape disturbance. Such climate forcing is forecast to occur more frequently and over larger areas in the coming years and decades [27]. Using EMSI to establish the rates of change in these ecosystems is an important area of potential future research. For example, forests that are drought stressed may in addition express greater severity for a given fire intensity [135,136]. Adapting EMSI to formulate a family of indices, such as NDWI, coupled with finer scale precipitation and temperature data, should contribute dimensional robustness to EMSI as an important analytical tool.     The years 1989 and 2008 were below average for the delta (which is typically green in the NDVI), but above average for the surrounding uplands. The three highlighted boxes are shown below in more detail in Figure 8. Pixels with no data due to masking are gray.  Figure 11. Pixels with no data due to masking are gray. Figure A6. EMSI for 1986-2020 for Landsat 5, 7, and 8 time series for the month of August in Acre state, Brazil. The years 1988,1989,1995,1996,2002,2009,2012, and 2014 had no clear sky scenes. The years 1993The years , 1993The years , 2001The years , 2005, and 2010 were historic drought years in Brazil and appear as highly stressed. The region undergoes rapidly expanding conversion to agriculture in the last decade of the time series. The three years with blue highlighted boxes are shown below in more detail in Figure 10. Pixels with no data due to masking are gray. Figure A7. Histograms of the raster pixel distributions for EMSI across the entire time series for Acre Brazil with the reference period set over 34-years (eight years are absent due to lack of clear sky scenes). Years with deeper departures from normal are exemplified by redder or bluer means. Note, the increasing agricultural area extents result in an increase in negative values over time which are evidenced by a left tail. Figure A8. EMSI for 1999-2020 for Landsat 5, 7, and 8 time series for July in Yakutia, Russian Federation. The year 2000 had no clear sky scene. The two years with blue highlighted boxes are shown below in more detail in Figure 11. Pixels with no data due to masking are gray. Multiple repeated wildfires burn across the time series. In the EMSI trained on a global reference period (1999-2020) the areas which burn later in the scene are exceptionally high EMSI before the fire, and then disturbed after.