Robust Trend Analysis in Environmental Remote Sensing: A Case Study of Cork Oak Forest Decline

Abstract

We introduce a novel methodological framework for robust trend analysis (RTA) using remote sensing data to enhance the accuracy and reliability of detecting significant environmental trends. Our approach sequentially integrates the Theil–Sen (TS) slope estimator, the Contextual Mann–Kendall (CMK) test, and the false discovery rate (FDR) control. This comprehensive method addresses common challenges in trend analysis, such as handling small, noisy datasets with outliers and issues related to spatial autocorrelation, cross-correlation, and multiple testing. We applied this RTA workflow to study tree cover trends in Los Alcornocales Natural Park (Southern Spain), Europe’s largest cork oak forest, analysing interannual changes in tree cover from 2000 to 2022 using Terra MODIS MOD44B data. Our results reveal that the TS estimator provides a robust measure of trend direction and magnitude, but its effectiveness is dramatically enhanced when combined with the CMK test. This combination highlights significant trends and effectively corrects for spatial autocorrelation and cross-correlation, ensuring that genuine environmental signals are distinguished from statistical noise. Unlike previous workflows, our approach incorporates the FDR control, which successfully filtered out 29.6% of false discoveries in the case study, resulting in a more stringent assessment of true environmental trends captured by multi-temporal remotely sensed data. In the case study, we found that approximately one-third of the area exhibits significant and statistically robust declines in tree cover, with these declines being geographically clustered. Importantly, these trends correspond with relevant changes in tree cover, emphasising the ability of RTA to detect relevant environmental changes. Overall, our findings underscore the crucial importance of combining these methods, as their synergy is essential for accurately identifying and confirming robust environmental trends. The proposed RTA framework has significant implications for environmental monitoring, modelling, and management.

1. Introduction

1.1. Trend Analysis Using Remote Sensing

Remote sensing provides a consistent and repeatable method for capturing data at regular intervals, making it an ideal tool for trend analysis in ecological and environmental studies [1]. In this sense, trend analysis in remote sensing involves monitoring and quantifying changes in environmental variables over time using satellite imagery [2]. This approach is crucial for detecting and understanding long-term and short-term patterns, such as vegetation cover changes, afforestation, and deforestation [3].

Remote sensing time series can vary in temporal resolution, including daily, weekly, monthly, and annual inputs [4]. While monthly data are valuable for capturing phenological and seasonal dynamics [5], yearly data are more suitable for assessing interannual variations, mainly when focusing on monotonic trends [6]. Monotonic trends, representing consistent increases or decreases over time without reversal, are crucial to understanding interannual changes over time [7]. This approach allows for more consistent measurement of parameters of interest over interannual and decadal periods, minimising the effects of short-term variability [8]. Additionally, it provides a more robust framework for analysing persistent environmental patterns when proper methods are used.

The Theil–Sen (TS) slope estimator is widely used to estimate the slopes and magnitude of monotonic trends [9,10]. This non-parametric method is particularly suited for environmental data, which often contain noise and outliers [11]. Less robust alternatives, such as ordinary least squares (OLS) linear regression [12], can be more sensitive to outliers, leading to potentially skewed results [13]. While OLS is commonly used for its simplicity, its vulnerability to extreme values makes TS a preferred choice for more accurate and reliable trend estimation in remote sensing applications [14].

The statistical significance of the estimated trends is typically assessed using tests such as the Mann–Kendall (MK) test [15,16]. Alternatives to the MK test include Spearman’s rank correlation test [17], which is non-parametric and evaluates the correlation between ranks. However, it is less commonly used in trend analysis due to its broader focus on rank correlation rather than trend detection. However, neither the MK test nor Spearman’s test accounts for spatial correlation or cross-correlation in spatiotemporal data, which can affect the reliability of trend detection in geographically distributed datasets [18].

One of the critical challenges in remote sensing trend analysis, particularly in statistical analysis of time series, is dealing with small sample sizes. In such cases, the robustness of methods like the TS slope estimator becomes crucial, as it reduces the influence of outliers and provides a reliable measure of trend magnitude. The TS technique is robust to outliers, enabling it to ignore extreme values without affecting the overall slope estimation. It can reject up to approximately 29% of the sample size as outliers (known as the breakdown bound: $1-{\displaystyle \frac{1}{\sqrt{2}}}\approx 29.3\%$) without compromising the results, as described by [19]. Additionally, the MK test is well-suited for small datasets as it is a non-parametric method that does not assume a specific data distribution, enhancing the reliability of trend detection in limited temporal series. These two methods are often employed together to provide a more robust and comprehensive analysis of trends in remote sensing time series [20].

Spatial autocorrelation is critical in remote sensing data analysis [21], as nearby locations often exhibit similar trends. While previous methods have incorporated spatial context into the analysis, such as the regional methods for trend detection [22,23], they usually need help with the complex interactions between spatial and cross-correlation, leading to potentially misleading trend detections. The Contextual Mann–Kendall (CMK) test addresses these limitations by accounting for both spatial and cross-correlation, enhancing the robustness of trend detection [24].

Another major challenge in spatiotemporal trend analysis is the potential for Type I errors [25], which occur when a false positive trend is detected [26]. This issue is exacerbated by the large number of statistical tests conducted when analysing gridded remote sensing data [27], leading to an increased risk of incorrectly identifying significant trends. To address this problem, corrections for multiple testing are applied [28], such as the family-wise error rate (FWER) control methods. However, these methods tend to be overly conservative. In contrast, the false discovery rate (FDR) correction is more suitable for controlling the proportion of false positives in large datasets. It is important to note that using such corrections is uncommon in remote sensing studies [29], despite their proven effectiveness in improving the reliability of results.

Although spatial autocorrelation, small sample sizes, and multiple testing have been addressed individually in previous environmental remote sensing studies, an integrated approach to tackle these challenges jointly has yet to be explored, particularly in tree cover analysis. Notably, applying multiple testing corrections is infrequent in environmental remote sensing [29], yet it is necessary to ensure the reliability of trend detection [30]. Therefore, there remains a significant gap in research for a comprehensive framework that simultaneously addresses these issues, providing a more robust method for detecting significant trends in tree cover dynamics over time.

1.2. The Decline in Mediterranean Cork Oak Forests

Quercus suber L., commonly known as cork oak, is a medium-sized broadleaved tree renowned for its slow growth and remarkable longevity [31]. Cork oak is commonly found in closed forests or open woodlands, often as the dominant tree species or alongside other Mediterranean trees [32,33]. This species, renowned for its corky bark, can be harvested in late spring or early summer every 9–12 years, promoting a sustainable cork-based economy [34,35].

Cork oak forests and open woodlands cover approximately 2.2 million hectares worldwide, with the vast majority situated in the western areas of the Mediterranean ecoregion [36,37,38,39]. These cork oak landscapes provide significant ecological, economic, and social importance, highlighting their biodiversity and sustainable forest production in these areas [40], mainly in the Iberian Peninsula, where its management is historical [41,42,43].

As a result of its evolutionary history [44], the cork oak has tended to form dense forest formations, further favoured by anthropic management, even when facing the challenge of summer drought in Mediterranean environments. However, over the past few decades, Mediterranean oak forests have become increasingly vulnerable, raising concerns about their ongoing decline [45,46]. Cork oak decline is widespread and may result from mismanagement [47,48], abiotic factors such as temperature increases and droughts [49,50,51], and biotic factors such as pest and pathogen attacks, which can trigger the decline [52,53,54,55,56]. Additionally, forest degradation impacts root dynamics, particularly the lifespan of fine roots, a process influenced by the species’ evolutionary history [44].

Since 1990, the Iberian Peninsula has reported rapid death and decline in cork oak and other oak species [57]. While some authors have proposed spatially explicit models to address the causes of decline at the plot scale [58,59], a significant research gap exists in addressing the tree cover loss using an extensive spatiotemporally explicit approach.

This approach maps both the geographical distribution and temporal progression of the decline, facilitating the identification of the most affected areas to offer potential solutions adapted to the specific conditions of each location, using spatiotemporally explicit data derived from remote sensing data. It applies robust time series analysis to estimate the interannual trends of tree cover in the largest cork oak-dominated mixed Mediterranean forest in Europe: Los Alcornocales Natural Park. Since cork oak decline has been documented at the plot scale in the same area (more details in the “Study Area” section), we hypothesised that trends would be identified, indicating widespread tree cover loss in the study area.

1.3. Research Objectives

Our research aim is primarily to introduce a robust trend analysis (RTA) that integrates advanced statistical techniques to address trends in spatiotemporal grid data, including the TS estimator, the CMK test for accounting for spatial and cross-correlation, and the FDR control to address multiple testing issues. Importantly, none of the leading software tools for trend analysis provide tools for controlling multiple testing, as evidenced by modules such as Kendall and Earth Trend Modeller in Terrset [60], LandTrendr Analysis in ArcGIS Pro 3.3 [46], and R packages such as Trends, Kendall, and Modifiedmk [61,62,63]. By integrating these methods, the study aims to overcome the limitations of traditional approaches, particularly in managing spatial data and enhancing the reliability of trend detection.

Additionally, we aim to provide reliable empirical evidence on the spatiotemporal trends of tree cover in the largest cork oak forest in Europe and quantify its extent in a statistically rigorous way. Specifically, we aim to apply the proposed RTA to assess the extent and patterns of previously reported cork oak decline in Los Alcornocales Natural Park and compare the results with those derived from previous approaches.

2. Materials and Methods

2.1. Study Area

Los Alcornocales Natural Park (~174,000 hectares) is a protected area southwest of the Iberian Peninsula, straddling the provinces of Málaga and Cádiz in Andalusia, Spain (Figure 1). The Spanish expression “Los Alcornocales” refers to the “cork oak forests”. Indeed, this protected area, near the Strait of Gibraltar, is home to Europe’s largest and most well-conserved cork oak forest and one of the most important in the world [64], situated within a significant biogeographical crossroads in the Western Palearctic [65].

The study area has a relatively humid Mediterranean climate with oceanic influence [66]. Winters are temperate with scarce frosts and rare snowfalls, while summers are warm. Precipitation averages around 1000 mm, with frequent fog near the Strait of Gibraltar. Mixed sclerophyll forests of Quercus suber dominate the tree forest cover in the Alcornocales Natural Park and coexist with the deciduous, shade-tolerant Quercus canariensis, forming closed forests; in contrast, it forms mixed open woodlands with the evergreen drought-tolerant Olea europaea var. sylvestris. In the drier areas of the natural park, the understorey of woodlands is dominated by shrub species such as Pistacia lentiscus L. In the wetter areas, the understorey of closed forests is more diverse, featuring species such as Arbutus unedo L., Phillyrea latifolia L., Erica arborea L., and Erica scoparia [67]. These forests have been relatively well conserved since their designation as a natural park in 1989, and their tree forest covers have remained relatively stable in the second half of the 20th century [68].

The Regional Environmental Administration of Andalusia has been monitoring the decline in oak forests and open woodlands in the region [69]. This decline is attributed to various factors, including abiotic stressors like drought and biotic factors like insect infestations and pathogenic infections [58,70,71]. One of the most significant impacts of cork oak forest decline on ecosystems is the ongoing defoliation and dieback [72,73], which can have devastating consequences for biodiversity and ecosystem services in Los Alcornocales Natural Park (Figure 2).

2.2. Research Workflow

This research is based on the yearly product of MOD44B Version 6.1 Vegetation Continuous Fields (VCF) [74]. The MODIStsp package was used in our workflow to download and create a raster time series derived from MODIS Land Products. It facilitates several preprocessing steps, such as downloading, mosaicking, reprojection, and resizing MODIS products within a selected period and area [75]. The MOD44B data product includes many layers, and we used the percent tree cover layer, which covers the entire temporal range available from 2000 to 2022.

In the study area, we assessed the interannual trends in tree cover (%) from 2000 to 2022. To achieve this, we processed 23 images representing the annual tree cover at a spatial resolution of 250 m using the ETRS89 H30 projected coordinate system (EPSG: 25830). Each image has dimensions of 195 columns and 330 rows, totalling 64,350 pixels; of these, 27,828 correspond to the study area, which covers around 173,925 hectares. In operational terms, this defines our study area. The time series was converted into a multidimensional GIS object consisting of 23 layers, each representing the tree cover percentage for a specific year from 2000 to 2022 at the pixel level. This approach captures the interannual dynamics in tree cover across the entire study area, allowing for a detailed analysis of changes over time.

We conducted a robust trend analysis (RTA) using a three-step approach to control trends’ magnitude, direction, and significance while addressing spatial autocorrelation, cross-correlation, and multiple testing (Figure 3). In the Level 1 analysis, we estimated the magnitude (annual rate, expressed as a percentage %) and direction (+/−) of interannual trends using the TS slope estimator [9,10]. Moving to the Level 2 analysis, we evaluated the significance of the interannual trends identified in Level 1 using the CMK test [76], only revealing statistically significant trends (alpha = 0.05). The CMK test modifies the traditional MK test [15,16] by integrating contextual information from adjacent pixels. In the subsequent Level 3 analysis, we addressed the potential for Type I errors (false positives) by applying FDR control to adjust the p-values [77] to the output derived from the Level 2 analysis.

The TS slope results were visually represented on maps and plots of the three-level trend analysis. Additionally, we categorised the TS slopes to facilitate their clarity and interpretation: (1) Sharp Decrease (>−0.5): Indicates a notably negative annual rate of change; (2) Moderate Decrease (−0.5 to −0.1): Indicates a moderate annual rate of change; (3) Stable (−0.1 to 0.1): Implies a negligible annual rate of change; (4) Moderate Increase (0.1 to 0.5): Indicates a moderate positive annual rate of change; (5) Sharp Increase (>0.5): Indicates a notably positive annual rate of change.

As a subsequent step to the RTA, Moran’s I was used to assess the spatial autocorrelation of the trends, focusing exclusively on the significant trends controlled by the FDR. The index ranges from −1 (indicating perfect dispersion) to +1 (indicating perfect clustering) [78,79].

We used the R environment [80] and the Terrset system [81] regarding the software tools and processed the gridded data using the Earth Trend Modeler module and the terra package [82]. Additionally, we used the Kendall module from TerrSet 2020 Geospatial Monitoring and Modelling Software to apply both the TS estimator and the CMK test [81], and the stats package in R to calculate the FDR correction. Finally, we used the spdep package in R to calculate and plot Moran’s I [83], enabling us to assess the spatial autocorrelation of our data. This allowed us to quantify and graphically represent spatial dependence, which is crucial for understanding the clustering patterns in our research.

2.3. Statistical Foundations of the Robust Trend Analysis (RTA)

2.3.1. Robust Trend Analysis Considering Spatial and Cross-Correlation

The Theil–Sen (TS) slope estimator calculates the median slope from all possible pairwise comparisons of observation values, resulting in a total of N slopes, where N = ${\displaystyle \frac{n\left(n-1\right)}{2}}$. This non-parametric technique for estimating the magnitude of trends in time series data was introduced by Theil (1950) [9] and later refined by Sen (1968) [10]. The equations used to estimate the TS slope and intercept are as follows:

$$\mathrm{TS} \mathrm{Slope}=\mathrm{Median}\left({\displaystyle \frac{X_j-X_i}{t_j-t_i}}\right)$$

$$\mathrm{Intercept}=X_i-\mathrm{MedianSlope} \ast t_i,\hspace{1em}\mathrm{for} i=1\dots n$$

where:

$X_j-X_i$ is the difference in the observation values between two points in time; $t_j-t_i$ is the difference in time between two observations, with $t_j$ being later than $t_i$; and, here, N represents the total number of non-zero differences $t_j-t_i$ for all pairs $\left(i,j\right)$ such that $1 \le i<j \le n$. Variables $i$ and $j$ are indices that iterate through the data series, ensuring that $i$ is always less than $j$ so that all possible pairs of observations are considered without repetition or reversing the order. This is necessary to calculate the slope differences between each pair of observations in the TS method.

This approach offers several advantages. The TS technique is robust to outliers, enabling it to ignore extreme values without affecting the overall slope estimation. It can reject up to approximately 29% of the sample size as outliers (known as the breakdown bound: $1-{\displaystyle \frac{1}{\sqrt{2}}}\approx 29.3\%$) without compromising the result, as described by [19].

A widely used method is the MK test to assess the significance of the TS slope [15,16]. Like the TS approach, the MK test evaluates the slopes between all possible pairs of data points. In the MK test, the data are ordered chronologically, with each data point serving as a reference for the subsequent data points in time. The statistic S, introduced by Kendall (1975) [16], is defined as follows:

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^n\mathrm{sign}\left(x_i-x_j\right)$$

where:

$$\mathrm{sign}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}\right)=\left\{\begin{array}{c}1\\ 0\\ -1\end{array}\begin{array}{c}if x_j-x_i>0\\ if x_j-x_i=0\\ if x_j-x_i<0\end{array}\right.$$

and:

$S$ is the MK test statistic that sums the signs of all differences between pairs of observations; $n$ represents the number of observations in the time series; and $x_i-x_j$ are the values of the time series at time points $i$ and $j$, respectively. In the MK test, the $\mathrm{sign}$ function is used to evaluate whether the difference between two values in a time series is positive, negative, or zero.

The variance ${\sigma }^2$ of $S$ is used to normalise the statistic and is calculated as follows:

$${\mathsf{\sigma }}^2={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)}{18}}-{\displaystyle \sum }_{t=1}^g{\displaystyle \frac{t\left(t-1\right)\left(2t+5\right)}{18}}$$

where:

${\sigma }^2$ is the variance of the statistic $S$; $g$ is the number of tied groups in the dataset; and $t$ represents the size of each tied group. The first term $({\sigma }^2={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)}{18}}$) accounts for the total variance, while the second term $({{\displaystyle \sum }}_{t=1}^g{\displaystyle \frac{t\left(t-1\right)\left(2t+5\right)}{18}})$ corrects for ties in the data.

The Z statistic is used to determine the significance of the trend as follows:

$$\mathrm{Z} =\left\{\begin{array}{cc}{\displaystyle \frac{\mathrm{S}-1}{\mathsf{\sigma }}}& if S>0\hfill \\ 0& if S=0\hfill \\ {\displaystyle \frac{\mathrm{S}+1}{\mathsf{\sigma }}}& if S<0\hfill \end{array}\right.$$

and:

$$p=2\times \mathsf{\Phi }\left(-\left|Z\right|\right)$$

where:

$Z$ is the standardised test statistic, which follows a normal distribution; $p$ is the p-value, the probability that measures the evidence against the null hypothesis; and $\mathsf{\Phi }$ represents the cumulative distribution function (CDF) of the standard normal distribution. The value of $Z$ indicates the strength and direction of the trend. If $Z$ is positive and significant, there is an increasing trend; if negative and significant, there is a decreasing trend.

We used the Contextual Mann Kendall (CMK) test, an enhanced version of the traditional MK test, to assess the statistical significance of interannual trends [24]. Unlike the original MK test [15,16], the CMK test sequentially incorporates contextual information from eight first-order neighbours to correct for spatial correlation and addresses cross-correlation. The equations used to estimate CMK significance are as follows:

$$Z_m={\displaystyle \frac{{\overline{S}}_m-E\left({\overline{S}}_m\right)}{\frac{\mathsf{\sigma }}{\sqrt{m}}}}$$

and:

$${\overline{S}}_m={\displaystyle \frac{1}{m}}{\displaystyle \sum }_{j=1}^m{S}_j$$

and:

$$Var\left({\overline{S}}_m\right)={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)}{18m}}={\displaystyle \frac{{\mathsf{\sigma }}^2}{m}}$$

where:

$S_j$ is Kendall’s for the j-th neighbour; $m$ = 9 pixels, which includes eight neighbours with a central pixel; and E(S) and σ are the mean (expected value) and standard deviation, respectively.

Adjustment for cross-correlation is performed by introducing a covariance term in the calculation of variance, as follows:

$$Var\left({\overline{S}}_m\right)={\displaystyle \frac{1}{m^2}}\left[{\displaystyle \sum }_{j=1}^{m}Var\left(S_j\right)+2{\displaystyle \sum }_{j=1}^m{\displaystyle \sum }_{l=1}^{m-j}Cov\left(S_j,S_{j+l}\right)\right]$$

According to Neeti and Ronald Eastman (2014) [24], the CMK method reduces the detection of spurious trends while increasing confidence in the presence of consistent ones. Furthermore, the application of the TS slope offers the distinct advantage of filtering out inter-annual variability shorter than 0.29 times the length of the series. Together, these techniques are non-parametric and robust against the influence of outliers.

2.3.2. False Discovery Rate Control

The false discovery rate (FDR) control procedure aims to control the proportion of spurious significances (false discoveries) concerning the total number of discoveries (i.e., the total number of tests declared significant) in a multiple-testing context.

Benjamini and Hochberg (BH) introduced the original FDR approach, which is the expected number of false discoveries within all discoveries [77].

The FDR is formally defined as follows:

$$\mathrm{FDR}=E\left[{\displaystyle \frac{V}{R}}\right]$$

where:

$E$ denotes the expected value; $V$ is the number of false positives (also known as false discoveries); and $R$ is the total number of rejected hypotheses (referred to as discoveries).

The BH procedure is applied as follows:

1.

Order the p-values of all the hypothesis tests in ascending order:

$$p_{\left(1\right)},p_{\left(2\right)},\dots ,p_{\left(m\right)}$$

where:

$p_{\left(1\right)}$ represents the smallest p-value from the set of hypothesis tests; $p_{\left(2\right)}$ represents the second-smallest p-value; $p_{\left(m\right)}$ represents the largest p-value; and $m$ is the total number of hypothesis tests performed.

2.

Determine the critical value k as the largest i such that:

$$p_{\left(i\right)}\le {\displaystyle \frac{i}{m}}\cdot \alpha$$

where:

$p_{\left(i\right)}$ denotes the smallest p-value from the set of hypothesis tests; $i$ is the rank of the p-value when ordered from smallest to largest; $i$ is the total number of hypothesis tests being performed; and α is the desired FDR level.

3.

Reject all null hypotheses:

$$H_{\left(1\right)},H_{\left(2\right)},\dots ,H_{\left(k\right)}$$

where:

$H_{\left(1\right)}$ represents the hypothesis corresponding to the smallest p-value; $H_{\left(2\right)}$ represents the hypothesis corresponding to the second-smallest p-value; $H_{\left(k\right)}$ represents the hypothesis corresponding to the k-smallest p-value; and k is the largest rank for which the p-value meets the BH criterion.

3. Results

3.1. The Effect of Robustness in Trend Analysis

Overall, according to the three-level trend analysis conducted in line with our workflow (Figure 4), we observed a prevalence of trends skewed towards negative values. Negative trends were prevalent when assessing their statistical significance and further accentuated when adjusting the p-values to control the false discovery rate. In contrast, positive trends, particularly those with slopes close to zero, were notably minor and almost entirely disappeared after the statistical significance tests were applied.

This apparent shift in trend direction can be attributed to the integration of the CMK test, which accounts for both spatial autocorrelation and cross-correlation. While the TS estimator is robust against outliers, it does not consider the spatial or cross-correlational relationships between data pixels. Applying the CMK test enhances the robustness of the analysis by recognising that neighbouring areas may exhibit similar trends due to underlying spatial dependencies and cross-correlations. In addition to capturing these spatial dynamics, this method discounts low slopes, likely filtering out less significant trends that could be driven by noise rather than true underlying patterns. Consequently, when these factors are incorporated, the analysis tends to reveal a more pronounced directional bias, particularly towards negative trends, which may reflect decline patterns that the TS estimator alone might not fully discriminate. This comprehensive approach underscores the importance of accounting for spatial structure, cross-correlation, and the significance of trend magnitudes to avoid overlooking or misinterpreting key spatially driven patterns.

The FDR control further refines the trend analysis by building upon the initial filtering completed by the CMK test. While the CMK test already accounts for spatial autocorrelation and cross-correlation, effectively filtering out trends that are not statistically significant, the FDR correction adds a layer of robustness. It specifically targets the control of false discoveries across multiple tests, ensuring that even fewer trends around zero, which might have been considered marginally significant, are identified as meaningful. Consequently, the analysis produces a cleaner and more focused graphical representation, with a sharper distinction between truly significant trends and those that could be attributed to noise or statistical artefacts. This results in a more stringent identification of robust patterns, particularly those indicative of broader decline trends.

Figure 5 illustrates the effect of applying the FDR control using the BH procedure on the p-values obtained from the CMK test. It demonstrates a transition where p-values shift from statistically significant (blue) to non-significant (red). Following the FDR control, this shift indicates that many p-values initially considered significant are now identified as non-significant, as will be detailed further below (Figure 6).

This finding supports the notion that the FDR control enhances the robustness of the CMK test by filtering out trends that might otherwise appear significant due to mere random variation. The presence of red dots representing p-values from the CMK test that no longer meet the significance criteria after FDR adjustment underscores the reduction in potential false positives. This reduction is crucial for ensuring that the trends identified as significant are statistically robust and not merely the result of multiple testing artefacts. By eliminating these less significant trends, the FDR control sharpens the distinction between significant and non-significant trends, resulting in a more robust and rigorous statistical analysis. Finally, as depicted in Figure 5, the FDR control facilitates a more focused examination of robust trends, particularly those indicative of decline patterns, by filtering out marginally significant trends. This outcome aligns with previous findings, where integrating spatial autocorrelation, cross-correlation, and FDR control collectively enhances the identification of the most reliable trends.

Figure 6 shows histograms comparing the distribution of original and adjusted p-values (FDR) derived from the CMK trend test across 27,828 pixels within our masked data analysis area. The left histogram details the original p-values, revealing a pronounced skew towards lower values and marked by a green dashed line at 8218 (denoting the first significant bin frequency in the second histogram). In contrast, the right histogram displays adjusted p-values post-FDR control, demonstrating a moderated frequency of low p-values, effectively curbing false discoveries. Initially, we observed that 11,677 p-values met the significance criterion. However, upon applying the false discovery rate (FDR) adjustment, which effectively controls for false discoveries, only 8218 p-values retained their significance. This adjustment revealed approximately 3459 p-values initially identified as significant that did not meet the stringent FDR criteria. Consequently, this suggests that about 29.6% of the initially significant p-values were false positives within the framework of our FDR control.

3.2. Tree Cover Loss in the Largest Cork Oak Forest in Europe

In this section, we present the results of the case study, which followed a three-step approach for robust trend analysis (RTA). The outcomes are structured across three levels and are presented sequentially to illustrate the comparison between each level. Figure 7 and Figure 8 in the manuscript visually represent the key findings of this RTA, highlighting the interannual trends in tree cover in the study area.

In the Level 1 analysis, our findings indicate that, of the total study area, 75% exhibited a negative trend in the TS slope. Further categorisation based on the TS slope revealed that sharp decreases accounted for an area equivalent to 42,831 hectares (24.63% of the study area). In comparison, moderate decreases accounted for an area equivalent to 65,700 hectares (37.77%). Stable areas, representing trends between −0.1 and 0.1 (including both slightly positive and slightly negative values), comprised 22.99% of the study area, totalling 39,993 hectares. Likewise, moderate increases constituted 13.45% of the study area (23,393 hectares), while sharp increases constituted 1.15% (2006 hectares).

In the Level 2 analysis, 41.96% of the study area exhibited statistically significant trends (α = 0.05) based on the CMK test and the TS slope estimator. Further categorisation based on the TS slope revealed that sharp decreases accounted for an area equivalent to 40,056 hectares (23.03% of the total study area and 54.89% of the area with statistically significant trends). In comparison, moderate decreases accounted for an area equivalent to 25,900 hectares (14.89% of the total study area and 35.49% of the area with statistically significant trends). Stable areas totalled 881 hectares (0.51% of the total study area and 1.21% of the area with statistically significant trends); moderate increases constituted 5000 hectares (2.87% of the total study area and 6.85% of the area with statistically significant trends); and sharp increases, totalling 1143 hectares (0.66% of the total study area and 1.57% of the area with statistically significant trends), represented a smaller portion.

In the Level 3 analysis, 29.53% of the study area exhibited statistically significant trends (α = 0.05) based on the CMK test and the TS slope estimator, with the significance of the CMK test adjusted using the FDR. Further categorisation based on the TS slope revealed that sharp decreases accounted for an area equivalent to 34,818 hectares (20.02% of the total study area and 67.79% of the area with statistically significant trends). In comparison, moderate decreases accounted for an area equivalent to 12,937 hectares (7.44% of the total study area and 25.19% of the area with statistically significant trends). Stable areas totalled 212 hectares (0.12% of the total study area and 0.41% of the area with statistically significant trends); moderate increases constituted 2688 hectares (1.53% of the total study area and 5.20% of the area with statistically significant trends); and sharp increases, totalling 725 hectares (0.42% of the total study area and 1.41% of the area with statistically significant trends), represented a smaller portion.

Figure 9 represents the relevant changes in tree cover within the identified significant trends. The analysis of actual changes in forest cover from 2000 to 2022 shows that 88.86% of the pixels with significant trends experienced a decrease, with a mean loss in tree cover of −18.24%. Of these pixels, 69.34% had less than 10% losses, highlighting a considerable decline in forest cover across the study area. This indicates that, over the 23 years, the areas with significant trends identified by the RTA correspond closely to those where tree cover losses were most relevant.

Figure 10 shows the spatial autocorrelation of significant TS slopes filtered with the CMK test and FDR control using Moran’s I. The high Moran’s I value of 0.812 indicates strong positive spatial autocorrelation, meaning that areas with similar rates of change in forest cover—increasing or decreasing—are likely to be geographically clustered. This strong clustering suggests that the processes driving changes in forest cover are not uniformly distributed across the landscape. Instead, certain areas are experiencing similar trends. Focusing on the areas where trends are most prevalent, most points fall into regions indicating significant decreases in forest cover, as seen in the lower left quadrant of the plot. These points represent geographic areas where the local trend and the neighbouring areas consistently decline in forest cover. This clustering of negative trends suggests that deforestation or forest degradation is occurring in a spatially concentrated manner.

4. Discussion

We introduce a robust methodological framework for trend analysis using remote sensing data, which effectively handles outliers in small, noisy datasets and challenges such as spatial autocorrelation and cross-correlation, as well as false discoveries. Applied to Los Alcornocales Natural Park, Europe’s largest cork oak forest, it reveals significant declines in tree cover from 2000 to 2022, underscoring the severity of oak decline in this unique protected area. The combined use of the CMK test and the FDR control further enhances the reliability of trend detection by filtering out spurious trends and isolating genuine patterns of decline.

While the TS slope estimator is highly robust against outliers, it alone does not differentiate between significant and non-significant trends. By integrating the CMK test, we enhance the analysis by clearly distinguishing statistically significant trends and accounting for spatial autocorrelation [76]. This integration allows for a more nuanced understanding of the trends observed, ensuring that only those trends that are both robust and statistically significant are highlighted. Furthermore, applying the FDR control mitigates the risk of false positives arising from multiple testing. Ventura et al. (2004) [84] demonstrated that the original FDR–BH procedure maintains strong performance even when dealing with spatially correlated data, performing comparably to more complex methods specifically designed for such data while being more straightforward to implement. Our findings underscore the importance of combining these statistical tools, as this approach significantly enhances the reliability and validity of the results, ensuring that interpretations are based on significant patterns rather than data artefacts [85].

Given that cork oak decline has previously been documented at the local scale in the same area, as detailed in the “Study Area” section, our results confirm our research hypotheses by demonstrating significant interannual tree cover loss. For the same study area, different researchers have highlighted biotic and abiotic factors as the main drivers of this forest decline [58,86,87]. Gómez-Aparicio et al. (2012) [58] found strong empirical support for spatial concordance between the distribution and health status of trees of different species and the abundance of soil-borne pathogens in the soil. Ibáñez et al. (2014) [86,87] investigated the consistency of neighbourhood effects on tree seedling performance over several consecutive years with contrasting climatic conditions. They noted that the negative impacts on oak seedling establishment due to the dieback of Q. suber and loss of tree cover could increase as these systems become drier and warmer with ongoing climate change. The elevated spatial autocorrelation we observed reinforces previous findings of cork oak decline in the study area, suggesting that the combined biotic and abiotic stressors driving this decline are spatially clustered and consistently impact specific regions, potentially linked to the concentration of pathogens and the intensity of drought [88,89].

A significant decrease in tree cover in the study area is critical from a conservation perspective. Several authors have investigated the impacts and feedback of cork oak decline in the same study area [59,70,71]. Avila et al. (2017) studied the interactions between Q. suber and decline and mortality induced by the exotic soil-borne pathogen Phytophthora cinnamomi. They suggest that a lack of coevolutionary history between Q. suber and P. cinnamomi likely contributes to the high vulnerability of this tree species to the pathogen. Their findings indicate that while Q. suber can ecophysiologically respond to water stress at local and landscape scales, the collective response of defoliated trees may not suffice to counteract the stress induced by P. cinnamomi. Homet et al. (2019) [71] investigated the potential interactions between pathogens and climate change, revealing that their findings suggest a higher frequency of extreme rain events, saturating the soil and inducing temporary waterlogging, will particularly benefit P. cinnamomi infections, amplifying its density beyond any potential response capacity of susceptible hosts. Domínguez-Begines et al. (2019) [59] found that the decline in Q. suber significantly altered the functional composition of soil communities. Their findings underscore the critical role of Q. suber in maintaining soil biodiversity and ecosystem functionality. Additionally, their results suggest potential broader implications for ecosystem health, emphasising the importance of understanding and preserving the role of Q. suber in terrestrial ecosystem dynamics.

This is the first research study to address the interannual trends in tree cover in a cork oak-dominated mixed Mediterranean forest by combining remote sensing with robust trend analysis (RTA). Despite the robust methodological framework used in this study, certain limitations must be acknowledged. One of the primary challenges lies in the reliance on MODIS remote sensing data, which, while powerful, can be subject to resolution limitations and potential inaccuracies in detecting subtle changes in forest cover, particularly in mixed Mediterranean forests where canopy structure is complex [90]. While MODIS data and robust trend analysis provide a solid foundation, the inclusion of higher-resolution data and field sampling methods, such as airborne local remote sensing with sample plots [91], would significantly enhance the precision and validation of the detected trends. This approach could help better capture the fine-scale variability crucial for understanding the full scope of oak decline in the study area.

Future research should focus on integrating these higher-resolution data (Landsat, Sentinel) and incorporating field-based observations to validate and refine the detected trends. Additionally, exploring the role of specific biotic factors, such as the spatial distribution and dynamics of pathogens like Phytophthora cinnamomi, in relation to the observed patterns of decline could provide deeper insights into the mechanisms driving these changes [92]. Extending the analysis to consider the interactions between a drier climate and pathogen dynamics over longer temporal scales would also be valuable in predicting future trends and informing conservation strategies [93]. Collaborative efforts combining ecological modelling with advanced remote sensing techniques and field validation could offer a more comprehensive understanding of the broader implications of cork oak decline on ecosystem health and resilience.

While our methodology has proven effective in detecting robust trends, there are several areas where future research could further refine and enhance these techniques. One potential improvement could be the integration of machine learning algorithms to complement the statistical methods, strengthening the ability to identify complex patterns and interactions that traditional approaches may not fully capture. Additionally, future studies could explore the application of multi-temporal analyses that consider variations on different timescales [94]. Assessing the impact of a prewhitening time series, a less commonly used technique in remote sensing of vegetation, could also provide insights into its relevance for improving trend detection [95,96]. Furthermore, developing and applying more sophisticated FDR control methods in scenarios with strong spatial dependence would be beneficial to ensure the robustness of results in such contexts.

5. Conclusions

This research introduces a novel methodological framework for robust trend analysis (RTA) in environmental remote sensing. RTA effectively and jointly addresses key challenges such as spatial autocorrelation, cross-correlation, handling small, noisy datasets with outliers, and filtering out potential false discoveries.

By integrating the Theil–Sen (TS) slope estimator and the Contextual Mann–Kendall (CMK) test, we offer a more reliable tool for detecting significant environmental trends using remote sensing data. The subsequent application of the false discovery rate (FDR) control further refines the analysis by controlling for false positives, ensuring that only the most robust trends are retained. As a case study, we apply this methodology to assess the decline in forests in Los Alcornocales Natural Park, Europe’s largest cork oak forest.

Our analysis revealed significant and statistically robust decreases in tree cover across a third of the study area, even after filtering out 29.6% of false discoveries using the FDR control. This underlines the effectiveness of our approach in accurately capturing these trends, discarding a high proportion of spurious trends, which were considered significant when applying previous methodologies. Moreover, the trends identified correspond with relevant changes in tree cover, emphasising the importance of trend analyses that truly signify meaningful changes.

Furthermore, this new approach has potentially relevant implications for environmental monitoring, modelling, and management, as it allows us to more accurately and robustly assess areas undergoing significant changes and focus attention and, where appropriate, resources on those key areas. Going forward, targeted management practices, supported by robust data analysis and continuous monitoring, will be crucial for mitigating the impacts on these valuable ecosystems facing environmental crises.