USING REMOTELY PILOTED AIRCRAFT (RPA) IMAGERY TO MAP THE PROFITABILITY OF COTTON CROPS

ABSTRACT The goal of this paper was to apply a method for the delineation of vegetation homogeneous zones (HZs) to map the profitability of cotton ( Gossypium sp.) using NDVI (Normalized Difference Vegetation Index) data from RPA imagery. Two irrigated cotton crops, at the FAU (103 ha) and FRP farms (106 ha), located in the western region of Bahia State, Brazil, were studied. The NDVI images were classified into three HZs: “high”, “medium” and “low” plant vigor; using the k-means clustering method. In each HZ the yield was measured, and the profitability estimated. In FAU we found that the lower the HZ’s NDVI values, the lower its profitability, because there was profit in the “high vigor” HZ, of US$ 829.40 ha-1, and loss in the “low vigor” HZ, of about US$ ‒1.256.09 ha-1 (‒251.44% compared with the “high vigor”). In the FRP plantation, the lower the NDVI values, the lower the loss of the plantation, as there were losses in all HZs: from “high” to “low vigor” HZ, the percentage difference in profitability was ‒343.43%. Thus, the use of a low-cost modified near-infrared Canon S100 on an RPA enabled the mapping of crop profitability, aiding the search for the factors behind yield variability.


INTRODUCTION
Optical vegetation indices (VIs) are mathematical models, applied to remote sensing measurements, that were developed to represent vegetation properties such as photosynthetic activity, pigment content, canopy structure, and others.After a systematic review of the scientific literature, Zeng et al. (2022) proposed a classification of VIs into three main ecological traits: biophysical, biochemical, and physiological.The VIs associated with biophysical plant parameters are those related to parameters such as the leaf area index and green biomass.In this regard, one of the most used VIs is the NDVI (Normalized Difference Vegetation Index), which is an index developed to represent the structure of canopies.
In precision agriculture, VIs such as the NDVI have already been used in agricultural applications, mainly in delineating management zones (MZs), as performed by Breunig et al. (2020), Oliveira Maia et al. (2022), Georgi et al. (2018), Kuiawski et al. (2017), Romano et al. (2021), and Scudiero et al. (2018).MZs are a basic concept for most approaches in precision agriculture (Breunig et al., 2020), and are characterized by areas of the field where plant productivity and soil attributes are similar in space and time.Among the advantages of delineating the MZs is the possibility of applying site-specific management per management unit, aiming to maximize profit and yield and minimize production costs and the impact of the activity on the environment.
The best variables for delineating MZs should be spatially correlated with crop yield (Gavioli et al., 2016).Yield maps, however, are not always readily available, so other methods for delineating MZs have been proposed.Such methods include the fusion of yield and soil attribute maps (such as soil apparent electrical conductivity) or Engenharia Agrícola, Jaboticabal, v.43, n.3, e20220218, 2023 delineating MZs using only soil attribute maps, while other methods are based on the application of VIs (Breunig et al., 2020), which are founded on the correlation of crop yield and biophysical, biochemical, and/or physiological properties of the plants.
MZs delineation is commonly applied with the use of data from more than one crop and when the other variables used in the delineation are considered stable in time and space (such as soil attributes), as can be observed in the works of Breunig et al. (2020), Georgi et al. (2018), and Maestrini & Basso (2018).When information from plants of the same harvest is used to delineate zones, these zones of similar plant behavior are commonly called homogeneous zones (HZs) (Marino & Alvino, 2018;Romano et al., 2021).
The use of VIs for delineation of MZs or HZs and their use for site-specific management configure a practical application of data surveyed by remote sensing platforms, from orbiting satellites to the aerial survey level using remotely piloted aircraft (RPAs).In fact, the use of RPAs for crop monitoring has become increasingly frequent (Chen et al., 2021;Marino & Alvino, 2018;Rasmussen et al., 2021), since they allow a greater level of detail than satellite images.
In crops like cotton, in which financial investments, agrochemical application, and mechanization are intensive, the use of technologies for continuous monitoring is fundamental to ensure the high yield and profitability of the crop.Thus, in this work, we applied a method for the delineation of HZs to estimate and map the profitability of two cotton crops (Gossypium sp.) using NDVI images obtained using an RPA.

Aerial survey of cotton fields using an RPA
Two central pivot-irrigated farms in the municipality of Barreiras, in the western region of Bahia State, Brazil, were used in this study (FIGURE 1).The first farm, Pivot 23 of the FAU farm, has approximately 103 ha of planted area, while the second farm, Pivot 01 of the FRP farm, has approximately 106 ha.Coupled to the ISIS RPA was a low-cost multispectral camera used to collect the images.A consumer-grade point-and-shot Canon S100 modified camera was used to collect light in the green, red and nearinfrared (NIR) bands.The modified Canon S100 collected the three green-red-NIR bands from 520 to 880 nm.The sensor was used because of its low acquisition cost, as well as being widely used in plant evaluation studies and vegetation vigor zoning (Barrows & Bulanon, 2017;Duarte-Carvajalino et al., 2018;Haghighattalab et al., 2017;Jewan et al., 2022;Ogliari et al., 2022;Van der Merwe & Price, 2015), including phenotyping (Haghighattalab et al., 2016).
The Canon S100, with a GPS receiver attached, captured frames at a resolution of 12.1 MP (effective pixels), which corresponds to 4000 × 3000 pixels per frame.With a focal length of 5.2 mm for all shots, it was possible to obtain a GSD (ground sample distance), or spatial resolution, of 5 cm (pixel size) at the flight height of 120 m.The flights were planned for image captures with an image overlap of 75% (overlap in the flight direction) and 65% (side overlap).

Delineation of homogeneous zones (HZs) with NDVI images
The images obtained in each plantation with the RPA were processed in Pix4Dmapper software.This step consisted of the generation of orthomosaics, from georeferencing, orthorectification, and mosaicking of the obtained images.Then the NDVI was calculated for each orthomosaic and the areas of interest (AOI) were cut out, a step in which only the pixels of the crops were kept.
The orthomosaics were generated with a spatial resolution of 0.11 m, which tends to be a very fine resolution for the purposes of applications in agriculture such as delineating HZs.Georgi et al. (2018) reached this conclusion by analyzing RapidEye images, of 6.5 m spatial resolution.Thus, the NDVI images underwent postprocessing in order to attenuate their detailing level.
The post-processing consisted of filling in missing data, i.e., possible pixels with no-data value had their values interpolated by averaging the values of neighboring pixels (in a 3 × 3 kernel, with a maximum search of 9 × 9 pixels); resampling, to reduce the spatial resolution, configuring an output resolution of 1.0 m and employing the averaging resampling method; then the images were smoothed with an averaging filter with a 51 × 51 pixel moving window (kernel).This post-processing was adapted from the methodology of Georgi et al. (2018), and although resampling and smoothing have a strong impact on the level of detail in the images, the focus of this work was on delineating geometric zones, and high levels of detail tend to produce small-scale detail, which may be too much for practical implementation in agricultural machinery.
After post-processing, the images were classified using the k-means clustering method, where three classes (number of clusters) of NDVI were segmented.An odd number of classes was chosen to ensure an "average" class of crop behavior (Georgi et al., 2018), and a smaller number of classes (such as three) facilitates an overall analysis of the crops (Albornoz et al., 2018;Song et al., 2009), in this case, an analysis of the crop's profitability.
The post-processing and classification of the NDVI images were performed with Python programming language resources (version 3.7).Subsequently, the delineation of homogeneous zones was performed using the QGIS software (version 3.22.10),where the NDVI classes were converted into polygons.All the processing for obtaining the NDVI-based HZs is presented in FIGURE 3.
FIGURE 3. Post-processing of the NDVI images and delineation of the homogeneous zones (HZs) in the cotton fields.
After defining the HZs, the NDVI distribution of each HZ was submitted to statistical similarity tests.For this, the non-parametric Kruskal-Wallis test was applied to test the overall statistical difference between the three groups of NDVI (coming from the three HZs), followed by Dunn's post-hoc test (also non-parametric), where pairs of NDVI distributions were tested.Both non-parametric tests were chosen because the NDVI distributions did not meet the requirement of homogeneity of variances to be submitted to the parametric test (ANOVA).

Estimating crop yield and profitability
Once the HZs were delineated, the yields (CY) in each of them were calculated.To obtain the CY, the harvest was carried out, in a sampled area, followed by weighing, to obtain the harvested production (HP) and the harvested area (HA) per HZ.With the CY per HZ, the profitability (P) of each HZ was computed using [eq.( 1 The first and second terms of [eq.( 1)] refer to the revenue obtained from the sale of the production harvested (selling of cotton lint and cottonseed, respectively), while the third term (CPC) refers to the unit cost of cotton crop production in the Western region of the State of Bahia (for the 2018/2019 crop season).The values adopted for CLP and CSP were R$ 5.79 kg -1 and R$ 500.18Mg -1 , respectively, while the CPC value adopted was about R$ 10,106.20 ha -1 .All data were obtained for the study region, from the 2018/2019 crop season, through the Brazilian Confederation of Agriculture and Livestock (CNA, 2019) and the Associação Brasileira dos Produtores de Algodão (the Brazilian Association of Cotton Producers) (ABRAPA, 2022).The CLP, CSP, and CPC values were converted from Engenharia Agrícola, Jaboticabal, v.43, n.3, e20220218, 2023 R$ to US$ using the 2019 average exchange rate of 3.9445 R$/US$, provided by the Banco Central do Brasil (IPEA, 2022).So, the values employed in [eq.( 1)] were:  = 1.47 $  , CSP = 126.80$  , and  = 2,562.10$ ℎ .
Considering that the studied crops did not receive any site-specific management, that is, they received the same treatment (of irrigation, fertilizers, etc.), the expectation was that the yield and profitability should perform uniformly.However, due to the interaction of production factors and the spatial variability of the yield, this homogeneity of productivity and profitability was not evident.Thus, in addition to estimating CY(HZ) and P(HZ), the loss of profitability was also calculated.The loss of profitability was calculated using [eq.( 2)], where the difference between the P of the high vigor HZ and the P of the other HZs is calculated.

RESULTS AND DISCUSSION
The NDVI of the crops, after image resampling from 0.11 to 1.0 m, is presented in FIGURE 4. In both cotton plots, we have identified visually regions with different cotton canopy biomass behavior.After image classification, the HZs were delineated and are displayed in FIGURE 4b) for FAU and FIGURE 4e) for FRP.For the FAU plantation, the "high vigor" HZ was delineated as 47.40 ha, the "medium vigor" HZ as 44.20 ha, and the "low vigor" HZ as 11.55 ha.For the FRP plantation, the "high vigor" HZ was delineated with 50.06 ha, the "medium vigor" HZ with 36.62 ha, and the "low vigor" HZ with 21.06 ha    et al., 2018;Gavioli et al., 2016;Ohana-Levi et al., 2021;Oldoni & Bassoi, 2016;Schenatto et al., 2017;Song et al., 2009).In these works, there is a consensus that the feasible limit of MZs per plot is approximately six classes.Methods for automatic definition of the optimum number of classes use indices such as the FPI (fuzziness performance index) to indicate the degree of dissimilarity between delineated classes (Schenatto et al., 2017).This type of methodology was not used in this work because we sought to represent the general behavior of crops concerning their economic return, and a larger number of classes would make it difficult to obtain the crop's cash yield in each HZ since even smaller areas would be delineated.Therefore, the number of HZs was set to three, ensuring the detection of an average behavior class.Indeed, employing methods such as the FPI, Song et al. (2009) likewise detected an optimal number of three HZs also using vegetation indices.
In each delineated HZ, the cotton harvest was done to obtain an estimate of the crop yield (CY), as presented in TABLE 2. Crossing the delineation of the HZs, presented in FIGURE 4, with the CY for each HZ, we have noticed that the lower the NDVI (see boxplots in FIGURE 4c) and FIGURE 4f)), the lower the CY.2017) developed a study to monitor the yield variability of soybean (Glycine max (L.) Merrill) and black oat (Avena strigosa Screb.) using NDVI images obtained by an RPA.The authors concluded that the NDVI images were efficient in monitoring crop yield since the crop cover (or canopy) is correlated with its yield.This correlation was also found by Breunig et al. (2020), who studied rapeseed (Raphanys sativus L.), white oat (Avena sativa L.), and rye (Secale cereale L.) crops.Although we did not apply correlation tests in the methodology from the data presented in FIGURE 4 and TABLE 2, we could infer that the NDVI of cotton also tends to decrease when its yield decreases, and this finding allowed us to map the profit and loss of profit of the crops based on the crop yield of each HZ.
The result of this economic profitability analysis, i.e., the calculation of profit and loss of profit in the crops, is presented in FIGURE 5. FIGURE 5a) and FIGURE 5b) display the profitability and loss of profitability, respectively, of the FAU plantation, while FIGURE 5c) and FIGURE 5d) show the profitability and loss of profitability, respectively, of the FRP plantation.In the FAU plantation, there are two HZs where there is profit from the activity, in the "high and medium vigor" HZs.However, of the 103-ha planted, there are losses in approximately 11 ha, in the "low vigor" HZ.In this region of the farm, the interaction of the genotype and environmental factors did not favor the maximum productivity of the plants.In the FRP plantation, there was no profit from the activity.In all the HZs, we calculated losses in cotton production; however, the degree of loss is different for each HZ.In FRP, the lower the NDVI (ranked), the greater the loss.
The loss of profitability (FIGURE 5b) and FIGURE 5d)) represents the negative difference between the highest profitability achieved in the plantation (P of the "high vigor" HZ) and the profitability of the other HZs.For the FAU plantation, two HZs showed a profit, and the third HZ showed a loss.Although the inputs and water were applied aiming for uniformity of crop yield throughout the area, obtaining different results, in this situation, configures a loss of yield and, consequently, profitability.
In these situations, if the patterns are repeated for successive crops, it is necessary to investigate the causal agent(s) of this limitation, be it the occurrence and interaction of pathogens with the plants, the carry-over of fertilizers along with soil losses, or even unevenness in the irrigation depth, and make a decision, which could be to reanalyze the agricultural suitability of the crop in the existing edaphic conditions, or to perform specific local management of the crop, seeking to optimize the production cost and cash yield.Unlike FAU, in FRP the whole crop would need site-specific management to obtain profits.
Throughout the literature, there is a consensus on the importance of the yield map for knowledge of the spatial variability of crops (Momin et al., 2019;Sanches et al., 2019), since it represents the response of plants to the conditions of the production factors during cultivation.However, yield monitors are not always available in existing harvesters on farms, and this may be a hindrance to the use of precision farming practices, since there is a close relationship between plant yield and vegetation indices (Breunig et al., 2020;Oliveira Maia et al., 2022;Ferrer et al., 2020;Georgi et al., 2018;Romano et al., 2021), which aim to represent their biophysical or biochemical parameters.
Operational and net profits tend to increase with the employment of yield maps (Schimmelpfennig, 2016), so we cannot neglect the existence of on-board technologies in the most recent harvesters available, such as yield monitors, which, after correct calibration of the system and appropriate data treatment (which requires personnel with expertise), are used to generate the so-called yield maps.However, not all farms, even in developed countries such as the United States of America, have access to or use these technologies, as shown in the USDA (U.S.Department of Agriculture) report published by Schimmelpfennig (2016).The data in that report also show that the adoption of these precision farming technologies is even lower in the cotton Engenharia Agrícola, Jaboticabal, v.43, n.3, e20220218, 2023 crop.Therefore, the use of technologies such as RPAs -and data processing techniques generated by these platformsis an alternative for these farms, since the producer does not necessarily need to own the equipment required for aerial mapping and expertise for data processing, since specialized companies can be hired.
Furthermore, both yield and vegetation index maps do not show the possible causes of the spatial variability of plant yield, which can be associated with (for example) spatial variability of soil attributes.In this case, the use of RPA images should be employed not as an end, but as a means for searching for these causes, which should be investigated by obtaining and evaluating other variables of the agricultural systems environment.This subject is widely explored in literature reviews such as those by Maes & Steppe (2019) and Aslan et al. (2022).

CONCLUSIONS
In the FAU plantation, we verified losses to the farmer in the "low vigor" homogeneous zone.As for the FRP plantation, we verified losses throughout the crop, although in different degrees: the lower the NDVI values, the greater the estimated loss.Thus, the methodology employed with the use of remotely piloted aircraft allows the farmer to map the areas with profitability and loss, and provides evidence for the adoption of site-specific management in the plots, such as the search for the causes of spatial variability of crop yield.

FIGURE 1 .
FIGURE 1. Location of the studied cotton fields.The field data survey was conducted in the 2018/2019 crop season: cotton was sown on January 25 and February 02, 2019, on the FAU and FRP farms, respectively.The aerial survey of both crops was conducted with the crop already

FIGURE 4 .
FIGURE 4. Study areas: FAU, in subplot a) the NDVI map, in b) the homogeneous zones (HZs), and in c) the distribution of NDVI values for each HZ; for FRP d) shows its NDVI map, b) its homogeneous zones, and c) the distribution of NDVI values in each HZ.
FIGURE 4c) and FIGURE4f) display the NDVI distribution for FAU and FRP, respectively, in each HZ.The NDVI distributions of HZs are statistically different (P = 0.05) in each crop, as can be seen in TABLE1.

FIGURE 5 .
FIGURE 5. Profitability (P) and loss of P from cotton planting in each homogeneous zone (HZ) in a) and b) for FAU plantation and in c) and d) for FRP plantation.

TABLE 1 .
Results of the Kruskal-Wallis and Dunn's tests displaying the statistical difference of NDVI values from each HZ in each crop area. is the chi-squared Kruskal-Wallis test statistics, df is the degrees of freedom, Z-stat is the statistics of Dunn's pairwise test, and * indicates significant statistical difference between groups (for the Kruskal-Wallis test) and between pairwise groups (for Dunn's test).

TABLE 2 .
Obtained harvested area (HA), harvested production (HP), and crop yield (CY) for each homogeneous zone (HZ) from FAU and FRP farms.