Pest Detection in Citrus Orchards Using Sentinel-2: A Case Study on Mealybug (Delottococcus aberiae) in Eastern Spain

Abstract

The Delottococcus aberiae is a mealybug pest known as Cotonet de les Valls in the province of Castellón (Spain). This tiny insect is causing large economic losses in the Spanish agricultural sector, especially in the citrus industry. The European Copernicus program encourages the progress of Earth observation (EO) in relation to the development of agricultural monitoring tools. In this context, this work is based on the analysis of the temporal evolution of spectral surface reflectance data from Sen2Like, analyzing healthy and fields affected by the mealybug. The study area is focused on the surroundings of Vall d’Uixó (Castellón, Spain), involving an approximate area of 25 ha distributed in a total of 21 fields of citrus trees with different mealybug incidence, classified as healthy or unhealthy, during the 2020–2021 season. The relationship between the mealybug infestation level and the Normalized Difference Vegetation Index (NDVI) and other optical bands (Red, NIR, SWIR, derived from Sen2Like) were analyzed by studying the time-series evolution of each parameter across the time period 2017–2022. In this study, we also demonstrate that evergreen fruit trees such as citrus, show a seasonality across the EO-based time series, which is linked to directional effects caused by the sensor–sun geometry. This can be mitigated by using a Bidirectional Reflectance Distribution Function (BRDF) model such as the High-Resolution Adjusted BRDF Algorithm (HABA). To study the infested fields separately from healthy ones and avoid mixing fields with very different spectral responses caused by field type, separation between rows, or age, we studied the evolution of each parcel separately using monthly linear regressions, considering the 2017–2018 seasons as a reference when the pest had not developed yet. The observations indicate the feasibility of the distinction between affected and healthy plots during a year utilizing specific spectral ranges, with SWIR proving a notably effective channel, enabling separability from mid-summer to the fall. Furthermore, the anomaly inspection demonstrates an increase in the effects of the pest from 2020 to 2022 in all spectral regions and enables a first approximation for identifying healthy and affected fields based on negative anomalies in the red and SWIR channels and positive anomalies in the NIR and NDVI. This work contributes to the development of new monitoring tools for efficient and sustainable action in pest control.

1. Introduction

The Delottococcus aberiae (De Lotto) (Hemiptera: Pseudococcidae) is a mealybug insect known locally as “Cotonet de les Valls”. It was detected for the first time in 2009 in Benifairó de les Valls (València) in citrus trees and spread to different regions in the north of the province of Valencia and the south of Castellón province [1]. This mealybug (Figure 1) feeds on tree sap, directly weakening the tree, as well causing the deformation of the fruits and a reduction in their size. Cotonet also causes indirect damage due to its molasses excretion, which leads to the development of a fungus on the tree which accumulates on the leaves and can negatively affect the photosynthetic process of the leaves and the quality of the fruits [1]. This pest has several generations a year, being active even in the winter. The population’ maximum may vary in its intensity and its moment of appearance depending on the area and year; it can be detected between 15 June and 30 July, so the critical period for pest control happens just before, from March to June. Sometimes, a second peak is observed between September and October, which can be very high. The most visible cotonet effects occur later, between October and December. With respect to the size of the population detected, the general trend is its increase since 2018, with a considerable increase in populations in 2020 compared to 2019, although in 2021 they were reduced in the areas with the greatest impact from the pest [2].

It should be noted that population level and fruit damage are not well correlated yet. However, the high level of damage at harvest (>20% of fruit attacked) is matched by high maximum annual catches [2].

Mealybug nymphs are usually located in protected areas like the calyxF or between fruits in contact, although they can also be found in branches and leaves. During the adult state, from March to June, some females migrate to the base of the trunk and to the ground to perform the laying of the cottony cover that protects the eggs from desiccation and natural enemies.

The permitted alternatives to chemical control are natural treatments based on sexual confusion, but there are a large number of fields in Comunitat Valenciana that are affected by mealybug populations in which the alternative measures do not show significant improvements [3]. The latest studies carried out by the “Instituto Valenciano de Investigaciones Agrarias” (IVIA) in terms of biological control have shown that the parasitoid Anagyrus aberiae, released by the Universitat Politecnica of Valencia in the spring of 2020, can be effective in controlling the cotonet, since the damage caused by the pest significantly decreased in the fields that were not treated with insecticides and where the parasite was present. However, the long-term effectiveness of that parasitoid has not been determined yet [4]. Given the critical consequences produced by the mealybug on the affected trees, its constant presence throughout the year in affected fields, and the complication in its population control, the mealybug is causing severe economic losses in the citrus industry in the Valencian Community [3]. Within this context, in the area of precision agriculture, the use of remote sensing can play a key role in the modernization of agricultural practices. There are different studies that have demonstrated the usefulness of reflectivity for the detection of diseases, stress, and pests in citrus trees. However, the majority use hyperspectral [5,6] and multispectral [7,8] images acquired with drones (UAVs). Therefore, it is interesting to consider the potential of high-resolution multispectral images from satellites.

The Sentinel satellites of the European Earth observation program Copernicus help control rapid changes in crop production and agriculture at national and field levels [9]. In particular, the Sentinel-2 mission is based on a constellation formed by two identical satellites in the same orbit (Sentinel2A launched in 2015 and Sentinel2B in 2017); every 5 days (since 2017) the satellites cover the entire Earth’s surface at a spatial resolution of 10 to 60 m, providing data in the visible range (400–700 nm), and near- (700–1300 nm) and short-wave infrared (1300–2500 nm). This mission focuses on providing useful information for agricultural and forestry practices, as well as for food security management [10]. Thus, Sentinel-2 provides frequent measurements of the spectral reflectivity of vegetation, which contain information about the pigmentation of leaves, and their cellular structure and moisture content, in addition to highlighting possible anomalies that plants present. In fact, Sentinel-2 was used by our research team in the monitoring of mealybug pest in citrus, showing high sensitivity to red and NIR (near-infrared) reflectivity in fields affected by the pest [11]. In other studies, such as that of M. Yones et al. 2019 [12], spectral differences are identified between several healthy and affected trees such as almond, olive, and lemon trees, demonstrating that it is possible to detect pests in crops using remote sensing. Although Sentinel-2 has a good revisit frequency, clouds can reduce it. As an alternative to improving temporal resolution, the ESA Sen2Like processor is considered, which generates a virtual constellation, merging Sentinel-2 and Landsat at the Sentinel-2 scale, increasing Sentinel-2’s temporal frequency of data by 30% [13].

Regarding plant health, this can be studied, on the one hand, through information on chlorophyll, since an increase in reflectivity in the blue (0.44 $\mathsf{\mu }$m) or red (0.65 $\mathsf{\mu }$m) channels could be a response to a decrease in the concentration of chlorophyll, which can occur due to certain diseases or the presence of invasive insects. On the other hand, the relationship of the NIR (0.7–1.3 $\mathsf{\mu }$m) with internal structure also provides information about the presence of possible stresses; so, if a lower-than-expected NIR reflectivity is obtained, the internal processes of the vegetation may not be working correctly. Furthermore, the SWIR (short-wave infrared, 1.3–2.5 $\mathsf{\mu }$m) presents higher reflectivity outside the water absorption regions (1.45, 1.95, and 2.5 $\mathsf{\mu }$m), when the water content of the leaves decreases, so a high response to SWIR may indicate the presence of stress due to lack of water.

Several studies have reported successful results in identifying different diseases using the NDVI (Normalized Difference Vegetation Index), such as the detection of “Twister” disease in onions [14], and the discrimination of areas with diseased coffee in plantations [15]. In the case of citrus, we also find studies that use the NDVI as a damage assessment tool, e.g., Vanella et al. 2020 [16] for citrus tristeza virus. Many other studies on pests use Sentinel-2 as a database, for instance, to achieve early detection of infestations in forests [17,18], to retrieve pest incidences [19,20], and to make predictions on defoliation [21,22]. Specifically for citrus fields, S2 has been used to study the impact of Huanglonbing disease [23,24]. However, there is an information gap concerning the study of Delottococcus aberiae through remote sensing.

The vegetal covers formed by fruit trees (such as orange trees) present heterogeneity due to their geometry at the tree level and the presence of soil (or vegetal cover), depending on the tree separation, the line orientations, and the trees’ age and trimming. All these contributions generally cause artificial artifacts, generating high (or low) surface reflectance values. Fewer shadows are caused by a low SZA (higher solar altitude) and the ground within the pixel size of Sentinel-2 data takes more prominence in the reflectivity; on the other hand, for a high SZA (illumination closer to the horizon), the generation of shadows increases, leading to lower surface reflectance values. The Bidirectional Reflectance Distribution Function (BRDF) is a function that describes how each surface interacts with incident radiation and reflects it back to a sensor in a given angular geometry [25]. Thus, it determines the geometric characteristics of the spectral reflectivity, being fundamental for its spectral and directional description. In the case of heterogeneous surfaces, the BRDF can present anisotropies, whose main directions are driven by the morphology of the cover, its optical properties, and the illumination zenith angle (SZA) [26]. In fact, Wang et al. demonstrated that there is a relationship between BRDF behavior and the canopy height and soil gaps [27]. Correcting these BRDF effects is necessary to normalize the sensor data, especially when analyzing the phenological evolution [28]. Several studies have demonstrated that BRDF correction improves their data quality, such as Colgan et al. with the classification of tree species in the Savanna [29] and Jia et al. in classification with hyperspectral imaging [30]. However, the correction of BRDF effects in medium- to high-resolution satellite images is really challenging given that their observation geometry is nearly constant. However, Franch et al. [31,32] presented the High-Resolution Adjusted BRDF Algorithm (HABA) method to normalize the BRDF effects for Landsat and Sentinel-2, and this method can be implemented in the Sen2Like product [33,34]. That is why the HABA algorithm [31] is used to minimize these angular effects. On the other hand, the literature does not include studies that apply this correction to citrus trees, an issue that may impact the development of EO-data-based algorithms that are transferable over time.

Given the limited number of studies that have analyzed the incidence of the mealybug pest in orange trees using remote sensing, the main goal of this study is to analyze how remote sensing can be used to detect it. Particularly, the main science questions that we want to address in this work are as follows: How can optical remote sensing data be used to detect the pest? How can we leverage each spectral region within the solar spectral range to detect it? When is the most significant period in the season for detecting it?

The present work has the aim of characterizing how the incidence of the mealybug affects different regions of the orange trees’ spectral response, as well as analyzing which range is optimal for the separation of healthy and affected fields.

2. Materials and Methods

2.1. Study Area

The study was carried out using 24.6 Ha distributed in 21 plots located in the vicinity of the municipality of Vall d’Uixó (Castelló, Valencian Community, Spain). This zone has a Mediterranean climate, with mild winter temperatures and hot, dry summers, with an annual average temperature of 22°. The hot season begins in June and generally lasts until September, with an average maximum temperature of 28°. The cold season runs from November to March, with an average minimum temperature of 7°.

Figure 2 shows the spatial distribution of the parcels. The plots are divided into two groups: 12 affected and 9 healthy (red and blue, respectively, on the map in Figure 2). The area of the fields affected by cotonet covers 14.58 ha, while the healthy ones occupy 10.04 ha. This classification was made for the 2021–2022 season, and it was based on the amount of harvested fruits which were damaged by the pest. The fields with more than 50% damaged fruits were considered as affected. This information was provided by the San Vicent Ferrer and Vallfruit cooperatives, located in the Vall d’Uixó. Furthermore, these local sources reported a substantial increase in the pest population in 2019 compared to previous years. Through the considered dataset, each parcel may have a different age as well as different variety: Clemenules, Lanelate, Villa, Murco, Hernandina, Clemenvilla, Clemenrubí. It is essential therefore to consider that the phenology of the studied fields may be different in each case, regardless of their condition. Furthermore, each plot is susceptible to different treatments, geometries, and different tree arrangements.

2.2. Sen2Like Processor and NBAR

For agriculture applications, when there is a need to work with fine spatial resolution and with a medium/high revisit time, the Sen2Like (version 4.3.0) processor is chosen, as it combines the Sentinel-2 MSI and Landsat-8 OLI L1C data with the aim of creating a mosaic-based multi-temporal surface reflectance dataset [13,34]; it includes geometric and atmospheric corrections. For the atmospheric corrections, the Sen2Cor 3.1.0 [35,36] algorithm is applied, and for dealing with the presence of clouds, the SCL mask is used.

Since Sen2Like combines images from different sensors, variations in observation geometry and lighting conditions, along with the heterogeneity of the analyzed vegetation cover, together induce variations in reflectance [31]. Gao [37] concluded that for Landsat-like narrow swath sensors, the major BRDF effect arises from the day-of-year effect, and can cause variations of 0.04–0.06 reflectance compared to mid-summer observations. That is why a correction of the reflectance’s angular effects is needed. The product normalized with the BRDF is called the Nadir BRDF-Adjusted Reflectance (NBAR) [13]; in this case the annual mean solar zenith angle of each tile is used as reference angle.

The processed tile was the 30TYK, using as a time period the years between 2017 and 2022, inclusive.

2.3. Preprocessing

In our case, given the background, the regions of the red, NIR, and SWIR are analyzed. In addition, we also implement the calculation of the NDVI since it relates the NIR and red bands and some studies have shown good results when using this index to study other pests [14,15,16] (Equation (1)).

$$\mathrm{NDVI}={\displaystyle \frac{\mathrm{NIR}-\mathrm{red}}{\mathrm{NIR}+\mathrm{red}}}$$

(1)

Given that the spatial unit where the data are collected is at parcel level, first of all we extracted the polygons that limit each parcel from SIGPAC [38].Then, a morphological erosion was applied to each parcel to subtract their edges and the possible outliers were filtered using a manual threshold method for the reflectivity of each channel.

In this work, for simplicity we focus the study on three different bands, considering each of them as representative of the different spectral regions of the solar spectral range:

• Band 4 (red): (0, 0.15) $\mathsf{\mu }$m, as representative of the visible spectral range;
• Band 8A (NIR): (0.2, 0.5) $\mathsf{\mu }$m;
• Band 11 (SWIR): (0.1, 0.3) $\mathsf{\mu }$m;
• NDVI: (0.2, 1).

As the ground data are provided at the parcel level, not pixel level, the entire region of a field cannot be considered affected with the same impact. Hence, in order to to homogenize the input data, the calculations of the average and standard deviation per band and NDVI are performed, first at a parcel level and then at a group level (healthy–affected).

2.4. Main Process

The average (${\mu }_p$) and standard deviation (${\sigma }_p$) of the NBAR (and NDVI) of each parcel in each available image (from 2017 to 2022) were obtained using Equations (2) and (3), respectively.

$${\mu }_p={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_i$$

(2)

$${\sigma }_p=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{({\rho }_i-{\mu }_p)}^2}$$

(3)

where ${\rho }_i$ refers to the reflectance per pixel, and i and n represent the amount of pixels in parcel p. Then, the average ($\mu$; Equation (4)) and standard deviation ($\sigma$; Equation (5)) of both parcel groups, the healthy and the affected, were calculated.

$$\mu ={\displaystyle \frac{1}{N_p}}\sum _{i=1}^{N_p}{\mu }_p$$

(4)

$$\sigma ={\displaystyle \frac{1}{N_p}}\sum _{i=1}^{N_p}{\sigma }_p$$

(5)

where $N_p$ is the number of parcels in each group.

The temporal evolution from 2017 to 2022 of all the healthy and affected averages, considering parcel-by-parcel means and group means (${\mu }_p$ and $\mu$) was determined. In the case of the groups, the differences between the means of the healthy and affected group averages were calculated. This was performed for all the channels mentioned. This method was motivated by the difficulty of intercomparing parcel signals, as each one has its own characteristics (variety, phenology, treatments…). Figure 3 shows the temporal evolution of the means by plot; the impossibility of finding clear separability by simple observation is evident, because the signals of healthy and affected parcels are superposed. For the red, NIR, and NDVI cases, mixed signals also occur. In this figure, a clear seasonality can also be perceived, where reflectivity increases (summer) and decreases (winter) periodically.

In order to avoid seasonality, all the time series were grouped by months, constructing the monthly means per year at a parcel level. In this way we were able to study the evolution of each month of the year separately, without taking account of these maximums and minimums that influence the analysis, since all years have the maximums and minimums located on the same dates. As local sources indicated, it can be assumed that in 2017–2018 the parcels were not significantly affected by the pest, so this period is considered as a reference. A linear regression over the monthly means (from 2017 to 2022) was performed, obtaining the slopes, intercepts, the Pearson correlation coefficients ($r_{value}$), $p_{values}$, and the standard error of the estimated slope, under the assumption of residual normality. This method focuses on the analysis of the slopes, and compares healthy and affected plots by using the monthly slope evolution at a parcel level.

Due to this obstacle, building an algorithm to separate healthy and affected signals becomes more complicated. The monthly regressions can be a helpful technique to take into account the parcel’s historical signal and analyze how each parcel has changed since it was healthy in 2017.

Furthermore, to help in the parcel by parcel evolution analysis, the differences with respect to 2018 were calculated, and are denoted as anomalies. Equation (6) shows the method for obtaining these anomalies ($\Delta$), where ${\rho }_{18m}$ refers to the reflectance mean of month “m” for the year 2018 and ${\rho }_{Xm}$ refers to the reflectance mean of month “m” for year X. In this case, 2018 is used as the reference year because in 2017 Sentinel-2 did not yet have its two platforms up and running.

$$\Delta ={\displaystyle \frac{{\rho }_{18m}-{\rho }_{Xm}}{{\rho }_{18m}}}$$

(6)

Figure 4 shows a summary of the methodology workflow.

3. Results

3.1. BRDF Correction and NBAR Time Series

Figure 5 shows the directional surface reflectance (SDR) (red) and the NBAR (blue) dependency with the SZA, focusing on the red, NIR, and SWIR bands, averaging the signal of all the studied fields and across 2022.

The NIR and SWIR, as shown in Figure 5b,c, present significant differences between the NBAR and SDR signals, of more than 0.3 reflectance points. When the sun is higher (SZA is lower), fewer shadows are generated, when the sun goes down, the generation of shadows increases, especially considering the morphology of the orange trees. The raw signal is initially higher than the corrected one in all three cases; so, the reflectance is overestimated. On the other hand, at higher SZAs (lower sun), the SDR signal is lower than the NBAR, since the reflectance is underestimated due to the presence of shadows.

Furthermore, Figure 6 shows the temporal evolution of the SDR (red) and NBAR (blue) for the $\rho NIR$ mean of healthy parcels (Figure 6a) and affected ones (Figure 6b), as well as the evolution of the SDR and NBAR signal differences for healthy (Figure 6c) and affected (Figure 6d) parcels too. In this comparison, it is notable, on one hand, that SDR and NBAR series present the same behavior, as both are provided by periodic signals which accompany the crop cycle, with peaks in summer and minimums in winter. But, on the other hand, there are no relevant variations when comparing the evolution of healthy and affected mean parcels $\rho$. As shown in Figure 6c,d, the evolutions of healthy and affected parcels are practically identical.

Thus, given the clear differences between the SDR and NBAR signals, given that these differences are independent of whether the plots are healthy or affected and the reduction in the seasonality when using the NBAR reflectances, the NBAR is used throughout all the following development.

3.2. Temporal Evolution

The evolution of the means is shown in Figure 7 and the differences between healthy and affected signals can be observed in Figure 8. From these we can extract various information. Regarding the red channel, its time evolution (Figure 7a) and difference (Figure 8a) evidence that there is a cyclic behavior, where the NBAR red presents maximums located in the months of June–July and minimums in December–January, coinciding with the end of the orange tree growing cycle. It can be appreciated how initially the response of the healthy parcels is greater than that of the affected parcels, but as time passes, the curves become closer little by little. The NIR evolution and difference (Figure 7b and Figure 8b) also manifests a periodic behavior, which moves according to the crop cycle. It is shown that the maximums and minimums occur at the same time as in the red case. The main period of interest (2020–2021) presents positive differences, which implies that the NIR NBAR of the healthy parcels is greater than the NIR NBAR of the affected ones. Furthermore, a periodic signal following the crop cycle can also be observed in the temporal evolution of the SWIR means (Figure 7c), with peaks in summer and troughs in winter. The differences (Figure 8c) between the healthy signal and the affected one become smaller as time progresses, even becoming negative in the 2020–2021 cycle. The temporal evolution of the NDVI (Figure 7d) also shows cyclical behavior, but with the minimums (June–July) and maximums (December–January) reversed with respect to the previous cases. Regarding the NDVI differences (Figure 8d) between healthy and diseased plots, it can be observed that the initial NDVI of the healthy areas was lower than that of the affected areas. As time progresses, the pest influences the affected plots and the differences between the NDVIs of the two groups are shortened.

3.3. Temporal Tendency Analysis

The visual inspection of the evolution of the time series is not precise enough to establish a clear separation between healthy and affected parcels. Seasonality and phenological differences between parcels may hinder the analysis; so, the slope method was implemented. This is considered as a way to determine the monthly evolution of the parcels without the need to intercompare them, in the case of analysis per parcel. A box plot representation was chosen to show the slopes of the monthly regressions (Figure 9). This helps in the separation when parcels’ time evolutions are regressed monthly one by one. In this case, the SWIR and the NDVI present the best visual separability. For SWIR (Figure 9c), the box plots of affected parcels present higher values than healthy ones. According to visual inspection, separation is seen in the second part of the year. Red (Figure 9a) and NDVI (Figure 9d) do not show a distinction for classification. The separability detected in May in the NIR range (Figure 9b) is also interesting, because the pest starts increasing its population in May–June.

All the results are summarized in

Table 1. Summary of results.

Spectral Range	Periodicity		Trend Sign	Separability
	Minimum	Maximum	Difference	Months
Red	Winter	Summer	Negative	8, 10, 11, 12
NIR	Winter	Summer	Positive	1, 5, 11
SWIR	Winter	Summer	Negative	From 6 to 12
NDVI	Summer	Winter	Positive	7, 8, 10, 11, 12

.

3.4. Anomaly Analysis

Figure 10, Figure 11, Figure 12 and Figure 13 show the anomalies from both healthy and affected parcels, comparing the monthly mean of each year’s signal (from 2019 to 2022) with 2018. The standard deviations of the affected parcel means are higher in all cases. The magnitude and sign of the anomalies change depending on the channel.

Regarding the red channel, in the case of healthy plots (Figure 10a), the anomalies are positive throughout the year, across all years. On the other hand, affected plots (Figure 10b) show negative anomalies at the beginning and end of the year, starting their negative values in August, except in 2019 which presents a maximum in November. The magnitudes of the anomalies are similar in both cases, generally between −0.2 and 0.2. For both healthy and affected plots, 2019 shows positive anomalies at the end of the year. In addition, 2022 stands out for its peak in June, which decreases to show negative anomalies from August onward. It is also the year with the largest anomaly. As for 2020 and 2021, both years show very similar behavior.

The NIR band shows lower anomalies than in the red channel. The healthy parcels (Figure 11a) present tiny anomaly values, with positive values lower than 0.1, while affected parcels (Figure 11b) show higher and positive values, excepting the end of 2019, which has negative values. In addition, the deviations of the anomalies of the affected plots are higher again. Regarding the annual intercomparability, this time all the years behave in a very similar way, except 2019.

In healthy parcels (Figure 12a), the SWIR shows positive values, except at the end of 2019. The affected parcels (Figure 12b anomalies) present larger negative values, especially in May of 2020, but show a large standard deviation, and from September to January of 2020, 2021, and 2022, with the latter showing the lowest value.

With respect to the NDVI, in general, the healthy parcels (Figure 13a) present small values in magnitude with a negative sign, excluding in 2019, when they display a positive maximum in October. The affected parcels (Figure 13b) show higher values, which are mainly positive, especially at the end of the year. The year 2022 shows positive anomalies through the year. Finally, 2019 shows a maximum anomaly both in healthy and affected parcels in October. The years 2020 and 2021 show similar behavior, with positive anomalies starting in July.

To give a spatial conception of the possible classification using the anomalies we constructed the map in Figure 14. It shows the anomalies per plot for the SWIR channel and the month of October 2022. As can be seen, in general, the affected plots (red border) show more negative anomaly values than the healthy plots (blue border), which show higher values.

4. Discussion

In this work, we analyzed the temporal evolution of the surface reflectance retrieved by Sentinel-2 and Landsat to monitor several orange tree fields affected by mealybugs and compared them with healthy parcels. First, citrus trees’ surface reflectance time series showed a strong seasonality, showing maximum signals during June–July and minimum signals during December–January across all spectral bands. In this work, we linked this seasonality to directional effects caused by the illumination and observation geometry. However, corrections of these directional effects are very challenging for high-resolution sensors such as Sentinel-2 or Landsat, whose sun-synchronous orbit provides near-nadir observations at the same time of the day. The HABA BRDF normalization method shows a generalized reduction in said seasonality. However, the NBAR signal still shows a remaining seasonality; this can be linked to limitations (i) of the BRDF algorithm (Ross Thick Li Sparse reciprocal model [39] to capture the particularities of evergreen fruit trees such as row orientations or the shape of the crowns (as stated by Lucht and Roujean, 2000 [40]; (ii) the inversion of the BRDF coarse-resolution pixels over heterogeneous areas such as the study area, which is characterized by small-size fields that may mix different crops, varieties, or row orientations; or (iii) the transition from coarse to medium resolution, that is, the HABA algorithm. Future studies will investigate all these possibilities to improve the performance of the HABA algorithm over these heterogeneous sites. Despite these limitations, in this work we highlight the need to consider these directional effects when working with agriculture applications over evergreen fruit trees when targeting robust algorithms that can be transferred across seasons.

Another possibility is that the remaining seasonality could be linked to any phenological pattern of the trees. Focusing individually on each spectral band, the NBAR red band shows a strong correlation with the SZA and barely reduces its angular dependency. It is interesting to see how the descending (from maximum to minimum) of the red NBAR’s temporal evolution coincides exactly with the growth of the fruit, which occurs from August to December. This phenomenon could be caused by the presence of fruit on the tree directly altering the reflectivity of orange trees. In the visible range, the reflectance spectrum is driven by the absorption derived from photosynthesis; therefore, the transmissivity is limited, causing the canopy elements that are localized underneath (oranges) to contribute minimally to reflectivity measured by the sensor, which only detects a signal from the upper parts of the tree [41]. As we mentioned in the Introduction, the red band is strongly related to the chlorophyll content; therefore, the observed cyclical changes may be influenced by variations in this. Furthermore, the summer maximums could be due to a lesser presence of leaves at this stage, since the leaves of the spring (April) budbreak have not yet grown enough. In addition, young leaves (3 to 6 months) have a greater photosynthetic capacity than those which are older than one year [42]; therefore, the decline to the December minimum could be a consequence of the increase in photosynthetic capacity of the leaves from the spring budbreak.

With respect to the monthly regressions, the red channel does not show very good results, since the slopes overlap in all months of the year, except in December, a month that presents the clearest difference between the healthy and affected groups. Unfortunately, at this time of year the damage caused by the pest is already very advanced.

Regarding the NIR, it is the spectral band with the largest reduction in the directional effects when applying HABA, decreasing the R2 and the slope significantly when evaluated against the SZA. However, the NBAR product still shows a remaining seasonality. The NBAR signal is initially similar for the two parcel groups, which could mean that the differences between these (variety, age, humidity, land conditions, etc.) are lower at the structural level than at the chlorophyll concentration level. In addition, the signal is cyclical with time, with minimums and maximums located around December and August, respectively, as in the previous case. In the NIR region, the radiation experiences so-called dispersion by volume, where an important part of the incident energy is transmitted through the upper leaves (from the top of the canopy to the bottom), and then it is reflected from the lower ones retransmitting (from below to the top), increasing the NIR NBAR of the tree. The presence of oranges can hinder this NIR transmission inside the tree, reducing this dispersion by volume, and therefore the measured NIR NBAR. This could be a good explanation for the fact that with the thickening of the fruit the NIR NBAR is reduced even though the spectral response of the fruit itself presents elevated NIR values [41].

It is also seen that as the parcels become affected, the NIR NBAR of the healthy plots increases with respect to that of the affected ones, with the trend of differences being positive with time, this may be the consequence of a reduction in the vigor of the plants directly induced by the disease [43]. The slopes method evidences separability in May, which is interesting because the population peak occurs in this month. In addition, from April to June, field sampling is performed to determine which plots it is necessary to treat. Thus, if it is possible to provide a remote indicator from satellites, field efforts could be partially avoided.

In regard to the SWIR, the NBAR product shows a significant reduction in the seasonality but after the BRDF normalization there is a remaining variability. As we mentioned in the Introduction, this is correlated with the water content of the leaves. Thus, a plant with elevated water stress will provide a higher SWIR signal (SWIR NBAR) than another plant with a higher water level. Therefore, the decline in SWIR NBAR as the fruit ripens (between September and December) could be related to the fact that oranges contain more water than leaves and the water’s presence increases the total absorption of the tree in the SWIR band, decreasing the canopy reflectivity in this spectral region [41,44]. This phenomenon is also reflected in the slopes method applied to the SWIR region. The fact that the SWIR response in affected parcels has been rising since 2017 shows the presence of certain stress due to the presence of mealybugs. There are several studies which have found a correlation between crop stresses and the SWIR region [45]. This spectral region has similarities to the thermal infrared, and in fact, studies such as the one in 2016 by Anderson et al. [46], conclude that stress can be detected earlier in the thermal than in the optical, since the first behavior of vegetation in the face of stress is to close its stomata, causing the temperature to rise. The parcel-by-parcel analysis reflects the behavior of the stress incidence, there are perceptible differences between healthy and affected groups. As mentioned in the Introduction, as the SWIR response rises, more signs of stress are presented by vegetation, in accordance with the obtained results.

Furthermore, the great separability observed in August is a valuable result, since for the farmers, the damage in the fields becomes visible in November–December (IVIA); therefore detecting the infestation in August could mean a warning of 3 months. In addition, this warning would come at the same time that the treatments are carried out (July, August, and September), meaning farmers could focus them on the most problematic plots. Likewise, the separability throughout the second half of the year is interesting, because it will be possible to see if the treatments applied have worked correctly without having to visit the field.

The NDVI differences present a positive trend (positive slope of the linear regression; with the highest one 5.65 $\times {10}^{-10}$), which may imply that the vigor and health of the affected plots has reduced over the years. A high dispersion is seen in the data, as evidenced by the higher RMSE (RMSE = 2.34 $\times {10}^{-2}$), but this value must be analyzed relatively, as the magnitude of the NDVI is greater than the spectral channels; therefore, its RMSE will also be higher.

On the other hand, the NDVI slopes (Figure 9d) from healthy parcels present larger values than from affected parcels, where they can even be negative. This fact indicates that the NDVI of the affected plots does not increase, and even decreases (with negative slopes), which may be a signal of stress in vegetation. Additionally, the NDVI shows limited separability and the differences between slope distributions can be noticed in August and December.

It is important to note that each parcel must be studied separately to perceive the healthy and affected parcels’ slope differences, given that the phenology, varieties, and treatments may be diverse, and also the incidence of the infestation, since it is variable between years. However, analyzing the trends may mask years with high incidence depending on when it occurred within the time series analyzed.

The analysis of the slopes based on the monthly regressions shows a first approach of a separability criterion, focusing mainly on the SWIR and NIR bands. The definition of a specific value is complicated in this preliminary study given the limited statistical significance of the sample, but it demonstrates that a criterion could be established around the separation between the medians of healthy and affected parcels. Also, given that we have presented the data as box plots, this criterion could be characterized by the groups’ differences between the first and the third quartile, which could be positive or negative, depending on the spectral range.

The anomaly analysis is proposed as a possible first way to obtain a quantitative method for distinguishing between affected parcels and healthy ones. It is also proposed as a way to make year-to-year comparisons. On the other hand, the analysis of anomalies, although starting from a pre-infestation year, allows us to study the evolution on a monthly level year by year, and to compare the infestation incidence in each of them.

Following the relationships between reflectance and the biophysical parameters of trees, as described in the Introduction, and in a similar way to the slope method and the temporal evolution, we can intuit how the evolution curves of the monthly anomalies of each band will theoretically behave based on their own definition. Theoretically, for the red channel, negative anomalies mean that the reflectance in the red has increased compared to 2018, meaning chlorophyll decreased. This phenomenon is indeed reproduced, since negative values are found only in the case of affected plots. In fact, as shown in Figure 10b, most of the negative values occur at the end of the year where the effects of the mealybug are visible. Something similar happens with SWIR, where negative anomalies can be associated with tree stress, because as SWIR increases, so does the hydric stress, and Figure 12b shows more negative values than Figure 12a, confirming the theory. In contrast, positive anomalies in both the NIR and NDVI channels may be associated with lower health given their relationship to plant volume and vigor, since the higher the anomalies, the lower the NIR (or NDVI) compared to 2018. Actually, anomalies in the NIR channel and NDVI (Figure 11b and Figure 13b) from affected parcels display higher values than the evolution of their respective healthy parcels (Figure 11a and Figure 13a).

The year 2019 shows similar behaviors in all cases, whether the parcels are healthy or not. This may indicate that the pest was not as influential in this year, or that the contextual conditions (especially in the second part of the year) were very important. Such external influences could be climatic parameters such as precipitation or temperature, which can affect the crop cycle.

The pattern of the red anomalies in 2022 is remarkable, since they decrease very clearly from August. Even though it is a year with high variability. The negative anomalies of the SWIR in the same year also stand out from August. On the other hand, the anomalies of 2020 and 2021 are similar, especially in the red, NIR, and NDVI channels, and in the second part of the year for SWIR. This phenomenon could be due to the mealybug populations increase in 2020 compared to 2019. In 2021, although there was an extension in the affected areas, those areas where there was already a high level of infestation sometimes reduced their population [2]. This could explain why the difference between the anomalies of 2020 and 2021 is not very high.

In general, differences between healthy and affected plots are more pronounced in the second part of the year, while there is little difference in the first part of the year. This is consistent with the timing of pest emergence, as pest populations begin to increase in June [2]. In addition, according to the hypotheses related to the biophysical parameters discussed above, interannual differences are evident in the second part of the year. In all spectral regions, there is a gradient of areas being increasingly affected from 2020 to 2022, which is consistent with the observation of a considerable increase in mealybug populations from 2020 onward [2].

It is also interesting to note the difference between the standard deviations of the healthy and affected anomalies. Anomalies from affected plots show higher standard deviations, which could indicate that the damage level is not homogeneous. As mentioned in the Introduction, mealybug population levels vary by zone and year.

Thus, the study of the temporal evolution of monthly anomalies could be a good tool to separate healthy and affected plots for each season. By establishing a boundary that will vary according to the channel, it may be possible to identify affected plots. We could summarize it in a first general criterion to characterize the affected plots. Regarding the red and SWIR channels, affected parcels are characterized by negative anomalies from July onward. On the other hand, NIR (greater than 0.05) and NDVI positive anomalies are related with a greater incidence of pests.

The used database is limited: there is very little statistical evidence, and it is all based on a single observation from the 2020–2021 season. The impact data were measured on the fruit, not on the trees, which may not correlate with the satellite information. Therefore, in future studies we will try to expand our database, both spatially and temporally.

5. Conclusions

The present work focuses on the study of the Delottococcus aberiae plague on orange trees by the utilization of remote sensing images at a 10 m spatial resolution. Images from the Sen2Like processor are used to analyze several parcels located around La Vall d’Uixo (Castello, Spain).

In the development of this research, some complementary methods are implemented: On one hand, a comparison analysis of the spectral reflectance (and NDVI) time series from healthy parcels and those affected by mealybug, focusing on the red, NIR, and SWIR channels, and also using the NDVI. On the other hand, a monthly linear regression is performed to reduce seasonality and be able to explore and isolate the tendency of the reflectance evolution of each month. Finally, an anomaly study is effected to provide alternatives in the search for quantitative outcomes and allow a year-on-year comparison of the evolution of the pest.

From the obtained results, several conclusions can be extracted:

• Reflectivity has a strong dependence on SZA in all the studied bands, but especially in the NIR and SWIR channels. The HABA BRDF normalization is applied in order to reduce that relationship. Furthermore, the behavior of cyclical reflectivity is related to the crop cycle, resulting in a clear seasonal pattern; the BRDF normalization also helped mitigate that seasonality.
• The separability method based on monthly regressions is able to differentiate between healthy and affected plots, especially in the SWIR range, at different times of the year, given that positive slopes are related to affected fields and negative slopes to healthy ones from August onward. This is especially interesting because of the huge differentiation observed and its coincidence in time with the treatments application, which are used from July to September [47]. The slight separability detected in May in the NIR is also remarkable, since in April, May, and June, control samplings are carried out to see if there is the pest is present; being able to detect problematic plots in May from satellite could facilitate this arduous task.
• The evolution of the differences between the daily averages for healthy and affected plots shows negative tendencies in the SWIR, which could imply that the SWIR reflectivity average is increasing in the affected fields, which could be related to an increase in water stress.
• Thanks to the inspection of the anomalies, it is possible to compare each year with respect to the others; this shows an increase in the affected level from 2020 to 2022 in all spectral regions. The second part of the year is again the period with the best separability between healthy and affected plots.
• A preliminary criterion of separability can be established according to the anomalies, since it is seen that the affected parcels stand out due to negative anomalies in the red channel (beginning and end of the year) and the SWIR region (at the end of the year), as well as positive anomalies in the NDVI and in the NIR band (greater than 0.05).

These results indicate that spectral differences due to the presence of mealybug pests can be found. The mentioned conclusions are a starting point. It will be necessary to improve the analysis by extending the database to other areas and increasing the temporal sampling, because one of the most conditioning factors is having infestation data from just a single time point (observation made in the 2020–2021 campaign). Furthermore, it would be optimal to have quantitative sampling of a parameter that can be related to the condition of the pest (such as the count of individuals), instead of having a global field classification based on just one campaign. It would also be interesting to extend the methodology to other crops (such as persimmons), and to explore the full spectrum (hyperspectral) to characterize the spectral signatures of healthy orange trees and those affected by the pest.