[go: up one dir, main page]
Next Article in Journal
A Novel Approach for Farmland Size Estimation in Small-Scale Agriculture Using Edge Counting and Remote Sensing
Next Article in Special Issue
Mapping Shrub Biomass at 10 m Resolution by Integrating Field Measurements, Unmanned Aerial Vehicles, and Multi-Source Satellite Observations
Previous Article in Journal
Evaluation of Sentinel-5P TROPOMI Methane Observations at Northern High Latitudes
Previous Article in Special Issue
Classifying Stand Compositions in Clover Grass Based on High-Resolution Multispectral UAV Images
 
 
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 16
10.3390/rs16162980
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

New Insights on the Information Content of the Normalized Difference Vegetation Index Sentinel-2 Time Series for Assessing Vegetation Dynamics

by
César Sáenz
1,2,3,
Víctor Cicuéndez
1,4,
Gabriel García
5,
Diego Madruga
1,
Laura Recuero
3,6,
Alfonso Bermejo-Saiz
1,
Javier Litago
6,
Ignacio de la Calle
2 and
Alicia Palacios-Orueta
1,3,*
1
Departamento de Ingeniería Agroforestal, ETSIAAB, Universidad Politécnica de Madrid, Av. Puerta de Hierro, nº 2—4, Ciudad Universitaria, 28040 Madrid, Spain
2
Quasar Science Resources, S.L., Camino de las Ceudas 2, Las Rozas de Madrid, 28232 Madrid, Spain
3
Centro de Estudios e Investigación para la Gestión de Riesgos Agrarios y Medioambientales (CEIGRAM), Universidad Politécnica de Madrid, C/Senda del Rey 13, 28040 Madrid, Spain
4
Departamento de Física de la Tierra y Astrofísica, Facultad de Ciencias Físicas, Universidad Complutense de Madrid (UCM), 28040 Madrid, Spain
5
Department of Computer Architecture and Automation, Complutense University of Madrid, 28040 Madrid, Spain
6
Departamento de Economía Agraria, Estadística y Gestión de Empresas, ETSIAAB, Universidad Politécnica de Madrid (UPM), Av. Puerta de Hierro, nº 2—4, Ciudad Universitaria, 28040 Madrid, Spain
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(16), 2980; https://doi.org/10.3390/rs16162980
Submission received: 18 June 2024 / Revised: 10 August 2024 / Accepted: 11 August 2024 / Published: 14 August 2024
(This article belongs to the Special Issue Crops and Vegetation Monitoring with Remote/Proximal Sensing II)
Browse Figures
Graphical abstract
"> Figure 1
<p>Image in true color from © Google Earth (Landsat/Copernicus). (<b>a</b>) Map of Spain and (<b>b</b>) study area.</p> "> Figure 2
<p>Sentinel-2 image processing flowchart of time series filtering.</p> "> Figure 3
<p>Spatial distribution of the vegetation species in the study area.</p> "> Figure 4
<p>Valid observations according to the SCL band of Sentinel-2. (<b>a</b>) Spatial distribution in percentage and (<b>b</b>) number of pixels for each level of valid observations.</p> "> Figure 5
<p>(<b>a</b>) Length of the longest gap within time series; (<b>b</b>) season in which the longest gap occurs.</p> "> Figure 6
<p>Interpolating Efficiency Indicator for Tile 30TUN: (<b>a</b>) spatial distribution and (<b>b</b>) frequency distribution.</p> "> Figure 7
<p>Time series of pure pixels for a wheat (<b>a</b>), barley (<b>c</b>), maize (<b>e</b>), and alfalfa (<b>g</b>) crop. Time series of the interpolated NDVI (<b><span style="color:#00B0F0">−</span></b>) and its four filters: SG (<b><span style="color:red">−</span></b>), FFT (<b><span style="color:#92D050">−</span></b>), WHI (<b>−</b>), and MVF (<b><span style="color:yellow">−</span></b>), and (<b>b</b>,<b>d</b>,<b>f</b>,<b>h</b>) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.</p> "> Figure 8
<p>Time series of pure pixels for a beech forest (<b>a</b>) and a pine forest (<b>c</b>). Time series of the interpolated NDVI (<b><span style="color:#00B0F0">−</span></b>) and its four filters: SG (<b><span style="color:red">−</span></b>), FFT (<b><span style="color:#92D050">−</span></b>), WHI (<b>−</b>), and MVF (<b><span style="color:yellow">−</span></b>), and (<b>b</b>,<b>d</b>) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.</p> "> Figure 9
<p>Representation of the evolution at the image level of the interpolation and filtering processes of the Sentinel-2 NDVI time series using the Whittaker filter. (<b>a</b>) NDVI image (10/06/20) without interpolation, (<b>b</b>) NDVI image (15/06/20) without interpolation, (<b>c</b>) NDVI image (20/06/20) without interpolation, (<b>d</b>) NDVI image (15/06/20) interpolated, and (<b>e</b>) NDVI image (15/06/20) filtered using the Whittaker filter.</p> "> Figure 10
<p>Values of the Q-test for the short term (lags 1, 2, 3, 4, 5, 6, and 7).</p> "> Figure 11
<p>Values of the Q-test at lag 36 (6 months), lag 73 (one year), and lag 146 (two years).</p> "> Figure 12
<p>Average Fk value for period 73 (one year) for each vegetation type (left axis) and percentage of increase after interpolation and filtering (right axis).</p> "> Figure 13
<p>Time series of a pixel declared as alfalfa implemented during the first year.</p> ">
Versions Notes

Abstract

:
The Sentinel-2 NDVI time series information content from 2017 to 2023 at a 10 m spatial resolution was evaluated based on the NDVI temporal dependency in five scenarios in central Spain. First, time series were interpolated and then filtered using the Savitzky–Golay, Fast Fourier Transform, Whittaker, and Maximum Value filters. Temporal dependency was assessed using the Q-Ljung-Box and Fisher’s Kappa tests, and similarity between raw and filtered time series was assessed using Correlation Coefficient and Root Mean Square Error. An Interpolating Efficiency Indicator (IEI) was proposed to summarize the number and temporal distribution of low-quality observations. Type of climate, atmospheric disturbances, land cover dynamics, and management were the main sources of variability in five scenarios: (1) rainfed wheat and barley presented high short-term variability due to clouds (lower IEI in winter and spring) during the growing cycle and high interannual variability due to precipitation; (2) maize showed stable summer cycles (high IEI) and low interannual variability due to irrigation; (3) irrigated alfalfa was cut five to six times during summer, resulting in specific intra-annual variability; (4) beech forest showed a strong and stable summer cycle, despite the short-term variability due to clouds (low IEI); and (5) evergreen pine forest had a highly variable growing cycle due to fast responses to temperature and precipitation through the year and medium IEI values. Interpolation after removing non-valid observations resulted in an increase in temporal dependency (Q-test), particularly a short term in areas with low IEI values. The information improvement made it possible to identify hidden periodicities and trends using the Fisher’s Kappa test. The SG filter showed high similarity values and weak influence on dynamics, while the MVF showed an overestimation of the NDVI values.

Graphical Abstract">

Graphical Abstract

1. Introduction

Satellite time series constitutes one of the most extensive sources of information about the environment. The large spatial coverage acquired at regular time intervals results in enormous data sets that go beyond the concept of big data. This type of information provides spatial and temporal coherence to uncover patterns of environmental processes yielding, sometimes, more robust information than field data, especially in heterogeneous areas.
In this context, Sentinel-2 provides the most informative set of satellite information because of the combination of high spatial and temporal resolutions. On the one hand, the minimum pixel size of ten meters enables the assessment of heterogeneous landscapes such as crop mosaics [1,2,3] with high field size variability as well as areas with strong environmental gradients such as mountainous regions [4,5,6]. On the other hand, the high frequency of data acquisition (5 days) makes it possible to monitor vegetation seasonality in high detail [7]. Several studies have demonstrated its utility when identifying specific phenological events (e.g., start and end of growing season), which is highly relevant when evaluating the influence on the vegetation of agricultural management practices and the impact of extreme weather events [8,9,10,11,12]. This is especially interesting in regions with a Mediterranean type of climate, where the interannual variability of meteorological factors results in different vegetation adaptation strategies [13,14] and agricultural practices, such as crop rotation and sowing and harvest dates [15,16,17]. In these areas, rainfed crops, dependent on rainfall, have their main growing season in spring [18,19], whereas irrigated crops grow during summer in specific areas with infrastructure [20]. In addition, economic circumstances and political regulations increase cropping dynamics variability and, hence, make monitoring more complex [21].
The experience acquired from working with low-spatial-resolution time series, such as MODIS, highlights the potential of Sentinel-2 time series. To achieve Sentinel-2’s full potential, there are two significant limitations: (1) the computational needs due to the increasingly longer time series with frequent observations and high spatial resolution, and (2) the presence of low-quality observations and gaps in time series [22] due to the presence of snow, sensor failure, and atmospheric disturbance (e.g., clouds) [23,24]. To overcome these issues, some authors used only cloud-free images, losing temporal resolution and limiting the dataset to some few images per year [25,26,27,28]. However, developing computationally efficient methods based on temporal analysis will make it possible to fully exploit the length of the Sentinel-2 time series while maintaining a high spatial resolution.
Compositing procedures such as Maximum Value Composite (MVC) [29] are frequently used before filtering time series from sensors with daily data acquisition (e.g., MODIS); although the effective temporal frequency is reduced, it does not significantly limit the information [30,31,32].
The impact of atmospheric disturbance becomes more limiting in sensors with longer revisiting periods such as Sentinel-2 and Landsat (5 and 16 days, respectively). In these cases, the decrease in temporal resolution due to compositing methods may [32] have a significant impact on the time series information contents, so other methods need to be explored.
The Copernicus System provides easily usable quality data that make it possible to accomplish the first assessment of the information content, including the number and length of gaps. The presence of long gaps in time series makes it more difficult to uncover vegetation dynamics patterns, and it seems appropriate to explore gap-filling and filtering strategies such as temporal filters and spectral domain procedures [33,34,35,36].
The methods used depend on the noise content in the series and the type of features to enhance. Usually, temporal filters are appropriate for preserving details and short-term patterns, while spectral domain procedures are better for highlighting long-term dynamics. In this sense, Zhou et al. [37] showed that different noise distributions require specific methodologies with temporal filters working better in the tropics and subtropics and spectral domain procedures in boreal latitudes.
Temporal filters such as the Savitzky–Golay (SG) [38] and Whittaker (WHI) [39,40] filters are based on time series smoothing algorithms using moving windows and are focused on the fact that the Normalized Difference Vegetation Index (NDVI) [41] time series are consistent with the vegetation growth curve [42]. The SG, based on a least square fitting, is the most widely used method for reducing noise in vegetation index time series [43]. The parameters used are the degree of the polynomial and the moving window size, which are defined based on the knowledge of the time series characteristics [37]. Htitiou et al. [44] and Stendardi et al. [45] applied the SG filter to the annual Sentinel-2 time series, preserving the NDVI profile and obtaining effective results in the identification of phenological metrics, especially the Start of Season (SOS). The basis of the WHI is to keep a balance between the fidelity of the data and the roughness of the filtered time series. Goldberg et al. [46] and Huang et al. [47] used it to smooth the Sentinel-2 NDVI time series to develop maps for different crops using classifiers, and a 90% accuracy was obtained. Nishanta Khanal et al. [48] compared the WHI, Fast Fourier Transform (FFT), and linear-fit averaging (LFA) smoothing methods to develop land cover classification maps in Nepal using Landsat 5, 7, and 8 images. They concluded that WHI worked better with stable classes while the FFT showed better results with the rest of the classes.
Conversely, procedures that work in the frequency domain, such as Fast Fourier Transform (FFT) [49] and Harmonic Analysis (HANTS) [37], are effective for assessing periodic time series [43].
Despite the potential of Sentinel-2 vegetation index time series for assessing vegetation dynamics with high spatial detail, only few studies testing filters and reconstruction strategies have been conducted [22,50]. Selecting the appropriate process is important to fill the gaps in time series due to the lack of images or the presence of low-quality values while maintaining the dynamics [51,52].
Evaluating the behavior of the filters and comparing their fit with the original Sentinel-2 time series data and dynamics is essential [53,54] for extracting relevant and reliable information, such as phenological metrics.
The main objective of this study was to assess the information content based on Sentinel-2 NDVI time series in terms of temporal dependency in an area with a significant climatic gradient and high variation in land cover dynamics. Specifically, five scenarios of NDVI time series patterns with different land covers and environmental conditions were assessed: (1) wheat and barley, (2) maize, (3) alfalfa, (4) beech forests, and (5) pine forests.
The specific objectives were as follows:
(1)
To assess the information content in terms of temporal dependency in the raw NDVI time series and evaluate the interactions between land cover dynamics and other sources of variability.
(2)
To assess the impact of interpolation and filtering procedures on the information content, specifically changes in the dynamics and the significance of the main periodicities.
(3)
To calculate a concise indicator to evaluate the number and distribution of invalid observations due to meteorological perturbances based on the image classification scene provided by Copernicus.

2. Materials and Methods

2.1. Study Area

The study area covered the Sentinel-2 tile 30TUN (212,056.05 km2), with the upper left corner at 43°19′34″N–5°28′0″W and lower right corner at 42°21′31″N–4°5′43″W. This tile includes the provinces of León, Palencia, Cantabria, Asturias, and Burgos (Figure 1). The study period was from January 2017 to August 2023.
This area has a temperate climate with dry summer in the southern area (Koppen Csb), and it is fully humid in the northern area (Koppen Cfb) [55,56]. Precipitation is usually concentrated in autumn, winter, and spring, with an accumulated yearly rainfall of 940 mm [55]. Soils are Fluvisols and Cambisols according to the FAO-UNESCO system [57]. Forested areas constitute approximately 35% of the Sentinel-2 tile, and they are distributed at higher altitudes than croplands. The area shows a strong topographical gradient that results in high variability in short distances.

2.2. Experimental Design and Working Scheme

The work was focused on assessing the information content of Sentinel-2’s complete NDVI time series and how several pre-processing methods can enhance the information content. To accomplish a comprehensive analysis, the following criteria were proposed: (1) the presence of distinct sources of variability, (2) the use of different pre-processing methods, and (3) the use of several evaluation methodologies.
To assure variability, a complete Sentinel-2 tile with a significant climatic gradient, representative land cover types, and vegetation dynamics was chosen; data were assessed and pre-processing levels were from 2017 to 2023. In the agricultural area, the crops selected were wheat, barley, and irrigated crops alfalfa and maize. For forest species, beech and pine forests were selected.
Figure 2 shows the flowchart from the building of the time series to the filtering process. The Jupyter Notebook application was used in Python 3.9 language.

2.3. Data Sources

2.3.1. Reference Data

In agricultural areas, crops were identified based on the 2018 to 2023 PAC farmers’ declarations for each field. This information, included in the “Sistema de Información Geográfica de Parcelas Agrícolas (SIGPAC)”, was downloaded from the “Fondo Español de Garantía Agraria (FEGA)” webpage (https://www.fega.gob.es/atom/es.fega.sigpac.xml, accessed on 14 June 2023). The surface evaluated was 455.74 km2, including rainfed and irrigated crops. Forest information was obtained from the forest map of Spain at a scale of 1:50,000 [58]; the area evaluated was 872.63 and 563.82 km2 for beech and pine forests, respectively. Figure 3 shows the spatial distribution of the land covers studied. Only fields with the same vegetation type throughout the study period were evaluated.
Non-irrigated wheat and barley evolution depends on meteorological conditions while irrigated maize shows a more stable cycle during summer. Alfalfa shows a distinct variability due to cuts during the growing period. Table 1 shows the growing period for different crops in our study area. Rainfed crops such as wheat and barley have a growing period, which is highly dependent on rainfall, during spring. Alfalfa, cultivated under irrigated conditions, is cut two to six times per year, and irrigated maize grows in more stable conditions during summer [59].
Beech and pine forests were selected for this study due to their distinct dynamics, being deciduous and evergreen, respectively. This choice was aimed to enhance our understanding of the impact of interpolation and filtering in mountainous areas with a higher cloud cover.

2.3.2. Sentinel-2 Data

The Sentinel-2 Multispectral Instrument (MSI) on board the satellites, 2A and 2B, from the Copernicus program provides information in 13 bands in the optical range at spatial resolutions of 10, 20, and 60 m and every five days [60].
In this research, we used the highest resolution images of 10 m to be able to work with pure pixels. A total of 485 images from January 2017 to August 2023 for the 30TUN tile were downloaded at Level-1C: Top-Of-Atmosphere (TOA) from https://dataspace.copernicus.eu/browser/ (accessed on 5 September 2023), and then they were atmospherically corrected using the Sen2Cor V.02.09 tool (https://step.esa.int/main/snap-supported-plugins/sen2cor/, accessed on 15 October 2023) to convert them to Level-2A: Bottom-Of-Atmosphere (BOA). This atmospheric correction process generated the Scene Classification Map (SCL) used to assess data quality; thus, an observation was considered invalid when the pixel’s class was “no data”, “saturated or defective”, “dark area pixels”, “cloud shadows”, “cloud medium probability”, “cloud high probability”, and “snow”. In this study, we did not use “thin cirrus” clouds because they are partially transparent and still allow for the capture of useful information [61]. Using this classification could significantly reduce the available data, especially in regions with frequent cloud cover, like our study area in the northern region.

2.4. Compilation of NDVI and Quality Data Time Series

The NDVI was calculated from the Level-2A images (b4-red: 646–685 nm and b8-nir: 767–908 nm), and the complete time series were generated. The SCL quality band of each date was used to identify the low-quality observations and assess their frequency and distribution in the time series. In the context of time series analysis, a “gap” refers to a period of time when data are scarce, missing, or incomplete. These gaps may arise due to various factors such as technical issues in data collection or simply a lack of data availability for certain periods. Identifying and understanding these gaps is crucial to ensure the accuracy and integrity of time series analysis.
The Interpolating Efficiency Indicator (IEI) was proposed for assessing the potential information content of a pixel time series based on the number of invalid observations and the presence of long gaps in terms of length and distribution within the year. While short gaps (i.e., 1 to 3 invalid consecutive observations) may be interpolated with high reliability, the interpolation process may be limited when gaps are longer, resulting in inaccurate information, especially when they occur during significant periods within the growing cycle.
The IEI (Equation (1)) values range between 0 and 100, with higher values indicating a higher capability of the original time series to retain significant information content to be used in analyses such as phenological studies, crop dynamics analysis, and time series analysis, among others.
∀ x ∈ 1 , 16 ∪ G t 16
F I E I x =   100 − 1 / N ∑ x = 1 m f x ∗ 100   i f   F I E I ≥ 0       0   i f   F I E I < 0  
f x =   100   i f   x = N   W t x ∗ n x   i f   1 ≤ x   ≤ 16   W t 16 + x ∗ 100 N   i f   x > 16  
where the function FIEI(x) ranges from 0 to 100 and is based on f(x) and N (time series length); f(x) is a piecewise function that depends on x (gap) and Wt (the weight of x), and n (the number of gaps that range from 1 to m) for the two intervals defined in Table 2: for x in Equation (2) and N. Gt16 represents the upper limit of the x values in the calculation of FIEI. It is a list of values calculated from the expression of f(x) for x > 16 and contains each gap longer than 16 in the original time series. The weights associated with the gaps in the IEI were defined based on the understanding that larger gaps have a detrimental impact on the overall series. This implies that larger gaps require greater weights to be applied.

2.5. Time Series Interpolation and Filtering

Linear interpolation was performed to complete the time series by replacing non-valid observation values as defined in the Sentinel-2 SCL quality band.
In the second step, four classical filters were applied: Savitzky–Golay, Whittaker, Fast Fourier Transform, and Maximum Value. These filters modify the original data values to achieve their respective objectives, such as noise reduction or smoothing. The changes are a necessary consequence of the filtering processes used to extract or emphasize features of the data.
The Savitzky–Golay (SG) filter [38] is based on the least squares method and is one of the most used methods for reducing noise and smoothing remote sensing time series [62,63]. In our research, the width of the smoothing window and the degree of the smoothing polynomial were set to 7 and 2, respectively.
The Whittaker (WHI) filter [39,40] is based on the balance between fidelity to the original data (S) and roughness (R). The smoothing parameter chosen in our research was 10% of the length of the time series.
The Fast Fourier Transform (FFT) filter [49] is an algorithm that computes the Discrete Fourier transform (DFT) and its inverse efficiently [64]. The time domain data are decomposed into the frequency domain using Fourier analysis.
The Maximum Value Filter (MVF) is based on the MVC [29]. We set as a threshold value the mean calculated using a window width of seven observations. NDVI values lower or equal to the mean are substituted by the maximum value within the window width of seven observations. Finally, each observation is averaged with the previous and subsequent values.

2.6. Assessment of Similarity

The similarity between valid observations in the raw and filtered time series was estimated by two indictors: first, the Pearson Correlation Coefficient (CC) to measure the strength of the adjustment, where high CC values represent a better adjustment to the original valid data [43]; second, the Root Mean Square Error (RMSE) based on the standard deviation of the residuals in the relationship between the original and filtered valid observations; lower RMSE values represent a better adjustment [65].

2.7. Assessment of Time Series Information Content

In order to evaluate the information content of the time series, both the intra-annual and interannual dynamics in the original and processed time series were considered.

Assessment of NDVI Dynamics

Time series analysis was used to study the dynamics of the original (including noise) and the interpolated and filtered time series, including short-term, seasonal phenological patterns, and interannual variability. Three tools of time series analysis were used:
(1)
The average year was calculated by applying Buys-Ballot tables [66] to study the intra-annual dynamics of the interpolated and the filtered time series.
(2)
The Q-Ljung-Box test (L-B Q) (Equation (3)) [67] was originally defined to test model adequacy in time series. In this work, this test was used to assess the amount of significant information content in a time series by the estimation of the autocorrelation function [17,68,69]. The test was carried out under the classical null hypothesis that the time series were white noise.
The statistic for the test is as follows:
L − B   Q M = N N + 2   ∑ k = 1 M ( N − k ) − 1 r ^ k 2   ~   x 2   M   d . f .
where
N is the number of time series observations,
r ^ k is the autocorrelation coefficient at lag k,
M is the number of lags up to which the sum of autocorrelation coefficients is calculated in the Q test and defines the number of degrees of freedom (d.f.) for the X2 probability distribution. The Q-test was calculated for several time horizons (M), short-term analysis (lags 1 to 7), seasonal analysis (lags 36 and 73), and biannual analysis (lag 146).
(3)
The periodogram analysis [70] was performed to identify the presence of periodic components of the NDVI time series [71]. The presence of a seasonal component and its significance were evaluated using Fisher’s Kappa test (Fk) [72] calculated as a ratio of the maximum periodogram ordinate to the mean of all ordinates. The null hypothesis is the lack of a periodic component.
The statistical analyses were performed using Python 3 and SAS software (SAS 9.4 Software, SAS Institute Inc., Cary, NC, USA).
For the calculation of the IEI and the assessment of similarity, the complete time series from January 2017 was used. However, for the assessment of time series information content, data from the years when Sentinel-2b was operational (starting in March 2017) were more suitable. Therefore, we used the series from 2018 onwards.

3. Results

3.1. Assessment of Pixels Quality

The invalid observations and gap frequency and length due to clouds for each pixel along the study period (SCL band) are shown in Figure 4 and Figure 5. Figure 4a shows the percentage of valid observations for each pixel in the study area. Figure 4b shows the number of pixels for each number of valid observations (the far extremes are not shown). The number of missing images was 21, 1, 2, 1, and 2 in the years 2017, 2018, 2020, 2022, and 2023, respectively. Low frequencies of valid observations in buildings and water bodies are not shown in the histogram.
The spatial distribution of the percentage of valid observations shows a significant spatial coherence with two clear regions and a significant gradient related to latitude and elevation (Figure 4a). Time series with a greater number of valid observations were mainly in the southern flat areas (Koppen Cs, temperate with dry summer) while in the northern mountainous areas (Koppen, temperate, fully humid Cfb), the percentage of valid observation was lower.
This pattern is shown clearly in the number of pixels for each number of valid observations (Figure 4b), with two distinct populations corresponding to the mountainous areas in the north and the flat areas in the south, with approximately 110 to 236 and 245 to 299 valid observations, respectively.
Figure 5a,b show the length of the longest gap in the time series and season when it occurs, both on a pixel basis. Table 3 shows the total area in which the longest gap in the time series occurs at a given season, including the statistics of the length of the gaps for each season. The maps show a consistent spatial distribution; the longest gaps with more than 16 consecutive observations make up around 25% of the tile, and most of these gaps are in the mountainous region situated in the northern part of the tile, occurring more frequently in winter and secondarily in spring. Most of the surface shows the longest gap between 7 and 9 observations (36% of the tile), and the range of 10 to 12 observations represents 26%; these are in the southern part of the tile.

3.2. Interpolating Efficiency Indicator (IEI)

Figure 6a shows the IEI for Tile 30TUN with high values in light yellow and low values in dark purple and black; Figure 6b shows the frequency distribution.
In the flat areas, values of up to 96% were observed, with the gap length ranging from 1 to 8 and having a smaller impact in the interpolation and filtering process than in the mountainous areas, where maximum values were around 78%. In 40.9% of the surface (4928.13 km2), IEI values ranged between 89% and 92%, indicating that the interpolation and filtering process will be able to retain more information content for subsequent analysis.
Table 4 shows the IEI values by season and land use, both annually and seasonally. IEI for beech forests, located in Cfb climate areas, presents the lowest annual values, whereas maize, in Cs areas, shows the highest values; the pine forests in the transitional zones show intermediate values. IEI is lowest in winter and highest in summer across all land covers.

3.3. Examples of Interpolation and Filtering Data

Examples of the interpolated and filtered NDVI time series are shown in Figure 7 and Figure 8 for crops and forests, respectively; time series and the average year corresponding to each land cover are shown. For all types of vegetation, the SG is the one that best fits the interpolated series, while the MVF seems to overestimate the NDVI values.
Figure 9 shows an example of the evolution of one image (10/06/20) after the interpolation and Whittaker filtering processes. Figure 9a (10/06/20), Figure 9b (15/06/20), and Figure 9c (20/06/20) show three consecutive dates of the raw NDVI time series, where invalid values according to the SCL band are indicated in purple; these low-quality values were replaced after linear interpolation (Figure 9d) and filtering (Figure 9e), showing an improvement in the spatial consistency of the scene.

3.4. Similarity between Original and Filtered Time Series in Pure Plots

Pure plots of different vegetation species were selected to evaluate the performance of the filters in terms of the Correlation Coefficient (CC) and the Root Mean Square Error (RMSE) (Table 5). In the crops, alfalfa showed the worst values for both indicators (CC = 0.83 and RMSE = 0.15) and the highest variability among filters. The SG filter showed the minimum RMSE and the maximum CC value, whereas the MVF showed the worst performance with the opposite pattern. For the forest species, beech forests have higher CC values than pines but similar values of RMSE.

3.5. Assessment of NDVI Time Series Information Content and Dynamics

The information content in the NDVI time series from the raw to the filtered data was estimated by the Q-Ljung-Box test, calculated for the short term (lags 1, 2, 3, 4, 5, 6, and 7) and at lags 36 (seasonal), 73 (annual), and 146 (interannual).
To clearly illustrate the results, Q-test values of raw NDVI time series (including low-quality observations) are shown in Table 6. Figure 10 (short term) and Figure 11 (seasonal and annual terms) show the Q-Ljung-Box test values at different lags (1 to 146) for the raw data and the interpolated and filtered time series. Q-values in all land covers were higher than the respective critical values at the 5% significance level for lags 1 to 146; consequently, the null hypothesis of white noise could be rejected in all cases. The largest Q-value increase was after the interpolation, while the filtering procedure resulted in less improvement; the maximum increase was observed when applying the MVF (i.e., wheat Q73 = 5043) and the lowest when using the SG (i.e., wheat Q73 = 4058). Beech forest showed the highest increase while pine forest the lowest at all lags.
Temporal dependency until lag 7 (35 days) could be divided in four levels: (1) maize, (2) wheat and barley, (3) alfalfa and beech forest and, (4) pine forest. The patterns of the seasonal and annual terms (Q36, Q73, and Q146) were similar to those of the short term, except for the beech forest, where NDVI time series presented a high information content compared to the short term (i.e., Q73/Q1 = 28).
Maize presented the highest values for the raw data in the short and long terms, while beech forests reached the highest values when filtered. The lowest values were for pine forests. It was observed that after interpolation, at seasonal and annual horizons, there were three groups with similar patterns: beech and maize with high values, wheat and barley with medium values, and alfalfa and pines with low values.
Table 7 shows the percentage of the surface for each land cover with the NDVI time series main periodicities at the periods 73 (annual component or phenological cycle), 36 (6 months), and 203 (representing the presence of a trend). The rest of the possible periodicities are not shown because none of them represented more than 1% of the surface. The periodicity at the annual period (period 73) was the most widely extended across all vegetation types but with different surface percentages. With the raw series, more than 90% of the surface had the main periodicity at one year in maize, barley, and beech forest, while in alfalfa and pine forest, the annual cycle was dominant only in 74% of the surface. In alfalfa, 13% of the time series with the maximum periodicity at 203 may indicate the presence of a trend while in the pine forest’s 5% of the surface shows the presence of a second intra-annual cycle.
In the interpolated and filtered time series, the percentage of the surface with an annual cycle decreased in all cases except in the beech forest case, where it reached 99%. In the rest of cases, the surface with the semi-annual cycle increased; in addition, in alfalfa and pine, there was a significant increase in the surface with a dominant cycle at a period of 203.
Figure 12 shows the average Fk values for the periodicity of 73 for each vegetation type. The Fk values of the raw time series are represented on the left axis, and the percentage of increase of the Fk after NDVI time series interpolation and filtering is shown on the right axis.
In the raw time series, maize showed the highest Fk values; however, when interpolating and filtering, the percentage of increase was lower than in the other land uses (wheat, barley, and beech forest). Beech forest also showed high Fk values and a high increase. In wheat and barley, a better identification of the cycle was obtained when applying the filters after interpolation. Alfalfa also showed a high increase but not as much as other crops; this was probably due to cuts during the vegetative cycle. Pine forest showed the lowest improvement. Generally, the MVF and WHI filters showed higher improvements, whereas the SG filter presented the lowest.

4. Discussion

The information content in terms of the temporal dependency of Sentinel-2 NDVI time series from 2018 to 2023 for a tile (1,200,000 ha) in central Spain was evaluated. Five scenarios with different climate and distinct phenological patterns were assessed: (1) wheat and barley, (2) maize, (3) alfalfa, (4) beech forest, and (5) pine forest. In these scenarios, the time series information content was quantified at three processing levels: raw, interpolated, and filtered. The IEI summarized information on atmospheric conditions based on the frequency of invalid observations and number and length of gaps due to clouds. For clarity purposes, general patterns are discussed first, and then specific issues for each land cover are addressed.
In the raw time series, Q-test values were higher than the test critical values in all scenarios, indicating that NDVI Sentinel-2 time series were not white noise and contained significant information related to vegetation dynamics. In more than 87% of the surface of wheat, barley, maize, and beech, and around 74% of alfalfa and pine forest, one seasonal cycle per year was found (Table 7, column Raw); this was expected, given that only fields with the same land cover in five years were evaluated.
After interpolation, the information content (i.e., Q and Fk values) increased in all land covers and all terms (Figure 10b, Figure 11b and Figure 12), indicating an increase in time series coherence; particularly, the increase in the Q1 to Q7 slope (i.e., 35 days, Figure 10b) demonstrated a higher consistency in the short term, which could be useful for phenological studies. While Chen et al. [73] used linear interpolation to reduce noise, maintaining the time series dynamics, in this study, the surface with one annual cycle decreased (Table 7, column Interpolated), except in beech forest. These results suggest that the filtering procedure emphasized certain dynamics not identified previously, such as a subtle second cycle due to weeds.
After filtering, crops showed better adjustment than forests (Table 5) located in cloudier areas with lower IEI values. The SG filter showed the highest similarity between the filtered values and the original high-quality ones, especially in maize (0.98), due to the high IEI value and information content in the interpolated time series. These results agree with those of Li et al. [43] who compared different filters in a forestry district in China with MODIS data in five years and obtained the highest CC and lowest RMSE values with the SG filter. In addition, Fontana et al. [74] showed excellent results with the SG filter in AVHRR images evaluating grassland phenology. The MVF and WHI filters showed a high smoothing capability (visual inspection, Figure 7 and Figure 8); however, the MVF overestimated NDVI values, which is in accordance with other studies [75,76,77]. The lowest CC value (0.68) was found when applying the MVF filter to the pine forest time series in agreement with the intermediate IEI values and lowest information content. The WHI filter showed slightly better results in adjustment than the MVF, which did not agree with Liu et al. [78] who found that the WHI filter did not have a good overall performance in filling the gaps of 8-day MODIS time series on a global scale. On the other hand, the WHI filter was faster than the others in computation time, as other authors have reported [79,80]. The FFT showed better performance in capturing small decreases in alfalfa NDVI time series due to cuts during summer (Figure 7).
The temporal dependency (i.e., Q values) of NDVI time series was slightly higher in all terms than in the interpolated data, especially when using the MVF and WHI, which showed similar Q values in the short term in all land covers, except in pine. This is consistent with their higher smoothing capability and the lower values in similarity measures. In the annual term (Figure 11c–f), they also showed the highest Q values, especially in beech forest.
When applying the MVF filter in pine forest and maize, the surface where the main periodicity was the annual cycle increased (Table 7), as well as the significance of the Fk test value, again indicating the highest smoothing capability of this filter. Figure S1 shows the surface where the maximum ordinate after the interpolating and filtering procedures changed. Patel and Oza [81] were able to identify double cycles in their time series through effective filtering techniques.
In the spatial context (Figure 9), NDVI images generated during the interpolation and filtering process represented the field structure with spatial consistency and detail due to the high spatial resolution of Sentinel-2 time series. In this regard, Reza et al. [82] also achieved good results when filtering contaminated images using the HANTS filter. Similar results were obtained by other researchers in areas with regular NDVI patterns, such as crops or urban areas [3,83,84].
  • Case 1. Wheat and Barley. Rainfed crops grown in winter–spring in a Mediterranean climate.
These rainfed crops presented significant random short-term variability in NDVI time series (i.e., lower Q1–Q7 values and slope) due to highly variable atmospheric conditions and frequent cloud presence (IEI values: 86% in winter and 91% in spring, see Table 4).
The temporal dependency under seasonal, annual, and interannual dynamics (i.e., Q36, Q73, and Q146 values) was lower than in maize and beech, with a low ratio between memory in one year and in the short term (Q73/Q1 = 11) suggesting the low stability of the annual cycle and the presence of high interannual variability. This was expected in Mediterranean regions where temperature and precipitation show a high interannual variability [85], which could be increased by farmers who decide on the date of sowing every year depending on precipitation and soil status [86].
Barley presented a predominant annual cycle in a larger surface than wheat, with 95% vs. 87%, also showing higher accumulated information in the annual terms (i.e., Q73, 787.05 vs. 718.86) and slightly higher Fk values (19.23 vs. 17.30). In the wheat area, 12.16% of the time series did not show any defined annual cycle (i.e., class “others”) versus 3.83% in barley fields, suggesting a higher stability of the growing cycle in barley than wheat [87].
Interpolated time series significantly increased Q values in the short term, and this result may be important in cases involving identifying phenological parameters, such as in Zhu et al. [88] who observed that time series dynamics and phenological stages remained constant after linear interpolation. In our case, after interpolation, there was a small percentage of time series with significant annual dynamics that were re-allocated to the classes “period 36” (two intra-annual cycles) and “others”. It is probable that there was an enhancement of subtle information content related to the presence of weeds or short-cycle crops used as forage that were previously masked. Figures S2 and S3 show examples of time series with different main periodicities.
  • Case 2. Maize. An irrigated crop grown in summer in a Mediterranean climate.
The raw time series of maize showed the highest accumulated information in all terms, with 91% of the surface displaying an annual cycle; the average maize Fk value (Fk = 62.10) was the highest among land covers, which agrees with the high Q value in one year (Q73 = 2050.96) and in two years (Q146) (Figure 11). These results are consistent with the stable and repetitive annual growing period in the irrigated areas during the Mediterranean summer, with optimal conditions of temperature and stable moisture availability. The high IEI values in the maize fields during summer (97.78%) explain the fast increase in Q values from 1 to 7 lags, indicating a low high-frequency variability due to clouds, which results in a regular NDVI evolution within the growing cycle.
After interpolation, the temporal dependency increase was the lowest among land covers (i.e., Qint/Qraw ≅ 2), likely because it was already high from the outset. Also, the surface with two intra-annual cycles increased slightly, suggesting that the noise elimination made it possible to uncover subtle variability sources such as the one due to the growth of weeds before preparing the soil between the growing seasons; this is shown in the exploratory analyses (Figure S4a). In addition, a significant surface percentage with other cycles suggested the presence of two crops in some years, such as vetch until the end of July and then, maize, resulting in a double greening cycle; this was not stated by the farmer, so it probably was not a common practice every year, and, therefore, the periodicity could not be defined.
After filtering, there was a slight increase in temporal dependency for all lags and in Fk values, with the lowest values for the SG filter. The MVF filter showed the highest increase at the interannual scale, together with an increase in the area with an annual cycle (90.08%), suggesting that the strong smoothing capability of this filter emphasized the annual cycle over seasonal patterns.
  • Case 3. Alfalfa. A permanent irrigated crop (5-6 years) in a Mediterranean climate.
Alfalfa showed low Q values, both in the short, seasonal, and annual terms, despite growing in summer with no water limitations and optimal atmospheric conditions (IEIsummer = 97.58%). Crop management including several cuts through the growing period may explain the short-term high variability (Figure 7g) that could hide the annual cycle information (Q73/Q1 = 10.16); this was corroborated by the fact that only 74.00% of the surface of alfalfa was shown in the annual cycle, with low significance (Fk value = 19.73). The other 26% of the surface was mainly distributed between the classes “others” (11.71%) and “period 203” (13.12%), suggesting the presence of a trend.
After interpolation, there was still a low amount of cumulated information in all terms, suggesting that noise due to atmospheric effects was not the main cause of high frequency variability, but rather the decrease in NDVI due to the cuts during the growing period. In addition, there was a significant increase in the surface percentage with a possible trend (33.56%); since the interpolation process tends to emphasize subtle patterns and decrease high-frequency variability, and the alfalfa crop usually lasts five to six years, this result suggested that these time series corresponded to fields where alfalfa was at the first growing stages in 2018 (Figure 13) or at the last in 2023. This hypothesis is consistent with Figure S5 in the Supplementary Material, which shows these time series located in homogenous fields, indicating that it was due to management.
The rest of the alfalfa surface (17%) is divided between the classes “period 36” and “others”, probably due to the high variability introduced by the cuts. The MVF and WHI significantly increased temporal dependency in the short term; however, the similarity indicators (Table 5) showed low values, which was probably due to the low adjustment of filtered NDVI time series to the original NDVI values when alfalfa is cut (Figure 7h).
  • Case 4. Beech forest. A deciduous forest in a temperate climate.
Beech forest showed significant but remarkably low values of temporal dependency in the short term (Q1–Q7), indicating an important source of high-frequency random variability. Annual information was highly relevant in relation to the short term (i.e., Q73/Q1 = 28); in addition, Q146 was higher than Q73 (i.e., Q146/Q73 = 1.81), suggesting the stability of inter-annual dynamics. Growing in the northern area with high cloudiness (i.e., IEI = 74.90%) resulted in high short-term variability. However, the indicators on the annual and inter-annual scales suggested that the distinct phenological pattern of the Fagus species growing in summer (IEIsummer = 92.74%) provided enough consistent information to override short-term variability. Similar results were found by Recuero et al. [89] who were able to map maize crop dynamics in Ecuador despite the high presence of clouds due to the significant annual cycle.
Ninety one percent of the beech surface presented one seasonal cycle with a high Fk value (Figure 12). On the other hand, 8.80% of the surface, specifically located in the highest latitudes, did not show a dominant periodicity; in this area, the average Q value of raw time series was significantly lower than in the whole beech forest areas (Q73 = 98 vs. Q73 = 677.62) probably due to the extremely high presence of clouds that resulted in an IEIwinter = 56.10%. Still, even in this area, the information content in the annual and bi-annual terms showed Q values greater than the critical.
After interpolation, temporal dependency increased, especially at seasonal, annual, and interannual lags (i.e., Q36, Q73, and Q146), reaching the highest values among land covers, and suggesting that these periodic components were highly significant and regular in intensity. In addition, the percentage of the area with one annual cycle increased to almost 99% (Table 7), and the significance of the annual cycle increased significantly (i.e., Fk increase percentage, Figure 12). This suggests that the interpolation process resulted in a relevant enhancement of significant annual information by diminishing the short-term random variability derived from the frequent cloud coverage.
  • Case 5. Pine forests. An evergreen forest in a transition area from mediterranean to temperate climate.
The lowest Q values (in the limit of the critical values) in the short term indicate the presence of high-frequency variability; this is consistent with the relatively high NDVI values in the year of evergreen species and the fast response to changes in environmental conditions; fast NDVI especially increases when optimal conditions of temperature are reached [90]. In addition, the pine forest was distributed in a wide area with a significant climate gradient and an average annual IEI value of 88.15% due to cloudy conditions. This, together with their irregular NDVI dynamics, entailed the lowest Fk values, suggesting a lack of a stable annual cycle.
In the interpolated time series, the Q values were still in the limit of significance, despite the increase in information content. In the beech forest, the short-term variability due to clouds was effectively ameliorated; in the pine forest, the variability introduced by physiological responses resulted in low-temporal-dependency information. In the annual term, the Q values remained low due to the unstable phenological dynamics, which has also been reported by other authors; specifically, Aragones et al. [91] found high interannual variability when working with a 17-year MODIS NDVI data set. Furthermore, in our case, the Fk values even decreased after interpolation and SG filtering, indicating an even lower significance in the annual cycle. In addition, in 34% of the surface associated with the annual cycle’s “period 73”, other seasonal components were uncovered (Table 7), such as the presence of a trend of “period 203” and “other” cycles. The NDVI time series showing these changes were distributed throughout all the pine forest areas and presented high local variability (see Figure S7). This is consistent with the significant differences found by Aragones et al. [91] among pine species in terms of phenological patterns. The MVF filter showed an increase in the surface with an annual cycle and the presence of trends, which is probably related to the high smoothing capability of this filter. Mantas et al. [92] demonstrated the capability of Sentinel-2 imagery and the importance of identifying trends in pine forests to analyze forest decline due to disease.

5. Conclusions

This work quantitatively demonstrates that Sentinel-2 NDVI time series contains relevant information on land cover dynamics based on temporal dependency (Q-test) in five land covers, even in areas with significant atmospheric disturbance. It was also possible to identify the main periodic components and trends by the classic Fk test.
Temporal dependency increased significantly in all terms after interpolation and filtering; in addition, the surface with an annual cycle decreased, making it possible to find fields with other cycles as main periodicity, which was not possible with the raw time series.
The Q and Fk tests, together with the similarity indicators CC and RMSE, showed that a better fit (i.e., high similarity values) does not always imply an improvement in the information content in terms of dynamics assessment.
The IEI adequately represented atmospheric disturbance in the region on the annual and seasonal scales, and it was highly related to short-term variability. Furthermore, it was demonstrated that the enhancement derived from interpolating and filtering was more notable in time series located in areas with lower IEI values. It is expected that the IEI could provide the first insight into the final information content based on a quality dataset.
Significant differences in information content patterns among scenarios were found due to interactions between vegetation dynamics, cloud disturbances, and specific management practices.
Specific insights into the information content and the impact of interpolation and filtering for the scenarios assessed can be summarized as follows:
In rainfed wheat and barley, highly variable cloud cover (low IEI) in winter and spring, together with high dependence on precipitation and temperature, resulted in significant high-frequency and interannual variabilities. Interpolation and filtering reduced short-term variability; however, the impact on annual and inter-annual dynamics was lower.
Irrigated crops, such as maize, growing in summer with a significant annual cycle and low cloud coverage expressed the largest amount of information content, so that interpolating and filtering did not result in a significant increase in information (Q-values); however, this slight increase and the use of the Fk test made it possible to find fields where the main periodicity corresponded to subtle cycles such as those caused by weeds. In addition, reducing high short-term variability could facilitate phenological analyses.
In irrigated alfalfa, a multiyear crop, all filters, especially the MVF, made it possible to identify fields with trends, caused by low NDVI values at the beginning or at the end of the study period, due to crop implantation or removal, respectively. On the other hand, the MVF and WHI filters masked the small NDVI decreases due to cuts during summer.
In beech forest, which is temperate deciduous, diminishing the cloud disturbance through interpolating and filtering combined with the strong and stable annual cycle could facilitate an accurate dynamics assessment. Thus, in these areas, the Sentinel-2’s interpolated time series could provide a significant amount of ecological information.
In the pine forest, which is evergreen coniferous in Mediterranean areas, the combination of high physiological short-term variability and variable cloud coverage resulted in complex time series. Despite these difficulties, after interpolation and filtering, the Fk test was able to identify areas with trends that would not appear in the raw time series, and this would be useful for decaying studies.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs16162980/s1.

Author Contributions

Conceptualization, C.S. and A.P.-O.; methodology, C.S., J.L., G.G., D.M., L.R., and A.P.-O.; software, C.S., G.G., D.M., L.R., and A.B.-S.; formal analysis, C.S., A.P.-O., V.C., and J.L.; investigation, C.S., V.C., G.G., D.M., L.R., A.B.-S., J.L., I.d.l.C., and A.P.-O.; resources, C.S., J.L., and A.P.-O.; writing—original draft preparation, C.S., V.C., L.R., A.B.-S., J.L., and A.P.-O.; writing—review and editing, C.S., V.C., A.P.-O., and J.L.; supervision, A.P.-O.; funding acquisition, A.P.-O. and C.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by a predoctoral scholarship awarded by the Community of Madrid and Quasar Science Resources, S.L. (No. IND2020/AMB-17747). It was also conducted in the framework of the I+D+i Spanish National Project PID2020-115509RB-I00 project (INFOLANDYN) funded by the Ministerio de Ciencia e Innovación of Spain MCIN/AEI/ 10.13039/501100011033. VC was supported by a postdoctoral Juan de la Cierva fellowship (FJC2021-046735-I) funded by the Spanish Ministerio de Ciencia e Innovación MCIN/AEI/ 10.13039/501100011033 and by the European Union ≪NextGenerationEU≫/PRTR. LR was supported by a UPM grant as part of the national research program Recualificación del Sistema Universitario Español, financed by the Recovery and Resilience Package—NextGenerationEU (European Commission).

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

We would like to thank the European Space Agency (ESA) for the images provided, which were necessary for this study. The authors would also like to acknowledge the participation of Quasar Science Resources, S.L., in this study.

Conflicts of Interest

César Sáenz and Ignacio de la Calle were employed by the company Quasar Science Resources, S.L. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Harmanny, K.S.; Malek, Ž. Adaptations in Irrigated Agriculture in the Mediterranean Region: An Overview and Spatial Analysis of Implemented Strategies. Reg. Environ. Chang. 2019, 19, 1401–1416. [Google Scholar] [CrossRef] [PubMed]
  2. Gumma, M.K.; Tummala, K.; Dixit, S.; Holecz, F.; Kolli, R.N.; Whitbread, A.M.; Krishna, M.; Tummala, K.; Dixit, S. Crop Type Identification and Spatial Mapping Using Sentinel-2 Satellite Data with Focus on Field-Level Information. Geocarto Int. 2020, 37, 1833–1849. [Google Scholar] [CrossRef]
  3. Tran, K.H.; Zhang, H.K.; Mcmaine, J.T.; Zhang, X.; Luo, D. 10 m Crop Type Mapping Using Sentinel-2 Reflectance and 30 m Cropland Data Layer Product. Int. J. Appl. Earth Obs. Geoinf. 2022, 107, 102692. [Google Scholar] [CrossRef]
  4. Rujoiu-Mare, M.-R.; Olariu, B.; Mihai, B.-A.; Nistor, C.; Săvulescu, I. Land Cover Classification in Romanian Carpathians and Subcarpathians Using Multi-Date Sentinel-2 Remote Sensing Imagery. Eur. J. Remote Sens. 2017, 50, 496–508. [Google Scholar] [CrossRef]
  5. Wakulinska, M.; Marcinkowska-Ochtyra, A. Multi-Temporal Sentinel-2 Data in Classification of Mountain Vegetation. Remote Sens. 2020, 12, 2696. [Google Scholar] [CrossRef]
  6. Wang, M.; Li, M.; Wang, F.; Ji, X. Exploring the Optimal Feature Combination of Tree Species Classification by Fusing Multi-Feature and Multi-Temporal Sentinel-2 Data in Changbai Mountain. Forests 2022, 13, 1058. [Google Scholar] [CrossRef]
  7. Solano-Correa, Y.T.; Bovolo, F.; Member, S.; Bruzzone, L.; Fernández-prieto, D. A Method for the Analysis of Small Crop Fields in Sentinel-2 Dense Time Series. IEEE Trans. Geosci. Remote Sens. 2020, 58, 2150–2164. [Google Scholar] [CrossRef]
  8. Pan, L.; Xia, H.; Yang, J.; Niu, W.; Wang, R.; Song, H.; Guo, Y.; Qin, Y. Mapping Cropping Intensity in Huaihe Basin Using Phenology Algorithm, All Sentinel-2 and Landsat Images in Google Earth Engine. Int. J. Appl. Earth Obs. Geoinf. 2021, 102, 102376. [Google Scholar] [CrossRef]
  9. Amin, E.; Belda, S.; Pipia, L.; Szantoi, Z.; El Baroudy, A.; Moreno, J.; Verrelst, J. Multi-Season Phenology Mapping of Nile Delta Croplands Using Time Series of Sentinel-2 and Landsat 8 Green LAI. Remote Sens. 2022, 14, 1812. [Google Scholar] [CrossRef]
  10. Dimov, D.; Uhl, J.H.; Löw, F.; Negash, G. Sugarcane Yield Estimation through Remote Sensing Time Series and Phenology Metrics. Smart Agric. Technol. 2022, 2, 100046. [Google Scholar] [CrossRef]
  11. Kosczor, E.; Forkel, M.; Hernández, J.; Kinalczyk, D.; Pirotti, F.; Kutchartt, E. Assessing Land Surface Phenology in Araucaria-Nothofagus Forests in Chile with Landsat 8/Sentinel-2 Time Series. Int. J. Appl. Earth Obs. Geoinf. 2022, 112, 102862. [Google Scholar] [CrossRef]
  12. Ye, J.; Bao, W.; Liao, C.; Chen, D.; Hu, H. Corn Phenology Detection Using the Derivative Dynamic Time Warping Method and Sentinel-2 Time Series. Remote Sens. 2023, 15, 3456. [Google Scholar] [CrossRef]
  13. Del Pozo, A.; Brunel-saldias, N.; Engler, A.; Ortega-farias, S.; Acevedo-opazo, C.; Lobos, G.A.; Jara-rojas, R.; Molina-montenegro, M.A. Climate Change Impacts and Adaptation Strategies of Agriculture in Mediterranean-Climate Regions (MCRs). Sustainability 2019, 11, 2769. [Google Scholar] [CrossRef]
  14. Aguilera, E.; Díaz-gaona, C.; García-laureano, R.; Reyes-palomo, C.; Guzmán, G.I.; Ortolani, L.; Sánchez-rodríguez, M.; Rodríguez-estévez, V. Agroecology for Adaptation to Climate Change and Resource Depletion in the Mediterranean Region. A Review. Agric. Syst. 2020, 181, 102809. [Google Scholar] [CrossRef]
  15. Ryan, J.; Singh, M.; Pala, M. Long-Term Cereal-Based Rotation Trials in the Mediterranean Region: Implications for Cropping Sustainability. Adv. Agron. 2008, 97, 273–319. [Google Scholar]
  16. Jacobsen, S.; Jensen, C.R.; Liu, F. Field Crops Research Improving Crop Production in the Arid Mediterranean Climate. Field Crops Res. 2012, 128, 34–47. [Google Scholar] [CrossRef]
  17. Recuero, L.; Wiese, K.; Huesca, M.; Cicuéndez, V.; Litago, J.; Tarquis, A.M.; Palacios-orueta, A. Fallowing Temporal Patterns Assessment in Rainfed Agricultural Areas Based on NDVI Time Series Autocorrelation Values. Int. J. Appl. Earth Obs. Geoinf. 2019, 82, 101890. [Google Scholar] [CrossRef]
  18. Yang, C.; Fraga, H.; Van Ieperen, W.; Trindade, H.; Santos, J.A. Effects of Climate Change and Adaptation Options on Winter Wheat Yield under Rainfed Mediterranean Conditions in Southern Portugal. Clim. Chang. 2019, 154, 159–178. [Google Scholar] [CrossRef]
  19. Oweis, T.; Zhang, H.; Pala, M. Water Use Efficiency of Rainfed and Irrigated Bread Wheat. Agron. J. 2000, 92, 231–238. [Google Scholar] [CrossRef]
  20. Katerji, N.; Mastrorilli, M.; Rana, G. Water Use Efficiency of Crops Cultivated in the Mediterranean Region: Review and Analysis. Eur. J. Agron. 2008, 28, 493–507. [Google Scholar] [CrossRef]
  21. Pradipta, A.; Soupios, P.; Kourgialas, N.; Doula, M.; Dokou, Z.; Makkawi, M.; Alfarhan, M.; Tawabini, B.; Kirmizakis, P. Remote Sensing, Geophysics, and Modeling to Support Precision Agriculture—Part 2: Irrigation Management. Water 2022, 14, 1157. [Google Scholar] [CrossRef]
  22. Liu, X.; Ji, L.; Zhang, C.; Liu, Y. A Method for Reconstructing NDVI Time-Series Based on Envelope Detection and the Savitzky- Golay Filter. Int. J. Digit. Earth 2022, 15, 553–584. [Google Scholar] [CrossRef]
  23. Pelletier, C.; Valero, S.; Inglada, J.; Champion, N.; Marais Sicre, C.; Dedieu, G. Effect of Training Class Label Noise on Classification Performances for Land Cover Mapping with Satellite Image Time Series. Remote Sens. 2017, 9, 173. [Google Scholar] [CrossRef]
  24. Kong, D.; Zhang, Y.; Gu, X.; Wang, D. A Robust Method for Reconstructing Global MODIS EVI Time Series on the Google Earth Engine. ISPRS J. Photogramm. Remote Sens. 2019, 155, 13–24. [Google Scholar] [CrossRef]
  25. Hemmerling, J.; Pflugmacher, D.; Hostert, P. Mapping Temperate Forest Tree Species Using Dense Sentinel-2 Time Series. Remote Sens. Environ. 2021, 267, 112743. [Google Scholar] [CrossRef]
  26. Abdollahi, A.; Liu, Y.; Pradhan, B.; Huete, A.; Dikshit, A.; Tran, N.N. Short-Time-Series Grassland Mapping Using Sentinel-2 Imagery and Deep Learning-Based Architecture. Egypt. J. Remote Sens. Space Sci. 2022, 25, 673–685. [Google Scholar] [CrossRef]
  27. Wang, Y.; Fang, S.; Zhao, L.; Huang, X.; Jiang, X. Parcel-Based Summer Maize Mapping and Phenology Estimation Combined Using Sentinel-2 and Time Series Sentinel-1 Data. Int. J. Appl. Earth Obs. Geoinf. 2022, 108, 102720. [Google Scholar] [CrossRef]
  28. Marino, S. Assessing the Agronomic Subfield Variability by Sentinel-2 NDVI Time-Series and Landscape Position. Agronomy 2023, 13, 44. [Google Scholar] [CrossRef]
  29. Holben, B.N. Characteristics of Maximum-Value Composite Images from Temporal AVHRR Data. Int. J. Remote Sens. 1986, 7, 1417–1434. [Google Scholar] [CrossRef]
  30. Narasimhan, R.; Stow, D. Daily MODIS Products for Analyzing Early Season Vegetation Dynamics across the North Slope of Alaska. Remote Sens. Environ. 2010, 114, 1251–1262. [Google Scholar] [CrossRef]
  31. Nagai, S.; Saitoh, T.M.; Suzuki, R.; Nasahara, K.N.; Lee, K.; Son, Y.; Muraoka, H. The Necessity and Availability of Noise-Free Daily Satellite-Observed NDVI during Rapid Phenological Changes in Terrestrial Ecosystems in East Asia. For. Sci. Technol. 2011, 7, 174–183. [Google Scholar] [CrossRef]
  32. Zeng, L.; Wardlow, B.; Hu, S.; Zhang, X.; Zhou, G.; Peng, G.; Xiang, D.; Wang, R.; Meng, R.; Wu, W. A Novel Strategy to Reconstruct NDVI Time-Series with High Temporal Resolution from MODIS Multi-Temporal Composite Products. Remote Sens. 2021, 13, 1397. [Google Scholar] [CrossRef]
  33. Yang, Z.; Li, J.; Shen, Y.; Miao, H.; Yan, X. A Denoising Method for Inter-Annual NDVI Time Series Derived from Landsat Images. Int. J. Remote Sens. 2018, 39, 3816–3827. [Google Scholar] [CrossRef]
  34. Yu, W.; Li, J.; Liu, Q.; Zhao, J.; Dong, Y.; Zhu, X.; Lin, S.; Zhang, H.; Zhang, Z. Gap Filling for Historical Landsat NDVI Time Series by Integrating Climate Data. Remote Sens. 2021, 13, 484. [Google Scholar] [CrossRef]
  35. Gao, F.; Masek, J.; Schwaller, M.; Hall, F. On the Blending of the Landsat and MODIS Surface Reflectance: Predicting Daily Landsat Surface Reflectance. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2207–2218. [Google Scholar] [CrossRef]
  36. Zhu, X.; Chen, J.; Gao, F.; Chen, X.; Masek, J.G. An Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model for Complex Heterogeneous Regions. Remote Sens. Environ. 2010, 114, 2610–2623. [Google Scholar] [CrossRef]
  37. Zhou, J.; Jia, L.; Menenti, M.; Gorte, B. On the Performance of Remote Sensing Time Series Reconstruction Methods—A Spatial Comparison. Remote Sens. Environ. 2016, 187, 367–384. [Google Scholar] [CrossRef]
  38. Savitzky, A.; Golay, M.J.E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 1964, 36, 1627–1639. [Google Scholar] [CrossRef]
  39. Whittaker, E. On a New Method of Graduation. Proc. Edinb. Math. Soc. 1923, 41, 63–75. [Google Scholar] [CrossRef]
  40. Eilers, P.H.C. A Perfect Smoother. Am. Chem. Soc. 2003, 75, 3631–3636. [Google Scholar] [CrossRef]
  41. Tucker, C.J. Red and Photographic Infrared Linear Combinations for Monitoring Vegetation. Remote Sens. Environ. 1979, 8, 127–150. [Google Scholar] [CrossRef]
  42. Hird, J.N.; Mcdermid, G.J. Noise Reduction of NDVI Time Series: An Empirical Comparison of Selected Techniques. Remote Sens. Environ. 2009, 113, 248–258. [Google Scholar] [CrossRef]
  43. Li, S.; Xu, L.; Jing, Y.; Yin, H.; Li, X.; Guan, X. High-Quality Vegetation Index Product Generation: A Review of NDVI Time Series Reconstruction Techniques. Int. J. Appl. Earth Obs. Geoinf. 2021, 105, 102640. [Google Scholar] [CrossRef]
  44. Htitiou, A.; Boudhar, A.; Lebrini, Y.; Hadria, R.; Lionboui, H.; Elmansouri, L.; Tychon, B.; Benabdelouahab, T. The Performance of Random Forest Classification Based on Phenological Metrics Derived from Sentinel-2 and Landsat 8 to Map Crop Cover in an Irrigated Semi-Arid Region. Remote Sens. Earth Syst. Sci. 2019, 2, 208–224. [Google Scholar] [CrossRef]
  45. Stendardi, L.; Karlsen, S.R.; Niedrist, G.; Gerdol, R.; Zebisch, M.; Rossi, M.; Notarnicola, C. Exploiting Time Series of Sentinel-1 and Sentinel-2 Imagery to Detect Meadow Phenology in Mountain Regions. Remote Sens. 2019, 11, 542. [Google Scholar] [CrossRef]
  46. Goldberg, K.; Herrmann, I.; Hochberg, U.; Rozenstein, O. Generating Up-to-Date Crop Maps Optimized for Sentinel-2 Imagery in Israel. Remote Sens. 2021, 13, 3488. [Google Scholar] [CrossRef]
  47. Huang, X.; Huang, J.; Li, X.; Shen, Q.; Chen, Z. Early Mapping of Winter Wheat in Henan Province of China Using Time Series of Sentinel-2 Data. GIsci. Remote Sens. 2022, 59, 1534–1549. [Google Scholar] [CrossRef]
  48. Khanal, N.; Matin, M.A.; Uddin, K.; Poortinga, A. A Comparison of Three Temporal Smoothing Algorithms to Improve Land Cover Classification: A Case Study from NEPAL. Remote Sens. 2020, 12, 2888. [Google Scholar] [CrossRef]
  49. Cooley, B.J.W.; Tukey, J.W. An Algorithm for the Machine Calculation Complex Fourier Series. Math. Comput. 1965, 19, 297–301. [Google Scholar] [CrossRef]
  50. Peng, N.; Yang, S.; Tao, Y.; Zhai, D.; Fan, W.; Liu, Q. Lai Time Series Reconstruction from Sentinel-2 Imagery Using Vegetation Growing Phenology Feature. In Proceedings of the IGARSS 2023–2023 IEEE International Geoscience and Remote Sensing Symposium, Pasadena, CA, USA, 16–21 July 2023; pp. 3046–3049. [Google Scholar] [CrossRef]
  51. Padhee, S.K.; Dutta, S. Spatio-Temporal Reconstruction of MODIS NDVI by Regional Land Surface Phenology and Harmonic Analysis of Time-Series. GIsci. Remote Sens. 2019, 56, 1261–1288. [Google Scholar] [CrossRef]
  52. Son, N.; Chen, C.; Chen, C.; Guo, H. Classification of Multitemporal Sentinel-2 Data for Field-Level Monitoring of Rice Cropping Practices in Taiwan. Adv. Space Res. 2020, 65, 1910–1921. [Google Scholar] [CrossRef]
  53. Beck, P.S.A.; Jönsson, P.; Høgda, K.A.; Karlsen, S.R. A Ground-Validated NDVI Dataset for Monitoring Vegetation Dynamics and Mapping Phenology in Fennoscandia and the Kola Peninsula. Int. J. Remote Sens. 2007, 28, 4311–4330. [Google Scholar] [CrossRef]
  54. Maignan, F.; Bréon, F.; Bacour, C.; Demarty, J.; Poirson, A. Interannual Vegetation Phenology Estimates from Global AVHRR Measurements Comparison with In Situ Data and Applications. Remote Sens. Environ. 2008, 112, 496–505. [Google Scholar] [CrossRef]
  55. AEMET/IMP. Atlas Climático Ibérico Iberian Climate Atlas; AEMET—Agencia Estatal de Meteorología, IMP—Instituto de Meteorología de Portugal, Ministerio de Medio Ambiente Rural y Marino, Eds.; Closas Orcoyen S.L.: Madrid, Spain, 2011; ISBN 978-84-7837-079-5.
  56. Beck, H.E.; Zimmermann, N.E.; McVicar, T.R.; Vergopolan, N.; Berg, A.; Wood, E.F. Present and Future Köppen-Geiger Climate Classification Maps at 1-Km Resolution. Sci. Data 2018, 5, 180214. [Google Scholar]
  57. FAO. FAO/Unesco Soil Map of the World, Revised Legend, with Corrections and Updates; World Soil Resources Report 60; FAO: Rome, Italy, 1988. [Google Scholar]
  58. MITECO Ministerio Para La Transición Ecológica y El Reto Demográfico. Mapa Forestal de España 50000; 2013. Available online: https://www.miteco.gob.es/es/biodiversidad/servicios/banco-datos-naturaleza/informacion-disponible/mfe50.html (accessed on 20 October 2023).
  59. Ministerio de Agricultura, Pesca y Alimentación. Calendario de Siembra, Recolección y Comercialización 2014–2016; Madrid, 2021. Available online: https://www.mapa.gob.es/es/estadistica/temas/estadisticas-agrarias/agricultura/calendarios-siembras-recoleccion/ (accessed on 20 October 2023).
  60. Drusch, M.; Del Bello, U.; Carlier, S.; Colin, O.; Fernandez, V.; Gascon, F.; Hoersch, B.; Isola, C.; Laberinti, P.; Martimort, P.; et al. Sentinel-2: ESA’s Optical High-Resolution Mission for GMES Operational Services. Remote Sens. Environ. 2012, 120, 25–36. [Google Scholar] [CrossRef]
  61. Zekoll, V.; Main-Knorn, M.; Louis, J.; Frantz, D.; Richter, R.; Pflug, B. Comparison of Masking Algorithms for Sentinel-2 Imagery. Remote Sens. 2021, 13, 137. [Google Scholar] [CrossRef]
  62. Cao, R.; Chen, Y.; Shen, M.; Chen, J.; Zhou, J.; Wang, C.; Yang, W. A Simple Method to Improve the Quality of NDVI Time-Series Data by Integrating Spatiotemporal Information with the Savitzky-Golay Filter. Remote Sens. Environ. 2018, 217, 244–257. [Google Scholar] [CrossRef]
  63. Cai, Y.; Lin, H.; Zhang, M. Mapping Paddy Rice by the Object-Based Random Forest Method Using Time Series Sentinel-1/Sentinel-2 Data. Adv. Space Res. 2019, 64, 2233–2244. [Google Scholar] [CrossRef]
  64. Menenti, M.; Azzali, S.; Verhoef, W.; Swol, R. Van Mapping Agroecological Zones and Time Lag in Vegetation Growth by Means of Fourier Analysis of Time Series of NDVI Images. Adv. Space Res. 1993, 13, 233–237. [Google Scholar] [CrossRef]
  65. Geng, L.; Ma, M.; Wang, X.; Yu, W.; Jia, S.; Wang, H. Comparison of Eight Techniques for Reconstructing Multi-Satellite Sensor Time-Series NDVI Data Sets in the Heihe River Basin, China. Remote Sens. 2014, 6, 2024–2049. [Google Scholar] [CrossRef]
  66. Buys-Ballot, C.H.D. Les Changements Périodiques de Température, Dépendants de La Nature Du Soleil et de La Lune, Mis En Rapport Avec Le Pronostic Du Temps, Déduits d’observations Néerlandaises de 1729 à 1846; Kemink: Utrecht, The Netherlands, 1847. [Google Scholar]
  67. Ljung, G.M.; Box, G.E.P. On a Measure of Lack of Fit in Time Series Models. Biometrika 1978, 65, 297–303. [Google Scholar] [CrossRef]
  68. Huesca, M.; Litago, J.; Merino-de-Miguel, S.; Cicuendez-López-Ocaña, V.; Palacios-Orueta, A. Modeling and Forecasting MODIS-Based Fire Potential Index on a Pixelbasis Using Time Series Models. Int. J. Appl. Earth Obs. Geoinf. 2014, 26, 363–376. [Google Scholar] [CrossRef]
  69. Huesca, M.; Merino-de-Miguel, S.; Eklundh, L.; Litago, J.; Cicuéndez, V.; Rodríguez-Rastrero, M.; Ustin, S.L.; Palacios-Orueta, A. Ecosystem Functional Assessment Based on the “Optical Type” Concept and Self-Similarity Patterns: An Application Using MODIS-NDVI Time Series Autocorrelation. Int. J. Appl. Earth Obs. Geoinf. 2015, 43, 132–148. [Google Scholar] [CrossRef]
  70. Hamilton, J.D. Time Series Analysis; Princeton University Press: Princeton, MA, USA, 1994. [Google Scholar]
  71. Recuero, L.; Litago, J.; Pinzón, J.E.; Huesca, M.; Moyano, M.C.; Palacios-Orueta, A. Mapping Periodic Patterns of Global Vegetation Based on Spectral Analysis of NDVI Time Series. Remote Sens. 2019, 11, 2497. [Google Scholar] [CrossRef]
  72. Fuller, W.A. Introduction to Statistical Time Series; Wiley: New York, NY, USA, 1976. [Google Scholar]
  73. Chen, Y.; Cao, R.; Chen, J.; Liu, L.; Matsushita, B. A Practical Approach to Reconstruct High-Quality Landsat NDVI Time-Series Data by Gap Filling and the Savitzky–Golay Filter. ISPRS J. Photogramm. Remote Sens. 2021, 180, 174–190. [Google Scholar] [CrossRef]
  74. Fontana, F.; Rixen, C.; Jonas, T.; Aberegg, G.; Wunderle, S. Alpine Grassland Phenology as Seen in AVHRR, VEGETATION, and MODIS NDVI Time Series—A Comparison with In Situ Measurements. Sensors 2008, 8, 2833–2853. [Google Scholar] [CrossRef] [PubMed]
  75. Van Leeuwen, W.J.D.; Huete, A.R.; Laing, T.W. MODIS Vegetation Index Compositing Approach: A Prototype with AVHRR Data. Remote Sens. Environ. 1999, 69, 264–280. [Google Scholar] [CrossRef]
  76. Maselli, F.; Battista, P.; Chiesi, M.; Rapi, B.; Angeli, L.; Fibbi, L.; Magno, R.; Gozzini, B. Int J Appl Earth Obs Geoinformation Use of Sentinel-2 MSI Data to Monitor Crop Irrigation in Mediterranean Areas. Int. J. Appl. Earth Obs. Geoinf. 2020, 93, 102216. [Google Scholar] [CrossRef]
  77. Tong, X.; Tian, F.; Brandt, M.; Liu, Y.; Zhang, W.; Fensholt, R. Trends of Land Surface Phenology Derived from Passive Microwave and Optical Remote Sensing Systems and Associated Drivers across the Dry Tropics. Remote Sens. Environ. 2019, 232, 111307. [Google Scholar] [CrossRef]
  78. Liu, R.; Shang, R.; Liu, Y.; Lu, X. Global Evaluation of Gap- Fi Lling Approaches for Seasonal NDVI with Considering Vegetation Growth Trajectory, Protection of Key Point, Noise Resistance and Curve Stability. Remote Sens. Environ. 2017, 189, 164–179. [Google Scholar] [CrossRef]
  79. Eilers, P.H.C.; Pesendorfer, V.; Bonifacio, R. Automatic Smoothing of Remote Sensing Data. In Proceedings of the 9th International Workshop on the Analysis of Multitemporal Remote Sensing Images (MultiTemp), Brugge, Belgium, 27–29 June 2017; pp. 1–3. [Google Scholar]
  80. Vuolo, F.; Atzberger, C. Exploiting the Classification Performance of Support Vector Machines with Multi-Temporal Moderate-Resolution Imaging Spectroradiometer (MODIS) Data in Areas of Agreement and Disagreement of Existing Land Cover Products. Remote Sens. 2012, 4, 3143–3167. [Google Scholar] [CrossRef]
  81. Patel, J.H.; Oza, M.P. Deriving Crop Calendar Using NDVI Time-Series. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2014, 40, 869–873. [Google Scholar] [CrossRef]
  82. Reza, H.; Malamiri, G.; Zare, H.; Rousta, I.; Olafsson, H.; Verdiguier, E.I.; Zhang, H.; Mushore, T.D. Comparison of Harmonic Analysis of Time Series (HANTS) and Multi-Singular Spectrum Analysis (M-SSA) in Reconstruction of Long-Gap Missing Data in NDVI Time Series. Remote Sens. 2020, 12, 2747. [Google Scholar] [CrossRef]
  83. Cetin, M.; Aksoy, T.; Nihan, S.; Anil, M.; Kurkcuoglu, S.; Cabuk, A. Employing Remote Sensing Technique to Monitor the Influence of Newly Established Universities in Creating an Urban Development Process on the Respective Cities. Land Use Policy 2021, 109, 105705. [Google Scholar] [CrossRef]
  84. Mehta, D.H.; Shukla, S.H.; Kalubarme, M.H.; Mehta, D.H. Winter Crop Growth Monitoring Using Multi-Temporal NDVI Profiles in Kapadvanj Taluka, Gujarat State. Int. J. Environ. Geoinform. 2021, 8, 33–38. [Google Scholar] [CrossRef]
  85. Cicuéndez, V.; Litago, J.; Sánchez-Girón, V.; Recuero, L.; Sáenz, C.; Palacios-Orueta, A. Identification and Modeling Carbon and Energy Fluxes from Eddy Covariance Time Series Measurements in Rice and Rainfed Crops. Eng. Proc. 2021, 9, 9. [Google Scholar] [CrossRef]
  86. Debaeke, P.; Aboudrare, A. Adaptation of Crop Management to Water-Limited Environments. Eur. J. Agron. 2004, 21, 433–446. [Google Scholar] [CrossRef]
  87. Cossani, C.M.; Slafer, G.A.; Savin, R. Do Barley and Wheat (Bread and Durum) Differ in Grain Weight Stability through Seasons and Water-Nitrogen Treatments in a Mediterranean Location? Field Crops Res. 2011, 121, 240–247. [Google Scholar] [CrossRef]
  88. Zhu, W.; Pan, Y.; He, H.; Wang, L.; Mou, M.; Liu, J. A Changing-Weight Filter Method for Reconstructing a High-Quality NDVI Time Series to Preserve the Integrity of Vegetation Phenology. IEEE Trans. Geosci. Remote Sens. 2012, 50, 1085–1094. [Google Scholar] [CrossRef]
  89. Recuero, L.; Maila, L.; Cicuéndez, V.; Sáenz, C.; Litago, J.; Tornos, L.; Merino-de-Miguel, S.; Palacios-Orueta, A. Mapping Cropland Intensification in Ecuador through Spectral Analysis of MODIS NDVI Time Series. Agronomy 2023, 13, 2329. [Google Scholar] [CrossRef]
  90. Cicuéndez, V.; Litago, J.; Sánchez-Girón, V.; Román-Cascón, C.; Recuero, L.; Saénz, C.; Yagüe, C.; Palacios-Orueta, A. Dynamic Relationships between Gross Primary Production and Energy Partitioning in Three Different Ecosystems Based on Eddy Covariance Time Series Analysis. Front. For. Glob. Chang. 2023, 6, 1017365. [Google Scholar] [CrossRef]
  91. Aragones, D.; Rodriguez-Galiano, V.F.; Caparros-Santiago, J.A.; Navarro-Cerrillo, R.M. Could Land Surface Phenology Be Used to Discriminate Mediterranean Pine Species? Int. J. Appl. Earth Obs. Geoinf. 2019, 78, 281–294. [Google Scholar] [CrossRef]
  92. Mantas, V.; Fonseca, L.; Baltazar, E.; Canhoto, J.; Abrantes, I. Detection of Tree Decline (Pinus Pinaster Aiton) in European Forests Using Sentinel-2 Data. Remote Sens. 2022, 14, 2028. [Google Scholar] [CrossRef]
Remotesensing 16 02980 g001
Figure 1. Image in true color from © Google Earth (Landsat/Copernicus). (a) Map of Spain and (b) study area.
Figure 1. Image in true color from © Google Earth (Landsat/Copernicus). (a) Map of Spain and (b) study area.
Remotesensing 16 02980 g001
Remotesensing 16 02980 g002
Figure 2. Sentinel-2 image processing flowchart of time series filtering.
Figure 2. Sentinel-2 image processing flowchart of time series filtering.
Remotesensing 16 02980 g002
Remotesensing 16 02980 g003
Figure 3. Spatial distribution of the vegetation species in the study area.
Figure 3. Spatial distribution of the vegetation species in the study area.
Remotesensing 16 02980 g003
Remotesensing 16 02980 g004
Figure 4. Valid observations according to the SCL band of Sentinel-2. (a) Spatial distribution in percentage and (b) number of pixels for each level of valid observations.
Figure 4. Valid observations according to the SCL band of Sentinel-2. (a) Spatial distribution in percentage and (b) number of pixels for each level of valid observations.
Remotesensing 16 02980 g004
Remotesensing 16 02980 g005
Figure 5. (a) Length of the longest gap within time series; (b) season in which the longest gap occurs.
Figure 5. (a) Length of the longest gap within time series; (b) season in which the longest gap occurs.
Remotesensing 16 02980 g005
Remotesensing 16 02980 g006
Figure 6. Interpolating Efficiency Indicator for Tile 30TUN: (a) spatial distribution and (b) frequency distribution.
Figure 6. Interpolating Efficiency Indicator for Tile 30TUN: (a) spatial distribution and (b) frequency distribution.
Remotesensing 16 02980 g006
Remotesensing 16 02980 g007
Figure 7. Time series of pure pixels for a wheat (a), barley (c), maize (e), and alfalfa (g) crop. Time series of the interpolated NDVI (−) and its four filters: SG (−), FFT (−), WHI (−), and MVF (−), and (b,d,f,h) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.
Figure 7. Time series of pure pixels for a wheat (a), barley (c), maize (e), and alfalfa (g) crop. Time series of the interpolated NDVI (−) and its four filters: SG (−), FFT (−), WHI (−), and MVF (−), and (b,d,f,h) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.
Remotesensing 16 02980 g007
Remotesensing 16 02980 g008
Figure 8. Time series of pure pixels for a beech forest (a) and a pine forest (c). Time series of the interpolated NDVI (−) and its four filters: SG (−), FFT (−), WHI (−), and MVF (−), and (b,d) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.
Figure 8. Time series of pure pixels for a beech forest (a) and a pine forest (c). Time series of the interpolated NDVI (−) and its four filters: SG (−), FFT (−), WHI (−), and MVF (−), and (b,d) the average year, i.e., Buys-Ballot. The colors of each filter and interpolation are the same as those of the entire time series.
Remotesensing 16 02980 g008
Remotesensing 16 02980 g009
Figure 9. Representation of the evolution at the image level of the interpolation and filtering processes of the Sentinel-2 NDVI time series using the Whittaker filter. (a) NDVI image (10/06/20) without interpolation, (b) NDVI image (15/06/20) without interpolation, (c) NDVI image (20/06/20) without interpolation, (d) NDVI image (15/06/20) interpolated, and (e) NDVI image (15/06/20) filtered using the Whittaker filter.
Figure 9. Representation of the evolution at the image level of the interpolation and filtering processes of the Sentinel-2 NDVI time series using the Whittaker filter. (a) NDVI image (10/06/20) without interpolation, (b) NDVI image (15/06/20) without interpolation, (c) NDVI image (20/06/20) without interpolation, (d) NDVI image (15/06/20) interpolated, and (e) NDVI image (15/06/20) filtered using the Whittaker filter.
Remotesensing 16 02980 g009
Remotesensing 16 02980 g010
Figure 10. Values of the Q-test for the short term (lags 1, 2, 3, 4, 5, 6, and 7).
Figure 10. Values of the Q-test for the short term (lags 1, 2, 3, 4, 5, 6, and 7).
Remotesensing 16 02980 g010
Remotesensing 16 02980 g011
Figure 11. Values of the Q-test at lag 36 (6 months), lag 73 (one year), and lag 146 (two years).
Figure 11. Values of the Q-test at lag 36 (6 months), lag 73 (one year), and lag 146 (two years).
Remotesensing 16 02980 g011
Remotesensing 16 02980 g012
Figure 12. Average Fk value for period 73 (one year) for each vegetation type (left axis) and percentage of increase after interpolation and filtering (right axis).
Figure 12. Average Fk value for period 73 (one year) for each vegetation type (left axis) and percentage of increase after interpolation and filtering (right axis).
Remotesensing 16 02980 g012
Remotesensing 16 02980 g013
Figure 13. Time series of a pixel declared as alfalfa implemented during the first year.
Figure 13. Time series of a pixel declared as alfalfa implemented during the first year.
Remotesensing 16 02980 g013
Table 1. Growing calendar of different crops in northern Spain.
Table 1. Growing calendar of different crops in northern Spain.
CropMonths
Jan.Feb.Mar.Apr.MayJun.Jul.Aug.Sep.Oct.Nov.Dec.
Wheat
Barley
Maize
Alfalfa Several cuts
Beech forest leaf fall
Pine forest
Light green represents the growing period for the different vegetation types, the dark green represents the several cuts in Alfalfa, and the orange represents the leaf fall in Beech Forest.
Table 2. Weights (Wt) for the different numbers of consecutive invalid observations (gaps).
Table 2. Weights (Wt) for the different numbers of consecutive invalid observations (gaps).
Gaps (x)Weight (Wt)Gaps (x)Weight (Wt)
10.198
20.3109
30.51110
40.71212
50.91314
621416
741518
861620
Table 3. Statistics of the longest gap length values and percentage of the total surface where the longest gap is in each season.
Table 3. Statistics of the longest gap length values and percentage of the total surface where the longest gap is in each season.
SeasonLength of the Longest Gap
MeanStandard DeviationPercentage (%)
Winter9.504.3446.78
Spring7.872.3334.31
Summer4.261.782.14
Autumn7.163.0616.78
Dates of each season: Winter (21 Dec.–19 Mar.), Spring (20 Mar.–19 Jun.), Summer (20 Jun.–21 Sep.), and Autumn (22 Sep.–20 Dec.).
Table 4. IEI values by vegetation species and season.
Table 4. IEI values by vegetation species and season.
Vegetation SpeciesAnnualWinterSpringSummerAutumn
Wheat91.2386.3691.7397.4991.51
Barley91.7887.9192.3097.5891.53
Maize92.8591.4091.5597.7892.23
Alfalfa92.2588.9391.9897.5892.05
Beech forest74.9056.1081.8992.7472.57
Pine forest88.1584.3487.1696.4888.05
Crops (Cs), Beech Forest (Cf), and Pine Forest (Transition).
Table 5. Average values of Correlation Coefficient and the Root Mean Square Error for the plots of each filter and vegetation species.
Table 5. Average values of Correlation Coefficient and the Root Mean Square Error for the plots of each filter and vegetation species.
CropsForest Species
WheatBarleyMaizeAlfalfaBeechPine
FiltersCCRMSECCRMSECCRMSECCRMSECCRMSECCRMSE
SG0.970.040.980.040.980.050.950.070.940.080.890.08
WHI0.940.060.950.060.960.070.870.100.850.130.760.11
FFT0.930.070.940.060.950.080.870.100.840.130.730.12
MVF0.910.100.920.100.930.120.830.150.800.170.660.15
SG: Savitzky–Golay (SG), Whittaker (WHI), Fast Fourier Transform (FFT), and Maximum Value Filter (MVF).
Table 6. Mean Q-test for time series of vegetation species plots (raw).
Table 6. Mean Q-test for time series of vegetation species plots (raw).
Vegetation SpeciesQ1Q2Q3Q4Q5Q6Q7Q36Q73Q146
Wheat65.71132.21171.03212.49243.69269.16291.74504.05718.86954.10
Barley70.50140.43178.93221.23251.81275.36295.50500.39787.051158.08
Maize155.56302.68428.80539.56629.75705.33763.331377.242650.964611.12
Alfalfa26.2248.8158.0368.6177.5684.2092.33174.69266.44424.26
Beech forest24.0647.1665.8080.9296.38109.45122.64342.62677.621227.30
Pine forest12.3822.2127.2630.5134.2238.5642.11100.06180.82336.58
The critical values of a χ2 (Chi-Square distribution) at 5% for the Q-Ljung-Box test with 1, 2, 3, 4, 5, 6, 7, 36, 73 and 146 degrees of freedom are 3.84, 5.99, 7.81, 9.49, 11.07, 12.59, 14.07, 51.00, 93.95, and 175.20, respectively.
Table 7. Percentage of surface in which the maximum ordinate of the NDVI periodogram corresponds to periods 36, 73, 203, and others for different filters and land covers.
Table 7. Percentage of surface in which the maximum ordinate of the NDVI periodogram corresponds to periods 36, 73, 203, and others for different filters and land covers.
Land CoversPeriodFraction of Surface (%)
RawInterpolatedSGWHIFFTMVF
Wheat360.251.331.331.041.331.13
7387.3184.6584.6084.9684.6083.70
2030.280.780.780.820.781.24
Others12.1613.2413.2913.1813.2913.93
Barley361.124.214.203.384.183.19
7394.9788.7688.7289.5088.7488.54
2030.080.310.320.330.320.46
Others3.836.726.766.796.767.81
Maize362.194.914.904.134.921.19
7390.9886.2886.2487.0786.2390.08
2030.110.780.780.800.770.96
Others6.728.038.088.008.087.77
Alfalfa361.178.378.387.298.425.92
7374.0049.4149.4350.3449.4550.34
20313.1233.5633.5634.2033.5237.11
Others11.718.668.638.178.616.63
Beech forest360.000.170.180.160.180.18
7391.1998.7698.7198.8398.7198.78
2030.010.240.240.260.240.36
Others8.800.830.870.750.870.68
Pine forest365.458.238.318.158.513.92
7373.2048.0248.3052.8548.3169.90
2030.8615.8215.8517.3215.5718.04
Others20.4927.9327.5421.6827.618.14
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

Sáenz, C.; Cicuéndez, V.; García, G.; Madruga, D.; Recuero, L.; Bermejo-Saiz, A.; Litago, J.; de la Calle, I.; Palacios-Orueta, A. New Insights on the Information Content of the Normalized Difference Vegetation Index Sentinel-2 Time Series for Assessing Vegetation Dynamics. Remote Sens. 2024, 16, 2980. https://doi.org/10.3390/rs16162980

AMA Style

Sáenz C, Cicuéndez V, García G, Madruga D, Recuero L, Bermejo-Saiz A, Litago J, de la Calle I, Palacios-Orueta A. New Insights on the Information Content of the Normalized Difference Vegetation Index Sentinel-2 Time Series for Assessing Vegetation Dynamics. Remote Sensing. 2024; 16(16):2980. https://doi.org/10.3390/rs16162980

Chicago/Turabian Style

Sáenz, César, Víctor Cicuéndez, Gabriel García, Diego Madruga, Laura Recuero, Alfonso Bermejo-Saiz, Javier Litago, Ignacio de la Calle, and Alicia Palacios-Orueta. 2024. "New Insights on the Information Content of the Normalized Difference Vegetation Index Sentinel-2 Time Series for Assessing Vegetation Dynamics" Remote Sensing 16, no. 16: 2980. https://doi.org/10.3390/rs16162980

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