Monitoring Forest Resilience: A Methodological Framework for Generating Cross-Calibrated Long-Term Landsat Datasets

1. Introduction

Global forests are increasingly threatened by climate and land-use changes across biomes (McDowell et al., 2020) that progressively erode their capacity to recover after disturbances (Li et al., 2018; Pimm, 1984). This loss of resilience manifests as critical slowing down (CSD), a phenomenon characterized by slower recovery rates, often measured through statistical indicators (Dakos et al., 2008). CSD-based metrics like variance, autocorrelation and skewness are expected to increase in proximity to a critical transition, or regime shift (Bathiany et al., 2025; Scheffer et al., 2015, 2009). These indicators, known as early warning signals (EWS), were originally derived from theoretical models of spatially homogeneous systems (Scheffer and Carpenter, 2003). Nevertheless, spatial and temporal EWS have been widely applied in climate, ecological, health, physical and social sciences for the determination of tipping point behaviour (Dakos et al., 2024); although they can also precede non-catastrophic shifts (Boettiger et al., 2013; Kéfi et al., 2013).

Assessments of forest resilience have primarily used field (Cailleret et al., 2019; Camarero et al., 2015; Colangelo et al., 2017), paleoecological (Stegner et al., 2019) and EO data (Carnicer et al., 2019) to identify CSD in timeseries. However, remote sensing has become crucial for monitoring and managing forest ecosystems (Lechner et al., 2020), given their vast spatial extent and inaccessibility. Despite covering relatively short time frames, satellite EO data are particularly important for understanding and detecting EWS across terrestrial ecosystems (Swingedouw et al., 2020), and have been used in dryland (Uden et al., 2019; Veldhuis et al., 2022; Wu et al., 2023), grassland (Hu et al., 2018), peatland (Alibakhshi et al., 2017; Lees et al., 2021a, 2021b) and forest (Boulton et al., 2022; Dwomoh and Wimberly, 2017; Liu et al., 2019; Majumder et al., 2019; Saatchi et al., 2021) biomes. This is because satellite EO datasets enable identification of biophysical proxies, such as vegetation indices (VIs), which are commonly used as state variables to quantify changes in vegetation resilience (Jiang et al., 2021; Lv et al., 2025; Rogers et al., 2018; Wang et al., 2025; Zhao et al., 2026). To date, many empirical EO assessments of forest resilience rely on Moderate Resolution Imaging Spectroradiometer (MODIS) and Advanced Very-High-Resolution Radiometer (AVHRR) satellite-derived products to detect EWS from timeseries of such variables at landscape scales (Error! Reference source not found., Error! Reference source not found.).

Delivering dense, continuous, and well-calibrated data for global vegetation monitoring, MODIS provides near-daily observations at 250-m to 1-km resolution from 2000 onward (Justice et al., 2002; Masuoka et al., 1998). Similarly, AVHRR has provided daily timeseries since 1978 at approximately 1-km resolution (Cracknell, 2024), commonly aggregated into global products (e.g., GIMMS at 8-km). Yet characterizing forest disturbance and recovery processes requires quantifying landscape patterns at ecologically meaningful scales (Levin, 1992; Turner et al., 1989; Wiens, 1989), meaning MODIS’ and AVHRR’s native spatial resolutions may be too coarse to resolve these dynamics (Cohen and Goward, 2004; Saura, 2004). Likewise, the spatial aggregation and multi-sensor merging underlying the Harmonized Landsat-Sentinel-2 (HLS) (Claverie et al., 2018; Ju et al., 2025) and the VODCA vegetation optical depth (VOD) datasets (Moesinger et al., 2020; Zotta et al., 2024) impact higher-order statistics with spurious variability (Dakos et al., 2012b; Smith and Boers, 2023). This is known to affect the reliable estimation of CSD-based indicators due to signal contamination (Smith et al., 2023; Tao et al., 2023), discrepancies between the spatial scale and species-specific trait responses (Liu et al., 2019) or mistakenly reflecting exogenous drivers (Miramontes, 2026).

In addition, monitoring forest resilience requires high resolution and long-term EO records capable of distinguishing gradual forcing-driven change against background variability (Williamson et al., 2016). This is particularly important for forest ecosystems with slower internal dynamics, where the decadal lifespan of trees and stand development result in longer recovery times following disturbances (Seidl and Turner, 2022) and lagged responses to environmental change (Seidl and Turner, 2022). Therefore, optical satellite-based approaches involve confronting a set of methodological constraints that are essential for evaluating whether potential critical transitions across such spatially heterogeneous landscapes can be reliably detected.

The Landsat archive provides 30-m spatial resolution with a nominal 16-day revisit interval dating back to 1972, offering one of the longest and most consistent records for assessing long-term forest dynamics (Wulder et al., 2019). However, only a quarter of the studies reviewed here employed Landsat data, and fewer exploited its full archive length (Error! Reference source not found.). This is mainly due to sensor inconsistencies, missing data and low temporal resolution (Bathiany et al., 2025; Lenton et al., 2024; Runge et al., 2025), along with anomalous features, or artifacts, that arise from sensor, acquisition, and processing errors rather than true surface conditions (Yang et al., 2026). Nevertheless, Landsat higher spatial resolution can capture processes at the finer stand and patch scales of forests (Banskota et al., 2014; Cohen and Goward, 2004; Townshend and Justice, 1988), including subtle trends in vegetation productivity (Rogers et al., 2018), abrupt disturbance events and slow recovery trends (Kennedy et al., 2010). It is theoretically ideal for estimating resilience based on spatial statistics at meaningful ecological scales (Lenton et al., 2024), thereby enabling additional mechanistic understandings of spatial patterns in ecosystems that traditional measures often miss by not accounting for spatial structure (Bastiaansen et al., 2020; Lennon, 2000). A Landsat-based spatiotemporal approach could then be directed to monitoring the resilience of forest biomes, as is more commonly done for non-forested terrestrial ecosystems (Error! Reference source not found.), enriching the insights gained from temporal EWS. However, no established framework has been developed for processing Landsat records tailored to the detection of critical transitions in forests.

We address this gap by outlining a reproducible end-to-end methodological framework that aims to mitigate current limitations to generate high-quality, long-term Landsat datasets for assessing forest resilience. This could facilitate the use of Landsat imagery in ecological applications of EWS, promoting its broader use for consistent spatiotemporal assessments of ecosystem resilience and support efforts to identify regions where the increasing probability of developing instabilities may lead to regime shifts. We first describe the main aspects of the proposed methodology, then demonstrate its application through a case study within a largely undisturbed area of the boreal forest in the province of Quebec in Eastern Canada.

2. Methods

The proposed methodological framework (Error! Reference source not found.) offers a robust processing workflow to tackle practical and multispectral limitations of Landsat imagery, and generate high-quality, analysis-ready composites suitable for the computation of EWS across large extents of forest ecosystems. Based on state-of-the-art tools such as the LandsatTS pipeline for automated cross-sensor harmonization, a multidimensional compositing algorithm implemented via cloud computing, and a scalable hexagonal discrete global grid system, our framework specifically addresses three major methodological limitations: i) CSD-based indicators requiring long observational records. This is addressed by maximising the length of the Landsat timeseries through biome-specific cross-calibration coefficients and a season specific mosaicking style. (ii) Sensor artifacts. These are reduced through robust imagery processing and compositing tailored to higher order statistics. (iii) The mismatch between the scale of observation and ecological processes. This is addressed via a multiscale hexagonal grid system to identify optimal spatial resolution(s) to measure an ecosystem. The same grid was then used to support the efficient computation of spatiotemporal EWS (STEWS) from higher spatial resolution imagery at multiple scales over extensive geographic domains. The JavaScript and RStudio scripts for this case study analysis are available on GitHub https://github.com/camuradr/LandsatCompositesGEE and can be modified with custom input parameters.

Figure 2. Methodological scheme for the generation of long-term high-quality Landsat composites and the computation of statistical indicators of resilience. The workflow and general criteria to address practical (i) and multispectral (ii) limitations (Panel A), along the multiscale grid approach to tackle the scale (iii) limitations (Panel B) are delineated. The corresponding examples from the case study are displayed as well. Nobs = number of observations per pixel. ROI = region of interest.

Figure 2. Methodological scheme for the generation of long-term high-quality Landsat composites and the computation of statistical indicators of resilience. The workflow and general criteria to address practical (i) and multispectral (ii) limitations (Panel A), along the multiscale grid approach to tackle the scale (iii) limitations (Panel B) are delineated. The corresponding examples from the case study are displayed as well. Nobs = number of observations per pixel. ROI = region of interest.

2.1. Cross-Calibration for Long-Term Records

The Landsat collections available in Google Earth Engine (GEE) provide high-resolution satellite imagery with global coverage for over four decades since 1984 (Appendix S2). Nonetheless, systematic differences arising from slight variations in wavelengths between Landsat sensors can significantly impact multidecadal timeseries interpretation or subtle ecological change detection analyses (Hemati et al., 2021). The most recent version of Landsat products, Collection 2, has undergone important corrections by improving data interoperability between sensors (Crawford et al., 2023); however, differences from sensor specific wavelength ranges between the spectral bands remain. Consequently, sensor cross-calibration is still recommended to avoid the introduction of false trends (Berner et al., 2020; Chen et al., 2024; Trevisiol et al., 2024), from known minor reflectance discrepancies influencing multi-sensor data compatibility (Ju and Masek, 2016).

Calibration coefficients developed for past collections (Roy et al., 2016a), and linear regressions (Wu et al., 2022) were previously used to correct sensor bias and ensure radiometric consistency between sensor observations (Appendix S3Error! Reference source not found.). Recent collections, however, require up-to-date coefficients to ensure their ongoing compatibility. A recently developed algorithm to cross-calibrate spectral bands among Landsat sensors was implemented using the LandsatTS package in RStudio (Berner et al., 2023). Rather than normalising to a reference scene (Wang et al., 2025), it employs a machine learning approach that uses representative pixel samples across a region to standardise data to Landsat 7 ETM+, due to its temporal overlap with Landsat 5 TM and Landsat 8 OLI (Error! Reference source not found.). A polynomial regression model then computes calibration coefficients derived from pixels within the study region, suggested to be good practice for change detection analysis (Richardson et al., 2025). The coefficients were applied to their corresponding Landsat collections in GEE, adjusting each surface reflectance (SR) band values to the fitted polynomial relationship (Error! Reference source not found. and Error! Reference source not found.). This bandpass adjustment offers a Landsat dataset specifically calibrated for an ecoregion or biome, minimizing the limitations of merging records from various sensors. Likewise, VIs derived from calibrated SR bands can be used as spatially explicit proxies of vegetation condition, including ecological responses to environmental change (Pettorelli et al., 2005), land monitoring and mapping (Wulder et al., 2008) at broad spatial extents.

2.2. Compositing for Resilience Indicators

The framework for generating high-quality, cloud-free images of the landscape was adapted from Li et al., (2019), implementing a pixel-based composite paradigm that selects the most representative reflectance value across all available Landsat scenes over the compositing period. First, imagery with less than 70% cloud cover was retained to avoid geometric correction errors, when reference points used to align each image to the ground are obscured by clouds (Goward et al., 2022). Second, the Quality Assessment (QA) bands were used to mask unreliable pixels, including those saturated, filled, or affected by clouds, cloud shadows, snow or ice features. Persistent residual contamination where QA masking omitted cloud or haze pixels (Hughes and Kennedy, 2019) was removed using post-compositing thresholds for unusually bright pixels (Sun et al., 2016). Water pixels were also masked, using the JRC global surface water dataset (Pekel et al., 2016), since they alter spatial metrics by decreasing spatial autocorrelation of surrounding land pixels due to their contrasting spatial patterns (Griffith and Chun, 2016). Additionally, edge cropping was applied to remove artifacts within ~0.5 km margins of the Landsat scenes, mostly associated with TM and ETM+ sensors (Robinson et al., 2017).

The preservation of spectral covariance structure across bands is also essential for the reliable computation of indicators (Smith et al., 2023; Smith and Boers, 2023). While algorithm performance is context-dependent with no one-size-fits-all method (Qiu et al., 2023), the geometric median (GeoMedian) based compositing has been suggested to provide a robust representation of a given period, even when the number of clear observations over large spatial extents is low (Roberts et al., 2017). In fact, annual high-quality GeoMedian composites, such as the GeoMAD Digital Earth Africa and Digital Earth Australia products, are highly suitable for spatial EWS analysis (Roberts et al., 2018) and classification (Wellington and Renzullo, 2021). Unlike pixel-selection compositing methods which select a single real observation such as the medoid and best available pixel (BAP), the GeoMedian computes a representative synthetic pixel by considering each clear observation as a point in a multi-band reflectance space and selecting the multivariate centre that minimises the summed Euclidean distance to all observations (Flood, 2013; Roberts et al., 2017). This avoids the mathematical biases introduced by temporal aggregation (Liu et al., 2026), reduces spatial noise, improves radiometric consistency, and better preserves the spectral relationships among bands (Francini et al., 2023).

By contrast, selecting single observations produces spatially discontinuous outputs (Meng et al., 2023), especially when neighbouring pixels are drawn from different acquisition dates (Roberts et al., 2017). For example, the maximum value composite (MVC) of the normalised difference vegetation index (NDVI) is a widely used approach for capturing peak vegetation growth, as it selects the greenest pixel values (Holben, 1986), which better discriminates against snow and cloud cover than mean and median compositing (Qiu et al., 2023). However, residual cloud shadows can mistakenly generate anomalously high NDVI values (Sousa et al., 2003). This is further compounded by an inherent bias towards phenological extremes, producing less radiometrically consistent image mosaics (Van doninck and Tuomisto, 2017) and reducing the contrast between vegetation types such as trees and grasslands (Flood, 2013). Moreover, MVC is sensitive to directional reflectance effects which, even though smaller for Landsat compared to wide-swath sensors (e.g. MODIS), remain non-negligible at high latitudes (Roy et al., 2016b). These effects are mitigated by the GeoMedian (Roberts et al., 2017), providing a particularly well-suited approach for peak of season (POS) compositing for temperate and tropical forest landscapes.

In addition, since conventional temporal EWS indicators are sensitive to residual non-stationarity and noise structure, rather than to the absolute magnitude of the underlying signal, compositing approaches that retain greater artifact-driven variability (e.g. MVC), risk generating spurious trends independent of genuine loss of resilience (Boettner and Boers, 2022; Dakos et al., 2012a). This is expected to be more pronounced for spatial EWS, rendering artifact removal critical for the utilisation of VIs for monitoring resilience (Zeng et al., 2022). Conversely, limitations of temporal aggregation are avoided by spatial EWS (Kéfi et al., 2014; MacLaren et al., 2025), because they do not necessarily require high-frequency timeseries to detect CSD trends (Butitta et al., 2017). The aggregation window length is ecosystem-specific and defined as the POS, when the dominant vegetation reaches peak phenology (Jeganathan et al., 2014). This is the summer months in boreal and temperate forests (Jarvis and Leverenz, 1983; Myneni et al., 1998) and the period of maximum solar radiation in tropical forests (Huete et al., 2006; Myneni et al., 2007).

2.3. Hexagonal Grids for Multiscale Indicators

To ensure that indicators computed from biophysical proxies are comparable across regions and scales of analysis, a reproducible approach for partitioning the landscape into consistent spatial units is necessary to reduce their sensitivity to arbitrary zoning and resolution choices (Jelinski and Wu, 1996). Squared grids are frequently used because of their congruency and alignment for grid partitioning (Birch et al., 2007). However, they become significantly affected by distortions when assessing large extents of the Earth (Sahr et al., 2003), particularly at higher latitudes (Fletcher, 2023). Avoiding this is crucial when assessing STEWS over extensive areas across different biomes.

Unlike conventional squared grids (Feng et al., 2024), Hexagonal Discrete Global Grids (HDGGs) result in compact equal-area tessellations that reduce directionality bias, provide uniform adjacency, and preserve equal area on the sphere (Kimerling et al., 1999a; Uher et al., 2019). Their applications have proven suitable for constructing globally consistent spatial units that minimize geometric and statistical artifacts (Sahr et al., 2003). Existing HDGG configurations provide balanced compromises between equal-area features, near-isotropic cell geometry, and a hierarchical index system (Kimerling et al., 1999b; Kmoch et al., 2025), which is well suited for environmental analyses at various scales (Birch et al., 2007). These features are advantageous when dealing with substantial geospatial datasets through cloud computing and parallel processing, supporting statistically robust spatial aggregation and map-based operations (Li and Stefanakis, 2020), while improving efficiency of STEWS computation (Appendix S3).

Multiple HDGGs were generated using the DGGRID open-source Python implementation (Sahr, 2025). Compared to other open-source tools (Kmoch et al., 2022a, 2022b), DGGRID provided a reproducible workflow which facilitated the construction of tessellations that can be seamlessly incorporated into standard GIS software as vector layers (Barnes, 2018). We propose HDGGs as a multiscale approach for unbiased cross-biome comparisons, expected to help in the identification of appropriate observation scales for system-specific ecological process, indicators or proxies.

2.4. The Boreal Forest: A Case Study in Eastern Canada

Boreal forests are currently experiencing some of the fastest and most significant warming on the planet (IPCC, 2022), and have been reported as a tipping element of the Earth’s biosphere (Armstrong McKay et al., 2022; Lenton et al., 2008). Despite diverging findings on boreal forest resilience, with reported declines (Smith et al., 2022) and increases (Forzieri et al., 2022; Southworth et al., 2023), evidence of an undergoing southern dieback and a northern expansion (Rotbarth et al., 2023) suggests it is potentially transitioning from the current bimodal distribution of tree cover, consisting mostly of either dense forests or sparse woodlands, towards a unimodal semi-open state with intermediate tree cover (Rotbarth et al., 2025). The observed trends of greening and densification in the northern interior forests (Berner and Goetz, 2022; Sulla-Menashe et al., 2018; Zhang et al., 2024) are consistent with the patterns of boreal forest change.

Although large disturbance footprints are captured by coarser spatial resolution imagery, small disturbances often remain under-detected (Bucha and Stibig, 2008; Jin and Sader, 2005; Potapov et al., 2008). Considering the broad spatiotemporal scales at which major disturbances typically trigger changes in the vegetation structure of boreal forests (Holling, 1992), finer-resolution datasets are required to resolve landscape changes (Wang and Friedl, 2019). The relatively high resolution of Landsat captures the stand and patch scales at which boreal disturbance and recovery processes occur (Frazier et al., 2015; White, 2024; White et al., 2017), closer to the key ecological scales that better capture boreal forest stability and change (Berner et al., 2020). In fact, Landsat has been able to characterize forest dynamics related to disturbances and climate-induced growth (Wang et al., 2020), including insect defoliation (Rahimzadeh-Bajgiran et al., 2018; Townsend et al., 2012), wildfires and timber harvest (Schroeder et al., 2011). These inherently spatial phenomena (Seidl et al., 2023) are often missed or aggregated away in coarse resolution sensors, such as MODIS or AVHRR, and resilience monitoring products (Zhang and Liu, 2024). In addition, Landsat’s native compatibility with the National Forest Inventory (NFI) and National Terrestrial Ecosystem Monitoring System (NTEMS) datasets, developed to annually quantify and characterise forest change (Hermosilla et al., 2016; Wulder et al., 2012), provides a further advantage of using Landsat to quantify forest resilience. The HDGG resolution levels considered in the tessellation analysis were then ecologically linked to the hierarchical structure of the boreal forest (Holling, 1992), comprising the stand, community and landscape scales across five different configurations (Error! Reference source not found.8).

To facilitate robust testing of our framework, we chose a large, undisturbed region of the boreal forest in northern Quebec as our case study (Error! Reference source not found.), based on the intact forest landscapes (IFL) dataset (Potapov et al., 2017). This selection minimises the influence of anthropogenic disturbances, allowing us to focus on natural disturbances and thereby simplifying the evaluation of our framework. The case study area is sufficiently large to accommodate the broad range of HDGG extents and configurations. Covering 325,514.11 km², approximately 2.5 times the size of England, it enables the demonstration of the utility of our proposed scale-aware study design (Chave, 2013). It is worth noting that while IFL areas have no evidence of anthropogenic disturbances, natural disturbances such as wildfires still occur within them. Nevertheless, compared to similar areas in western Canada, the Eastern Canadian shield taiga remained relatively undisturbed until the wildfires of unprecedented scale and severity in 2023 (Byrne et al., 2024). Following the spatial extent delimitation, steps 1 to 3 of the framework (Error! Reference source not found., Panel A) were implemented.

Figure 3. Spatial extent of the region of interest (ROI in red) corresponding to an undisturbed boreal forest area in Quebec, Eastern Canada.

Figure 3. Spatial extent of the region of interest (ROI in red) corresponding to an undisturbed boreal forest area in Quebec, Eastern Canada.

The quality of the composites was evaluated for data gaps, noise and radiometric consistency according to the criteria for pixel-based Landsat composite assessment delineated by Francini et al., (2023). The revised composites were then exported and mosaicked in RStudio to build VIs raster stacks from which indicators (i.e. spatial variance, autocorrelation and skewness) at five different HDGGs scales were computed, following steps 5 and 6 (Error! Reference source not found., Panel B). The processing parameters chosen specifically for the case study ecoregion are depicted in Error! Reference source not found..

3. Results

3.1. Reducing Multispectral Limitations

Spatially and radiometrically consistent cloud-free composites were generated for all SR bands. For the chosen region in Quebec, a two-year aggregation window was determined to be adequate to achieve near-complete land coverage (≥ 96%) for all composites, significantly reducing the spatial gaps (Error! Reference source not found.). This was determined based on the total number of observations available for annual composites being notably less when compared to biennial composites, especially for the periods when Landsat 5 provided the only available imagery (Error! Reference source not found., Panel A). The incorporation of Landsat 7 ETM+ was also required to achieve constant total (Error! Reference source not found., Panel B) and per month (Error! Reference source not found.Error! Reference source not found.) numbers of observations between composites, as the additional observations it provided significantly reduced the spatial gaps across the composites.

The resulting biennial composites were consequently distributed into fourteen composites between 1984 and 2011 and six composites between 2013 and 2024, constituting a dataset with a total of 20 composites. It is worth noting that a one-year gap corresponding to the summer of 2012 remained, since there is no dataset overlap between the Landsat 5 TM and Landsat 8 OLI collections, negatively impacting the number of observations available. Nonetheless, the incorporation of the available imagery from 2012 was tested to determine whether this had any significant impact on the composites and trends observed, with no significant change on the overall trend observed (Error! Reference source not found.). Incorporating the Landsat 7 ETM+ data to specifically account for 2012 in the timeseries was then deemed unnecessary, as the available observations for the selected years were already satisfactory and its inclusion would result in inconsistent aggregation windows.

Figure 4. Valid number of observations available for (A) land coverage between one- and two-year composites, and (B) for the total area when applying a one- and two-year aggregation window. The temporal coverage of each Landsat sensor is displayed at the bottom of panel (B). The shaded grey area indicates the gap year (2012).

Figure 4. Valid number of observations available for (A) land coverage between one- and two-year composites, and (B) for the total area when applying a one- and two-year aggregation window. The temporal coverage of each Landsat sensor is displayed at the bottom of panel (B). The shaded grey area indicates the gap year (2012).

Both the MVC and GeoMedian approaches were compared to provide insights into the suitability of incorporating the latter in the proposed scheme. The GeoMedian was particularly effective at suppressing artifacts commonly found in optical satellite imagery (Figure 1). It successfully reduced the presence of ETM+ stripes (SCL-off), along with eliminating residual haze and shadow noise. Residual cloud contamination was observed for the GeoMedian composites in cases with fewer available observations (Figure 1, Panel E). These were removed using a saturation threshold derived from the mean reflectance value across all bands, selected empirically as 0.2 from visual inspection for this study area. Altogether, the spatial and radiometric consistency of the composites was improved, supporting the reliable computation of STEWS from Landsat imagery and the use of the GeoMedian approach for doing so.

Figure 1. Landsat composite comparison between the (A) MVC NDVI and (B) GeoMedian (GM) approaches, displaying improved artifact removal for spatial consistency for the latter. Samples from the biennial composite show differences in atmospheric haze (C), residual clouds (C), SCL-off stripes (E) and residual shadows (F). Images in false colour infrared (SR bands 5/4/3) band combination without water mask to highlight vegetation cover (in red hues), contrasting water bodies (dark blue or black).

Figure 1. Landsat composite comparison between the (A) MVC NDVI and (B) GeoMedian (GM) approaches, displaying improved artifact removal for spatial consistency for the latter. Samples from the biennial composite show differences in atmospheric haze (C), residual clouds (C), SCL-off stripes (E) and residual shadows (F). Images in false colour infrared (SR bands 5/4/3) band combination without water mask to highlight vegetation cover (in red hues), contrasting water bodies (dark blue or black).

3.2. Merging the Landsat Collections

Region-specific calibration coefficients were derived from representative pixel samples across the province of Quebec (Error! Reference source not found.A). The obtained coefficients were similar to the ones provided by Berner et al. (2023) from a high-latitude dataset of randomly sampled terrestrial locations within the Arctic tundra and boreal forest biomes globally (Error! Reference source not found.). Small variations in coefficient values arise partly because the collections of pixel samples extracted from the imagery differ between calibrations. In addition, the calibration functions have some degree of random sub-setting of those samples, which causes the model coefficients to vary even if run multiple times on the same dataset. This resulted in slightly different NDVI timeseries (Table S7), but with the overall trends remaining quite consistent across networks of sites for both the MVC (Error! Reference source not found.) and the GM (Error! Reference source not found.). As expected for coefficients tailored to the same biome, the results were similar when applying them to the Landsat imagery as both sample sizes were large.

Temporal trends in mean NDVI from 1984 to 2024 evidenced the greening of the boreal forest biome of Quebec (Figure 6). The NDVI was taken as an exemplar to evaluate the composite timeseries as it is primarily associated with the reported greening trends. An analysis of the greening trend in NDVI illustrates the value of the compositing and calibration pipeline we propose, as it adjusts the magnitude and rate of the overall NDVI trend. Even though an increasing trend in NDVI was consistently observed, the rate of change varies significantly. Uncalibrated versions of the Landsat dataset were particularly different from calibrated ones, inflating the trends (0.0035 ± 0.0007 NDVI/yr), compared with calibrated datasets (0.0017 ± 0.0008 NDVI/yr). This reflects the systematic offset introduced by uncalibrated sensor transitions, which is observed in the abrupt discontinuity of the NDVI increase when Landsat sensors shift (Figure 6), particularly after 2012, when Landsat 5 imagery was no longer available, and Landsat 8 was introduced. The influence of sensor-shift induced leaps in the timeseries was further reduced by the incorporation of Landsat 7, although trends remained comparable within both the calibrated and uncalibrated versions.

Similarly, the compositing strategy impacted the magnitude of the trend, where MVC trends were steeper compared to the GeoMedian calibrated trends, inflated by 40.2% for the calibrated MVC (Figure S6), and up to 130.7% for the uncalibrated MVC. Cross-calibrating was clearly beneficial to reduce the influence of sensor shift leaps on the timeseries. This highlights the importance of carefully distinguishing an increase in NDVI as evidence of canopy densification, a shift toward broadleaf dominance, or recovery from disturbance. The GeoMedian compositing strategy, along with the regional cross-calibration coefficients, provided a higher spatial and temporal consistency across the composites for a more reliable computation of STEWS from Landsat imagery.

3.3. Beyond Pixels with Multiscale Grids

The HDGGs enabled the re-sampling of the Landsat data into ecologically meaningful units. This allowed for STEWS to be calculated at multiple scales of observation. The effect of scale of observation and resolution could then be assessed using varying grid configurations and resolutions. Integrating a multiscale grid system to analyse biomes at relevant, ecologically grounded scales is crucial for ecological applications of EWS. While processing higher resolution imagery is computationally expensive, our HDGG approach also addresses this via parallelisation. It greatly improved the computation times of indicators from raster data (Error! Reference source not found.Error! Reference source not found.), reducing overall processing times by 69.6% compared to sequential attempts.

4. Discussion

Our workflow provided consistent high-quality Landsat composites for a large region of the boreal forest in eastern Canada utilising the full Landsat Collection-2 archive. It therefore represents an important advance for existing studies that use Landsat-derived EWS to assess forest resilience. These have generally covered shorter time windows at local or regional scales (Liu et al., 2019; Vargas et al., 2019; Wu et al., 2022). Hui Wang et al. (2025) combined all three sensors (TM, ETM+ and OLI) over a ~16-year window to evaluate local vegetation resilience. Rogers et al. (2018) detected EWS of tree mortality across boreal North America over a ~28-year window but including two sensors (TM and ETM+). The employed Landsat-derived NDVI timeseries addressed multi-sensor differences with a scaling factor (Ju and Masek, 2016), primarily calibrating at the NDVI-level. However, index-level cross-calibration does not eliminate per-band reflectance discontinuities, which are known to introduce artificial long-term trends in combined NDVI records (Sulla-Menashe et al., 2016), to which statistical indicators are particularly sensitive (Dakos et al., 2012b; Lenton et al., 2012). Our study clearly show that the region-specific calibration coefficients allowed for a consistent multi-sensor Landsat record suitable for reliable resilience analysis. Global scale monitoring of the boreal forest could rely on existing, already-developed biome-specific coefficients. Therefore, we recommend deriving broader biome-level calibration coefficients for further terrestrial ecosystems (e.g. rainforests, temperate forests and drylands) and, if necessary, compare with regionally specific coefficients.

Our study used a relatively long aggregation window of two years, as for forest ecosystems, longer aggregation windows may be acceptable due to their higher stability and slower vegetation dynamics compared to grasslands and drylands (Han et al., 2025). However, the length of the aggregation window can be dynamically adjusted within our framework based on the local density of cloud-free observations, or the characteristic ecological recovery timescale of the specific biome under study. As such, shorter aggregation windows are possible for less extensive and/or less cloudy regions, allowing for annual composites to increase temporal resolution.

Through our case study it was also confirmed that without a built-in bias toward greenness, the GeoMedian obtains a physically representative state over the entire compositing period, capturing abrupt changes between compositing periods. While it smooths out rapid within-period events (Roberts et al., 2017), the GeoMedian excels at detecting cross-period long-term changes rather than transient fluctuations (Gupta et al., 2024), making it ideally suited for spatial EWS. By contrast, the MVC algorithm will actively avoid these observations since disturbances usually represent a state of low productivity (Flood, 2013). With at least one pixel observation from a pre-disturbance and/or a recovered post-disturbance state, disturbances within the aggregation window are effectively removed from the MVC composite. The associated timeseries would then only reflect the most productive state of the landscape rather than its actual average condition.

Furthermore, our framework can be applied for any VI. This is important, as the choice of VIs should be tailored to the target ecosystem or biome, as different indices vary in their sensitivity and applicability across vegetation types (Lenton et al., 2022; Zeng et al., 2022). We therefore recommend employing multiple VIs when drawing ecological interpretations to improve robustness and reduce the risk of biome-specific bias. For instance, in the boreal forest, an increase in broadleaf species relative to conifer-dominated stands following disturbance is expected to be reflected in the increase of vegetation productivity indices, such as the Enhanced Vegetation Index (EVI) and NDVI (Fiore et al., 2020; White et al., 2023). By contrast, other VIs are spatially correlated to natural disturbances such as wildfires and insect outbreaks (Mackey et al., 2024). The Normalized Burn Ratio (NBR), for example, is more sensitive to structural recovery of the canopy (Pickell et al., 2016). Similarly, the Shortwave Infrared Water Stress Index (SIWSI) captures water-stressed canopies (Shestakova et al., 2024), whereas the Normalized Difference Moisture Index (NDMI) decreases under drought stress and pre-fire fuel availability (Sivrikaya et al., 2024), distinguishing canopy water content from structural or compositional change. Indices derived from the near-infrared (NIR) and shortwave infrared (SWIR) bands (i.e. NBR, NDMI and SIWSI) tend to be more strongly correlated with changes in above-ground biomass than greenness-based indices are (Meimand et al., 2026). Therefore, to better identify processes driving a given trend (e.g. stand succession, biome-margin shift, shrubification, tree mortality), complementary indices serve specific disambiguation roles. Our framework is also flexible enough to enable future work which could incorporate biophysical or structural variables such as fraction of photosynthetically active radiation (fAPAR) or leaf area index (LAI), which are less confounded by ground surface reflectance, atmospheric conditions, and sensor characteristics (Myneni and Williams, 1994).

Future work should consider integrating Landsat 9, given the mission continuity it provides (Masek et al., 2020) and its robust instrument similarity to Landsat 8 (Gross et al., 2022), although sensor-specific cross-calibration or transformation coefficients (Zhang et al., 2025) remain necessary to ensure inter-sensor consistency. On the other hand, geospatial foundation models that generate satellite embeddings from multi-sensor imagery, such as AlphaEarth Foundations (Brown et al., 2025) and TESSERA (Feng et al., 2026), represent a different paradigm in which sensor harmonization is learned as part of model training rather than applied as a separate calibration step. Beyond this conceptual difference, both products are currently restricted to the Sentinel-era archive (2017 onwards) and cannot yet support the multi-decadal analysis required to detect ecological tipping points. Their ecological interpretability would also likely prove more challenging compared to VIs based on physical reflectance properties. AlphaEarth native availability in GEE could nonetheless make it useful as a complementary validation layer for recent periods, since it could integrate easily with existing repositories.

Consistent, cross-calibrated long-term archives such as Landsat remain essential for monitoring forest resilience at large scales, especially as a global decline in vegetation resilience since the early 2000s has been reported (Forzieri et al., 2022; Rocha, 2022; Seddon et al., 2016; Smith et al., 2022; Zhang et al., 2024). Tropical rainforests have shown slower recovery rates under drier conditions associated with mean annual precipitation declines (Verbesselt et al., 2016). This is particularly evident in the Amazon rainforest where loss of resilience driven by deforestation and climate change is increasing the likelihood of rainforest dieback (Boulton et al., 2022; Van Passel et al., 2024; Wang et al., 2024). Similarly, significant declines in the resilience of the northwestern American boreal forests have been detected in regions affected by forest losses and wildfires along the southern margins of warm, dry and densely vegetated areas (Rotbarth et al., 2025, 2023; Zhang et al., 2024). In contrast, other studies report opposite patterns, with increasing resilience in boreal forests accompanied by declining resilience in tropical rainforests (Forzieri et al., 2022; Southworth et al., 2023) or, conversely, increasing resilience in the tropics coupled with declining resilience at higher latitudes (Smith et al., 2022). While forest vegetation responses are highly variable when measured by season, biome and land use (Saatchi et al., 2021; Southworth et al., 2023), these divergent findings may be partly attributable to variable practices for measuring resilience in forest systems, along with the spatiotemporal scales at which forests are being monitored. Resolving these discrepancies requires consistent methodological frameworks, such as the one developed here for Landsat data, to support the application of EWS across forest ecoregions that reliably characterize the changing resilience of Earth tipping elements.

5. Conclusion

The utilisation of satellite EO for measuring ecosystem resilience through statistical indicators has been widely explored across terrestrial biomes. Long-term remote sensing data should be carefully treated when applied to measure loss of ecosystem resilience because current EO limitations for detecting EWS need to be considered. Reliable indicators of resilience often require (i) long-term timeseries, (ii) consistent spatial data and (iii) scale-aware approaches. We have presented a methodological framework that minimizes these limitations by leveraging the full length of the higher-resolution Landsat archive for resilience assessments. While no single compositing or calibration approach is unambiguously optimal for all applications, we outline the approach we consider best suited for STEWS computation, on the basis that it minimises the artifact-driven non-stationarity to which indicators are particularly sensitive. The regionally cross-calibrated GeoMedian composites provided a consistent multi-sensor Landsat record suitable for a more reliable forest resilience assessment, though further testing would be needed for non-boreal biomes. Critically, our framework’s multiscale HDGG approach provides a highly efficient way to process EO data through parallel processing that is scale-comparable because of the indexing system it offers, and which can be easily integrated with additional geospatial datasets, thereby, providing a consistent method for quantifying STEWS in any terrestrial ecosystem. In combination, application of our framework could shed light on the discrepancies among the trends observed from different indicators and approaches applied to ecosystems worldwide. This is central to understanding how forest biomes are changing, and to assessing whether they are approaching critical transitions.