Detecting semi-arid forest decline using time series of Landsat data

ABSTRACT

Detecting forest decline is crucial for effective forest management in arid and semi-arid regions. Remote sensing using satellite image time series is useful for identifying reduced photosynthetic activity caused by defoliation. However, current studies face limitations in detecting forest decline in sparse semi-arid forests. In this study, three Landsat time-series-based approaches were used to distinguish non-declining and declining forest patches in the Zagros forests. The random forest was the most accurate approach, followed by anomaly detection and the Sen’s slope approach, with an overall accuracy of 0.75 (kappa = 0.50), 0.65 (kappa = 0.30), and 0.64 (kappa = 0.30), respectively. The classification results were unaffected by the Landsat acquisition times, indicating that rather, environmental variables may have contributed to the separation of declining and non-declining areas and not the remotely sensed spectral signal of the trees. We conclude that identifying declining forest patches in semi-arid regions using Landsat data is challenging. This difficulty arises from weak vegetation signals caused by limited canopy cover before a bright soil background, which makes it challenging to detect modest degradation signals. Additional environmental variables may be necessary to compensate for these limitations.

KEYWORDS:

• forest decline
• Landsat time series
• random forest
• anomaly
• Sen’s slope
• semi-arid

Introduction

Recent research highlights the rapid decline of (semi-)arid forests worldwide, emphasizing the need for a better understanding of the causes of forest decline (e.g. Anderegg et al., 2015). While many studies have investigated deforestation in these regions (e.g. Hoyos et al., 2013; Le Polain de Waroux & Lambin, 2012), forest decline has received relatively little attention (Le Polain de Waroux & Lambin, 2012; Li et al., 2009). The decline of forests is often the result of degradation, which negatively impacts their functional and structural properties (Vásquez-Grandón et al., 2018; Wang et al., 2020). Forest decline can also decrease the regeneration capacity of a forest by reducing the number of available reproductive trees and altering micro-environmental conditions, in addition to the direct loss of trees (Cailleret et al., 2019). Moreover, forest decline can have long-term consequences for ecosystem services and regulatory services (Le Polain de Waroux & Lambin, 2012).

Forest decline is often caused by interacting biotic and abiotic factors (Hosseini et al., 2017). Among abiotic factors, droughts can have a particularly damaging effect on forests (Sulla-Menashe et al., 2014, Le Polain de Waroux & Lambin, 2012. (Semi-)arid ecosystems must deal with more frequent and severe droughts, which, in combination with other biotic stressors, can lead to significant forest decline (Sánchez-Pinillos et al., 2021). While forest decline occurs globally, (semi-)arid ecosystems are particularly vulnerable, with an increasing tree mortality rate (David et al., 2022).

The Zagros forests in western Iran have been suffering from significant forest decline due to climatic extremes, wildfires, and overexploitation by local populations. Their low tree species diversity, dominated by a few oak species, in particular Brant’s Oak (Quercus brantii Lindl.), makes them susceptible to various forest pathogens (Goodarzi et al., 2016; Moradi et al., 2021). Earlier studies conducted in the region often focused on deforestation but not forest decline (e.g. Goodarzi et al., 2016; Hosseini et al., 2017; Jahanbazy Goujani et al., 2020; Moradi et al., 2021). Understanding forest decline in the region is crucial for effective forest management and adaptation initiatives but is still lacking for larger continuous areas (Moradi et al., 2021).

Field data on forest decline is often limited due to cost, particularly in less accessible and extensive areas (Diao et al., 2020; Lausch et al., 2016; Shafeian et al., 2021) such as the Zagros forests. Studies performed by, for example, Wang et al. (2020) and Senf et al. (2020), indicate that remote sensing (RS) is an alternative to assess forest decline in the field and over large regions. Comparing multi-temporal satellite images enables the detection of changes in vegetation conditions in a particular area (Wang et al., 2020). Various studies have used vegetation indices (VIs) like normalized vegetation index (NDVI) (Jin et al., 2016) and normalized burn ratio (NBR) (Francini & Chirici, 2022) to detect forest decline. However, multispectral VIs are known to have limitations in monitoring biophysical and biochemical vegetation properties in (semi-)arid ecosystems (e.g. David et al., 2022). For instance, NDVI is sensitive to green components, not to woody components, and factors like moisture content and species composition can influence the relationship between biomass and NDVI. The Enhanced Vegetation Index (EVI) was developed to overcome the drawbacks associated with variable background reflectance and atmospheric interference in the NDVI (David et al., 2022). The EVI provides more comprehensive information about changes in vegetation over time and space and is particularly suitable for assessing vegetation in semi- and arid ecosystems. Meanwhile, Bae et al. (2022) and Li et al. (2022) note that studies on forest decline commonly utilize very high-resolution remote sensing data, which can be costly and challenging to acquire for extensive areas.

As mentioned above, many RS studies in the context of forest monitoring have focused on deforestation rather than forest decline. The identification of deforestation (using optical RS) is based on the distinction of forest and non-forest land-cover types, such as open landscapes, bare soil, crops, and settlements (Hoekman et al., 2020). Among several types of RS data, multi-spectral optical images are the most cost-effective option for monitoring large forest areas. Time series of multispectral optical satellite data, such as Landsat and Sentinel-2, have facilitated the development of various methodologies for mapping forest changes (Giannetti et al., 2020). In particular, the Landsat time series has been recognized as a valuable data source for tracking forest decline and disturbances due to their free accessibility, relatively high spatial resolution (30 m), long and consistent acquisition record, and accessibility (Diao et al., 2020; Dutrieux et al., 2015; Giannetti et al., 2020; Rodman et al., 2021; Senf et al., 2020; Zhu et al., 2020).

A variety of algorithms have been introduced for identifying forest disturbances throughout time, such as Continuous Change Detection and Classification (CCDC; Zhu & Woodcock, 2014), Landsat-based Trends in Disturbance and Recovery Detection (LandTrendr; Kennedy et al., 2010), and Breaks for Additive Season and Trend (BFAST; Verbesselt et al., 2010). Nevertheless, the majority of these algorithms are more advantageous for monitoring forest changes or deforestation than assessing forest decline or minor changes in forest status (Zhu et al., 2020). The difficulty of detecting subtle signals of degradation is a challenge for all algorithms, as spectral signals associated with tree mortality depend on various factors, including the number of canopy layers and forest cover (FC) (Hoekman et al., 2020).

This study seeks to assess the effectiveness of well-established methods like random forest (RF), anomaly detection, and Sen’s slope analysis for identifying and mapping forest decline in the Zagros semi-arid forests.

Study area

The study area is located in the Chaharmahal and Bakhtiari provinces of Iran, in the Zagros semi-arid forests, with average rainfall of 250–800 mm and temperatures ranging from 9–25°C (Attarod et al., 2016). The study area is located in the southwestern part of the province (Figure 1). The Zagros forests play an important role in preserving the unique vegetation and habitats indigenous to the region (Sagheb-Talebi et al., 2014; Shafeian et al., 2021). The most prominent and widespread tree species found in the area include Brant’s oak (Quercus brantii Lindl.), which is partially mixed with Quercus infectoria G.Olivier, Quercus libani G.Olivier, wild pistachio (Pistacia atlantica Desf.), and Acer monspessulanum L. In addition, one can commonly find Crataegus spp., Amygdalus spp., and Pyrus spp. in the area (Daneshmand Parsa et al., 2016; Erfanifard et al., 2014).

Figure 1. (A) location of Iran; (B) digital elevation model of Chaharmahal and Bakhtiari province; (C) sampled area covered in the green square (AOI); (D) an example UAV image of the study region.

Methodology

The workflow for the study is outlined in Figure 2 and further described below.

Figure 2. The workflow of the present study (steps 1, 2, and 3 were performed in GEE, and step 4 was performed in R.).

Remote sensing analysis

In this study, we looked at a variety of methodical approaches developed for analyzing time series and image stacks of optical Landsat images to see if we could distinguish between declining and non-declining forest regions. In the following, we will first describe the pre-processing steps and then introduce the tested methodical approaches.

Landsat time series and image processing

The GEE platform was used to process the Landsat image time series. The analysis employed the Landsat Surface Reflectance Tier 1 data from Landsat 4 Enhanced Thematic Mapper (ETM), Landsat 5 ETM, Landsat 7 ETM+, and Landsat 8 Operational Land Imager (OLI) (IDs, respectively: LANDSAT/LT04/C01/T1; LANDSAT/LT05/C01/T1; LANDSAT/LE07/C01/T1; and LANDSAT/LC08/C02/T1_L2). We filtered the image collection to include all images that were available for the summer periods of 1986 to the summer periods of 2021. We examined three definitions for the summer periods to see if there was any effect on the results: (1) Landsat data from the months between July and the end of September, (2) July and the end of August, and (3) only August, i.e. the dry season, with the highest spectral signal differences between trees and other vegetation, crops, and grass, as suggested by Symeonakis et al. (2018). Then we calculated the annual median reflectance value for each pixel and band from the remaining cloud-free image stack, which resulted in a cloud-free, high-quality mosaic. In order to detect forest decline, a set of VIs was computed using these annual composites (Martin-Ortega et al. 2020). The examined VIs are summarized in Table 1. In the case of the RF classifications (see below), we did not use the annual mosaics but calculated a single mosaic using all Landsat data from the corresponding time periods.

Table 1. Equations of the used vegetation indices (VIs).

VI	Name	Equations	Reference
NDVI	Normalized Difference Vegetation Index	(NIR – RED)/(NIR + RED)	Rouse et al. (1973)
GNDVI	green NDVI	(NIR- GREEN)/(NIR+ GREEN)	Gitelson et al. (1996)
kNDVI	Kernel NDVI	tanh ((NIR-RED)^2/2 sigma) sigma = NIR + RED (0.5)	Camps-Valls et al. (2021)
NRGI	Normalized Green-Red Vegetation index	(GREEN – RED)/(GREEN + RED)	Buce Saleh et al. (2019)
EVI	Enhanced Vegetation Index	2.5 ((NIR – RED)/(NIR + C1 × RED − 7.5 × BLUE + 1))	Liu and Huete (1995)
SR	Simple Ration	(NIR)/(RED)	Jordan (1969)
SLAVI	Specific Leaf Area Vegetation Index	(NIR)/(RED+SWIR)	Lymburner et al. (2000)
NDWI	Normalized Difference Water Index	(GREEN-NIR)/(GREEN+NIR)	Hardisky et al. (1983)
NDMI	Normalized Difference Moisture Index	(NIR – SWIR)/(NIR + SWIR)	United States Geological Survey
NDTI	Normalized Difference Tillage Index	(SWIR1 - SWIR2)/(SWIR1 + SWIR2)	van Deventer et al. (1997)
NBR	Normalized Burn Ratio	(NIR- SWIR2)/(NIR + SWIR2)	United States Geological Survey
ARVI	Atmospherically Resistant Vegetation Index	(NIR-(2×RED)-BLUE)/(NIR+(2×RED) +BLUE)	Kaufman and Tanre (1992)
GCI	Green Chlorophyll Index	NIR/GREEN − 1	Gitelson et al. (2005)
GLI	Green Leaf Index	(GREEN-RED-BLUE)/(2× RED+GREEN+BLUE)	Gitelson et al. (2002)

Table 1. Equations of the used vegetation indices (VIs)

Identifying forest decline

The datasets described in Section 4.1 were used as inputs to the following approaches:

Random Forest classification

The RF classifier (Breiman, 2001) was used to classify the study area into declining and non-declining classes. Previous studies (e.g. Belgiu & Dragut, 2016) successfully utilized the RF classifier for satellite data classification and forest disturbance analysis. In this study, several RF classifiers were trained with Landsat data covering various year intervals, and each classifier was validated using iterative splitting of the reference data (see below) into training and validation sets (with 100 repetitions) to distinguish between declining and non-declining areas throughout the study area. This was done using the randomForest package in R with the number of trees (ntree) set to 500 and the second parameter, mtry, set to default settings. In order to prevent the issue of having two different classes in one pixel (a few of the available field plots were located too close to each other), we built a function during the RF classification training to eliminate any duplicated reference pixels, resulting in the removal of 164 points. As a result, a total of 719 reference points remained. As inputs, we initially used the Landsat time series with images from summer 1986 to summer 2021 as well as from 2000–2021 and from 2010–2021 and examined the three input options described above: (1) images from July to the end of September; (2) images from July to the end of August; and (3) images from August. We examined several time periods since we did not have a clear idea of when the field-observed decline processes had started. According to local experts and previous literature, the decline accelerated after the year 2000 (e.g. Ghanbari Motlagh & Kiadaliri, 2021), but no detailed information exists. Since the RF was showing very stable results, we then examined even more time periods reaching further in the past (1986–1990, 1986–1999, 1991–1995, 1996–2000, 2001–2005, 2006–2010, 2011–2015, and 2016–20201).

We used all of the VIs listed in Table 1 as predictors, as well as the seven original bands of the Landsat data. Then, in addition to using all the predictors separately for each input period, we also applied the feature selection algorithm VSURF in R (Genuer et al., 2015). The first step in the variable selection process with VSURF is to remove all non-essential, i.e. highly interrelated variables, from the dataset. The second step is to select all variables that are related to the response and help the model perform well. The third step refines the collection for prediction purposes by removing redundancy in the range of variables chosen in the second step (Genuer et al., 2015). Table 2 shows the predictor variables that were selected for each of the three input options. To obtain the final maps, we used all available reference data for classification. The result was a binary map of binary raster layers showing two classes (non-declining, and declining). To validate the RF model, the reference data set was repeatedly divided into 70% training samples and 30% validation samples (Shafeian et al., 2021). We obtained kappa and overall accuracy from the validation samples.

Table 2. The results of VSURF on RF for each season. The selected predictors of Landsat 8 images are written with a suffix (_LS8).

Landsat data for summer periods	Selected predictors
August 1986–2021	B1, B2, B3_LS8, NDVI, B2_LS8, B3, NDWI, B5, NRGI, NDTI
August 2000–2021	B2, B1, B3_LS8, B3, NDMI, NDWI, B6, B1_LS8, NDTI
August 2010–2021	B2, B1, B3_LS8, B3, B2_LS8, NDWI, B1_LS8, B6, B5
July- August 1986–2021	B1, B2, B3_LS8, B3, NDVI, NDWI, NRGI, B5, B6, NDTI
July- August 2000–2021	B1, B2, B3, NRGI, NDTI, NDWI
July- August 2010–2021	B2, B1, B2_LS8, NDWI, B6, B5, NRGI, B1_LS8
July- September 1986–2021	B1, NDMI, B3, B2_LS8, NDWI, B6, NDTI
July- September 2000–2021	B2, B1, NDMI, B3, B2_LS8, NRGI, NDWI, NDTI
July- September 2010–2021	B2, B1, NDMI, B3_LS8, B2_LS8, B3, B6, NDWI

In addition to the Landsat series, we also applied an RF classification to Sentinel-2 imagery captured between 2015 and 2021 using the same procedure (see Appendix 1). The spatial resolution of Sentinel-2 imagery is higher than that of the Landsat series at 10 × 10 m, making it a good option for forest decline detection. However, these images have only been available since 2015, limiting the ability to perform time series analysis such as anomaly detection over a longer period.

Anomaly analysis

As a second approach, we applied a pixel-wise anomaly analysis using the vegetation indices. For this purpose, we used the annual Landsat mosaics from 1986–2021. Again, we examined the three input options to define the summer period. Following the calculation of the vegetation indices for each image collection as described above, a median value filter of a temporal moving window with a three-year interval was applied to the annual time series of the indices, that is, 1986–1989, 1987–1990,… 2018–2021. Then all these images were stacked into one single image time-series stack. This step was conducted to smooth the datasets and avoid the potential influence of particularly high or low values in an individual year, which may have been related to extraordinary weather conditions or a limited number of cloud-free sample pixels. For each pixel of this time series stack, the “all-time high” and “all-time low” values, defined as the 90^th and 10^th percentiles, respectively, were calculated. Then the 10^th percentile was subtracted from the 90^th percentile to check how widely the VI value of a given pixel varied over the whole time period. This range was then used to define an anomaly threshold. For example, the range was multiplied by 0.5, and the resulting value was added to the value of the 10^th percentile of the VI time series. We then considered an area to be declining if the current VI value (see below) was below the 10^th percentile plus the 0.5 × range value. The corresponding areas were identified by applying the pixel-specific thresholds, resulting in a binary raster layer (zero or black = non-declining, and one or orange = declining) for each VI. We varied the factor with which the range was multiplied and examined the results for all factors from 0.15 to 0.9 with steps of 0.05. The current VI value for each pixel was obtained from a cloud-free mosaic image of the current status of the area calculated from Landsat images of the years 2018–2021. These calculations were performed for all indices given in Table 1.

Sen’s slope analysis

Sen’s slope, or Theil-Sen’s slope (Sen, 1968), is a non-parametric trend analysis that is resistant to outliers and rejects odd values without changing the slope (Correa‐Díaz et al., 2019). We used this approach because it is useful to identify whether a regular time series has a statistically significant positive or negative trend. The Sen’s slope can account for the magnitude of changes (Reygadas et al., 2019). Thus, after importing a cloud-free Landsat time series and calculating the aforementioned VIs (e.g. KNDVI, ARVI, etc.), the Sen’s slope was calculated using the ee.Reducer.sensSlope() function available in GEE. With this algorithm, we could estimate the trend of changes in the region over the previous 36 years using mosaics of Landsat imagery. We again examined all three input options with respect to the summer period. The results were then exported as a raster layer for each summer period and VI (14 VIs and three summer periods gave a total of 42 single-band raster layers). To translate the Sen’s slope results into a binary map that we could compare with our reference data, we defined a function to calculate 500 different threshold values for each VI (based on the range of Sen’s slope values obtained for the VI) and determined the best threshold for detecting the declining areas (according to our reference data). We calculated the thresholds for each VI individually since each index had a different range of values. We ran this for each index separately over three summer seasons. And finally, the threshold with the highest overall accuracy was recorded (the results of the selected threshold for each index are summarized in Table S3 in Appendix 3).

Visualization of vegetation indices time series

To guide the interpretation of the results, particularly with respect to the anomaly and Sen’s slope analysis, we additionally visualized the VI time series from 1984–2022 for the two most accurate (according to overall accuracy) VIs using the pixels representing the reference data locations for declining and non-declining areas from the spring 2022 campaign. We calculated the median and standard deviation of all pixels coinciding with the corresponding reference plots.

Reference data

So far, no official and regular forest inventories exist for the majority of the Zagros region. Therefore, we tried to maximize the number of field samples from which the state of the forest is known by merging various datasets collected between 2021 and 2022 in the region.

The first part of our reference data was collected during a field campaign in the spring of 2022, during which declining and non-declining areas were inventoried. After consultation with local experts to define the study area, a stratified random sampling approach was used to select areas for sampling. Two factors were considered to stratify the study area: Illumination Condition (IC) and Forest Cover (FC). The IC for each pixel in the study area was calculated using NASA SRTM Digital Elevation imagery available on Google Earth Engine (ID: USGS/SRTMGL1_003). Although IC is not directly related to forest decline, it was assumed that it could affect satellite signals. In addition, an existing forest cover map for the area provided by Shafeian et al. in 2021 was used. FC is known to significantly affect remote sensing signals. The study area was divided into four classes each for IC and FC, resulting in a total of sixteen layers (combinations of IC and FC classes). To further account for accessibility, a 400-meter buffer was applied around roads, and random points were generated within these buffers with a minimum distance of 100 meters between each point. Nine random points were selected for each stratum, for a total of 144 candidate positions for field verification. Approximately half of these points were eventually visited and field recorded. During the field campaign, information on the condition of the forests was collected and classified into two general categories: “declining” and “not declining”. Overall, this process allowed the collection of reference data that can be used to validate and calibrate remote sensing data to accurately assess forest decline in the study area.

The second source of reference data was a field campaign conducted in September 2021. The locations of the field surveys were determined by first stratifying the study area into two zones of declining and non-declining areas according to an EVI time series as observed for the last 36 years using Landsat data. To identify these zones, we applied an EVI anomaly approach based on data between 1986 and 2021. The resulting binary image indicated declining (orange color) and non-declining (black color) areas. The field plot locations were chosen at random within the defined declining area. Plots in inaccessible areas were dropped and replaced with plots that were located closer to the roads. During the field campaign, we visited and recorded 43 sites. Any signs of deterioration, such as foliar deficiencies, branch mortality, or unusual leaf coloration, were documented (Figure 3 shows some examples of the decline signs).

Figure 3. Example of declining trees in the study area. Images are taken during the field campaign in September 2021. The exact reasons for the decline were often unknown. In some cases signs of the presence of damaging insects or pathogens were visible, but it was impossible to determine whether these were a consequence or the reason for the decline.

The third source of data is a set of additional GPS locations for declining and non-declining tree groups collected by local forest service officers in September 2022. All the reference points used in this study are jointly plotted in Figure 4. In total, we had positions for 461 declining and 422 non-declining plots.

Figure 4. The reference points, declining areas in red and non-declining in green.

Results

Figure 5 displays binary decline maps from all three approaches (Figure 5a: Sen´s slope, Figure 5b: anomaly, and Figure 5c: RF). To ensure that the analysis focused only on relevant land cover (LC) types, we used an LC map to mask out irrelevant classes such as water, bare soil, and agriculture, using the LC map of Shafeian et al. (2021) (Figure 5d). This left only rangeland, forest, and plantations for further analysis, allowing for a more accurate and targeted assessment of the study area.

Figure 5. A) Sen’s slope (NDMI index with threshold = 107.93 (August)), B) anomaly (NRGI 20th percentile (July-August)), and C) RF (July-August) for 1986 to 2021. Declining forest areas are marked in orange, while non-declining areas are depicted in black. D) LC map (masked classes in white and the relevant classes in green), E) google Earth view of the study area.

Visual inspection revealed that certain areas in the southwest were consistently detected as declining zones in all approaches, but each approach contained some errors. RF performed better than the other two approaches. However, overall, the examined approaches performed moderately to poor, as also reflected in the noticeable differences between the maps.

Validation and accuracy of forest decline detection

Random forest

The RF accuracy for different summer seasons was stable, with the median overall accuracy out of 100 iterations for splitting training and testing reference points ranging between 0.72 and 0.76 (Figure 6). The best accuracy was achieved using input datasets containing images taken between July 1 and the end of August for the years between 2000 and 2021 (overall accuracy = 0.76, kappa = 0.50). The user’s accuracy for the non-declining class was approximately 0.80, whereas it was approximately 0.72 for the declining class. The producer’s accuracies for the two classes were 0.78 and 0.72, respectively (see Figure 7). Images for the period between 1986 and 2021 had the second-highest median values for overall accuracy and kappa with 0.75 and 0.48, respectively (see Figure 6). The user´s and producer´s accuracies for the non-declining class of this analysis were around 0.76 and 0.77, respectively, and for the declining class, 0.70 and 0.69, respectively (see Figure 7).

Figure 6. The upper row of the plot shows the overall accuracy and the kappa values of the RF forest classifications for the years 1986–2021 (86–21), 2000–2021 (00–21), and 2010–2021 (10–21) for the three examined summer periods shown in different colors (August, July–August, and July–September). The lower row of the plot shows the same analysis for the years 1986–1999 (86–99), 1991–1995 (91–95), and 1996–2000 (96–00). The fine dashed lines indicate the overall mean across all bean plots.

Figure 7. The upper row of the plot shows the user and producer accuracies for declining and non-declining classes of RF classifications for the years 1986–2021 (86–21), 2000–2021 (00–21), and 2010–2021 (10–21) for the July–August summer period. The lower row of the plot shows the same analysis for the years 1986–1999 (86–99), 1991–1995 (91–95), and 1996–2000 (96–00). The fine dashed lines show the overall mean across all bean plots.

More detailed RF results for the years 1986–2021, 2000–2021, and 2010–2021 for all the summer seasons are shown in Figure S3 in Appendix 1 (i.e. specificity, F1 Score, etc.). The RF results for the years further in the past, including 1986–1999, 1990–1995, and 1996–2000, for the three examined summer periods and their comparison with other year interval results are shown in Figure 6. The highest accuracy was obtained for the input dataset with images from July 1 to the end of August for the years between 1996 and 2000 (overall accuracy = 0.76, kappa = 0.50). Also, the user´s accuracy of the non-declining class and the declining class for this analysis was almost 0.78 and 0.72, respectively, and the producer´s accuracy of the non-declining and declining classes was almost 0.80 and 0.71, respectively (Figure 7). The user´s accuracy of the non-declining class and the declining class for this analysis for other summer seasons (August and July-September) can be found in Figures S1 and S2 in Appendix 1.

The RF results based on VSURF-selected predictors hardly improved accuracy and kappa values. The selected predictors are summarized in Table 2.

Table 2. The results of VSURF on RF for each season. The selected predictors of Landsat 8 images are written with a suffix (_LS8).

Anomaly analysis

The best results for anomaly analysis were obtained when using the annual composite time series comprising images from July to August. In this configuration, the GLI (30^th percentile) achieved a kappa value of 0.30 and an overall accuracy of 0.65. The user accuracies for the non-declining and declining classes were 0.61 and 0.69, respectively, while the producer accuracies for these classes were 0.71 and 0.59, respectively. In contrast, when using images from July and September, the NDWI (15^th percentile) and the NDWI (20^th percentile) had the lowest accuracy values, with kappa values and overall accuracy of −0.13 and 0.45, respectively. The Kappa values of anomaly analyses for three summer periods at different percentiles for all the VIs are shown in Figure 8, and the overall accuracies are shown in Figure S9, Appendix 2. Further details on the results can be found in Table 3. The detailed results of this analysis for user´s and producer´s accuracies for declining and non-declining classes can be found in Table S2, Appendix 2.

Figure 8. Kappa values of anomaly analysis at different percentiles for all the VIs of August (left panel), July to August (center panel), and July to September (right panel).

Table 3. A summary of the results (overall accuracy (in bold) and kappa value) of the anomaly approach for different VIs at different threshold percentiles (15–90) for July to August. The percentile (P) with the highest overall accuracy is marked with an asterisk, and the percentile with the highest overall accuracy is underlined.

P VI	15	20	25	30	40	50	60	70	80	90
NDVI	0.54 0.11	0.54 0.11	0.56 0.15	0.56 0.15	0.59 0.20	0.59 0.21	0.60 0.22	0.60 0.21	0.61* 0.24	0.60 0.20
GNDVI	0.47 −0.002	0.47 −0.002	0.47 −0.002	0.47 −0.002	0.47 0	0.48 0.02	0.53 0.09	0.55 0.14	0.58 0.19	0.58 0.18
KNDVI	0.50 0.04	0.52 0.09	0.53 0.10	0.53 0.10	0.54 0.13	0.56 0.15	0.57 0.17	0.57 0.17	0.59 0.20	0.58 0.18
NRGI	0.62* 0.25	0.63* 0.26	0.62* 0.23	0.60 0.18	0.57 0.12	0.57 0.11	0.55 0.07	0.54 0.04	0.53 0.02	0.52 0.009
EVI	0.53 0.10	0.53 0.10	0.54 0.11	0.55 0.12	0.55 0.13	0.55 0.14	0.57 0.16	0.59 0.19	0.59 0.20	0.61* 0.23
SR	0.50 0.04	0.52 0.07	0.53 0.10	0.53 0.10	0.54 0.12	0.56 0.15	0.57 0.16	0.57 0.16	0.59 0.20	0.58 0.18
SLAVI	0.52 0.07	0.54 0.12	0.55 0.14	0.55 0.14	0.56 0.15	0.58 0.18	0.58 0.19	0.59 0.20	0.60 0.22	0.60 0.22
NDWI	0.47 −0.007	0.49 −0.005	0.50 −0.003	0.50 −0.002	0.52 0.007	0.52 0.004	0.52 0.002	0.52 0.002	0.52 0.002	0.52 0
NDMI	0.52 0.08	0.55 0.14	0.56 0.15	0.56 0.15	0.59 0.20	0.60 0.22	0.61* 0.24	0.61 0.23	0.60 0.22	0.59 0.18
NDTI	0.47 −0.002	0.47 −0.001	0.47 −0.002	0.47 −0.006	0.47 −0.009	0.47 −0.007	0.47 −0.007	0.47 −0.006	0.47 −0.002	0.48 0
NBR	0.47 −0.006	0.47 −0.002	0.47 −0.01	0.47 −0.01	0.46 −0.002	0.47 −0.008	0.47 −0.005	0.46 −0.02	0.47 −0.01	0.50 0.03
ARVI	0.50 0.04	0.50 0.04	0.50 0.05	0.52 0.08	0.53 0.10	0.54 0.12	0.55 0.14	0.56 0.15	0.57 0.16	0.56 0.14
GCI	0.47 0	0.47 −0.002	0.47 −0.002	0.47 −0.002	0.47 −0.002	0.47 −0.006	0.47 −0.002	0.49 0.02	0.50 0.04	0.53 0.09
GLI	0.56 0.14	0.60 0.22	0.62* 0.26	0.65* 0.30	0.61 0.22	0.60 0.19	0.60 0.19	0.60 0.12	0.57 0.07	0.55 0.05

Table 3. A summary of the results (overall accuracy (in bold) and kappa value) of the anomaly approach for different VIs at different threshold percentiles (15–90) for July to August. The percentile (P) with the highest overall accuracy is marked with an asterisk, and the percentile with the highest overall accuracy is underlined.

Sen’s slope

The highest overall accuracy for Sen´s slope was achieved by GLI and NDMI for August, with approximately 0.64 (kappa = 0.53) and 0.63 (kappa = 0.56), respectively. The Kappa values of the Sen´s slope for all the VIs of three summer periods at different thresholds are shown in Figure 9; the overall accuracies are shown in Figure S10, Appendix 3. The lowest overall accuracy (0.39) was obtained by GCI and GNDVI for both August and July – August. Detailed information for the median overall accuracy and kappa values for each VI is shown in Table 4. The corresponding results for user´s and producer´s accuracies for declining and non-declining classes can be found in Table S4, Appendix 3.

Figure 9. Kappa values of the Sen´s slope for all the VIs of August (left panel), July to August (center panel), and July to September (right panel).

Table 4. A summary of the most accurate results of Sen’s slope for different VIs for August after trying 500 different VI-specific thresholds; the highest overall accuracies among all the VIs are marked with an asterisk sign.

VI	Overall Accuracy	Kappa
NDVI	0.51	.01
GNDVI	0.47	−.02
KNDVI	0.51	−.004
NRGI	0.52	.04
EVI	0.51	.004
SR	0.51	.48
SLAVI	0.52	.06
NDWI	0.52	.02
NDMI	0.55*	.13
NDTI	0.49	.002
NBR	0.52	.009
ARVI	0.51	−.01
GCI	0.47	−.014
GLI	0.53	.065

Table 4. A summary of the most accurate results of Sen’s slope for different VIs for August after trying 500 different VI-specific thresholds; the highest overall accuracies among all the VIs are marked with an asterisk sign.

Time series of NDMI and GLI

NDMI and GLI were the best VI indicators of forest decline in the semi-arid forests studied. Figure 10 shows the time series (1984–2022) of Landsat-based NDMI and GLI values with median and standard deviation values for non-declining and declining reference points. Differences in VI over time for declining and non-declining regions were subtle, with no clear trend visible. GLI showed a downward trend after 2000, but it was visible in both declining and non-declining areas, making it hard to interpret. The trend also coincided with the launch of Landsat-7 and may relate to a technical artifact.

Figure 10. Time series of NDMI and GLI for August: the median (solid line) and standard deviation (dashed line) of declining pixels (red) and non-declining pixels (blue) are depicted.

Discussion

In this study, we investigated whether well-established methods such as RF, anomaly detection, and Sen’s slope analysis can detect and map forest decline in sparse and patchy semi-arid forests using the Landsat time series. In the following, we will discuss the technical results for the applied VIs and examined approaches, followed by a general discussion exploring challenges in detecting and mapping vegetation decline in mountainous semi-arid regions.

Spectral bands and vegetation indices

Choosing appropriate RS indicators is critical for successfully detecting forest decline. In our study, the performance of VI-based approaches varied significantly depending on the VI used. This finding is consistent with previous studies on forest decline (Yu et al., 2021). Certain VIs, such as NDMI and GLI, appear to be more effective at detecting subtle spectral changes in (semi-)arid forests than other VIs. Such indices, which include a SWIR band, are more closely related to physiological variables such as hydraulic conductivity and water potential than greenness indices, particularly in (semi-)arid regions (e.g. Marusig et al., 2020; Moreno-Fernández et al., 2021). However, the usefulness of NDMI has also been demonstrated in studies of ecosystem types other than (semi-)arid regions. For example, Li et al. (2022) successfully used the NDMI to detect forest disturbance in a subtropical forest.

The NDMI is calculated using both NIR and SWIR bands. Both bands are known to be directly related to key vegetation properties, including the high NIR reflectance of parenchyma cells and the high sensitivity of the SWIR band to leaf water content. Hence, NIR and SWIR bands are frequently reported to be important predictors in studies examining vegetation dynamics using satellite data. For example, Meyer et al. (2019) found the NIR and SWIR bands as important variables for tracking changes in leaf area in satellite data time series. Similarly, Moreno-Fernández et al. (2021) found for Landsat time series data that forest decline was linked to the trend component of the spectral index series, with the wetness index NDMI exhibiting declines earlier than the greenness indices (EVI, and NDVI). Negative trends occurred earlier for wetness indices like NDMI than for greenness indices, indicating that the former may be better suited for detecting declines.

Gu et al. (2008) and Li et al. (2022) emphasized the significance of SWIR bands in identifying forest decline. Li et al. (2022) observed a decline in accuracy of up to 26% when excluding SWIR bands from their modified continuous monitoring of land disturbance approach. In our study, the SWIR bands and VIs that included SWIR bands were frequently chosen by the VSURF variable selection for various RF classifications, underscoring the importance of SWIR spectral bands in detecting semi-arid tree decline. Furthermore, the NDWI index, calculated from the green and NIR bands, was frequently selected but performed poorly in detecting anomalies during the July to August and July to September periods. Previous research, such as that by Das et al. (2023), has reported the relevance of NDWI in identifying declining and stressed forests.

In addition to NDMI and NDWI, GLI showed promising results in our study area, despite being rarely used in prior studies to detect forest decline. GLI has been utilized successfully for other vegetation monitoring purposes, such as detecting vegetation using aerial images (e.g. Eng et al., 2019). The GLI uses blue, green, and red bands. It was one of the best indices in anomaly and Sen´s slope approaches but was not selected by VSURF as an important variable for the RF classifications. However, all the bands that are required to calculate GLI were frequently selected by VSURF for the RF classifications in all summer periods. The selection of the green band partly agrees with Higginbottom et al. (2018), who reported the green band as the second-best individual band for monitoring canopy properties in sparse woody vegetation.

NDVI was only considered a medium-important index in anomaly and Sen’s slope approaches and was not selected in any of the VSURF runs. Therefore, NDVI does not appear to be a reliable indicator of forest decline in semi-arid vegetation, which is consistent with Madonsela et al. (2018) but contrasts with Camarero et al. (2015), who suggested that NDVI is a dependable indicator of drought-induced productivity decline.

It should be noted that leaf structure and vegetation type are also important to be considered when choosing the appropriate VIs. Only a few studies have investigated VIs for different leaf structures. Croft et al. (2014) evaluated several VIs to estimate leaf chlorophyll content in different leaf and canopy structures, including broadleaf and coniferous trees across a Canadian test site. They found that the canopy-scale relationship between spectral indices and foliar chlorophyll content was particularly strong for broadleaf samples compared to coniferous trees. However, in our study area, the tree species are mostly broad-leaved Brant oaks, so this was not an issue to consider. But, in heterogeneous forests with both types of leaf structure, this could be an important factor in selecting appropriate VIs.

The importance of acquisition time windows of satellite images

This study found that the best time window for detecting forest decline differed slightly depending on the approach used. The RF and anomaly approaches were more accurate in detecting declines in images captured from July to August, while Sen’s slope was most accurate for images captured in August only. This may be due to the particularly dry conditions during these months, which lead to greater spectral differences between trees and other vegetation. This finding is supported by a study by Symeonakis et al. (2018). We observed that RF performed better for the predictors based on July and August images compared to Sen’s slope, which worked slightly better for only August images. However, we did not find any explanation for this observation, and the performance differences were subtle.

Differences in performances between the examined approaches

In our study, RF outperformed other methods in terms of accuracy, which may be due to its use of multiple predictor variables. Although the overall performance was moderate, the balanced class-specific accuracy suggests that RF was able to capture spectral information related to declining forest areas. Interestingly, we observed no notable differences in performance between RF models trained with Landsat data from different time periods (1986–2021, 2000–2021, and 2010–2021), despite the reference data being collected in 2021 and 2022. This observation was unexpected, as we expected that using more recent data would lead to better results. Results for different RF models with various additional year intervals are provided in Appendix 1. According to local experts, the widespread severe oak decline in the study area only emerged after 2010, which may indicate that the subtle spectral differences enabling to differentiate declining from non-declining areas in this study may not necessarily relate to changes in the spectral signal due to decline of the trees but rather to general differences in the site conditions on which declining and non-declining forest areas are located. This may very well relate to other environmental parameters such as differing soil conditions, co-occurring vegetation, topography, or a combination of both.

The anomaly detection and Sen’s slope analysis approaches performed notably worse than RF in the study. A direct comparison between RF and the other two approaches (anomaly and Sen’s slope) is not feasible due to the distinct methodologies employed in these analyses, in particular concerning the input data and the validation. RF, as a supervised classification method, a portion of the reference data (equal to 70% of all reference data) is directly utilized for model training. Contrarily, in the anomaly analysis and Sen’s slope approach the reference data is only used for the validation and accuracy assessment; even though determining the optimal threshold based on the reference data is similar to applying the reference data as training data. Nevertheless, the better performance of RF is somewhat expected due to the notably increased number of predictors as compared to the other two tested approaches relying on a single VI.

Several earlier studies found Sen’s slope approach useful for trend detection in vegetation degradation. For instance, Zhang et al. (2022) used three trend analyses, including Sen’s slope and Mann-Kendall tests, to assess the changing trend and degradation of the environmental quality of the Loess Plateau in China. Sen’s slope and Mann-Kendall analysis were observed to be robust and insensitive to outliers, and their combination was suggested to offer advantages over a simple linear regression. The Loess Plateau has experienced serious erosion processes, resulting in clearer spectral change signals in satellite imagery, which may explain the successful application of trend-analysis techniques in this area.

Our study found that anomaly detection and Sen’s slope approach, which are based on temporal analysis, perform poorly. This suggests that we may not be able to detect a clear change in the spectral signal caused by forest decline. This though is also backed up by the visualization of the time series of the two most accurate VIs which showed similar signals for declining and non-declining areas, which makes a reliable separation with only one single VI’s information difficult (Figure 10). The downward trend observed for GLI after 2000 may relate to the start of a general vegetation decline but could also be an artifact related to the launch of Landsat 7 in 1999 with slightly deviating spectral properties. This finding supports the idea that spectral differences due to differing environmental conditions, rather than spectral changes caused by forest decline, enable us to at least partly distinguish declining from non-declining areas in our semi-arid study area.

Environmental variables with a potential influence on forest decline

Earlier studies showed that topography is occasionally related to forest decline in the Zagros forests (Goodarzi et al., 2016). Adding topographic information to the remote sensing predictors may be one way to improve the detection of declining areas. However, we utilized elevation, slope, and aspect data derived from a 30-meter DEM (see Hawker et al., 2022 for more data details) for one of the RF classifications (1986–2021 July 2001to August), and the results were not significantly improved (Appendix 1).

Differences in regional climatic conditions could be further important variables affecting forest decline. Multiple studies revealed an increase in forest decline worldwide associated with global warming and a more pronounced drought. Moreno-Fernández et al. (2019) found that warmer conditions are directly linked to worse vegetation health in open, evergreen oak woodlands. Ahmadi et al. (2014) showed that climatic variables such as yearly precipitation, temperature, and moisture were strongly correlated with the risk of forest decline within Zagros forests. Similarly, Kooh Soltani et al. (2018) showed that these three elements were most crucial in the emergence of forest decline. Hence, one way to further improve the detection of declining oak stands may be to integrate climatic data. However, it is not fully clear whether our study area is large enough to span a sufficiently large climate gradient (independent from topography) to improve our results. The highly clustered occurrence pattern of declining areas in our reference dataset (Figure 4) suggests that the use of standard climatic datasets such as Bioclim at a spatial resolution of 1 km may lead to notably improved accuracy (Appendix 1), which is, however, mostly due to the spatial autocorrelation in our reference data. Thinning out the reference dataset could be accomplished but would likely lead to a too-small dataset to draw meaningful conclusions, whereas field inventory of larger sample sizes would entail a high level of logistics and financial means. Uncrewed aerial vehicles may contribute to making the corresponding effort more efficient (Latifi, 2023).

The influence of different soils and co-occurring vegetation communities on the observed decline processes remains unclear and requires further investigation. These variables could potentially lead to changes in the observed spectral signal. While we currently have limited information on these two variables for our study region, we plan to explore their influence in future work.

Technical challenges of detecting forest decline in semi-arid areas

Future work on detecting and mapping forest decline in arid and semi-arid areas could also benefit from technical improvements addressing some of the challenges that may have impacted the results of the present study. For example, in (semi-)arid regions, detecting subtle decline processes is complicated by mixed pixels with bright soil backgrounds, as noted by Wang et al. (2022) and Maier et al. (2022). To address these challenges, we utilized blue bands and VIs that incorporated the soil background, such as GLI and ARVI. However, it should be noted that even with this approach, the correction for the soil background is not entirely accurate, and more sophisticated approaches may be needed.

A further, often-neglected influential factor on the spectral signal is the influence of the cast shadows of trees. Depending on spatial resolution and acquisition geometries, tree shadows can negatively impact the quality of RS data and negatively affect vegetation monitoring through mixed pixels (Wang et al., 2022), particularly in low-density forests with potentially large shadow fractions. To address this issue, Larsen and Salberg (2010) proposed tree shadow removal approaches. Although not the focus of this study, developing an algorithm to quantify and correct the influence of shadows could be useful for monitoring Zagros’ semi-arid forests via RS.

Other technical challenges that could have affected the identification and mapping of declining forest areas include the lack of correction factors to calculate VIs across different Landsat sensors. Landsat 4–7 and Landsat 8 have slight spectral band offsets, which can result in small changes in the spectral signal that are not related to underlying changes in the forest. Additionally, some spectral bands used in this study (red, NIR, and SWIR) are narrower for Landsat 8 than for Landsat 4–7. These challenges should be considered in future work to improve the accuracy of forest decline detection and mapping.

In an ideal scenario, satellite images would be acquired using identical technology and acquisition conditions, but this is not possible to achieve retrospectively. To address the differences between images in a time series, correction factors could be applied. However, in regions with limited information on the stability of the spectral signal over the last few decades, reference data or prior knowledge of spectrally stable areas would be necessary to successfully implement such corrections. Unfortunately, this was not feasible for our study due to the lack of corresponding information.

Conclusion

We evaluated three approaches in order to map declining forest areas in a semi-arid region using the Landsat time series from 1986 to 2021. These approaches were RF, anomaly detection, and Sen’s slope analysis. RF outperformed the other two methods. Sen’s slope and anomaly approaches produced unsatisfactory results, which did not result in maps that could serve as valuable input to forest managers. We found that VIs, including NIR and particularly SWIR, were crucial for all methods. Interestingly, the RF classification was not sensitive to the Landsat data’s acquisition time periods, suggesting that the ability to distinguish declining from non-declining forest areas was likely due to general environmental differences rather than due to spectral changes caused by the decline process. Further research is required to identify the environmental factors that contribute to these variances. In conclusion, mapping forest decline in semi-arid sparse cover forests using medium-resolution Landsat data is challenging, and technical improvements are needed to address issues with bright soil backgrounds and tree-cast shadows.

Authorship statement

Elham Shafeian: Conceptualization, Methodology, Fieldworks, Software, Formal analysis, Investigation, Visualization, Writing – original draft, Writing – review, and editing.

Fabian Ewald Fassnacht: Conceptualization, Methodology, Supervision, Software, Formal analysis, Visualization, Writing – review, and editing.

Hooman Latifi: Conceptualization, Fieldworks, Supervision, Writing – review and editing.

Acknowledgments

The German Academic Exchange Service (DAAD) is acknowledged by the first author for granting a Ph.D. scholarship. We would like to acknowledge the support received from the KIT-Publication Fund of the Karlsruhe Institute of Technology. Additionally, the Research Institute of Forests and Rangelands of Iran, the Agricultural and Natural Resource Research Center, and the Provincial Office of Natural Resources and Watershed Management of the Chaharmahal and Bakhtiari are appreciated for their aid in organizing site visits and assisting with fieldwork. Special thanks to Prof. Yaghoub Iranmanesh and the local foresters in the province for their valuable contributions. This study was conducted within the FORZA project at the Karlsruhe Institute of Technology and the K. N. Toosi University of Technology and was co-funded by the Iranian Ministry of Science, Research and Technology and the German National Space Agency on behalf of the German Federal Ministry of Education.

Disclosure statement

No known competing financial interests or personal relationships that could influence the study were declared by the authors.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/22797254.2023.2260549.

Additional information

Funding

This work was supported by the Deutscher Akademischer Austauschdienst; KIT-publication fund (library).