Detection of Longhorned Borer Attack and Assessment in Eucalyptus Plantations Using UAV Imagery

: Eucalyptus Longhorned Borers (ELB) are some of the most destructive pests in regions with Mediterranean climate. Low rainfall and extended dry summers cause stress in eucalyptus trees and facilitate ELB infestation. Due to the di ﬃ culty of monitoring the stands by traditional methods, remote sensing arises as an invaluable tool. The main goal of this study was to demonstrate the accuracy of unmanned aerial vehicle (UAV) multispectral imagery for detection and quantiﬁcation of ELB damages in eucalyptus stands. To detect spatial damage, Otsu thresholding analysis was conducted with ﬁve imagery-derived vegetation indices (VIs) and classiﬁcation accuracy was assessed. Treetops were calculated using the local maxima ﬁlter of a sliding window algorithm. Subsequently, large-scale mean-shift segmentation was performed to extract the crowns, and these were classiﬁed with random forest (RF). Forest density maps were produced with data obtained from RF classiﬁcation. The normalized di ﬀ erence vegetation index (NDVI) presented the highest overall accuracy at 98.2% and 0.96 Kappa value. Random forest classiﬁcation resulted in 98.5% accuracy and 0.94 Kappa value. The Otsu thresholding and random forest classiﬁcation can be used by forest managers to assess the infestation. The aggregation of data o ﬀ ered by forest density maps can be a simple tool for supporting pest management.


Introduction
Eucalyptus Longhorned Borers (ELB), Phoracantha semipunctata (Fabricius), and P. recurva Newman (Coleoptera: Cerambycidae), are among the most destructive eucalypt pests in regions with Mediterranean climate [1,2]. Eucalyptus globulus is one of the most planted eucalypt species in these regions, and it is known to have low resistance to ELB [2,3].
ELB activity generally starts in late spring when adults emerge and start laying eggs. Upon hatching, ELB larvae bore galleries along the phloem and cambium of trees, eventually preventing sap from flowing, which leads to rapid tree death during summer and fall [1,[4][5][6]. The ability of larvae to successfully colonize the host plant depends on low bark moisture content, leaving water-stressed trees particularly susceptible to attack by ELB [1,5,7].
Under the Mediterranean climate, with low rainfall and extended dry summers periods, the attack by ELB often result in significant tree mortality, despite long-lasting efforts to select for more resistant E. globulus genotypes [8]. With droughts expected to increase due to climate change, ELB outbreaks will likely become more frequent and more severe [9].
ELB control methods include biological control, selection of more resistant eucalypts, and various cultural practices aimed at increasing tree adaptation and resilience, but the most effective curative measure is felling all attacked trees and removing them from the stands [6,10]. Damage caused by ELB is usually not detected using traditional surveillance techniques until significant mortality has occurred [11]. Early detection is recognized as the first step to reduce the impact of ELB [1], hence new approaches to monitoring, particularly through remote sensing, can provide an invaluable tool for forest managers in terms of planning and executing control actions.
Traditional survey techniques are restricted by small area coverage and subjectivity [12] but combined with remote-sensing technology can lead to expanded spatial coverage, minimize the response time, and reduce the costs of monitoring forested areas [13]. Following Wulder et al. [14], the appropriate sensor and resolution should be adopted according to the spatial scale which better adjusts to the situation. To date, several studies have demonstrated successful pest and diseases detection and monitoring in forests using different types of sensor and platform [15][16][17]. For instance, in terms of spaceborne optical sensors for large areas, Meddens et al. [18] detected multiple levels of coniferous tree mortality using multi-date and single-date Landsat imagery. Mortality in ash trees was assessed by Waser et al. [19] using multispectral WorldView-2 imagery. In the eucalypt forest context, to predict bronze bug damage Oumar and Mutanga [20] tested WorldView-2 imagery. The airborne optical sensor Compact Air-Borne Spectrographic Imager 2 (CASI-2) has been tested by Stone et al. [21] who conducted a study to assess damage caused by herbivorous insects in Australian eucalypt forests and pine plantations.
In recent years the use of unmanned aerial vehicle (UAV) platforms has become widely employed for pests and disease detection and monitoring [12,[22][23][24][25][26][27][28]. The main attributes are very high resolution, suitability for multitemporal analysis, lower operational costs compared with airplanes and satellites, independence of cloud cover, and the ability to operate in specific phenological phases of plants or pest/disease outbreaks [12]. In addition, a large number of passive and active sensors can be assembled, such as RGB (red, green, blue), multispectral and hyperspectral cameras, LiDAR (laser imaging detection and ranging), and RADAR (radio detection and ranging) [29,30]. On the other hand, as stressed by Pádua et al. [29], the disadvantages include small area-coverage when compared with satellites, sensitivity to bad weather, increasingly strict regulations that may restrict operations, and the high volume of data generated.
The imagery classification approaches used in the more recently published UAV studies to detect pests and diseases in forests are diverse. Lehmann et al. [22] used a modified normalized difference vegetation index (NDVI) to discriminate between five classes of infestation by the oak splendor beetle through OBIA (objected-based imagery analysis) classification with an overall Kappa index of agreement of 0.81-0.77. Näsi et al. [23] investigated bark beetle damage at the tree level in South Finland using a combination of RGB, near infrared (NIR), red-edge band and NDVI. Object-based K-NN (K-nearest neighbor) classification was applied, and overall accuracy was 75%. In New Zealand, Dash et al. [12] used NDVI and red-edge NDVI to study the discoloration classes of Pinus radiata through the OBIA classification approach with a random forest (RF) algorithm. In Catalonia (Spain), Otsu et al. [26] detected defoliation of pine trees affected by pine processionary and distinguished pine species at a pixel level using four spectral indices and NIR band. The authors estimated the threshold values using histogram analysis. For comparisons, another classification used was the OBIA with a random forest algorithm and an overall accuracy of 93%. Finally, Iordache et al. [28] studied pine wild disease in Portugal using OBIA and machine-learning random forest algorithm for multispectral and hyperspectral imagery. The overall accuracy of both classifications was 95% and 91%, respectively. To date, no published work has been developed with UAV imagery for ELB detection. This work arises from the necessity to find monitoring tools that could support pest management decisions regarding ELB attacks in eucalypt stands. In the past, this has been hampered since the identification of dead trees in the field is bothcostly and time-consuming. Furthermore, it is critical that ELB attacks can be identified at early stages of infestation. Current surveying methods only detect problems when large number of trees are infested.
The main objectives of this experimental work were: (1) to detect ELB attacks through UAV imagery by using selected spectral indices and Otsu thresholding; (2) to map trees crown status using large-scale mean-shift (LSMS) segmentation, as well as the machine-learning classification approach with a RF algorithm; and (3) to map tree density using the hexagonal tessellations technique to support phytosanitary interventions.

Study Area
A 30.56 ha E. globulus stand located in central Portugal close to Gavião locality (39 • 27.909 N, 7 • 56.124 W) was selected for this study ( Figure 1). The study area is managed by The Navigator Company (NVG), a Portuguese pulp and paper company, and was selected based on reports of severe pest attacks provided by forest managers of the company.
Remote Sens. 2020, 12, x FOR PEER REVIEW  3 of 21 This work arises from the necessity to find monitoring tools that could support pest management decisions regarding ELB attacks in eucalypt stands. In the past, this has been hampered since the identification of dead trees in the field is bothcostly and time-consuming. Furthermore, it is critical that ELB attacks can be identified at early stages of infestation. Current surveying methods only detect problems when large number of trees are infested.
The main objectives of this experimental work were: (1) to detect ELB attacks through UAV imagery by using selected spectral indices and Otsu thresholding; (2) to map trees crown status using large-scale mean-shift (LSMS) segmentation, as well as the machine-learning classification approach with a RF algorithm; and (3) to map tree density using the hexagonal tessellations technique to support phytosanitary interventions.

Study Area
A 30.56 ha E. globulus stand located in central Portugal close to Gavião locality (39°27.909′ N, 7°56.124′ W) was selected for this study ( Figure 1). The study area is managed by The Navigator Company (NVG), a Portuguese pulp and paper company, and was selected based on reports of severe pest attacks provided by forest managers of the company. The eucalypt plantation was established in 2014 with a mean stand density of 1250 plants per hectare. Initial mortality caused by ELB was detected in early 2018. Due to the severe drought that occurred in 2017, ELB attack increased over the course of a few months. This is supported by the occurrence of low rainfall in 2017, especially during the summer months, as shown in Figure 2.
Annual average rainfall is low in this region (700-800 mm/year), so if there is a year drier than the average, the following year's mortality is likely to increase due to water stress and ELB attacks. The eucalypt plantation was established in 2014 with a mean stand density of 1250 plants per hectare. Initial mortality caused by ELB was detected in early 2018. Due to the severe drought that occurred in 2017, ELB attack increased over the course of a few months. This is supported by the occurrence of low rainfall in 2017, especially during the summer months, as shown in Figure 2. This area is characterized by sandstones both in valleys and hillsides. The soils are generally poor and shallow. According to the Forest and Agriculture Organization (FAO)classification system [31], soils are Leptosols or Plinthic Arenosols. The relief is slightly inclined, and the altitude of the site varies between 160 and 250 m a.s.l. (above sea level). Southern aspects prevail in the study site.

Data Acquisition
Multispectral images were acquired through four flights on 21 January 2019 using a Parrot SEQUOIA camera (Parrot S.A., Paris, France) on a remotely piloted fixed wing eBee SenseFly drone (Parrot S.A., Paris, France). This single propeller drone allows more stable and longer flights, which in turn reduces the cost and length of missions. Considering the frequency of cloudy days and strong winds during winter season, this aerial vehicle has greater speed in the acquisition of data and efficiency in the operation. The details of flight planning parameters are presented in Table 1. The Parrot SEQUOIA is a compact bundle with multispectral and sunshine sensors. The camera collects four discrete and separated spectral bands with 1.2 megapixels (MP) resolution: green (530-570 nm), red (640-680 nm), red-edge (730-740 nm), and near-infrared (770-810 nm). Additionally, the RGB camera captures 16 MP of resolution [32,33].
To improve equal resolution and less likelihood of holes at higher elevation, terrain awareness was used (Figure 3a,b). The mission planning application used was SenseFly eMotion (Parrot S.A., Paris, France). To refine vertical and horizontal accuracy, nine ground control points (GCP's) were installed and measured with real-time kinematic (RTK) global navigation satellite systems (GNSS). Annual average rainfall is low in this region (700-800 mm/year), so if there is a year drier than the average, the following year's mortality is likely to increase due to water stress and ELB attacks. This area is characterized by sandstones both in valleys and hillsides. The soils are generally poor and shallow. According to the Forest and Agriculture Organization (FAO)classification system [31], soils are Leptosols or Plinthic Arenosols. The relief is slightly inclined, and the altitude of the site varies between 160 and 250 m a.s.l. (above sea level). Southern aspects prevail in the study site.

Data Acquisition
Multispectral images were acquired through four flights on 21 January 2019 using a Parrot SEQUOIA camera (Parrot S.A., Paris, France) on a remotely piloted fixed wing eBee SenseFly drone (Parrot S.A., Paris, France). This single propeller drone allows more stable and longer flights, which in turn reduces the cost and length of missions. Considering the frequency of cloudy days and strong winds during winter season, this aerial vehicle has greater speed in the acquisition of data and efficiency in the operation. The details of flight planning parameters are presented in Table 1. The Parrot SEQUOIA is a compact bundle with multispectral and sunshine sensors. The camera collects four discrete and separated spectral bands with 1.2 megapixels (MP) resolution: green (530-570 nm), red (640-680 nm), red-edge (730-740 nm), and near-infrared (770-810 nm). Additionally, the RGB camera captures 16 MP of resolution [32,33].
To improve equal resolution and less likelihood of holes at higher elevation, terrain awareness was used (Figure 3a,b). The mission planning application used was SenseFly eMotion (Parrot S.A., Paris, France). To refine vertical and horizontal accuracy, nine ground control points (GCP's) were Remote Sens. 2020, 12, 3153 5 of 20 installed and measured with real-time kinematic (RTK) global navigation satellite systems (GNSS). Then, the multispectral camera was calibrated in loco, to adjust the sunshine sensors to the local lighting conditions at solar noon.
Then, the multispectral camera was calibrated in loco, to adjust the sunshine sensors to the local lighting conditions at solar noon.
Regarding the imagery processing workflow, all discrete bands imagery and GCPs were imported to Pix4Dmapper Pro software (Version 4.2, Pix4D S.A., Prilly, Switzerland), a photogrammetry and drone mapping software broadly used for UAV imagery processing [34]. The agricultural (multispectral photogrammetry workflow was chosen, which consists in alignment, geometric calibration based on GCP's, reflectance calibration, point cloud generation, and classification, raster digital surface model (DSM) and orthomosaic generation based on the digital terrain model (DTM). The next step was to analyze the vitality of the stand in the field. UAV multispectral imagery was used to elaborate color composites in which matching of the tree canopies pixel were compared with the real field status of trees. For this end, two different tree vitality status were assigned: healthy trees ( Figure 4a) and dead trees-crown completely dead ( Figure 4b). Field data collection was carried out through random sampling. Data collection was performed using Arrow Gold Antenna (Eos Positioning Systems, Inc., Terrebonne, QC, Canada) with RTK sensor. The cause of death was verified by removing the bark in order to find signs of ELB attack, namely galleries caused by larvae ( Figure 4c). Regarding the imagery processing workflow, all discrete bands imagery and GCPs were imported to Pix4Dmapper Pro software (Version 4.2, Pix4D S.A., Prilly, Switzerland), a photogrammetry and drone mapping software broadly used for UAV imagery processing [34]. The agricultural (multispectral photogrammetry workflow was chosen, which consists in alignment, geometric calibration based on GCP's, reflectance calibration, point cloud generation, and classification, raster digital surface model (DSM) and orthomosaic generation based on the digital terrain model (DTM).
The next step was to analyze the vitality of the stand in the field. UAV multispectral imagery was used to elaborate color composites in which matching of the tree canopies pixel were compared with the real field status of trees. For this end, two different tree vitality status were assigned: healthy trees ( Figure 4a) and dead trees-crown completely dead ( Figure 4b). Field data collection was carried out through random sampling. Data collection was performed using Arrow Gold Antenna (Eos Positioning Systems, Inc., Terrebonne, QC, Canada) with RTK sensor. The cause of death was verified by removing the bark in order to find signs of ELB attack, namely galleries caused by larvae ( Figure 4c).
Once the status of the trees was assigned in the field, spectral values were extracted for the training and validation datasets.

Selected Vegetation Indices
Among the considerable number of vegetation indices (VIs) usually used as indicators of vegetation status [35], there are a set of VIs that are generally applied to assess symptoms of stress [12,36,37]. Stressed eucalypt trees reveal chlorophyll content reduction, red discoloration produced by secondary metabolites accumulation such as anthocyanins and carotenoids, or loss of photosynthetic tissues due to defoliation or necrosis [38,39]. Commonly, NIR, red, red edge, and green bands are used to assess vegetation status due to high sensitivity to changes in moisture content, pigment indices, and vegetation health, respectively [20].
Based on their ability to reveal stress in eucalypts, a total of five vegetation indices were calculated using the available spectral reflectance bands from UAV imagery ( Table 2). The difference vegetation index (DVI) and NDVI were estimated based on the evidence of correlation with plant stress [40,41]. The green normalized difference vegetation index (GNDVI) and normalized difference red-edge (NDRE) were also estimated, as both have shown high sensitivity to changes in chlorophyll concentrations [20,[42][43][44]. In addition, some studies have shown that GNDVI has higher sensitivity to Remote Sens. 2020, 12, 3153 6 of 20 pigment changes than NDVI [45,46]. The soil-adjusted vegetation index (SAVI) was also estimated to reduce soil brightness influence [47].
The next step was to analyze the vitality of the stand in the field. UAV multispectral imagery was used to elaborate color composites in which matching of the tree canopies pixel were compared with the real field status of trees. For this end, two different tree vitality status were assigned: healthy trees ( Figure 4a) and dead trees-crown completely dead (Figure 4b). Field data collection was carried out through random sampling. Data collection was performed using Arrow Gold Antenna (Eos Positioning Systems, Inc., Terrebonne, QC, Canada) with RTK sensor. The cause of death was verified by removing the bark in order to find signs of ELB attack, namely galleries caused by larvae (Figure 4c). The values of the spectral indices were extracted for all features collected in the field using the zonal statistics tool provided by QGIS (Version 3.10) [48]. Table 2. Selected vegetation indices derived from unmanned aerial vehicle (UAV) imagery.

Otsu Thresholding Method
The Otsu thresholding method is an unsupervised method used to select a threshold automatically from a gray level histogram [52]. This approach assumes that the image contains two classes of pixels following a bimodal histogram, and estimates a threshold value which splits the foreground and the background of an image [26,[52][53][54]. The threshold value is calculated aiming to minimize intra-class intensity variance [52]. Although this method is usually applied to images with bimodal histograms, it might be also used for unimodal and multimodal histograms if the precision of the objects is not a requirement [53]. As stressed by Otsu [52], when the histogram valley is flat and broad, it is difficult to detect the threshold with precision. This phenomenon frequently occurs with real pictures and remote sensing imagery. Generally, as eucalypt stands are aligned and understory vegetation is managed, tree and bare soil reflectance can be easily discriminated. Hence, by splitting soil and vegetation, gaps of eucalypt trees would be more visible. In this study, this method was applied in vegetation indices (VIs)to separate healthy and dead eucalypt. Using the Scikit-image library [55] in Python (Version 3.6.9) [56], the histogram threshold values were calculated.

Local Maxima of a Sliding Window
The extraction of tree location can be performed on high-resolution imagery by applying local maxima filter, as stressed by Wulder et al. [57][58][59] and Wang et al. [60]. Therefore, this algorithm enables the estimation of tree density, volume, basal area, and other important parameters used in forest management and planning [57]. The case of pest management is of special interest since the algorithm can be combined with other techniques in order to determine how many trees are infested or dead. To find the number of trees in the multispectral image, the QGIS Tree density plugin (Ghent University, Lieven P.C. Verbeke, Belgium) [61] was used on a brightness image calculated through the mean of the four bands. According to Crabbé et al. [61] this algorithm uses a sliding window to move over the image. The pixel is marked as a local maximum when the central pixel is the brightest of the window. From Figure 5, highlighting of local maxima filtering using a sliding window can be observed.

Local Maxima of a Sliding Window
The extraction of tree location can be performed on high-resolution imagery by applying local maxima filter, as stressed by Wulder et al. [57][58][59] and Wang et al. [60]. Therefore, this algorithm enables the estimation of tree density, volume, basal area, and other important parameters used in forest management and planning [57]. The case of pest management is of special interest since the algorithm can be combined with other techniques in order to determine how many trees are infested or dead. To find the number of trees in the multispectral image, the QGIS Tree density plugin (Ghent University, Lieven P.C. Verbeke, Belgium) [61] was used on a brightness image calculated through the mean of the four bands. According to Crabbé et al. [61] this algorithm uses a sliding window to move over the image. The pixel is marked as a local maximum when the central pixel is the brightest of the window. From Figure 5, highlighting of local maxima filtering using a sliding window can be observed. In the next step, all false positives (Figure 5b) were removed through the interpretation of the false-color composite image (Figure 5c). The position of the trees will be used to estimate the number of individual dead and healthy trees. Finally, local maxima points selected will be also used to extract the crowns from the segmented image using spatial tools.

Object-Based Analysis and Classification
Geographic object-based image analysis (GEOBIA), as an extension of object-based analysis, consists on dividing imagery into different meaningful image-objects that group neighboring pixels with similar characteristics or semantic meanings such as spectral, shape (geometric), color, intensity, texture, or contextual measures [62][63][64]. As a result of the GEOBIA procedure, a vector file with several polygons or segments with similar properties is obtained instead of individual pixels. This approach is convenient for individual tree classification [65,66], reducing the intra-class spectral variability caused by crown textures, gaps, and shadows [67,68]. On the other hand, the main limitation of this method is related to the over-and under-segmentation, which affect the subsequent classification process. These errors occur when the generated segments do not represent the real shape and area of the image objects [69]. In the next step, all false positives (Figure 5b) were removed through the interpretation of the false-color composite image (Figure 5c). The position of the trees will be used to estimate the number of individual dead and healthy trees. Finally, local maxima points selected will be also used to extract the crowns from the segmented image using spatial tools.

Object-Based Analysis and Classification
Geographic object-based image analysis (GEOBIA), as an extension of object-based analysis, consists on dividing imagery into different meaningful image-objects that group neighboring pixels with similar characteristics or semantic meanings such as spectral, shape (geometric), color, intensity, texture, or contextual measures [62][63][64]. As a result of the GEOBIA procedure, a vector file with several polygons or segments with similar properties is obtained instead of individual pixels. This approach is convenient for individual tree classification [65,66], reducing the intra-class spectral variability caused by crown textures, gaps, and shadows [67,68]. On the other hand, the main limitation of this method is related to the over-and under-segmentation, which affect the subsequent classification process. These errors occur when the generated segments do not represent the real shape and area of the image objects [69].
Image segmentation was performed by applying the LSMS algorithm with four spectral bands stacked. The LSMS algorithm, developed by Fukunaga and Hostetler in 1975, is a non-parametric and iterative process that groups pixels with similar meanings [70] by using the radiometric mean and variance of each band [71]. The selection of the LSMS algorithm was motivated by three main reasons: (1) LSMS can be used by the ORFEO ToolBox (OTB) (CNES, Paris, France), which is an easy to use open-source project for remote-sensing imagery processing and analyses [72]; (2) LSMS has been especially developed to be applied in large very high resolution (VHR) images processing [70]; and (3) previous work has shown good results when using UAV imagery [68].
The segmentation technique is not a fully automatized process since its procedure may not create the desired "meaningful" segments. Hence, segmentation parameters are commonly "optimized" to achieve the intended entities [63]. For this aim, different segmentation parameters were tested using OTB version 7.1.0 and evaluated by visual observation. The chosen parameters were spatial radius of 5 (default value), range radius of 10, and minimum segment size of 30. Subsequently the means of bands and VI values in each segment were addressed. From Figure 6, visual comparison of the segmentation parameters tested can be ascertained.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 21 Image segmentation was performed by applying the LSMS algorithm with four spectral bands stacked. The LSMS algorithm, developed by Fukunaga and Hostetler in 1975, is a non-parametric and iterative process that groups pixels with similar meanings [70] by using the radiometric mean and variance of each band [71]. The selection of the LSMS algorithm was motivated by three main reasons: (1) LSMS can be used by the ORFEO ToolBox (OTB) (CNES, Paris, France), which is an easy to use open-source project for remote-sensing imagery processing and analyses [72]; (2) LSMS has been especially developed to be applied in large very high resolution (VHR) images processing [70]; and (3) previous work has shown good results when using UAV imagery [68].
The segmentation technique is not a fully automatized process since its procedure may not create the desired "meaningful" segments. Hence, segmentation parameters are commonly "optimized" to achieve the intended entities [63]. For this aim, different segmentation parameters were tested using OTB version 7.1.0 and evaluated by visual observation. The chosen parameters were spatial radius of 5 (default value), range radius of 10, and minimum segment size of 30. Subsequently the means of bands and VI values in each segment were addressed. From Figure 6, visual comparison of the segmentation parameters tested can be ascertained. Given the most suitable segmentation result, one more step was undertaken to improve the segments that correspond to tree canopies. This optimization consisted of obtaining tree position by applying the local maxima algorithm included in the Tree density plugin of QGIS [61]. By using this tool, bare soil and shadows were eliminated, as they might interfere with the identification of dead trees, which had mostly already lost their leaves.
In order to classify tree canopies into two different classes, healthy and dead trees, supervised machine-learning (ML) classification was used. Among the different object-oriented ML classifiers, the RF algorithm was applied as its performance is one of the most accurate ML algorithms when supervised classification for GEOBIA is conducted [27,28,68,73]. Presented by Breiman [74], RF is an automatic ensemble method based on decision trees where each tree depends on a collection of random variables [75]. RF consists of a large number of independent individual decision trees working together and the final output is estimated based on the outputs of all decision trees involved. Thereby, the return classification is the one that has been the most recurrent. Random selection of the variables used in each decision tree is a keystone to avoid correlation among decision trees. For this purpose, the so-called bootstrap is performed and a significant percentage of samples of the training areas with replacements are used to construct each decision tree. The remaining training areas, called out-of-bag (OOB) data, are used to validate the classification accuracy of RF [74,76]. Therefore, all the decision trees are always constructed using the same parameters but on different training sets. Given the most suitable segmentation result, one more step was undertaken to improve the segments that correspond to tree canopies. This optimization consisted of obtaining tree position by applying the local maxima algorithm included in the Tree density plugin of QGIS [61]. By using this tool, bare soil and shadows were eliminated, as they might interfere with the identification of dead trees, which had mostly already lost their leaves.
In order to classify tree canopies into two different classes, healthy and dead trees, supervised machine-learning (ML) classification was used. Among the different object-oriented ML classifiers, the RF algorithm was applied as its performance is one of the most accurate ML algorithms when supervised classification for GEOBIA is conducted [27,28,68,73]. Presented by Breiman [74], RF is an automatic ensemble method based on decision trees where each tree depends on a collection of random variables [75]. RF consists of a large number of independent individual decision trees working together and the final output is estimated based on the outputs of all decision trees involved. Thereby, the return classification is the one that has been the most recurrent. Random selection of the variables used in each decision tree is a keystone to avoid correlation among decision trees. For this purpose, the so-called bootstrap is performed and a significant percentage of samples of the training areas with replacements are used to construct each decision tree. The remaining training areas, called out-of-bag (OOB) data, are used to validate the classification accuracy of RF [74,76]. Therefore, all the decision trees are always constructed using the same parameters but on different training sets.
To supervise RF classification, 1241 out of 25,911 segments were manually selected as training areas based on field data and on-screen interpretation of different image color composites. A stratified division of the training areas was carried out according to the preset classes being divided as 80% for training and 20% for the model validation. As shown in previous studies [68,77], default values were set in OTB when RF was performed, since OTB parameters for training and classification processes worked optimally.
The summary of the performance of the machine-learning RF training model can be observed in Table 3. For each canopy class, high precision and global model performance were obtained, as shown by the Kappa value of 0.96.

Accuracy Assessment
Since it is a very effective way to represent both global accuracy and individual class accuracy, an error matrix was performed to determine the agreement between the classified image and the ground truth data with regard to tree health status [78]. Field data collection through stratified random sampling was conducted for validation purposes. From the error matrix, the overall accuracy, commission error (CE), omission error (OE), producer's accuracy (PA), and user's accuracy (UA) were established. CE describes the error associated with the aerial surveyors that classify objects into a different category from the one they belong to, whereas OE represents the error associated with the aerial surveyors that place objects outside of the category that they belong to [79]. In addition, the Kappa (к) statistic was estimated to quantitatively assess the agreement in the error matrix and the chance of an agreement, by determining the degree of association between the remotely sensed classification and the reference data [78,80].
The accuracy assessment of the classifications obtained with Otsu's method was validated with the data collected in the field. The supervised classification with the random forest was validated with 2% (517 samples) of randomly stratified selected segments using the research tools provided by QGIS [48].

Density Forest Maps
The results of the classification were joined to the tree's position layer due to the slightly touching trees phenomenon described by Wang et al. [60]. Counting the number of trees using the classification results could induce error because the tree crows are not always clearly separated.
Density maps were created by grouping the position of the trees with the hexagonal binning technique for 0.1 ha. This strategy serves to minimize the irregular shape of the stand limits and to group the trees into less than discrete groups. More importantly, hexagons are more similar to the circle than squares or diamonds [81]. This technique was first described by Carr et al. [82]. According to Carr [83] maps based on hexagon tessellations offer the opportunity to summarize point data and show the patterns for a single variable, that in this case was made into two types of the canopy.

Spectra Details and Vegetation Indexes
Spectra details of tree canopies are shown in Figure 7a, with the mean reflectance values for each band. NIR spectral region revealed the largest differences between the two tree classes, and therefore the highest discriminating power when compared with red-edge and other bands. As shown in Figure 7b, all indices have discriminative capability, although the NDRE index showed the smallest interval difference.

Pixel-Based Analysis
Otsu Thresholding Analysis Imagery interval values and histogram threshold values obtained from the Otsu thresholding method, as shown in Table 4 and Figure A1 (Appendix A). The DVI index imagery interval ranged between 0.020 and 1.150, setting the threshold, value at 0.394. The GNDVI and NDVI indices showed a similar histogram threshold with global minima of 0.703 and 0.712, respectively. Their imagery intervals ranged between 0.150 and 0.930 for the NDVI index and between 0.248 and 0.890 for the GNDVI index. The other two indexes, NDRE and SAVI, had imagery intervals between −0.339 and 0.660 and between 0.057 and 0.983, respectively, whereas their threshold values were 0.215 and 0.526, respectively. Figure 8 highlights the threshold values of all spectral indices by applying the Otsu thresholding method. If compared with the false-color composite in Figure 8a, where red corresponds to healthy trees and blue-grey to dead trees mixed with bare soil, all indices showed good performance to discriminate between healthy and dead trees. Hence, the application of Otsu methodology seemed to improve the differentiation between classes, at least in planted stands where there is active silvicultural management.

Otsu Thresholding Analysis
Imagery interval values and histogram threshold values obtained from the Otsu thresholding method, as shown in Table 4 and Figure A1 (Appendix A). The DVI index imagery interval ranged between 0.020 and 1.150, setting the threshold, value at 0.394. The GNDVI and NDVI indices showed a similar histogram threshold with global minima of 0.703 and 0.712, respectively. Their imagery intervals ranged between 0.150 and 0.930 for the NDVI index and between 0.248 and 0.890 for the GNDVI index. The other two indexes, NDRE and SAVI, had imagery intervals between −0.339 and 0.660 and between 0.057 and 0.983, respectively, whereas their threshold values were 0.215 and 0.526, respectively.  Figure 8 highlights the threshold values of all spectral indices by applying the Otsu thresholding method. If compared with the false-color composite in Figure 8a, where red corresponds to healthy trees and blue-grey to dead trees mixed with bare soil, all indices showed good performance to discriminate between healthy and dead trees. Hence, the application of Otsu methodology seemed to improve the differentiation between classes, at least in planted stands where there is active silvicultural management.

Object-Based Analysis and Classification
Discrimination of tree canopy status was performed through the classification of crown segments extracted from the points generated with local maxima filtering by applying a machinelearning random forest algorithm ( Figure 9). The spectral indices and band information of each segment allowed the identification of tree canopy status (Figure 9b). Large-scale mean-shift segmentation algorithm performed through OTB had good performance, improving the ability to obtain meaningful segments, which in turn resulted in improvements in RF classification capability (Figure 9c).

Object-Based Analysis and Classification
Discrimination of tree canopy status was performed through the classification of crown segments extracted from the points generated with local maxima filtering by applying a machine-learning random forest algorithm ( Figure 9). The spectral indices and band information of each segment allowed the identification of tree canopy status (Figure 9b). Large-scale mean-shift segmentation algorithm performed through OTB had good performance, improving the ability to obtain meaningful segments, which in turn resulted in improvements in RF classification capability (Figure 9c).

Otsu Thresholding
Histogram thresholding analysis and RF classifier were validated by estimating the confusion matrix and Kappa statistic (Tables 5 and 6). The confusion matrix to classify tree status with NDVI histogram analysis had the highest overall accuracy (98.2%) and Kappa value (0.96). On the other hand, NDRE histogram analysis had the lowest performance, with 86.4% overall accuracy and a Kappa value of 0.73. The other spectral indices histogram analyses had a similar performance.

Otsu Thresholding
Histogram thresholding analysis and RF classifier were validated by estimating the confusion matrix and Kappa statistic (Tables 5 and 6). The confusion matrix to classify tree status with NDVI histogram analysis had the highest overall accuracy (98.2%) and Kappa value (0.96). On the other hand, NDRE histogram analysis had the lowest performance, with 86.4% overall accuracy and a Kappa value of 0.73. The other spectral indices histogram analyses had a similar performance.

Object-Based Analysis and Classification
To perform the error matrix for machine-learning RF, 517 predicted observations were used, 431 for healthy trees and 86 for dead trees. Predicted observations were selected by stratified random sampling using QGIS tools ( Table 6). Overall accuracy was 98.5%, which means that this model had the highest accuracy of all the models tested, while the Kappa value was 0.94.

Density Forest Maps
After RF classification, tree canopy condition data (healthy and dead) was validated and two different outputs were created by applying a hexagonal binning. Both products are representations of the classified data in order to support management decisions. Prior to the density forest map creation, the data was already validated so as no new data was generated, the validation procedures may not be required.
Hexagonal binning was created to allow the aggregation of trees in small areas and, therefore, identify hotspots and density of healthy trees ( Figure 10). A total of 422 hexagons were produced with an area of 0.1 ha each, which represents one-tenth of the hectare which is the standard measure. The hexagon with the highest number of healthy trees had 130, while there was only one hexagon where all the trees were dead. Considering the total number of trees classified, 16 hexagons did not show any dead tree whereas a total of 18 hexagons had more than 40 dead trees ( Figure 11). On average, each hexagon was expected to have an initial density of 125 trees, since planting density per hectare was 1250 trees.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 21 Forest mortality maps for the plantation area ( Figure 11) allow the identification of the most critical spots and are important in the planning of dedicated surveying or sanitary fellings to remove infested trees thus reducing ELB local populations and minimizing future outbreaks.

Discussion
This study explored the capability of multispectral images captured with a parrot sequoia camera to detect the mortality caused by ELB on E. globulus plantations. Before establishing a flight plan, the phenology of the trees and the life cycle of the insect must be known, in order to extract as much information as possible. The spatial resolution of the images obtained was 17 cm, which allowed the reduction of spatial scale to the individual tree level.
NIR bands had the greatest discriminating power of the all spectra analyzed and provided more information about different tree canopies. As stressed by Otsu et al. [26], NIR band was more sensitive to defoliation than the red-edge band, as observed in Figure 7b. Regarding the remaining bands, the reflectance of the green band was generally higher than that of red band because of the high The total number of trees recorded in the studied stand was 28,271 of which 4507 were dead. Since total stand area is 30.56 ha, initial stocking would be about 38,000 trees at plantation, representing a 25% difference. This discrepancy may be due the leafless canopies or missing trees from early mortality events, which are common shortly after establishment.
Forest mortality maps for the plantation area ( Figure 11) allow the identification of the most critical spots and are important in the planning of dedicated surveying or sanitary fellings to remove infested trees thus reducing ELB local populations and minimizing future outbreaks.

Discussion
This study explored the capability of multispectral images captured with a parrot sequoia camera to detect the mortality caused by ELB on E. globulus plantations. Before establishing a flight plan, the phenology of the trees and the life cycle of the insect must be known, in order to extract as much information as possible. The spatial resolution of the images obtained was 17 cm, which allowed the reduction of spatial scale to the individual tree level.
NIR bands had the greatest discriminating power of the all spectra analyzed and provided more information about different tree canopies. As stressed by Otsu et al. [26], NIR band was more sensitive to defoliation than the red-edge band, as observed in Figure 7b. Regarding the remaining bands, the reflectance of the green band was generally higher than that of red band because of the high absorption related to chlorophyll content [84]. Spectral indices, whose calculation includes bands such as NDVI, also presented highly discriminating behavior. This VI is particularly used to assess defoliation, and damage at high spatial resolutions using multispectral and RGB cameras mounted on UAV [12,22,23,26,85,86].
Histogram analysis allowed the discrimination between healthy and dead trees. Although the detection of the histogram valley in spectral indices was difficult to determine, this method may become a valid alternative to split the two canopy types. NDVI was the most accurate index at 98.2% and Kappa value of 0.96. The second most accurate was the GNDVI index, which had an overall accuracy of 97.4% and Kappa value of 0.95. NDRE was the less accurate index at 84.4% overall accuracy and Kappa value of 0.73. Similar overall accuracy values were obtained by Otsu et al. [26] to detect defoliation of pine needles by pine processionary in four different locations. However, in our study, the NIR band was not used to remove shadows as Miura and Midokawa et al. [87] applied. To minimize the shadow effect, the flight was undertaken at solar noon, so treetops could be easily captured, with shadows present mostly along the plantation lines.
The combination of segmentation and the maximum filtering location allowed the extraction of the tops of trees and carrying out their classification. Applying this method, shadows were removed from other vegetation and bare soil, which has a reflectance very close to that of dead trees. The RF learning machine obtained 98.5% global accuracy and the Kappa value was 0.94. This precision is explained by the great distance between reflectance values and the differences between the two classes. Iordache et al. [28] applied RF classifier on the classification of Pinus pinaster canopy types (infected, suspicious, and healthy) affected by pine wild and obtained an overall accuracy of 95%. Pourazar et al. [27], obtain as overall accuracy 95.58% using five spectral bands and five indices to detect dead and diseased trees.
The forest density maps produced through hexagon tessellations aimed to group the position of classified trees. The ease of reading allows the identification of the most critical areas with tree mortality and to extract important metrics for forest management such as the total number of trees, the standard deviation, and other landscape metrics. Barreto et al. [88] set a hexagonal grid to represent classes of natural Cerrado vegetation in the the Northeast of Brazil, in an area of about 25,590 km 2 , to study the remaining habitats through quantitative indices. More recently, Amaral et al. [89], used a hexagonal grid subdivided into 1000 ha units to study the restinga (a type of coastal vegetation) in Northern Brazil.
The current strategy to control ELB relies mostly on the identification of dead trees and their removal from plantations before a new generation of adults emerges in late spring. Traditional field surveys are extremely labor-intensive, as they require visual assessment of large eucalypt plantations. As a result, outbreaks are often only detected when large numbers of trees are already dead and ELB populations are high. By allowing the identification of individual infested trees in a much more efficient way, widespread application of the methods described in this study will allow forest managers to detect ELB attacks at an early stage, thus reducing the cost and efficiency of sanitary fellings.
Future work will focus on adding additional classes to the survey in order to improve the discrimination of tree health status. The integration of a UAS-derived canopy model at a 3D tree model could be performed to automatically outline individual tree crowns. Another improvement is the ambition to implement periodical flights at different attack stages to provide a multitemporal analysis. Further research might focus on other remote-sensing platforms at different scales and considering different bioclimatic and geographical settings.

Conclusions
In this experimental research, ELB attack detection and assessment were proposed at pixel and object-based level using UAV multispectral imagery. The NIR band showed discriminative power for two canopy classes. In general, all selected spectral indices were effective in discriminating healthy from dead trees. Otsu thresholding analysis resulted in high overall accuracy and Kappa value for most of the vegetation indices used, with the exception of NDRE. The RF machine-learning algorithm applied also resulted in high overall accuracy (98.5%) and Kappa value (0.94). From an operational point of view, these results can have important implications for forest managers and practitioners. In particular, forest mortality and density maps using hexagonal tessellations allow a complete survey and an accurate measure of the scale of the problem and where the most critical hotspots of mortality are located. By carefully planning the time of the year when these surveys are carried out, timely quantification of mortality rates is an important metric for ensuring phytosanitary fellings are conducted as early and as thoroughly as possible.