Prediction of the Kiwifruit Decline Syndrome in Diseased Orchards by Remote Sensing

Eight years after the first record in Italy, Kiwifruit Decline (KD), a destructive disease causing root rot, has already affected more than 25% of the area under kiwifruit cultivation in Italy. Diseased plants are characterised by severe decay of the fine roots and sudden wilting of the canopy, which is only visible after the season’s first period of heat (July–August). The swiftness of symptom appearance prevents correct timing and positioning for sampling of the disease, and is therefore a barrier to aetiological studies. The aim of this study is to test the feasibility of thermal and multispectral imaging for the detection of KD using an unsupervised classifier. Thus, RGB, multispectral and thermal data from a kiwifruit orchard, with healthy and diseased plants, were acquired simultaneously during two consecutive growing seasons (2017–2018) using an Unmanned Aerial Vehicle (UAV) platform. Data reduction was applied to the clipped areas of the multispectral and thermal data from the 2017 survey. Reduced data were then classified with two unsupervised algorithms, a K-means and a hierarchical method. The plant vigour (canopy size and presence/absence of wilted leaves) and the health shifts exhibited by asymptomatic plants between 2017 and 2018 were evaluated from RGB data via expert assessment and used as the ground truth for cluster interpretation. Multispectral data showed a high correlation with plant vigour, while temperature data demonstrated a good potential use in predicting health shifts, especially in highly vigorous plants that were asymptomatic in 2017 and became symptomatic in 2018. The accuracy of plant vigour assessment was above 73% when using multispectral data, while clustering of the temperature data allowed the prediction of disease outbreak one year in advance, with an accuracy of 71%. Based on our results, the unsupervised clustering of remote sensing data could be a reliable tool for the identification of sampling areas, and can greatly improve aetiological studies of this new disease in kiwifruit.


Introduction
Italy is, beside New Zealand, one of the two major kiwifruit producers and exporters worldwide, accounting for 15% of the total world demand on its own [1]. Since 2012, a major new disease known as Kiwifruit Decline (KD) has been threatening the future of kiwifruit production in Italy [2]. Up to now, KD has already destroyed more than 25% (~6600 ha) of the Italian kiwifruit growing area, and it has become a serious threat to the future of the crop in this country (Tacconi, personal communication which would be otherwise undetectable by the naked eye. Since leaves emit radiation according to their temperature in the thermal infrared band (TIR; 3 to 15 µm) [23,30], this wavelength range has been widely used to estimate the water status of several plants [31][32][33], and to detect a number of plant diseases that influence leaf stomatal conductance and transpiration rates directly or indirectly [34][35][36]. Regarding kiwifruit, thermal imaging has been previously adopted for the early detection of another vascular disease caused by the bacterium Pseudomonas syringae pv. actinidiae [37].
Among the remote sensing technologies, drone-based surveys provide several advantages in KD studies: (i) they allow us to generate high resolution images that can be registered to image the whole field and quickly evaluate the health status of the plants; (ii) they can be easily equipped with different sensors, collecting data in non-visible wavelengths such as the near infrared region (NIR) and thermal infrared; (iii) they can fly at lower altitudes, reducing the influence of the atmosphere; (iv) they are cost efficient; and (v) they have a high scheduling flexibility, which is crucial for surveying plant diseases that are strictly related to climatic condition [38,39]. Considering KD, the latter aspect is probably the most important, since the survey must be performed not only with favourable weather conditions (radiation, clouds, wind), but also with optimal soil water content, since both waterlogging and drought stress can affect the transpiration rates [26,27,29].
Images acquired with these systems can be classified using either supervised or unsupervised machine learning algorithms. However, the application of supervised algorithms is limited, since it is difficult to produce sufficient high-quality labelled data to create a robust supervised classifier that is able to account for all the variability associated with KD and kiwifruit orchards. The best method for labelling the data is the observation of the disease spread between two consecutive growing seasons, since protocols to rapidly and objectively detect diseased plant are still missing [2]. However, the swiftness of the disease spread make the use of supervised algorithms impractical. Indeed, by the time the classifier has been trained, the field could have already been irreversibly compromised. Unsupervised learning may prove to be an effective alternative, since it can reveal the underlying patterns and structure of the data without the need for labelling [40].
Thus, this study aims to understand the feasibility and reliability of thermal and multispectral imaging in predicting and assessing disease outbreaks, via an unsupervised classification approach. To the best of our knowledge, the present study is the first one that introduces remote sensing as a reliable tool for addressing aetiological and epidemiological studies of Kiwifruit Decline.

Experimental Site
UAV-based surveys were performed over one hectare of Actinidia deliciosa (A. Chev.) CF Liang et AR Ferguson cultivar 'Hayward' in San Giorgio della Richinvelda (Friuli Venezia Giulia, NE Italy, 46.039 • N 12.864 • E). The orchard was 20 years old, grown in a 5 × 5-m T-bar plantation without anti-hail netting, with a micro-irrigation system (micro-sprinkler, 20 mm/hour) and with the following soil properties: 8% gravel, 36% sand, 59% silt, 5% clay and 1.04% organic matter (w/w). The first symptoms of KD were recorded in 2015 after an irrigation pipe broke in the central rows of the orchard. Later in 2016, the disease rapidly spread towards the border of the field, compromising approximately 30% of the field. From 2017 the field was regularly inspected to check the pattern of disease spread. In that year, the disease affected almost the entire orchard, compromising approximately 70% of the field, and by 2018 only a few plant rows on the eastern boundary were still apparently healthy. At the end of 2018, the field was definitively explanted. From the first appearance of the disease, the orchard was irrigated by following the usual agronomic practices and avoiding both drought stress and excessive irrigation.

Image Acquisition
The UAV surveys were carried out on 8 August 2017 and 18 July 2018, respectively at 11.00 a.m. and 10.00 a.m. on hot sunny days. The timing of the flight was based on a previous report suggesting Remote Sens. 2020, 12, 2194 4 of 19 this time frame as one that would maximise the differences between well-fed and drought-stressed plants [28]. Weather condition were favourable for the acquisition: no clouds, wind speed <2 km/h, relative air humidity~44%, radiation~2900 MJ/m 2 and temperature~29 • C. The acquisitions were performed with a single flight lasting less than 20 min, including take-off and landing. The flights were undertaken, respectively, two days after an irrigation event in 2017 and three days after a heavy rainfall (20 mm) in 2018, to exclude water deficit from the variables influencing the canopy temperature and to make sure that the soil water content was not a limiting factor. The gravimetric soil moisture was measured the day before the flight at nine locations in the surveyed area. Three sampling areas were selected inside an apparently healthy area, three in a heavily compromised area and three in a transition area between these two. The surveys were performed if the soil humidity was between 70% and 80% of the field capacity, representing an optimal condition for kiwifruit growth [41]. The field capacity value was derived from the soil texture using soil-plant-atmosphere-water field and pond hydrology (SPAW) models [42]. SPAW models are pedo-transfer functions, which use soil texture (percentage of sand, silt and cay) and percentages of soil organic matter to estimate the soil hydrodynamics properties relevant to agronomical practices, such as volumetric field capacity and wilting point [43].
The UAV flights were performed with a hexa-copter using Real Time Kinematic-Global Navigation Satellite System (RTK-GNSS) and equipped with a gimbal system on two axes (Adorn-technologies, Italy). Thermal, Multispectral and RGB images were acquired simultaneously at a flying height of 35 m, at a flying speed of 5 m/s. Frame rate was adjusted for each sensor in order to achieve a forward image overlap of 80%, while sidelap overlap was 80%, 88% and 90% between thermal, RGB and multispectral images, respectively. In 2017, the acquisition combined the three sensors: GoPro Hero 4 (GoPro, San Mateo, CA, USA) as the RGB camera, Tetracam ADC snap (Tetracam, Chatsworth, LA, USA) as multispectral sensor and FLIR TAU2 640 for thermal imaging (FLIR System, Wilsonville, OR, USA). In 2018 the RGB and multispectral sensors were exchanged for Sequoia+ sensors (Parrot, Paris, France). For detailed information about the sensors used see Table 1. Tetracam ADC snap is a low-cost multispectral sensor, modified from an RGB sensor equipped with a blue absorbing glass filter to eliminate the blue sensitivity, and to use the blue pixels in the sensor to measure NIR bands. FLIR TAU2 640 is a thermal camera capable of detecting temperature differences of ±0.05 • C with a sensibility of 2 • C. It has the advantage of performing automatic flat-field corrections (FFC) while recording, to account for microbolometer variation and remove artefacts from the 2D images. For this study, the FFC was performed every 100 frames. FLIR TAU2 640 also measures the internal temperature of the camera to frequently update the non-uniformity coefficients used to convert raw data into radiometric values [44]. Air temperature and relative humidity were measured before the flight and inserted in the camera settings, together with the plant emissivity value, which was set to 0.98 as suggested in Maes et al. (2014) [37]. Sequoia+ is a camera equipped with a sunshine sensor for image calibration, an RGB camera and four single-band cameras measuring reflectance in the green, red, red-edge and near-infrared (NIR) bands.

Pre-Processing: Creation of Orthomosaics and Alignment Checking
An overall visualisation of the workflow is presented in Figure 1. The analysis was performed using all data (RGB, multispectral and thermal) acquired in 2017, while in 2018, only RGB data were used to evaluate the disease spread, as detailed below. After the acquisition, raw RGB and multispectral data were converted to reflectance and thus normalised from 0 to 1 using Pix4D software (Pix4D., Lausanne, Switzerland). In 2017, calibration was performed using a white Teflon plate included in the Tetracam ADC snap packaging as white reference, and as dark reference a picture acquired with the shutter closed. White and dark references were acquired at the beginning of the flight. Raw measurements were then calibrated as reported in Elmasry et al. (2012) [45]. In 2018, RGB images were instead normalised following the manufacture-recommended methods for the Sequoia+ sensor, using the one-point calibration plus sunshine sensor method as reported in Poncet et al. (2019) [46]. Thermal data was not calibrated; however, rather than absolute leaf temperature, it was of higher Remote Sens. 2020, 12, 2194 5 of 19 interest in the relative difference occurring between plants with different health statuses. Temperature values ( • C) were derived from the TIR wavelengths using FLIR-studio (FLIR System Inc., USA). Stitching, georeferencing and orthorectification were performed independently for each set of data (RGB, multispectral and thermal data from 2017 and RGB from 2018) with Pix4D. Image processing is based on structure from motion (SfM) algorithms to reconstruct the three-dimensional scene based on shared features detected across the images. The four orthomosaics (multispectral and thermal from 2017; RGB from 2017 and 2018) were combined in Qgis [47], and planting lines and control points were drawn separately over each orthomosaic to check alignment of the datasets. A total of 11 control points were selected as follows: four points in both the northern and southern part of the field, and three in the orchard. Errors between planting line slope were minimal (on average 0.01 • ), while the mean error between control points was 4.97 cm horizontally and 4.51 cm vertically, suggesting that no major deformation occurred during the stitching. Alignment was deemed good enough for the analysis purposes, therefore co-registration of the orthomosaics was not performed. Errors linked to temperature drifts were checked at two pairs of points: two points diagonally collected on a grass-free compacted portion of the earth roads delimiting the field on the northeast and southwest sides; and two points were collected on the water contained in an irrigation canal which flanks the north side of the field. Temperature values for each point were extracted from the thermal orthomosaic and respectively compared to check for temperature drifts, which were 0.2 • C on the bare soil and 0.04 • C on the water.

Post-Processing: Clustering and Reference Data
Single-band orthomosaics were extracted from the multispectral orthomosaic and tested for band correlation. Green and red bands were de-correlated from the NIR band using the formula R adj = R ini − R NIR * 0.8, where R adj is the corrected reflectance value, R ini is the initial value of the red or green bands and R NIR is the reflectance value in the NIR spectrum. This correction reduced the R 2 correlation from 0.89 to 0.64 between the red and NIR band, and from 0.92 to 0.81 between the green and NIR band. Temperature and RGB values were not modified after the stitching. Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 19 Figure 1. Overall processing pipeline. In (a) processes are represented with a darker background, while data have a light background. The green steps represent data and processes, used for the unsupervised classification of multispectral and thermal data; the blue steps represent factors used for reference data estimation based on RGB clips; the yellow steps are common between the two processing pipelines. In (b) details of the segmentation process are depicted: starting with mask positioning along planting lines, proceeding with rectangular shape clipping and concluding with circular mask extraction. In (c) the classes derived by the expert assessment are shown for each vigour and plant health shift class. On the left, the four vigour classes are represented from highly vigorous and asymptomatic (V4) to dead (V1). On the right, the classification of the health shifts that occurred between 2017 and 2018 in plants that were initially highly vigorous and asymptomatic (V4) is reported: from those that remain highly vigorous and asymptomatic in the second year (S1), to the ones that in 2018 were weakened (S2), diseased (S3) or dead (S4). Abbreviations: ROI, region of interest; GeoRef, georeferenced.

Post-Processing: Clustering and Reference Data
Single-band orthomosaics were extracted from the multispectral orthomosaic and tested for band correlation. Green and red bands were de-correlated from the NIR band using the formula = − * 0.8, where is the corrected reflectance value, is the initial value of the red or green bands and is the reflectance value in the NIR spectrum. This correction reduced the R 2 correlation from 0.89 to 0.64 between the red and NIR band, and from 0.92 to 0.81 between the green and NIR band. Temperature and RGB values were not modified after the stitching.
Identification of the stump positions and determination of the canopy boundaries is very difficult to achieve when performing in-field UAV surveys over kiwi vines because, firstly, T-bar or pergola training systems create a flat and dense canopy preventing stump identification, and secondly, the canopies of two neighbouring plants usually overlap due to the disordered growth of shoots. For these reasons, a rough segmentation algorithm was applied to the RGB, thermal, and each single-band multispectral orthomosaic. Firstly, the orthomosaics were clipped along the planting Figure 1. Overall processing pipeline. In (a) processes are represented with a darker background, while data have a light background. The green steps represent data and processes, used for the unsupervised classification of multispectral and thermal data; the blue steps represent factors used for reference data estimation based on RGB clips; the yellow steps are common between the two processing pipelines. In (b) details of the segmentation process are depicted: starting with mask positioning along planting lines, proceeding with rectangular shape clipping and concluding with circular mask extraction. In (c) the classes derived by the expert assessment are shown for each vigour and plant health shift class. On the left, the four vigour classes are represented from highly vigorous and asymptomatic (V4) to dead (V1). On the right, the classification of the health shifts that occurred between 2017 and 2018 in plants that were initially highly vigorous and asymptomatic (V4) is reported: from those that remain highly vigorous and asymptomatic in the second year (S1), to the ones that in 2018 were weakened (S2), diseased (S3) or dead (S4). Abbreviations: ROI, region of interest; GeoRef, georeferenced.
Identification of the stump positions and determination of the canopy boundaries is very difficult to achieve when performing in-field UAV surveys over kiwi vines because, firstly, T-bar or pergola training systems create a flat and dense canopy preventing stump identification, and secondly, the canopies of two neighbouring plants usually overlap due to the disordered growth of shoots. For these reasons, a rough segmentation algorithm was applied to the RGB, thermal, and each single-band multispectral orthomosaic. Firstly, the orthomosaics were clipped along the planting lines with rectangular masks. These digital masks were 5 m wide and 1.5 m long and were placed with the longer side perpendicular to the direction of planting (Figure 1b). The rectangles were placed 1.5 m apart from each other to avoid overlapping. Secondly, a circular shape with a radius of 1 m was applied to the centre of each rectangle clip, retaining only the area within the circle. The shape resulting from the rectangle's and circle's intersection was selected to reduce the border effect, while the mask dimensions and the distances between the mask centres were set to minimise the interference from background pixels of the inter-row grass, and to reduce the canopy overlap within the same clip. Segmentation of RGB data was stopped at the rectangular shape clipping, while multispectral and thermal data were subjected to the complete segmentation workflow. The RGB clips were used to estimate reference values for plant vigour and plant health shifts via expert assessment, while multispectral and thermal data were used for unsupervised clustering. Maintaining the entire rectangular clips for RGB data improved expert awareness of the canopy dimensions, and consequently the precision of plant vigour classification. On the other hand, multispectral and thermal data was reduced to the assumed canopy extent to reduce background bias, since the central portion of the canopy is usually flat and orthogonal to the acquisition of the images (Figure 1b).
Two types of reference data were evaluated from the rectangular RGB clips via expert assessment: plant vigour and plant health shifts (Figure 1c). Plant vigour is a good indicator to distinguish asymptomatic and symptomatic plants, but it is not directly associated with the status of the root system [2]. Therefore, the a posteriori evaluations of vigour shifts and the death rates between years were the only available means to perform realistic evaluations of plant health status in 2017 without uprooting plants.
Reference data for canopy vigour was estimated based on RGB clips derived from the segmentation process of the 2017 orthomosaics. Each clip was evaluated individually and independently by two experts and ranked in four classes based on canopy vigour, considering: (i) the presence/absence of symptoms on the leaves (yellowing, scorching or wilting), (ii) disease severity (percentage of canopy with wilting symptoms), and (iii) fractional vegetation coverage (Table 2 and Figure 1c). Reference data for plant health shifts were derived by comparing plant vigour classes in 2017 with those observed in 2018. Vigour classes for 2018 were evaluated via the same procedure described above. Shifts that occurred in highly vigorous plants from 2017 (V4) were the only ones considered due to their relevance for aetiological studies (Table 2 and Figure 1c). Indeed, plants with canopies resembling the diseased classes (V2 and V1) usually possessed a heavily deteriorated root system that may have already been colonised by secondary invaders, while the vigour reduction of plants in class V3 could have been caused by several aspects unrelated to the disease being studied.
Data reduction was applied to multispectral and thermal data by averaging the pixel values within each masked area and storing the results in a matrix associated with the coordinates of the masking centres. K-means [48] and Ward's hierarchical [49] clustering were then performed separately for multispectral (NIR and red bands) and thermal data, using two, three and four clusters. The highest number of clusters was selected based on the capability of the experts to correctly classify the reference Remote Sens. 2020, 12, 2194 8 of 19 data, while lower numbers were tested to better understand the clustering prediction capability (see below for details).

Cluster Interpretation
Finally, clustering results were compared with the reference classes of plant vigour in 2017, and the plant health shifts between 2017 and 2018, via a two-step analysis: firstly, mosaic plots were used to visualise and explore possible correlations between clusters and reference data; secondly, predictive statistics were evaluated based on confusion matrices considering the result of the correlation analysis.
Mosaic plots can display the relationship between categorical variables using rectangles whose areas represent the proportion of cases for any given combination of the multivariate categorical data [50]. Paired with residual analysis, such as Pearson Xˆ2, the significance of such a correlation can be estimated [51]. A cluster can be considered representative of a reference class if it is positively correlated with only one class and negatively correlated with all the others.
Predictive statistics (sensitivity, accuracy, precision and F1 score) were then evaluated for combinations (reference x cluster) whose correlation was significant (p > 0.05) and meaningful ( Figure 4). Therefore, vigour classes were compared to the clusters obtained from multispectral data, while thermal data were used for disease outbreak prediction. If the number of clusters was inferior to the number of reference classes, the latter was reduced to match the former, aggregating classes of the reference data ( Figure 4). Reclassification was performed by merging the neighbour classes reported in Table 2 only if they fit the following criteria: (i) they were correlated (or at least highly represented) with the same cluster, (ii) the overall accuracy of the confusion matrix increased after class aggregation, and (iii) the new aggregated classes preserved an internal logic. Therefore, the four reference classes for plant vigour (Table 2 and Figure 4) were reduced to three and two by respectively merging the V3 and V2 classes (weakened plants) and V4, V3 and V2 classes (plants with leaves). The new aggregated reference data for vigour were then compared via multispectral clustering using three or two clusters respectively (Figure 4).
The same criteria adopted for the aggregation of vigour classes were applied also for the plant health shifts, which were compared with two clusters derived from temperature data. The plants that remained asymptomatic without changing their vigour in 2018 (S1) were merged with those displaying slightly reduced vigour without showing symptoms on the leaves (S2), creating a new class characterised by plants that remained unharmed, whereas the plants that died in 2018 (S4) or were heavily compromised (S3) were grouped in a new class composed of plants that most likely were already diseased in 2018 ( Figure 4).
All the analyses regarding decorrelation, segmentation, clustering and statistics of the clusters were performed with RStudio [52].

Multispectral and Thermal Data as an Indicator of Plant Vigour and Disease Outbreaks
For both methods tested, multispectral data were mostly correlated with plant vigour (Figure 2a-f), but it was found to not be a good predictor for disease outbreaks (Figure 2g-l). Indeed, the spatial distribution of multispectral clusters resembles that of the reference data for vigour classes (Figures 3c-h  and 3a, respectively). With four clusters, it was possible to estimate all the vigour classes needed for the evaluation of the disease. With three, the weakened plants (V3 and V2) were separated from the asymptomatic and the dead ones, while with two clusters, plants with leaves were distinguished from those already dead. Overall, K-means performed better than the hierarchical method, since the differences between the middle vigour classes V3 (weakened) and V2 (diseased) could be better identified. However, hierarchical clustering showed higher correlations towards the highly vigorous plants (V4). (v-x) report those from the hierarchical method. White represents non-significant correlation. Positive and negative correlations are highlighted in blue and red, respectively, and shaded in relation to their confidence level (Pearson test): 95% confidence in brighter colours, 99% in darker tones. The bases of the rectangles are proportional to the total number of images within each reference class (numbers at the bottom of the graph), whereas heights are the proportion of reference data within the specific cluster.  centres. In (b) and (i-n) white represents plants that were not highly vigorous in 2017 and were therefore removed from the plant health shift analysis. Yellow rectangles in (b) and (f-h) identify the hottest area within the orchard that happened to have the highest disease incidence in 2018.
When K-means was applied to four clusters to classify multispectral orthomosaics (Figure 2a), each cluster was positively correlated with only one vigour class, and negatively correlated with all the others. Almost all the correlations were significant or highly significant, suggesting that cluster KM4 can be used for the detection of apparently asymptomatic plants with high vigour (V4), cluster KM3 for asymptomatic plants with low vigour (V3), cluster KM2 for plants that are heavily compromised but still alive (V2), and KM1 for the dead plants (V1). Misclassification errors occurred mostly between the narrower cluster classes KM3 and KM4. Hierarchical clustering showed a higher correlation between the V4 class and the HM4 cluster, and performed similarly regarding the HM1-V1 and the HM2-V2 associations (Figure 2d). However, the cluster HM3 was not clearly correlated with the V3 class, since it also contained some diseased plants belonging to the V2 class.
Using K-means with three clusters on the multispectral data (Figure 2b), the extreme vigour classes V4 and V1 were still well discriminated by KM3 and KM1, respectively, while cluster KM2 was mostly correlated with classes V2 and V3, corresponding to plants with a weakened status. Considering this association, errors can occur when assuming cluster KM3 to be solely associated with the highly vigorous plants in V4. By using hierarchical clustering, the association HM3-V4 and HM2 (V2 and V3) was maintained (Figure 2e), but most of the misclassifications occurred within the cluster HM1, which was not highly correlated solely with dead plants (V1) and was also correlated with the V2 class.
Finally, with two clusters only, both methods performed equally (Figure 2c,f). The plants with an asymptomatic canopy (V4 and V3) were clearly separated from the dead plants (V1), and highly correlated with clusters KM2 or HM2, and KM1 or HM1, respectively. Major errors occurred in the classification of plants that were diseased but still alive (V3), which seemed to be better correlated with clusters KM1 and HM1.
For both methods tested, thermal data were mostly correlated with health shifts (Figure 2s-v), but could not be associated with any vigour classes (Figure 2m-r). The hottest area in 2017 within the orchard happened to be the one with the highest disease incidence in 2018 (Figure 3b,j-n). Canopy temperature ranged between 22 • C and 30 • C, with the clusters averaging between 23 • C and 26 • C, depending on the number of clusters used. K-means centres were on average slightly higher (+0.3 • C) than the means evaluated for each hierarchical cluster. The mosaic plots suggested that the temperature data were able to predict disease outbreaks by clustering the plants into two groups: one where disease symptoms appeared (S3 and S4) and another where the plants remained asymptomatic (S1 and S2). Indeed, among the tested combinations, the two-cluster model showed the best discrimination between the plant health shift classes. Both algorithms performed similarly, with K-means manifesting slightly better correlations between the plants that remained highly vigorous in 2018 (S1).
Nevertheless, the plants that died in 2018 (S4) were always positively associated with the hottest clusters, regardless of the number of clusters used (Figure 2s-v). Furthermore, plants that were heavily compromised, but still alive, in 2018 (S3) were usually abundant in the hottest cluster, although not always with a significant correlation for all numbers of clusters tested (Figure 2s-v). Lastly, the coldest clusters, KT1 and HT1, were consistently positively associated with plants that also maintained high vigour in 2018 (S1), and negatively correlated with plants that became diseased (S3 and S4) (Figure 2s-v).

Clustering Interpretation
Based on the mosaic plot results, the accuracy in estimating the vigour classes or plant health shifts was tested via a confusion matrix (Figure 4). The clustering of multispectral data performed similarly for both methods, with better statistics for K-means compared to hierarchical clustering. Dead plants (V1) were the class with the highest predictive statistics for all the numbers of clusters tested in both methods, obtaining the best score set when four clusters were used. By reducing the number of clusters, the precision dropped to a minimum of 60%, but the sensitivity increased greatly (98%). The same behaviour was observed for plants with high vigour (V4) that were correctly classified into clusters KM4 or HM4 and KM3 or HM3, when the data were grouped into four or three clusters, respectively. between weakened plants (V3+V2) and the middle cluster HM2 still existed (precision 81%), but the accuracy was reduced by the association of some plants of class V3 (diseased) with the HM1 cluster that was highly correlated with dead plants as well (V1) (Figures 2e and Figure 4f). Finally, the predictive statistics obtained with two clusters were almost equal to those obtained by the K-means method (Figure 4g).
The predictive capability of thermal data was tested for health shifts, which occurred in plants that were highly vigorous in 2017 (Figure 4d,h). The best predictive capability for disease spread was observed when clustering was performed with only two clusters, which allowed discrimination between plants that showed wilting and those that did not. The performances of the two methods were similar, with K-means being slightly better than hierarchical clustering (Figure 4d,h). Figure 4. Confusion matrices of unsupervised clustered data for the assessment of disease spreading (plant vigour) and the prediction of future outbreak (plant health shifts). K-means (a-c) and hierarchical clustering (e-g) performance in assessing disease spreading using 4, 3 or 2 clusters, and using multispectral data. (d) and (e), performance of thermal clusters in predicting future outbreaks using 2 clusters, for K-means and hierarchical clustering, respectively. In (b-d) and (f-h) reference data were merged if the number of clusters was inferior to the number of reference classes. In (b) and (f), a class was created aggregating weakened and diseased plants (V3 and V2). In (c) and (g), the three classes with leaves (V4, V3 and V2) were grouped together. In (d) and (h), plants that remained asymptomatic until 2018 (left) were split apart from those that showed wilting symptoms in 2018 (right). Abbreviation: Sesn, Sensitivity; Prec, Precision; F1, F1-score; Acc, accuracy; KM1-KM4 and KT1, KT2 multispectral and thermal clusters.
For the K-means method, using four clusters was the approach that showed the best association with plant vigour. Each cluster was associated with a single class with an accuracy above 73% (Figure 4a). In particular, the extreme classes (V1 and V4) were precisely identified (precision above 78%), with a low number of false positives. Misclassification errors mostly occurred between neighbour classes, mostly due to false positives that occurred in the association of cluster KM3 with the vigour class V3. Reducing the number of clusters to three, the central cluster KM2 was associated with weakened plants (V3 and V2) with an accuracy of 73% and a precision above 82% (Figure 4b). Finally, using only two clusters, it was possible to distinguish between dead plants with no canopy and plants with leaves (Figure 4c).
Hierarchical clustering was more sensible than K-means clustering for highly vigorous plants (V4), especially when four clusters were used (Figure 4e). The classification of the middle classes V3 and V2 was the major reason for errors. Indeed, with four clusters, this caused a reduction of the accuracy by as much as 0.68% and 0.67%, when V3 and V2 classes were respectively associated with HM3 and HM2 (Figure 4e). By reducing the number of clusters to three (Figure 4f), the association between weakened plants (V3+V2) and the middle cluster HM2 still existed (precision 81%), but the accuracy was reduced by the association of some plants of class V3 (diseased) with the HM1 cluster that was highly correlated with dead plants as well (V1) (Figures 2e and 4f). Finally, the predictive statistics obtained with two clusters were almost equal to those obtained by the K-means method (Figure 4g). The predictive capability of thermal data was tested for health shifts, which occurred in plants that were highly vigorous in 2017 (Figure 4d,h). The best predictive capability for disease spread was observed when clustering was performed with only two clusters, which allowed discrimination between plants that showed wilting and those that did not. The performances of the two methods were similar, with K-means being slightly better than hierarchical clustering (Figure 4d,h).

Discussion
Our results suggest that UAVs combined with multispectral and thermal imaging can be useful tools for scouting activities at a field scale in kiwifruit orchards, and can greatly improve monitoring activities. For instance, the sole observation of an RGB map is highly valuable, because it provides experts with an overview of the entire field, improving their awareness of the disease spread ( Figure A1 in Appendix A).
Multispectral data can be used to speed up the assessment of the symptomatic plants due to their strong correlation with plant vigour. By identifying dead or highly compromised plants, it was certainly possible to precisely delimit and monitor the spread of the disease, while the detection of very vigorous plants made it possible to identify homogeneous areas, where other analyses should be focused for the prediction of disease outbreaks. Our findings are in line with those of other studies that used NIR (850-1700 nm), green (495−570 nm) and red (620-750 nm) wavelengths to develop vegetation indices for the estimation of plant biophysical traits [14,[53][54][55]. Indeed, healthy plants are usually characterised by a higher reflectance in the NIR bands [56]. Studies of multispectral data regarding soil-borne diseases identified NIR bands as an important factor for the detection of plant vigour and the estimation of disease severity indices [57][58][59][60][61]. However, as observed in our study, multispectral data seem to have a poor predictive capability for this kind of diseases, and they are usually useful once symptoms are visible [34,36,[57][58][59][60][61].
The lack of predictive capability within the multispectral data can be offset by thermal data. The canopy temperature was found to be a reliable predictor of plant health status, being able to indicate the spread of the disease one year before the actual outbreak. It is interesting to note that the best prediction metrics were obtained precisely for plants that were completely asymptomatic in 2017, far exceeding the predictive capability of any expert. Indeed, these plants were still apparently healthy in both October 2017 and in the first half of the vegetative season in 2018. Until June 2018, they produced a wide canopy, but suddenly after the first heat waves of July, they showed scorching, leaf drop and, in the worst cases, complete defoliation. Results obtained from the clustering of thermal data confirmed that plant responses to KD might be similar to those induced by drought stress, as also observed for other soil-borne diseases [23,34,36]. Indeed, just like abiotic factors, pathogens may affect the stomatal response by influencing the temperature gradient between the plant tissue and the air [62]. Soil-borne diseases cause the reduction of water absorption translocation and transpiration functions as a net result of their infection [63]. These alterations induce a closure of the stomata, and consequently the increase of leaf temperature as a consequence of the reduced evaporative cooling effect [24,63,64]. Indeed, for soil-borne diseases, canopy temperature has been shown to be particularly useful in detecting compromised plants, even during the early stages of infection when visual symptoms were invisible [23,34,36,63,64]. It can be speculated that the multispectral data probably failed to predict disease outbreaks because KD does not cause internal structure modifications until a critical point is reached, at which time the root system can no longer support the plant's transpiration rate [65]. It seems that the kiwifruit has a very high root/canopy ratio [25,66], and thus can cope with a substantial (80%) loss of the root system before the growth of shoots and leaves is affected [25]. Conversely, a justification for the predictive capability of canopy temperature data may reside in the physiological response of kiwi vines to drought stress or root loss, which induce a rapid reduction in gas exchange fluxes and consequently an increase in leaf temperature [24,25,67]. Finally, it should also be noted that even flooding conditions can induce a reduction in the transpiration rates of kiwi vines [27], and thus the water content of the soil cannot be overlooked during temperature data acquisition for KD detection.
Based on our results, the use of unsupervised clustering can be a reliable method for quickly exploring and identifying sampling areas that are suitable for aetiological studies. By using multispectral data, symptomatic plants can be easily identified, simplifying the disease assessment process. The clustering of thermal data is even more useful, since the early detection of diseased plants allows us not only to identify the best areas where samples for laboratory analysis should be collected, but also to observe the evolution of symptoms over time, and improve the understanding of root degradation dynamics. However, we believe that the unsupervised clustering of remote sensing data might also be implemented for other diseases. To the best of our knowledge, unsupervised clustering has not been used previously for forecasting purposes, although it has been used for the assessment of cotton root rot [68], for structuring the aggressiveness of fungal infections on peas [69], and for the segmentation of plant backgrounds for weed detection [70]. With our experimental set up, K-means gave better results than hierarchical clustering in both assessing and predicting the disease spread. However, it cannot be ignored that under other conditions the performances of these two methods might differ. The biggest drawback of the unsupervised classification methods is the need for an expert to assign meaning to each cluster. Nevertheless, for KD, cluster interpretation can be easily achieved by overlapping cluster results with an RGB orthomosaic. The use of supervised machine learning algorithms might be a solution to the problem of directly classifying plant health status, but when the aetiological background is not clear or the labelling is impractical (such as in the KD syndrome), the development of a reliable training dataset is difficult because rapid tests to discriminate between diseased and healthy plants are currently unavailable. Therefore, unsupervised clustering is an appropriate method of seeing or finding groups within unknown data when labels are not available, especially when we do not know what kind of sensor data is needed to find patterns and groups [71].
The study was limited by low cost technologies, and it was performed by only analysing shifts that occurred between two seasons, so several future improvements can be envisaged with this approach. A better segmentation algorithm, based on sensor fusion processes and improved geometric accuracy of the orthomosaics, is needed to remove background noise and to improve the quality of the analysis. In this regard, a similar approach was successfully adopted to detect olive plants infected by Verticillium wilt, using airborne data captured simultaneously with thermal, fluorescence and hyperspectral detectors [34]. Moreover, the advantages derived from sensor fusion are already evident in fruit safety and quality control studies [72][73][74]. Segmentation could also be improved by using automatic tree segmentation with object-based image analysis [75]. In our study, the low-cost multispectral sensors available precluded the possibility of correctly segmenting the inter-row grass from the kiwifruit canopy. Nevertheless, testing flexible and low-cost methods is important to ensure the real application of remote sensing technologies to the field scale. A better understanding of the relationship between temperature and plant health status is also needed. Models that are able to predict leaf temperature using data from weather stations (e.g., air temperature, air relative humidity, radiation, soil water content) should be developed to set more objective thresholds for discriminating healthy from diseased plants, and thus use temperature as a more reliable predictor. Temporal analysis of simple vegetation indices might also improve our understanding of the disease dynamics and increase the detection accuracy of KD. Studies on a wider range of wavelengths are also needed to identify the best bands for disease detection and prediction. In particular, it would be interesting to study the possible correlations between KD appearance and short-wave infrared wavelengths (1100 to 2500 nm), as these bands are influenced by leaf chemical composition and water content [76,77].

Conclusions
This study is the very first to demonstrate the prediction of the spread of Kiwifruit Decline using remote sensing technologies. Our results indicate that unsupervised clustering can be a reliable algorithm for characterising the disease in its early stages, in order to identify homogenous areas out in the field. Multispectral data can be used to discriminate symptomatic from asymptomatic plants, allowing a quick estimation of disease spread. On the other hand, thermal data have been shown to be (Italy)-for his collaboration and support in field-scouting activities. Luca Zuliani, from Adron Technology srl, for his help in the pre-processing of images and generation of the orthomosaics. Gianni Tacconi from CREAthe Council for Agricultural Research and Agricultural Economics Analysis-for his insights regarding the spread of Kiwifruit Decline in Italy.