**1. Introduction**

As a result of the environmental impact factors and the growing demands on food production and consumption, combined with the global market demand to keep merchandise prices low, the modern agricultural industry faces a major challenge [1,2]. There is greater urgency than ever before for producers, farmers, and agronomists around the world

**Citation:** Alexandris, S.; Psomiadis, E.; Proutsos, N.; Philippopoulos, P.; Charalampopoulos, I.; Kakaletris, G.; Papoutsi, E.-M.; Vassilakis, S.; Paraskevopoulos, A. Integrating Drone Technology into an Innovative Agrometeorological Methodology for the Precise and Real-Time Estimation of Crop Water Requirements. *Hydrology* **2021**, *8*, 131. https:// doi.org/10.3390/hydrology8030131

Academic Editors: Aristoteles Tegos and Nikolaos Malamos

Received: 18 June 2021 Accepted: 25 August 2021 Published: 1 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to improve the management of their agricultural practices on farms in response to the reduction in their budgets through the optimization in water-saving and farm inputs and, at the same time, maintain high product quality [3]. Future water requirements, combined with limited or bad quality water resources and requirements to adequately meet a growing population's increasing nutritional needs, require improved crop yields and productivity and increased irrigated land, even in the planet's arid areas.

Globally irrigated agriculture delivers 40% of the food production and consumes 70% of the available water [1,4]. Therefore, irrigation scheduling is an important tool for improving the efficient use of irrigation water. Recently, several new sensors and emerging technologies, such as the Internet of things (IoT), have been developed and applied in precision agriculture for the management of available water resources and the biometeorological monitoring of plants and soil [2,5]. These provide significant potential in precision agriculture (PA) and smart farming since it has a direct impact on improving the management of irrigation systems and enabling a long-term increase in productivity [6]. Nowadays, the scientific community focuses not only on estimating the duration of irrigation but also on the accurate calculation of actual needs (actual evapotranspiration (Eta)) regarding the water required by each crop [7]. This trend is contrary to the general tendency of determining the maximum rate of evapotranspiration of the crop (ETC) in large estimation time steps. Additionally, the effect of the vegetation surface (temperature, physical characteristics, etc.) on the precise estimation of actual evapotranspiration values is also one of the main issues of research [7].

Earth observation data with variations in spectral, spatial, and temporal characteristics have been widely used for vegetation mapping and crop water stress monitoring [8–11]. Traditional remote sensing methods place remote sensors over crop fields or use aircraft and satellites where the appropriate fixed position or the temporal and spatial resolution significantly bounds their utility for PA [12–14]. However, freely available optical-sensorsbased imagery is often unsuitable for monitoring crop water stress at the farm level due to the poor revisiting times and coarse spatial resolutions. At the same time, the higher accuracy data are too expensive [15,16]. Furthermore, these data are typically unavailable or not useful on cloudy days [2,17]. Nevertheless, for decades, a lot of studies were undertaken to monitor the crop water stress index using satellite- and aerial-based instruments in several crops, such as those of Rud et al. [18] and Cucho-Padin et al. [19] in potato fields, Veysi et al. [20] and Lebourgeois et al. [21] regarding sugarcane, da Silva et al. [22] regarding melon, Sepulcre-Cantó et al. [23] in olive orchards, and Gutiérrez et al. [24] regarding vineyards.

The synergistic utilization of innovative UAV advanced high-precision thermal cameras and other infrared sensors have enhanced the usefulness of these systems to monitor water statuses [25,26]. As a result, numerous studies have researched several different crops to evaluate the crop water stress conditions, such as those of Zhang et al. [16,27], Yang et al. [28], Gago et al. [12], Gonzalo-Dugo et al. [29], Bellvert et al. [30], Berni et al. [31], Matese et al. [32], and Santesteban et al. [33]. However, the process of calibrating and processing thermal images takes a long time to give the final decision on whether to irrigate. In addition, this process involves more empiricism. Therefore, even a high-precision thermal camera could not be described as a useful tool for everyday use on extensive arable land. The CWSI can be determined using at least two different methodologies: the theoretical model proposed by Jackson et al. [34], which is based on meteorological models, and the empirical model proposed by Idso et al. [35], which has achieved considerably more popularity due to the limited data requirements [8,16]. Due to its simplicity and requiring only three variables to be measured, this practical approach has received much attention in the literature. However, it received criticism concerning its inability to account for temperature changes due to radiation and wind speed [36]. The theoretical method is more complicated because it requires these two additional variables to be measured to evaluate the aerodynamic resistance. Moreover, given the assumed net radiation and wind speed, the theoretical approach to calculating the CWSI promises to improve the

estimation of plant water trends [36]. Considering the above, the methodology for estimating the CWSI coefficients is based mainly on the theoretical approach. For this purpose, all the necessary sensor equipment was installed for the CWSI calibration stage for three crops independently.

The GreenWaterDrone (GWD) project proposed a pilot implementation of an innovative and autonomous system (onboard micrometeorological station) that identifies a farmer's real-time irrigation needs through the CWSI to estimate each crop's water status. The direct estimation of the CWSI is based on spatial data that is remote and directly from the field and is based mainly on the infrared temperature of the crop canopy. Even in the recent (2021) literature, but also in previous years, the calculation of the CWSI with direct spatial measurements of the temperature of the crop canopy has never been proposed or adapted [37–43]. Such measurements are constantly achieved indirectly through thermal images. The raw data and thermal images from the infrared camera were also used for comparison and calibration purposes. The use of the spatial measurements of canopy temperature, air temperature, and relative humidity from sensors incorporated into an unmanned aerial vehicle (UAV) may be a global novelty. This study aimed to present the new methodology and the equipment used in the assessment of crop water stress by presenting preliminary results from the initial stages of a pilot implementation in a potato cultivation field at the farm scale.

#### **2. Description of the GWD Concept**

#### *2.1. CWSI—Empirical and Theoretical Approaches*

CWSI is a measure of the relative transpiration rate that occurs from a crop and is found by using the canopy temperature and the vapor pressure deficit as relevant variables. The latter measure is related to the atmospheric dryness over the crop. The CWSI approach [34] utilizes the energy balance theory that separates the net radiation over the canopy from the sensible heat (thermal content of the air) and latent heat that is consumed for transpiration. The energy balance considerations show how the difference between the canopy and air temperatures (*T<sup>c</sup>* − *Ta*) is related to the vapor pressure deficit (VPD), and the flux density of the net radiation (*Rnet*) presents a theoretical basis for the CWSI.

When the plant canopy is fully transpiring, the leaf temperature is some degrees below the overlying air layer temperature and the CWSI is equal to 0 (stomata are closed). Conversely, as the transpiration decreases, the leaf temperature rises and can reach some degrees above the overlying air temperature. When the canopy is no longer transpiring, the CWSI is equal to 1. Because of the scatter in the measured canopy minus air temperature vs. the vapor pressure deficit, the crop does not need to be watered until the CWSI reaches 0.1 to 0.15. For example, for eggplant crops, it could be possible to use values of CWSI between 0.18–0.20 for high and good quality yields. Irmak et al.'s [44] work showed that when the CWSI value exceeds more than 0.22, this resulted in a decreased corn grain yield. Determining the temperature upper line (maximum stress) is not a simple process for any crop. However, this empirical approach has attracted much interest due to its simplicity. It requires measuring the canopy temperature (*Tc*), air temperature (*Ta*), and the vapor pressure of the atmosphere's deficit (VPD).

After the empirical approach suggested by Idso et al. [35], a theoretical method for calculating the CWSI was presented. The index's practical estimation is calculated by determining the relative distance between the lower baseline representing the soil water adequacy conditions in the rhizosphere (no stress) and the upper baseline representing the crop's maximum water stress (no transpiration by plants; Figure 1).

**Figure 1.** A schematic representation of the CWSI, which equals the relative distance of the measurement between the lower and the upper baselines. **Figure 1.** A schematic representation of the CWSI, which equals the relative distance of the measurement between the lower and the upper baselines.

The two lines (blue and red) presented in Figure 1 are suitable for assessing a crop's water status during the warm seasons of the year when irrigation is required. These lines are the result of local calibration and concern the specific crop that is studied each time. The lower baseline is determined for the crop by plotting (Tc − Ta) against the local VPD after a significant watering. A sufficient number of measurements (10 min time step) around midday were taken for the precise estimation of the canopy temperature (Tc). The canopy's temperature sampling must be repeated several times, including measurements under dry and wet atmospheric conditions. This results in a full range of measures representing the lower regression line for a wide range of value pairs (Τc − Ta vs. VPD). The upper baseline could be determined after cutting the plant and taking the canopy temperature the next day, as the plant no longer transpires. The theoretical procedure requires a measurement of the net flux density of the radiation (Rn) and an aerodynamic resistance factor, in addition to the air temperature and relative humidity of the air that is required by the empirical approach [36]. The two lines (blue and red) presented in Figure 1 are suitable for assessing a crop's water status during the warm seasons of the year when irrigation is required. These lines are the result of local calibration and concern the specific crop that is studied each time. The lower baseline is determined for the crop by plotting (*T<sup>c</sup>* − *Ta*) against the local VPD after a significant watering. A sufficient number of measurements (10 min time step) around midday were taken for the precise estimation of the canopy temperature (*Tc*). The canopy's temperature sampling must be repeated several times, including measurements under dryand wet atmospheric conditions. This results in a full range of measures representing the lower regression line for a wide range of value pairs (Tc − *T<sup>a</sup>* vs. VPD). The upper baseline could be determined after cutting the plant and taking the canopy temperature the next day, as the plant no longer transpires. The theoretical procedure requires a measurement of the net flux density of the radiation (*Rn*) and an aerodynamic resistance factor, in addition to the air temperature and relative humidity of the air that is required by the empirical approach [36].

The theoretical development of the CWSI is based on the well-known equation of energy balance of the surface ( *Rn* −= + *GH E*λ ), where *Rn* : net radiation, *G*: soil heat flux, *H*: sensible heat flux, and λ*E* : latent heat flux of evaporation. All the terms are in Wmିଶ. The terms *H* and λ*E* are expressed in Equations (1) and (2), which are based on two assumptions (see [36,45–47]): The theoretical development of the CWSI is based on the well-known equation of energy balance of the surface (*R<sup>n</sup>* − *G* = *H* + *λE*), where *Rn*: net radiation, *G*: soil heat flux, *H*: sensible heat flux, and *λE*: latent heat flux of evaporation. All the terms are in Wm−<sup>2</sup> . The terms *H* and *λE* are expressed in Equations (1) and (2), which are based on two assumptions (see [36,45–47]):

$$H = \rho \mathbb{C}\_p (T\_\mathfrak{c} - T\_\mathfrak{a}) / r\_\mathfrak{a} \tag{1}$$

$$
\lambda E = \left[\rho \mathbf{C}\_p (e\_\mathbf{s} - e\_\mathbf{d})\right] / \left[\gamma (r\_\mathbf{a} + r\_\mathbf{c})\right] \tag{2}
$$

where : density of air (kg m−3), *Cp*: specific heat moist air (1.013 kJ kg−1 °C−1), : canopy temperature (°C), : air temperature (°C), *<sup>s</sup> e* : saturated vapor pressure (kPa), *ae* : actual vapor pressure (kPa), : psychrometric constant (kPa °C−1), *ar* : aerodynamic resistance where *ρ* : density of air (kg m−<sup>3</sup> ), *Cp*: specific heat moist air (1.013 kJ kg−<sup>1</sup> ◦C −1 ), *T<sup>c</sup>* : canopy temperature (◦C), *Ta*: air temperature (◦C), *e<sup>s</sup>* : saturated vapor pressure (kPa), *ea*: actual vapor pressure (kPa), *γ* : psychrometric constant (kPa ◦C −1 ), *ra*: aerodynamic resistance (sm−<sup>1</sup> ), and *rc*: bulk surface resistance (s m−<sup>1</sup> ).

(sm−1), and : bulk surface resistance (s m−1). Clothier et al. [48] showed that the soil heat flux estimation can be calculated from *Rn* measurements with reasonable accuracy. The ratio *G/Rn* for full canopies (more than Clothier et al. [48] showed that the soil heat flux estimation can be calculated from *R<sup>n</sup>* measurements with reasonable accuracy. The ratio *G*/*R<sup>n</sup>* for full canopies (more than 45 cm) is nearly constant at 0.1. Thus, the soil heat flux does not contribute significantly to the

45 cm) is nearly constant at 0.1. Thus, the soil heat flux does not contribute significantly to the energy balance above the canopy. Taking 0.1 *G R* = *net* the equation becomes

*IR H E c n* = +λ , where *cI* = 0.9 is an interception coefficient.

energy balance above the canopy. Taking *G* = 0.1*Rnet* the equation becomes *IcR<sup>n</sup>* = *H* + *λE*, where *I<sup>c</sup>* = 0.9 is an interception coefficient.

$$T\_c - T\_a = \frac{r\_a I\_c R\_n}{\rho \mathcal{C}\_p} \cdot \frac{\gamma^\*}{\Delta + \gamma^\*} - \frac{VPD}{\Delta + \gamma^\*} \tag{3}$$

where *γ\** (kPa ◦C −1 ) is the modified psychrometric constant and VPD is the vapor pressure deficit.

Equation (4) represents the slope of the saturated vapor pressure relation:

$$
\Delta = (e\_{\rm s\_C} - e\_{\rm s}) / (T\_{\rm c} - T\_{\rm a}) \tag{4}
$$

where *es<sup>C</sup>* and *e<sup>s</sup>* represent the saturated vapor pressure of the canopy temperature and the air temperature, respectively.

The upper limit *T<sup>c</sup>* − *T<sup>a</sup>* is determined when the canopy resistance tends to infinity (*T<sup>c</sup>* − *T<sup>a</sup>* → ∞) and Equation (3) becomes

$$[T\_c - T\_a]\_{ul} = \frac{r\_a I\_{cu} R\_n}{\rho \mathcal{C}\_p} \tag{5}$$

The lower limit, determined by setting the bulk surface resistance *r<sup>c</sup>* = 0, is

$$T\_c - T\_a = \frac{r\_a I\_{cl} R\_n}{\rho \mathcal{C}\_p} \cdot \frac{\gamma}{\Delta + \gamma} - \frac{VPD}{\Delta + \gamma} \tag{6}$$

Thom and Oliver [49] proposed an effective aerodynamic resistance *rae* that represents a semi-empirical equation for *ra*:

$$r\_{a\varepsilon} = 4.72 \text{l} \{ \ln[(z - d) / z\_o] \}^2 / (1 + 0.54u) \tag{7}$$

The crop water stress index can be estimated using the equation below:

$$\text{CWSI} = \frac{(T\_c - T\_a) - (T\_c - T\_a)\_{ll}}{(T\_c - T\_a)\_{ul} - (T\_c - T\_a)\_{ll}} \tag{8}$$

#### *2.2. Structure and Architecture of the GWD Project*

The GreenWaterDrone system is functionally and physically divided into four subsystems, as follows (Figure 2):


multimedia, end-user preferences, and interfaces with external services (satellite imagery, photogrammetry applications). In addition, it interconnects and supports all other subsystems and is responsible for providing the services of the system (alerting and multimedia content) to all types of supported end users. other subsystems and is responsible for providing the services of the system (alerting and multimedia content) to all types of supported end users.

• The service provision I/Fs (FrontEnd) includes appropriate web interfaces of the system to predefined types of GWD users, such as plain (farmers/agronomists), group (partnerships), and strategic (local/regional authorities) end users, with graded access to the three supported applications through different devices (PCs, smartphones, etc.) and relevant GUIs. • The service provision I/Fs (FrontEnd) includes appropriate web interfaces of the system to predefined types of GWD users, such as plain (farmers/agronomists), group (partnerships), and strategic (local/regional authorities) end users, with graded access to the three supported applications through different devices (PCs, smartphones, etc.) and relevant GUIs.

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 6 of 28

**Figure 2.** The concise physical architecture of the system with the basic physical entities, implementing the respective groups of functional entities and their placement in the above subsystems, as well as the relevant interfaces between entities and subsystems. **Figure 2.** The concise physical architecture of the system with the basic physical entities, implementing the respective groups of functional entities and their placement in the above subsystems, as well as the relevant interfaces between entities and subsystems.

> The GWD system supports the following applications to provide a comprehensive irrigation management and field surveillance framework to both the plain and strategic end users that are targeted: The GWD system supports the following applications to provide a comprehensive irrigation management and field surveillance framework to both the plain and strategic end users that are targeted:


**3. Case Study—Materials and Methodology** 

*3.1. Study Area* 

## **3. Case Study—Materials and Methodology**

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 7 of 28

#### *3.1. Study Area*

The study area was situated in the prefecture of Messinia (western Peloponnese), close to the city of Kyparissia in the Municipality of Trifilia (Figure 3). Based on long-term data obtained from the nearby station of Kalamata (37.07◦ N, 22.10◦ E, alt. 8 m a.s.l.) over the last century, the broader area has a humid climate according to Thornthwaite's climate classification system (UNEP 1992), which has, however, become drier nowadays compared to the past, as seen via the decreasing aridity index [50] values from 0.88–0.89 in the previous climatic periods (1900–1930 and 1930–1960) to 0.85 in the period 1960–1997 [51]. The study area was situated in the prefecture of Messinia (western Peloponnese), close to the city of Kyparissia in the Municipality of Trifilia (Figure 3). Based on long-term data obtained from the nearby station of Kalamata (37.07° N, 22.10° E, alt. 8 m a.s.l.) over the last century, the broader area has a humid climate according to Thornthwaite's climate classification system (UNEP 1992), which has, however, become drier nowadays compared to the past, as seen via the decreasing aridity index [50] values from 0.88–0.89 in the previous climatic periods (1900–1930 and 1930–1960) to 0.85 in the period 1960–1997 [51].

**Figure 3.** Location of the GWD project's experimental fields in Trifilia. **Figure 3.** Location of the GWD project's experimental fields in Trifilia.

The experimentation was conducted in three fields with three different crops that are representative of the area, i.e., potato, watermelon, and tomato, which were utilized for the validation procedure of the GWD project methodology (Figure 3). Indicatively, in this work, potato field 2 (*Solanum tuberosum*; 1.25 hectares; 37°13′20.65″ N, 21°36′41.51″ E, alt. 1 m a.s.l.) was used throughout the next sections for demonstrating the preliminary results of the different imaging used in the project and GWD methodology validation during the 2019 growing season. The potato crop was selected by considering its relatively high water requirements, with its maximum in the hot and dry summer period, and the critical susceptibility of potato plants to water stress, which is attributed mainly to their relatively shallow root system (Figure 4a) [18]. The field had a small slope of 0–3% and a clay loam surface soil layer (0–20 cm depth), becoming sandy clay loam further down (20–60 cm depth; Figure 4b). The experimentation was conducted in three fields with three different crops that are representative of the area, i.e., potato, watermelon, and tomato, which were utilized for the validation procedure of the GWD project methodology (Figure 3). Indicatively, in this work, potato field 2 (*Solanum tuberosum*; 1.25 hectares; 37◦13020.6500 N, 21◦36041.5100 E, alt. 1 m a.s.l.) was used throughout the next sections for demonstrating the preliminary results of the different imaging used in the project and GWD methodology validation during the 2019 growing season. The potato crop was selected by considering its relatively high water requirements, with its maximum in the hot and dry summer period, and the critical susceptibility of potato plants to water stress, which is attributed mainly to their relatively shallow root system (Figure 4a) [18]. The field had a small slope of 0–3% and a clay loam surface soil layer (0–20 cm depth), becoming sandy clay loam further down (20–60 cm depth; Figure 4b).

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 8 of 28

**Figure 4.** (**a**) Experimental field 2 with the potato cultivation; (**b**) categorization of the soils of the three experimental fields at depths of 0–20, 20–40, and 40–60 cm in the soil triangle. The soil texture classification according to the USDA (U.S. Department of Agriculture). **Figure 4.** (**a**) Experimental field 2 with the potato cultivation; (**b**) categorization of the soils of the three experimental fields at depths of 0–20, 20–40, and 40–60 cm in the soil triangle. The soil texture classification according to the USDA (U.S. Department of Agriculture). **Figure 4.** (**a**) Experimental field 2 with the potato cultivation; (**b**) categorization of the soils of the three experimental fields at depths of 0–20, 20–40, and 40–60 cm in the soil triangle. The soil texture classification according to the USDA (U.S. Department of Agriculture).

#### *3.2. Ground Micrometeorological Station 3.2. Ground Micrometeorological Station 3.2. Ground Micrometeorological Station*

Three automatic micrometeorological stations were installed in the middle of the fields (Figure 5) for observations of the upper atmospheric layer of the crops and measurements within the soil profile (at the root zone). This was an indispensable condition for the accurate estimations of CWSI and the determination of all the biophysical attributes of the cultivations during all stages of the growing season. This ground-based measuring system composed the core of the calibration for estimating the CWSI of each crop, which was the basis for determining the real-time irrigation needs according to the methodology proposed by the GWD project. Three automatic micrometeorological stations were installed in the middle of the fields (Figure 5) for observations of the upper atmospheric layer of the crops and measurements within the soil profile (at the root zone). This was an indispensable condition for the accurate estimations of CWSI and the determination of all the biophysical attributes of the cultivations during all stages of the growing season. This ground-based measuring system composed the core of the calibration for estimating the CWSI of each crop, which was the basis for determining the real-time irrigation needs according to the methodology proposed by the GWD project. Three automatic micrometeorological stations were installed in the middle of the fields (Figure 5) for observations of the upper atmospheric layer of the crops and measurements within the soil profile (at the root zone). This was an indispensable condition for the accurate estimations of CWSI and the determination of all the biophysical attributes of the cultivations during all stages of the growing season. This ground-based measuring system composed the core of the calibration for estimating the CWSI of each crop, which was the basis for determining the real-time irrigation needs according to the methodology proposed by the GWD project.

in field 1, with early watermelon cultivation with a plastic cover for protection against the spring frost; (**b**) the second station in field 2 with the cultivation of spring potato in the middle growing stage; (**c**) the third station installed in field 3 with young tomato plants. **Figure 5.** Three stations were installed in the middle of the experimental fields: (**a**) the first station in field 1, with early watermelon cultivation with a plastic cover for protection against the spring frost; (**b**) the second station in field 2 with the cultivation of spring potato in the middle growing stage; (**c**) the third station installed in field 3 with young tomato plants. **Figure 5.** Three stations were installed in the middle of the experimental fields: (**a**) the first station in field 1, with early watermelon cultivation with a plastic cover for protection against the spring frost; (**b**) the second station in field 2 with the cultivation of spring potato in the middle growing stage; (**c**) the third station installed in field 3 with young tomato plants.

**Figure 5.** Three stations were installed in the middle of the experimental fields: (**a**) the first station

#### 3.2.1. Measurements above Crops

The three ground micrometeorological stations (GMMS) were equipped with the following sensors:


Given the importance of measuring the temperature of the plant foliage (canopy) remotely, reliable infrared radiometer sensors were crucial for achieving the program's objectives. The model SI-111-SS (Apogee, Berkeley, CA, USA), which is an unamplified analog sensor with a standard field of view, was selected to measure the canopy temperature. The response time of the sensor is 0.6 s, with a measurement repeatability of less than 0.05 ◦C, and the calibration uncertainty is 0.2 ◦C. The spectral range of the measurement is between 8 and 14 µm, that is, the atmospheric window. The sensor has a 22◦ half-angle field of view. The placement height and the angle formed between the mast with the longitudinal axis of the infrared thermometer determine the elliptical area of measurement (in m<sup>2</sup> ) of the surface of the dense canopy.

Moreover, portable infrared thermometers (Model MI-210 Apogee, USA) were used for conducting spatial measurements in all fields' canopies during the flights of the UAV for calibration and comparative data analysis.

#### 3.2.2. Measurements in the Root Zone

The soil moisture (SMois) (*v*/*v*%) and soil temperature (Tsoil) were measured continuously at six soil layers (from 0–10 cm to 50–60 cm in 10 cm intervals) using the Drill and Drop probe (model 00620 Sentek Technologies, Australia). In this way, the actual changes in soil moisture were constantly obtained throughout the depth of the root zone at any time. Furthermore, a soil heat flux plate (HFP01 Hukseflux Thermal Sensor, Netherlands) was installed in the topsoil layer (at a 5 cm depth). Figure 6 graphically shows the soil moisture profile estimated using a geostatistical gridding method (Kriging), which was proven to be useful in many fields. The Kriging method estimates the surface at successive nodes in the grid using only a selection of the nearly closed data points and produces surfaces from irregularly spaced data.

surfaces from irregularly spaced data.

nodes in the grid using only a selection of the nearly closed data points and produces

The method efficiently and naturally incorporates anisotropy and underlying trends from a data set by specifying the appropriate variogram model. To create this graph, 1080 daily soil moisture measurements were used from six depths that were monitored in the potato cultivation field for 51 days. The method efficiently and naturally incorporates anisotropy and underlying trends from a data set by specifying the appropriate variogram model. To create this graph, 1080 daily soil moisture measurements were used from six depths that were monitored in the potato cultivation field for 51 days.

For the detailed fluctuations of the soil moisture of the potato field at depths of 0–10 and 10–20 cm (the two most variable layers), we used 17,284 average values of ten-minute continuous soil moisture measurements. These measurements are shown in Figure 6 and refer to the interval from 10 May to 9 July. The smooth, continuous, and detailed datacollection process, which took place throughout the growing season in the soil profiles, allowed us to accurately estimate the actual evapotranspiration (ETa) from the soil moisture profile in hourly time steps. Moreover, from the observations and characteristics of the plants in all stages (height, LAI, aerodynamic, and surface resistance), in combination with the micrometeorological observations just above the crop (Rs, Rn, albedo, soil heat flux, T, RH), we estimated the crop potential evapotranspiration (ETc) with satisfactory accuracy. All the above estimates contributed to estimating the crop water stress index (CWSI), both practically and theoretically. For the detailed fluctuations of the soil moisture of the potato field at depths of 0–10 and 10–20 cm (the two most variable layers), we used 17,284 average values of ten-minute continuous soil moisture measurements. These measurements are shown in Figure 6 and refer to the interval from 10 May to 9 July. The smooth, continuous, and detailed datacollection process, which took place throughout the growing season in the soil profiles, allowed us to accurately estimate the actual evapotranspiration (ETa) from the soil moisture profile in hourly time steps. Moreover, from the observations and characteristics of the plants in all stages (height, LAI, aerodynamic, and surface resistance), in combination with the micrometeorological observations just above the crop (Rs, Rn, albedo, soil heat flux, T, RH), we estimated the crop potential evapotranspiration (E*Tc*) with satisfactory accuracy. All the above estimates contributed to estimating the crop water stress index (CWSI), both practically and theoretically.

#### 3.2.3. Data Acquisition System and Telemetry 3.2.3. Data Acquisition System and Telemetry

The telemetry unit (model A753 addWAVE GPRS RTU, ADCON (Business Unit of OTT HydroMet GmbH)) provided the possibility of local storage in EPROM Memory for at least half a million measurements. The telemetry unit has an internal rechargeable battery that is charged using a solar panel, which is independent of the one used for the energy needs of the sensors. The duration of operation in normal operation is up to 14 days. According to the GSM standard, the frequency band is 850/900/1800/1900 MHz, and the maximum transmission distance is 36 km. To achieve the maximum accuracy (connecting sensors to a central station unit), especially in the connection of thermocouple sensors (pyranometers, soil heat flux plate) to a central station a unit HD978TR3 (Delta OHM), The telemetry unit (model A753 addWAVE GPRS RTU, ADCON (Business Unit of OTT HydroMet GmbH)) provided the possibility of local storage in EPROM Memory for at least half a million measurements. The telemetry unit has an internal rechargeable battery that is charged using a solar panel, which is independent of the one used for the energy needs of the sensors. The duration of operation in normal operation is up to 14 days. According to the GSM standard, the frequency band is 850/900/1800/1900 MHz, and the maximum transmission distance is 36 km. To achieve the maximum accuracy (connecting sensors to a central station unit), especially in the connection of thermocouple sensors (pyranometers, soil heat flux plate) to a central station a unit HD978TR3 (Delta OHM), signal amplifiers were used.

signal amplifiers were used. The time step for storing the average data values for each sensor is ten minutes for all stations (144 average values per 24 h). From the ten-minute average values, we obtained and recorded hourly and daily average values for all sensors. All average values from each sensor were found using 4320 reads per 24 h. This means that almost all sensors were automatically excited to give a reading sample every 20 s (3 readings per minute). Of course, this programming considers the response time of each sensor, i.e., the ability of the sensor to respond to two different successive values in time.

#### *3.3. Experimental Design—Unmanned Aerial Vehicle, Cameras, and IR-TH*

In the present study, the quadcopter DJI Matrice 200 (DJI Technology Co., Ltd., Shenzhen, China) was used for the measurements of the multispectral and IR-TH data, and the fixed-wing Q100 Datahawk was employed for the photogrammetric mapping of the fields. The quadcopter was selected due to the compatibility of incorporating both multispectral and thermal cameras, and because the flights with infrared sensors require a low flight altitude, low-speed navigation, and high stability, which cannot be achieved with a fixed-wing platform [52]. The accuracy of the system is 10 cm in the vertical and 30 cm in the horizontal, which is considered suitable for PA applications, while at the same time, it offers the possibilities of low, slow, and stable flight. It has significant autonomy and a maximum take-off mass with a large payload (6.14 and 1.61 kg, respectively) and can withstand strong winds (up to 12 m/s), which ensures its smooth operation in demanding flight environments. In addition, the system incorporates an advanced obstacle recognition and avoidance subsystem, which consists of a combination of optical, infrared, and ultrasonic sensors, making the platform one of the safest for flight operations in existence, and incorporates inertial measurement units (IMU) sensors, supporting both automatic, semi-automatic, and manual flights.

The flight planning of the photogrammetric, multispectral, and thermal cameras was conducted with the ground control station software, GS PRO (for the IR), Pix4D (for the thermal camera), Field Agent (for the multispectral camera), and the SkyCircuits Plan and Flight (for the photogrammetric camera), which allow the user to generate a route over a grid path as a function of the field of view of the sensor, flight altitude, flight speed, direction, degree of overlap between images, and ground resolution.

#### 3.3.1. UAV for IR Spatial Row Data Measurements

The DJI Matrice quadcopter was selected to install an autonomous aerial micrometeorological system (MicroStation) consisting of a high-performance data logger (Symetron type Stylitis-12), as well as four sensors: IR radiometer (SI-111-SS Apogee), miniature Thermo-Hygrometer EE08, and an accurate GPS (Figure 7). Through the "Opton 4" Symetron software, the data logger operation was managed and controlled. Furthermore, all the appropriate equations and precalibrated factors for any crop were incorporated into the CWSI estimation software.

The radiation shield plates (which protected the temperature and RH sensors), the watertight box (which protected the data logger), and all the sensors' mounts were made from lightweight synthetic materials. All these components were designed in a 3D CAD program and printed using a 3D printer in the lab. Particular attention was paid to the placement of the thermo-hygrometer so that it was not affected by the flow of the drone propellers (Figure 7). The vertical positioning of the radiation shield on top of the UAV ensured mild air vortexing and mixing, avoiding the violent airflow from the downside of the drone due to the propellers.

One of the most important parts of the IR instrument sampling workflow was the calculation and design of the most suitable flight planning in order to have the best point samples coverage and complementarity to ideally cover the whole cultivated area. The flight height depends on the sensor type and the characteristics of the crop. For a specific IRT sensor, the visible surface area is a function of height. For the sensor SI-111 that was placed on the UAV in such a way to be vertical to the foliage for a flight height of 4 m (h), the surface area was estimated to be 8.2 m<sup>2</sup> (E), and the diameter of the circular surface was about 3.2 m (Figure 8).

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 12 of 28

**Figure 7.** (**a**) Watertight box (data logger Stylitis-12 inside); (**b**) temperature and relative humidity (thermo-hygrometer) sensor in the radiation shield; (**c**) infrared sensor and the mounting base; (**d**) GPS sensor. **Figure 7.** (**a**) Watertight box (data logger Stylitis-12 inside); (**b**) temperature and relative humidity (thermo-hygrometer) sensor in the radiation shield; (**c**) infrared sensor and the mounting base; (**d**) GPS sensor. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 13 of 28

The radiation shield plates (which protected the temperature and RH sensors), the

**Figure 8.** Graphic illustration of the measurement process with an infrared radiometer adapted to a UAV over dense cultivation foliage under a clear sky. The flight height depends on the sensor type and the characteristics of the crop. For a specific IRT sensor, the visible surface area is a function of height. For the sensor SI-111 that was placed vertically to the foliage with a height of flight h = 4 m, the surface area was E = 8.2 m2, and the diameter of the circular surface was about 3.2 m. **Figure 8.** Graphic illustration of the measurement process with an infrared radiometer adapted to a UAV over dense cultivation foliage under a clear sky. The flight height depends on the sensor type and the characteristics of the crop. For a specific IRT sensor, the visible surface area is a function of height. For the sensor SI-111 that was placed vertically to the foliage with a height of flight h = 4 m, the surface area was E = 8.2 m<sup>2</sup> , and the diameter of the circular surface was about 3.2 m.

The photogrammetric camera that was used for the present study was the SONY rx 100, with a 20.2 megapixels CMOS lens that has a ground resolution of 3.2 cm/pixel. The camera was adjusted on a fixed-wing Q100 Datahawk UAV (Figure 9a). By processing the

plore the surface slope inclination and, indirectly, the conditions of drainage and water flow in the soil, as well as the growth conditions of the crop (canopy height and density).

3.3.2. UAVs for Photogrammetric, Multispectral, and Thermal Measurements

#### 3.3.2. UAVs for Photogrammetric, Multispectral, and Thermal Measurements

The photogrammetric camera that was used for the present study was the SONY rx 100, with a 20.2 megapixels CMOS lens that has a ground resolution of 3.2 cm/pixel. The camera was adjusted on a fixed-wing Q100 Datahawk UAV (Figure 9a). By processing the images that were produced, an accurate description of the canopy height and a digital surface model (DSM) were created. This information is significant in endeavoring to explore the surface slope inclination and, indirectly, the conditions of drainage and water flow in the soil, as well as the growth conditions of the crop (canopy height and density). *Hydrology* **2021**, *8*, x FOR PEER REVIEW 14 of 28

**Figure 9.** (**a**) The fixed-wing Q100 Datahawk; (**b**) the multispectral Sentera camera; (**c**) the thermal infrared camera Zenmuse XT2. **Figure 9.** (**a**) The fixed-wing Q100 Datahawk; (**b**) the multispectral Sentera camera; (**c**) the thermal infrared camera Zenmuse XT2.

The DSM from UAV image stereo pairs used the structure from motion (SfM) algorithm, which is a photogrammetric technique that was initially developed for archaeological sites [53–55]. The algorithm assesses the unknown camera orientations through a comparison of multiple detected image feature points in multiple images. It subsequently produces the 3D structural model of a scene from the overlapping two-dimensional (2D) image sequences that are taken from various spots and orientations [56]. The DSM from UAV image stereo pairs used the structure from motion (SfM) algorithm, which is a photogrammetric technique that was initially developed for archaeological sites [53–55]. The algorithm assesses the unknown camera orientations through a comparison of multiple detected image feature points in multiple images. It subsequently produces the 3D structural model of a scene from the overlapping two-dimensional (2D) image sequences that are taken from various spots and orientations [56].

The multispectral camera that was used in the present study was the Sentera AGX710 (12.3 megapixels) multi-spectral sensor with five bands (Sentera Inc., Minneapolis, MN, USA), which cover the blue, green, red, red edge, and near-infrared parts of the electromagnetic spectrum (Figure 9b). Their central wavelengths are 446 nm, 548 nm, 650 nm, 720 nm, and 840 nm, respectively. A 15 cm × 15 cm white reference panel (MicaSense Inc., Seattle, DC, USA) with a 60% nominal reflectance was used for the radiometric correction. UAV imaging was conducted in sunny weather and the period between 11 a.m. and 12 p.m. was chosen for imaging to minimize sunshade. To certify correct image acquisition during the flight, the FieldAgent Mobile App from Sentera was used to plan the flight path and automate the operation of the UAV. A constant flight height was maintained at 30 m and the flight speed was set to 3 m/s. The ground sample distance of the imagery was approximately 1 cm/pixel and a 90% overlap between two images both for side-lap and front-lap was implemented. In total, about 820 images were taken for the experimental field and stacked into a single image for analysis. The multispectral camera that was used in the present study was the Sentera AGX710 (12.3 megapixels) multi-spectral sensor with five bands (Sentera Inc., Minneapolis, MN, USA), which cover the blue, green, red, red edge, and near-infrared parts of the electromagnetic spectrum (Figure 9b). Their central wavelengths are 446 nm, 548 nm, 650 nm, 720 nm, and 840 nm, respectively. A 15 cm × 15 cm white reference panel (MicaSense Inc., Seattle, DC, USA) with a 60% nominal reflectance was used for the radiometric correction. UAV imaging was conducted in sunny weather and the period between 11 a.m. and 12 p.m. was chosen for imaging to minimize sunshade. To certify correct image acquisition during the flight, the FieldAgent Mobile App from Sentera was used to plan the flight path and automate the operation of the UAV. A constant flight height was maintained at 30 m and the flight speed was set to 3 m/s. The ground sample distance of the imagery was approximately 1 cm/pixel and a 90% overlap between two images both for side-lap and front-lap was implemented. In total, about 820 images were taken for the experimental field and stacked into a single image for analysis.

Agisoft Metashape Professional (v 1.5.5) was preferred as the Structure for Motion Multi-View Stereo (SfM-MVS) processing software to generate the digital surface model (DSM) and the correspondent normalized difference red-edge index (NDRE) vegetation Agisoft Metashape Professional (v 1.5.5) was preferred as the Structure for Motion Multi-View Stereo (SfM-MVS) processing software to generate the digital surface model (DSM) and the correspondent normalized difference red-edge index (NDRE) vegetation index.

overall amount and quality of photosynthetic material in vegetation, which is essential for understanding the state of vegetation (photosynthetic capacity and canopy chlorophyll content, structure, and plant health status). NDRE (Equation (9)) was measured by using the NIR and red-edge band (RE; 705 nm), where RE replaces the red band in the widely used normalized difference vegetation index (NDVI) equation. Compared with the NDVI, the NDRE index as a widely used red-edge-based VI was shown to be more resistant to

index.

NDRE is a narrowband greenness VI that was designed to provide a measure of the overall amount and quality of photosynthetic material in vegetation, which is essential for understanding the state of vegetation (photosynthetic capacity and canopy chlorophyll content, structure, and plant health status). NDRE (Equation (9)) was measured by using the NIR and red-edge band (RE; 705 nm), where RE replaces the red band in the widely used normalized difference vegetation index (NDVI) equation. Compared with the NDVI, the NDRE index as a widely used red-edge-based VI was shown to be more resistant to the saturation problem and is more sensitive than the NDVI to chlorophyll content in vegetation [57].

$$\text{NDRE} = \frac{\left(\text{B}\_{\text{NIR}} - \text{B}\_{\text{RedEdge}}\right)}{\left(\text{B}\_{\text{NIR}} + \text{B}\_{\text{RedEdge}}\right)}\tag{9}$$

The final image nominal spatial resolution, expressed as a ground sample distance (GSD), was about 1.2 cm [58].

The final NDRE image ideally demonstrates the plant health status and the existing gaps in the canopy cover and along the cultivation lines, which are important to distinguish and exclude since the reflection from the soil significantly affects the calculation results of the CWSI. Thus, the results of the NDRE image were used to create a binary mask of the soil. Sunlit soil parts were extracted from the image by utilizing a threshold value on a local spatial statistic that characterized the illumination levels [18,59]. This final mask was applied on the co-registered TIR image to produce the thermal product with pixels belonging only to sunlit leaves.

The thermal camera used in the present study was the DJI Zenmuse XT2 (12 MP), which incorporates a high-resolution forward-looking infrared FLIR Tau 2 thermal sensor (FLIR Systems, Wilsonville, OR, USA; resolution 640 × 512 pixels, lens 9 mm; thermal sensitivity of <50 mK) and a 4 K visual camera (1/1.7" active-pixel CMOS) with a leading stabilization and machine intelligence technology (Figure 9c) [60]. In general, the Zenmuse XT2 camera used in the present study is a very high technological and sophisticated camera that performs radiometric calibration of the acquired thermal images through its advanced digital system, which records the absolute temperature. The calibration of the FLIR camera was performed annually according to the factory specifications, which is carried out by the official supplier by measuring targets with known temperatures and comparing the known and the measured temperatures (e.g., boiling water and melting ice) [16,61]. The thermal lens has a focal length of 19 mm, which avoids previous deformations and obtains rectilinear images. The second lens has a focal length of 8 mm and it also retrieves highresolution RGB images [62]. The ground sample distance of the imagery was approximately 2 cm/pixel and an 85% overlap between two images for both side-lap and front-lap was implemented. In total, about 382 images were acquired for the experimental field and stacked into a single image for analysis.

UAV thermal imaging was performed in sunny weather from 12:00 to 1:00 p.m. to minimize sunshade and because it has general unique image-acquiring conditions since the maximum temperature differences occur between the soil and vegetation [63,64].

The thermal images acquired by UAV were mosaicked using Pix4D mapper (version 4.5 from PIX4D, Prilly, Switzerland). This software offers an automated SfM procedure, taking a set of images as input and routinely going through the steps of feature identification, matching, and bundle adjustment. The process aligned the images that were acquired by the thermal camera. A polygon mesh was calculated from the dense 3D point, and the pixel values of each image were then projected onto the mesh to create an orthomosaic. When combined with the GPS locations, this procedure allowed for the establishment of a high-resolution orthophoto and a digital elevation model (DEM) of the crop field [33,65].

The final steps considered the exclusion of the soil cover and the gaps among the canopy by using the very high resolution multispectral and true color orthomosaic. Then, the interpolation of the selected points for the selected area comprised the most homogeneous and appropriate parts of the canopy for the following measurements.

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 16 of 28 and were surveyed using an RTK GNSS with a range of <1 cm in the horizontal plane and 1.7 cm in the vertical axis (Figure 11).

For the better calibration of the thermal image and IR data derived from the UAV, field samples and continuous measurements of the meteorological station using the same IR instrument were utilized and evaluated (Figure 10). and were surveyed using an RTK GNSS with a range of <1 cm in the horizontal plane and 1.7 cm in the vertical axis (Figure 11).

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 16 of 28

**Figure 10.** The field measurements that were taken using the IR sensors on (**a**) the meteorological station and (**b**) a portable device. **Figure 10.** The field measurements that were taken using the IR sensors on (**a**) the meteorological station and (**b**) a portable device.

The information on GPS locations for each image was obtained from the UAV during collection; however, to achieve higher accuracy, eight rectangular aluminum plates were used as ground control points (GCPs), which were equally spread around the potato field and were surveyed using an RTK GNSS with a range of <1 cm in the horizontal plane and 1.7 cm in the vertical axis (Figure 11). **Figure 10.** The field measurements that were taken using the IR sensors on (**a**) the meteorological station and (**b**) a portable device.

Thus, a linear regression model using ground truth temperatures (measured with the handheld IR) and those obtained from the UAV was formed following analogous results that were in the recent literature, which found that the linear regression model based on **Figure 11.** (**a**) The RTK GNSS that was used for the collection of the GCPs; (**b**) the rectangular aluminum plates that were used as ground control points on the field **Figure 11.** (**a**) The RTK GNSS that was used for the collection of the GCPs; (**b**) the rectangular aluminum plates that were used as ground control points on the field.

data acquired at a 1 m height had a slope of 1.0, while that based on data acquired at a 50 m height had a slope of 1.4 [16]. **4. Preliminary Results and Discussion**  Finally, for the temperature calibration, a linear regression model for the precise acquisition height (30 m) of the UAV thermal data was applied since the temperature obtained using the IR sensors considerably decreased with increasing flight height [28,66,67]. Thus, a linear regression model using ground truth temperatures (measured with the handheld IR) and those obtained from the UAV was formed following analogous results Finally, for the temperature calibration, a linear regression model for the precise acquisition height (30 m) of the UAV thermal data was applied since the temperature obtained using the IR sensors considerably decreased with increasing flight height [28,66,67]. Thus, a linear regression model using ground truth temperatures (measured with the handheld

that were in the recent literature, which found that the linear regression model based on

The results listed below came from flights that took place over an autumn potato field (25 October 2019). The pilot flights were conducted around midday. The meteorological

The results listed below came from flights that took place over an autumn potato field (25 October 2019). The pilot flights were conducted around midday. The meteorological

m height had a slope of 1.4 [16].

**4. Preliminary Results and Discussion** 

IR) and those obtained from the UAV was formed following analogous results that were in the recent literature, which found that the linear regression model based on data acquired at a 1 m height had a slope of 1.0, while that based on data acquired at a 50 m height had a slope of 1.4 [16].

#### **4. Preliminary Results and Discussion** *Hydrology* **2021**, *8*, x FOR PEER REVIEW 17 of 28

The results listed below came from flights that took place over an autumn potato field (25 October 2019). The pilot flights were conducted around midday. The meteorological conditions prevailing during the day were characterized by a clear sky with maximum global solar and net radiation flux densities of 670 and 395 Wm−<sup>2</sup> , respectively, at 13:00 (Figure 12a). The air temperature attributes were normal for the season with a daily Tavg = 19.1 ◦C, Tmax = 24.8 ◦C, and Tmin = 14.4 ◦C (Figure 12b). The VPD, E*Tc*, and wind speed hourly averages depicted in Figure 12c indicate that during the flights, the wind speed was about 1.5 m/s and the E*T<sup>c</sup>* maximum rate was 0.43 mm/h. For the specific day, the total E*T<sup>c</sup>* reached 2.7 mm/day. The relatively low crop evapotranspiration rate was in line with the small changes in the soil moisture profile presented in Figure 12d. conditions prevailing during the day were characterized by a clear sky with maximum global solar and net radiation flux densities of 670 and 395 Wm−2, respectively, at 13:00 (Figure 12a). The air temperature attributes were normal for the season with a daily Tavg = 19.1 °C, Tmax = 24.8 °C, and Tmin = 14.4 °C (Figure 12b). The VPD, ETc, and wind speed hourly averages depicted in Figure 12c indicate that during the flights, the wind speed was about 1.5 m/s and the ETc maximum rate was 0.43 mm/h. For the specific day, the total ETc reached 2.7 mm/day. The relatively low crop evapotranspiration rate was in line with the small changes in the soil moisture profile presented in Figure 12d.

**Figure 12.** Hourly values for the date 25/10/2019 for the (**a**) global solar radiation Rs, net radiation Rnet, photosynthetically active radiation PAR, and soil heat flux SHF; (**b**) air and canopy temperatures; (**c**) vapor pressure deficit VPD, wind speed WS, crop evapotranspiration; (**d**) soil moisture profile for six depths measured in the potato field. **Figure 12.** Hourly values for the date 25/10/2019 for the (**a**) global solar radiation Rs, net radiation Rnet, photosynthetically active radiation PAR, and soil heat flux SHF; (**b**) air and canopy temperatures; (**c**) vapor pressure deficit VPD, wind speed WS, crop evapotranspiration; (**d**) soil moisture profile for six depths measured in the potato field.

*4.1. Photogrammetry—Multispectral and Thermal Imagery* 

revealed all the aspects of the canopy attributes and its surrounding area.

The photogrammetric analysis provided very high-resolution information about the canopy conditions by creating the RGB true-color imagery of the field (Figure 13a), which

#### *4.1. Photogrammetry—Multispectral and Thermal Imagery*

The photogrammetric analysis provided very high-resolution information about the canopy conditions by creating the RGB true-color imagery of the field (Figure 13a), which revealed all the aspects of the canopy attributes and its surrounding area. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 18 of 28

**Figure 13.** (**a**) The very high-resolution true color orthomosaic and (**b**) the very high-resolution digital surface model of field 2. **Figure 13.** (**a**) The very high-resolution true color orthomosaic and (**b**) the very high-resolution digital surface model of field 2.

The DSM indicated a gentle inclination of the field, which was ideal for the appropriate drainage of the field and avoiding overlogging situations, where the elevation varied from 1.1 m at the south-southeastern part to 0.4 m at the north-northwestern part of the field (Figure 13b). The DSM indicated a gentle inclination of the field, which was ideal for the appropriate drainage of the field and avoiding overlogging situations, where the elevation varied from 1.1 m at the south-southeastern part to 0.4 m at the north-northwestern part of the field (Figure 13b).

The result of the NDRE equation provided an image with continuum pixel values that ranged from −1 to 1, where negative values corresponded to non-vegetated surfaces, in this case, soil cover, while positive values corresponded to the vegetation reflectance (Figure 14). The values between 0.25 and 0.64 were related to healthy photosynthetic vegetation of potato canopy, while lower values, especially those <0.2, were related to stressed vegetation or bare soil. Values higher than 0.64 represented the dense tree cover in the northwestern and southeastern parts of the field. As it was mentioned before, the NDRE image was utilized to mask the soil cover from the co-registered TIR image (XT2) and the interpolation image derived from the IR sensor point values to produce the thermal products The result of the NDRE equation provided an image with continuum pixel values that ranged from −1 to 1, where negative values corresponded to non-vegetated surfaces, in this case, soil cover, while positive values corresponded to the vegetation reflectance (Figure 14). The values between 0.25 and 0.64 were related to healthy photosynthetic vegetation of potato canopy, while lower values, especially those <0.2, were related to stressed vegetation or bare soil. Values higher than 0.64 represented the dense tree cover in the northwestern and southeastern parts of the field. As it was mentioned before, the NDRE image was utilized to mask the soil cover from the co-registered TIR image (XT2) and the interpolation image derived from the IR sensor point values to produce the thermal products with pixels belonging only to sunlit leaves.

with pixels belonging only to sunlit leaves. The thermal image (Figure 15) that was taken by the Zenmuse XT2 on the UAV was used to calculate the canopy temperature (*Tc*) of the experimental field. The image was classified (using custom classification), keeping only the pixels of the crop canopy area to separate and present the reliable foliage temperature values. The canopy temperatures (from the thermal image) ranged between approximately 21.2 and 24.3 ◦C (range of about 3 ◦C), which was very close to the station average IRT temperature (23.9 ◦C) and the average spatial IRT (21.8 ◦C), respectively. These measurements were taken during a day with clear sky conditions and relatively low temperatures. The crop was in the middle stage of growth (almost full canopy cover), where the soil and atmospheric conditions are presented in the following charts.

**Figure 14.** The red-edge normalized difference vegetation index of field 2, which was calculated using the multispectral camera data. **Figure 14.** The red-edge normalized difference vegetation index of field 2, which was calculated using the multispectral camera data. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 20 of 28

**Figure 15.** Thermal image acquired by the Zenmuse XT2 camera showing only the crop vegetation's temperature values by masking all the non-vegetation pixels. **Figure 15.** Thermal image acquired by the Zenmuse XT2 camera showing only the crop vegetation's temperature values by masking all the non-vegetation pixels.

Although thermal RS has various potential advantages in crop monitoring, there are several practical difficulties in its use, including atmospheric attenuation and absorption, calibration, climatic conditions, crop growth stages, and complex soil and plant interaction that have to be adequately addressed. Although thermal RS has various potential advantages in crop monitoring, there are several practical difficulties in its use, including atmospheric attenuation and absorption, calibration, climatic conditions, crop growth stages, and complex soil and plant interaction that have to be adequately addressed.

The high-resolution multispectral and thermal imagery acquired using the UAV plat-

The first pilot flight with the UAV that brought the aMMS (aerial micrometeorological system) equipment took place on 25 October 2019 over field 2 during late potato cultivation. It was essentially the first comprehensive process for estimating the CWSI spatial mapping. In the same field, during the early cultivation of potatoes in May 2019, a previous estimation of CWSI calibration coefficients from specific days with ideal weather conditions was undertaken. We used all ground station data and parameters to determine the lower and upper baselines using theoretical and practical approaches. The day of the UAV flight was a clear sunny day with relatively low temperature levels, as the period was in the middle of autumn. The crop was in the middle stage of development. The UAV flight for the spatial measurements of the canopy temperature, air temperature, and relative humidity started at 12:14 p.m., which corresponded to the local solar time of the ground

opy. Furthermore, they demonstrate great capability for the recognition and monitoring of drought stress and were used in several different crops, such as potato, cotton, soybean, cotton, maize, vineyards, and orchards, over the last two decades [16,18,26,27,68–71]. Additionally, several studies showed that cloud cover and other haze or dew conditions can influence the quality of the thermal data obtained (due to lower contrast and poor radiometric resolution) [72–74]. Hence, these variables, along with relative humidity, altitude, and viewing angle, should be thoroughly examined while acquiring thermal data [71,73]. Moreover, the acquisition time during midday was demonstrated by several scientists that

were the optimal time for thermal image acquisition [75].

*4.2. Direct Canopy IRT Measurements and Spatial CWSI Estimation* 

The high-resolution multispectral and thermal imagery acquired using the UAV platform ideally enables the extraction and identification of crop conditions (vegetation vigor, spatial distribution, canopy cover, and bare ground gaps) and the temperature of the canopy. Furthermore, they demonstrate great capability for the recognition and monitoring of drought stress and were used in several different crops, such as potato, cotton, soybean, cotton, maize, vineyards, and orchards, over the last two decades [16,18,26,27,68–71]. Additionally, several studies showed that cloud cover and other haze or dew conditions can influence the quality of the thermal data obtained (due to lower contrast and poor radiometric resolution) [72–74]. Hence, these variables, along with relative humidity, altitude, and viewing angle, should be thoroughly examined while acquiring thermal data [71,73]. Moreover, the acquisition time during midday was demonstrated by several scientists that were the optimal time for thermal image acquisition [75].

#### *4.2. Direct Canopy IRT Measurements and Spatial CWSI Estimation*

The first pilot flight with the UAV that brought the aMMS (aerial micrometeorological system) equipment took place on 25 October 2019 over field 2 during late potato cultivation. It was essentially the first comprehensive process for estimating the CWSI spatial mapping. In the same field, during the early cultivation of potatoes in May 2019, a previous estimation of CWSI calibration coefficients from specific days with ideal weather conditions was undertaken. We used all ground station data and parameters to determine the lower and upper baselines using theoretical and practical approaches. The day of the UAV flight was a clear sunny day with relatively low temperature levels, as the period was in the middle of autumn. The crop was in the middle stage of development. The UAV flight for the spatial measurements of the canopy temperature, air temperature, and relative humidity started at 12:14 p.m., which corresponded to the local solar time of the ground station (lat. 37.222194 N, long. 21.611806 E in degrees). The flight duration was 21 min, and 217 primary spatial measurements of infrared temperature from the crop canopy, air temperature, and relative humidity were collected from a height of 4 m above the ground (about 3.5 m from the top of the foliage). The flight was autonomous and scheduled following the standard flight plan described in Section 3.3.1 (Figure 8). The speed of the quadcopter was 1 m/s, and the excitation of the sensors for measurement was every 6 s (measurement every 6 m along the flight line). All microclimatic data obtained from the ground station (GMMS) above the crop, as well as the mean spatial values (N = 217) from the sensors of the aerial micrometeorological system, are shown in Table 1.


**Table 1.** Ten-minute average climatic data from the ground micrometeorological station (GMMS) from field 2 and the averages values from 217 measurements of the air temperature, canopy infrared temperature, and relative humidity from the aerial micrometeorological system (AMMS).

\* Mean values of measurements.

12:35 24.1 24.2 62.2 320 1.2 665 393 GMMS \* 24.0 23.9 64.3 318 1.4 669 395

AMMS \* 24.3 21.8 56.0

Moreover, Table 2 presents the moisture conditions in the six layers (from 0–60 cm with 10 cm intervals), and the soil heat flux density at a depth of 8 cm. As shown in Figure 16, the surface layer (0–10 cm) had lost 67% of its moisture from the total height of the available water between two irrigation events. Usually, the irrigation frequency for this period was every 3–4 days, as shown in Figure 16 for the last three irrigation events (16, 20, and 23 October 19).


**Table 2.** Soil water moisture (vol%) for the depth of 60 cm of the soil profile in the potato field.

**Figure 16.** The fluctuation of the soil moisture in the surface layer of the soil (0–10 cm) after three **Figure 16.** The fluctuation of the soil moisture in the surface layer of the soil (0–10 cm) after three consecutive irrigations.

consecutive irrigations.

Figure 17 presents the isothermal lines from a spatial-observation-derived analysis (N = 217) using the geostatistical gridding Kriging method, which produces maps from Figure 17 presents the isothermal lines from a spatial-observation-derived analysis (N = 217) using the geostatistical gridding Kriging method, which produces maps from atypical, spaced data.

atypical, spaced data. The method is very flexible and produces an accurate data grid by a smoothing interpolator depending on the user-specified parameters. The spacing nodes grid size used 86 rows × 100 columns (total nodes: 8600). Some univariate grid statistics for the canopy IR temperature for the grid data (in ◦C) were min: 21.603, max: 22.133, mean: 21.886, root mean square: 21.8869, standard deviation: 0.1486, average abs. deviation: 0.1252, and relative mean diff.: 0.0078.

> The exact area (enclosed area by the boundary line) of field 2 was 9940.3 m<sup>2</sup> . Each observation was taken from about 4 m above the crop and corresponded to the mean value of the infrared temperature of the circular surface of the cone base radius (target) of about 3 m (Figure 8).

> The evaluation of the CWSI index during the clear sky and dry hours of the day offered a real and more direct method to determine the ideal time for irrigation, as it effectively found and normalized the short-term microclimatic changes of the environment in which the crop grew up and responded reliably to maintaining the optimal water content of the soil porosity into the rhizosphere zone. This significantly improved the irrigation efficiency and water-saving at all stages of plant growth. The spatial estimates of the CWSI are shown in Figure 18.

The method is very flexible and produces an accurate data grid by a smoothing interpolator depending on the user-specified parameters. The spacing nodes grid size used

**Figure 17.** The temperature values that were derived from the IR sensor measurements.

**Figure 16.** The fluctuation of the soil moisture in the surface layer of the soil (0–10 cm) after three

Figure 17 presents the isothermal lines from a spatial-observation-derived analysis (N = 217) using the geostatistical gridding Kriging method, which produces maps from

consecutive irrigations.

atypical, spaced data.

are shown in Figure 18.

**Figure 17.** The temperature values that were derived from the IR sensor measurements. **Figure 17.** The temperature values that were derived from the IR sensor measurements. efficiency and water-saving at all stages of plant growth. The spatial estimates of the CWSI

**Figure 18.** Spatial estimates (isolines of CWSI) from the first pilot flight with AMMS. **Figure 18.** Spatial estimates (isolines of CWSI) from the first pilot flight with AMMS.

The calibration resulting from the spring potato cultivation (non-stressed lower base-

the first pilot approach of the spatial CWSI in the autumn potato cultivation. The high relative humidity that prevailed during the calibration season, in combination with the low temperatures, gave a limited range of vapor pressure deficit of the atmosphere (VPD), leading to a high slope of the lower non-stressed baseline bottom line with a relatively high dispersion of the observed measurements (r2 = 0.613). However, this calibration provided satisfactory values for the spatial CWSI estimates in conjunction with the actual water loss from the ground profile just before the next irrigation (Figure 16 and Table 2).

The calibration resulting from the spring potato cultivation (non-stressed lower baseline: *T<sup>c</sup>* − *T<sup>a</sup>* = −1.74 × VPD − 1.23, stressed upper baseline *T<sup>c</sup>* − *T<sup>a</sup>* = 2.32 ◦C) was used in the first pilot approach of the spatial CWSI in the autumn potato cultivation. The high relative humidity that prevailed during the calibration season, in combination with the low temperatures, gave a limited range of vapor pressure deficit of the atmosphere (VPD), leading to a high slope of the lower non-stressed baseline bottom line with a relatively high dispersion of the observed measurements (r<sup>2</sup> = 0.613). However, this calibration provided satisfactory values for the spatial CWSI estimates in conjunction with the actual water loss from the ground profile just before the next irrigation (Figure 16 and Table 2).

#### **5. Conclusions**

The integration of new, improved methods and tools in modern agriculture, combined with the limitations of the water availability for irrigation, pushes research into new and smart approaches. A major advantage in this context is that, nowadays, it is both feasible and affordable to use proven and reliable methods that required expensive instruments and high installation costs in the past. Furthermore, in agricultural production, the involvement of irrigation water is crucial in maintaining the quality and conservation of it as a natural resource.

Nowadays, the use of UAVs in agriculture is universal. The continuous evolution of UAV technology with lighter materials and the evolution of battery technology allows for long flight periods with heavier equipment loads. However, spatial data from crops are mostly acquired only through multispectral and thermal cameras to date.

Taking advantage of the UAV evolution, GreenWaterDrone (GWD) proposes an innovative system, integrating aerial measurements of canopy temperatures for the near-realtime calculation of crop irrigation needs (crop water stress detection) and dynamic crop surveillance, achieving measurements with high temporal and spatial resolution, while at the same time, minimizing the intervention in crop activities and system maintenance costs. The GWD solution is especially suited for relatively small fields with increased crop diversity that cannot be efficiently addressed by satellite remote sensing. An onboard aerial micrometeorological system (AMMS) allows for calculating the known CWSI (crop water stress index) in spatial estimation for many crops in the area. Using the thermal images in parallel with spatial measurements by the infrared radiometers is an advantage that allows us to compare and calibrate the methodology at an accurate level. Furthermore, it must be noted that the estimation of the CWSI at the implementation level is achieved by the spatial estimates of the infrared thermometers from UAVs (AMMS). Moreover, it must be mentioned that the multispectral and photogrammetric cameras significantly assist the results by providing information concerning the canopy cover and field conditions useful for crop surveillance.

In this work, we present the GWD methodology, along with the equipment specifications used. All processes at all steps of applying the methodology on a pilot potato field are described, and preliminary results from the application are presented. From the pilot study, the functionality of the proposed system (GWD) was certified (accuracy of the UAV path and flight altitude, reliability of the aerial data acquisition system, communication stability between UAVs and ground base). Our findings indicated that the canopy temperatures derived from the ground meteorological station (GMMS), the AMMS, and the portable IRT radiometers produced a suitable thermal image from the surface of the crop. However, it should be noted that the scheduled flight took place in the middle of the autumn season when the crop evapotranspiration rate was very low (E*T<sup>c</sup>* was 2.7 mm/d and the soil moisture profile remained almost unchanged).

The project subsystems can be useful for supporting applications that are significant for irrigation water management and programming, such as irrigation alerting and scheduling, crop surveillance, and irrigation water management. However, more efforts are necessary to make these technologies more user-friendly and available for all end users, covering different advantages for a precise crop water stress evaluation. Still, the stakeholders, farmers, and industry responded positively to this effort's augmented interest. Furthermore, similar pilot campaigns need to be conducted in other regions on different crops to enhance the applicability of the GWD system. There is ongoing work by the authors to collect data for calibrating the methodology for other crops.

**Author Contributions:** Conceptualization, S.A., E.P., N.P., P.P., I.C. and G.K.; methodology, S.A., E.P., N.P., I.C. and P.P.; software, S.A., E.P., N.P. and I.C.; validation, S.A. and E.P.; formal analysis, S.A., E.P. and N.P.; resources, S.A., E.P., N.P. and I.C.; data curation, S.A., E.P., S.V., E.-M.P., P.P. and A.P.; writing—original draft preparation, S.A., E.P., N.P. P.P., I.C., G.K., S.V. and E.-M.P.; writing—review and editing, S.A., E.P., N.P. and P.P.; supervision, S.A. and E.P. All authors read and agreed to the published version of the manuscript.

**Funding:** This research was co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship, and Innovation under the call RESEARCH—CREATE—INNOVATE (project code: T1EDK-04021).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable for this article.

**Acknowledgments:** We want to thank all the Department of Agriculture, Economics, and Veterinary Medicine of Trifylia, and especially the director of the Department, Antonios Paraskevopoulos, for their tremendous support throughout the project GreenWaterDrone. Furthermore, we want to acknowledge the farmers Nikos Xirokostas, Charalabos Papadopoulos, and George Kalianis for their collaboration in our experiments in their field.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

