Nothing Special   »   [go: up one dir, main page]
Next Article in Journal
Landslide Hazard Prediction Based on UAV Remote Sensing and Discrete Element Model Simulation—Case from the Zhuangguoyu Landslide in Northern China
Next Article in Special Issue
Intelligent Monitoring Applications of Landslide Disaster Knowledge Graphs Based on ChatGLM2
Previous Article in Journal
The Zenith Total Delay Combination of International GNSS Service Repro3 and the Analysis of Its Precision
Previous Article in Special Issue
Spatio-Temporal Knowledge Graph-Based Research on Agro-Meteorological Disaster Monitoring
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 16
Issue 20
10.3390/rs16203886
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Oliver Gutiérrez-Hernández
1,* and
Luis V. García
2
1
Department of Geography, University of Málaga (UMA), 29071 Málaga, Spain
2
Institute of Natural Resources and Agrobiology of Seville (IRNAS), Spanish National Research Council (CSIC), 41012 Seville, Spain
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(20), 3886; https://doi.org/10.3390/rs16203886
Submission received: 27 August 2024 / Revised: 7 October 2024 / Accepted: 11 October 2024 / Published: 19 October 2024
(This article belongs to the Special Issue Temporal Resolution, a Key Factor in Environmental Risk Assessment II - Integrating Data from Multiple Data Sources)
Browse Figures
Figure 1
<p>Study area: Los Alcornocales Natural Park (≈174,000 ha), a representative cork oak-dominated mixed Mediterranean forest in the Southwest Iberian Peninsula (Andalusia, Spain). Source: Author’s work. Data from the Department of the Environment of the Regional Government of Andalusia, specifically from the Andalusian Network of Environmental Information (REDIAM).</p> "> Figure 2
<p>Evidence of cork oak forest decline in Los Alcornocales Natural Park. Source: Lorena Gómez-Aparicio.</p> "> Figure 3
<p>Research workflow diagram for robust trend analysis (RTA). Source: Author’s work.</p> "> Figure 4
<p>Density of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: TS slope (dashed green line); combined TS slope and CMK significance (thin solid green line); combined TS slope and CMK significance with FDR control (thick solid green line). Axes details: X-axis—TS slope; Y-axis—spatial distribution density.</p> "> Figure 5
<p>Controlling the FDR for interannual tree cover trends (2000–2022) in Los Alcornocales Natural Park: FDR–BH <span class="html-italic">p</span>-value adjustment for the CMK test. Source: Author’s work. Note: The <span class="html-italic">x</span>-axis represents the rank (<span class="html-italic">k</span>) of 27,828 <span class="html-italic">p</span>-values sorted in ascending order, while the <span class="html-italic">y</span>-axis displays the observed <span class="html-italic">p</span>-values. The celestial blue line symbolises the threshold defined by the FDR, which is calculated using the following expression: <math display="inline"><semantics> <mrow> <mstyle scriptlevel="0" displaystyle="true"> <mfrac> <mrow> <mi>k</mi> <mo>·</mo> <mi>α</mi> </mrow> <mi>m</mi> </mfrac> </mstyle> </mrow> </semantics></math>. Herein, <span class="html-italic">k</span> denotes the rank, α is the significance level, and m is the total number of tests conducted. The part of the line coloured blue indicates the <span class="html-italic">p</span>-values that are deemed significant following the application of the FDR control, using the FDR–BH procedure. Conversely, the part of the line coloured red represents the <span class="html-italic">p</span>-values that are not considered significant according to this adjustment. Altogether, this bi-coloured line comprises 27,828 points (pixels), corresponding to the number of tests performed (<span class="html-italic">m</span>).</p> "> Figure 6
<p>Histogram of original and adjusted <span class="html-italic">p</span>-values (FDR). Both histograms share the same scale axes for easy comparison. Source: Author’s work. Note: The <span class="html-italic">x</span>-axis shows <span class="html-italic">p</span>-values ranging from 0 to 1, and the <span class="html-italic">y</span>-axis shows frequencies up to 12,000 (pixels). The first bin in each histogram represents pixels with significant trends (α = 0.05). The left histogram displays the original <span class="html-italic">p</span>-values, heavily skewed towards 0, with a green dashed line at 8218 (pixels), which serves as a reference for comparison with the adjusted <span class="html-italic">p</span>-values at the alpha level of 0.05. The right histogram shows the adjusted <span class="html-italic">p</span>-values, illustrating how the FDR–BH procedure reduces the frequency of small <span class="html-italic">p</span>-values to control the FDR at the desired level (5%).</p> "> Figure 7
<p>Spatial distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The choropleth map classifies trends into five intervals based on TS slope values: “Sharp Decrease” (&lt;−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (&gt;0.5) indicates a notably positive annual rate of change.</p> "> Figure 8
<p>Distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The bar chart classifies trends into five intervals based on the TS slope values: “Sharp Decrease” (&gt;−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (&gt;0.5) indicates a notably positive annual rate of change.</p> "> Figure 9
<p>Spatial distribution of tree cover and tree cover change (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: The first map represents tree cover derived from MOD44B Version 6.1 Vegetation Continuous Fields (VCF) data for 2000, while the second map shows tree cover for 2022. The third map illustrates the percentage changes in tree cover between 2000 and 2022. However, only relevant changes are presented, specifically those that coincide with the statistically significant trends identified by the RTA.</p> "> Figure 10
<p>Moran’s I scatter plot of spatial autocorrelation for significant TS trend slopes filtered with the CMK test and FDR control. Source: Author’s work. Note: The scatter plot presents the standardised Theil–Sen slope values on the x-axis against their spatially lagged counterparts on the y-axis. Each point corresponds to a geographic unit exhibiting a significant trend, filtered through the CMK test and controlled with FDR. The red regression line depicts the overall spatial autocorrelation, with Moran’s I value quantifying the strength and direction of this spatial relationship.</p> ">
Versions Notes

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 − 1 2 ≈ 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 = n n − 1 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:
TS   Slope = Median X j − X i t j − t i
Intercept = X i − MedianSlope   ∗   t i ,   for   i = 1 … 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 i , j such that 1   ≤ i < j   ≤ 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 − 1 2 ≈ 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 = ∑ i = 1 n − 1 ∑ j = i + 1 n sign x i − x j
where:
sign x j − x i = 1 0 − 1 i f   x j − x i > 0 i f   x j − x i = 0 i f   x j − x i < 0
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 sign function is used to evaluate whether the difference between two values in a time series is positive, negative, or zero.
The variance σ 2 of S is used to normalise the statistic and is calculated as follows:
σ 2 = n n − 1 2 n + 5 18 − ∑ t = 1 g t t − 1 2 t + 5 18
where:
σ 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 ( σ 2 = n n − 1 2 n + 5 18 ) accounts for the total variance, while the second term ( ∑ t = 1 g t t − 1 2 t + 5 18 ) corrects for ties in the data.
The Z statistic is used to determine the significance of the trend as follows:
Z   = S − 1 σ i f   S > 0 0 i f   S = 0 S + 1 σ i f   S < 0
and:
p = 2 × Φ − Z
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 Φ 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 = S ¯ m − E S ¯ m σ m
and:
S ¯ m = 1 m ∑ j = 1 m S j
and:
V a r S ¯ m = n n − 1 2 n + 5 18 m = σ 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:
V a r S ¯ m = 1 m 2 ∑ j = 1 m V a r S j + 2 ∑ j = 1 m ∑ l = 1 m − j C o v S j , S j + l
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:
FDR = E V R
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 1 , p 2 , … , p m
where:
p 1 represents the smallest p-value from the set of hypothesis tests; p 2 represents the second-smallest p-value; p m 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 i ≤ i m ⋅ α
where:
p i 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 1 , H 2 , … , H k
where:
H 1 represents the hypothesis corresponding to the smallest p-value; H 2 represents the hypothesis corresponding to the second-smallest p-value; H k 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.

Author Contributions

Conceptualization: O.G.-H.; methods: O.G.-H. and L.V.G.; software: O.G.-H.; formal analysis: O.G.-H. and L.V.G.; investigation: O.G.-H. and L.V.G.; writing—original draft preparation: O.G.-H.; writing—review and editing: O.G.-H. and L.V.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The complete dataset supporting this study’s findings is available from the corresponding author, O.G.-H.

Acknowledgments

The authors would like to express their gratitude to several people. To Lorena Gómez Aparicio for providing photographic material and her research projects in the Parque Natural de Los Alcornocales. In this regard, we thank the colleagues, technicians, and students from IRNAS-CSIC, with whom we worked for several years on field and laboratory tasks in this cork oak-dominated forest in the south of the Iberian Peninsula. Thanks to the reviewers for their contribution to improving the original manuscript. The review process has been challenging but fair, with an excellent outcome. We also thank the editorial team for the waiver of the publication fee.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Kuenzer, C.; Dech, S.; Wagner, W. (Eds.) Remote Sensing Time Series; Springer International Publishing: Cham, Switzerland, 2015; Volume 22, ISBN 978-3-319-15966-9. [Google Scholar]
  2. Hipel, K.; McLeod, A. Nonparametric Tests for Trend Detection. In Time Series Modelling of Water Resources and Environmental Systems; Elsevier: Amsterdam, The Netherlands, 1994; pp. 853–938. [Google Scholar]
  3. Mitchell, A.L.; Rosenqvist, A.; Mora, B. Current Remote Sensing Approaches to Monitoring Forest Degradation in Support of Countries Measurement, Reporting and Verification (MRV) Systems for REDD+. Carbon. Balance Manag. 2017, 12, 9. [Google Scholar] [CrossRef] [PubMed]
  4. Kuenzer, C.; Dech, S.; Wagner, W. Remote Sensing Time Series Revealing Land Surface Dynamics: Status Quo and the Pathway Ahead; Springer: Cham, Switzerland, 2015; pp. 1–24. [Google Scholar]
  5. Eklundh, L.; Jönsson, P. TIMESAT: A Software Package for Time-Series Processing and Assessment of Vegetation Dynamics; Springer: Cham, Switzerland, 2015; pp. 141–158. [Google Scholar]
  6. Hostert, P.; Griffiths, P.; van der Linden, S.; Pflugmacher, D. Time Series Analyses in a New Era of Optical Satellite Data; Springer: Cham, Switzerland, 2015; pp. 25–41. [Google Scholar]
  7. McLeod, A.I.; Hipel, K.W. Tests for Monotonic Trend; Springer: Dordrecht, The Netherlands, 1994; pp. 245–270. [Google Scholar]
  8. Privalsky, V. Time and Frequency Domain Models of Scalar Time Series; Springer: Cham, Switzerland, 2021; pp. 27–44. [Google Scholar]
  9. Theil, H. A Rank-Invariant Method of Linear and Polynomial Regression Analysis I, II and III. In Proceedings of the Section of Sciences, Koninklijke Academie van Wetenschappen te; Springer: Amsterdam, The Netherlands; 1950; pp. 386–392. [Google Scholar]
  10. Sen, P. Estimates of the Regression Coefficient Based on Kendall’s Tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  11. Peng, H.; Wang, S.; Wang, X. Consistency and Asymptotic Distribution of the Theil–Sen Estimator. J. Stat. Plan. Inference 2008, 138, 1836–1850. [Google Scholar] [CrossRef]
  12. Freedman, D. Statistical Models. Theory and Practice; University of California: Berkeley, CA, USA, 2009. [Google Scholar]
  13. Bal, A. Improving the Robustness of the Theil-Sen Estimator Using a Simple Heuristic-Based Modification. Symmetry 2024, 16, 698. [Google Scholar] [CrossRef]
  14. Fernandes, R.; Leblanc, S.G. Parametric (Modified Least Squares) and Non-Parametric (Theil–Sen) Linear Regressions for Predicting Biophysical Parameters in the Presence of Measurement Errors. Remote Sens. Environ. 2005, 95, 303–316. [Google Scholar] [CrossRef]
  15. Mann, H. Nonparametric Tests against Trend. Econometrica 1945, 13, 245–259. [Google Scholar] [CrossRef]
  16. Kendall, M. Rank Correlation Methods; Charles Griffin: London, UK, 1975. [Google Scholar]
  17. Spearman, C. The Proof and Measurement of Association between Two Things. Am. J. Psychol. 1904, 15, 72. [Google Scholar] [CrossRef]
  18. Helsel, D.R.; Frans, L.M. Regional Kendall Test for Trend. Environ. Sci. Technol. 2006, 40, 4066–4073. [Google Scholar] [CrossRef]
  19. Hoaglin, D.; Mosteller, F.; Tukey, J. Understanding Robust and Exploratory Data Analysis; John Wiley and Sons: New York, NY, USA, 2000. [Google Scholar]
  20. Eastman, J. ClarkLabs. TerrSet: Geospatial Monitoring and Modeling. Tutorial; Clark University: Worcester, MA, USA, 2016. [Google Scholar]
  21. Cliff, A.; Ord, J. Spatial Autocorrelation; Pion: London, UK, 1973. [Google Scholar]
  22. Douglas, E.M.; Vogel, R.M.; Kroll, C.N. Trends in Floods and Low Flows in the United States: Impact of Spatial Correlation. J. Hydrol. 2000, 240, 90–105. [Google Scholar] [CrossRef]
  23. Renard, B.; Lang, M.; Bois, P.; Dupeyrat, A.; Mestre, O.; Niel, H.; Sauquet, E.; Prudhomme, C.; Parey, S.; Paquet, E.; et al. Regional Methods for Trend Detection: Assessing Field Significance and Regional Consistency. Water Resour. Res. 2008, 44, W08419. [Google Scholar] [CrossRef]
  24. Neeti, N.; Eastman, J.R. A Contextual Mann-Kendall Approach for the Assessment of Trend Significance in Image Time Series. Trans. GIS 2011, 15, 599–611. [Google Scholar] [CrossRef]
  25. Piedmont, R.L. Type I Errors. In Encyclopedia of Quality of Life and Well-Being Research; Springer: Dordrecht, The Netherlands, 2014; p. 6758. [Google Scholar]
  26. Wang, F.; Shao, W.; Yu, H.; Kan, G.; He, X.; Zhang, D.; Ren, M.; Wang, G. Re-Evaluation of the Power of the Mann-Kendall Test for Detecting Monotonic Trends in Hydrometeorological Time Series. Front. Earth Sci. 2020, 8, 14. [Google Scholar] [CrossRef]
  27. Clements, N.; Sarkar, S.K.; Zhao, Z.; Kim, D.Y. Applying Multiple Testing Procedures to Detect Change in East African Vegetation. Ann. Appl. Stat. 2014, 8, 286–308. [Google Scholar] [CrossRef]
  28. James, G.; Witten, D.; Hastie, T.; Tibshirani, R.; Taylor, J. Multiple Testing. In An Introduction to Statistical Learning; Springer Texts in Statistics; Springer: Cham, Switzerland, 2023. [Google Scholar] [CrossRef]
  29. Heumann, B.W. The Multiple Comparison Problem in Empirical Remote Sensing. Photogramm. Eng. Remote Sens. 2015, 81, 921–926. [Google Scholar] [CrossRef]
  30. Cortés, J.; Mahecha, M.; Reichstein, M.; Brenning, A. Accounting for Multiple Testing in the Analysis of Spatio-Temporal Environmental Data. Environ. Ecol. Stat. 2020, 27, 293–318. [Google Scholar] [CrossRef]
  31. Gil, L.; Varela, M. EUFORGEN Technical Guidelines for Genetic Conservation and Use for Cork Oak (Quercus Suber); American Society for Photogrammetry and Remote Sensing: Rome, Italy, 2008. [Google Scholar]
  32. Agrillo, E.; Alessi, N.; Jiménez-Alfaro, B.; Casella, L.; Angelini, P.; Argagnon, O.; Crespo, G.; Fernández-González, F.; Monteiro-Henriques, T.; Neto, C.S.; et al. The Use of Large Databases to Characterize Habitat Types: The Case of Quercus Suber Woodlands in Europe. Rend. Lincei Sci. Fis. Nat. 2018, 29, 283–293. [Google Scholar] [CrossRef]
  33. Aronson, J.; Pereira, J.; Pausas, J. Cork Oak Woodlands on the Edge Ecology, Adaptive Management, and Restoration; Island Press: Washington, DC, USA, 2009. [Google Scholar]
  34. Pereira, H. Cork. Biology Production and Uses; Elsevier: Amsterdam, The Netherlands, 2007; ISBN 9780444529671. [Google Scholar]
  35. Aroso, I.M.; Araújo, A.R.; Pires, R.A.; Reis, R.L. Cork: Current Technological Developments and Future Perspectives for This Natural, Renewable, and Sustainable Material. ACS Sustain. Chem. Eng. 2017, 5, 11130–11146. [Google Scholar] [CrossRef]
  36. San-Miguel-Ayanz, J.; de Rigo, D.; Caudullo, G.; Houston Durrant, T.; Mauri, A. European Atlas of Forest Tree Species; Publication Office of the European Union: Luxembourg, 2016. [Google Scholar]
  37. Houston Durrant, T.; de Rigo, D.; Caudullo, G. Quercus Suber in Europe: Distribution, Habitat, Usage and Threats. In European Atlas of Forest Tree Species; Publication Office of the European Union: Luxembourg, 2016; pp. 164–165. [Google Scholar]
  38. Olson, D.M.; Dinerstein, E.; Wikramanayake, E.D.; Burgess, N.D.; Powell, G.V.N.; Underwood, E.C.; D’amico, J.A.; Itoua, I.; Strand, H.E.; Morrison, J.C.; et al. Terrestrial Ecoregions of the World: A New Map of Life on EarthA New Global Map of Terrestrial Ecoregions Provides an Innovative Tool for Conserving Biodiversity. Bioscience 2001, 51, 933–938. [Google Scholar] [CrossRef]
  39. Serrasolses, I.; Pérez-devesa, M.; Vilagrosa, A.; Pausas, J.G.; Sauras, T.; Cortina, J.; Vallejo, V.R. Cork Oak Distribution. In Cork Oak Woodlands on the Edge: Conservation, Adaptive Management, and Restoration; Island Press: Washington, DC, USA, 2009; pp. 89–99. [Google Scholar]
  40. Sørensen, I.H.; Torralba, M.; Quintas-Soriano, C.; Muñoz-Rojas, J.; Plieninger, T. Linking Cork to Cork Oak Landscapes: Mapping the Value Chain of Cork Production in Portugal. Front. Sustain. Food Syst. 2021, 5, 787045. [Google Scholar] [CrossRef]
  41. Faísca, C.M. Economy and Ecology in the Iberian Cork Oak Forests: Land Use in the Second Half of the 19th Century. Configurações 2020, 26, 83–105. [Google Scholar] [CrossRef]
  42. Parsons, J.J. The Cork Oak Forests and the Evolution of the Cork Industry in Southern Spain and Portugal. Econ. Geogr. 1962, 38, 195. [Google Scholar] [CrossRef]
  43. Arosa, M.L.; Bastos, R.; Cabral, J.A.; Freitas, H.; Costa, S.R.; Santos, M. Long-Term Sustainability of Cork Oak Agro-Forests in the Iberian Peninsula: A Model-Based Approach Aimed at Supporting the Best Management Options for the Montado Conservation. Ecol. Modell. 2017, 343, 68–79. [Google Scholar] [CrossRef]
  44. Li, F.; Qian, H.; Sardans, J.; Amishev, D.Y.; Wang, Z.; Zhang, C.; Wu, T.; Xu, X.; Tao, X.; Huang, X. Evolutionary History Shapes Variation of Wood Density of Tree Species across the World. Plant Divers. 2024, 46, 283–293. [Google Scholar] [CrossRef]
  45. Kim, H.N.; Jin, H.Y.; Kwak, M.J.; Khaine, I.; You, H.N.; Lee, T.Y.; Ahn, T.H.; Woo, S.Y. Why Does Quercus Suber Species Decline in Mediterranean Areas? J. Asia Pac. Biodivers. 2017, 10, 337–341. [Google Scholar] [CrossRef]
  46. Vagniluca, S.; Moricca, S.; Vigniluca Source, S. European Expansion of Oak Decline: Involved Microrganisms and Methodological Approaches. Phytopathol. Mediterr. 1995, 34, 207–226. [Google Scholar]
  47. Acha, A.; Newing, H.S. Cork Oak Landscapes, Promised or Compromised Lands? A Case Study of a Traditional Cultural Landscape in Southern Spain. Hum. Ecol. 2015, 43, 601–611. [Google Scholar] [CrossRef]
  48. Costa, A.; Pereira, H.; Madeira, M. Analysis of Spatial Patterns of Oak Decline in Cork Oak Woodlands in Mediterranean Conditions. Ann. For. Sci. 2010, 67, 204. [Google Scholar] [CrossRef]
  49. Camarero, J.J.; Gazol, A.; Valeriano, C.; Colangelo, M.; Rubio-Cuadrado, Á. Growth Responses to Climate and Drought in Relict Cork Oak Populations as a Benchmark of the Species Tolerance. Forests 2024, 15, 72. [Google Scholar] [CrossRef]
  50. Matías, L.; Pérez-Ramos, I.M.; Gómez-Aparicio, L. Are Northern-Edge Populations of Cork Oak More Sensitive to Drought than Those of the Southern Edge? Environ. Exp. Bot. 2019, 163, 78–85. [Google Scholar] [CrossRef]
  51. Gentilesca, T.; Camarero, J.; Colangelo, M.; Nolè, A.; Ripullone, F. Drought-Induced Oak Decline in the Western Mediterranean Region: An Overview on Current Evidences, Mechanisms and Management Options to Improve Forest Resilience. iForest 2017, 10, 796–806. [Google Scholar] [CrossRef]
  52. Tiberi, R.; Branco, M.; Bracalini, M.; Croci, F.; Panzavolta, T. Cork Oak Pests: A Review of Insect Damage and Management. Ann. For. Sci. 2016, 73, 219–232. [Google Scholar] [CrossRef]
  53. Moricca, S.; Linaldeddu, B.T.; Ginetti, B.; Scanu, B.; Franceschini, A.; Ragazzi, A. Endemic and Emerging Pathogens Threatening Cork Oak Trees: Management Options for Conserving a Unique Forest Ecosystem. Plant Dis. 2016, 100, 2184–2193. [Google Scholar] [CrossRef]
  54. Brasier, C. Phytophthora Cinnamomi and Oak Decline in Southern Europe. Environmental Constraints Including Climate Change. Ann. Des Sci. For. 1996, 53, 347–358. [Google Scholar] [CrossRef]
  55. Moralejo, E.; García-Muñoz, J.A.; Descals, E. Susceptibility of Iberian Trees to Phytophthora Ramorum and P. Cinnamomi. Plant Pathol. 2009, 58, 271–283. [Google Scholar] [CrossRef]
  56. de Sampaio e Paiva Camilo-Alves, C.; da Clara, M.I.E.; de Almeida Ribeiro, N.M.C. Decline of Mediterranean Oak Trees and Its Association with Phytophthora Cinnamomi: A Review. Eur. J. For. Res. 2013, 132, 411–432. [Google Scholar] [CrossRef]
  57. Brasier, C. Oak Tree Mortality in Iberia. Nature 1992, 360, 539. [Google Scholar] [CrossRef]
  58. Gómez-Aparicio, L.; Ibáñez, B.; Serrano, M.S.; De Vita, P.; Ávila, J.M.; Pérez-Ramos, I.M.; García, L.V.; Esperanza Sánchez, M.; Marañón, T. Spatial Patterns of Soil Pathogens in Declining Mediterranean Forests: Implications for Tree Species Regeneration. New Phytol. 2012, 194, 1014–1024. [Google Scholar] [CrossRef]
  59. Domínguez-Begines, J.; De Deyn, G.B.; García, L.V.; Eisenhauer, N.; Gómez-Aparicio, L. Cascading Spatial and Trophic Impacts of Oak Decline on the Soil Food Web. J. Ecol. 2019, 107, 1199–1214. [Google Scholar] [CrossRef]
  60. Eastman, J. ClarkLabs. Earth Trends Modeler in TerrSet 2020; ClarkLabs; Clark University: Worcester, MA, USA, 2021. [Google Scholar]
  61. Pohlert, T. Trend: Non-Parametric Trend Tests and Change-Point Detection; R Package Version 1.1.6; 2023; Available online: https://cran.r-project.org/web/packages/trend/vignettes/trend.pdf (accessed on 7 October 2024).
  62. McLeod, A. Kendall: Kendall Rank Correlation and Mann-Kendall Trend Test. R Package Version 2.2.1; 2022; Available online: https://cran.r-project.org/web/packages/Kendall/Kendall.pdf (accessed on 7 October 2024).
  63. Patakamuri, S.; O’Brien, N. Modifiedmk: Modified Versions of Mann Kendall and Spearman’s Rho Trend Tests. R Package Version 1.6; 2021; Available online: https://cran.r-project.org/web/packages/modifiedmk/modifiedmk.pdf (accessed on 7 October 2024).
  64. Jacques, B.; Aronson, J. Biology and Wildlife of the Mediterranean; Oxford University Press: Oxford, UK, 1999. [Google Scholar]
  65. Rodríguez-Sánchez, F.; Pérez-Barrales, R.; Ojeda, F.; Vargas, P.; Arroyo, J. The Strait of Gibraltar as a Melting Pot for Plant Biodiversity. Quat. Sci. Rev. 2008, 27, 2100–2117. [Google Scholar] [CrossRef]
  66. Gómez-Zotano, J.; Alcántara-Manzanares, J.; Olmedo-Cobo, J.A.; Martínez-Ibarra, E. La Sistematización Del Clima Mediterráneo: Identificación, Clasificación y Caracterización Climática de Andalucía (España). Rev. Geogr. Norte Gd. 2015, 61, 161–180. [Google Scholar] [CrossRef]
  67. Ojeda, F.; Marañón, T.; Arroyo, J. Plant Diversity Patterns in the Aljibe Mountains (S. Spain): A Comprehensive Account. Biodivers. Conserv. 2000, 9, 1323–1343. [Google Scholar] [CrossRef]
  68. Gutiérrez-Hernández, O.; Senciales-González, J.M.; García, L.V. Evolución de La Superficie Forestal En Andalucía. Procesos y Factores. Rev. Estud. Andal. 2016, 33, 111–148. [Google Scholar] [CrossRef]
  69. Regional Ministry for the Environment of Andalusia. Procesos de Decaimiento Forestal (La Seca) En Andalucía. Situación Del Conocimiento; Consejería de Medio Ambiente; Junta de Andalucía: Córdoba, Spain, 2009; Available online: https://www.juntadeandalucia.es/medioambiente/portal/documents/20151/38640721/procesos_decaimiento_forestal_1.pdf/f1b89e3d-901c-9e02-e6ee-72161c595121?t=1645789298123 (accessed on 7 October 2024).
  70. Avila, J.M.; Linares, J.C.; García-Nogales, A.; Sánchez, M.E.; Gómez-Aparicio, L. Across-Scale Patterning of Plant–Soil–Pathogen Interactions in Quercus Suber Decline. Eur. J. For. Res. 2017, 136, 677–688. [Google Scholar] [CrossRef]
  71. Homet, P.; González, M.; Matías, L.; Godoy, O.; Pérez-Ramos, I.M.; García, L.V.; Gómez-Aparicio, L. Exploring Interactive Effects of Climate Change and Exotic Pathogens on Quercus Suber Performance: Damage Caused by Phytophthora Cinnamomi Varies across Contrasting Scenarios of Soil Moisture. Agric. For. Meteorol. 2019, 276–277, 107605. [Google Scholar] [CrossRef]
  72. Ibáñez, B.; Gómez-Aparicio, L.; Stoll, P.; Ávila, J.M.; Pérez-Ramos, I.M.; Marañón, T. A Neighborhood Analysis of the Consequences of Quercus Suber Decline for Regeneration Dynamics in Mediterranean Forests. PLoS ONE 2015, 10, e0117827. [Google Scholar] [CrossRef]
  73. Gómez-Aparicio, L.; Domínguez-Begines, J.; Kardol, P.; Ávila, J.M.; Ibáñez, B.; García, L.V. Plant-Soil Feedbacks in Declining Forests: Implications for Species Coexistence. Ecology 2017, 98, 1908–1921. [Google Scholar] [CrossRef]
  74. DiMiceli, C.; Sohlberg, R.; Townshend, J. MODIS/Terra Vegetation Continuous Fields Yearly L3 Global 250m SIN Grid V061 [Data Set]; NASA EOSDIS Land Processes Distributed Active Archive Center: Sioux Falls, SD, USA, 2022.
  75. Busetto, L.; Ranghetti, L. MODIStsp: An R Package for Automatic Preprocessing of MODIS Land Products Time Series. Comput. Geosci. 2016, 97, 40–48. [Google Scholar] [CrossRef]
  76. Neeti, N.; Ronald Eastman, J. Novel Approaches in Extended Principal Component Analysis to Compare Spatio-Temporal Patterns among Multiple Image Time Series. Remote Sens. Environ. 2014, 148, 84–96. [Google Scholar] [CrossRef]
  77. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B (Methodol.) 1995, 57, 89–300. [Google Scholar] [CrossRef]
  78. Moran, P.A.P. Notes on Continuous Stochastic Phenomena. Biometrika 1950, 37, 17. [Google Scholar] [CrossRef]
  79. Goodchild, M. Spatial Autocorrelation; Geo Books: Norwich, UK, 1986. [Google Scholar]
  80. R Core Team. R: A Language and Environment for Statistical Computing; R Core Team: Vienna, Austria, 2023. [Google Scholar]
  81. Eastman, J.; ClarkLabs. TerrSet: Geospatial Monitoring and Modeling Software; Version 19.08; 2023; Available online: https://clarklabs.org/terrset/ (accessed on 7 October 2024)Version 19.08.
  82. Hijmans, R. Terra: Spatial Data Analysis 2024. Available online: https://cran.r-project.org/web/packages/terra/terra.pdf (accessed on 7 October 2024).
  83. Bivand, R. R Packages for Analyzing Spatial Data: A Comparative Case Study with Areal Data. Geogr. Anal. 2022, 54, 488–518. [Google Scholar] [CrossRef]
  84. Ventura, V.; Paciorek, C.J.; Risbey, J.S. Controlling the Proportion of Falsely Rejected Hypotheses When Conducting Multiple Tests with Climatological Data. J. Clim. 2004, 17, 4343–4356. [Google Scholar] [CrossRef]
  85. García, L.V. Controlling the False Discovery Rate in Ecological Research. Trends Ecol. Evol. 2003, 18, 553–554. [Google Scholar] [CrossRef]
  86. Ibáñez, B.; Gómez-Aparicio, L.; Ávila, J.M.; Pérez-Ramos, I.M.; Marañón, T. Effects of Quercus Suber Decline on Woody Plant Regeneration: Potential Implications for Successional Dynamics in Mediterranean Forests. Ecosystems 2017, 20, 630–644. [Google Scholar] [CrossRef]
  87. Ibáñez, B.; Ibáñez, I.; Gómez-Aparicio, L.; Ruiz-Benito, P.; García, L.V.; Marañón, T. Contrasting Effects of Climate Change along Life Stages of a Dominant Tree Species: The Importance of Soil-Climate Interactions. Divers. Distrib. 2014, 20, 872–883. [Google Scholar] [CrossRef]
  88. Gutiérrez-Hernández, O.; Sánchez, M.E.; Ramo, C.; Sánchez-Solana, J.E.; García, L.V. The Occurrence of Phytophthora Cinnamomi in Southern Spain: Presence—Absence Records and Potential Distribution Area. IOBC-WPRS Bull. 2017, 127, 105–109. [Google Scholar]
  89. Rebollo, P.; Moreno-Fernández, D.; Cruz-Alonso, V.; Gazol, A.; Rodríguez-Rey, M.; Astigarraga, J.; Zavala, M.A.; Gómez-Aparicio, L.; Andivia, E.; Miguel-Romero, S.; et al. Recent Increase in Tree Damage and Mortality and Their Spatial Dependence on Drought Intensity in Mediterranean Forests. Landsc. Ecol. 2024, 39, 38. [Google Scholar] [CrossRef]
  90. Owen, H.J.F.; Lines, E.R. Common Field Measures and Geometric Assumptions of Tree Shape Produce Consistently Biased Estimates of Tree and Canopy Structure in Mixed Mediterranean Forests. Ecol. Indic. 2024, 165, 112219. [Google Scholar] [CrossRef]
  91. Mu, X.; Hu, M.; Song, W.; Ruan, G.; Ge, Y.; Wang, J.; Huang, S.; Yan, G. Evaluation of Sampling Methods for Validation of Remotely Sensed Fractional Vegetation Cover. Remote Sens. 2015, 7, 16164–16182. [Google Scholar] [CrossRef]
  92. Gómez-Aparicio, L.; Domínguez-Begines, J.; Villa-Sanabria, E.; García, L.V.; Muñoz-Pajares, A.J. Tree Decline and Mortality Following Pathogen Invasion Alters the Diversity, Composition and Network Structure of the Soil Microbiome. Soil Biol. Biochem. 2022, 166, 108560. [Google Scholar] [CrossRef]
  93. Serrano, M.S.; Villa-Sanabria, E.; Homet, P.; Gutiérrez, E.; Gómez-Aparicio, L. Impact of a Drier Climate on the Exotic Pathogen Phytophthora Cinnamomi in Mediterranean Forests Differing in Soil Properties and Species Composition. For. Ecol. Manag. 2024, 556, 121721. [Google Scholar] [CrossRef]
  94. Güçlü, Y.S. Multiple Şen-Innovative Trend Analyses and Partial Mann-Kendall Test. J. Hydrol. 2018, 566, 685–704. [Google Scholar] [CrossRef]
  95. Bayazit, M.; Önöz, B. To Prewhiten or Not to Prewhiten in Trend Analysis? Hydrol. Sci. J. 2007, 52, 611–624. [Google Scholar] [CrossRef]
  96. Yue, S.; Wang, C.Y. Applicability of Prewhitening to Eliminate the Influence of Serial Correlation on the Mann-Kendall Test. Water Resour. Res. 2002, 38, 4-1–4-6. [Google Scholar] [CrossRef]
Remotesensing 16 03886 g001
Figure 1. Study area: Los Alcornocales Natural Park (≈174,000 ha), a representative cork oak-dominated mixed Mediterranean forest in the Southwest Iberian Peninsula (Andalusia, Spain). Source: Author’s work. Data from the Department of the Environment of the Regional Government of Andalusia, specifically from the Andalusian Network of Environmental Information (REDIAM).
Figure 1. Study area: Los Alcornocales Natural Park (≈174,000 ha), a representative cork oak-dominated mixed Mediterranean forest in the Southwest Iberian Peninsula (Andalusia, Spain). Source: Author’s work. Data from the Department of the Environment of the Regional Government of Andalusia, specifically from the Andalusian Network of Environmental Information (REDIAM).
Remotesensing 16 03886 g001
Remotesensing 16 03886 g002
Figure 2. Evidence of cork oak forest decline in Los Alcornocales Natural Park. Source: Lorena Gómez-Aparicio.
Figure 2. Evidence of cork oak forest decline in Los Alcornocales Natural Park. Source: Lorena Gómez-Aparicio.
Remotesensing 16 03886 g002
Remotesensing 16 03886 g003
Figure 3. Research workflow diagram for robust trend analysis (RTA). Source: Author’s work.
Figure 3. Research workflow diagram for robust trend analysis (RTA). Source: Author’s work.
Remotesensing 16 03886 g003
Remotesensing 16 03886 g004
Figure 4. Density of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: TS slope (dashed green line); combined TS slope and CMK significance (thin solid green line); combined TS slope and CMK significance with FDR control (thick solid green line). Axes details: X-axis—TS slope; Y-axis—spatial distribution density.
Figure 4. Density of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: TS slope (dashed green line); combined TS slope and CMK significance (thin solid green line); combined TS slope and CMK significance with FDR control (thick solid green line). Axes details: X-axis—TS slope; Y-axis—spatial distribution density.
Remotesensing 16 03886 g004
Remotesensing 16 03886 g005
Figure 5. Controlling the FDR for interannual tree cover trends (2000–2022) in Los Alcornocales Natural Park: FDR–BH p-value adjustment for the CMK test. Source: Author’s work. Note: The x-axis represents the rank (k) of 27,828 p-values sorted in ascending order, while the y-axis displays the observed p-values. The celestial blue line symbolises the threshold defined by the FDR, which is calculated using the following expression: k · α m . Herein, k denotes the rank, α is the significance level, and m is the total number of tests conducted. The part of the line coloured blue indicates the p-values that are deemed significant following the application of the FDR control, using the FDR–BH procedure. Conversely, the part of the line coloured red represents the p-values that are not considered significant according to this adjustment. Altogether, this bi-coloured line comprises 27,828 points (pixels), corresponding to the number of tests performed (m).
Figure 5. Controlling the FDR for interannual tree cover trends (2000–2022) in Los Alcornocales Natural Park: FDR–BH p-value adjustment for the CMK test. Source: Author’s work. Note: The x-axis represents the rank (k) of 27,828 p-values sorted in ascending order, while the y-axis displays the observed p-values. The celestial blue line symbolises the threshold defined by the FDR, which is calculated using the following expression: k · α m . Herein, k denotes the rank, α is the significance level, and m is the total number of tests conducted. The part of the line coloured blue indicates the p-values that are deemed significant following the application of the FDR control, using the FDR–BH procedure. Conversely, the part of the line coloured red represents the p-values that are not considered significant according to this adjustment. Altogether, this bi-coloured line comprises 27,828 points (pixels), corresponding to the number of tests performed (m).
Remotesensing 16 03886 g005
Remotesensing 16 03886 g006
Figure 6. Histogram of original and adjusted p-values (FDR). Both histograms share the same scale axes for easy comparison. Source: Author’s work. Note: The x-axis shows p-values ranging from 0 to 1, and the y-axis shows frequencies up to 12,000 (pixels). The first bin in each histogram represents pixels with significant trends (α = 0.05). The left histogram displays the original p-values, heavily skewed towards 0, with a green dashed line at 8218 (pixels), which serves as a reference for comparison with the adjusted p-values at the alpha level of 0.05. The right histogram shows the adjusted p-values, illustrating how the FDR–BH procedure reduces the frequency of small p-values to control the FDR at the desired level (5%).
Figure 6. Histogram of original and adjusted p-values (FDR). Both histograms share the same scale axes for easy comparison. Source: Author’s work. Note: The x-axis shows p-values ranging from 0 to 1, and the y-axis shows frequencies up to 12,000 (pixels). The first bin in each histogram represents pixels with significant trends (α = 0.05). The left histogram displays the original p-values, heavily skewed towards 0, with a green dashed line at 8218 (pixels), which serves as a reference for comparison with the adjusted p-values at the alpha level of 0.05. The right histogram shows the adjusted p-values, illustrating how the FDR–BH procedure reduces the frequency of small p-values to control the FDR at the desired level (5%).
Remotesensing 16 03886 g006
Remotesensing 16 03886 g007
Figure 7. Spatial distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The choropleth map classifies trends into five intervals based on TS slope values: “Sharp Decrease” (<−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (>0.5) indicates a notably positive annual rate of change.
Figure 7. Spatial distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The choropleth map classifies trends into five intervals based on TS slope values: “Sharp Decrease” (<−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (>0.5) indicates a notably positive annual rate of change.
Remotesensing 16 03886 g007
Remotesensing 16 03886 g008
Figure 8. Distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The bar chart classifies trends into five intervals based on the TS slope values: “Sharp Decrease” (>−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (>0.5) indicates a notably positive annual rate of change.
Figure 8. Distribution of interannual trends in tree cover (2000–2022) in Los Alcornocales Natural Park using RTA. Source: Author’s work. Note: The bar chart classifies trends into five intervals based on the TS slope values: “Sharp Decrease” (>−0.5) indicates a notably negative annual rate of change; “Moderate Decrease” (−0.5 to −0.1) indicates a moderate annual rate of change; “Stable” (−0.1 to 0.1) implies a negligible annual rate of change; “Moderate Increase” (0.1 to 0.5) indicates a notably positive annual rate of change; and “Sharp Increase” (>0.5) indicates a notably positive annual rate of change.
Remotesensing 16 03886 g008
Remotesensing 16 03886 g009
Figure 9. Spatial distribution of tree cover and tree cover change (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: The first map represents tree cover derived from MOD44B Version 6.1 Vegetation Continuous Fields (VCF) data for 2000, while the second map shows tree cover for 2022. The third map illustrates the percentage changes in tree cover between 2000 and 2022. However, only relevant changes are presented, specifically those that coincide with the statistically significant trends identified by the RTA.
Figure 9. Spatial distribution of tree cover and tree cover change (2000–2022) in Los Alcornocales Natural Park. Source: Author’s work. Note: The first map represents tree cover derived from MOD44B Version 6.1 Vegetation Continuous Fields (VCF) data for 2000, while the second map shows tree cover for 2022. The third map illustrates the percentage changes in tree cover between 2000 and 2022. However, only relevant changes are presented, specifically those that coincide with the statistically significant trends identified by the RTA.
Remotesensing 16 03886 g009
Remotesensing 16 03886 g010
Figure 10. Moran’s I scatter plot of spatial autocorrelation for significant TS trend slopes filtered with the CMK test and FDR control. Source: Author’s work. Note: The scatter plot presents the standardised Theil–Sen slope values on the x-axis against their spatially lagged counterparts on the y-axis. Each point corresponds to a geographic unit exhibiting a significant trend, filtered through the CMK test and controlled with FDR. The red regression line depicts the overall spatial autocorrelation, with Moran’s I value quantifying the strength and direction of this spatial relationship.
Figure 10. Moran’s I scatter plot of spatial autocorrelation for significant TS trend slopes filtered with the CMK test and FDR control. Source: Author’s work. Note: The scatter plot presents the standardised Theil–Sen slope values on the x-axis against their spatially lagged counterparts on the y-axis. Each point corresponds to a geographic unit exhibiting a significant trend, filtered through the CMK test and controlled with FDR. The red regression line depicts the overall spatial autocorrelation, with Moran’s I value quantifying the strength and direction of this spatial relationship.
Remotesensing 16 03886 g010
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Gutiérrez-Hernández, O.; García, L.V. Robust Trend Analysis in Environmental Remote Sensing: A Case Study of Cork Oak Forest Decline. Remote Sens. 2024, 16, 3886. https://doi.org/10.3390/rs16203886

AMA Style

Gutiérrez-Hernández O, García LV. Robust Trend Analysis in Environmental Remote Sensing: A Case Study of Cork Oak Forest Decline. Remote Sensing. 2024; 16(20):3886. https://doi.org/10.3390/rs16203886

Chicago/Turabian Style

Gutiérrez-Hernández, Oliver, and Luis V. García. 2024. "Robust Trend Analysis in Environmental Remote Sensing: A Case Study of Cork Oak Forest Decline" Remote Sensing 16, no. 20: 3886. https://doi.org/10.3390/rs16203886

APA Style

Gutiérrez-Hernández, O., & García, L. V. (2024). Robust Trend Analysis in Environmental Remote Sensing: A Case Study of Cork Oak Forest Decline. Remote Sensing, 16(20), 3886. https://doi.org/10.3390/rs16203886

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop