Suitability of the MODIS-NDVI Time-Series for a Posteriori Evaluation of the Citrus Tristeza Virus Epidemic

The technological advances of remote sensing (RS) have allowed its use in a number of fields of application including plant disease depiction. In this study, an RS approach based on an 18-year (i.e., 2001–2018) time-series analysis of Normalized Difference Vegetation Index (NDVI) data, derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) and processed with TIMESAT free software, was applied in Sicily (insular Italy). The RS approach was carried out in four orchards infected by Citrus tristeza virus (CTV) at different temporal stages and characterized by heterogeneous conditions (e.g., elevation, location, plant age). The temporal analysis allowed the identification of specific metrics of the NDVI time-series at the selected sites during the study period. The most reliable parameter which was able to identify the temporal evolution of CTV syndrome and the impact of operational management practices was the “Base value” (i.e., average NDVI during the growing seasons, which reached R2 values up to 0.88), showing good relationships with “Peak value”, “Small integrated value” and “Amplitude”, with R2 values of 0.63, 0.70 and 0.75, respectively. The approach herein developed is valid to be transferred to regional agencies involved in and/or in charge of the management of plant diseases, especially if it is integrated with ground-based early detection methods or high-resolution RS approaches, in the case of quarantine plant pathogens requiring control measures at large-scale level.


Introduction
The adoption of remote sensing (RS) approaches is gaining special relevance to monitor, quantify and map vegetation dynamics resulting from life-cycle patterns, climatic conditions, photosynthetic activity and plant diseases [1,2]. Plant pathogens cause damage to crops, quantifiable in direct and indirect costs. In addition to the direct costs, quantified as yield losses, in the case of virus diseases and consequent quarantine plan, the indirect costs (i.e., plant protection treatments, environmental impact of pesticides, replacement of plants and loss of biodiversity) find the major item being prevention costs with a strong economic and social impact for the community. The most modern challenges of plant pathology are to combine the classic systems of prevention and quarantine with modern automated methods of detection, which, being characterized by precision and reliability, allow for a more efficient management of phytosanitary emergencies. In situ monitoring approaches for plant health have made tremendous progress, but they are intensive and often integrate subjective indicators. RS bridges the fragmentation of citriculture in this region and late legislative action. Effective management of the disease would have required the assessment of epidemic evolution followed by the application of medium and long-term containment programs on a territorial scale, which also included the eradication and replanting of tolerant rootstocks.
In this scenario, the approach carried out in this study is based on the application of the asymmetric function CF method on NDVI time-series products derived by MODIS for depicting the Citrus tristeza virus (CTV) epidemic effects. In light of the aforementioned scientific context, the aim of the study is to identify the dynamics of the vegetation of citrus trees in relation to the occurrence and/or evolution of the signs of tristeza disease in the agricultural context of Eastern Sicily (southern Italy). The study aims to develop an effective, rapid and exportable tool for the monitoring and management of cultivated areas affected by plant pathogens similar to CTV.

Case Studies Selection Criteria
Four case studies (CS), named from CS1 to CS4, have been selected in Eastern Sicily (insular Italy) for the application of the proposed RS methodological approach (Figure 1). The study area (Catania Plain, located at the foothills of Etna Volcano) is the most suitable production area for blood orange (pigmented cultivars) and is characterized by semi-arid Mediterranean climatic conditions. Remote Sens. 2020, 12, x 3 of 21 management of the disease would have required the assessment of epidemic evolution followed by the application of medium and long-term containment programs on a territorial scale, which also included the eradication and replanting of tolerant rootstocks. In this scenario, the approach carried out in this study is based on the application of the asymmetric function CF method on NDVI time-series products derived by MODIS for depicting the Citrus tristeza virus (CTV) epidemic effects. In light of the aforementioned scientific context, the aim of the study is to identify the dynamics of the vegetation of citrus trees in relation to the occurrence and/or evolution of the signs of tristeza disease in the agricultural context of Eastern Sicily (southern Italy). The study aims to develop an effective, rapid and exportable tool for the monitoring and management of cultivated areas affected by plant pathogens similar to CTV.

Case Studies Selection Criteria
Four case studies (CS), named from CS1 to CS4, have been selected in Eastern Sicily (insular Italy) for the application of the proposed RS methodological approach (Figure 1). The study area (Catania Plain, located at the foothills of Etna Volcano) is the most suitable production area for blood orange (pigmented cultivars) and is characterized by semi-arid Mediterranean climatic conditions. The CSs were selected on the basis of the following common features (Table 1): (i) presence of the rootstock susceptible to CTV (Citrus aurantium L.); (ii) an extent identifiable with the spatial resolution of the MODIS pixel (250 × 250 m); (iii) characterized by a quite uniform canopy ground cover (i.e., higher than 80%) within the entire MODIS pixel dimension (i.e., between rows, there is only the natural ground cover).
All the CSs were supplied by full (i.e., potential) irrigation, in order to avoid water deficit and potential damage to crops due to water stress conditions e.g., [45,46]. In particular, CS1 and CS2 were supplied by micro-irrigation systems, while CS3 and CS4 were subjected to drip irrigation. The soil type at the CSs is clay loam (USDA classification) with mean % (± standard deviation) of 42.1 (± 4.8), 28.0 (± 4.6), 29.9 (± 3.9) for silt, sand and clay, respectively, and bulk density varies from 0.92 to 1.19 gr cm −3 [47,48]. The CSs were selected on the basis of the following common features (Table 1): (i) presence of the rootstock susceptible to CTV (Citrus aurantium L.); (ii) an extent identifiable with the spatial resolution of the MODIS pixel (250 × 250 m); (iii) characterized by a quite uniform canopy ground cover (i.e., higher than 80%) within the entire MODIS pixel dimension (i.e., between rows, there is only the natural ground cover).
All the CSs were supplied by full (i.e., potential) irrigation, in order to avoid water deficit and potential damage to crops due to water stress conditions e.g., [45,46]. In particular, CS1 and CS2 were supplied by micro-irrigation systems, while CS3 and CS4 were subjected to drip irrigation. The soil type at the CSs is clay loam (USDA classification) with mean % (± standard deviation) of 42.1 (± 4.8), 28.0 (± 4.6), 29.9 (± 3.9) for silt, sand and clay, respectively, and bulk density varies from 0.92 to 1.19 gr cm −3 [47,48]. The preliminary evaluation of the evolution of the CTV and the timing of the related corrective actions (e.g., decline of the green, eradication of plants) during the MODIS 2001-2018 monitoring is based on the photointerpretation of the satellite images of Google Earth (https://earth.google.com/web/) and on interviews with farmers. To that purpose, the spatial resolution of the Google images was scaled at the MODIS pixel resolution ( Figure 1).

MODIS Data
Free RS data were retrieved from MODIS sensor on board of the Terra satellite (https://modis. gsfc.nasa.gov/). Specifically, Terra-MODIS vegetation index (VI) product (MOD13Q1, Version 6 Level 3, [49]), generated every 16 days at 250 m of spatial resolution, was used for the time-series analyses. This product provides two primary vegetation layers: NDVI and Enhanced Vegetation Index (EVI). The MODIS algorithm chooses the best available pixel value from all the acquisitions from the 16-day period, on the basis of low cloud coverage, low view angle and the highest NDVI/EVI value. The main difference between the two VI layers is that EVI has improved sensitivity over high biomass areas, avoiding the problem of NDVI saturation [50]. However, in this study, NDVI was selected as the representative index of the vegetation dynamics because of the absence of NDVI saturation at the selected CSs (i.e., orange orchards).
In general, NDVI is representative of vegetation vigor [33], offering a valid tool for identifying/ monitoring the vegetation status, mainly related to biotic and abiotic site conditions. It is calculated as follows: where NIR and RED are the reflectance values in the red and near-infrared bands of the electromagnetic spectrum, respectively. In this study, NDVI time-series, referring to CS1-CS4, were analyzed in the period from 1 January 2001 to 31 December 2018. Data of the vegetation index were extracted and downloaded from the Global Subsets Tool [51]. This tool provides a summary of selected MODIS and VIIRS (Visible Infrared Imaging Radiometer Suite) products that can be freely used by the community for several applications (e.g., for validating models and remote-sensing products and for characterizing field sites) (https://modis.ornl.gov/cgi-bin/MODIS/global/subset.pl). The pixel selection (ID pixel, Table 1 and Figure 1) was defined using the geographic coordinates (latitude and longitude) of the center of the four selected CS. Eighteen years of NDVI time-series (23 images per year) were retrieved by Terra-MODIS for a total of 414 layers for each CS. Data quality was assessed by checking the VI Quality and Pixel Reliability indicators, that come together with the MOD13Q1 product.

Meteorological Data Clustering
Daily meteorological data, including solar radiation (R s , MJ m −2 ), relative humidity (RH, minimum and maximum; %), air temperature (T air , minimum, maximum, mean; • C), wind speed and direction (u 2 , Remote Sens. 2020, 12,1965 5 of 20 m s −1 and WD, • ) and rainfall (P, mm), were obtained by the Sistema Informativo Agrometeorologico Siciliano (SIAS). Specifically, the meteorological data were measured at two SIAS stations (in operation since 2002) located 6 km from CS1 (S234, 37.51 • N, 14.85 • E, 100 m a.s.l.) and about 2-8 km from CS2-CS4 (S292, 37.35 • N, 14.91 • E, 50 m a.s.l.). Daily meteorological data were averaged at 16-day intervals as the Terra-MODIS products used for deriving the NDVI time-series (see Section 2.2). Figure 2 reports the 16-day meteorological data recorded at the weather stations S234 and S292 during the 2002-2018 period. No significant changes were observed in the main meteorological patterns during the reference period. In particular, the daily average values (± standard deviation) at 16-day intervals were 17.12 ± 6.56 MJ m −2 for R s , 63.32 ± 9.63% for RH, 18.18 ± 6.25 • C T air , 1.46 ± 0.38 m s −1 for u 2 Figure 2 reports the 16-day meteorological data recorded at the weather stations S234 and S292 during the 2002-2018 period. No significant changes were observed in the main meteorological patterns during the reference period. In particular, the daily average values (± standard deviation) at 16-day intervals were 17.12 ± 6.56 MJ m −2 for Rs, 63.32 ± 9.63% for RH, 18.18 ± 6.25 °C Tair, 1.46 ± 0.38 m s −1 for u2 and 1.54 ± 2.21 mm for P.

TIMESAT Curving Fitting Method
The TIMESAT software (version 3.3), developed by [41], was used for generating smoothed NDVI time-series in the period 2001-2018 from Terra-MODIS satellite spectral measurements (see Section 2.2). The TIMESAT software is freely available at http://web.nateko.lu.se/timesat/timesat.asp and has a user-friendly graphical interface. It implements three CF methods: (i) the adaptive Savitzky-Golay filter that uses local polynomial functions for data fitting; and (ii) the asymmetric Gaussian and (iii) the double logistic model functions, both based on the least-squares methods, where data are fitted to non-linear model functions of different complexity. For these latter methods, the model functions are fitted to the data in intervals between maxima and minima of the time-series (t). The general form of the model functions is defined as: where linear c parameters determine the base level (c1) and the amplitude (c2) for the seasons. The non-linear parameter x determines the shape of the basis function g (t; x) that in this study was considered as a double logistic filter [52], as reported in Equation (3):

TIMESAT Curving Fitting Method
The TIMESAT software (version 3.3), developed by [41], was used for generating smoothed NDVI time-series in the period 2001-2018 from Terra-MODIS satellite spectral measurements (see Section 2.2). The TIMESAT software is freely available at http://web.nateko.lu.se/timesat/timesat.asp and has a user-friendly graphical interface. It implements three CF methods: (i) the adaptive Savitzky-Golay filter that uses local polynomial functions for data fitting; and (ii) the asymmetric Gaussian and (iii) the double logistic model functions, both based on the least-squares methods, where data are fitted to non-linear model functions of different complexity. For these latter methods, the model functions are fitted to the data in intervals between maxima and minima of the time-series (t). The general form of the model functions is defined as: Remote Sens. 2020, 12,1965 6 of 20 where linear c parameters determine the base level (c 1 ) and the amplitude (c 2 ) for the seasons. The non-linear parameter x determines the shape of the basis function g (t; x) that in this study was considered as a double logistic filter [52], as reported in Equation (3): where the non-linear parameters x 1 and x 3 determine the position of the left and right inflection points for the season, respectively, whereas x 2 and x 4 determine the time period of increase and decrease (i.e., rate of change), respectively. Firstly, a pre-processing phase was performed. It consists of the seasonality definition (e.g., fixed unimodal), which considers the growing seasons' timings for the CSs under study (i.e., annual for citrus). No spikes and outliers were removed from the Terra-MODIS NDVI time-series at the CSs due to the high quality of the RS products used. Then, the season "start" and "end" were determined using the seasonal amplitude method, defined as the difference between the base level and the maximum NDVI value for each individual season. As reported in [53], the "start" occurs when the left part of the fitted curve has reached a specified fraction of the amplitude, counted from the base level. The "end" of the season is defined similarly, but for the right side of the curve. The details on the setting parameters used in this study are reported in Table 2. In order to exploit the information of the NDVI time-series, the derived phenological TIMESAT metrics were extracted and analyzed at the selected CSs during the reference period. Figure 3 shows the typical representation of the seasonality parameters [41] (as reported in Table 3).
Remote Sens. 2020, 12, x 6 of 21 where the non-linear parameters x1 and x3 determine the position of the left and right inflection points for the season, respectively, whereas x2 and x4 determine the time period of increase and decrease (i.e., rate of change), respectively. Firstly, a pre-processing phase was performed. It consists of the seasonality definition (e.g., fixed unimodal), which considers the growing seasons' timings for the CSs under study (i.e., annual for citrus). No spikes and outliers were removed from the Terra-MODIS NDVI time-series at the CSs due to the high quality of the RS products used. Then, the season "start" and "end" were determined using the seasonal amplitude method, defined as the difference between the base level and the maximum NDVI value for each individual season. As reported in [53], the "start" occurs when the left part of the fitted curve has reached a specified fraction of the amplitude, counted from the base level. The "end" of the season is defined similarly, but for the right side of the curve. The details on the setting parameters used in this study are reported in Table 2. In order to exploit the information of the NDVI time-series, the derived phenological TIMESAT metrics were extracted and analyzed at the selected CSs during the reference period. Figure 3 shows the typical representation of the seasonality parameters [41] (as reported in Table 3).   (h) small integrated value, (h + i) large integrated value. Graph provided by [41]. Table 3. Description of the seasonal parameters analyzed in the study.

Seasonal Parameters Description
Length of the season Time from the "start" to the "end" of the season Peak value Maximum NDVI for the fitted function during the season Base level Average of the left and the right minimum NDVI values Seasonal amplitude Difference between maximum NDVI and the base level Small seasonal integral Small integrated NDVI value for the fitted function during the season Large seasonal integral Integral of the function describing the season from the "start" to the "end"

Statistical Analysis
Linear regression models (Statistix analytical software, v. 9.0, Tallahassee, FL, USA) were adopted for analyzing the seasonal parameter trends. These trends were identified on the basis of the increasing/decreasing patterns (regression slope) and goodness-of-fit (coefficient of determination, R 2 ). The significance of the trends was assessed for each seasonal parameter (Table 3) using the t-test at significance levels (p-value) of 0.05, 0.01 and 0.001, respectively.

Case Studies Selection
The preliminary analysis of the Google Earth imagery provides useful insight for visualizing the CTV syndrome effects at the CSs ( Figure 4). However, the limited temporal availability of the images does not permit a clear picture of the on-site CTV-related dynamics. In particular, CS1 (located 4 km in the southern-eastern direction from the first CTV epidemic focus) showed the presence of an irregular distribution of CTV symptomatic trees since 2002 ( Figure 1). CS2-CS4 showed quite similar conditions, with the distinct identification of the following phases, corresponding to the epidemic evolution of CTV as a function of the distance from the first area of focus of the virus (Figures 1 and 4): (i) growing phase (e.g., more evident in CS4 since trees were younger than in CS2 and CS3) together with the appearance of the symptomatic trees, and (ii) the adoption of corrective operations for the containment of CTV. In particular, CS2 (located 16 km in the south-western direction from the first CTV epidemic focus) exhibited the presence of the first symptomatic trees since 2007. At CS3 and CS4, the beginning of the CTV plant decline was observed starting from 2010 and 2013, respectively. These CSs were located 18 and 23 km from the CTV epidemic focus in the southern and southern-eastern directions, respectively.

Terra-MODIS NDVI Data
The analysis of the MOD13Q1 VI pixel quality and reliability indicators showed that the overall quality of the input data used here is classified as "good data" for 99% of the NDVI products. Figure 5 shows the temporal patterns of the Terra-MODIS NDVI data referring to CS1-CS4 during the 2001-2018 period. NDVI values ranged from 0.28 to 0.81 in CS1, from 0.43 to 0.89 in CS2, from 0.47 to 0.90 Remote Sens. 2020, 12,1965 9 of 20 in CS3 and from 0.33 to 0.86 in CS4 (Figure 4). The highest mean (± standard deviation) NDVI values were observed at CS2 and CS3 (0.74 ± 0.09 and 0.76 ± 0.08, respectively), whereas lower NDVI values were detected at CS1 (0.59 ± 0.14) and CS4 (0.64 ± 0.11). Figure 6 shows the monthly NDVI box diagram for CS1-CS4 sites in the period of 2001-2018; this graph shows the minimum and maximum NDVI values found, respectively, in Summer (June-August) and in Autumn-Winter (January-February/October-December).

Terra-MODIS NDVI Data
The analysis of the MOD13Q1 VI pixel quality and reliability indicators showed that the overall quality of the input data used here is classified as "good data" for 99% of the NDVI products. Figure  5 shows the temporal patterns of the Terra-MODIS NDVI data referring to CS1-CS4 during the 2001-2018 period. NDVI values ranged from 0.28 to 0.81 in CS1, from 0.43 to 0.89 in CS2, from 0.47 to 0.90 in CS3 and from 0.33 to 0.86 in CS4 (Figure 4). The highest mean (± standard deviation) NDVI values were observed at CS2 and CS3 (0.74 ± 0.09 and 0.76 ± 0.08, respectively), whereas lower NDVI values were detected at CS1 (0.59 ± 0.14) and CS4 (0.64 ± 0.11). Figure 6 shows the monthly NDVI box diagram for CS1-CS4 sites in the period of 2001-2018; this graph shows the minimum and maximum NDVI values found, respectively, in Summer (June-August) and in Autumn-Winter (January-February/October-December).

Terra-MODIS NDVI Data
The analysis of the MOD13Q1 VI pixel quality and reliability indicators showed that the overall quality of the input data used here is classified as "good data" for 99% of the NDVI products. Figure  5 shows the temporal patterns of the Terra-MODIS NDVI data referring to CS1-CS4 during the 2001-2018 period. NDVI values ranged from 0.28 to 0.81 in CS1, from 0.43 to 0.89 in CS2, from 0.47 to 0.90 in CS3 and from 0.33 to 0.86 in CS4 (Figure 4). The highest mean (± standard deviation) NDVI values were observed at CS2 and CS3 (0.74 ± 0.09 and 0.76 ± 0.08, respectively), whereas lower NDVI values were detected at CS1 (0.59 ± 0.14) and CS4 (0.64 ± 0.11). Figure 6 shows the monthly NDVI box diagram for CS1-CS4 sites in the period of 2001-2018; this graph shows the minimum and maximum NDVI values found, respectively, in Summer (June-August) and in Autumn-Winter (January-February/October-December).

TIMESAT Fitting Curves
The time-series analysis was applied successfully to Terra-MODIS NDVI data of the different CSs. The fitted functions obtained by applying the double logistic filter CF method in TIMESAT (in red) to the original NDVI data (in blue) are shown in Figure 7. The "start" and the "end" of the growing seasons are marked with the red filled circles. Generally, the fitting curves reproduced the NDVI data quite well (Figure 7) due to their high quality, typical of clear sky conditions (i.e., absence of spikes of outliers).
The time-series analysis was applied successfully to Terra-MODIS NDVI data of the different CSs. The fitted functions obtained by applying the double logistic filter CF method in TIMESAT (in red) to the original NDVI data (in blue) are shown in Figure 7. The "start" and the "end" of the growing seasons are marked with the red filled circles. Generally, the fitting curves reproduced the NDVI data quite well (Figure 7) due to their high quality, typical of clear sky conditions (i.e., absence of spikes of outliers).

Seasonal Parameters from NDVI Time-Series
The TIMESAT seasonal parameters extracted in CS1-CS4 were useful for depicting the site-specific dynamics that occurred during the reference period (2001-2018) (Figures 8-11).

Seasonal Parameters from NDVI Time-Series
The TIMESAT seasonal parameters extracted in CS1-CS4 were useful for depicting the site-specific dynamics that occurred during the reference period (2001-2018) (Figures 8-11).        In particular, in CS1, near-stable trends of the seasonal parameters were observed i 2001-2018, with slopes ranging from −0.03 to 0.05 and low values of R 2 for all parameters (R 2 between 0 and 0.35, even with the significant trend for "Small integrated value", p < 0.05) except for the "Base value", which reached an R 2 of 0.62 (p < 0.001) (Figure 8).
In contrast to CS1, two different phases were identified in CS2-CS4 on the basis of the photo-interpretation analysis (Figures 9-11). Specifically, in CS2, the parameter that better identified the first phase was the "Base value" (R 2 of 0.48, p < 0.01), whereas in CS3 and CS4, this phase was well represented by the "Peak value" (R 2 of 0.37 and 0.63, respectively, p < 0.05). At CS3, "Amplitude" and "Small integrated value" also contributed in representing the first phase (R 2 of 0.33 and 0.36, respectively, p < 0.05).
Regarding the second phase, the "Base value" (R 2 of 0.58) and the "Peak value" (R 2 of 0.54) represented better this trend in CS2, even with no significance trends (p < 0.05). In CS3, this second phase was more clearly observed than in CS2, with the "Base value" and "Peak value" (R 2 of 0.88 and 0.93, respectively) being the most representative parameters (p < 0.05 and 0.01). The performance obtained from TIMESAT in CS4 was different, with the parameters explaining the trend being the "Length" (R 2 of 0.74), "Base value" (R 2 of 0.81), "Small and Large seasonal integral values" (R 2 of 0.86 and 0.83, respectively); those parameters provided the highest R 2 values, even with no significant trends.
When analyzing individually each seasonal parameter, it was observed that, in general, the trends observed in all the CSs were quite similar (Figures 8-11). The "Length" parameter presented a positive slope term that ranged from 0.10 to 0.19 and from 0.04 to 0.89 in the first and second phases, respectively. The "Base value" exhibited a positive slope term in the first phase (ranging from 0.00 to 0.01 in all CS), whereas the patterns are inverted in the second phase, showing negative slope terms (−0.03 to −0.02). "Peak value" showed a trend similar to "Base value", with positive slopes in the first phase (0.00-0.01) and negative slopes in the second phase (−0.01 to 0.00). Regarding the "Amplitude" and "Small integrated value", in CS3 and CS4 they exhibited positive slope values in the first phase (ranging from 0.00 to 0.01 and from 0.05 to 0.10, for "Amplitude" and "Small integrated value", respectively).
Conversely, in the first phase in CS2, these parameters had negatives slopes (ranging from −0.02 to −0.00). Nevertheless, the behavior of "Amplitude" and "Small integrated value" in the second phase was similar in CS2-CS4, with slopes ranging from 0.01-0.02 and from 0.12-0.44 for "Amplitude" and "Small integrated value", respectively. The "Large integrated value" presented positive slopes in the first phase (0.15-0.23), whereas it changed the sign in the second phase (slope terms of −0.14 and −0.23 at CS2 and CS3, respectively), except for CS4, where a positive slope value was observed (0.50; Figure 11).
By analyzing the inter-seasonality parameter relationship in terms of R 2 , it can be observed that, in general, good relationships were observed between "Length" and "Large integrated value" (R 2 = 0.52) and between "Amplitude" and "Small integrated value" (R 2 = 0.83) (Figure 12). In addition, the "Base value" showed good relationships with "Peak value", "Small integrated value" and "Amplitude", with R 2 values of 0.63, 0.70 and 0.75, respectively. The relationships between the other parameters were quite poor with R 2 values lower than 0.41. positive slope term that ranged from 0.10 to 0.19 and from 0.04 to 0.89 in the first and second phases, respectively. The "Base value" exhibited a positive slope term in the first phase (ranging from 0.00 to 0.01 in all CS), whereas the patterns are inverted in the second phase, showing negative slope terms (−0.03 to −0.02). "Peak value" showed a trend similar to "Base value", with positive slopes in the first phase (0.00-0.01) and negative slopes in the second phase (−0.01 to 0.00). Regarding the "Amplitude" and "Small integrated value", in CS3 and CS4 they exhibited positive slope values in the first phase (ranging from 0.00 to 0.01 and from 0.05 to 0.10, for "Amplitude" and "Small integrated value", respectively). Conversely, in the first phase in CS2, these parameters had negatives slopes (ranging from −0.02 to −0.00). Nevertheless, the behavior of "Amplitude" and "Small integrated value" in the second phase was similar in CS2-CS4, with slopes ranging from 0.01-0.02 and from 0.12-0.44 for "Amplitude" and "Small integrated value", respectively. The "Large integrated value" presented positive slopes in the first phase (0.15-0.23), whereas it changed the sign in the second phase (slope terms of −0.14 and −0.23 at CS2 and CS3, respectively), except for CS4, where a positive slope value was observed (0.50; Figure 11).
By analyzing the inter-seasonality parameter relationship in terms of R 2 , it can be observed that, in general, good relationships were observed between "Length" and "Large integrated value" (R 2 = 0.52) and between "Amplitude" and "Small integrated value" (R 2 = 0.83) (Figure 12). In addition, the "Base value" showed good relationships with "Peak value", "Small integrated value" and "Amplitude", with R 2 values of 0.63, 0.70 and 0.75, respectively. The relationships between the other parameters were quite poor with R 2 values lower than 0.41.

Discussion
TIMESAT software and MODIS products have been widely used in the agriculture context. Specifically, [54] explored their potential for accurately deriving crop phenological metrics in comparison with ground-observed vegetation green-up dates. Additionally, the suitability of TIMESAT with different MODIS-derived VI time series has been assessed for large area vegetation dynamic monitoring over arid and semi-arid lands [55]. TIMESAT has been also used to evaluate the scaling effects on spring phenology detections using MODIS data at multiple spatial resolutions [56]. In addition, [57][58][59] highlighted the strengths of the TIMESAT approach for mapping crop phenological stages and crop-calendar events, and for crop classification. Despite these promising applications, the main weakness of TIMESAT implementation is related to the correct choice of the CF method capable of providing the most robust description of the seasonal dynamics [58][59][60]. Other limitations are due to the influence of noise level of data input on CF technique performance [60,61] and the presence of sample impurity and landscape heterogeneities that can largely affect the classification accuracy [59]. Nowadays, although several authors have compared different CF methods for identifying vegetation features changes through VI metrics, no consensus has been reached in the election of the best CF method (e.g., [22,24,56,[60][61][62]). In general, authors corroborated that the choice of the most suitable CF method (and filtering parameters) depends on the quality of the input data and has to be assessed case-by-case by inspecting how well the CF functions match the row data for preserving the signal integrity [63].
In this study, the double-filter CF method was selected for analyzing the main seasonal parameters of areas affected by CTV in Sicily, with this choice being supported by the excellent quality of the Terra-MODIS NDVI data used as inputs (Figures 5-7). Specifically, no spikes or outliers were detected in the NDVI time-series, which were quite robust, consisting of 18 years of continuous data with 16-day temporal resolution. Even though multiple definitions of seasonality parameters and ways for extracting these parameters are reported in RS literature [53], the importance of the seasonal parameters lies in the possibility to describe the temporal changes in the vegetation cover resulting from abiotic or biotic changes (climatic or land use or disease).
The time-series of NDVI at the different CSs ( Figures 5 and 6) have shown a typical annual pattern deriving from the sum of two contributions: (i) the fairly stable vegetation dynamic of citrus trees and (ii) the presence of typical conditions of grassed soil [64][65][66][67]. It can be assumed, given that the ground cover plays a role of constant "background" every year, that the dynamics of the seasonal NDVI parameters are mainly linked to the presence of citrus groves. In addition, some of the seasonal parameters analyzed (e.g., "Base value", defined as the average of the minimum left and right NDVI values within a season) are independent of the presence of ground cover since they refer to summer periods, when the soil is generally bare due to weed control. Furthermore, the age of trees can play a crucial role in specific seasonal parameter trends. This effect is more evident when the trees are not subject to any management of the size of the canopy (for example, some rain-fed olive groves with large planting structures). However, the size of the crown of citrus groves is commonly managed by pruning, especially in the Sicilian context, where the citrus systems are maintained in a fairly stable size. Therefore, despite the wide age range of the citrus groves examined in this study (which vary from 25 to 60 years, Table 1), and given that the canopy size is kept almost constant, the main differences in the seasonal parameters of NDVI depend on the citrus groves pathology. This is also supported by similar time patterns for the main meteorological variables observed in the two stations located close to the CSs during the reference period (Figures 1 and 2). This allows us to imagine that the observed trends for NDVI are quite independent of the variability of meteorological parameters.
Results derived from TIMESAT analysis reflected a stable pattern of the seasonal parameters in CS1, which can be related to the in situ conditions influenced both by the CTV effects on plants and by the adoption of spotted mitigation actions such as eradication of infected plants and replanting with trees grafted on tolerant rootstocks in order to prevent economic damages since 2002 ( Figure 4). This stability also reflected the low rate of plant growing in CS1 as a consequence of the tree age ( Table 1), indicating that the plant growing stage has been reduced. Nevertheless, the temporal trends of the seasonal parameters in CS2-CS4 allowed recognizing two different phases (Figures 9-11), accordingly identified by the Google Earth imagery analysis (Figure 4). The first phase (observed in the period 2001-2014 in CS2 and CS3 and in the period 2001-2015 in CS4) corresponded to the citrus growing period in combination with the occurrence of the CTV epidemic effects. Therefore, this phase was characterized by a slight increase of all seasonality parameters with the exception of "Amplitude" and "Small integrated value" in CS2, which experienced a decrease, probably because in this site CTV has been present for longer without being countered by corrective measures. The second phase (observed in the period of 2014-2018 in C2 and CS3 and in the period of 2015-2018 in CS4) was characterized by the corrective measures applied by the growers for facing CTV (such as eradication and trees replanting on CTV-tolerant rootstocks) (see the grey areas in Figures 9-11). In this phase "Length", "Amplitude" and "Small integrated value" experienced increasing patterns in CS2-CS4, whereas "Base value", "Peak value" and "Large integrated value" showed a decreasing trend with the exception of "Large integrated value" in CS4. At this site, the intensive corrective actions taken since 2015 by the grower strongly modified the NDVI signals, with an abrupt increase of the "Length" and "Amplitude" (positive slopes), resulting in a more evident increase of "Small integrated value" and "Large integrated value".
In the first phase, a number of parameters were able to accurately describe the vegetation dynamics at the different CS. Specifically, in CS1 (Figures 4 and 8), the contemporary occurrence of CTV with mitigations actions was observed in the "Base value", whereas in CS2-CS4, site-specific patterns in the seasonality parameters were distinguished within the first phase due to tristeza severity. These patterns are more related to the range of variations (minimum and maximum values) of the seasonal parameters (e.g., "Base value", "Peak value", "Amplitude" and "Small integrated value"), than to their "Length".
For example, the CTV epidemic effects were better identified once more by "Base value" in CS2 (i.e., which experienced CTV syndrome since 2007) than in CS3 and CS4 (i.e., which experienced CTV since 2010 and 2013 ), by the "Peak value" in CS3 and CS4, and by this latter parameter plus "Amplitude" and "Small integrated value" in CS3.
Regarding the second phase, the "Base value" and "Peak value" were recognized as the most representative seasonal parameters of the site-specific conditions in CS3, giving a picture of the timing of CTV-corrective actions and the epidemic evolution.
However, care should be taken in analyzing the linear relationship of the seasonal parameters, which in some cases showed a low magnitude of the R 2 values. In this sense, future research must tend to couple linear and non-linear trends and to analyze seasonal variations, cyclical variations, random or irregular patterns.
In addition, the "Base value" showed a well-defined increasing trend in CS2-CS4 (Figures 9-11), as a function of the changes in NDVI due to the adoption of CTV-corrective actions, even without a significant trend due to the limited numbers of years considered in the second phase. More specifically, these actions were more intense in CS4, followed by CS3 and CS2, as can be observed by the slope terms in Figures 9-11. Therefore, as observed in Figures 8-11, the "Base value" helps to provide information on the evolution of the areas affected by the CTV.
Furthermore, the matrix in Figure 12 shows high R 2 values between the "Base value" and most of the other parameters, suggesting that the analysis of the seasonal parameters for monitoring the areas affected by CTV can be limited to the evaluation of this parameter.
This study contributes to corroborate the significant use of RS technologies in the phytopathology field of application. In fact, RS applications have recently already demonstrated their promising ability for describing and monitoring plant diseases with respect to field-and laboratory-based analyses (e.g., [68,69]). The main limitation of these latter approaches (e.g., including visual inspection for symptoms at early stages and spectroscopic and imaging techniques) is that they are time-consuming, i.e., needing a large number of samples, expensive [70] and not standardized because they are influenced by the observer's expertise. The main advantage of using RS technologies is their cheapness and the ability to provide a large-scale detector for pathogen control [38]. The potential use of time-series RS procedures can be exploited for being adapted in phytopathology studies and being used as a monitoring system of long-term changes [71].
CTV is one of the most complex pathosystems known in plant virology [72]. This complexity is determined by some key factors: the quasispecies nature of viruses, which makes them particularly prone to the continuous evolution of populations of variants under selective pressure operated by the genome of the hosts and the climatic conditions; the vastness of the genetic patrimony of rutaceous hosts; the use of rootstocks with different susceptibility, which determines another important factor of deep interaction with the virus; the quite vast climatic range as well as the soil conditions of citrus growing areas. The number and diversity of these involved factors leads to a very large plethora of disease phenotypes being related to the pathogen. In this context, the effects of CTV on yield are difficult to estimate also due to the fact that, in most phenotypes, the plants died within a few years, as in the case of the CTV epidemic in Sicily, where the syndrome had a rapid evolution with the death of the plants in about 2 years after the infection. In 2002, it was estimated that CTV had caused the death of about 80 million citrus trees worldwide, about half of which were in Spain [73]. When the main symptoms of yellowing foliage or plant decline are visible, the pathogenic process is already at an advanced stage; for this reason, the development of stress detection systems based on RS or proximal sensing techniques is desirable [74,75] because it would allow the early interception of viral infection and the application of more effective containment methods (e.g., elimination of inoculation sources).
The approach developed in this study is valid to be transferred to regional agencies involved in and/or in charge of the management of devastating plant diseases, especially if it is integrated with ground-based early detection methods or high-resolution thermal and hyperspectral imagery acquired by airborne and/or unmanned aerial vehicles (UAV), for the early detection and quantification of plant diseases, (i.e., Xylella fastidiosa [76], Verticillium dahliae infections in Olea europaea L. [77], red leaf blotch in Prunus amigdalus [78]), especially in the case of quarantine plant pathogens requiring complex control schemes that must be readily adaptable to the evolution of epidemics. Furthermore, the application carried out under these heterogeneous conditions of citrus crops can allow the methodology for large-scale interventions and epidemic management programs to be exported.

Conclusions
The adoption of the RS approach described for the characterization of long-term changes in NDVI Terra-MODIS data was performed with the aim of bridging the gap between the temporal distribution of CTV diseases and the recognition of corrective mitigation actions implemented autonomously by growers.
This study outlined the history of the tristeza epidemic in Sicily and showed how the management choices were ineffective because they were not coordinated at the regional level. In fact, the legislative gap of about 14 years between the first version of the Italian Ministerial Decree for compulsory control of the CTV (1996) and the current version, implemented only in 2014, has reinforced the difficulty of farmers in implementing effective operational choices. As a result, fragmentation of interventions is observed, which, if analyzed from an economic point of view, outlines a picture of avoidable economic losses.
The main results obtained from the study are the following: • The Terra-MODIS products represented a valuable data source for the implementation of long-term time series approaches and for the support of phytopathological studies, the limitation of the spatial resolution being largely compensated by their high temporal resolution and the availability of images for long time intervals; • The Terra-MODIS time series approach has proven to be reliable for identifying the specific phenological phases of citrus groves linked to the evolution of CTV, with reference to conditions prior to and subsequent to the implementation of corrective measures by farmers; • Considering TIMESAT statistics analyzed, the "Base value" was identified as a representative proxy for identifying the timing of corrective actions to contain the CTV.
To conclude, it is necessary to underline that the management of severe epidemic events cannot be left solely to individual farmers, but instead, it needs coordinated short-, medium-and long-term interventions. The remote-sensed-based time series analysis proposed in this study is economically advantageous and, if used at the level of territories or regions, it would therefore be useful both in the forecasting phase of the epidemic trend and in the phase of the definition of management interventions.