Generating the Baseline in the Early Detection of Bud Rot and Red Ring Disease in Oil Palms by Geospatial Technologies

: Oil palm cultivation in Ecuador is important for the agricultural sector. As a result of it, the country generates sources of employment in some of the most vulnerable zones; it contributes 0.89% of the gross domestic product and 4.35% of the agricultural gross domestic product. In 2017, a value of USD $252 million was generated by exports, and palm contributed 4.53% of the agricultural gross domestic product (GDP). It is estimated that 125,000 hectares of palm were lost in the Republic of Ecuador due to Red Ring Disease (RRD) and speciﬁcally Bud Rot (BR). The current study aimed to generate an early detection of BR and RRD in oil palm. Image acquisition has been performed using Remotely Piloted Aircraft System (RPAS) with Red, Green, and Blue (RGB) cannons, and multispectral cameras, in study areas with and without the presence of the given disease. Hereby, we proposed two phases. In phase A, a drone ﬂight has been conducted for processing and georeferencing. This allowed to obtain an orthomosaic that serves as input for obtaining several vegetation indices of the healthy crop. The data and products obtained from this phase served as a baseline to perform comparisons with plantations a ﬀ ected by BR and RRD and to di ﬀ erentiate the palm varieties that are used by palm growers. In phase B, the same process has been applied three times with an interval of 15 days in an a ﬀ ected plot, in order to identify the symptoms and the progress of them. A validation for the diseases detection has been performed in the ﬁeld, by taking Global Positioning System (GPS) points of the palms that presented symptoms of BR and RRD, through direct observation by ﬁeld experts. The inputs obtained in each monitoring allowed to analyze the spatial behavior of the diseases. The values of the vegetation indices obtained from Phase A and B aimed to establish the di ﬀ erences between healthy and diseased palms, with the purpose of generating the baseline of early responses of BR and RRD conditions. However, the best vegetation index to detect the BR was the Visible Atmospherically Resistant Index (VARI).


Introduction
Diseases and pests that affect plants can harm a wide range of cash crops, resulting in significant yield loss [1][2][3]. As a solution to this situation, pesticides have been developed to protect crops from diseases and pests, which increase the cost of production, and the danger, due to the presence of toxic residues in agricultural products. A further approach has been proposed by the timely identification of patches of diseases and pests within the fields for local treatments [4,5]. Detection of pests and diseases The use of remote sensors allows evaluating the state of the palms to obtain an early diagnosis of diseases or infestation of pests based on the symptoms shown in certain places [5,34,35]. Reference [36] has proposed the common hypothesis that oil palm infected with Ganoderma boninense shows observable symptoms at an early stage, and conducted a study to discriminate the oil palm affected by this disease. Among these, a portable field hyperspectral instrument was used to differentiate healthy and infected oil palm [37]. In this respect, there are several studies carried out to determine the spatial behavior of certain affectations in crops, specifically in the study of oil palm, research has been carried out with the application of remote sensors in order to carry out monitoring and detection of diseases [34,[37][38][39]. Research has also been conducted with the aim of establishing predictive models of the spatiotemporal behavior of phenomena associated with the management of plant health, such as the damage caused by defoliators or the incidence of the bud rot disease (BR), through the use of digital processing of images acquired by remote satellite sensors [40]. It has also been considered to use multitemporal analysis methodologies in which areas that have undergone alteration can be verified [41,42]. In this case, the decrease in vegetation indexes is considered, which allows to assume the possible involvement of plants.
The spatial distribution pattern of lethal wilt and Red Ring disease (RRD) has been determined, in addition to their relationships with environmental and nutritional factors of a study area in Colombia [38]. Regarding the analysis of spatial distribution, it has been stated that it is possible to appreciate that there are foci of palms with affectation, and the progression of diseases occurs around these foci (which may be considered to have occurred randomly) [43].
From the inputs purchased by RPAS, spectral reflectance combinations of two or more wavelengths are able to be developed, which are called spectral indices. In this specific case, vegetation indices were used, which are designed to maximize sensitivity to vegetation characteristics and minimize confounding factors, such as ground bottom reflectance, directional, or atmospheric effects [44]. Plants in general reflect a considerably small portion of the incident light in the range of the visible spectrum, while in the blue and red ranges, there is a very small amount (approximately 10%), whereas in the green it is 15%. That's why human eyes see predominantly green in plants. The remaining 90% of the energy in these wavelengths is absorbed by plant pigments to be used by the photosynthesis process, of which absorption occurs 65% by chlorophylls; therefore, the activity is concentrated in the areas of blue and red, while the β-carotene is centered in blue, and phycoerythrin in green and yellow. Finally, the phycocyanin is concentrated in the yellow-red, and a certain part in the ultraviolet. On the contrary, the near infrared region, for the most part (40 to 60%) is scattered upwards or transmitted downwards (50%), and only 10% is absorbed, since the pigments and cellulose of plant cells are transparent to this radiation [45]. The most widely used indices apply the information contained in the reflectance or radiance of the red canopy and near infrared (NIR) [46].
Typical general behavior of vegetation in the electromagnetic spectrum is conditioned by different factors, such as the structure of the leaves and the presence of certain pigments. This is because as the plant ages, its leaves lose chlorophyll, which results in the increase of the reflectivity of yellow and red. In addition to the water content in the tissues, nutritional deficiencies are linked to the concentration of pigments. It is worth mentioning the presence of pests that modify the amount of pigments and water exchange, among others. It is also fundamental to consider the phenological state of the plant, since it will determine the type of information of the vegetation that will be acquired from remote sensing and indices, since slightly different spectral responses are presented [45].
By considering the interaction of phenological states and the presence of pests with the use of spectral indices, the main objective of the current study is to determine significant differences between varieties of palms, according to their phenological state, in order to establish a baseline of plants (unaffectedly for the early detection of BR and RRD diseases in oil palm crops). This occurs by using the statistical analysis of their differences in terms of palm varieties and affectations. Additionally, we performed a multitemporal evaluation of the behavior of BR and RRD in the study palms in order to detect changes in the growth of existing diseases.

Study Area
The study is composed of two plots in Ecuador (Figure 1a,b). The first one (CIPAL plot) is located at (GPS data of the WGS84 reference system being −0.0421, −79.3662), with an area of 52 hectares (Figure 1c). The second one is called Arteaga plot; it is located at (0.0749, −79.41329), having 40 hectares (Figure 1d).

Study Area
The study is composed of two plots in Ecuador (Figure 1a,b). The first one (CIPAL plot) is located at (GPS data of the WGS84 reference system being −0.0421, −79.3662), with an area of 52 hectares (Figure 1c). The second one is called Arteaga plot; it is located at (0.0749, −79.41329), having 40 hectares (Figure 1d). In the Centro de Investigación en Palma Aceitera, Oil Palm Investigation Center (CIPAL) plot ( Figure 2), the incidence of BR was 0.5% in the varieties of Instituto Nacional de Investigaciones Agropecuarias, Agricultural Research National Institute (INIAP), Centre de coopération internationale en recherche agronomique pour le développement, Center for International Cooperation in Agricultural Research for Development(CIRAD), and Agricultural Services and Development from Costa Rica (ASD), of approximately 15 years old, belonging to the phenological state 305 [47], while the Oleifera x Guineensis (OxG) interspecific hybrids, called Taisha, Unipalma, and Amazon, of 2 years old (phenological state 302), have no BR affectation. These plots were considered healthy plots due to a low incidence and the practices of management of the disease [48]. The Arteaga plot has the majority of the INIAP variety, which was BR's incidence of 30% at the beginning of the current study. This plot was the sick one. The lots were considered in a way in which a lot lacks affectation in order to establish the baseline (CIPAL lot) and one that has affectations (Arteaga lot), as this already had affectations of all kinds of degrees. In the Centro de Investigación en Palma Aceitera, Oil Palm Investigation Center (CIPAL) plot ( Figure 2), the incidence of BR was 0.5% in the varieties of Instituto Nacional de Investigaciones Agropecuarias, Agricultural Research National Institute (INIAP), Centre de coopération internationale en recherche agronomique pour le développement, Center for International Cooperation in Agricultural Research for Development (CIRAD), and Agricultural Services and Development from Costa Rica (ASD), of approximately 15 years old, belonging to the phenological state 305 [47], while the Oleifera x Guineensis (OxG) interspecific hybrids, called Taisha, Unipalma, and Amazon, of 2 years old (phenological state 302), have no BR affectation. These plots were considered healthy plots due to a low incidence and the practices of management of the disease [48]. The Arteaga plot has the majority of the INIAP variety, which was BR's incidence of 30% at the beginning of the current study. This plot was the sick one. The lots were considered in a way in which a lot lacks affectation in order to establish the baseline (CIPAL lot) and one that has affectations (Arteaga lot), as this already had affectations of all kinds of degrees. Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 22

Material
Drone flights has been conducted for image acquisition with the DJI Mavic Pro equipment with the MAPIR Survey3W camera Red, Green, and Near Infrared (RGNIR), which includes a calibration card (MAPIR Calibration Ground Target Package V2), and the DJI Matrice 100 drone, with two different Zenmuse X3 cameras, Blue, Green, and Near Infrared (BGNIR), and a Red, Green, Blue (RGB) in order to have a wider data bank and statistical redundancy in the obtained data. The specifications of the cameras are listed in Table 1. All of the cameras are of converted types, which has as a positive aspect due to their low costs, compared to pure NIR cameras (Table 1). A handheld Global Positioning System (GPS) Garmin Ventura with an accuracy of 2 m was used in order to obtain coordinates of samples of oil palm. For the flight planning (for the displaying of the data), we used Pix4D Capture, while we used Agisoft PhotoScan Professional 1.2.6.2834 for orthomosaic processing and atmospheric correction (QUAC ® ), automatic plant count, and the generation of the vegetation indexes ENVI 5.5.

Material
Drone flights has been conducted for image acquisition with the DJI Mavic Pro equipment with the MAPIR Survey3W camera Red, Green, and Near Infrared (RGNIR), which includes a calibration card (MAPIR Calibration Ground Target Package V2), and the DJI Matrice 100 drone, with two different Zenmuse X3 cameras, Blue, Green, and Near Infrared (BGNIR), and a Red, Green, Blue (RGB) in order to have a wider data bank and statistical redundancy in the obtained data. The specifications of the cameras are listed in Table 1. All of the cameras are of converted types, which has as a positive aspect due to their low costs, compared to pure NIR cameras (Table 1). A handheld Global Positioning System (GPS) Garmin Ventura with an accuracy of 2 m was used in order to obtain coordinates of samples of oil palm. For the flight planning (for the displaying of the data), we used Pix4D Capture, while we used Agisoft PhotoScan Professional 1.2.6.2834 for orthomosaic processing and atmospheric correction (QUAC ® ), automatic plant count, and the generation of the vegetation indexes ENVI 5.5.

Methodology
Two phases have been developed to execute this project, which are described (in a general way) in Figure 3. For the photogrammetric flying and planning, 120-m flights with 70% longitudinal and transverse overlaps were carried out in both phases, which is sufficient to generate orthomosaics of high quality. The Ground Sample Distance (GSD) was approximately 5 cm for the BGNIR sensor and approximately 5.5 cm for the RGNIR. Immediately afterwards, the orthomosaic generation process was conducted (internal and external orientation, creation of DEM, creation of the orthomosaic). It should be mentioned that, in both phases, from the first orthomosaic obtained, the others have been georeferenced. The GSD provided by the cameras used varied from 5 to 5.5 cm at a flight height of 120 m. If an effective resolution of three pixels is taken into account, and the discrimination factor is that established by the Ecuadorian Military Geographical Institute (IGM) of 0.3 mm, a maximum scale of 1:500 is obtained, and a minimum mapping unit of 4 to 4.84 m 2 , so that a palm is fully identifiable within this range, since the average palm surface in variety INIAP was 300 m 2 , Taisha 45 m 2 , Unipalma 42 m 2 y, Amazon 21 m 2 [48].

Methodology
Two phases have been developed to execute this project, which are described (in a general way) in Figure 3. For the photogrammetric flying and planning, 120-m flights with 70% longitudinal and transverse overlaps were carried out in both phases, which is sufficient to generate orthomosaics of high quality. The Ground Sample Distance (GSD) was approximately 5 cm for the BGNIR sensor and approximately 5.5 cm for the RGNIR. Immediately afterwards, the orthomosaic generation process was conducted (internal and external orientation, creation of DEM, creation of the orthomosaic). It should be mentioned that, in both phases, from the first orthomosaic obtained, the others have been georeferenced. The GSD provided by the cameras used varied from 5 to 5.5 cm at a flight height of 120 m. If an effective resolution of three pixels is taken into account, and the discrimination factor is that established by the Ecuadorian Military Geographical Institute (IGM) of 0.3 mm, a maximum scale of 1:500 is obtained, and a minimum mapping unit of 4 to 4.84 m 2 , so that a palm is fully identifiable within this range, since the average palm surface in variety INIAP was 300 m 2 , Taisha 45 m 2 , Unipalma 42 m 2 y, Amazon 21 m 2 [48].

Phase A
In this phase, two flights have been made for each sensor (with BGNIR on 18 April 2019 and 20 May 2019, and with RGNIR on 20 May and 22 August) in the CIPAL plot, in order to obtain an extensive database to contrast the results. The orthomosaics obtained allowed to generate vegetation indexes in order to generate vegetation indexes that allow characterizing and obtaining the baseline of affected palms. Once these data are obtained, a statistical analysis is carried out from the ANOVA-Fisher tests to differentiate the different palm varieties. The main objective of using the analysis of variances (ANOVA) is to compare the groups or treatments with each other [49]. Hereby, they consist of the lots, according to each species (INIAP, CIRAD, ASD) present in the CIPAL lot, and as repetitions, the number of lots according to species. The null hypothesis may be stated in such a way

Phase A
In this phase, two flights have been made for each sensor (with BGNIR on 18 April 2019 and 20 May 2019, and with RGNIR on 20 May and 22 August) in the CIPAL plot, in order to obtain an extensive database to contrast the results. The orthomosaics obtained allowed to generate vegetation indexes in order to generate vegetation indexes that allow characterizing and obtaining the baseline of affected palms. Once these data are obtained, a statistical analysis is carried out from the ANOVA-Fisher tests to differentiate the different palm varieties. The main objective of using the analysis of variances (ANOVA) is to compare the groups or treatments with each other [49]. Hereby, they consist of the lots, according to each species (INIAP, CIRAD, ASD) present in the CIPAL lot, and as repetitions, the number of lots according to species. The null hypothesis may be stated in such a way that there are no significant differences between the population means of all the treatments and, as an alternative hypothesis, it is stated that there are significant differences between the treatments, that is, at least two groups differ. Then, in order to verify the hypothesis of equality between means, a Fisher F statistic is calculated, which allows to obtain the degree of similarity that exists between the means that are being compared. In the current study, the Least Significant Difference (LSD) test of Fisher has been used, as it is considered the best option in order to conduct an exploratory analysis of differences between treatments, since it has simplicity and consistency due to the known and constant type I error [50]. In order to corroborate Fisher's LSD test, the Bonferroni method was also used, which is based on the Student's f distribution and on the Bonferroni inequality. This method attempts to solve the problem of applying numerous Student tests by reducing the probability of making a type I error in each comparison. This occurs by the control of the error rate by dividing the level of significance α by the number of performed comparisons k [51]. The data obtained from this phase served as baseline for comparison with plantations affected by BR and DRR.

Phase B
In this phase, flights were made on 30 May 2019, 13 June 2019, and 26 August 2019 in the Arteaga plot, in which can be denoted the presence of the affectation of BR and RRB in its different degrees of severity. Three monitorings were applied in order to analyze the spread of the present diseases.
Based on the description of the methodology (Figure 4), we have to start from the generation of orthomosaics; in both phases, an atmospheric correction was applied using the QUAC ® method in ENVI. This algorithm determines the atmospheric correction parameters directly from the pixel spectra observed in a scene, without auxiliary information. Therefore, it performs a more approximate atmospheric correction than Fast Line-of-sight Atmospheric Analysis of Hypecubes (FLAASH), or other physics-based first principles methods, generally producing spectra reflectance within the range of about 10% of the ground truth. It is based on the empirical finding that the average reflectance of various material spectra (excluding highly structured materials, such as vegetation, shallow water, and mud) does not depend on each scene. Therefore, the processing is much faster compared to the first principles methods. It also allows any viewing angle or solar elevation angle. If a sensor lacks the proper wavelength or radiometric calibration, or if the intensity of the solar illumination is unknown (with cloud cover, among others), reasonably accurate reflectance spectra are still able to be recovered, as long as the following described conditions are met. There are at least 10 different materials in a scene (in this case there is the highway, the road signs, gutters, stony road, six different structures at the entrance, three structures in the location of the laboratories, a structure in the directory, drums).
There are sufficiently dark pixels in a scene in order to allow a good estimate of the baseline spectrum [52]. Next, the Normalized Difference Vegetation Index (NDVI) first Equation (1), Green Normalized Difference Vegetation Index (GNDVI) second Equation (2), Green Vegetation Index (GVI) third Equation (3), and Visual Atmospheric Resistance Index (VARI) fourth Equation (4) indices were calculated with the corresponding channel for the two cameras, depending on the variety of palm being studied. ( An automatic count was applied only in phase A, using the indices that allowed to identify the individual palms with respect to the vegetation or soil elements that surround them, which could affect the response of the calculated index. That is why we studied which index is better to identify the different species. This procedure was performed with the "Crop Science from ENVI" package through the "count crops" section. The results obtained from this type of tools demonstrate a high average precision (greater than 93%) between the counts according to the proposed method and the manual counts (based on RPAS images) [53]. As a result of the count, a polygon-type vector layer was obtained, in which each covers the entire surface of the canopy of a palm. Once the automatic count Remote Sens. 2020, 12, 3229 8 of 21 was performed, a zonal statistic was applied, which allowed calculating the mean of the value of the pixels that fall within each polygon.
This procedure served as an input to spectral differentiation, through statistical analysis of variance tests (ANOVA and Fisher), the different palm varieties present in the study area, and also to obtain the baseline of palms with or without affectation in each case. The assumptions to apply ANOVA (normality, independency, and homoscedasticity) were demonstrated with Shapiro-Wilks modification, Pearson correlation, and Levene's test. There are sufficiently dark pixels in a scene in order to allow a good estimate of the baseline spectrum [52]. Next, the Normalized Difference Vegetation Index (NDVI) first equation (1), Green Normalized Difference Vegetation Index (GNDVI) second equation (2), Green Vegetation Index (GVI) third equation (3), and Visual Atmospheric Resistance Index (VARI) fourth equation (4) indices were calculated with the corresponding channel for the two cameras, depending on the variety of palm being studied. (3) An automatic count was applied only in phase A, using the indices that allowed to identify the individual palms with respect to the vegetation or soil elements that surround them, which could affect the response of the calculated index. That is why we studied which index is better to identify the different species. This procedure was performed with the "Crop Science from ENVI" package through the "count crops" section. The results obtained from this type of tools demonstrate a high average precision (greater than 93%) between the counts according to the proposed method and the manual counts (based on RPAS images) [53]. As a result of the count, a polygon-type vector layer was obtained, in which each covers the entire surface of the canopy of a palm. Once the automatic

Detection of Diseases and Control
A sampling was performed with a navigating GPS at the same time as the first monitoring in order to understand the condition of certain palms with respect to diseases through the expertise of people from National Association of Oil Palm Cultivators (ANCUPA). This allowed to analyze their multi-temporal changes and has led to a validation of the detection of aforementioned affectations when conducting the first processing of the obtained orthomosaics.

Multitemporal Analysis
In the phase B for the change detection process, orthomosaics have been compared, taken at different times, in order to identify their differences. The difference has been calculated with the comparison of a specified input a selected index. A thresholding is applied, which allows establishing parameters that help the process to determine the areas that have large changes. In this case, to select an appropriate gray level threshold to extract objects from their background, the Otsu method [54] was applied, which employs the shape of the histogram. It is based on a discriminatory analysis and uses the cumulative zero and first order moments of the histogram to calculate the threshold level value We proceeded to integrate through a color-composition of multitemporal indexes that allowed to discriminate the progress of the affectations from the decrease in the values of the indexes.
Remote Sens. 2020, 12, 3229 9 of 21 Vegetation indices generated in the first monitoring were located in the red canyon and for the green canyon, the calculated index corresponding to the second monitoring. Thus, the areas where the spectral values have not presented considerable variations will remain uniform with respect to the previous scene. Finally, in the blue canyon, the change detection, has been entered, where the affected palms and the areas that have presented a decrease in the index, indicates a possible alteration and will stand out compared to the plants that do not have any affectation.

Ground-Truth
Once the analysis of the detection of the changes and multitemporal indices had been performed, a validation of the areas that tended to decrease the values of the vegetation indices, that is, areas of possible spread of the disease, was elaborated, with the aim to detect a distribution pattern. In this way, a random sampling was performed in these areas. The coordinates of these points were exported immediately, to be imported into the GPS navigator, together with the CIPAL staff, checking the palms within those areas and validating the growth of affected areas.

Phase A
For the first phase, four orthomosaics of different dates have been generated in order to obtain the baseline of palm crops without affectations. The indices that made it possible to differentiate the palm from its environment, in the case of the materials of the phenological state 305, INIAP, CIRAD, and ASD were the Difference Vegetation Index (DVI) and Modified Triangular Vegetation Index (MTVI), for the phenological state 301, the Taisha material was used for the NDVI, for Amazon and Unipalma the Burn area index (BAI). The automatic palm count ( Figure 5) allowed to analyze each palm belonging to each lot by species in order to generate zonal statistics for obtaining the baseline of healthy palms. Based on the zonal statistics, the characteristic indices of a healthy palm could be defined according to each type of species (Table 2). Hereby, we observed within this display that the vegetation indices calculated from green and infrared are affected by haze, weather conditions, and sunlight at the time of monitoring. Therefore, two conditions are specified, which in certain cases establish a baseline. Hereby, the first condition belongs to 30 May 2019 with greater luminosity and less cloudiness, while the second corresponds to 18 April 2019, with less presence of luminosity and greater cloudiness.  Based on the zonal statistics, the characteristic indices of a healthy palm could be defined according to each type of species (Table 2). Hereby, we observed within this display that the vegetation indices calculated from green and infrared are affected by haze, weather conditions, and sunlight at the time of monitoring. Therefore, two conditions are specified, which in certain cases establish a baseline. Hereby, the first condition belongs to 30 May 2019 with greater luminosity and less cloudiness, while the second corresponds to 18 April 2019, with less presence of luminosity and greater cloudiness. After conducting the analysis of variance, we obtained that for the phenological state 305, at a confidence level of 95%, there is statistical significance regarding the analysis of variance of the varieties INIAP, CIRAD, and ASD when using all of the indexes (ρ ≤ 0.05). Therefore, the alternative hypothesis test is accepted in all cases, which means that there are differences between the several varieties. This trend is reiterated in the two performed flights in phase A. Subsequently, Fisher's LSD and the Bonferroni test of each of the indexes used resulted in statistically differences between the INIAP, CIRAD, and ASD varieties, only for the RGNIR sensor and the two monitors. As for the BGNIR sensor, only the INIAP variety differs from CIRAD and ASD, but not CIRAD and ASD from each other. That means that all of the indices made with the RGNIR sensor is capable of detecting more rigorously the differences between species of the phenological state 305. On the other hand, the index of the BGNIR sensor allowed to differentiate between the Ecuadorian (INIAP, CIRAD, UNIPALMA, TAISHA) and Costa Rican (ASD and Amazon) species, but not the Costa Rican species. In this way, the feasibility of using the proposed vegetation indexes to differentiate palm species in the phenological state 305 was confirmed from the MAPIR RGNIR sensor.
Subsequently, the analysis of variance for the phenological status 301, corresponding to the approximate age of two years of the species (Taisha, Unipalma, and Amazon), was performed. In the case of the RGNIR, it was obtained when applying the vegetation indices to differentiate hybrid species. The use of geospatial technology is viable, since statistical significance (5%) is appreciated, related with the chlorophyll content of the plant. In the case of BGNIR, significant differences could be obtained between the UNIPALMA species with respect to the Taisha and Amazon species, but not between them.

Phase B
For the second phase, aerial surveys were conducted in the affected batch, which resulted in two orthomosaics (RGB and BGNIR) by date, which served as input for the analysis of the affectations of BR and RRD in the batches. In addition, VARI, GNDVI, and GVI vegetation indexes were generated for three monitoring ( Figure 6). In the three stages it was observed that the same range of index values is not always maintained. Despite the conducted atmospheric correction, the GVI and GNDVI indices are affected by other conditions, such as soil moisture and different environmental conditioners, while the VARI, for being an index that reduces atmospheric effects, responds better to the expected spectral response in a sick crop as it evolves. They did not remain constant, but it can be seen that the trend of the minimum indices representing affected palms maintained the trend, mainly in the areas of affected palms. orthomosaics (RGB and BGNIR) by date, which served as input for the analysis of the affectations of BR and RRD in the batches. In addition, VARI, GNDVI, and GVI vegetation indexes were generated for three monitoring (Figure 6). In the three stages it was observed that the same range of index values is not always maintained. Despite the conducted atmospheric correction, the GVI and GNDVI indices are affected by other conditions, such as soil moisture and different environmental conditioners, while the VARI, for being an index that reduces atmospheric effects, responds better to the expected spectral response in a sick crop as it evolves. They did not remain constant, but it can be seen that the trend of the minimum indices representing affected palms maintained the trend, mainly in the areas of affected palms.

Detection of Diseases and Control
Plants with disease detection are shown in Table 3. It can be verified that the monitoring points in the field coincide with palms that have low vegetation indices.

Detection of Diseases and Control
Plants with disease detection are shown in Table 3. It can be verified that the monitoring points in the field coincide with palms that have low vegetation indices. The points monitored in the field allowed to know that, indeed, the foci of palms affected by BR and RRD have indexes with low values in relation to the range managed by each of the calculated indices, which reveal affectations in the vigor of the plant. Next, the contrast between the RGB orthomosaic and the generated indices is detailed. Firstly, BR involvement has been considered, which demonstrated greater feasibility to be identified in the different states that are considered foci of infection (Figure 7a-f).

10
BR The points monitored in the field allowed to know that, indeed, the foci of palms affected by BR and RRD have indexes with low values in relation to the range managed by each of the calculated indices, which reveal affectations in the vigor of the plant. Next, the contrast between the RGB orthomosaic and the generated indices is detailed. Firstly, BR involvement has been considered, which demonstrated greater feasibility to be identified in the different states that are considered foci of infection (Figure 7a-f) (Figure 7e). With the GVI, the baseline is 5,124, grade 1 has 2.8 (Figure 7a), grade 2 2.5 (Figure 7b), grade 3 2 (Figure 7c), pc crater 1.8 (Figure 7d) and dead BR 2.1 (Figure 7e).
It is noticeable that in the case of BR crater severity, the index is much lower than a plant that is already dead from the damage. This type of problem can be caused by wilting, rotting and subsequent drying of the leaves, since the biomass of the plant is lost so that at the time of calculating the index values that correspond to weeds that have grown on the crop are interfering or from the ground itself. Regarding the involvement of RR, it is shown in all the indices generated, that the values are lower than the baseline of healthy palms. In the case of VARI, there are average values of 0.18, with GNDVI of 0.52 and with GVI of 3.4 (Figure 7f).
The results obtained in terms of disease detection indicate that the indices calculated from the RGB and BGNIR sensors of the MATRICE drone with the Zenmuse camera have the ability to differentiate between the degrees of severity, coinciding with previous studies where the ability of RGB and NIR sensors with RPAS technology have been able to differentiate the severity of the Ganoderma infection [55]. Nonetheless, the analyzes conducted in the current study were only carried out at the canopy level, without taking into account a study and comparison with the effects of the soil history of the study area. Therefore the need has been reached to determine these effects with sufficient precision in order to discriminate between the different categories of diseases and to develop methods of control and mitigation of the problems.

Change Detection
For the detection of changes between the first and second monitoring through the application of the VARI index. As a result we have the areas represented by color red in Figure 8, as the palms that presented a reduction in index values despite the fact that as mentioned above. The second monitoring showed an increase in value in the indexes due to atmospheric conditions. Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22

Change Detection
For the detection of changes between the first and second monitoring through the application of the VARI index. As a result we have the areas represented by color red in Figure 8, as the palms that presented a reduction in index values despite the fact that as mentioned above. The second monitoring showed an increase in value in the indexes due to atmospheric conditions. For the detection of changes between the second and third monitoring through the application of the VARI index+. Hereby, it is visible that the areas that demonstrate a change in decreasing indexes are more numerous. (Figure 9), which is consistent with the progress of the affectations present in the lot. The foci of palms with affectations present a greater coloration red, therefore it can For the detection of changes between the second and third monitoring through the application of the VARI index+. Hereby, it is visible that the areas that demonstrate a change in decreasing indexes are more numerous. (Figure 9), which is consistent with the progress of the affectations present in the lot. The foci of palms with affectations present a greater coloration red, therefore it can be deduced that the spatial distribution of the diseases of BR and RRD are distributed randomly around mentioned infection points. Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 Figure 9. Change detection VARI second and third monitoring.

Multitemporal Indexes
Regarding the interpretation of the results according to the methodology of multitemporal indexes, in the comparison of the first two monitoring the areas or crops that have a blue and violet coloration are sites where there has been a decrease in the values of the indexes, so in this case it coincides with the result obtained and described with the detection of changes for this pair of generated indices (Figure 10).

Multitemporal Indexes
Regarding the interpretation of the results according to the methodology of multitemporal indexes, in the comparison of the first two monitoring the areas or crops that have a blue and violet coloration are sites where there has been a decrease in the values of the indexes, so in this case it coincides with the result obtained and described with the detection of changes for this pair of generated indices ( Figure 10).  For the pair of multi-temporal indexes VARI 2 and VARI 3, a similar result was obtained to that from the detection of changes, by which the analysis carried out is contrasted. In this image the areas with a dark coloration are the ones that present changes in terms of decreasing values of the examined index ( Figure 11).
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 22 For the pair of multi-temporal indexes VARI 2 and VARI 3, a similar result was obtained to that from the detection of changes, by which the analysis carried out is contrasted. In this image the areas with a dark coloration are the ones that present changes in terms of decreasing values of the examined index ( Figure 11). Figure 11. VARI 2 and VARI 3 multitemporal indices.

Ground-Truth
In order to validate the multitemporal analysis, a total of 47 samples, which were the number of validation points that could be verified on the ground due to difficulties with the topography to access all, were examined from the results generated from areas with decreasing indexes. Based on the aforementioned, the VARI index has been considered as the basis to generate the sampling points. The coordinates of the contemplated palms were obtained and a field inspection was carried out in August 2019, as indicated with a summary in Table 4. All the monitoring presented indicated a mean below the established baseline. Of the 47 analyzed points, only two exposed another type of affectation that need to be considered when analyzing the vegetation indices. Such similar behavior, like the known Boron problem, may occur in the other points where around the correctly identified foci there are more plants distributed at random with affectation with all kinds of severity. Therefore, by means of the applied methodology it has not been possible to appreciate a specific spatial distribution pattern. On the contrary, a random advance of the problems was encountered, which is a dilemma that cannot yet be clarified, due to the lack of identification of the factors that give rise to the BR, and therefore the relationship it has with the RRD.
With the VARI index, it is evident that as the degree of severity of the affectation by BR progresses, the index decreases in each one of the conducted monitors, in addition to the fact that the RRD tends to have a lower value than the degrees of affectation two and three of the BR. This trend is confirmed in a certain way after the application of the GNDVI index. The base values for the prompt identification of foci of involvement of BR and RRD and of the degrees of severity that were possible to determine are detailed below, from all the generated indices (Table 5).

Ground-Truth
In order to validate the multitemporal analysis, a total of 47 samples, which were the number of validation points that could be verified on the ground due to difficulties with the topography to access all, were examined from the results generated from areas with decreasing indexes. Based on the aforementioned, the VARI index has been considered as the basis to generate the sampling points. The coordinates of the contemplated palms were obtained and a field inspection was carried out in August 2019, as indicated with a summary in Table 4. All the monitoring presented indicated a mean below the established baseline. Of the 47 analyzed points, only two exposed another type of affectation that need to be considered when analyzing the vegetation indices. Such similar behavior, like the known Boron problem, may occur in the other points where around the correctly identified foci there are more plants distributed at random with affectation with all kinds of severity. Therefore, by means of the applied methodology it has not been possible to appreciate a specific spatial distribution pattern. On the contrary, a random advance of the problems was encountered, which is a dilemma that cannot yet be clarified, due to the lack of identification of the factors that give rise to the BR, and therefore the relationship it has with the RRD.
With the VARI index, it is evident that as the degree of severity of the affectation by BR progresses, the index decreases in each one of the conducted monitors, in addition to the fact that the RRD tends to have a lower value than the degrees of affectation two and three of the BR. This trend is confirmed in a certain way after the application of the GNDVI index. The base values for the prompt identification of foci of involvement of BR and RRD and of the degrees of severity that were possible to determine are detailed below, from all the generated indices (Table 5).

Discussion
It has been able to determine the values for the baseline of palms without affection according to each proposed vegetation index. These values were obtained for each index. In the case of, the NDVI values for each material present at CIPAL plot vary between 0.737 and 0.683, the lowest values correspond to the phenological state 301. For the GNDVI index, two lines were obtained instead base of values according to the environmental conditions in which the monitoring was performed, since when grouping the data obtained, a normal distribution was not obtained. In the same way, the trend in the GVI index is maintained. Hereby the CIRAD species and ASD have two baselines. On the other hand, with the BGNIR sensor, two baselines were obtained according to the conditions for the GNDVI and GVI indices. In the same way that with what happened with the RGNIR sensor, a normal distribution could not be obtained by grouping the values of the two monitors an exception of the INIAP material, which did present distribution. When using the VARI index, similar results were obtained in both monitoring and the data could be grouped.
We used advanced methods of classification and discrimination of these categories in order to have as a result a significant difference between each one of them [56]. Statistical tests concluded that there are significant differences between species for the two sensors. With the help of Fisher's statistical LSD test and Bonferroni adjustment, it was possible to differentiate classes between palm varieties, which may be a great support to the palm regarding the alienation of INIAP species in the local market. Based on the statistical results, we encountered that with the BGNIR sensor it is feasible to differentiate the INIAP material from the others if we use the VARI index and with the RGNIR sensor instead, it was possible to draw differences between all species.
Regarding the phenological state 301, only UNIPALMA differs from the other materials with GNDVI and GVI, while with VARI it only differs from the Amazon material, in the first monitoring. On the second flight, UNIPALMA is differentiated from Taisha with all the proposed indices, therefore the feasibility to differentiate UNIPALMA is concluded.
Certain investigations confirm the feasibility of the proposal of this project when using the canyons (green, red and NIR) for the detection of diseases in plants [57]. In addition, the statistical test of the Lowest Significant Difference (LSD) was used to describe the difference between vigorous, affected and dead oaks. Hereby, it was possible to demonstrate that green, red and NIR have a statistically significant difference between all the analyzed categories. Finally, it is verified when comparing with the studies that the wavelengths of the green, red and infrared canyons have the aptitude to distinguish healthy and infected palms [58].
This result may be corroborated with previous studies, where an infected adult plant can intervene as a source of infection for three to five neighboring palms, which cannot necessarily be contiguous, therefore, they need to be strictly eliminated [19]. As for BR, the spatial distribution of the disease, the severity of the symptoms may vary considerably, even among palms that grow in proximity, without establishing an established spatial distribution pattern [11]. Therefore it is corroborated that the trend of detecting palms with problems is feasible.
The early detection of diseases according to the symptoms they present and the analysis of the spatial distribution of the diseases by means of a multi-temporal study were demonstrated. In the present case study, it has been established that in the boundaries of the lands with affectations there are a greater number of sources of affectations that harm the contiguous lands, as in the case of the Arteaga lot. In the particular case of the CIPAL lot, there are affectations close to the hybrid lots that are in early phenological stages, for which it is recommended to conduct a study in which it is disclosed that a sanitary strip of 200 m need to be established [39]. Thus, based on the Ecuadorian experience and measured from the hybrids palm cultivation edge, it may be applied to the local conditions of the study area.
Some multi-temporal studies in oil palm were previously performed, but used mostly satellite images and not RPAS, which provides data with higher spatial resolution [59,60]. In contrast, multi-temporal drone studies are generalized in forests [61] and in cultivation as floodplain [62]. On the other hand, in the present study, multitemporal studies were conducted in palm crops with RPASs, being innovative as far as known.
As mentioned in the results of phase A and B, there are technical limitations that affect the results obtained in the monitoring. In this case, the values of the indices do not indicate a certain uniformity. The limitations in the spectral characterization projects for the detection of the problem of Dry Arrow in oil palm in Costa Rica are technical, logistical and environmental, since in this case the digital levels of the orthomosaics obtained may have the noise of field observations, such as the water absorption peaks [63].

Conclusions
From the NDVI, GNDVI, GVI indices, statistically significant differences could be determined and the palm species corresponding to the phenological status of approximately fifteen years (INIAP, CIRAD, ASD) and those of the phenological state corresponding to two years (Taisha, UNIPALMA, Amazon).
In the detection of diseases, it is concluded that it was possible to appreciate that the palms have low values compared to the values of the established baseline. In addition to that visually with the rasters obtained from the proposed vegetation indices, palms with possible affectation can be detected. In the detection of diseases, it was possible to verify the feasibility of applying geospatial technologies by being able to determine foci of palms affected by BR and DRR, since these plants presented vegetation indices with lower levels than the established baseline.
In the present study, the methodology allows finding outbreaks of palms affected by RRD with higher incidence, since this disease does not present visible symptomatology at the beginning. This is a complication for its early detection, which is why it is recommended to implement methodologies that focus on the study of biomass lost by plants in all degrees of severity, in order to establish a baseline for early detection.
A methodology was proposed in order to conduct a multitemporal evaluation of the behavior of the affectations of BR and RRD in the study area, which were admissible to detect changes in the spread of the diseases studied. In the case of the detection of changes, areas were generated with decreasing values that allow knowing palms that have lowered the level of vigor compared to the previous monitoring. It was determined that the index that best fits to perform this method is the VARI. In order to corroborate the obtained results, the multitemporal index methodology was applied, which allowed appreciating similar results and strengthening the areas that were determined as progress of the problem. Therefore it was deduced that the distribution of the studied problem is of a random type around the identified sources of infection.
One of the main qualities of the VARI index is that it is based on the red, green, and blue bands, which are the typical bands of any RGB camera. Therefore, it is not necessary to use very expensive equipment in order to detect diseases. In case of using other sensors with infrared bands (RGNIR and BGNIR), the trend of detection of the disease is repeated.
Once the methodological proposal is established, it is recommended to add the characterization of spectral signatures of the affectations studied, based on multispectral information acquired in the field and laboratory. This serves as an input in order to estimate with greater solvency the biophysical characteristics of the culture, when adjusting the parameters of reflectance that best expresses the behavior of the plant health problem. Hereby, it was possible to characterize the spectral signature of the oil palm from satellite images [40]. These spectral signatures were not taken into account since they belong to a completely different study area than the one studied in this project. Furthermore, it is recommended that the adjustment of the model that characterizes the studied phenomena is supported with monitoring data in the field, in order to have time series that provide information for the forecasting models approach.
For the development of the methodological proposal, for the early detection of affectations of oil palm, it is recommended to conduct a phytosanitary census simultaneously with the application of geospatial technologies, in order to strengthen the results obtained for decision-making. In addition to the fact that it is proposed to perform the methodology in a more continuous way, such as weekly monitoring in order to analyze the results obtained with the vegetation indices in different climatic conditions, it should be possible to be able to better analyze the outliers obtained from the geodatabase generated and updated with each performed monitoring.
It is necessary to complement the proposed methodology with the execution of phytosanitary censuses in order to corroborate the results obtained from geospatial technology and direct decision-making towards the early detection of damages. Hereby, it enhances a great feasibility to detect foci. Therefore, it is necessary to develop a continuous monitoring system with palm growers, so that the proposed methodology has greater efficiency in early detection in palms precisely asymptomatic.