Maize Crop Coefﬁcient Estimation Based on Spectral Vegetation Indices and Vegetation Cover Fraction Derived from UAV-Based Multispectral Images

: Remote sensing-based crop monitoring has evolved unprecedentedly to supply multispectral imagery with high spatial-temporal resolution for the assessment of crop evapotranspiration (ETc). Several methodologies have shown a high correlation between the Vegetation Indices (VIs) and the crop coefﬁcient (Kc). This work analyzes the estimation of the crop coefﬁcient (Kc) as a spectral function of the product of two variables: VIs and green vegetation cover fraction ( f v ). Multispectral images from experimental maize plots were classiﬁed to separate pixels into three classes (vegetation, shade and soil) using the OBIA (Object Based Image Analysis) approach. Only vegetation pixels were used to estimate the VIs and f v variables. The spectral Kc fv :VI models were compared with Kc based on Cumulative Growing Degree Days (CGDD) (Kc-c GDD ). The maximum average values of Normalized Difference Vegetation Index (NDVI), WDRVI, amd EVI2 indices during the growing season were 0.77, 0.21, and 1.63, respectively. The results showed that the spectral Kc fv :VI model showed a strong linear correlation with Kc-c GDD (R 2 > 0.80). The model precision increases with plant densities, and the Kc fv :NDVI with 80,000 plants/ha had the best ﬁtting performance (R 2 = 0.94 and RMSE = 0.055). The results indicate that the use of spectral models to estimate Kc based on high spatial and temporal resolution UAV-images, using only green pixels to compute VI and f v crop variables, offers a powerful and simple tool for ETc assessment to support irrigation scheduling in agricultural areas.


Introduction
Crop water requirement (CWR) is the water needed to compensate for water vapour released into the atmosphere by evapotranspiration (ETc) and is one of the most critical variables in irrigated agriculture to better match water demand with irrigation supply [1]. There are several practical methodologies to estimate ETc. The most used method in the irrigation practice is well described in the FAO-56 manual [2], which is based on the knowledge of two variables: the reference evapotranspiration (ETo), representing the atmospheric evaporative demand, and the basal crop coefficient (Kcb), representing the water crop factors.
The ETo values represent the evaporative effects of environmental conditions on the evapotranspiration process and are commonly estimated by equations using weather data [3]. The Kcb coefficient is species-specific and varies according to the phenology, Green Vegetation Cover Fraction (f v ) has also been used for studying the phenological and physiological status of vegetation [28], estimating crop yields and monitoring crop developmental stages [29]. Marcial-Pablo et al. [30] showed that f v could be estimated with precision from VI indices from UAV-based RGB and multispectral images. Since VIs and f v are correlated with crop phenology and crop canopy density, it is necessary to study their correlation with Kc. Johnson and Trout [31] reported that Kc can be estimated based on f v indices with the support of vegetation indexes. Several authors reported a strong relation of f v and VIs with biophysical variables such as Fraction of Absorbed Photosynthetically Active Radiation (fAPAR), an important variable in the agricultural monitoring [32,33].
ETc estimation has focused on standardized crop conditions and full irrigation conditions [34]. There are very robust methodologies to estimate Kcb based on thermal time expressed as a continuous function of Cumulative Growing Degree Days (CGDD) reported by Ojeda-Bustamante et al. [6]. CGDD-based functions are considered a suitable option to account for the climate variability on crop growth and phenological development. However, Kc under non-standard conditions is important to know for crops with highly-variable water demands. For better water planning and management, and given the limited meteorological information in many agricultural zones, it is necessary to use geospatial technologies to indirectly estimate crop water demands. Remote sensing data have been used to assess water stress in non-standard conditions of various crops using low resolution such as satellites images and high resolution (UAV) imagery [18].
The existing research has been concentrated on using spectral data to estimate crop coefficient (Kc) based mainly on a linear relationship either with the vegetation index or with vegetation cover fraction. However, with the availability of high spatial and temporal resolution images and better image analysis algorithms, it is now better to split and classify pixels of vegetation, shadow, and soil and isolate their effects and result in better spectral models to fit experimental Kc data.
Therefore, the purpose of this paper was to estimate the maize crop coefficient (Kc) as a function of the product of spectral vegetation indices (VIs) and Vegetation Cover Fraction (f v ), using high-resolution multispectral UAV-images and an efficient OBIA-based approach to classify image pixels. The proposed methodology was validated using Kc, estimated with a Cumulative Growing Degree Days (CGDD) approach, and gathered experimental data with two planting densities.

Study Area and Site Data
The field study was conducted in the experimental station ( Figure 1) of the National Institute of Forestry, Agriculture and Livestock Research (INIFAP) in Zacatepec, in the Southern State of Morelos, Mexico (18 • 39 6.45" N and 99 • 11 59.63" W). The climate in the study area is hot sub-humid, with average annual temperature and precipitation of 24.3 • C and 892 mm, respectively; the rainy season is from mid-May to October, and the dry season, November to mid-May.
The experimental cropped area was 1 ha, sowed with maize hybrid H-515 on 5 July 2016. Twenty-four lots of 6 m × 25 m were planted with a furrows' separation of 0.8 m at two plant densities (60,000 and 80,000 plants/ha). The volumetric field capacity (FC) and permanent wilting point (PWP) soil water content were 41% and 22.5%, respectively. Weather daily data were available from a meteorological station located at 160 m of the study area (18 • 39 11.19 N and 99 • 12 00.47 S).

Crop Evapotranspiration (ETc)
Crop Evapotranspiration (ETc) is rarely measured in the field and is usually estimated by multiplying the reference evapotranspiration (ETo) with a crop coefficient (Kc) [1,2,6].

ETc = Kc × ETo
(1) Many equations have been developed to estimate ETo and extensively evaluated and compared with measured lysimeter ET under different climatic conditions. However, the FAO Penman--Monteith method is reliable and physically-based and adopted as the standard method for ETo computation [1,2]. With data registered in weather stations, the PM method estimates ETo as a function of net radiation, soil heat flux, wind speed, air temperatures, vapour pressure, and other coefficients as described in Allen et al. [1].
The crop coefficient Kc is a crop-specific coefficient representing the physical and physiological differences between a specific crop and the reference crop, considering the relationship between the atmosphere, crop physiology, and agricultural practices [1,2,15]. Kc values for annual crops increase from a minimum value at planting until reaching a maximum value at full canopy cover. One of the main limitations of this approach is to generate specific Kc values instead of using the general FAO Kc curves supplied by Allen et al. [1].

Estimation of GDD-Based Crop Coefficient (Kc)
Under water stress conditions and neglecting soil water evaporation, the crop coefficient (Kc) can be estimated using the following equation [1]: Kcb is the basal crop coefficient obtained under optimal crop development conditions, and Ks is the soil's water stress coefficient. Kcb was estimated using a GDD-based

Crop Evapotranspiration (ETc)
Crop Evapotranspiration (ETc) is rarely measured in the field and is usually estimated by multiplying the reference evapotranspiration (ETo) with a crop coefficient (Kc) [1,2,6].
Many equations have been developed to estimate ETo and extensively evaluated and compared with measured lysimeter ET under different climatic conditions. However, the FAO Penman-Monteith method is reliable and physically-based and adopted as the standard method for ETo computation [1,2]. With data registered in weather stations, the PM method estimates ETo as a function of net radiation, soil heat flux, wind speed, air temperatures, vapour pressure, and other coefficients as described in Allen et al. [1].
The crop coefficient Kc is a crop-specific coefficient representing the physical and physiological differences between a specific crop and the reference crop, considering the relationship between the atmosphere, crop physiology, and agricultural practices [1,2,15]. Kc values for annual crops increase from a minimum value at planting until reaching a maximum value at full canopy cover. One of the main limitations of this approach is to generate specific Kc values instead of using the general FAO Kc curves supplied by Allen et al. [1].

Estimation of GDD-Based Crop Coefficient (Kc)
Under water stress conditions and neglecting soil water evaporation, the crop coefficient (Kc) can be estimated using the following equation [1]: Kcb is the basal crop coefficient obtained under optimal crop development conditions, and Ks is the soil's water stress coefficient. Kcb was estimated using a GDD-based approach with Equation (3) proposed by Ojeda-Bustamante et al. [6], as a function of the cumulative Growing Degree Days (x n ) to day i.
Kc 0 is the basal crop coefficient for the first phenological stage that depends essentially on soil evaporation; K max is the maximum value of the basal crop coefficient (Kcb) during the growing period; erfc is the complement error function; x i is the cumulative Growing Degree Days (CGDD i ) elapsed until the day i from sowing or crop emergence, normalized to the maximum CGDD value required to complete the phenological cycle (CGDD x ), i.e., x i = CGDD i /CGDD x ; x Kmax is the normalized dimensionless value x when the maximum value K max is reached; α 1 is the regression coefficient obtained from fitting the model to the experimental data.
The Cumulative Growing Degree Days (CGDD) are calculated as the accumulated daily Degree Days (GDD) according to the following Equation (4) used by Ojeda-Bustamante et al. [6].
T a is the average daily air temperature, T c−min and T c−max are the basal (10 • C) and maximum (30 • C) growing temperatures for maize [35]. The stress adjusted crop coefficient of the day i (Ks i ) was determined as a function of the available soil water for day i (AW i ) as used by Jensen et al. [36]: Ojeda-Bustamante et al. [37] calibrated and validated the above equations for fully irrigated maize with a planting density of 95,000 plants/ha, obtaining the following parameters: K max = 1.25, Kc 0 = 0.2, X Kmax = 0.59. In our study, for plant densities of 60,000 and 80,000, the K max values were estimated at 1.03 and 1.10, respectively.
Daily values of Kc were estimated using the methodology and information validated for the study area as reported by Ojeda-Bustamante et al. [6] and adjusted under field conditions considering the experiment's planting densities.

UAV Multispectral Image Acquisition and Orthomosaic Generation
UAV multispectral images were acquired during the growing season as described as follows. Twelve Ground Control Points (GCPs) evenly distributed over the study area were placed and georeferenced with a TopCon model GR-5 GPS RTK with a vertical and horizontal accuracy of <1 cm. The hexacopter UAV used was a DJI A2 model (equipped with IMU) with a maximum load capacity of 2.5 kg and top ascent and descent speed of 6 m/s. The UAV was equipped with Tetracam ADC Snap (Tetracam Inc. Chatsworth, CA, USA) three-band multispectral camera (Table 1).
Two flights ( Figure 1c) were made on six different dates (from August to October, 2016, between 11:00 and 13:00 local time). The effective flight time for each date was 11 min and the average flight height was 50 m above the ground (GSD = 2.1 cm/pixel). Image overlap was 75% (front and side overlap) ( Table 2). The flights were planned and executed during the clear sky condition with the UgCS software (SPH Engineering, Riga, Latvia). RAW images acquired by the sensor are then radiometrically calibrated in postprocessing using PixelWrench2 software, thanks to a calibration target (spectralon) acquired before each flight. The digital elevation model and orthomosaic were generated using Pix4D software (Pix4D SA, Prilly, Switzerland) as shown in the Figure 2. Figure 3 shows as example, two multispectral orthomosaics obtained from photogrammetric processing of UAV images at the phenological stages. Two flights ( Figure 1c) were made on six different dates (from August to October, 2016, between 11:00 and 13:00 local time). The effective flight time for each date was 11 min and the average flight height was 50 m above the ground (GSD = 2.1 cm/pixel). Image overlap was 75% (front and side overlap) ( Table 2). The flights were planned and executed during the clear sky condition with the UgCS software (SPH Engineering, Riga, Latvia). RAW images acquired by the sensor are then radiometrically calibrated in post-processing using PixelWrench2 software, thanks to a calibration target (spectralon) acquired before each flight. The digital elevation model and orthomosaic were generated using Pix4D software (Pix4D SA, Prilly, Switzerland) as shown in the Figure 2. Figure 3 shows as example, two multispectral orthomosaics obtained from photogrammetric processing of UAV images at the phenological stages.

Image Classification and Estimation of GreEN Vegetation Cover Fraction (fv)
The classification of orthomosaics was performed using the OBIA approach. The OBIA detectors were implemented using eCognition developer 9.0 (Trimble GeoSpatial, Munich, Germany). OBIA combines spectral, contextual, and morphological information of the objects created in the segmentation [39] and minimizes the average heterogeneity of the objects in the image [40]. Several papers show the successful use of this methodology to analyze high-resolution spatial images taken with UAV [41][42][43].
Two stages are involved in an OBIA approach: image segmentation and classification. Orthomosaics was segmented using the multi-resolution segmentation (MRS) algorithm that use three criteria (scale, shape, and compactness) to generate the segments. The values of the final configuration for the scale parameter, shape, and compactness were: 25, 0.3, and 0.5, respectively, as proposed by Drǎguţ et al. [44] and Drǎguţ et al. [45]. The output of this step produced homogeneous segments of different sizes.

Image Classification and Estimation of GreEN Vegetation Cover Fraction (f v )
The classification of orthomosaics was performed using the OBIA approach. The OBIA detectors were implemented using eCognition developer 9.0 (Trimble GeoSpatial, Munich, Germany). OBIA combines spectral, contextual, and morphological information of the objects created in the segmentation [39] and minimizes the average heterogeneity of the objects in the image [40]. Several papers show the successful use of this methodology to analyze high-resolution spatial images taken with UAV [41][42][43].
Two stages are involved in an OBIA approach: image segmentation and classification. Orthomosaics was segmented using the multi-resolution segmentation (MRS) algorithm that use three criteria (scale, shape, and compactness) to generate the segments. The values of the final configuration for the scale parameter, shape, and compactness were: 25, 0.3, and 0.5, respectively, as proposed by Drǎguţ et al. [44] and Drǎguţ et al. [45]. The output of this step produced homogeneous segments of different sizes.
Finally, the classifications with three classes (vegetation, shadow, and soil) were based on the segments generated ( Figure 4). The method described by Marcial-Pablo et al. [30] was adapted to classify the orthomosaic into three classes.
Finally, the classifications with three classes (vegetation, shadow, and soil) were based on the segments generated ( Figure 4). The method described by Marcial-Pablo et al. [30] was adapted to classify the orthomosaic into three classes.

Estimation Vegetation Indices
There is a vast repository of vegetation indices used to explain vegetation characteristics related to ET. NDVI is probably the most widely used spectral index. However, in studies with satellite imagery, it can be affected by changes in soil background brightness and atmospheric effects caused by scattering and absorption by aerosols. In this study the use of UAV imagery can overcome these limitations. VIs such as EVI2 and WDRVI, improved versions of EVI and RDVI, offer better results by overcoming the absence of a blue band or the non-linear relationship with LAI or VF in satellite images. However, so far, these indices have been little studied in UAV images. Table 3 shows the VIs calculated for the six segmented multispectral orthomosaics using the QGIS software. VIs was estimated using two cases: using all pixels and using vegetation pixels only. A Tukey's mean test (p ≤ 0.10) was used to analyze if there were significant differences between VIs value obtained from Case 1 and Case 2. Note: α is weighting coefficient of 0.2 according to [47].
Increasing positive NDVI values indicate increasing amounts of green vegetation [49] because healthy vegetation has a low reflectance of red light and high reflectance of f v = Number of vegetation pixels/Number of (vegetation + non-vegetation) pixels (6)

Estimation Vegetation Indices
There is a vast repository of vegetation indices used to explain vegetation characteristics related to ET. NDVI is probably the most widely used spectral index. However, in studies with satellite imagery, it can be affected by changes in soil background brightness and atmospheric effects caused by scattering and absorption by aerosols. In this study the use of UAV imagery can overcome these limitations. VIs such as EVI2 and WDRVI, improved versions of EVI and RDVI, offer better results by overcoming the absence of a blue band or the non-linear relationship with LAI or VF in satellite images. However, so far, these indices have been little studied in UAV images. Table 3 shows the VIs calculated for the six segmented multispectral orthomosaics using the QGIS software. VIs was estimated using two cases: using all pixels and using vegetation pixels only. A Tukey's mean test (p ≤ 0.10) was used to analyze if there were significant differences between VIs value obtained from Case 1 and Case 2.  [48] Note: α is weighting coefficient of 0.2 according to [47].
Increasing positive NDVI values indicate increasing amounts of green vegetation [49] because healthy vegetation has a low reflectance of red light and high reflectance of near-infrared, which results in high NDVI values. NDVI is closely related to biophysical vegetation parameters, such as LAI, but has the drawback that it loses sensitivity to change when it reaches a threshold and becomes saturated. This saturation is because chlorophyll is a very efficient absorber of red radiation, and therefore, as plant biomass increases, red reflectance will not change. On the other hand, several spectral indices such as WDRVI and EVI2 are suitable in UAV imagery without the saturation limitation [50].
WDRVI shows greater sensitivity to changes in leaf area index (LAI) values from moderate to high (between 2 and 6) than the NDVI and maintains a linear relationship with the LAI of maize and soybeans [47,51].
EVI2 considers two bands of the enhanced vegetation index (EVI); it has been developed for sensors without a blue band. Besides, similar to EVI, it maintains improved sensitivity and linearity in high biomass and LAI regions [52]. EVI2 captures spatial-temporal variation in photosynthetically active vegetation [48] and is less prone to saturation in dense vegetation cover and less sensitive to different soil reflectances [53]. EVI2 is better than NDVI at discriminating land cover diversity [54] and capturing vegetation condition and structure changes [53]. However, the seasonal variation of EVI2 in arid and semi-arid zones is not as strong as in the NDVI [55]. The three spectral vegetation indices' values (VIs) were normalized to a scale of 0 to 1 for equally and better comparison.

Estimation Crop Coefficient Based f v :VIs (Kc fv:VI )
Several studies demonstrate that the use of spectral vegetation indices from nearinfrared (NIR) and visible regions has a linear correlation with the interception of sunlight from the canopy [56,57]. Kc obtained from vegetation indices are more robust than assuming standard conditions (Kcb) since they estimate current crop conditions [34]. Since high VI values indicate plant areas under optimal conditions, it was assumed that the model Kc = f(f v × VIs) is better than the simpler model Kc = f(VIs) used in most studies [11,34]. Therefore, Kc has been estimated as a function of two independent relationships: Kc = f(VIs) and Kc = f(f v ), as reported by Johnson and Trout [31]. Moreover, direct estimation of Kcb has been documented by Pôças et al. [15] using sophisticated equations as a function of f v and VIs. Vegetation Cover Fraction (f v ) describes a vertical projection of the areal proportion of a landscape occupied by green vegetation [58]. Grattan et al. [59] and Pôças et al. [15] have reported a strong correlation between f v and Kc.
Kc based on VIs and f v was estimated using the following equation: (green pixels only, using VIs and f v ): Kc fv:VIs = a + b * (f v * VIs) Vegetation pixels (7) where a and b are fitting parameters between CGDD-based Kc data and predicted data by the three models.

Evaluation of Model Performance
Predicted (Kc fv:VIs ) and observed Kc-c GDD values were compared using correlation analysis and linear regression (including 95% confidence prediction intervals). Several statistical parameters were calculated, mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), efficiency (E), and coefficient of variation (CV), to assess the performance of the regression model, using the equations (8 to 12) [60,61].
O i are the observed values (Kc-c GDD ), P i are the predicted values Kc fv:VI , O m is the mean observed value, and n is the total number of samples.

Maize Crop Coefficient CGDD-Based (Kc-c GDD )
The maize phenological duration was 119 days (from planting to physiological maturity), equivalent to 1655 GDDs. The Kc values are detailed in Table 4 and Figure 5, which shows that the peak Kc values are in the blister-milk stages.

Vegetation Indices (VIs)-CGDD Relationship
A Tukey HSD test (p ≤ 0.10) reflected that there are statistically significant differences between VIs values obtained from case 1 (considering three classes) and for case 2 (considering one class) (RMSE = 0.66; p = 0.00) mainly due to the effect when considering shade and soil pixels [62]. This difference is not detectable at larger spatial resolutions, as in satellite images, where pixel disaggregating is not possible. In that sense, case 2 was used in this work because it has better results and better represents the characteristics of the vegetation in the field. The simulated Kc values ( Figure 5), for both plan densities, show good agreement between the optimal (Kcb) and stressed (Kc) curve, indicating a good match between crop water demands and rainfall, with small differences during blister and milk stages. As expected, Kc values are higher for higher plant density (80,000 plants/ha) because of the higher vegetation coverage due to high crop biomass.

Vegetation Indices (VIs)-CGDD Relationship
A Tukey HSD test (p ≤ 0.10) reflected that there are statistically significant differences between VIs values obtained from case 1 (considering three classes) and for case 2 (considering one class) (RMSE = 0.66; p = 0.00) mainly due to the effect when considering shade and soil pixels [62]. This difference is not detectable at larger spatial resolutions, as in satellite images, where pixel disaggregating is not possible. In that sense, case 2 was used in this work because it has better results and better represents the characteristics of the vegetation in the field.
The The results are comparable to those reported by Calera et al. [64], who obtained NDVI values between 0.10 and 0.85 (with a peak at flowering) using Landsat 7 images. Also by Gitelson et al. [65], who found NDVI values from 0.28 to 0.85 for vegetative and reproductive stages, respectively, using MODIS images, and by Raun et al. [66] with peak NDVI values (0.77 to 0.73) at flowering, estimated with a Greenseeker sensor.
Most remarkable differences between NDVI, EVI2, and WDRVI occurred in the early and late stages, where EVI2 and WDRVI are more sensitive to vegetation changes than NDVI. In this regard, Guzinski [67] reports that although the three indices were useful for estimating crop phenology, there are differences between them when the crop canopy density changes drastically in the early and late stages.

fv:VIs Approach for Kc Estimation
The high spatial resolution of UAV imagery, as calculated for case 2, improves the The results are comparable to those reported by Calera et al. [64], who obtained NDVI values between 0.10 and 0.85 (with a peak at flowering) using Landsat 7 images. Also by Gitelson et al. [65], who found NDVI values from 0.28 to 0.85 for vegetative and reproductive stages, respectively, using MODIS images, and by Raun et al. [66] with peak NDVI values (0.77 to 0.73) at flowering, estimated with a Greenseeker sensor.
Most remarkable differences between NDVI, EVI2, and WDRVI occurred in the early and late stages, where EVI2 and WDRVI are more sensitive to vegetation changes than NDVI. In this regard, Guzinski [67] reports that although the three indices were useful for estimating crop phenology, there are differences between them when the crop canopy density changes drastically in the early and late stages.

f v :VIs Approach for Kc Estimation
The high spatial resolution of UAV imagery, as calculated for case 2, improves the accuracy of clustering, detection, and classification of pixel classes and thus increases the ability to recognize vegetation characteristics from vegetation indices. A strong linear correlation (R 2 > 0.8) between VIs from green vegetation pixels and Kc-c GDD values was found in the study. Therefore, considering both f v and NDVI improved the Kc(f v :NDVI) prediction, since f v is an expression of "only-green" biomass coverage and NDVI estimates the greenness "quality". Both variables are related with canopy photosynthetic activity, and non-green pixels should be not considered, since they are associated with shadow soil and decay leaves.
It is interesting to note that the fitting model accuracy (model vs observed Kc-c GDD data) increases with plant density due to high canopy density, improve the correlation for the case of NDVI (R 2 > 0.91), while they do not change for EVI2 and decrease for WDRVI (R 2 > 0.8). Furthermore, all VIs were less sensitive to detect Kc changes because of reduced canopy density in the final stages (maturity in 1598 CGDD). On the other hand, these results show that the f v : NDVI and f v : EVI2 models disaggregate the vegetation canopy cover from the soil more effectively than f v WDRVI at vegetative stages (12 leaves). Table 5 and Figure 7 shows the statistical parameters and the linear regression between Kc-c GDD values and VIs multiplied by f v value (f v :NDVI, f v :EVI2, and f v :WGRDI). A strong linear correlation with an R 2 coefficient greater than 0. 80 for all cases is observed. The model precision is high R 2 values about 0.94 (Figure 6), indicating that VIs values are sensitive to canopy density, associated with planting densities [68]. Besides, in early stages with low levels of canopy cover, the relationship of reflectance levels with Kc values is not very clear, due to higher variation on soil spectral properties [9]. On the other hand, the Kc fv:NDVI model (R 2 = 0.94) is more sensitive to vegetation cover change than Kc fv:EVI2 (R 2 = 0.80) and Kc fv:WDRVI (R 2 = 0.93) due to a reduction of the saturation of the NDVI when leaf cover reaches high levels.   Table 4.
Furthermore, all VIs were less sensitive to detect Kc changes because of reduced canopy density in the final stages (Maturity in 1598 CGDD). On the other hand, these results show that the fv: NDVI and fv: EVI2 models disaggregate the vegetation canopy cover from the soil more effectively than fv WDRVI at vegetative stages (12 leaves).
In this study, there are no images of the early growth stages of the crop because the canopy cover of the crop has a minimal effect on evapotranspiration; however, it could be interesting in future studies to monitor the crop periodically from the beginning of plant development.

Validation of Kcfv:VIs Models
The validation of Kcfv:VI models is based on statistical analysis of performance indicators shown on Figure 8 and Table 6, respectively. Our results reveal that the model that best estimates Kc is fv:NDVI since it showed lower RMSE (0.055-0.061), higher efficiency (E = 92-96%), and lower CV (7.78-6.55); fv:WRDVI and fv:EVI2 models presented RMSE > 0.064, E < 95% and CV > 6.97.   [4] calculated Kc values from MODIS imagery using a simple linear regression model Kc -NDVI with an R 2 of 0.83. Gitelson [47] recommended WDRVI to estimate canopy cover from Landsat images with higher accuracy than NDVI.
Furthermore, all VIs were less sensitive to detect Kc changes because of reduced canopy density in the final stages (Maturity in 1598 CGDD). On the other hand, these results show that the f v : NDVI and f v : EVI2 models disaggregate the vegetation canopy cover from the soil more effectively than f v WDRVI at vegetative stages (12 leaves).
In this study, there are no images of the early growth stages of the crop because the canopy cover of the crop has a minimal effect on evapotranspiration; however, it could be interesting in future studies to monitor the crop periodically from the beginning of plant development.

Validation of Kc fv:VIs Models
The validation of Kc fv:VI models is based on statistical analysis of performance indicators shown on Figure 8 and Table 6, respectively. Our results reveal that the model that best estimates Kc is f v :NDVI since it showed lower RMSE (0.055-0.061), higher efficiency (E = 92-96%), and lower CV (7.78-6.55); f v :WRDVI and f v :EVI2 models presented RMSE > 0.064, E < 95% and CV > 6.97. Agronomy 2021, 11, x FOR PEER REVIEW 15 of 20   All vegetation indices show a better performance in 80,000 plants/ha plots than 60,000 plants/ha plots. These results indicate that estimation of Kc in low-planting densities, using high spatial resolution images obtained from UAV, is good as long as f v is high and soil evaporation is negligible, as reported by Er-Raki et al. [70] and Farg et al. [14]. Furthermore, multispectral vegetation indices (NDVI, WRDVI, and EVI2) obtained from UAV images, and Kc values calculated (Kc-c GDD ) and predicted (Kc fv:VI ) follow the same trend through the different growth stages. f v :NDVI and f v :EVI2 are similar in performance parameters and linear regression validation, although f v :NDVI is slightly more accurate. These results indicate that when vegetation is isolated from the ground at a high spatial resolution, NDVI can overcome problems with high canopy density and lead to better results than EVI2 and WDRVI. Using generated models, the spatial distribution of crop Kc can be determined. Figure 9 shows, as an example, the maps of Kc for high density estimated for all pixels (case 1) and for only green pixels (case 2).  All vegetation indices show a better performance in 80,000 plants/ha plots than 60,000 plants/ha plots. These results indicate that estimation of Kc in low-planting densities, using high spatial resolution images obtained from UAV, is good as long as fv is high and soil evaporation is negligible, as reported by Er-Raki et al. [70] and Farg et al. [14]. Furthermore, multispectral vegetation indices (NDVI, WRDVI, and EVI2) obtained from UAV images, and Kc values calculated (Kc-cGDD) and predicted (Kcfv:VI) follow the same trend through the different growth stages. fv:NDVI and fv:EVI2 are similar in performance parameters and linear regression validation, although fv:NDVI is slightly more accurate. These results indicate that when vegetation is isolated from the ground at a high spatial resolution, NDVI can overcome problems with high canopy density and lead to better results than EVI2 and WDRVI. Using generated models, the spatial distribution of crop Kc can be determined. Figure 9 shows, as an example, the maps of Kc for high density estimated for all pixels (case 1) and for only green pixels (case 2).  Table 5 for the Kcfv:NDVI model. The model is the most accurate model, since the Kcfv:NDVI value is almost equal to the observed KcGDD = 1.069. In contrast, the average Kc values without considering fv for cases 1 and 2 are 1.157 and 1.049, respectively.
One additional advantage of the Kcfv:NDVI model is the it can be used to estimate areal Kc that can be applied to compute water needs in irrigation management units such as  Table 5 for the Kc fv:NDVI model. The model is the most accurate model, since the Kc fv:NDVI value is almost equal to the observed Kc GDD = 1.069. In contrast, the average Kc values without considering f v for cases 1 and 2 are 1.157 and 1.049, respectively.
One additional advantage of the Kc fv:NDVI model is the it can be used to estimate areal Kc that can be applied to compute water needs in irrigation management units such as furrow, border or basin in surface irrigation, as well a tree or row-crop in micro-irrigation.
The results indicate that the use of spectral models to estimate Kc based on high spatial and temporal resolution UAV-images, using only green pixels to compute VI and f v crop variables, offers an accurate and simple tool for ETc assessment to support irrigation scheduling in agricultural areas.

Conclusions
The results indicate that Kc prediction accuracy increases as Kc is expressed as a function of the VF multiplied by VI, both obtained from UAV multispectral images. Splitting and classifying pixels of vegetation, shadow, and soil, and isolating their effects were crucial to improve the model performance.
UAV-based images provided detailed spectral crop information to estimate Kc values throughout the phenological season. On the other hand, the GDD approach is a powerful Kc estimation tool to better adjust to phenology and local-scale environmental conditions. Kc values obtained from spectral data were validated with accurate GDD-based Kc values using experimental crop data with two planting densities.
This study indicates that three spectral-based Kc models showed good performance (R2 > 0.80). The model that best fit the Kc values was f v :NDVI (R2 = 0.91 to 0.94). f v :EVI2 and f v :WDRVI models showed lower performance, indicating that UAV images when vegetation is isolated from soil in high spatial resolution overcome the weaknesses of the NDVI index. It can reach saturation spectral values due to high canopy density or by changes in background conditions through soil brightness or surface soil moisture content. Since the Kc fv:NDVI model is not based in one pixel but grouped pixels in a cropped area, it can be used to estimate areal Kc which is more powerful for irrigation management purposes.
This study concludes that the use of spectral models for Kc estimation based on high spatial and temporal resolution reflectance values, once they are validated locally with experimental data, offers a powerful and simple tool for ETc estimation to support irrigation scheduling in agricultural areas.