Use of Geostatistics for Multi-Scale Spatial Modeling of Xylella fastidiosa subsp. pauca (Xfp) Infection with Unmanned Aerial Vehicle Image

Abstract

In recent years, the use of Unmanned Aerial Vehicles (UAVs) has been spreading widely, as in plant pest control. The collection of huge amounts of spatial data raises various issues including that of scale. Data from UAVs generally explore multiple scales, so the problem arises in determining which one(s) may be relevant for a given application. The objective of this work was to investigate the potential of UAV images in the fight against the Xylella pest for olive trees. The data were a multiband UAV image collected on one date in an olive grove affected by Xylella. A multivariate geostatistics approach was applied, consisting firstly of estimating the linear coregionalization model to detect the scales from the data; and secondly, of using multiple factor kriging to extract the sets of scale-dependent regionalized factors. One factor was retained for each of the two selected scales. The short-range factor could be used in controlling the bacterium infection while the longer-range factor could be used in partitioning the field into three management zones. The work has shown the UAV data potential in Xylella control, but many problems still need to be solved for the automatic detection of infected plants in the early stages.

1. Introduction

Remote sensing technologies have evolved over the years, and the agricultural sector now has a multitude of options both in terms of platforms (from satellite to manned aircraft to UAS) and sensors (e.g., visible, multispectral, hyperspectral, and thermal) with which to collect a variety of useful data to support farmers’ informed decision making. Despite the wide applications of satellite-based RS in agriculture [1,2], its use still remains limited, especially in precision farming due to the coarse spectral resolution and the still high cost of services. Furthermore, the use of satellite data is severely limited by weather conditions (e.g., cloud and snow cover).

UAS platforms are a promising alternative solution for obtaining high-frequency data with appropriate spatiotemporal resolution, so that plant health can be monitored according to the farmer’s needs [3]. However, even UAS applications may be limited by weather conditions (e.g., a maximum wind resistance of 10 m/s precipitation), limited spatial coverage due to limited battery life and/or regulatory restrictions, as well as the maximum payload, which limits the simultaneous use of multiple sensors [4].

In any case, the use of UAVs is expected to increase in Precision Agriculture because there is a growing need by farmers to increase spatial resolution and to reduce the time interval between repeated measurements to monitor the health status of the plant, especially during the development of an infection [5]. Drone data are particularly suitable in this regard, as most experiments indicate that higher spatial resolutions lead to better results, since more meaningful information can be extracted when plants and/or leaves can be analyzed individually [6,7].

Xylella fastidiosa (Xf) is a well-known plant pathogenic bacterium that is transmitted through insects and causes serious infections affecting more than 350 plant species. In olive trees, it causes a rapid decline syndrome (OQDS) that currently poses a serious threat to the survival of more than 360 kha of olive groves in south-eastern Italy (Apulia region). Although controls and countermeasures are currently implemented by the regional administration to limit the rapid spread of the infection, the most effective measure for such containment is the early detection of the presence of the bacterium even on visually asymptomatic plants.

To date methods for detecting the presence of the bacterium rely on complex laboratory measurements with medium to long waiting times for results. Recently, alternative techniques have been proposed that involve the use of a UAV equipped with a multispectral radiometer [8,9,10,11,12].

The use of such modern technology has led to a proliferation of data and the development of methods for their analysis. While this has brought undisputed advantages, it has also raised both theoretical and practical problems that need to be resolved. Some theoretical issues include scale, which is of crucial importance in geospatial research and relevant applications even in Precision Farming (PF). Scale refers to the spatial extent to which certain processes operate in the environment and it is then quite important to identify this scale in view, also, of the application of the collected spatial data. Fine resolution is associated with more detailed information about the processes under investigation. However, this is not always the optimal solution in terms of costs–benefits, since the optimal scale should be defined according to the objectives of the survey and its application. Spatial research suggests matching the scale of spatial analysis to that of the process dynamics under study [13,14]. Anyway, a given process can evolve across multiple scales, with complex dependencies between scales [15,16]. It is then worth emphasizing the need for introducing the scale within the spatial data modeling.

From a theoretical point of view, one can state that scale is related to the concept of spatial dependence, which is an intrinsic property of many spatial variables. Spatial dependence refers to the property of neighboring observations to be similar and to influence each other, according to Tobler’s (1970) [17] definition: ‘Although everything is related to everything else, neighboring things tend to be more related than distant things.’ The consequence is that neighboring observations cannot be considered independent and randomly distributed but are characterized by varying degrees of spatial association. Therefore, the pixels of a high-resolution UAV image (of the order of centimeters) can no longer be treated as independent and, consequently, the techniques of traditional statistics should be suitably modified to account for spatial correlation.

Furthermore, when spatial phenomena occur over a wide range of scales with a nested or hierarchical structure [14], such as those explored by an RS image, it is necessary to identify the scales to decide which one(s) to analyze and represent for a specific application. There are several statistical methods for assessing spatial autocorrelation, among which is geostatistics, which is a robust system for quantifying and modeling spatial dependence [18]. Geostatistics most commonly uses variograms to measure the dissimilarity between two sample locations as a function of the distance between them. From the mathematical form of the variogram model and its parameters, in particular, the range, it is possible to make useful inferences regarding the scale(s) present in the spatial data. More specifically, the range represents the maximum scale of spatial autocorrelation, beyond which observations can be assumed spatially independent as in traditional statistics. However, other smaller scales could be included in a nested variogram model, formed by the sum of various mathematical models, each with its own range representing a different scale of spatial variation [19]. An important consequence is that the variogram can be used as a useful geostatistical tool to determine the multiple scales of spatial variation in a set of spatial data.

From all the above it is therefore evident that there is a relationship between spatial dependence, as expressed by a variogram, and a set of spatial scales relevant to spatial modeling. Thus, the variogram(s) of any variable or set of variables in a multivariate approach can be used as an analytical tool to inform multiscale spatial modeling in geostatistics. A key question then arises whether in spatial modeling it is more appropriate to use a single relevant spatial scale, or rather one should provide multiple representations of the same scene at the different scales corresponding to those identified from the spatial data. We believe that it is more reliable to produce multiscale descriptions, since the spatial relationships among variables vary as a function of scale [20]. Once a decomposition of the overall spatial variance into its various spatial components has been obtained, the one specific to the phenomenon under study and/or concerning the application of interest can be selected.

The objective of the paper is to show how it is possible to treat multi-scale spatial modeling using multivariate geostatistics and to test the potential of the UAV in combination with a multispectral radiometer in detecting Xylella fastidiosa subsp. pauca (Xfp) at a very early stage over an olive grove in the region of Apulia (Italy).

2. Materials and Methods

2.1. Study Site

The olive grove consists of 50–60-year-old trees of the “Cellina di Nardò” cultivar arranged in a regular sequence (Figure 1). The trees have a mainly vertical development bearing with large crown. The trunk is regular, fairly linear, and columnar in shape and the canopy is thick and shrubby. Periodic plowing was carried out on the plot to contain weeds and the plants were not irrigated.

2.2. Visual Inspection of Plants for Infection Detection and Assessment

Each olive tree was surveyed to assess the level of severity of Xfp infection. The visual inspection of symptoms with related photographic documentation was carried out by a well-trained operator who recorded the percentage of foliage with wilt symptoms in each of the four cardinal positions. Plants were grouped according to an empirical scale of symptom severity with increasing values from 1 (asymptomatic plants) to 6 (plants characterized by the absence of green parts and therefore dead) [9].

The assessment of the level of severity of desiccation, therefore, did not refer to the whole plant, but to an angular portion of the canopy with an opening of 90° according to the 4 cardinal directions. Since the present work was intended to identify asymptomatic plants or those at a very early stage of the infection, only two classes were considered: Class 0, which included asymptomatic plants or those with a degree of desiccation in any angular portion of the canopy of no more than 5%; Class 1, which included all remaining plants. The choice of the threshold value was based on the experience of people well practiced in visual assessment.

2.3. Acquisition of Reflected Electromagnetic Radiation with UAV

The aerial survey was conducted with a very light quadcopter (DJI Mavic Pro drone), equipped with an on-board Global Navigation Satellite System (GNSS), and able to perform programmed missions over waypoints routes. The UAV had an unloaded endurance flight time of about 12/13 min and a signal maximum effective distance of 7 km. The dataset consisted of 1152 images captured on the site, at 70 m height, with a theoretical resulting Ground Sampling Distance (GSD) of 0.066 m/pixel. A custom payload tray was designed to carry the multispectral sensor Parrot Sequoia camera, which captures four narrowband spectral images in green (550 ± 20 nm), red (660 ± 20 nm), red-edge (735 ± 5 nm), and near-infrared (790 ± 20 nm) at a resolution of 1.2 Megapixel.

The choice of these bands, instead of the usual RGB, is due to the well-known utility of the red-edge reflectance as a general indicator of plant stress, and the NIR reflectance as an indicator of plant canopy structure and leaf area [21,22], which are the features of the olive tree most sensitive to Xfp infection. In fact, it is well known that the disease causes rapid desiccation of the canopy of the olive trees with a consequent macroscopic change in the coloring of the leaves and the shape of the crown.

This type of sensor, although limited in the number of bands and spectral resolution, represents nevertheless a good cost-effective solution, and is widely used in precision agriculture applications.

2.4. Selection of Plant Pixel

The semi-automatic extraction of canopy polygons from UAV data was achieved with the following steps (see Figure 2).

• Pre-processing UAV images

UAV surveys acquired a dataset consisting of 288 aerial shoots for each spectral band, for a total of 1152 images, with different geometric and spectral characteristics, which required the application of appropriate techniques and methods of analysis. The pre-processing step on data, including geometrical correction and calibration, was applied using version 4.4.12 of Pix4d Mapper [23]. After that, the software automatically managed the layers of the different bands by applying photogrammetry, using at least two overlapping images of the objects themselves, taken from different points of view to generate orthophoto images.

In order to have data with the same extent and pixel number for each band, some GCPs-Ground Control Points (target on the ground with well-known coordinates) were used by Pix4d Mapper. For the experimental site, 10 GCPs in cross-format (about 0.50 m × 0.50 m size) were distributed uniformly.

Starting from some points detected on multiple frames, a quoted point cloud was obtained. Then the software applied the SfT (Structure from Motion) technique to generate a final orthophoto image.

The further elaborations on images, to get out the delimitation of olive trees crowns, were carried out using the commercial ENVI (ENvironment for Visualizing Images) 5.1 software (Harris Geospatial) and the open source QGIS software (Quantum Geographic Information System) version 3.22, Bialowieza, Poland.

• Supervised classification

Supervised classification algorithm of Maximum Likelihood (ML) was applied to discriminate the crowns from background. The algorithm is a conventional classification technique that determines a known class of distributions as the one corresponding to the maximum of a given statistic. For the training samples, a normality hypothesis is made. ML determines for each class the probability density function, and each pixel is assigned to the class whose a posterior probability is maximized [24].

ML is recognized as a stable and robust classifier with high precision and accuracy [25]. As a pixel-based method, it is computer time-consuming method in classifying remotely sensed data [26,27].

In our case, only three labels were identified: soil, canopy, and shadow.

To apply this classification some prior knowledge of pixels is required. The user has to extract the training samples from the image that represent the prototype pixels assigned to the selected classes. To avoid transition areas between the classes with mixed pixels, it was necessary to select a sufficient number of training pixels (minimum 150 pixels per class) [28]. In this way, a properly representative sample of each class was obtained and the pixels were clustered according to their radiometric characteristics.

To modify the geometry features of the signal, a morphological filter, based on closing operator [29], was employed. This kind of operator edits the geometry features of the signal, deleting small holes inside the polygon [30,31,32].

The last step using ENVI consisted of the conversion from smoothed raster image into vector data format (ESRI polygon shapefile format), which was necessary for the successive geostatistical analysis.

• GIS environment

Each closed polygon was imported into QGIS environment. An editing procedure was also applied, operating on the attribute table of shapefile and selecting polygons, to improve the extraction of the crown of each olive tree. This was necessary when a polygon enclosed two crowns, and so had to be split. Therefore, the polygons were marked and used for all subsequent analyses.

2.5. Geostatistical Analysis

The data processing was carried out according to the standard formulation of multivariate geostatistics, as exhaustively described in several geostatistical manuals [18,19,33,34] and the interested reader can refer to them. Here, the only emphasis will be placed on those procedures most closely related to the identification of spatial scales, and to mapping of a multivariate spatial dataset using synthetic scale-dependent indices.

To calculate the experimental variograms of the UAV data, it was necessary to reduce the large number of pixels by superimposing a regular grid with 1 m mesh on the study area and choosing randomly one point falling within each cell contained within the shapefile of an olive tree canopy.

The spatial correlation structure of the multi-band UAV data was analyzed by adapting a linear model of co-regionalization (LMC) to the reflectance data of the four bands. LMC [33,35] considers all of the study variables to be generated by the same independent physical processes acting at N_s different spatial scales. Although this assumption is rather crude, for multi-band images from the same radiometer, it may be reasonable. All (both direct and cross-) variogram models are expressed as linear combinations of the same basic structures for each spatial scale (range), represented by authorized variogram models standardized to unit sill, and with the coefficients equal to partial sills [35]. In the direct variogram, the partial sill reflects the influence of the specific spatial scale on the total spatial variation of the study variable, whereas in the cross-variogram the partial sill expresses the influence of the specific spatial-scale on the covariance between the two involved variables. Therefore, in the latter case, partial sill is an indicator of the level of interaction between two variables. The set of partial sills of the direct and cross variograms for each identified scale constitutes a matrix known as the co-regionalization matrix for that scale.

After jointly fitting a nested variogram model to all (both direct and cross) variograms using the same basic spatial structures, the next step consists in applying the multivariate factor kriging [36]. This procedure of multivariate geostatistics is quite similar to the traditional Principal Component Analysis (PCA), from which it differs, in that, instead of applying to the whole variance-covariance matrix, it applies to its individual spatial components (i.e., co-regionalization matrices), each corresponding to a specific spatial scale. PCA decomposes each co-regionalization matrix into two additional matrices: the eigenvalue and eigenvector matrices [33] through the matrix of transformation coefficients:

$$B^u=Q^u{\mathsf{\Lambda }}^u{\left(Q^u\right)}^T=A^u{\left(A^u\right)}^T$$

(1)

where T is the transposed matrix symbol; $B^u$ coregionalization matrix; $Q^u$ matrix of eigenvectors; ${\mathsf{\Lambda }}^u$ matrix of eigenvalues; $A$ matrix of transformation coefficient; and $u$ the scale.

The transformation (loading) coefficients correspond to the covariances between the variables and the principal component, here called the regionalized factor, relative to a given spatial scale. They express the influence of each variable on the factor at a given spatial scale and are then used to assign a physical meaning to the factor. Furthermore, the eigenvalue associated with each factor is proportional to the variance explained by the factor at the specific spatial scale.

Therefore, by summing up all the eigenvalues corresponding to a given scale, the spatial component of variance relative to that scale can be obtained. This decomposition of the total spatial variance of a multivariate data set is of great utility because it allows the identification of the most relevant scales, conditional to the observation data.

Mapping the retained scale-dependent regionalized factors, which explain most variance at the specific scale, provides a very powerful and synthetic way to show the spatial relationships among several variables at the different spatial scales. Each map is a sort of snapshot of the study scene at the selected scale disregarding all others [37]. The approach outlined can be advantageously used in Precision Agriculture to produce a partition of the field in homogeneous zones at the specific spatial scales as derived from spatial analysis.

All geostatistical analyses were performed using the software ISATIS 2018 (Geovariances, France, 2018).

3. Results and Discussion

Figure 3 shows the images of the olive tree canopy reflectance at the different bands. From the visual inspection, an extreme within-plant variability is revealed, registered at the radiometer scale on board the drone of about 0.07 m. Specifically, the NIR and red-edge (Figure 3) maps show that the plants in the western sector exhibit less vegetative vigor than in the eastern sector. Notably, the plants in the northwest corner of the field appear particularly stressed. Most probably, the vector insect of the bacterium entered from this corner and then spread to the whole field, as was verified by the subsequent drone flights.

An LMC was jointly fitted to the ensemble of direct and cross-experimental variograms of the four variables including the following three basic spatial structures: a nugget effect, a cubic model with a range equal to 6.69 m, and a sincard model [38] ($\gamma \left(h\right)=C_0+C\left[1-\frac{\sin \left(\frac{\delta h}{a}\right)}{\frac{\delta h}{a}}\right]\delta =20.371$, where $h$ is lag, a range, $C$ partial sill, $C_0$ nugget effect, $\gamma$ variogram) with a range equal to 32.64 m (Figure 4). The fitting was quite good as tested in cross-validation by mean errors and standardized error variances, which were quite close to zero and one, respectively, for all variables [39]. The total spatial variance was therefore divided into three components: the one of spatially uncorrelated error accounting for 14%, and two spatially correlated components, the former at a shorter range (6.69 m) accounting for only 4%, and the latter at longer range (32.64 m) accounting for most of the total spatial variance (82%) (

Table 1. Co-regionalization Matrix for each spatial scale.

Nugget Effect
	Variable 1 *	Variable 2	Variable 3	Variable 4
Variable 1	0.0000	0.0002	0.0000	0.0001
Variable 2	0.0002	0.0012	0.0001	0.0008
Variable 3	0.0000	0.0001	0.0001	0.0001
Variable 4	0.0001	0.0008	0.0001	0.0005
Cubic − Range = 6.69 m
Variable 1	0.0001	0.0001	0.0001	0.0001
Variable 2	0.0001	0.0001	0.0001	0.0001
Variable 3	0.0001	0.0001	0.0001	0.0001
Variable 4	0.0001	0.0001	0.0001	0.0002
Sincard − Range = 32.64 m
Variable 1	0.0003	0.0012	0.0004	0.0009
Variable 2	0.0012	0.0064	0.0009	0.0044
Variable 3	0.0004	0.0009	0.0006	0.0008
Variable 4	0.0009	0.0044	0.0008	0.0032

* Variable 1 = green, Variable 2 = NIR, Variable 3 = red, Variable 4 = red-edge.

).

Table 1 shows the co-regionalization matrices for each spatial component (scale), representing the partial sills of the direct (on the main diagonal) and cross-(off-diagonal) variograms.

Omitting the nugget effect as it is most affected by errors; it is evident that at a small (intra-canopy) scale the contribution in terms of variance (main diagonal) is very low for all bands. In contrast, at a larger scale, the greatest contributions (main diagonal) can be attributed to the variances of the NIR and red-edge [38]. It is well known how the former can be taken as an indicator of canopy architecture and the latter as an integral indicator of the physiological status of the plant, which can be determined by multiple types of stress.

Experimental (dot) and model variograms in the two main directions: N–S (N0) and E–W (N90) (bold line) are shown in Figure 4. No clear anisotropy is evident except moderately for red and red-edge. The nugget effect plus the short-range spatial structure represents the variability within plants, as 6 m is the mean canopy diameter of these olive trees. The short-range component, including that at a scale smaller than the sampling scale (nugget effect), which was detected thanks to the high spatial resolution of the drone survey, actually accounts for a small proportion (18%) of the total variability associated with multispectral images, of which only 4% is spatially structured. This might be due to the extreme randomness of the variation of the foliage at this scale, attributable to both natural (intrinsic, genetic) properties of the olive tree and anthropogenic interventions such as pruning [9].

Conversely, the greatest proportion of total variability is found at a scale of around 32 m, comprising four plants on average, with a characteristic sinusoidal pattern (Figure 4), due to the discrete nature of the tree canopy. This range or scale of variation can be determined by various factors, such as, undoubtedly, the particular arrangement of the plants but also the dynamics of the spread of the infection, which was the focus of the research.

From the proximity of the cross-variogram model to the dotted curve (Figure 4), representing the maximum spatial correlation between the two variables (intrinsic correlation) [33], it is possible to assess the degree of spatial association between two variables. The significant spatial correlations of red-edge with all other variables, especially NIR, with the exception of red, are quite evident. This further justifies the assumption of red-edge as an integral measure of plant health. It is worth underlining also the strong spatial correlation between red and green, related to chlorophyll function, which was positive against expectation based on the standard reflectance spectrum (signature) of a green leaf [38]. This can be justified because the multi-band radiation reflected by a highly irregular olive tree canopy is the result of complex interactions between the radiations reflected by individual leaves.

Among the two sets consisting of four regionalized factors, each set corresponding to a given spatial scale (6.69 m and 32.64 m), only the first factor for each scale was retained because explaining a proportion equal to 96.79% of the variance at short range and 94.04% of the variance at longer range. Furthermore, their composition determined by the loading coefficients shows that on the first factor at the scale of 6.69 m all four variables weigh significantly although less the green. It can then be assumed as a general indicator of plant health because the red-edge reflectance weighs mainly and positively on it. Differently, the first factor at long range is mainly affected by NIR and red-edge, and can therefore be taken as a spatial index of canopy structure and plant vigor (

Table 2. Composition of regionalized factors at the different spatial scales.

Nugget Effect (Component Variance = 0.0018)
	Variable 1 *	Variable 2	Variable 3	Variable 4	Var. Perc. (%)
Factor 1	0.13	0.83	0.06	0.54	93.37
Factor 2	−0.39	0.27	−0.85	−0.23	5.42
Factor 3	−0.09	−0.49	−0.33	0.80	1.07
Factor 4	−0.91	0.05	0.41	0.09	0.14
Cubic − Range = 6.69 m (Component Variance = 0.0005)
Factor 1	0.32	0.48	0.47	0.67	96.79
Factor 2	0.23	−0.53	−0.51	0.63	2.72
Factor 3	−0.90	−0.10	0.17	0.38	0.49
Factor 4	0.17	−0.69	0.70	−0.08	0.00
Sincard − Range = 32.64 m (Component Variance = 0.0105)
Factor 1	0.16	0.80	0.14	0.56	94.04
Factor 2	0.39	−0.34	0.84	0.17	5.06
Factor 3	−0.05	−0.49	−0.34	0.80	0.81
Factor 4	0.90	−0.02	−0.41	−0.12	0.10

* Variable 1 = green, Variable 2 = NIR, Variable 3 = red, Variable 4 = red-edge.

).

To obtain an overview of the spatial association of the olive trees at the two identified scales from the UAV multi-band image, the first scale-dependent regionalized factors were interpolated on the grid of the drone images and shown in Figure 5.

The two factors show a quite different pattern: for the factor at the scale of 6.69 m (Figure 5a), the intra-canopy somewhat randomly distributed variability prevails, due to the previously mentioned irregularity in the shape of the olive crown, which makes it rather difficult to identify spatial regularities. The situation is not so with the factor at the scale of 32.64 m, because from its map (Figure 5b), it seems that in the eastern sector, the olive tree crowns appear more vigorous in their foliage, whereas the plants in the north-western corner are most stressed. Furthermore, at this scale, the field could be roughly divided into three zones: the southern zone in which most plants appear to be in good health with vigorous crowns; a central zone in which the plants may be suffering some type of stress, and finally the northern zone characterized by larger variability with a tendency towards less vigor in the western sector, which might be due to the desiccation produced by the infection.

At this point, it is worth emphasizing some important features of Xfp infection in the early stages that are punctually localized due to the action of the transmitter insect, so parts of the crown with different vigor may coexist in the same infected plant. This explains the extreme variability in canopy reflectance as detected by UAV images and the great difficulty in visually identifying infected plants, especially in the early stages of infection. The above reasons may explain why it was not possible to delineate clear and compact zones of different infection levels, also due to the discrete nature of the study area (set of spaced trees in an olive grove).

3.1. Potential of UAV Imaging

In order to assess the potential of UAV images to detect plant stress conditions in the primordial stages, the operator’s findings of the status of infection were superimposed on the map of the intra-canopy scale factor (black dots in Figure 6a). For the graphical display of the map, three classes of equal frequency were used to which three levels of vegetational vigor could be attributed: low, medium, high, or inversely three degrees of stress: high, medium, and low. It should be underlined that the two types of observations not only differ profoundly in their nature but also refer to a different measurement support [40]. The spatial resolution of the UAV image is a pixel of 4.9 × 10⁻³ m², while the support of the visual survey by sector varies between 3 and 7 m² assuming a crown radius varying between 2 and 3 m, with differences therefore of about three orders of magnitude. The drone image can therefore detect very fine differences in plant vigor.

Examination of Figure 6a shows that only six plants were classified as infected by the operator in the northwest and southeast corners. However, based on the interpretation given to the short-range factor as an index of plant health, it can be hypothesized that other plants, indeed the majority, may already be suffering from some type of stress. Clearly, the limitations previously highlighted in Figure 5 also apply to Figure 6.

Since, from UAV images, it is not possible to discriminate plant desiccation due to Xfp infection from that attributable to, e.g., water stress, those plants with a greater extent of low level of the factor should be marked. Destructive samples of leaves would be taken and brought to the laboratory for direct analysis aimed at detecting the presence of the bacterium [8,41].

To verify the plausibility of our hypotheses in Figure 7, the maps of red-edge, assumed as an integral index of the physiological state of the plant, at the four dates (20 September 2017 (a), 9 March 2018 (b), 20 June 2018 (c), and 2 August 2019 (d)) of the successive drone flights are shown [9]. To represent the spatial variability of the red-edge, two types of scales were used: (1) A numerical one, reporting the reflectance values in the red edge, using not a linear but an isofrequency scale. The entire range of variation was divided into three classes of equal numerosity (33%), each representing a different level of reflectance. (2) A qualitative one. Given the general direct relationship between red-edge and canopy leaf vigor, each of the above classes was assigned a qualitative value, starting from the lowest to the highest value of reflectance, of low, medium, and high plant vigor. This isofrequency scale was preferred because it enhances the differences between classes and was used also for Figure 6. Obviously, a plant in the low class is not necessarily infected, but it is very likely that it is suffering some kind of stress. For the reasons expressed earlier, the same numerical scale was not used for all the maps as the aim was not to compare them, but to highlight the plants under stress over time. It is evident how the stress localized practically in the northwest corner is spreading, within a few months from September to March, to all plants, albeit with a different level of intensity and extension and preferably in the northwest sector direction.

The photos shown in Figure 8 illustrate the rapid desiccation produced by Xfp infection in a plant, classified as infected on the date of the first flight, over time (photographic material provided by the operator). Details on the temporal evolution of the infection can be found in [9].

3.2. Differences in the Use of the Maps at the Different Scales

Regarding the different use of the two maps related to different scales, for the short-range one (Figure 6a), its application in pest control has been discussed extensively in the previous sections. As far as the longer scale map (32.64 m) is concerned (Figure 6b), having filtered out the variance associated with the error and that related to the intra-canopy scale, it better highlights structural variations in the olive grove. These could be due to both intrinsic characteristics linked to the cultivar and the plant arrangement as well as the prevailing directions of infection spread resulting from the insect’s movements. At this scale of variation, it has already been mentioned above how the olive grove could roughly be divided into three zones of different vegetative vigor, with the central one clearly more stressed. Such a map could support agronomic operations working on a larger scale, such as irrigation, fertilization, or mechanical tillage. It should be emphasized that for effective decision support, UAV images alone would not be sufficient, but should be supplemented with remote or proximal sensing sensors operating in the electromagnetic range of SWIR (1000–2500 nm) [41,42] and/or IR (12–16 µ), the latter, especially for precision irrigation [43].

Moreover, at this scale, it is quite evident that the plants are less vigorous in the western than in the eastern sector. The entry of the vector insect probably took place in the northwest corner and the direction of spread was probably northwest, as confirmed by the successive UAV images (Figure 6). A further use of this map could therefore be to mark plants with different risks of infection [8] and consequently with different control priorities.

4. Conclusions

In this paper, the application of an approach based on multivariate geostatistics is illustrated, which allows the identification of the scales of spatial variation observed in a high spatial resolution multi-band image from a drone flying over an olive grove affected by Xfp infection. The method is sufficiently flexible and can also be applied to other RS types (from aircraft or satellite) with any wavelength range, and spectral and spatial resolutions. The use of a validated nested variogram model is of fundamental importance here. It also makes it possible to relate the variation observed in the data to both the scales of the phenomenon under study and of the potential agronomic applications. Multiple factor kriging has provided two different views of the same scene but at two distinct spatial scales. The maps obtained clearly show how the spatial patterns displayed are a function of scale. The variogram is therefore of fundamental importance when making decisions based on the most relevant scale for the phenomenon under study and/or for possible operational applications.

Although there are also other techniques for analyzing UAV data in the fight against infection, such as modern artificial intelligence and machine learning methods [44], we believe that these have two limitations. The first is related to the fact that to be efficient and accurate, they need to be used on large data sets (big data), and the second is that they are based on the assumption of spatial independence of the data [45]. This limitation is particularly critical in the case of the use of UAV data for pest control because of the very fine spatial (centimeter) resolution, whereby a high spatial correlation is expected. Moreover, as one of the objectives of the study was also to discover any directions of higher spatial correlation, where infection was most likely to spread, we believe that geostatistics can provide answers to this.

The proposed approach should have a broader appeal, not limited to precision agriculture alone, as it can also be used to describe the appropriate scales for modeling other environmental phenomena, for example for controlling water resources in land management or describing animal and plant habitat relationships in ecology.

Although this work has underlined the potential of using UAV images in plant infection control, it is still a long way from automatically identifying the disease in the plant, so further research is needed. Several points that are relevant to UAV-acquired images should be addressed. Isolating plants from the background (mostly soil) can be a difficult problem depending on the spatial resolution of the images. Mixed pixels (plant + soil, plant + shadows) will inevitably be present reducing the accuracy of segmentation. Add to this the presence of weeds, which complicates the delineation of regions of interest and then affects the correct detection and quantification of infection.

Most plant infections produce physiological changes that can be detected in only a few specific bands, which depend not only on the infection but also on the cultivar, the presence of other stresses, and environmental conditions. In this regard, the use of hyperspectral radiometers with high spectral resolution is claimed over broadband radiometers such as the one used in this work.

In practice, there may be many factors, internal and external to the plant, which produce in it a response like that produced by the infection. In addition, various stresses may be simultaneously present. This leads to the assertion that it is currently too limited to base an alarm solely on abnormal coloring of the drone image and that scouting in the field is always necessary to define the actual cause of such an abnormality.

However, it is evident that, while using different approaches and sensors, combining the information contained in a UAV image at different scales improves stress detection, monitoring capability, and efficacy of application.