A Novel ArcGIS Toolbox for Estimating Crop Water Demands by Integrating the Dual Crop Coefficient Approach with Multi-Satellite Imagery

Advances in information and communication technologies facilitate the application of complex models for optimizing agricultural water management. This paper presents an easy-to-use tool for determining crop water demands using the dual crop coefficient approach and remote sensing imagery. The model was developed using Python as a programming language and integrated into an ArcGIS (geographic information system) toolbox. Inputs consist of images from satellites Landsat 7 and 8, and Sentinel 2A, along with data for defining crop, weather, soil type, and irrigation system. The tool produces a spatial distribution map of the crop evapotranspiration estimates, assuming no water stress, which allows quantifying the water demand and its variability within an agricultural field with a spatial resolution of either 10 m (for Sentinel) or 30 m (for Landsat). The model was validated by comparing the estimated basal crop coefficients (Kcb) of lettuce and peach during an irrigation season with those tabulated as a reference for these crops. Good agreements between Kcb derived from both methods were obtained with a root mean squared error ranging from 0.01 to 0.02 for both crops, although certain underestimations were observed resulting from the uneven crop development in the field (percent bias of −4.74% and −1.80% for lettuce and peach, respectively). The developed tool can be incorporated into commercial decision support systems for irrigation scheduling and other applications that account for the water balance in agro-ecosystems. This tool is freely available upon request to the corresponding author.


Introduction
The important technological advances that have taken place in recent years allow for obtaining a greater predictive capacity and a more accurate assessment of different processes within agricultural systems [1].These advances have enabled the generation of models for the estimation of plant chemical and biophysical parameters for improving pest and disease management [2,3], nutrient assessment [4,5], irrigation scheduling [6,7], and yield predictions [8,9].Nevertheless, they frequently result in models that are difficult to apply to generate real solutions for the agricultural sector, as a consequence of limitations in public access to data, privacy and confidentiality issues, lack of funding for long-term collection and curation of data, and technical difficulties in managing and sharing data [10].
For the specific case of crop evapotranspiration (ET c ) estimation, several models have been developed following different approaches; the most common among them are those based on surface Water 2019, 11, 38 2 of 17 energy balance (SEB) and soil water balance (SWB).Most of these models are too complex for non-specialized users, as in the case of many farmers, and even by more qualified users such as technicians from irrigation communities or agricultural cooperatives.This low level of practical applicability is caused by the fact that these models are usually composed of a large number of variables and parameters, requiring, in many cases, long execution times [11,12].In addition, ET c is usually computed on a punctual location that may not be representative of the entire irrigated area.Using remote sensing technologies appears to be an efficient and economic approach to provide evapotranspiration (ET) estimations on large areas, allowing the assessment of the ET spatial variability [13].
The most common methodologies for ET estimation from remotely sensed imagery are those based on SEB and those based on vegetation indexes (VI).Models based on SEB take advantage of the relationship between surface reflectance data and energy balance components to derive ET c .Latent heat flux, and thus ET, is computed as the residual component of the SEB at the time of satellite overpass by subtracting the soil heat flux and sensible heat flux from the net radiation at the surface.In all SEB models, surface temperature data play a critical role in the estimation of crop water requirements.Nevertheless, spatial resolution improvements of satellites have occurred mainly in the visible and near infrared (VNIR) regions, whereas the thermal domain has remained at the same (or even lower) spatial resolution as that provided in the early 1980s [14].Thus, the resulting pixel size limits the use of SEB models, especially for heterogeneous crops.In addition, SEB approaches need to extrapolate the instantaneous ET (ET inst ) obtained at the time of the satellite pass to daily values, which can be usually done by assuming the conservation of the evaporative fraction or the ratio of ET to reference evapotranspiration, ET o [15,16].A recent study in Mexico showed that using the MOD16 algorithm for ET estimation, through an SEB model, on wheat crops and shrub lands had an average underestimation of 0.32 mm•day −1 (up to 1.48 mm•day −1 in certain cases) with respect to ET measurements in the field using eddy covariance equipment [17].
An alternative to SEB models is using VI-based approaches, which compute ET from empirical, site-specific relationships between the fraction of ground covered by vegetation and VI (for example, the normalized difference vegetation index (NDVI) or the soil-adjusted vegetation index (SAVI)).However, VI reflect the vegetation cover conditions of a given date for estimating the actual basal crop coefficient (K cb ), but require an additional approach, normally a water balance model, to estimate the evaporation coefficient (K e ).In addition, most VI-based methods assume baseline conditions (e.g., no soil water limitations) as being unable to consider the impact of soil water shortage, salinity, or disease on the estimation of the crop coefficient (K c ) [7,18].In this sense, previous reports indicated that, despite some drawbacks and an average error of 0.75 mm•day −1 , combining VI with the dual crop coefficient approach for estimating ET was a valid method for monitoring the spatial distribution of crop water requirements in several species including garlic, cotton, olive, mandarin, and peach [18].
Among all these approaches, the methodology proposed in the Food and Agricultural Organization of the United Nations (FAO) Irrigation and Drainage Paper No. 56 (FAO-56) [19] acquires special relevance for estimating water consumption by irrigated crops because of its worldwide use, especially in scientific research [20][21][22].This method has been criticized because of several sources of error that it may include, such as the consideration of a constant canopy resistance, the way in which some parameters are estimated, and the requirement of adapting the crop coefficient values to the local area in which the model is going to be applied, among others [23].Despite these drawbacks, in general, ET estimated using this procedure agrees relatively well with that measured with lysimeters for periods of bare soil and partial and full vegetation cover [24], as well as with ET measurements obtained through Bowen ratio procedures (a disagreement of 1.1% on average) [25].
Under this scenario, decision support systems (DSS) arise as tools intended for bridging the gap between an accurate ET estimation and its final application by farmers and technicians.Therefore, DSS facilitate using crop models to help decision-makers by reducing the time and human resources required for analyzing complex alternative decisions [26].Proof of the increasing importance of DSS is the number of them developed in the last decade, although most of them provide numerical results, neglecting the spatial variability of the considered variable [27][28][29].
In order to solve the above-mentioned limitations, the objective of this work was to develop a simple tool for assessing crop water demands and for providing advice in irrigation scheduling, thus generating a real solution for the agricultural sector.The tool was based on the joint application of the dual crop coefficient approach proposed in the FAO-56 document [19] and multi-satellite images that permit incorporating the spatial component into the estimations using a geographic information system (GIS) tool.

Model Theorical Framework
The dual K c approach estimates the evapotranspiration of a given crop (ET c ) using the basal coefficient (K cb ), the evaporation coefficient (K e ), and the evapotranspiration of a reference surface (ET o ), defined as the water consumed by an extensive surface of green, well-watered grass of uniform height, actively growing and completely shading the ground [19].
In this framework, K cb is the ratio of ET c over ET o when the soil surface is dry but transpiration is occurring at a potential rate, that is, water is not limiting transpiration; whereas K e refers to the soil evaporative component.The sum of both coefficients is known as the crop coefficient (K c ).
General K cb values for the different seasonal stages of the main crops are tabulated [19], but refer to sub-humid climates with moderate wind speeds (e.g., minimum relative humidity of 45% and wind speeds of 2 m•s −1 , approximately).For those regions that do not meet these conditions, tabulated K cb must be adjusted using the meteorological variables observed at those locations for correctly representing the local conditions [19].
In order to incorporate the spatial component into the estimation of K c , the approach proposed in this study includes data derived from remote sensing imagery.More specifically, the model calculates two VI from bottom of atmosphere reflectivity values derived from satellite images: NDVI and SAVI.The first index, NDVI, allows for estimating the soil fraction covered by vegetation (F c ) through the linear, or approximately linear, relation between some VI and F c in the range from bare soil to near full ground cover [30,31].In addition, F c in combination with SAVI, allows for correcting K cb when non-homogeneous crops are considered [32].
where the max and min subscripts refer to a high leaf area index (LAI) and a bare surface, respectively; and F c,max is the vegetated covered fraction for which K cb reaches its maximum value (K cb,max ).In this study, K cb,max values were the FAO-56 tabulated values adjusted (K cb,corr ) for the local relative humidity and wind speed conditions (K cb,corr in Figure 1) [19].
For the estimation of K cb , the days when satellite images are not available (e.g., cloudy conditions or days without scheduled pass), the model extrapolates from the last available satellite image.Thus, for initial and mid stages, K cb is assumed constant and, therefore, K cb of a given day is considered equal to the K cb value of the last available image.For development and final stages, the slope between the K cb of the last available image and K cb tabulated for mid and final stages, respectively, was determined.
In addition, K e depends on the water available in the soil surface layer, which is computed as Water 2019, 11, 38 4 of 17 where K c,max represents the maximum K c value after an irrigation or precipitation event, and is calculated using the wind speed (u 2 , m•s −1 ), minimum relative humidity (RH min , %), and crop height (h, m), as follows: K r = (TEW − D e )/(TEW − REW) (5) where TEW and REW are the total and the readily evaporable water, respectively (mm); and D e represents the cumulative depletion from the soil surface layer (mm) estimated from a daily soil water balance.
Kr = (TEW − De)/(TEW − REW) where TEW and REW are the total and the readily evaporable water, respectively (mm); and De represents the cumulative depletion from the soil surface layer (mm) estimated from a daily soil water balance.

Model Implementation
The dual Kc model was implemented using ArcPy [33] because of its capacity for geographical data analysis, conversion and management, and map automation.This environment utilizes Python as a programming language.The model was split into two sub-models.The first one is intended for VIs, Fc, and Kcb computation, whereas the second one is focused on Ke, Kc, and ETc estimation (Figure 1).Finally, the sub-model formulations were incorporated into an ArcGIS (v10.2;Esri, Redlands, CA, USA [33]) toolbox, which provides a friendly interface to execute the model, where the user only needs to select the values for the required inputs.The created tool was programmed for using images acquired from multiple satellites, identifying the source by the file name.This tool can be freely distributed to interested users by contacting the corresponding author.

Model Validation
The described model was validated in two sites: (i) a 1.32 ha lettuce orchard (Lactuca sativa L. cv.Little gem) located in Los Alcáceres, Murcia, southeast Spain (37°46′7′′ N, 0°50′13′′ W) during 2016; and (ii) a 3.56 ha peach orchard (Prunus persica L. Bastch cv.Amandina) located in Molina de Segura, Murcia, Southeast Spain (38°06′52′′ N, 1°12′46′′ W) during 2017.The lettuce planting date was 24 September and harvest took place on 11 November; therefore, the length of the crop growing cycle was 48 days.On the other hand, the peach orchard was planted in 2013 and the fruit harvest for 2017 occurred in May, followed by a long post-harvest period (namely, the peach cultivar is early maturing).Cloudless images from satellites Sentinel 2A, Landsat 7 and Landsat 8 collected during the studied growing seasons were considered as inputs.The spatial resolution of the images in the VNIR region was 10 m for Sentinel 2 and 30 m for Landsat 7 and 8, whereas the temporal resolution was 10 days for Sentinel 2A and 16 days for Landsat 7 and Landsat 8. Details about the different images used in this study are shown in Tables 1 and 2. A total of 8 and 23 cloudless satellite images were used for characterizing the lettuce and peach Kcb curves, respectively.They corresponded with the 73% (lettuce) and 45% (peach) of all satellite images acquired during the considered period.Only pure pixels, namely, those totally included in the plots, were considered for determining crop water Inputs Intermediate Results

Model Implementation
The dual K c model was implemented using ArcPy [33] because of its capacity for geographical data analysis, conversion and management, and map automation.This environment utilizes Python as a programming language.The model was split into two sub-models.The first one is intended for VIs, F c , and K cb computation, whereas the second one is focused on K e , K c , and ET c estimation (Figure 1).Finally, the sub-model formulations were incorporated into an ArcGIS (v10.2;Esri, Redlands, CA, USA [33]) toolbox, which provides a friendly interface to execute the model, where the user only needs to select the values for the required inputs.The created tool was programmed for using images acquired from multiple satellites, identifying the source by the file name.This tool can be freely distributed to interested users by contacting the corresponding author.

Model Validation
The described model was validated in two sites: (i) a 1.32 ha lettuce orchard (Lactuca sativa L. cv.Little gem) located in Los Alcáceres, Murcia, southeast Spain (37 • 46 7 N, 0 • 50 13 W) during 2016; and (ii) a 3.56 ha peach orchard (Prunus persica L. Bastch cv.Amandina) located in Molina de Segura, Murcia, Southeast Spain (38 • 06 52 N, 1 • 12 46 W) during 2017.The lettuce planting date was 24 September and harvest took place on 11 November; therefore, the length of the crop growing cycle was 48 days.On the other hand, the peach orchard was planted in 2013 and the fruit harvest for 2017 occurred in May, followed by a long post-harvest period (namely, the peach cultivar is early maturing).Cloudless images from satellites Sentinel 2A, Landsat 7 and Landsat 8 collected during the studied growing seasons were considered as inputs.The spatial resolution of the images in the VNIR region was 10 m for Sentinel 2 and 30 m for Landsat 7 and 8, whereas the temporal resolution was 10 days for Sentinel 2A and 16 days for Landsat 7 and Landsat 8. Details about the different images used in this study  1 and 2. A total of 8 and 23 cloudless satellite images were used for characterizing the lettuce and peach K cb curves, respectively.They corresponded with the 73% (lettuce) and 45% (peach) of all satellite images acquired during the considered period.Only pure pixels, namely, those totally included in the plots, were considered for determining crop water demands in order to avoid edge-pixel contamination.In this sense, for the lettuce plot, the number of pure pixels was 6 and 99 for Landsat and Sentinel, respectively; whereas the number of pure pixels for the peach field was 23 and 291 for Landsat and Sentinel, respectively.
During the growing seasons, the plants were cultivated using standard practices in the area with a cumulative irrigation application of 120 mm and 293 mm for lettuce and peach, respectively.The total amounts of N, P 2 O 5 , and K 2 O applied to fulfill the estimated crop nutrient demand were 56, 30, and 76 kg•ha −1 , respectively, for lettuce and 105, 85 and 192kg•ha −1 , respectively, for peach.
Lettuce crop height was measured weekly in 60 plants randomly distributed within the plot (Table 3), whereas the average peach height measured in the field was 2.40 m at the beginning of the year and was considered to be constant during the season.This information, as well as tabulated K cb values [19], was incorporated into the model (h and K cb (Tab) in Figure 1).The tabulated K cb values for initial, mid, and end stages were 0.15, 0.90, and 0.90, respectively, for lettuce and 0.45, 0.85, and 0.60, respectively, for peach.Daily and hourly weather data, including air temperature, relative humidity, wind speed, precipitation, reference evapotranspiration (calculated using the Penman-Monteith equation [19]), and solar radiation (Figures 2 and 3), were obtained from the nearest weather stations.The main utility of using hourly data is that this allows for a precise weather characterization at the satellite acquisition time.For the lettuce study site, the weather station considered was located 3 km from the experimental site at San Javier, Murcia (code TP22; http://siam.imida.es),whereas that for the peach site was located 1.5 km from the study site at Molina de Segura, Murcia (code MO22; http://siam.imida.es).Mean air temperature and relative humidity for the study period were 20.7 • C and 76.2%, and 18.57 • C and 54.23% for lettuce and peach, respectively; with a cumulative precipitation and ET o of 30.2 and 120.9 mm, respectively, for lettuce, and 140.9 and 1002.2 mm, respectively, for peach (Figures 2 and 3).Noticeably, only a few rainfall events greater than 5 mm occurred over both growing seasons and accounted for more than 80% of the total rainfall within the studied periods (Figures 2 and 3).The values of TEW and REW were selected according to those defined in the FAO-56 document [19] for the different soils in both sites.More specifically, lettuce was grown on a loamy soil (sand 42%, silt 37; clay 21%), whereas the peach orchard was planted on a soil that was clay-loamy (sand 40%, silt 17%, clay 43%); corresponding with REW and TEW values of 9.0 and 19.0 mm, respectively, for lettuce and 9.5 and 24.5 mm, respectively, for peach.In addition, soil water balance was initialized at 15 September (day of year, DOY, 259) for lettuce and 15 January (DOY 15) for peach, because no significant precipitation events occurred during the previous weeks, allowing the assumption that all evaporable water had been depleted from the soil surface layer at the beginning of calculations, that is, D e, 0 = TEW.

Statistical Analyses
The comparison between Kcb values derived from the ArcGIS approach proposed in this study and those values tabulated in FAO-56 [19] was performed over the validation sites using the root mean square error (RMSE).

RMSE
where Kcb, FAO56 is the Kcb value derived from FAO-56; Kcb, estimated is the Kcb value estimated with the created tool, and L is the length of the crop season.
Considering that no other direct ETc determination was available in the field sites, the tabulated FAO-56 values can be considered as a reasonable proxy for ET determinations because of their widespread use, particularly in herbaceous crops.The spatial variability within each plot was quantified by means of the coefficient of variation (CV) among the values of all pixels included in the studied fields:

Statistical Analyses
The comparison between K cb values derived from the ArcGIS approach proposed in this study and those values tabulated in FAO-56 [19] was performed over the validation sites using the root mean square error (RMSE).
where K cb, FAO56 is the K cb value derived from FAO-56; K cb, estimated is the K cb value estimated with the created tool, and L is the length of the crop season.
Considering that no other direct ET c determination was available in the field sites, the tabulated FAO-56 values can be considered as a reasonable proxy for ET determinations because of their widespread use, particularly in herbaceous crops.The spatial variability within each plot was quantified by means of the coefficient of variation (CV) among the values of all pixels included in the studied fields: CV = SD × 100 K cb, mean (7) where SD is the standard deviation within the plot and K cb, mean is the average K cb of all pixels included in the plot.
Additionally, the percent bias, PBIAS, was calculated to measure the average trend of the simulated data to be larger or smaller than the FAO-56 derived K cb values.
The optimal value of PBIAS is 0.0, with low-magnitude values indicating accurate model simulations.Positive values indicate overestimation, whereas negative values indicate model underestimation.

Model Interface
Figure 4 shows the final interfaces for the different sub-models developed.On the left-hand side of each interface, the inputs required for running the model are listed.Thus, the required inputs for sub-model 1 and sub-model 2 are shown in Table 4.To help the user, a detailed description for each input is provided on the right-hand side of the interface when the mousepad icon is placed on each input box.
Water 2018, 10, x FOR PEER REVIEW 9 of 17 where SD is the standard deviation within the plot and Kcb, mean is the average Kcb of all pixels included in the plot.Additionally, the percent bias, PBIAS, was calculated to measure the average trend of the simulated data to be larger or smaller than the FAO-56 derived Kcb values.
The optimal value of PBIAS is 0.0, with low-magnitude values indicating accurate model simulations.Positive values indicate overestimation, whereas negative values indicate model underestimation.

Model Interface
Figure 4 shows the final interfaces for the different sub-models developed.On the left-hand side of each interface, the inputs required for running the model are listed.Thus, the required inputs for sub-model 1 and sub-model 2 are shown in Table 4.To help the user, a detailed description for each input is provided on the right-hand side of the interface when the mousepad icon is placed on each input box.

Practical Validation
Regarding the model validation, Figure 5 shows the temporal evolution of the K cb values derived from the ArcGIS approach proposed in this study, as well as those values tabulated in FAO-56 [19], over the lettuce (Figure 5A) and peach growing seasons (Figure 5B).A good agreement between both K cb values was obtained, as the RMSE (0.01-0.02) was low for both crops.The K cb values estimated for lettuce were 0.15 and 0.85 for the initial and the mid-season stage, respectively, and values intermediate between these two for the crop-development stage.However, a slight underestimation (2.8% on average for the entire growing season) was observed in lettuce for the K cb values obtained from the ArcGIS approach when compared with those derived from FAO-56, especially during the crop development stage (5.7% on average for the development stage; Figure 5A).For the peach orchard, K cb values estimated were 0.43 and 0.84 for the initial and the mid-season stage, respectively, and values intermediate between these two for the crop-development stage.A slight underestimation was observed in peach for the K cb values obtained from the ArcGIS approach when compared with those derived from FAO-56, with this being more evident during the initial and mid-season crop development stages (2.6% on average; Figure 5B).The reason for the lower underestimation observed in peach when compared with lettuce could be the higher frequency of available images used to build the K cb curve (Figure 5).The PBIAS values obtained for lettuce and peach were −4.74% and −1.80%, respectively, confirming the underestimation pattern observed in Figure 5.
The K cb variability within the lettuce field, assessed through CV, was low at the beginning of the growing season (CV ≈ 7.4%), increased during crop-development (CV ≈ 12.2%), and decreased again at the mid-season (CV ≈ 7.7%).For the peach orchard, CV was relatively high at the initial stage (CV ≈ 11.0%), and decreased during crop-development and mid-season stages (CV ≈ 4.7%).This variability could have been induced by the fact that our approach takes into account the actual soil surface covered by the crop (as determined by satellite images) for calculating K cb , reflecting the uneven distribution of plant size for lettuce and the initial bud break for peach within the studied plots at the crop development and mid-season stages.
Examples of the spatially distributed K cb output resulting from the application of the model for lettuce and peach at the development stage and at the end of the growing cycle (day of year (DOY) 282 and 312, respectively, for lettuce; and 56 and 156, respectively, for peach) are shown in Figure 5.The mean K cb values on the lettuce plot were around 0.20 and 0.85 for DOY 282 and 312, respectively; whereas the mean K cb values on the peach orchard were 0.50 and 0.84 for DOY 56 and 156, respectively.The reductions of K cb near the field edges in both crops are the result of the pixel contamination due to the surrounding bare soil.This suggests that a larger pixel cut off could have reduced the errors.
The Kcb variability within the lettuce field, assessed through CV, was low at the beginning of the growing season (CV ≈ 7.4%), increased during crop-development (CV ≈ 12.2%), and decreased again at the mid-season (CV ≈ 7.7%).For the peach orchard, CV was relatively high at the initial stage (CV ≈ 11.0%), and decreased during crop-development and mid-season stages (CV ≈ 4.7%).This variability could have been induced by the fact that our approach takes into account the actual soil surface covered by the crop (as determined by satellite images) for calculating Kcb, reflecting the uneven distribution of plant size for lettuce and the initial bud break for peach within the studied plots at the crop development and mid-season stages.
Examples of the spatially distributed Kcb output resulting from the application of the model for lettuce and peach at the development stage and at the end of the growing cycle (day of year (DOY) 282 and 312, respectively, for lettuce; and 56 and 156, respectively, for peach) are shown in Figure 5.The mean Kcb values on the lettuce plot were around 0.20 and 0.85 for DOY 282 and 312, respectively; whereas the mean Kcb values on the peach orchard were 0.50 and 0.84 for DOY 56 and 156, respectively.The reductions of Kcb near the field edges in both crops are the result of the pixel contamination due to the surrounding bare soil.This suggests that a larger pixel cut off could have reduced the errors.

Discussion
One main advantage of the tool developed here is that it does not require many inputs, and those compulsory for running the model are readily available, making its use easier.This solves the over-parameterization issue, which is one of the main limitations for using models for practical applications in agriculture [36,37].In fact, after a certain level, an increase in model structure complexity does not necessarily lead to better estimations [38,39].In this sense, the tool proposed here was able to provide reliable estimations of crop water demands using a set of easily available parameters; this capability, along with its friendly interface, might favor the use of the tool by farmers and technicians.In its present form, the toolbox is implemented using ArcGIS because of the flexibility offered by this commercial software.In the future, other software developers might take advantage of the algorithm details described here to migrate the toolbox into other open-software systems such as Quantum Geographical Information System (QGIS), expanding the operational possibilities of the tool developed here.
In addition, the integration of remote sensing technology into the tool permits obtaining a spatial representation of the crop water demand variability over a given plot or orchard and, therefore, it will be possible to apply a different irrigation scheduling to a given sector within the field [40][41][42][43].This is the main improvement of our tool when compared with other approaches.Indeed, remote sensing technology has been incorporated into some DSS platforms (System of Participatory Information, Decision-support, and Expert knowledge for River-basin management, SPIDER [27]; Automated Radiative Transfer Models Operator, ARTMO [44]; AquaGIS [45]; Online Professional Irrigation Scheduling; Flexible and PrecIse IrriGation PlAtform to Improve FaRm Scale Water PrOductivity, FIGARO [46]).Some of these platforms are still ongoing developments, while some others provide information about other aspects of plant physiology and not directly on crop water requirements.To date, AquaGIS has been employed for simulating multi-year wheat yields under several irrigation strategies [45], proving the applicability of this approach for managing irrigation in large areas.However, AquaGIS is based on AquaCrop, a model that requires a higher number of inputs and parameters than the one implemented in the toolbox described in the current study.Therefore, our approach reduced model complexity when compared with AquaGIS.Other DSS, such as ARTMO [44], simulate and estimate vegetation structure or radiative transfer from vegetal canopies, with no regards to irrigation scheduling and, therefore, they have to be used by experts in order to gather the information provided by these platforms and translate it into irrigation strategies to be applied in the field.Another example of DSS with a module for irrigation scheduling is FIGARO [46], which also uses AquaCrop for determining crop water requirements, and has been applied only to field crops such as cotton [47].In summary, our approach is easier to use than current platforms because it is more intuitive and requires a lower number of parameters and inputs; moreover, it is easily adapted to a great number of crop species, as long as their crop coefficients are available, as exemplified for lettuce (herbaceous plant) and peach (woody perennial species) in the current manuscript.
Furthermore, most of the available applications intended for crop management do not consider the spatial component, providing only numerical outputs and neglecting the spatial variability within the field [8,28,29,[48][49][50].Additionally, the tool proposed in this study allows for integrating information obtained from different satellites; improving the reliability and consistency of vegetation state variables; and, consequently, enhancing the estimations of crop water requirements [51][52][53][54].Furthermore, the inclusion of the tool into a geographic information system (GIS) allows for task automation, which reduces computational time and permits modelling large areas at once, which is an important limitation in other models [11,12].
The validation with data for growing seasons of lettuce and peach in southeast Spain indicates that the proposed tool provides accurate crop water demand estimates when the FAO-56 methodology is taken as reference [19].In this context, K cb for lettuce is estimated at the crop development and mid-season stages in the current study as being from 0.6 to 0.85, which leads to K c values similar to those measured in southeast Spain using Bowen ratio procedures, although slightly lower at the initial stage of crop development [25,55].Regarding the peach orchard, K cb values ranged from 0.45 (at the initial stage) to 0.85 (at the mid-season).The values at the initial stage were considerably higher than those proposed by other authors [56,57], whereas estimated K cb values for mid stage were intermediate when compared with those previously reported [56,57].This highlights the importance of the K cb database used as input, which can clearly affect the model efficiency in estimating the actual ET c .The period considered here covered the initial and mid stages, being necessary to enlarge the study until the leaf fall period to reach the K cb value of the end season stage.It is necessary to develop and implement further new model routines, other than semi-empiric models, for computing ET, to be able to compute plant transpiration more specifically based on specific canopy resistance models.Additionally, the application of the proposed methodology allows the detection of the timing for bud break, complete canopy growth, and leaf senescence and fall.Nevertheless, a limitation with this type of crop is given by the spatial resolution of the satellite images, which does not allow distinguishing among orchard floor management practices.Moreover, outputs generated from the proposed tool collect K cb spatial variations within the field, and this information can allow for precise irrigation management [58][59][60]; for instance, if the field has several irrigation sectors, they can be differently managed according to a possible uneven distribution of crop water requirements.In this sense, the uneven development of the crops within the field was reflected by the slight underestimation observed in the K cb values derived from the tool, mainly during the crop development stage, which could have modified the SAVI and, therefore, the K cb value.In contrast, the FAO-56 method employs a fixed value for fraction cover at a given developmental stage of the crop, irrespective of the SAVI, likely causing the underestimation of our approach when compared with the FAO-56 as a reference.Moreover, the FAO-56 K cb values used compute crop water requirements under conditions of complete absence of biotic or abiotic stress and we cannot rule out the possibility that the experimental plants had temporarily suffered some slight limitations in their growth and development.In order to consider this issue, the next version of the developed tool will include the water stress coefficient, K s [19].Additionally, the temporal evolution of CV and RMSE showed that the variation observed was probably caused by the uneven distribution of plant sizes that left a variable surface of bare soil within the field.Once the crop started to develop, CV and RMSE increased, reflecting different growing rates of lettuce and peach canopies throughout the field.Finally, at mid-season, the variability decreased, probably because the crops reached a high-vegetated fraction cover and a more homogeneous distribution (for the case of lettuce) over the entire field.Furthermore, the selection of the pixels within the study area and the way the satellite images are cut off may cause additional errors for the estimation of crop water requirements at the edges of the field.
The implemented tool could be incorporated into a cloud-based DSS that allows for integrating data from various sources and, therefore, for assessing different scales ranging from small, single farms to large catchments; for centralizing the control of field devices; and for transforming the data into useful knowledge for decision making [61,62].In this sense, potential applications might consider the estimation of the soil water content, which could be, for instance, of interest in humid areas where high soil moisture levels often hamper soil management using heavy tractors and machinery.The tool created can be also used in combination with weather forecasts as inputs instead of past weather data, which may allow for estimating future irrigation requirements and predicting the occurrence of extreme events [63][64][65].Additionally, the satellite image download, pre-processing (e.g., radiometric, geometric, and atmospheric image corrections; cloud detection), and post-processing tasks are susceptible to being automatized, so that the results could be delivered to end-users (farmers, stakeholders, and policy makers) at near-real time, which facilitates the adoption of quick responses to agricultural issues [66].
Indeed, the proposed tool provides a user-friendly approach for the spatial estimation of crop water demands that can be incorporated into a DSS to generate real solutions for the agricultural sector in terms of irrigation scheduling and other applications related to the agro-ecosystem water balance.

Figure 1 .
Figure 1.Flowchart of the dual Kc model implemented into ArcPy, showing the inputs, intermediate results, and outputs.ρNIR and ρRED are the reflectances corresponding to the infrared and red bands, respectively; fw is the fraction of soil surface wetted by irrigation; and I is the irrigation depth on each day that infiltrates the soil (mm).For definitions of the rest of abbreviations used in the figure, readers are referred to the main text.

Figure 1 .
Figure 1.Flowchart of the dual K c model implemented into ArcPy, showing the inputs, intermediate results, and outputs.ρ NIR and ρ RED are the reflectances corresponding to the infrared and red bands, respectively; f w is the fraction of soil surface wetted by irrigation; and I is the irrigation depth on each day that infiltrates the soil (mm).For definitions of the rest of abbreviations used in the figure, readers are referred to the main text.

Figure 2 .
Figure 2. Daily weather data obtained from the weather station located at San Javier, Murcia (code TP22; http://siam.imida.es)for the period including the lettuce growing season in 2016.Abbreviations: Prec = precipitation; T mean = mean temperature; T max = maximum temperature; T min = minimum temperature; RH = relative humidity; u 2 = wind speed; ET o = reference evapotranspiration; SR = solar radiation.

Figure 3 .
Figure 3. Daily weather data obtained from the weather station located at Molina de Segura, Murcia (code MO22; http://siam.imida.es)for the period including the peach growing season in 2017.Abbreviations: Prec = precipitation; T mean = mean temperature; T max = maximum temperature; T min = minimum temperature; RH = relative humidity; u 2 = wind speed; ET o = reference evapotranspiration; SR = solar radiation.

Figure 4 .Figure 4 .
Figure 4. Toolbox interfaces of sub-model 1 (left) and sub-model 2 (right) implemented in ArcGIS (geographic information system), showing the different input parameters (Table4provides explanations for each item).RefET-reference ET; DSA-digital elevation model, slope, and terrain aspect; SAVI-soil-adjusted vegetation index.

Figure 5 .
Figure 5. Temporal evolution of the K cb values obtained from the developed ArcGIS tool and those derived from the Food and Agricultural Organization of the United Nations (FAO)-56 document for lettuce ((A): Coordinates 37 • 46 7 N, 0 • 50 13 W) and peach (B): Coordinates 38 • 06 52 N, 1 • 12 46 W); and examples of the resultant K cb images for different crop stages.Error bars represent the standard error (SE) within the field.)

Table 1 .
Satellite imagery used for model validation during the lettuce growing season.

Table 2 .
Satellite imagery used for model validation during the peach growing season.

Table 3 .
Temporal evolution of lettuce height at the study site (average from 60 measurements per date and standard error, SE).

Table 4 .
Description of the different inputs considered in each sub-model.FAO-Food and Agricultural Organization of the United Nations; RefET-reference evapotranspiration; SAVI-soil-adjusted vegetation index; LAI-leaf area index.