Timing recovery of ecosystems in sequential remotely sensed and simulated data

doi:10.1038/protex.2017.038

Method Article

Timing recovery of ecosystems in sequential remotely sensed and simulated data

https://doi.org/10.1038/protex.2017.038

This work is licensed under a CC BY 4.0 License

This protocol has been posted on Protocol Exchange, an open repository of community-contributed protocols sponsored by Nature Portfolio. These protocols are posted directly on the Protocol Exchange by authors and are made freely available to the scientific community for use and comment.

Version 1

posted

You are reading this latest protocol version

The time needed for ecosystems to recover from a disturbance has been proposed as a generic indicator of ecosystem resilience. The lengthening of the recovery time with increasing stress is referred to as “Critical Slowing Down” and has been proposed as an early warning of a nearing tipping point. Hence, methodologies for measuring recovery rates and critical slowing down in remotely sensed data might provide a powerful way to synoptically asses ecosystem resilience. Here, we present a protocol using an algorithm to measure the recovery time after a disturbance from sequential spatial data. The algorithm can be applied to both empirical, e.g. remotely sensed, and simulated spatial data.

Computational biology and bioinformatics

Mathematics and computing

disturbance

ecosystem dynamics

recovery

resilience

stability

tipping points

Changes in environmental conditions can have important ramifications for the resilience of ecosystems and understanding how they respond to subsequent perturbations is an urgent need in the context of global change¹. Theoretical studies have proposed that the time needed to recover from a small perturbation can be used as an indicator for resilience and can signal an upcoming collapse of the ecosystem state^2,3. This phenomenon (Fig. 1) of an increasing recovery time T_r with decreasing resilience (i.e. increasing stress) is referred to as “critical slowing down” (CSD). Although the idea is simple, straightforward implementations of the concept in real biological systems remain scarce. Proof of principle has been provided by controlled laboratory experiments with cultures of algae⁴ and yeast⁵. Yet, direct observations of CSD in real ecosystems are to our knowledge absent from the literature. This might be related to the fact that the methods to test the concept of CSD in (spatial) datasets of real ecosystems need further development⁶. Hence, in this protocol we propose a methodology for measuring recovery times in sequential spatial data that can e.g. be applied to remotely sensed data sets. Therefore, this method can provide a powerful way to synoptically asses ecosystem resilience and test for CSD.

The approach

Here a methodology is proposed to use changes in presence and absence of the key habitat-forming species (e.g. tidal marshes vegetation, as in the associated ref. [7]) in sequential spatial data as a natural disturbance-recovery experiment. Disturbance of e.g. vegetation can be observed if its presence changes into absence on the subsequent snapshot (Fig. 2). The time needed for a grid cell (i.e. pixel) to recover once disturbed, can be measured as the recovery time T_r. Pixels that do not recover within the sequence are flagged as censored observations to prevent biases when estimating the recovery rate λ from the recovery time T_r (e.g. pixel 4 in Fig. 2). Pixels that are not disturbed are considered stable and not used in the analysis, even if it was not present from the first snapshot (e.g. pixel 5 in Fig. 2). Although this methodology uses discrete data (presence/absence) instead of continuous data, it is analogous to the proposed procedure for measuring the recovery time (Fig. 1). Finally, this approach results in mapping the recovery time in space on a pixel-by-pixel basis.

Next, the recovery rate λ can be reliably estimated from the observed recovery times T_r. When assuming an exponential trajectory towards the predisturbed state, the recovery rate λ can be estimated as 1/T_r.

The algorithm

The algorithm (Fig. 3) computes the recovery time on a pixel-by-pixel basis, i.e. disturbance and subsequent recovery in each grid cell is tracked separately. Its code is provided as a function implemented in Matlab (fTimeRecov.m), which can be found in Supplementary File 1 (see the Source code section in Equipment for more information). As input, the algorithm needs a stacked grid g of L layers. The layers are ordered sequentially and hold a snapshot of the ecosystem containing information about the presence and absence (i.e. binary data [0, 1]) of the key habitat forming species at a particular time step. In order to calculate the time between disturbance and recovery T_r, an additional vector t with the timestamps of the snapshots in each layer is required. Although a sequence of snapshots taken at regular time intervals is prefered, the algorithm can also handle sequences taken at irregular time intervals.

Starting from the oldest snapshot, the algorithm compares a snapshot with the subsequent snapshot to find the disturbances. Next, the remaining snapshots are scanned in ascending order for recovery to occur. If no recovery occurs within the sequence the grid cell is flagged as censored. This procedure is repeated for the next successive pair of snapshots until all (L-1) pairs in the sequence have been visited.

After running the algorithm, all observed recovery times (T_r), censored recovery grid cells (Cens) and all disturbances (Dist) are outputted as a stacked grid with L-1 layers. To calculate recovery rates the T_r and Cens outputs are needed. Additionally, the average recovery time is outputted as a single grid, which enable easy visualization of the results.

The protocol is developed specifically for disturbance-recovery processes in tidal marshes where vegetation is classified on the basis of Normalized Difference Vegetation Index (NDVI) and inundation by seawater is the main stressor. However, this methodology can be applied to any vegetated or non-vegetated ecosystem in which the key habitat-forming species can be identified. It can, furthermore, be applied to compute recovery times in the sequential output of spatial models.

Sequential spatial data:

Sequential spatial data is needed in which disturbances and recovery can be observed. The algorithm requires a stacked grid as input, in which each layer is a snapshot in time containing information on the presence or absence of key habitat forming species (e.g. vegetation) of the ecosystem of interest. This can be simulated (e.g. the output of a spatial model) or empirical data (e.g. a sequence of remotely sensed images).

Source code:

The algorithm (Fig. 3) is provided as function (fTimeRecov.m) implemented in Matlab. The function is provided in the ZIP package that can be found as Supplementary File 1. The ZIP file contains, next to the function, two worked-out examples to illustrate the use of the functions and the workflow. The first example (Example1.m) shows how the five pixels in the 7-by-7 grids in Fig. 2 are traced and analyzed. In the second example (Example2.m) the full workflow, as described in Procedure from 1.4 onward including calculating the recovery rates along a stress gradient, is illustrated (see also Fig. 4). This example (Fig. 4) is based on a sequence of simulated spatial data (t0.grid to t39.grid) along a gradient (gradient.grid).

The procedure needs a sequence of binary (i.e. presence-absence) grids to compute disturbance and recovery. The algorithm can be applied both on empirical, e.g. remotely sensed data, as on simulated data. Step 1.1-1.3 explain steps to take in the data acquisition and pre-processing in case empirical data is used. When starting from a sequence of grids as the result of simulation the procedure can start from step 1.4.

1) Data and pre-processing:

1.1) Collect sequential spatial data:

First, a sequence of spatial data needs to be acquired. For instance, a time series of false color aerial images (i.e. false color composite (FCC), usually uses near infrared (NIR), red (R) and green (G) as color bands) work well for monitoring the development of vegetation. However, dependent on the ecosystem of interest and the contrast between ecosystem components, colour images collected with off the shelf (digital) camera’s (having red (R), green (G), and blue (B) color bands) can also be used. If collected analog, images need to be digitized. The sequence of snapshots should be sufficiently long containing sufficient disturbance and recovery events (see Troubleshooting for guidelines)

1.2) Georectify grids:

If the data is not yet a grid, rasterize the data. The snapshots need to be georectified to make sure that each pixel (i.e. grid cells) corresponds with the same pixel in the subsequent snapshot.

1.3) Classify grid cells:

The function requires binary data (i.e. absence = 0 and presence = 1). If the data is continuous, e.g. biomass values, the data can be classified based on a threshold. In case the data sequence consists of images, they need to be converted to a biomass proxy first before the pixels can be classified. Dependent on the image type used converting the information in the images to a proxy for biomass can be computed, among others, based on e.g.:“Normalized Difference Vegetation Index” (NDVI) for FCC or follow ref. [8] “Green excess Index” (2G-RBi).

For FCC images

As a proxy of the vegetation biomass the Normalized Difference Vegetation Index is calculated as:

NDVI = (NIR - R) / (NIR + R)

Here NIR, R, and G are the digital numbers (DN) of a 3 channels false color image, containing respectively the near infrared, red and green information.

For RGB images

The Green Excess Index can be used as proxy for biomass in RGB images and is calculated as⁸:

2G-RBi = 2*G - R - B

Here R, G, and B are the digital numbers (DN) of the 3 channels (red, green and blue respectively) in a color image.

Once the data is a continuous biomass (proxy) value presence/absence is determined on a threshold. Above a biomass threshold biomass the ecosystem component is present, otherwise it is absent.

1.4) Stack grids:

The grid of each classified snapshot needs to be stacked into one input file (Seq). The layers should be stacked in sequential order. The first layer in the stacked grid needs to contain the oldest snapshot. The stacking and the correct order of the stack is essential for the function to work.

1.5) Mapping environmental stressor:

Finally, a separate grid is needed in which information is present on the environmental stressor(s) that is/are likely to drive a critical slowing down response of the ecosystem. This grid is required as input for the analysis in step 3 (Recovery rates along a stress gradient). The grid should have the same dimensions as a single snapshot grid.

2) Timing recovery in sequential spatial data:

To compute recovery time after a disturbance the algorithm (see Fig. 2 and Fig. 3 implemented as the Matlab function fTimeRecov.m) is applied to the sequential spatial data (i.e. the stacked grid). The Matlab function provides four output files. For the final step (3) only the stacked grids ‘Tr_all’ and ‘Cens_all’ are required.

3) Recovery rates along a stress gradient:

3.1) Binning recovery times:

After the recovery times have been computed at a pixel-by-pixel basis the data can be used to test if recovery rates decline with increasing stress, that is if critical slowing down occurs in the ecosystem along a stress gradient. To test this hypothesis the observed recovery times need to be grouped (i.e. binned) in categories of stress levels. Therefore, the gradient is categorized in equally sized bins. As a rule of thumb at least five categories containing observations of recovery time should be defined. Based on the categories, the recovery times and censoring flags (i.e. the output variables ‘Tr_all’ and ‘Cens_all’) are binned.

3.2) Estimating recovery rate for each bin:

For each category the recovery rate is estimated based on the binned recovery times. Recovery rates are estimated using the maximum likelihood estimator for the exponential model. To avoid overestimation of the recovery rates, the censored recovery times should be taken into account when using the maximum likelihood estimator.

3.3) Relating recovery rates to stress gradient:

Finally, the relation between the estimated recovery rates and the stress gradient is investigated, e.g. by plotting the relation (Fig. 4). Moreover, the correlation between the stressor (i.e. bin centres of the categories) and associated recovery rate is calculated and tested for significance.

Note: An example of the full workflow (Example2.m) is provided in the ZIP package that can be found as Supplementary File 1. In the example it is assumed a data sequence of snapshots is already prepared (steps 1.1 to 1.3) and the Procedure is followed from 1.4 onward, including calculating the recovery rates along a stress gradient (step 3). This example, of which the results are shown in Fig. 4, is based on simulated spatial data (t0.grid to t39.grid) along a gradient (gradient.grid).

Most time is taken by the collection and preprocessing of the dataset itself (step 1). Once all data is collected data analysis can be performed using the algorithm. The time it takes the algorithm to run depends on the specifications of the computer used and the size of the dataset.

For the algorithm to work properly it is important to have sufficient disturbance-recovery events occurring in the sequential data set. Therefore, the ecosystem should be dynamic with disturbances occurring regularly and the sequence of snapshots should be long enough to allow for recovery to befall after a disturbance event. To guide if the sequence is sufficiently long or additional snapshots are needed in the sequence, the fraction of censored recovery times can be used as an indication. If the sequence is not long enough too much recovery times will be censored (i.e. within the sequence no recovery is observed after a disturbance). However, as the fraction of censored recovery times increases by progressing through the sequence (i.e. shorter remaining sequence), especially the fraction of censored recovery times in the oldest snapshots are relevant. To warrant reliable results, this fraction should not be too high (e.g. <50%).

If too little data is available and expanding the sequence is not possible, other methods which infer resilience indirectly from statistical properties in time or space, such as spatial variance and correlation between neighboring grid cells (i.e. spatial autocorrelation), might be considered (see e.g. ref. [1, 6]). However, these indirect methods might not work for all datasets (see ref. [7]).

Some of the observed disturbance-recovery events can be the result of the acquisition and processing of the data. Variations in the way the data is acquired (i.e. different instruments, illumination conditions due to weather and sun angle, etc.) and rectified can result in variations in the biomass indices (e.g. NDVI, 2G-RBi). Even though the data can be calibrated, some variations are unavoidable and will result in classification errors. This can, in part, be solved by removing unrealistically short recovery times from the dataset. However, if enough real disturbance-recovery events are observed the estimation of the recovery rate should still be reliable. This can be judged based on the confidence intervals of the estimate.

It is possible that no clear relation between the recovery time T_r or recovery rate λ and the stress gradient is found. Again, it is important that sufficient temporal data is used to be able to pick up a CSD signal. Check if sufficient (e.g. n > 50) disturbance-recovery events (i.e. recovery time T_r) have taken place within the stress category (i.e. bin). Increase the size of the category or the sequence of snapshots to obtain more observations, if required. The absence of a CSD signal might also be related to the use of an incorrect stress gradient and another stress is more important in driving the CSD response. But, ultimately the absence of CSD in the ecosystem can also point out that the hypothesis should be rejected and CSD is not a relevant signal of resilience in that particular ecosystem.

Overall, it is expected that the recovery rates decline with increasing stress^1-4 (Fig. 4). In the case of tidal marshes, the expected trend is that vegetation recovery rates decrease with increasing inundation time (i.e. lower intertidal elevation).

Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53-59 (2009)
Wissel, C. A universal law of the characteristic return time near thresholds. Oecologia 65, 101-107 (1984)
Van Nes, E. H. & Scheffer, M. Slow recovery from perturbations as a generic indicator of nearby catastrophic shift. Am. Nat. 169, 738-747 (2007).
Veraart, A. J. et al. Recovery rates reflect distance to a tipping point in a living system. Nature 481, 357-359 (2012)
Dai, L., Vorselen, D., Korolev, K. S. & Gore, J. Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336, 1175–1177 (2012).
Dakos, V., Kéfi, S., Rietkerk, M., van Nes, E. H., & Scheffer, M. Slowing down in spatially patterned ecosystems at the brink of collapse. Am. Nat. 177, E153-E166 (2011)
van Belzen, J., van de Koppel,, J., Kirwan, M.L., van der Wal, D., Herman, P.M.J., Dakos, V., Kéfi, S., Scheffer, M., Guntenspergen, G.R., Bouma, T.J. (2017). Vegetation recovery in tidal marshes reveals critical slowing down under increased inundation. Nature Communications, xxx
Richardson, A.D, Braswell, B.H, Hollinger, D.Y, Jenkins, J.P, Ollinger, S.V. Near-surface remote sensing of spatial and temporal variation in canopy phenology. Ecological Applications, 19(6), 1477-1428 (2009)

We thank VFA. de Witte for useful comments on earlier versions of the protocol. The work of JvB and TJB is funded by the European Commission through FP7 ENV2009-1, Contract 244104-THESEUS. JvB was further supported by the VNSC project "Vegetation modelling HPP" (contract 3109 1805) and the Schure-Beijerinck-Popping fund of the Royal Dutch Academy of Science to visit SK and VD. TJB was further supported by the NWO funded BE-SAFE project grant 850.13.011.

The authors declare no competing financial interests.

supplement0.zip
Supplementary File 1 Matlab source code Matlab sources code for the timing of recovery in sequential spatial data. The ZIP file contains the sources code as Matlab functions and two worked-out examples of the workflow to illustrate the use of the functions and the workflow (linked to Fig. 2 and 4).

Download PDF

Version 1

posted

You are reading this latest protocol version

Timing recovery of ecosystems in sequential remotely sensed and simulated data

Status:

Version 1

Abstract

Figures

Introduction

Equipment

Procedure

Timing

Troubleshooting

Anticipated Results

References

Acknowledgements

Additional Declarations

Supplementary Files

Associated Publications

Status:

Version 1

Privacy Policy

Terms of Service

Timing recovery of ecosystems in sequential remotely sensed and simulated data

Status:

Version 1

Abstract

Figures

Introduction

Equipment

Procedure

Timing

Troubleshooting

Anticipated Results

References

Acknowledgements

Additional Declarations

Supplementary Files

Associated Publications

Status:

Version 1

Privacy Policy

Terms of Service

Manage Cookie Preferences