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

Abstract

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 significance (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) wildfire 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 findable, accessible, interoperable, reusable (FAIR) data principles.

1. Introduction

Human land-use and the rapidly changing climate are creating widespread plant water stress at global scale [1]. Consequently, monitoring the viability of Earth’s ecosystems is a key objective of the global Earth Observation Systems (EOS) network [2,3]. Satellite remote sensing technology enables the assessment and measurement of variation in ecosystem health [4,5,6,7] with increasingly fine spatial resolution [8,9]. Ecosystem applications of EOS include measuring and monitoring of the global carbon cycle [3,10,11], wildfire severity [12,13], insect and disease outbreaks [14,15,16], and severity of drought-induced forest morbidity and mortality [17,18,19,20].

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 real-time, 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.

1.1. 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].

1.2. 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 co-hosted 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] FMSI and Peters et al. [42] SVI, and the EMSI z-score calculations are conceptually identical (defined in Methods, Section 2.3). Meroni et al. [43] utilized standardized indices of SVI for vegetation anomaly monitoring, extending the Peters et al. drought framework with Moderate Resolution Imaging Spectroradiometer (MODIS) data.

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.

1.3. 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.

2. Materials and Methods

2.1. 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. Available Earth Observation Systems (EOS) platform radiometric bands (in nanometers (nm)) used to calculate Normalized Difference Vegetation Index (NDVI), spatial resolution (meters (m²)), and temporal resolution repeat [days], and dates of availability.

Satellite 1	Red (nm) (Band)	NIR (nm) (Band)	Spatial (m2)	Temporal (day−1)	Dates (month/years)
Landsat 1 MSS	600–700 (B5)	700–800 (B6)	60	16	07/1972–01/1978
Landsat 2 MSS	600–700 (B5)	700–800 (B6)	60	16	01/1975–02/1982
Landsat 3 MSS	600–700 (B5)	700–800 (B6)	60	16	03/1978–03/1983
Landsat 4 TM	630–690 (B3)	770–900 (B4)	30	16	08/1982–12/1993
Landsat 5 TM	630–690 (B3)	770–900 (B4)	30	16	04/1984–11/2011
Landsat 7 ETM+	631–692 (B3)	772–898 (B4)	30	16	06/1999–Present
Landsat 8 OLI	636–673 (B4)	851–879 (B2)	30	16	04/2013–Present
MODIS	620–670 (1)	841–876 (B2)	250/500	8/16	02/2000–Present
AVHRR	580–680 (B1)	725–1100 (B2)	1090	1	06/1981–Present
VIIRS	600–680 (B I1)	846–885 (B I2)	375/500	1	01/2012–Present
PlanetScope 2	590–670 (B3)	780–860 (B4)	3	1	07/2014–Present
Sentinel-2	645–684 (B4)	763–907 (B8)	10	5	07/2015–Present

¹MSS—Multispectral Scanner System, TM—Thematic Mapper, ETM—Enhanced Thematic Mapper, OLI—Operational Land Imager, VIIRS—Visible Infrared Imaging Radiometer Suite. ² Commercial data require a license and be hosted in private Google Earth Engine (GEE) collections.

). 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], WorldClim2 [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.

2.2. 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:

NDVI = (ρNIR − ρRED)/(ρNIR + ρRED),

(1)

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.

2.3. 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:

EMSI_xyt = (ρ x_xyt − µ_xyT)/σ_xyT,

(2)

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:

z_xyt = NDVI_xyt − µ_xyT/σ_xyT,

(3)

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. Proposed EMSI ranking system for judging ecosystem health based on local standardized scoring.

Index	Ecological State	Characteristics
z ≤ −2.0	In Collapse, Disturbed, or Misclassified	Spatial extents that have recently undergone a major disturbance, such as land clearing, overgrazing, or wildfire, exceeding 95% of all historical observations used in the reference
−2.0 < z < −1.0	Degrading or highly stressed	Spatial extents exhibit significant departure from normal conditions, likely due to extreme moisture stress, severe drought, or herbivory.
−1.0 ≤ z ≤ 1.0	Stable or Normal conditions	Variation falls within one standard deviation, meaning conditions are most common for that date or period.
1.0 < z < 2.0	Improving or Robust	Spatial extent is undergoing a surplus or pluvial-like condition with enhanced greenness or new growth.
2.0 ≤ z	Exceptional or Type Conversion	Exceeds 95% of all historic observations in the reference period. May indicate recovery from collapsed state, or conversion to an altered use.

).

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].

Statistical analyses of pixel level extraction were done with the base R [83] as well as: boot [84], dlm [85], devtools [86], ggfortify [87], ggplot2 [88,89], gridExtra [90], foreign [91], knitr [92], leafletR [93], lubridate [94], PerformanceAnalytics [95], plotly [96], raster [97], rgdal [98], rmarkdown [99], sp [100], and truncreg [98]. The complete linux package dependency list along with Rmd notebooks are hosted on Github and CyVerse.

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.

2.4. 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. Links to JavaScript (.js) scripts on the GEE for exporting pixel level NDVI extractions as time series, the surface reflectance or Top of Atmosphere (TOA) data products sampled, their processing algorithms for deriving surface reflectance, and supporting references.

.js Explorer	Reflectance Product SR/TOA/Raw	Processing Algorithm	Reference
Landsat 4	SR	USGS LEDAPS	[65]
Landsat 5	SR
Landsat 7	SR
Landsat 8	SR	USGS LaSRC	[66]
Sentinel-2	TOA	L1-C
MODIS	SR	MOD13Q1 v6
AVHRR	SR	NOAA CDR	[54]
VIIRS	SR

) 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 and Appendix B 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 Appendix A and Appendix 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.

2.5. Study Areas

We selected four model ecosystems (Figure 4,

Table 4. Locations mean annual temperature (MAT, ° Celsius), mean annual rainfall (MAT, millimeters per year), Koppen Climate Classification [103], and vegetation type descriptions.

Site	Long., Lat. [Decimal Degree°]	MAT [C°]	MAP [mm/yr−1]	Koppen Climate	Vegetation/Land Use
Savannah-Delta, Botswana	22.5601° E, −19.0901° S	20°	450	Bsh	Grassland, thornscrub, swamp. Natural, agriculture, rangelands
Taiga Yakutia	120.3300° E, 63.4300° N	−14°	240	Dfd	Conifer forest, Alaas meadow, pastoral, riverNatural, wildfires
Desert Grassland USA-Mexico	−110.5901° W, 31.3400° N	19°	330	Csa	Desert grassland and oak woodlands Natural, rangelands, agriculture
Rainforest Brazil	−71.9610° W, −9.1901° S	25°	2270	Af	Rainforest & Agriculture Natural, farming, logging

) 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].

3. Results

Figures 5–7 are presented in the text for the Okavango in Section 3.2 of the main text, and Figures 8 and 9 for the US-Mexico border region in Section 3.3 of the main text, these same figures are replicated for all other example sample locations can be (re)generated using the .Rmd notebooks that are supported in Appendix B. Time series for the reference periods of each model ecosystem location are presented in the Appendix A figure series (Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7 and Figure A8).

3.1. 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 (Figure 5 and Figure 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.

3.2. 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 (Figure 5 and Figure 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 (Figure 5b and Figure 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 (Figure 6 and Figure 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,b). Grass fires are captured in many of the scenes, particularly in the south eastern portion of the delta during the years 1992, 1995, 2005, 2006, and 2019 (Figure A4).

3.3. 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 (Figure 8 and Figure 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.

3.4. 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 (Figure A1 and Figure 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, 2001, 2005, and again in 2010, Brazil underwent severe droughts [104,105] all of these are captured in the EMSI (Figure 10 and Figure 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).

3.5. 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 cloud-free 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.

4. 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).

4.1. 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.

4.2. 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 (Figure 5b and Figure 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 Figure 10, Figure A6 and Figure 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).

4.3. 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 Appendix A and Appendix B).

4.4. 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 site-specific 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.

5. 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 and Appendix B). 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 dis-turbances 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.