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.  “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 as8:
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).