Remote Sensing Analysis of Surface Temperature from Heterogeneous Data in a Maize Field and Related Water Stress

: Precision agriculture aims at optimizing crop production by adapting management actions to real needs and requires that a reliable and extensive description of soil and crop conditions is available, that multispectral satellite images can provide. The purpose of the present study, based on activities carried out in 2019 on an agricultural area north of Ravenna (Italy) within the project LIFE AGROWETLANDS II, is to evaluate the potentials and limitations of freely available satellite thermal images for the identiﬁcation of water stress conditions and the optimization of irrigation management practices, especially in agricultural areas and wetlands a ﬀ ected by saline soils and salt water capillary rise. Point ﬁeld surveys and a very-high resolution thermal survey (5 cm) by an unmanned aerial vehicle (UAV) supported thermal camera were performed on a maize ﬁeld tentatively at every Landsat-8 passage to check land surface temperature (LST) and canopy cover (CC) estimated from satellite. Temperature measured in the soil near ground surface and from UAV


Introduction
Precision agriculture, which designates a combination of technologies and practices aiming to optimize the production by accounting for variability and uncertainties within agricultural systems [1], is emerging as a possible answer to the challenge of increasing food production in a sustainable way [2]. In practice, precision farming tries to supply variable rates of farming inputs (e.g., water and nutrients), depending on local soil and plant conditions [3]. Therefore, a precise knowledge of crop deficiency and quality of soils is required for the delineation of management zones, the optimization of farming practices and, possibly, the prediction of crop yield within-field variability [4]. For this purpose, a near real-time monitoring is needed with cost-effective technologies.
Remote sensing is certainly the most cost-effective method for large-scale monitoring [5] while, to perform analyses at the scale of a single crop, a proper combination of satellite and unmanned aerial vehicle (UAV) technologies is advisable [6]. On one hand, the use of both Landsat-8 and Sentinel-2 satellite platforms offers the opportunity to acquire multispectral images with appropriate resolution between 10 and 30 m and a revisiting time of a few days [7]. On the other hand, UAV technologies provide very high-resolution capabilities with the maximum timing flexibility, even though they are suitable for surveys on smaller extent areas.
Remote sensing can provide accurate field mapping for plant stress monitoring, biomass or nutrient estimation, weed detection and inventory [6]. Retrieving variables from remote sensing data requires the inversion of radiative transfer models. Furthermore, transfer outputs require a not straightforward interpretation in order to provide useful information for crop management [2]. As a consequence, obvious and widely recognized standard solutions suitable for all cases are not available. Methods to retrieve the variables of interest can be purely empirical (generally consisting in a regression between the measured signal and the variable of interest), mechanistic (based on physical model inversion) or contextual (e.g., classification techniques) [2].
When monitoring plant stress, a key aspect is the variation of the crop evapotranspiration [8]. Several indices have been proposed in literature that can be based on optical band ratios or on estimates of surface temperature. While the optical domain is considerably exploited, the use of thermal sensors in agriculture is less frequent [9], even though it is proven to be effective in the assessment of water stress [10]. In fact, the surface temperature of vegetation canopy is closely related to the transpiration rate, which is in turn controlled by the environment evaporative demand and the availability of water in the soil. Furthermore, surface temperature is a rapid response variable, if compared to indices in the visible and near-infrared domain [5] reflecting canopy status.
Thermal remote sensing from satellite platforms is still limited in spatial resolution and revisiting time. The highest resolution is provided by Landsat-7 thermal band, with a ground sample distance (GSD) of 60 m, whilst the more recent Landsat-8 provides two thermal infrared (TIR) bands at 100 m resolution; in both cases the revisiting time is about 16 days. Unfortunately, the presence of data gaps (black stripes) in the area under consideration due to a failure of the Scan Line Corrector, represents a major obstacle to the use of Landsat-7 Enhanced Thematic Mapper Plus (ETM+) images. Due to these limitations, the use of thermal sensors installed on board UAV platforms becomes particularly interesting for precision farming, even though their accurate calibration and the correction of atmospheric effects are still a challenge [11].
Monitoring crop growth and water stress was a major objective of the supporting project LIFE AGROWETLANDS II (LIFE15 ENV/IT/000423); it was originally (2015) proposed to monitor water stress by UAV over a couple of fields, but was later decided to perform this monitoring action over extended areas in extensive form through satellite images combined with UAV survey over a restricted Remote Sens. 2020, 12, 2506 3 of 31 area. Carlson [12] and Liou and Kar [13] present reviews of methods for estimating land surface temperature (LST), soil moisture, and actual evapotranspiration (ETa) from satellite imagery; the project decided to use in particular the Surface Energy Balance Algorithm for Land (SEBAL, [14]) whose main output is the estimate of ETa, obtained from the surface energy balance, in which surface temperature plays an essential role causing long-wave outgoing radiation that is measured by Landsat Thermal Infrared Sensor (TIRS). LST is both an intermediate result of the ETa estimate procedure and a variable from which crop water stress can be derived. A check of the resulting temperature was therefore planned that motivated the UAV thermal survey. It was considered useful to plan at the same time a chromatic optical survey, since it provides complementary information with higher resolution.
The present paper aims at describing the UAV survey characteristics and calibration, the comparison between the UAV survey and satellite-derived temperatures over a maize field north of Ravenna (Italy), providing a cross-validation of the two approaches and evidence of their primary merits, as well as an assessment of the potential of freely available satellite imagery to estimate evapotranspiration and crop water stress in the farmer perspective of precise agriculture implementation.
The present paper is structured as follows. Section 2 describes the study site, the crop and field surveys, UAV and satellite imagery acquisition as well as ground truth information. Section 3 presents methods followed in the analysis of UAV and satellite images in order to obtain temperature, evapotranspiration, and water stress values. Section 4 presents the results of the analysis and comparisons for cross-validation. Section 5 presents a brief discussion about the utilization of results for irrigation management and perspectives. Section 6 summarizes the conclusions.

Study Site
The study site is located in the Adriatic coastal plain, between the outlets of the Reno (44 • 35 N) and Lamone (44 • 32 N) rivers, approximately 18 km north of the city of Ravenna in the Emilia-Romagna region, Northern Italy ( Figure 1).
The low-lying Ravenna coastal area has been historically affected by a widespread land subsidence process due to the superposition of natural processes, mainly the compaction of Quaternary deposits, and anthropogenic activities related to aquifer overexploitation, mechanical drainage, and gas withdrawal from a number of deep reservoirs scattered inland and offshore [15,16]. A recent survey of subsidence in the Emilia-Romagna plain carried out by the Regional Agency for Prevention, Environment and Energy (ARPAE) through interferometric data processing has quantified the subsidence of the study area over the period 2011-2016 being around 2-3 mm/year [17]. Combined with sea level rise, land subsidence significantly increases the exposure to flood risk of coastal areas and favors saltwater intrusion into the shallow coastal aquifer [18], seriously threatening groundwater quality, fertility and productivity of agricultural soils, and wetland biodiversity.
A field was selected for the paper investigation in the Biomarcabò unit which is managed by Agrisfera Cooperative according to the principles of organic farming. The field is located about 2.5 km to the west of the Adriatic coastline and is bordered to the north by the artificial embankment of the Reno river. The field is subdivided by a dirt road into two land parcels: the north one has an extension of about 16.5 ha and the south one of 12.2 ha. The mean elevation of the field is approximately 0.30-0.40 m above mean sea level (m a.m.s.l.). In the summer of 2018 the ditches of the field were replaced by a subsurface drainage system consisting of perforated corrugated polyvinyl chloride (PVC) pipes placed 0.80 m below the ground surface at 10 m interval, with slope 1:1000 towards the interceptor ditch, located at the southern edge of the two parcels. By lowering the water table, removing the excess water and leaching out the salt solutions from the upper horizons of the soil, the drainage pipes actually separate the upper agricultural soil from lower layers, where brackish or saline water is normally present [19]. During the crop growing season soil and groundwater characteristics were monitored at the node P09 ( Figure 1) by the wireless sensor network installed within the project. The low-lying Ravenna coastal area has been historically affected by a widespread land subsidence process due to the superposition of natural processes, mainly the compaction of Quaternary deposits, and anthropogenic activities related to aquifer overexploitation, mechanical drainage, and gas withdrawal from a number of deep reservoirs scattered inland and offshore [15,16]. A recent survey of subsidence in the Emilia-Romagna plain carried out by the Regional Agency for Prevention, Environment and Energy (ARPAE) through interferometric data processing has quantified the subsidence of the study area over the period 2011-2016 being around 2-3 mm/year [17]. Combined with sea level rise, land subsidence significantly increases the exposure to flood risk of coastal areas and favors saltwater intrusion into the shallow coastal aquifer [18], seriously threatening groundwater quality, fertility and productivity of agricultural soils, and wetland biodiversity.
A field was selected for the paper investigation in the Biomarcabò unit which is managed by Agrisfera Cooperative according to the principles of organic farming. The field is located about 2.5 km to the west of the Adriatic coastline and is bordered to the north by the artificial embankment of the Reno river. The field is subdivided by a dirt road into two land parcels: the north one has an extension of about 16.5 ha and the south one of 12.2 ha. The mean elevation of the field is approximately 0.30-0.40 meters above mean sea level (m a.m.s.l.). In the summer of 2018 the ditches of the field were replaced by a subsurface drainage system consisting of perforated corrugated polyvinyl chloride (PVC) pipes placed 0.80 m below the ground surface at 10 m interval, with slope 1:1000 towards the interceptor ditch, located at the southern edge of the two parcels. By lowering the water table, removing the excess water and leaching out the salt solutions from the upper horizons of the soil, the drainage pipes actually separate the upper agricultural soil from lower layers, where Soil properties in the area are extremely variable reflecting the sedimentary response of coastal progradation in the form of alternating littoral dunes and interdunal depressions, filled with recent fine sediments, mainly silt and clay, deposited in fresh swamp or lagoon and salt marsh environments [20].
The climate of the study area is Mediterranean, characterized by mild wet winters and warm to hot dry summers [21]. The climate in the area features an average annual precipitation of 696 mm/year (estimated by ARPAE [22] over the period 1991-2015) with a minimum in the summer season and two weak peaks, the main one in autumn and the secondary one in spring. The mean annual temperature is 14.1 • C with the coldest period in January, when the monthly average temperature is 3.4 • C, and the hottest period in July with 24.2 • C. Rainfall in summer is episodic and highly variable from year to year, so that a few irrigation events are normally necessary to obtain high yield from summer crops.
The area was the object of surveys and monitoring during the project period, autumn 2016 to spring 2020. The monitoring action referred in this paper was carried out during the summer crops growing season in 2019. It is highlighted that, after a very dry first quarter and a normal rainfall in April, May was characterized by exceptionally high precipitations and by temperatures much lower than normal [23]. Furthermore, exceptionally high temperatures, with maximum value close to 40 • C, occurred in the last week of June due to a high pressure system located over Europe, accompanied by hot, stable, and dry air intrusion of Saharan origin [24]. Local weather data are recorded by the two agro-meteorological stations installed in the study area at nodes P02 and P07 ( Figure 1).
Due to the aforementioned unfavorable weather conditions, sowing of maize in the field under investigation was delayed until May 23. Maize (Pioneer P1517W, FAO class 600, 128 days growing period to physiological maturity) was sown with row spacing of 0.75 m and seed spacing of 0.15 m Remote Sens. 2020, 12, 2506 5 of 31 along the north-south oriented rows. The maize plants started flowering in the penultimate week of July and were harvested on August 24 at the dough stage, 94 days after sowing. The plants were cut at a height of about 20-30 cm from the ground and chopped to produce maize silage.
Field observations were performed about once a week to monitor maize phenological development. Maize growth stages were classified according to the BBCH (Biologische Bundesanstalt, Bundessortenamt und Chemische Industrie) scale [25] with its associated decimal coding system.
The field plot measurements were integrated with satellite images to explore within-field spatial variability in crop growth and final yield. Furthermore, optical and thermal images obtained from UAV for a sample strip were used for field validation of satellite remote sensing.

Imagery Surveyed from UAV and Satellites
The reference area surveyed from sensors mounted on the UAV consists of a 64 m wide and 780 m long strip oriented from S-W to N-E ( Figure 2). The total surveyed area is therefore equal to 5 ha and covers two different kind of crops: maize in the N-E field and soybean in the circular S-W field. This paper aims to study the first crop. The area was in fact carefully selected as an optimal compromise between the UAV's flight autonomy and the coverage of a large portion of the studied crops and habitats. Interferences of the UAV flight with other activities in the field area were also taken into account, besides the compliance of the flight and the necessary authorizations. The UAV reference area is the portion of territory where image cover and a planimetric precision of 5 cm are guaranteed thanks to methods described in the following section.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 32 Reference area for UAV surveys is highlighted with a red polygon. The Landsat-8 grid of 30 × 30 m pixels at ground level in path 192 is overlaid in light blue. The 12 selected tiles discussed in this paper are highlighted. The orthophoto is produced from images acquired on July 19 accurately geo-referenced with 10 ground control points (GCPs) thermal target (a). Grid coordinates are UTM-ETRS89 Zone 33T. The Reno river is visible at the NE end of the strip.
By using a dual-gimbal UAV, it was possible to acquire imagery in one single flight with both sensors necessary for the experimental study over the surveyed area. It was therefore possible to optimize the use of the UAV battery and to perform the flight in the shortest time window.
In fact, the defined area was surveyed from a UAV in precise dates during the whole agricultural season from late spring (16 May 2019, first survey S01) to end of summer 2019 (4 September 2019, last survey S09; Figure 3). The dates were selected in order to be as close as possible The UAV that met the described requirements is an industrial grade quadri-copter, model Matrice 210 from manufacturer DJI. This UAV has a dual-gimbal configuration that can handle the simultaneous use of two sensors' payloads.
By using a dual-gimbal UAV, it was possible to acquire imagery in one single flight with both sensors necessary for the experimental study over the surveyed area. It was therefore possible to optimize the use of the UAV battery and to perform the flight in the shortest time window.
In fact, the defined area was surveyed from a UAV in precise dates during the whole agricultural season from late spring (16 May 2019, first survey S01) to end of summer 2019 (4 September 2019, last survey S09; Figure 3). The dates were selected in order to be as close as possible to Landsat-8 overpasses. Thanks to the availability of information on the satellite orbit, it was possible to anticipate each satellite overpass and plan the UAV surveys accordingly. The surveys were then confirmed accounting for weather forecasts and carried out only in clear-sky conditions. × 30 m pixels at ground level in path 192 is overlaid in light blue. The 12 selected tiles discussed in this paper are highlighted. The orthophoto is produced from images acquired on July 19 accurately geo-referenced with 10 ground control points (GCPs) thermal target (a). Grid coordinates are UTM-ETRS89 Zone 33T. The Reno river is visible at the NE end of the strip. By using a dual-gimbal UAV, it was possible to acquire imagery in one single flight with both sensors necessary for the experimental study over the surveyed area. It was therefore possible to optimize the use of the UAV battery and to perform the flight in the shortest time window.
In fact, the defined area was surveyed from a UAV in precise dates during the whole agricultural season from late spring (16 May 2019, first survey S01) to end of summer 2019 (4 September 2019, last survey S09; Figure 3 The dates were selected in order to be as close as possible to Landsat-8 overpasses. Thanks to the availability of information on the satellite orbit, it was possible to anticipate each satellite overpass and plan the UAV surveys accordingly. The surveys were then confirmed accounting for weather forecasts and carried out only in clear-sky conditions. In order to obtain an accurately geo-referenced survey, 10 topographic ground control points (GCPs) were properly displaced covering the whole study area. Appropriate positioning has been crucial for consistent results, in particular since the study area is elongated [26]. GCPs placement was also challenging due to the search for undisturbed locations during the whole season without interfering with farming activities carried out in the fields. Each GCP location was then identified, marked, and precisely located with GNSS NRTK instrument.

UAV Thermal and Optical Surveys
The first payload of the dual-gimbal UAV platform is a thermal camera, with technical specifications properly chosen for the experimental study. The thermal sensor acquires brightness temperature measurements within the spectral band from 7.5 to 13.5 μm. The camera is a DJI Zenmuse XT with resolution equal to 0.3 MP, pixel size of 17 μm, delivering thermal imagery with 640 × 512 pixels. The chosen focal length for the survey is equal to 13 mm.
The flight was planned with a GSD of approximately 12.5 cm, flying at an average height of 95 m above the surveyed surface.
Each flight was planned with a minimum overlap of 80% between subsequent images and a sidelap of at least 65% between each track, with heading perpendicular to the main strip dimension. In order to avoid blur movements in the acquired imagery, the UAV moved at a reduced speed and each thermal image was captured only when the UAV was hovering, avoiding any undesired movement.
In order to be visible from the UAV platform, each GCP on the ground consisted of a custom-made target, created using highly reflecting aluminum foils. Therefore, each target offered enough contrast in the thermal band, compared to the ground temperature, becoming highly visible in the surveyed imagery. Furthermore, the target overall dimensions (30 × 30 cm) were appropriately In order to obtain an accurately geo-referenced survey, 10 topographic ground control points (GCPs) were properly displaced covering the whole study area. Appropriate positioning has been crucial for consistent results, in particular since the study area is elongated [26]. GCPs placement was also challenging due to the search for undisturbed locations during the whole season without interfering with farming activities carried out in the fields. Each GCP location was then identified, marked, and precisely located with GNSS NRTK instrument.

UAV Thermal and Optical Surveys
The first payload of the dual-gimbal UAV platform is a thermal camera, with technical specifications properly chosen for the experimental study. The thermal sensor acquires brightness temperature measurements within the spectral band from 7.5 to 13.5 µm. The camera is a DJI Zenmuse XT with resolution equal to 0.3 MP, pixel size of 17 µm, delivering thermal imagery with 640 × 512 pixels. The chosen focal length for the survey is equal to 13 mm.
The flight was planned with a GSD of approximately 12.5 cm, flying at an average height of 95 m above the surveyed surface.
Each flight was planned with a minimum overlap of 80% between subsequent images and a sidelap of at least 65% between each track, with heading perpendicular to the main strip dimension. In order to avoid blur movements in the acquired imagery, the UAV moved at a reduced speed and each thermal image was captured only when the UAV was hovering, avoiding any undesired movement.
In order to be visible from the UAV platform, each GCP on the ground consisted of a custom-made target, created using highly reflecting aluminum foils. Therefore, each target offered enough contrast in the thermal band, compared to the ground temperature, becoming highly visible in the surveyed imagery. Furthermore, the target overall dimensions (30 × 30 cm) were appropriately set accordingly to the GSD of the planned survey. A white plastic square board with installed thermocouples was placed aside the target for thermal camera calibration.
The second payload of the dual-gimbal configuration was a DJI optical camera with resolution equal to 12 MP (4000 × 3000 pixels), with a pixel size of 1.6 µm. The focal length expressed in standard 35 mm equivalent sensor is 26 mm. Flying at an average height of 95 m, the resulting GSD for the optical camera is equal to 3 cm.

Acquisition of Satellite Data
The time series of images acquired by the Landsat-8 in the period from 16 May to 5 September 2019 were used, testing their consistency for monitoring of crop growth and with thermal UAV data.
In the present analysis Landsat-8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) Collection 1 level-1 images were used for the study. Landsat Collection 1 products are made publicly available by the U.S. Geological Survey (USGS) website: https://earthexplorer.usgs.gov/. The images are already processed to a level-1 precision and terrain corrected product (L1TP), that means radiometrically calibrated and orthorectified using GCPs and digital elevation model (DEM). The Landsat products are delivered in GeoTIFF format, using World Geodetic System (WGS) 84 datum and Universal Transverse Mercator (UTM) map projection (32 or 33 in the case under examination). The two thermal bands (band 10 and 11 with central wavelength 10,895 and 12,005 nm in the order) capture data with 100 m resolution but are delivered after resampling by cubic convolution to 30 m to match OLI multispectral bands. Landsat-8 level-2 data, i.e., atmospherically corrected surface reflectance data using the Land Surface Reflectance Code (LaSRC) [27] and Normalized Difference Vegetation Index (NDVI), provided by USGS, were used to evaluate surface emissivity for methods employing atmosphere effects corrections.
Technical features can be derived from official sensors handbooks [28,29], while image acquisition dates are presented in Figure 3, along with the crop growth stage, expressed in the BBCH scale, detected during the field surveys. Landsat-8 images acquired under clear-sky conditions were selected for the analysis.
Sentinel-2 images are marginally used in this paper, but due to the higher passage frequency, they gave extensive contextual information as field surveys provided detailed point-wise data.
A digital terrain model (DTM) is required by SEBAL. The DTM with spatial resolution of 5 m produced for the AGROWETLANDS II project was used. It covers the area of interest; it was derived from the National Cartographic Portal (http://www.pcn.minambiente.it/mattm/) and was integrated with the topo-bathymetric survey carried out in 2016 by ARPAE [30].

Ground Temperature Acquisition for Validation
A Sentek Drill & Drop probe 60 cm long was installed at the northern edge of the southern parcel for measuring soil properties during the growing season (P09 in Figure 1). Soil moisture, temperature, and electrical conductivity sensors were placed at depths of 0.05, 0.15, 0.25, 0.35, 0.45, and 0.55 m and measurements were recorded approximately every 15 min. Accounting for solar radiation, air conditions, and precipitation measured at the meteorological station P07 (Figure 1), temperature was hindcast at each sensor depth and extrapolated up to the soil surface with the numerical model HYDRUS-1D [31], representing 1D evolution of water phases, liquid and vapor, and of temperature through mass and heat balance equations.
Temperature was also monitored in several near water bodies: Valli di Comacchio at Bellocchio outlet (data kindly delivered by ARPAE), Valle della Canna outlet, and Canale Destra Reno in proximity of the sea outlet (data collected by AGROWETLANDS II project).

Irrigation Water Management
Water needed to fulfill the crop requirements for its full growth was supplied by above-canopy sprinkler irrigation, scheduled by using the Smart AGROWETLANDS Decision Support System Remote Sens. 2020, 12, 2506 8 of 31 (SA-DSS) implemented within the project. The SA-DSS is mainly constituted by a wireless sensor network for real-time monitoring of weather, soil, channel-and ground-water conditions [32] and by a forecast model of crop growth, for which the AquaCrop model developed by the Land and Water Division of FAO (Food and Agriculture Organization of the United Nations) [33,34] was adopted.
The field was irrigated twice; the first time on 3-7 July with 45 mm of water and the second time on 30 July-3 August with 50 mm. During the growing period rainfall contributed with 161 mm, the major amount concentrated in three events (42 mm around 27 May, 49 mm on 22 June, and 40 mm around 13 July).

UAV Imagery Pre-Processing
Both thermal and optical imagery datasets were further processed using Structure from Motion (SfM) techniques. In the first part of the SfM processing Global Positioning System (GPS) locations recorded by the UAV platform were used as approximate coordinates for the camera positions. Then, using coordinates obtained with greatest accuracy from GCP survey, the whole photogrammetric block was georeferenced and mosaicked. For each planned survey in the crop season both thermal and optical orthophotos were produced.
Having a good coverage and overlap between the whole dataset, at the end of the SfM processing, the final mosaicked output products were delivered with 0.05 m for thermal imagery and 0.02 m for the optical one, ready to be further processed in the experimental workflow. Figure 4 presents the thermal orthomosaic acquired on 19 July by the end of stem elongation stage.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 32 network for real-time monitoring of weather, soil, channel-and ground-water conditions [32] and by a forecast model of crop growth, for which the AquaCrop model developed by the Land and Water Division of FAO (Food and Agriculture Organization of the United Nations) [33,34] was adopted. The field was irrigated twice; the first time on July 3-7 with 45 mm of water and the second time on July 30-August 3 with 50 mm. During the growing period rainfall contributed with 161 mm, the major amount concentrated in three events (42 mm around May 27, 49 mm on June 22, and 40 mm around July 13).

UAV Imagery Pre-Processing
Both thermal and optical imagery datasets were further processed using Structure from Motion (SfM) techniques. In the first part of the SfM processing Global Positioning System (GPS) locations recorded by the UAV platform were used as approximate coordinates for the camera positions. Then, using coordinates obtained with greatest accuracy from GCP survey, the whole photogrammetric block was georeferenced and mosaicked. For each planned survey in the crop season both thermal and optical orthophotos were produced.
Having a good coverage and overlap between the whole dataset, at the end of the SfM processing, the final mosaicked output products were delivered with 0.05 m for thermal imagery and 0.02 m for the optical one, ready to be further processed in the experimental workflow. Figure 4 presents the thermal orthomosaic acquired on 19 July by the end of stem elongation stage. From the radiometric point of view, the images recorded by the thermal camera were corrected for the atmospheric effects using the model implemented in the FLIR ResearchIR software (version From the radiometric point of view, the images recorded by the thermal camera were corrected for the atmospheric effects using the model implemented in the FLIR ResearchIR software (version 4.40). At the same time, each image file format was converted from radiometric JPEG (RJPEG), as recorded in camera, to a Tagged Image File Format (TIFF). Air temperature and humidity were derived from the project meteorological stations. At a first step, the surface emissivity was assumed constant and equal to 0.95. Analyzing the histogram of these first surface temperature values, it is possible to recognize a bimodal distribution, roughly corresponding to vegetated and non-vegetated surfaces. In order to better separate the two populations, the mixed pixels were cancelled using the edge detection Canny method implemented in MATLAB environment.
The thermal camera was calibrated accounting for temperature measurements by thermocouples on white plastic surfaces used as thermal reference points (TRPs) on ground surface, having measured emissivity (0.84-0.86). The thermal camera calibration factor resulted equal to 1.049, leading to a significant temperature overestimation in the absence of compensation. The two modal temperatures were further corrected accounting for more realistic emissivity values of the two surface classes (i.e., 0.98 for vegetation and 0.95 for bare soil). Calibration produces a temperature reduction of about 4.5 • C, of which 3 • C are due to thermal camera and the rest to soil and vegetation emissivity.
For the UAV surveys on bare soil (16 May and 4 September), the obtained surface temperatures are in good agreement with those derived from field sensors in the soil and interpreted with HYDRUS model. Figure 5 shows the simulation results starting from 11 May. The detail presented for Landsat and UAV passage on 16 May shows that, despite the short delay, the bare soil surface temperature increased 2.5 • C between the two passages under the effect of strong sunshine. The air temperature over the same time interval varied only 0.5 • C and was 14 • C below ground surface temperature.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 32 4.40). At the same time, each image file format was converted from radiometric JPEG (RJPEG), as recorded in camera, to a Tagged Image File Format (TIFF). Air temperature and humidity were derived from the project meteorological stations. At a first step, the surface emissivity was assumed constant and equal to 0.95. Analyzing the histogram of these first surface temperature values, it is possible to recognize a bimodal distribution, roughly corresponding to vegetated and non-vegetated surfaces. In order to better separate the two populations, the mixed pixels were cancelled using the edge detection Canny method implemented in MATLAB environment. The thermal camera was calibrated accounting for temperature measurements by thermocouples on white plastic surfaces used as thermal reference points (TRPs) on ground surface, having measured emissivity (0.84-0.86). The thermal camera calibration factor resulted equal to 1.049, leading to a significant temperature overestimation in the absence of compensation. The two modal temperatures were further corrected accounting for more realistic emissivity values of the two surface classes (i.e., 0.98 for vegetation and 0.95 for bare soil). Calibration produces a temperature reduction of about 4.5 °C, of which 3 °C are due to thermal camera and the rest to soil and vegetation emissivity.
For the UAV surveys on bare soil (16 May and 4 September), the obtained surface temperatures are in good agreement with those derived from field sensors in the soil and interpreted with HYDRUS model. Figure 5 shows the simulation results starting from 11 May. The detail presented for Landsat and UAV passage on 16 May shows that, despite the short delay, the bare soil surface temperature increased 2.5 °C between the two passages under the effect of strong sunshine. The air temperature over the same time interval varied only 0.5 °C and was 14 °C below ground surface temperature. Despite efforts to reduce UAV flight in the shortest time window, UAV surveys and Landsat-8 passage are not perfectly synchronized. Therefore, UAV measurements are corrected accounting for soil surface and air temperatures time variation. Corrections are presented in Table 1. The resulting expected root mean square (r.m.s.) error for synchronous UAV values is estimated to be around 1.5 °C: 1.2 °C for the calibrated thermo-camera residual error and 0.8 °C for synchronization correction, assumed to be statistically independent.

Relation Between UAV Measured LST and Landsat Derived One
The relation is assessed by comparing surface temperature value at the Landsat pixels with the UAV derived temperature averaged over the corresponding 30 × 30 m tiles (see Figure 2 for tiles Despite efforts to reduce UAV flight in the shortest time window, UAV surveys and Landsat-8 passage are not perfectly synchronized. Therefore, UAV measurements are corrected accounting for soil surface and air temperatures time variation. Corrections are presented in Table 1. The resulting expected root mean square (r.m.s.) error for synchronous UAV values is estimated to be around 1.5 • C: 1.2 • C for the calibrated thermo-camera residual error and 0.8 • C for synchronization correction, assumed to be statistically independent. Table 1. Temperature statistics over tiles derived from UAV and Landsat-8 using single band method uncorrected for atmospheric effects (SB), radiative transfer equation (RTE), and split window method (SW). Statistical results of thermal surveys of the study area. Average (up row) and standard deviation (low) over 12 tiles.

Relation between UAV Measured LST and Landsat Derived One
The relation is assessed by comparing surface temperature value at the Landsat pixels with the UAV derived temperature averaged over the corresponding 30 × 30 m tiles (see Figure 2 for tiles identification). Comparison is restricted to 12 tiles corresponding to the central area of the maize field ( Figure 2) where irrigation scheduling under saline conditions was implemented by using the DA-DSS. Parts of the UAV derived orthomosaic contained in the selected tiles are identified and statistics over the tiles are evaluated.
Landsat TIRS data are distributed with 30 m resolution, even if original measurements are obtained with 100 m resolution; temperatures in adjacent pixels are therefore not independent. Selected tiles are scattered in the area to reduce the effect of non-independence and statistics over the set of tiles disregard any effect of dependence.
The temperature analysis of UAV thermal survey is carried out to obtain vegetation cover and the typical temperatures of vegetation and soil. The average surface temperature of the tile (T LS ) is evaluated by definition (arithmetic average of all pixel values). Typical temperatures of vegetation and soil were obtained by cancelling the mixed pixels and obtaining in such way a clearly bimodal distribution of temperatures over the tile (see Figure 6) from which the optimal threshold separating cold vegetation from hot soil can be easily identified, 45 • C in Figure 6. This temperature value is then used to count pixels below and above threshold, and hence to derive fractional vegetation cover (F V ) and bare soil fraction (F S ), with F V + F S = 1. Typical temperatures of vegetation and soil are preliminarily assumed equal to the two modal or peak values; average vegetation and soil temperatures, T V and T S , are then derived so that they conform to the measured average surface temperature, Equation (1) where T LS = LST, with minimal surface weighted discrepancy from modal values T VP and T SP , i.e., by minimizing F S (T S -T SP ) 2 + F V (T V -T VP ) 2 , resulting in the correction Equations (2) and (3): values TVP and TSP, i.e. by minimizing FS(TS -TSP) 2 + FV(TV -TVP) 2 , resulting in the correction Equations (2) and (3):

TV = TLS + FS (TVP -TSP).
(3) Figure 6. Image of the tile (left) and corresponding temperature distribution before calibration and after edge removal.
It is worth to remark that the tile image ( Figure 6, left) contains much more information than the temperature distribution, information relative to the spatial structure of the crop. Rows are easily identified, as well as a row oriented structure of dense/sparse canopy cover that suggests that sowing may be the cause of yield deficit. This information is lost at the tile scale.
Temperatures  It is worth to remark that the tile image ( Figure 6, left) contains much more information than the temperature distribution, information relative to the spatial structure of the crop. Rows are easily identified, as well as a row oriented structure of dense/sparse canopy cover that suggests that sowing may be the cause of yield deficit. This information is lost at the tile scale.
Temperatures presented in Figures 4 and 6 and used in Equations (1) In this study land surface temperature is retrieved from Landsat-8 satellite data according to three different methods: 1) by using Artis and Carnahan equation [35] as implemented in the Surface Energy Balance Algorithm for Land (SEBAL) described in Section 3.3.1; 2) the radiative transfer equation based method; and 3) the split window algorithm.
The simplest method 1) was implemented first with good results on water bodies' temperatures, but when applied to land surfaces a significant gap remained between satellite and UAV derived temperatures that has to be related to atmosphere effects. The complete and enlightening review [36] published at Landsat-8 launch time, and the specific paper on the subject [37] guided us to the choice of the two atmosphere correction methods better described in Sections 3.2.2 and 3.2.3.

LST from Artis and Carnahan Equation
In the first method surface temperature is derived from Landsat-8 TIRS band 10 by using Artis and Carnahan equation [35] as implemented in the SEBAL model which we used for estimation of actual evapotranspiration [38,39], not accounting for atmosphere correction parameters but only for surface emissivity: where T B is the satellite brightness temperature, λ is the thermal band central wavelength, and ε NB is the narrow band surface emissivity, representing the surface emissivity within the band range of the satellite thermal sensor. c 2 = hc/k B where h is the Planck constant, c is the velocity of light in vacuum and k B is the Boltzmann constant. Following Landsat-8 Handbook [28], determination of top of atmosphere brightness temperature is a two-step process. In the first step level-1 digital number (DN) values of Landsat-8 band 10 are converted to at-sensor spectral radiance values through Equation (5): where M L is the radiance multiplicative scaling factor for the given band provided in the image metadata file, A L is the radiance additive scaling term for the given band, available in the same file, and Q cal refers to level-1 pixel value in DN. Then TIRS data can be converted from spectral radiance to brightness temperature, which is the effective temperature viewed by the satellite under the assumption of unity emissivity, by using the formula where L sen λ is the top of atmosphere spectral radiance and K 1 and K 2 are band specific thermal conversion constants available from the image metadata file.
The following empirical equations are used in SEBAL model to evaluate ε NB from OLI reflectance data as a function of NDVI and Leaf Area Index (LAI) [40]. where NDVI > 0: ε NB = min(0.97 + 0.0033 · LAI, 0.98).
For water, where NDVI ≤ 0, ε NB = 0.99. NDVI is an indicator of the amount and condition of green vegetation. In SEBAL model NDVI is normally calculated from the top of atmosphere reflectance ρ of the red and near-infrared bands; for Landsat-8 bands 4 and 5 in the order: LAI is defined as the ratio of the total area of all leaves on a plant to the ground area represented by the plant. It is an indicator of biomass and canopy resistance to vapor flux and is computed using the following empirical equation [40]: where the Soil Adjusted Vegetation Index (SAVI) is calculated as follows: and L SAVI is a constant (a value of 0.5 is used). This single band method for retrieving LST from Landsat-8 TIRS data uncorrected for atmospheric effects is denoted as SB in what follows and T SB is the derived LST (T s in Equation (4)).

LST from the Radiative Transfer Equation Based Method
The second method used to retrieve LST from the Landsat-8 TIRS band 10 is based on the inversion of the radiative transfer equation. According to Ottle and Stoll [41], a simplified radiative transfer equation that indicates the radiance measured by a remote sensor at satellite level for a given wavelength λ can be expressed as: where λ is the given wavelength; L sen λ is the at-sensor radiance; B λ (T s ) is the emitted radiance by a black body at temperature T s ; τ λ is the atmosphere transmittance from the surface to the top; ε λ is the land surface emissivity; L ↓ λ is the downwelling atmospheric radiance and L ↑ λ is the upwelling atmospheric radiance. According to Planck law, the emitted radiance for a black body at temperature T s can be expressed in the form Remote Sens. 2020, 12, 2506 13 of 31 with c 1 = 2hc 2 (1.19104 × 10 8 W µm 4 m −2 sr −1 ) and c 2 = hc/k B (1.43877 × 10 4 µm K). Therefore, T s can be obtained from B λ (T s ) as in Equation (13): At-sensor spectral radiance L sen λ is obtained for Landsat-8 band 10 as indicated in Equation (5). In this study the atmospheric parameters needed for retrieving land surface temperatures, τ λ , L ↓ λ , and L ↑ λ , are calculated for each Landsat-8 scene by using the National Aeronautics and Space Administration (NASA) Atmospheric Correction Parameter Calculator developed by Barsi et al. [42,43], and available online at https://atmcorr.gsfc.nasa.gov/. A NDVI threshold method is used in this context to estimate land surface emissivity from Landsat-8 satellite data, based on the following equation, as suggested by Sobrino et al. [44]: where a λ and b λ are channel dependent regression coefficients and ε Vλ and ε Sλ are respectively the vegetation and soil emissivities at wavelength λ, estimated by using the values proposed by Skoković et al. [45]. Landsat-8 level-2 data and surface reflectance derived NDVI were used in ε λ evaluation. NDVI V and NDVI S represent the NDVI values corresponding to full vegetation and bare soil and were set at 0.65 and 0.15, in the order. C λ is a term which takes into account the cavity effect due to surface roughness. It is evaluated by using the expression suggested in [46]: where F is a geometrical factor ranging between 0 and 1, depending on the geometrical distribution of the surface; a mean value of 0.55 was assumed according to Cristóbal et al. [46] and Yu et al. [47]. P V refers to the vegetation fraction and is calculated as follows: P V is then set to 0 for pixels with NDVI < NDVI S and set to 1 for pixels with NDVI > NDVI V . This method is denoted as RTE in what follows and T RTE the derived LST (T s in Equation (13)).

LST from the Split Window Algorithm
The third used method is based on the split window algorithm proposed by Jiménez-Muñoz et al. [48]. It uses data of both TIR bands and derives from the difference of brightness temperature in the two bands an optimized corrective term to be applied to temperature derived from band 10. The split window algorithm used in this study has the following mathematical structure: where T i and T j are top of atmosphere brightness temperatures for the band i and j (in K) obtained through Equations (5) and (6). For Landsat-8, these two bands are band 10 and 11 respectively. ε m is the mean land surface emissivity of the two bands and is calculated as ε m = 0.5(ε i + ε j ); ∆ε is the emissivity difference given by ∆ε = (ε i − ε j ). w is the total atmospheric water vapor content (in g cm −2 ). c 0 to c 6 are coefficients derived by statistical linear regression performed on a simulated dataset of brightness temperatures over different surface and atmospheric conditions. In the split window algorithm given by Equation (17), the coefficients c 0 to c 6 take the following values as indicated in [48]: The NDVI threshold method described in the previous subsection for RTE method is also used in this case to derive the land surface emissivity from satellite data.
The total atmospheric water vapor content required for retrieving LST by this method was obtained for each Landsat-8 scene from the Moderate Resolution Imaging Spectroradiometer (MODIS)/Terra level-2 Total Precipitable Water Vapor product (MOD05_L2, freely available at https://ladsweb.modaps. eosdis.nasa.gov/search/), since the time of this satellite overpass over the study area is close enough to the Landsat one.
This method is denoted as SW in what follows and T SW the derived LST (T s in Equation (17)). The difference of temperature obtained from Landsat-8 data by using different methods is widely discussed in different papers [49,50] and is just documented in the present one. Our contribution is an analysis of the effect of survey resolution, as it is clear from Figure 6 that ground temperature is far from uniform on a 30 or 100 m grid element.

Actual Evapotranspiration and Crop Water Stress
The interest on temperature in agriculture is related to the fact that all biochemical processes depend on temperature and, specifically for irrigation management, that it is a variable necessary to evaluate evapotranspiration and crop water stress.
We have therefore finally made use of the Landsat derived temperatures to obtain actual evapotranspiration (ETa) and water stress (WS) maps by using SEBAL. ETa is obtained from the surface energetic balance with some information on the subdivision between latent and sensible heat fluxes. The innovative component of SEBAL is the use of a near-surface temperature difference dT related to the radiometric surface temperature T s by a linear relation determined by two extreme conditions: • a hot or dry pixel, representing zero ETa conditions, and • a cold or wet pixel, representing areas where ETa equals the potential value as LST equals air temperature canceling sensible heat flux convected to the air.
SEBAL is not autocalibrating and ETa maps depend on absolute LST maps, in the sense that if surface temperature is varied the energy balance items and in particular evapotranspiration are changed consequently.
METRIC [51,52], otherwise, is provided with an internal calibration fixing two anchor points where ETa can be evaluated, for instance, with the established and well-calibrated alfalfa-based formula for reference ET (ETr; ASCE-EWRI equation [53]; FAO paper 56 [54]). In any case maps will depend on how precisely anchor points satisfy conditions assumed for them, so that both approaches can be the best depending on contextual conditions. We have finally used the SEBAL model with the different methods described for retrieving LST from Landsat-8 data (Equations (4), (13) and (17)) to obtain actual evapotranspiration and water stress maps and verify how much maps depend on the temperature correction method to account for atmosphere effects.

Surface Energy Balance Algorithm for Land (SEBAL)
SEBAL is a satellite-based image processing model with the main objective to provide estimates of ETa at relatively high resolution from remotely sensed land surface temperatures. The major principles of the SEBAL model are described in detail in [14,55,56].
An instantaneous value of the ETa flux (kg m −2 s −1 or equivalently mm s −1 of liquid water) is calculated for every pixel of the satellite image at the overpass time as the residual of the surface energy balance equation: where λE is the latent heat flux and H is the sensible heat flux convected to the air, R n is the net radiation balance, and G is the sensible heat flux conducted into the ground. The latent heat flux represents the amount of energy absorbed to maintain a certain ET rate.

Net Radiation
The net radiation at the surface R n is estimated quantifying all the incoming and outgoing shortwave and longwave radiation components where R S↓ is the incoming shortwave radiation, R L↑ is the outgoing longwave radiation, R L↓ is the incoming longwave radiation, α is the surface albedo, and ε 0 is the broadband surface thermal emissivity.
The incoming shortwave radiation R S↓ representing the direct and diffuse solar radiation flux that actually reaches the Earth's surface, is the principal energy source for ET and can be measured or estimated assuming clear sky conditions, which is a prerequisite for a reliable evaluation from satellite images, by using the following equation: where G sc is the solar constant value (1367 W/m 2 ), β is the sun elevation angle contained in the metadata file of Landsat imagery data, and d r is the inverse squared relative distance between the Earth and the Sun depending on the day of the year. The albedo at the top of the atmosphere α TOA is evaluated as the weighted average of spectral reflectance over bands 1 to 7 of the Landsat-8 OLI: where ρ λ is the spectral reflectance of the band and ω λ is the fraction of the total solar irradiance in the considered band. The albedo at the top of the atmosphere α TOA is then adjusted for the broadband shortwave atmospheric transmissivity τ sw , as described in SEBAL manual [40], to obtain the surface albedo. τ sw is derived assuming clear sky conditions, in our case where land is at sea level τ sw = 0.75. The outgoing longwave radiation R L↑ , that is the thermal radiation flux emitted from the Earth's surface to the atmosphere, is calculated by using the Stefan-Boltzmann equation where ε 0 is the broadband surface emissivity, σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W/m 2 /K 4 ), and T s is the surface temperature in Kelvin. For surface temperature evaluation, in addition to the SB method, we tested the use of Equations (13) and (17) which provide T s values corrected for atmospheric effects.
The following empirical equation is used in SEBAL to evaluate the broadband surface thermal emissivity [40] on land, where NDVI > 0: For water, with NDVI ≤ 0, ε 0 = 0.985. The incoming longwave radiation is the downward thermal radiation flux from the atmosphere, which is calculated using the Stefan-Boltzmann equation where ε a is the effective atmospheric emissivity and T a is the near-surface air temperature. The following empirical equation is used for ε a [40]: ε a = 0.85(− ln(τ sw )) 0.09 (25) where τ sw is the broadband atmospheric transmissivity given above. Allen et al. [51], for calibration of METRIC, recommend to estimate T a in Equation (24) through the temperature of a "cold pixel" which is a pixel selected from a wet, well-irrigated crop surface with a complete vegetation cover (LAI > 4), since at that condition, H is generally small and near-surface air temperature and surface temperature can be assumed to be similar.

2.
Soil heat flux Soil heat flux (G) is the rate of heat storage in the soil and vegetation due to conduction. In SEBAL it is calculated through the ratio G/R n using the empirical equation proposed in [51,55] representing values near midday: where T s is the surface temperature (K) and α is the surface albedo. For water bodies (NDVI < 0), a mean value of 0.5 is assumed for the ratio G/R n .

Sensible heat flux
Sensible heat flux (H) is the flux of heat from the Earth's surface to the atmosphere and is estimated through similarity of turbulent diffusion of momentum and heat, i.e., by the aerodynamic resistance equation transposed to heat, as proposed in [51,55]: where ρ air is air density (kg m −3 ), C p is specific heat of air at constant pressure (1004 J kg −1 K −1 ), and r ah (s m −1 ) is the aerodynamic resistance to heat and/or momentum transfer between two near-surface heights, z 1 and z 2 (usually 0.1 and 2.0 m above the displacement plane), between which temperature difference is dT. The assumed wind velocity vertical distribution u(z) in the elevation range between surface and blending height z B , i.e., an elevation (usually 200 m) sufficient to retain wind velocity constant over the area and independent from local surface roughness is where k is the von Karman constant (0.41), u * is the friction velocity (m s −1 ), and z om is the roughness length for momentum transfer (m), that is calculated over the map as proposed in [57]: ψ m(z) is a correction term of the logarithmic velocity distribution accounting for stability conditions of the low atmosphere that is 0 for neutral stability and in general depends on the ratio between z and Monin-Obukhov length L [52]. The characteristic length scale L is defined as L = − ρ air C p u 3 * T s kgH (30) where g is the gravitational acceleration. At daytime L is usually negative and the boundary layer is unstable; −L is in this case the height at which the buoyant production of turbulence kinetic energy is equal to that produced by the shearing action of the wind; near the ground shearing action is dominant, and for neutrally stable conditions L is not finite.
The blending wind velocity is evaluated first using Equation (28) twice: the first time to evaluate friction velocity at the meteorological station from the measured wind and the second time to obtain blending velocity u z B from friction velocity at the station.
Then, since ψ(z) through L depends on u * and Equation (28) is not explicit, an iterative repeated evaluation procedure is activated to obtain the friction velocity map over the area. The procedure starts assuming neutral equilibrium conditions (ψ = 0), and evaluates u * with Equation (30), for given blending velocity and surface roughness z om . During iterations, adjusted values are computed also for H, L, ψ and r ah : where ψ h(z 1 ) and ψ h(z 2 ) are the stability corrections for heat transport at heights z 1 and z 2 . Empirical formulae for ψ(z) are given by Allen et al. [51], as Equations (41)- (45). The variable dT (K) is the vertical air temperature difference between the heights z 1 and z 2 . In SEBAL a linear relationship is assumed between the near-surface air temperature difference dT and the surface temperature T s : dT = a + bT s (32) where the coefficients a and b are determined empirically for a given satellite image through the so-called anchor pixels, the cold and hot pixels with assumed values of H. The cold pixel selected in well-irrigated fields is used to anchor dT = 0 which implies that LST equals air temperature canceling sensible heat flux convected to the air, H = 0. The hot pixel found in dry, bare agricultural soil, where evapotranspiration is assumed to be zero, is associated with a value of dT such that λE = 0 and H =R n − G. By comparing the conditions (Equation (32)) for the two anchor points, a and b can be easily obtained: and hence the following maps are obtained: • dT from T s and Equation (32); • H from dT and Equation (27); • finally, λE from H and the energy balance (Equation (18)). This is repeated until convergence is reached. In this study the cold and hot pixels were manually selected for each image in agricultural fields near the study area satisfying the indicated extreme conditions.
Once the instantaneous latent heat flux λE is calculated, it is used, by dividing for the latent heat of vaporization λ (J kg −1 ), to evaluate the instantaneous evapotranspiration E i (mm s −1 ) at the satellite overpass time.
In order to transform the instantaneous ET to daily average ET, that is the form that normally farmers use, for the specific day of satellite passage, ETr is evaluated at overpass time providing E r,i , as well as at several moments (every hour) from which, by integration over the day, E r,day is obtained. ETr estimates were performed by using the Penman-Monteith formula for alfalfa through the Ref-ET Reference Evapotranspiration Calculator [58].
The actual daily ET (E day ) is evaluated by the following proportion and is usually expressed in mm d −1 : assuming that the ratio between instantaneous ET and ETr is constant throughout the day.

Crop Water Stress Index (CWSI)
Crop water stress can be assessed from surface temperature in different ways (see [59,60], for instance). We used the Crop Water Stress Index (CWSI) introduced by Veysi et al. [61], that is defined so as to make use only of surface temperature data: where T V is the crop temperature, T hot is crop temperature under dry conditions and T cold under wet conditions, i.e., the two extreme temperatures used in the SEBAL model. By its definition CWSI is not sensitive to any linear transformation of temperatures, since the additive term is cancelled in Equation (36) differences, and the multiplicative factor in the quotient: it has for instance the same value when temperatures are measured in the Kelvin, Celsius or Fahrenheit scale. Since T V is not observable from satellite, T s is used, but care should be paid when canopy cover is low because this substitution could overestimate water stress.

Tile and Pixel Data
Fractional vegetation cover (F V ) and average temperatures (T LS ) over 30 m tiles, as the one presented in Figure 6, are obtained from the UAV survey by the averaging procedure described in Section 3.2. At the same time, almost contemporary surface temperature maps were obtained with methods presented in the same section from Landsat-8 images, distributed at 30 m resolution after a resampling of the 100 m TIRS bands by cubic convolution (L1T processing level). Table 1 presents the environmental conditions at UAV survey and Landsat-8 passage (i.e., air temperature Ta and vapor content in the atmosphere w), UAV derived temperatures over the 12 investigated tiles (T V , T S and T LS ), the temperatures of the corresponding 12 Landsat pixels (T SB , T RTE , and T SW ), the temperature correction term DT (which must be subtracted from UAV observed temperatures to compensate for the time lag between UAV and satellite acquisitions), and NDVI obtained from Landsat-8 level-1 data (top of atmosphere). For each survey (S01 to S09) the average value among the 12 investigated tiles/pixels is presented in the upper line and standard deviation of the population (not of averages) in the lower line. Table 2 presents average differences and average correlation coefficients among the sets of 12 corresponding values. The upper triangle shows correlation coefficients among pairs. The lower triangle shows differences (column-row).
The single tile/pixel values for UAV measured and satellite derived LSTs with the three methods are presented in Figure 7 (left), as well as the combined fractional vegetation cover derived from UAV thermal images and top of atmosphere NDVI retrieved from the satellite.
Effects of comparing data derived with different aggregation/disaggregation procedures are: (1) the variability of temperature among tiles/pixels is significantly lower in the Landsat pixels than in UAV tile averaged data; compare in Table 1 standard deviations in the second line of each survey; (2) contemporary tile temperatures derived from UAV and from Landsat-8 with any single method are very poorly correlated ( Table 2, first line), whereas correlation between pixel temperatures retrieved from satellite with the different methods is at least moderate.
temperatures being around 12 °C and difference between class averages around 9 °C. This allows to distinguish vegetation from soil and to evaluate FV in all tiles and surveys.
The good correlation between NDVI, derived from Landsat-8 OLI sensors at 30 m resolution, and FV evaluated from UAV survey over corresponding tiles is remarkable (see Figure 7, right) and supports the argument that the poor correlation of temperature is due to the difference in acquisition and distributed image pixel. The highest temperatures were measured by UAV; temperatures estimated by the RTE, SW, and SB methods are progressively lower, as shown in Table 2. The greatest differences occur in May and June, whereas at the end of the season difference are lower.
If we restrict comparison to UAV, RTE, and SW the differences can be however interpreted as due to systematic and random errors, separately of the order of 1 °C. Discrepancies between RTE and SW, based on the same information source, are normally lower than 1 °C.
Contemporary temperature has been measured also in the water bodies mentioned in Section 2.5. Since all these temperatures are measured in channels narrower than TIRS pixel, Landsat derived values are quantified for a wide upstream water body and compared with measured spilled water temperatures. Eleven reliable comparisons are obtained, and their statistics presented in Table  3. We interpret the absence of any correlation among UAV and satellite derived pixel temperature as the effect of aggregating actual variations of surface temperature, on one hand, and on the other as the effect of the convolution algorithm used to transform data from 100 to 30 m resolution; one responds to actual variation in prototype, the other is a mathematical artifact with poor connection with processes occurring at 30 m resolution.
The moderate correlation among SW temperatures ( Table 2, last column) and those derived from the other two methods is presumably due to the combination of small variability in pixel temperatures and different methods of corrections. SB and RTE methods, which use either no correction or a constant correction over the image related to environmental conditions of the moment, show unitary correlation and almost constant difference.
Temperature derived from Landsat TIRS data with SB, RTE, and SW methods (T SB , T RTE , and T SW ) are compared among each other and with LST derived from the UAV thermal camera (T LS ), after correction for non-synchronous observation with term DT in Table 1.
Since there is no correlation between contemporary temperatures derived at 30 m tile/pixel scale from UAV and Landsat data, only temperatures averaged over the 12 tiles/pixels are considered significant and compared.
Differences of several pairs remain significant (see Table 2, lower triangle). We can separate the comparison of different methods to derive temperatures from Landsat-8 TIRS measurements from the comparison with UAV measurement averages. SW method gives the highest temperatures, RTE the central ones, and SB the lowest. UAV average measurements still exceed the highest satellite derived temperatures by more than measurements uncertainty. We estimate that this is partly due to some residual variance in the average among surveys, but a difference is also due to emissivity assumed for the surface. In the analysis of UAV data, emissivity was 0.95 for soil and 0.98 for vegetation; in Equation (14) the two parameters representing the same variables are 0.971 and 0.987, plus an additive term and two variable bounds. In any case the difference is around 0.02, that according to Equation (4) for a brightness temperature around 305 K (30 • C) causes a surface temperature variation of 1.4 • C, systematically reducing the observed differences. Therefore, the difference between UAV and SW can be interpreted, as mainly due to the systematic difference of the assumed surface emissivities (1.4 • C of 1.6 • C).
The correlation coefficient of average values among tiles/pixels along the season (i.e., among columns of Table 1) is high, never descending below 0.85, or below 0.94 if only Landsat derived data are considered, and this confirms that an average over 12 tiles/pixel can be considered significant and used for further considerations.
Soil and vegetation temperatures result highly different, the average difference of peak temperatures being around 12 • C and difference between class averages around 9 • C. This allows to distinguish vegetation from soil and to evaluate F V in all tiles and surveys.
The good correlation between NDVI, derived from Landsat-8 OLI sensors at 30 m resolution, and F V evaluated from UAV survey over corresponding tiles is remarkable (see Figure 7, right) and supports the argument that the poor correlation of temperature is due to the difference in acquisition and distributed image pixel.
The highest temperatures were measured by UAV; temperatures estimated by the SW, RTE, and SB methods are progressively lower, as shown in Table 2. The greatest differences occur in May and June, whereas at the end of the season difference are lower.
If we restrict comparison to UAV, RTE, and SW the differences can be however interpreted as due to systematic and random errors, separately of the order of 1 • C. Discrepancies between RTE and SW, based on the same information source, are normally lower than 1 • C.
Contemporary temperature has been measured also in the water bodies mentioned in Section 2.5. Since all these temperatures are measured in channels narrower than TIRS pixel, Landsat derived values are quantified for a wide upstream water body and compared with measured spilled water temperatures. Eleven reliable comparisons are obtained, and their statistics presented in Table 3. Table 3. Relation among contemporary water bodies temperatures T W obtained on various dates by different methods. The upper triangle shows correlation coefficients among pairs. The lower triangle shows differences (column-row, • C). For these conditions only SB (uncorrected) estimates are lower than observed temperatures; RTE and SW corrected temperatures result around 2 • C higher than observed, but it must be remarked that, in summertime late morning, breeze and induced mixing are low, and some thermal stratification of the water body must be expected. Water is spilled from the water body almost uniformly along the vertical; average water body temperature 2 • C lower than surface temperature observed by the satellite might be real. Differences induced by the atmosphere correction methods are in this case weaker, ranking among methods remaining however the same: SB gives the lowest values, RTE the median, and SW the highest. Since satellite temperatures are evaluated over wide uniform surfaces all temperatures are mutually well correlated. Figure 8 shows the temperature distribution derived from SB method (the simplest one) at the experimental field and two dates. This image is for the end-user, who is presumably more interested in retrieving some information about his field than discussing precision of measurements. It represents also visually the resolution and the features that can be monitored from Landsat images. The S and W parts of the area are characterized by more permeable sandy soils that dry up in shorter time and become rapidly hot in dry conditions (3 July). After an irrigation event (4 August) the same area is uniformly cool. Based on these images and other considerations the farmer may decide to irrigate with variable intensity the portions of the field. This was not done in 2019. The farmer should, however, be informed that presented values should be incremented around 5 • C to represent the true surface temperature at noon. The S and W parts of the area are characterized by more permeable sandy soils that dry up in shorter time and become rapidly hot in dry conditions (3 July). After an irrigation event (4 August) the same area is uniformly cool. Based on these images and other considerations the farmer may decide to irrigate with variable intensity the portions of the field. This was not done in 2019. The farmer should, however, be informed that presented values should be incremented around 5 °C to represent the true surface temperature at noon. The correction methods represent radiative transfer of the atmosphere in different ways, but on a clear sky day the effect is almost constant over distances of dozens of kilometers, so that the correlation of the different satellite derived LST maps is very high and, for most uses, as for instance water stress identification, they result almost equivalent. The correction methods represent radiative transfer of the atmosphere in different ways, but on a clear sky day the effect is almost constant over distances of dozens of kilometers, so that the correlation of the different satellite derived LST maps is very high and, for most uses, as for instance water stress identification, they result almost equivalent. Figure 9 shows temperature maps on 3 July derived with RTE and SW methods over a 10.86 × 10.38 km wide area. The relation between temperatures in any two maps of this area, including the not presented SB derived map, is described in Table 4. Average differences do not deviate substantially from what is observed for dozens of tiles or water bodies; correlation is higher because results are derived by different analytical models based on the same image and no measurement uncertainty is introduced. In synthesis Landsat-8 derived temperatures uncorrected for atmosphere effects (SB method) underestimate UAV thermal camera measurements, whereas (see Tables 2-4), even if some differences remain for temperatures corrected according to methods RTE and SW, they do not differ between each other by more than 1 • C, RTE providing 0.6-1.0 • C lower average temperatures, and the major systematic recognized difference between our UAV measurements and SW derived estimates is the difference of surface emissivity assumed in the conversion radiance-to-temperature. Figure 9 shows that Landsat-8 derived maps are, however, highly correlated; they identify the same hot or cool areas at the patterns might be judged equal; the effects of these differences on most indicators used for detection of water stress conditions might be disregarded, as explained in the next section.

Surface Temperature Maps
Paying attention to small-scale effects, RTE images are more regular whereas some noise in the form of stripes is evident in SW images, that derives from the utilization of both TIRS bands; this evident feature balances the lower SW average discrepancy from UAV measurements.

Evapotranspiration and Water Stress Maps
Uncorrected SB temperatures are obviously the most simple to retrieve since the parameters obtained from Landsat-8 data are sufficient, whereas RTE and SW methods require more information (atmospheric transmittance, atmosphere upwelling and downwelling radiance; total atmospheric water content and surface emissivity difference between bands 10 and 11), so that, if the temperature pattern is the objective of the analysis, this method can be preferred. If ETa values are requested one should pay attention to the systematic underestimation of temperatures of this method and to the associated overestimation of ETa.
Daily ETa estimates obtained from SEBAL using SB temperatures are presented in the left column of Figure 10  Uncorrected SB temperatures are obviously the most simple to retrieve since the parameters obtained from Landsat-8 data are sufficient, whereas RTE and SW methods require more information (atmospheric transmittance, atmosphere upwelling and downwelling radiance; total atmospheric water content and surface emissivity difference between bands 10 and 11), so that, if the temperature pattern is the objective of the analysis, this method can be preferred. If ETa values are Figure 10. Evapotranspiration and water stress maps at and around the survey field on date 3 July, The first date, 3 July 2019, is short before an irrigation event and represents a dry condition of the experimental field when the estimated maize root depth was 0.76 m; the second date is short after an irrigation event and represents a wet condition; the last date represents conditions just before harvesting of maize at dough stage, when root depth exceeded 1 m, the field was not irrigated since 20 days and soil humidity was not far from the lower readily available water (RAW) limit. The right column of Figure 10 presents temperature maps at the same dates in the form of the Crop Water Stress Index that is defined so to make use only of surface temperature data. The hot and cold points are the same for both columns.
The round field on the left of the maize field is an irrigated soybean field, that remains transpiring and cool all along the growing season. The other cool areas are cane thickets, water bodies or woodland areas. Figure 11 shows the maps of ETa estimated with SEBAL and RTE (left) and SW (right) corrected temperatures. For comparison the correspondent map obtained from SB method is the top-left one in Figure 10. Patterns are quite similar, and correlation is very high, but a careful comparison will show that evapotranspiration amount is decreasing from SB down to SW, in the reverse order of temperatures. In fact, the higher is the temperature, the higher is the sensible heat flux from land surface to the atmosphere and the lower is latent heat flux necessary to balance net radiation minus ground heat flux. The daily ETa on 3 July 2019, averaged over the area presented in Figures 10 and 11, ranges from 5.68 mm/d (SW) to 6.09 mm/d (SB) and detailed comparison is presented in Table 5.
Remote Sens. 2020, 12, x FOR PEER REVIEW 25 of 32 requested one should pay attention to the systematic underestimation of temperatures of this method and to the associated overestimation of ETa. Daily ETa estimates obtained from SEBAL using SB temperatures are presented in the left column of Figure 10 for dates 3 July, 4 August, and 20 August 2019. The first date is short before an irrigation event and represents a dry condition of the experimental field when the estimated maize root depth was 0.76 m, the second date is short after an irrigation event and represents a wet condition; the last date represents conditions just before harvesting of maize at dough stage, when root depth exceeded 1 m, the field was not irrigated since 20 days and soil humidity was not far from the lower readily available water (RAW) limit. The right column presents temperature maps at the same dates in the form of the Crop Water Stress Index that is defined so to make use only of surface temperature data. The hot and cold points are the same for both columns.
The round field on the left of the maize field is an irrigated soybean field, that remains transpiring and cool all along the growing season. The other cool areas are cane thickets, water bodies or woodland areas. Figure 11. Evapotranspiration around the survey field on date 3 July obtained from SEBAL with RTE (left) and SW (right) correction methods of atmosphere effects on LST. Figure 11 shows the maps of ETa estimated with SEBAL and RTE (left) and SW (right) corrected temperatures. For comparison the correspondent map obtained from SB method is the top-left one in Figure 10. Patterns are quite similar, and correlation is very high, but a careful comparison will show that evapotranspiration amount is decreasing from SB down to SW, in the reverse order of temperatures. In fact, the higher is the temperature, the higher is the sensible heat flux from land surface to the atmosphere and the lower is latent heat flux necessary to balance net radiation minus ground heat flux. The daily ETa on 3 July 2019, averaged over the area presented in Figures 10 and  11, ranges from 5.68 mm/d (SW) to 6.09 mm/d (SB) and detailed comparison is presented in Table 5. Table 5. Comparison of actual evapotranspiration (ETa) estimated over the area presented in Figures  10 and 11 with SB, RTE, and SW methods based on 3 July 2019 Landsat scene. Correlation coefficient is presented above the diagonal and mean differences below it (mm/d, column minus row entries). CWSI is a linear transformation of temperature, it is therefore perfectly correlated with temperature, whichever is its distribution. The methods used to correct top of atmosphere temperature to account for atmosphere effects produce an almost linear transformation of the Figure 11. Evapotranspiration around the survey field on date 3 July obtained from SEBAL with RTE (left) and SW (right) correction methods of atmosphere effects on LST. CWSI is a linear transformation of temperature, it is therefore perfectly correlated with temperature, whichever is its distribution. The methods used to correct top of atmosphere temperature to account for atmosphere effects produce an almost linear transformation of the uncorrected temperatures. Therefore, CWSI or temperature pattern for different correction methods are practically identical and are not presented.

ESB
The relation between temperature and evapotranspiration is more complex, even if evident from Figure 10. Over a small and homogenous area, the correlation is strict, whereas, extending the area over which correlation is calculated or when the vegetation status/cover is variable, the correlation coefficient declines toward zero (Table 6). Table 6. Correlation coefficients between SB temperature and ETa over the experimental field and areas represented in Figures 9 and 10.

Date
Northern Parcel Southern Parcel Entire Figure 10 Entire Figure 9 03 Negative correlation (intrinsically dry areas are hot and wet areas cold) is very high only in the case of uniform conditions, as obtained with single crop under uniform growth conditions (i.e., in the northern parcel). The southern parcel suffered during the end of June heat wave, crop development was very irregular at the beginning of July and recovered partially in the two following months. Under different surface conditions, high temperature alone is not a good quantitative indicator of low evapotranspiration.

Discussion and Perspectives
Having verified that temperatures can be retrieved from satellite images with sufficient reliability, the scheduled goal of the experiment is achieved; moreover, derived actual evapotranspiration seems reliable; we can therefore forecast a non-academic implementation of the technique. Regarding the perspectives we would like to discuss briefly how the results can help the end-users (i.e., the farmers).
According to the principles of precision agriculture, a farmer may use any water stress indicator, as high CWSI value and/or low ETa, to plan irrigation according to crop need within the limit of soil storage capacity, and/or accounting for the excess water flow. Two problems are evident: • the long revisit interval of Landsat-8 (see Figure 3) and the frequent presence of cloud cover or atmospheric haze, which reduce dramatically the quality of satellite observations, limit the temporal resolution of useful images (around two useful passages per month) compared to irrigation frequency (one to three events per month); as a consequence the probability of not having the necessary information, when one needs it, is relevant so that the present technique, when used for irrigation management, may lead to disappointment; • the fact that soil storage capacity, groundwater conditions and groundwater dynamics should be known with a balanced reliability; soil in particular is not changing quickly so that information about its properties can be cumulated from a sequence of maps under several hydrologic conditions.
The relevance of the first problem can be significantly reduced when Sentinel-2 images can be used. In fact, as a consequence of shorter revisit interval (10 days), the presence of two satellites moving along the same orbit in phase opposition and of some sidelap between tracks, the time interval between passages is 2-3 days and the number of useful images is 7-8 per month during dry periods.
A useful index is the Normalized Difference Moisture Index (NDMI) [62], which is based on the effect of water content in canopy tissues on reflectance in bands near-infrared (NIR) and shortwave infrared (SWIR1) and is evaluated similarly to NDVI but with bands 5 and 6 for Landsat-8 and bands 8 and 11 for Sentinel-2. Figure 12 presents the distribution derived from almost contemporary passage of the two satellites on 3 July 2019. Landsat-8 map is obtained on demand from USGS. Sentinel-2 SWIR image with 10 m resolution was obtained after applying the super-resolution method described in [63] and integrated into the Sentinel Application Platform (SNAP) developed by the European Space Agency to increase the 20 m resolution of band 11 (SWIR1) to the highest resolution of 10 m.  NDMI interpretation requires some experience, since it is a significant moisture index only when vegetation is dominant, whereas it provides negative values for bare soil and water bodies so that when canopy cover is low the result is a weighted combination of soil and vegetation proper values. The effect of different resolution is evident. The correlation coefficient of the two images in Figure 12 is 0.8, because of the resolution difference and of a different evaluation of the index in the river Reno, that does not depend on the different central wavelength and bandwidth of Landsat band 5 and Sentinel band 8 (having 10 m resolution), because it remains unaltered if Sentinel band 8a is used (having 20 m resolution), which is practically coincident with Landsat band 5. The correlation coefficient between USGS NDMI and ESW is high (0.90 in the northern parcel) and similar to the correlation between CWSI and ESW.
The two water stress indexes, based on totally different physical bases and information, provide for crop water stress the same information (compare the upper right image in Figure 10 and the left image of Figure 12).
In the preliminary application the maps of water stress indicators (at least a dozen per year) will be easily used to recognize recurrent water stress patterns and under which meteorological conditions they appear, so that the farmer can analyze which are the most plausible causes of water stress, and which is the best solution for the specific case. In this preliminary phase the water stress spatial and temporal pattern is the most relevant information, whereas the exact intensity is related to the specific meteorological history and not easily utilized in the near future. The simplest solution will be probably the best one in this phase.
In order to upgrade the utilization of the thermal and water stress data they have to be related to properties that can be derived from OLI bands data, optical and infrared, that are obtained with greater spatial and temporal resolution. Sharpening of thermal data by fusion with short wave data has excellent perspectives. Improved observation capabilities can be provided also by recent and upcoming mini-and nano-satellite constellations, which are able to provide more frequent overpasses. In particular, the ConstellR technology seems promising for the monitoring of land surface temperature (http://www.constellr.space).
Higher spatial resolution will help the characterization of soils, and one map per week combined with pointwise monitoring of weather and soil conditions will be optimal for irrigation decisions. For operational irrigation scheduling purposes, satellite derived distributed estimates can be used in conjunction with point-wise monitoring of soil and meteorological conditions and with NDMI interpretation requires some experience, since it is a significant moisture index only when vegetation is dominant, whereas it provides negative values for bare soil and water bodies so that when canopy cover is low the result is a weighted combination of soil and vegetation proper values. The effect of different resolution is evident. The correlation coefficient of the two images in Figure 12 is 0.8, because of the resolution difference and of a different evaluation of the index in the river Reno, that does not depend on the different central wavelength and bandwidth of Landsat band 5 and Sentinel band 8 (having 10 m resolution), because it remains unaltered if Sentinel band 8a is used (having 20 m resolution), which is practically coincident with Landsat band 5. The correlation coefficient between USGS NDMI and E SW is high (0.90 in the northern parcel) and similar to the correlation between CWSI and E SW .
The two water stress indexes, based on totally different physical bases and information, provide for crop water stress the same information (compare the upper right image in Figure 10 and the left image of Figure 12).
In the preliminary application the maps of water stress indicators (at least a dozen per year) will be easily used to recognize recurrent water stress patterns and under which meteorological conditions they appear, so that the farmer can analyze which are the most plausible causes of water stress, and which is the best solution for the specific case. In this preliminary phase the water stress spatial and temporal pattern is the most relevant information, whereas the exact intensity is related to the specific meteorological history and not easily utilized in the near future. The simplest solution will be probably the best one in this phase.
In order to upgrade the utilization of the thermal and water stress data they have to be related to properties that can be derived from OLI bands data, optical and infrared, that are obtained with greater spatial and temporal resolution. Sharpening of thermal data by fusion with short wave data has excellent perspectives. Improved observation capabilities can be provided also by recent and upcoming mini-and nano-satellite constellations, which are able to provide more frequent overpasses. In particular, the ConstellR technology seems promising for the monitoring of land surface temperature (http://www.constellr.space).
Higher spatial resolution will help the characterization of soils, and one map per week combined with pointwise monitoring of weather and soil conditions will be optimal for irrigation decisions. For operational irrigation scheduling purposes, satellite derived distributed estimates can be used in conjunction with point-wise monitoring of soil and meteorological conditions and with agro-hydrological model, through data assimilation aiming to improve prediction of distributed soil water content.

Conclusions
The UAV thermal images were calibrated by comparison with temperature measured at TRPs and validated by comparison of a) temperature extrapolated at soil surface from measurement in the soil and b) temperature measured at near water bodies. Resulting temperatures are affected by a measured r.m.s. error around 1.5 • C; maps with very high resolution (5 cm) can be obtained.
Soil and vegetation temperatures result highly different, with the average difference being around 9 • C and the soil being warmer; this allows to distinguish vegetation from soil and to evaluate vegetation cover. The vegetation-air temperature difference is also high in the late mornings of late spring or early summer (up to 9 • C), making uncertain which is the average land surface temperature satellites do observe.
Temperatures derived from Landsat TIRS data, provided with resolution 30 m, according to SB, RTE, and SW methods were compared among each other and with LST derived from the thermal camera data averaged over tiles corresponding to Landsat pixels. Some differences are significant. In particular, neglecting atmosphere scattering of surface thermal radiation, SB method underestimates systematically the real LST.
There is no correlation between contemporary temperatures derived at 30 m tile/pixel scale from thermal camera and Landsat derived temperature, since the real resolution of TIRS images is 100 m. Therefore, only the comparison of temperatures averaged over several corresponding tiles and pixels is significant.
The good correlation between NDVI, derived from OLI sensors at 30 m resolution, and F V assessed over corresponding tiles is remarkable, and proves that averages over ground tiles corresponding to satellite acquisition pixels are significant and that the correspondence NDVI-F V is good and almost linear.
Over the summer season the highest temperatures were measured at land surface level by the UAV mounted thermal camera; temperatures estimated by the SW method were 1.7 • C below measurements on average; those estimated by the RTE method were 2.7 • C below measurements, and those estimated by SB method (without atmosphere correction) are 5-6 • C below measurements.
Different surface emissivity values assumed in UAV and Landsat-8 radiance to temperature transformation explain an average difference of 1.4 • C (i.e., most part of the observed difference between UAV and SW derived temperatures).
SW temperatures show some noise in the form of stripes, which is not accounted in averages. RTE estimate may be considered within measurement error range and finally the two correction methods produce reliable results. Differences are essentially due to different corrections of the atmosphere effects that are not significantly variable at scale of 10 km. Correlation among satellite derived temperature estimates over such areas are evaluated and are variable between 0.984 and 0.994, showing that almost exact linear relations hold among them; the relations are not far from a pure translation.
The pattern of CWSI is practically independent from the correction method. LST pattern as well as the related evapotranspiration are therefore affected by minor errors and are reliable and useful for irrigation management. In order to obtain quantitatively correct temperatures and evapotranspiration, methods correcting for atmosphere effects, as RTE or SW, should be used.
Similar information regarding crop water stress can be obtained from NDMI, based on bands that are provided by both Landsat-8 and Sentinel-2, which can be retrieved usually every 4-5 days, and can be used either directly or combined with thermal indicators to suggest irrigation practice.
Combination and assimilation of point-wise monitoring and agro-hydrological modeling with extensive satellite monitoring provide excellent perspectives for the near future.