Mapping Very-High-Resolution Evapotranspiration from Unmanned Aerial Vehicle (UAV) Imagery

: There is a growing concern about water scarcity and the associated decline in Australia’s agricultural production. Efﬁcient water use as a natural resource requires more precise and adequate monitoring of crop water use and irrigation scheduling. Therefore, accurate estimations of evapotranspiration (ET) at proper spatial–temporal scales are critical to understand the crop water demand and uptake and to enable optimal irrigation scheduling. Remote sensing (RS)-based ET estimation has been adopted as a method for large-scale applications when the detailed spatial representation of ET is required. This research aimed to estimate instantaneous ET using very-high-resolution (VHR) multispectral and thermal imagery (GSD < 8 cm) collected using a single ﬂight of a UAV over a high-density peach orchard with a discontinuous canopy. The energy balance component estimation was based on the high-resolution mapping of evapotranspiration (HRMET) model. A tree-by-tree ET map was produced using the canopy surface temperature and the leaf area index (LAI) resampled at the corresponding scale via a systematic feature segmentation method based on pure canopy extraction. Results showed a strong linear relationship between the estimated ET and the leaf transpiration ( n = 42) measured using a gas exchange sensor, with a coefﬁcient of determination (R 2 ) of 0.89. Daily ET (5.5 mm d − 1 ) derived from the instantaneous ET map was comparable with daily crop ET (6.4 mm d − 1 ) determined by the meteorological approach over the study site. The proposed approach has important implications for mapping tree-by-tree ET over horticultural ﬁelds using VHR imagery.


Introduction
Water scarcity has long been an important issue in Australia, and has resulted in declined production in irrigated agriculture at a national scale [1]. Efficient water use in agriculture requires improved irrigation management informed by accurate monitoring of crop water requirements. In the case of high-value horticultural crops, e.g., pear, peach, grape, and olive, various studies have investigated the effect of deficit irrigation, including regulated deficit irrigation (RDI) [2,3], to increase fruit quality while maintaining yield [4,5]. The accurate measurement of evapotranspiration (ET) plays a critical role in determining optimal irrigation water volume and understanding crop water uptake and loss. In particular, crop ET estimations inform farmers of water consumption and plant available water in the soil, which are key inputs for irrigation management and scheduling.
Various methods have been used to derive field-based ET; weather station based approaches, surface energy balance methods (SEBM), and soil water balance have been investigated and utilized in point measurements for several decades [6]. However, fieldbased ET has a limited utility when expanding the ET values over a large area. In particular, it is not feasible to estimate intra-field ET variations within a crop area [6]. With the consideration of spatial variations within land features (a mixture of vegetated and soil area), remote sensing (RS)-based ET has been adopted as a reliable and feasible method for estimating ET at regional scales [7]. Specifically, remotely sensed ET maps could reproduce spatially varying ET over heterogeneous crop fields at the resolution of the sensor. Satellite images have been used as a major input for RS-based ET estimation for several decades [8]. In general, four categories were proposed as RS-based ET methods [6,9]: (1) the direct empirical method using the direct relationship between ET and the combined RS and meteorological data; (2) the residual method of surface energy balance model; (3) the vegetation index method where a potential or reference ET is calculated from ground measurements and RS data is used to estimate the crop specific factor (crop coefficient) and; (4) deterministic methods based on the soil-vegetation atmospheric transfer (SVAT) model. Among these approaches, the residual method has been widely applied to estimate ET (or latent heat flux), where ET is obtained as a residual in SEBM and other energy balance components (net radiation, sensible heat flux and soil heat flux) are estimated based on the combined empirical and physical relationships [10]. The residual method in SEBM employs remote sensing data as surface temperature input and canopy structural characteristics derived from vegetation indices (e.g., normalized difference vegetation index (NDVI) or leaf area index (LAI)).
Some of the commonly used techniques to estimate ET are the surface energy balance algorithm (SEBAL) [11], the simple surface energy balance index (S-SEBI) [12], the mapping evapotranspiration with internalized calibration (METRIC) [13], and the two-source energy balance model (TSEB) [14,15]. SEBAL, SEBI, and METRIC require the cold and hot reference pixels within the image to estimate ET. The two extreme reference temperatures are allocated to set the limit of the evaporative fraction (0 and 1) where the latent heat flux (λET) is set to 0 on a dry surface and sensible heat energy (H) is set to 0 on a wet surface, respectively [11]. On the other hand, the TSEB does not require the extreme references within the image and partitions the sensible and latent energy fluxes separately from the canopy and soil in one pixel (or unit area). TSEB has been proven to be a reliable method to estimate E and T for uniform and row-crop fields.
The RS-based methods were designed to estimate ET over a relatively large area with primary input images from satellites. Satellites can provide the visible, near-infrared (VISNIR), and thermal (TIR) bands with a spatial resolution of 30 m-1.0 km, e.g., Landsat 7 with VISNIR 30 m and TIR 60 m, Landsat 5 with 30 m and 120 m, and MODIS with 0.5 km and 1 km, respectively. However, when an ET estimate is required at the resolution of small fields or trees for precision agriculture, the resolution of typical satellite data is too coarse to provide the fine information required for plant water use. In addition, the temporal resolution is rather low to capture a target area in a projected time on demand since a satellite's revisit cycle is usually several days [8,16,17]. On the other hand, highresolution TIR and VISNIR imagery acquired by an unmanned aerial vehicle (UAV) can provide surface information at a resolution higher than 10 cm depending on the flight altitude and camera sensor specifications [18][19][20][21]. The high-resolution imagery from a UAV enables capturing tree-level variabilities of crop water use and status, such as ET and water stress symptoms. It also facilitates more precise delineation of canopy pixels from the soil background. In the case of heterogeneous fields, e.g., fruit tree crops, if low-resolution imagery is used, it is difficult to show intra-field characteristics separately from tree rows and inter-rows as the imagery includes many mixed pixels of trees and soil [22]. The TSEB has been widely applied to partition the energy fluxes from the tree canopy and soil in the mixture pixels on this subject. As one of the ET methods for high-resolution imagery, Zipper and Loheide Ii [16] introduced a combined model of one-source and two-source schemes using a pixel-based ET estimate, called high-resolution mapping of evapotranspiration (HRMET). HRMET is designed to partition the available energy (net radiation) to the canopy and soil in a pixel and then calculates the sensible and latent heat flux iteratively. HRMET relies on simplified inputs that can be practically obtained from remotely sensed data and basic meteorological data. The other advantage of the model is that it does not require wet and dry reference features to calculate turbulent fluxes in the imagery.
In this study, tree-level energy balance components were estimated based on HRMET using UAV-borne multispectral and thermal data. The study area is a small-size peach orchard showing heterogeneous field characteristics due to tree crops with tree rows and inter-rows. Therefore, typical coarse resolution from satellite imagery has limitations in the intra-field ET characterization. Thus, this study included very-high-resolution imagery to represent the spatial variability of ET more accurately within the field by dealing with heterogeneous land cover, discontinuous canopy cover, and mixed features of vegetation and soil. The main objective is to assess the sub-field-scale variability of ET over a dripirrigated peach orchard.

Study Site
The experiment was carried out on 19 January 2017 at peach and nectarine orchards  The study site has two tree crops: a peach (Prunus persica (L.) Batsch cv. August Flame); and a nectarine (Prunus persica (L.) Batsch cv. Autumn Bright). The nectarine and peach trees were planted in 2013, presenting a late fruit maturity phenological stage at full leaf-up at the moment of measurments. Each cultivar is trained to two canopy structures: a vertical leader (VL) with a canopy size of 2.0 m high × 0.8 m wide, and a Tatura trellis (TT) with a canopy size of 1.9 m high × 1.0 m wide trellis arms (Figure 1b). TT is in a Y shape with two arm branches since TT's architecture is designed to maximize the leaf surface exposed to capture the sun to enhance productivity [23,24]. The trees were planted in a north-south direction in fine, sandy, loamy, Shepparton soil. The tree spacing is 1.0 m and the inter-row distance is 4.5 m. The orchard was drip-irrigated daily, and the irrigation amount was calculated by a weather-based ETc (crop evapotranspiration) model [25].

Net Radiation and Micrometeorological Data
A tower was installed to measure both net radiation and micrometeorological data near the center of the study area. A net radiometer with four components (NR01, Hukseflux Thermal Sensors B.V., Delft, The Netherlands) measured net radiation, incoming/outgoing short wave radiation, and longwave radiation 2.6 m above ground level (AGL) and 0.6 m above the trees. Temperature and humidity (HMP155, Vaisala Corporation, Helsinki, Finland) sensors were installed above the tree. A wind monitor (Wind Monitor-AQ, R.M. Young Company, Traverse City, MI, USA) was placed at the top of the tower at 4 m AGL and measured wind speed and direction.

Crop Physiology and Solar Radiation Interception Data
The measurement of crop physiological data was carried out at four adjacent trees in the same three-row where the tower was installed. Leaf temperature was measured in both sunlit and shaded leaves using a handheld infrared thermometer (TN410LCE, ZyTemp, Radiant Innovation Inc., HsinChu, Taiwan). Gas exchange measurements were performed using a photosynthesis system (LI-6400, LI-COR Inc., Lincoln, NE, USA) over three east-sided and three west-sided mature and fully expanded leaves to obtain stomatal conductance (g s , mol m −2 s −1 ) and transpiration rate (E, mmol m −2 s −1 ) per tree was measured (n = 6 leaves × 7 trees = 42 observations ). The measurements were averaged to a single transpiration rate to represent a tree's status. Photosynthetically active radiation (PAR) was measured above and under the canopy for light interception using a ceptometer (Sunfleck, model SF-80, Decagon Devices Inc., Washington, DC, USA). All measurements were carried out near solar noon.

UAV Field Campaign
The UAV campaign was conducted on the day of clear sky with moderate wind (0.6 m s −1 ), air temperature (30.6 • C), and relative humidity (26.7%). A thermal infrared (TIR) camera (A65, FLIR Systems, Inc., Wilsonville, OR, USA) and multispectral (MS) camera (RedEdge, MicaSense, Seattle, WA, USA) were integrated with a GPS in an onboard CPU for geo-tagging. All the instruments were mounted to a UAV platform (S1000, DJI, Shenzhen, China). The UAV used was an octocopter aircraft, and the maximum capacity of the payload is 11 kg.
All aerial images from TIR and MS sensors were taken by a single flight within the short time window (<15 min) to capture the homogenous features at midday over the site. TIR images were captured in the spectrum of 7.5-13 µm, a spatial resolution of 640 × 512 pixels, a focal length of 25 mm, and a FOV of 25 • (H) × 20 • (V). MS images were acquired concurrently with TIR sensing. The MS camera consisted of five discrete bands in the blue (475 nm), green (560 nm), red (668 nm), red edge (717 nm), and near-infrared (NIR) (840 nm) spectra. Spatial resolution was 1280 × 960 pixels with a focal length of 5.5 mm.
Both cameras were mounted to an active gimbal to enable the images to be acquired at nadir view from the UAV. The UAV sensing was conducted at solar noon to capture the period of high ET and minimize the shadows cast by the trees. The UAV flew at an altitude of 90 m AGL to capture images with over 80% forward and 40% side overlap by an autonomous flight plan. In total, 204 snapshots of TIR and 199 multi-band images from the MS camera were captured. The footprint of the TIR image was 39 m × 31 m with a ground sample distance (GSD) of 6 cm. The MS image had a footprint of 108 m × 81 m and a GSD of 8 cm.
During the UAV sensing, two types of ground targets were deployed at the site: (1) radiometric calibration target for TIR and MS images, and (2) ground control point (GCP) ground artificial feature (GAF) for image processing. The water body and rubber plates were used to calibrate TIR images as cold and hot features, respectively. Tarps (3 m × 3 m) in three different reflectance values (6%, 12%, and 33%) were placed for the calibration of MS images ( Figure 2). In total, 24 GCPs were distributed over the field and surveyed by a DGPS with 3 cm positional accuracy. The specific GCP and GAF targets were designed and used for the bundle block adjustment. These targets were made of aluminum sheets considering the spatial resolution and spectral characteristics of both TIR and MS images.

Remote Sensing Data Processing and ET Modelling
The data processing and modeling consisted of four steps ( Figure 3): (1) RS data processing to produce surface temperature and canopy characteristics variables of ET modeling; (2) tree-by-tree modeling to identify individual plants and allocate the representative spatial variable to each tree; (3) calculation of energy balance components in the ET model and mapping instantaneous ET rate and; (4) evaluation of the estimated ET with ground measurements.

Aerial Imagery Processing
When processing the very high-resolution imagery, a post-calibration process was carried out to produce reliable inputs, which are key factors for the accurate ET estimate. TIR images were recorded in a raw format, consisting of the signal-based values from the surface temperatures. A one-point calibration method [26] was applied to convert the raw digital numbers to temperature values using customized code written in Matlab (Mathworks Inc., Matick, MA, USA). When correlating the raw values to the temperature, calibration targets made of rubber sheets and water bodies were used to refer to the actual temperatures as hot and cold features. The target temperatures were measured with TIR snapshots by a handheld thermal imaging camera (T640, FLIR Systems, Inc., Wilsonville, OR, USA) and a handheld thermometer concurrently at the time of the UAV flight. In the case of MS imagery, radiometric calibration was performed to retrieve the reflectance values in each band. The reflectance measured from three different uniform tarps was correlated to the digital number (DN) values from the MS imagery, and the linear regression relationship was applied to convert all DN values to the reflectance values. All consecutive TIR and MS images were aligned with initial camera positions from onboard GPS in a photogrammetric software (PhotoScan, Agisoft LLC, St. Petersburg, Russia). Then, image processing of point cloud and build mesh/texture/DEM/orthomosaic was performed with the highest accuracy setup. RMSEs of georeferenced MS and TIR were horizontally 0.084 m and 0.137 m and vertically 0.229 m and 0.319 m, respectively.

Leaf Area Index (LAI) Estimation
The HRMET model used in this research requires LAI to determine canopy transmittance. LAI can be parameterized as a function of multispectral vegetation indices such as NDVI, simple ratio, and reduced simple ratio [27,28]. For example, NDVI is highly correlated with LAI, and the relationship can be adopted in various forms of a regression model (e.g., general exponential, simple linear, or higher curve model) [28,29]. In this research, photosynthetically active radiation (PAR)-based LAI is calculated and correlated to NDVI-based LAI, to obtain a more realistic structure of the canopy per tree as an input of the ET model. The exponential regression method [27] was employed to retrieve LAI from remote sensed NDVI.
The reference LAI values were determined using PAR measurements above and under the canopy used to calculate light interception by the canopy in the four sample trees. The measured PAR of each tree consists of transmitted and scattered radiation through the canopy and within the canopy, respectively. A radiative model of transmitting and scattering was introduced by Norman and Jarvis [30] and simplified [31] as follows: where τ: canopy transmittance, A: canopy absorption, f b : fraction of direct radiation, L: LAI, K: coefficient of canopy extinction, and θ: sun zenith angle. The canopy transmittance τ is the ratio of the transmitted PAR to the incident PAR, and f b is the fraction of radiation from the solar beam to the total incident radiation, which is also known as beam fraction. K is the canopy extinction coefficient and can be simplified with the solar zenith angle (θ). The PAR-based LAI in the reference tree was retrieved by inverting the radiative model (Equation (1)). Then, the LAI map over the study site was estimated by the exponential regression model with the NDVI map generated from MS imagery obtained from UAV sensing.

Tree Segmentation
The surface temperature and LAI maps described in the previous sections were generated into the co-registered raster imagery within an approximately 1 pixel error with geo-referenced coordinates. The tree-level analysis was conducted to systematically produce individual plant water use, automated by the informative field descriptions (plant spacing, row spacing, etc.). For the tree-by-tree analysis, the feature segmentation technique was used in this research. First of all, the non-canopy background pixels were excluded in both temperature and LAI imagery using a histogram method. A set of threshold values for surface temperature and LAI was determined to separate the canopy and non-canopy distributions in the histogram. Although we used very-high-resolution imagery over the experimental site, image pixels along the boundaries between the canopy and background soil can present biased surface temperature and NDVI due to the mixed soil-canopy features for the edge pixels and possible leakage of heat from adjacent pixels. To improve the canopy temperature estimation accuracy, the edge pixels along the boundary between canopy and background were removed using the edge detection method. The details of the edge detection method used in the research were presented in the previous study [32]. Second, pixels belonging to individual trees were classified and grouped systematically by the informative field descriptions, i.e., the plant spacing and row spacing distances. The grouped pixels were aggregated and averaged as a representative value of each tree using customized code written in Matlab (Mathworks Inc., Matick, MA, USA).

RSEB Algorithm in ET Modeling
The high-resolution mapping of evapotranspiration (HRMET) model was implemented to estimate energy balance components, and subsequently ET in this research. HRMET is a pixel-based surface energy balance model that was developed to incorporate high-resolution RS data. The model's most prominent advantage is that it does not depend on wet and dry reference features to calculate turbulent fluxes in the imagery. Most of the model's input variables can be obtained from RS data, and the model does not require a heavy process and is executed with minimized inputs as a feasible tool to adapt an applicable irrigation scheduling. The algorithm of the model follows the residual method of the general surface energy balance model as in Equation (3): where λET: latent heat flux (W m −2 ); R n : net radiation at the surface (W m −2 ); H: sensible heat flux to the air (W m −2 ) and; G: soil heat flux (W m −2 ). The instantaneous net radiation was calculated to separate canopy and soil components based on the two-source model [14] as follows: where S in , S out , L in , L out is incoming and outgoing shortwave and incoming and outgoing longwave (W m −2 ), respectively. The fractional cover f c is derived from LAI. The subscripts c and s stand for canopy and soil. α and ε represent the surface albedo and emissivity. σ is the Stefan-Boltzmann constant, and T is the surface temperature in Kelvin. The incoming and outgoing shortwave radiation was calculated separately for the canopy and soil using the radiation extinction model in Equations (4)-(7) and the Stefan-Boltzmann law equation (Equation (8)). Then, soil heat flux was approximated to be 35% of the net radiation reaching the soil, as follows [14]: The sensible heat flux, H, was computed as follows: where ρ a : the molar density of air (mol m −3 ), c p : the specific heat of air (J mol −1 C −1 ), г Ha : the aerodynamic resistance to heat transport (s m −1 ), and г ex : the excess resistance to heat transport (s m −1 ). The model employs an iterative convergence method to calculate H, since both H and г Ha are codependent variables. The level of H affects the aerodynamic stability and the aerodynamic resistance to heat transport [11]. The iterative method of H is described in a flow chart in Figure 4. This method does not require the wet and hot extreme pixels due to the turbulent heat fluxes' iterative calculation.
When the convergence criterion is met, which is a less than 0.1% change between the consecutive iterations, the latent heat flux of each pixel was computed using the updated sensible heat flux in Equation (3).

Site-Specific Crop ET
The crop ET (ET c ) value of the study field was calculated by using reference crop ET (ET o ) and the basal crop coefficient (K cb ) adjusted for tree size as follows: where K c is the crop coefficient and can be divided into two coefficients: the soil evaporation coefficient (K e ) and the basal crop coefficient (K cb ) refer to the contribution of soil evaporation and crop transpiration in modeling crop ET, respectively [25]. K cb has been considered to improve the estimation of daily crop water use in irrigated row crops. Thereby, ET c is expressed as follows: Previous research has determined a K cb of peach whereby a site-specific crop coefficient adjusted for tree size (1.52 × midday light interception) permits the estimation of orchard ET after Goodwin, et al. [33].
As ET c indicates the daily evapotranspiration, the estimated hourly ET from the UAV sensing was upscaled to the daily value using the extrapolation method [34] as follows: where LE i and R n,i are instantaneous latent heat flux and net radiation and obtained from the estimated energy components. R n,d is the mean net radiation daily (24 h) and measured by the net radiometer, c f is a unit conversion factor of time, and λ v and ρ w are the latent heat of vaporization and water density, respectively.

Canopy Temperatures and NDVI
The study area consisted of two canopy systems for both nectarine and peach trees: vertical leader (VL) and Tatura trellis (TT) as shown in Figure 1b. At the time of UAV flight over the study area, the canopy temperature, captured from TIR sensing, showed spatial variations regarding each canopy system and cultivar. The first three rows were planted with nectarine of VL canopy and had lower temperatures (29.8 • C) compared to the average temperature of the other tree rows (31.7 • C), showing a temperature difference of approximately 2 • C. Figure 5a shows the intra-field spatial variability of surface temperatures without the soil background in an enhanced contrast colormap. The distributed canopy temperatures over the field were visualized; the temperature of nectarine with VL was cooler than that of the others. The entire study field was under a full-irrigation regime. The coefficient of variation (CV) of canopy temperature was 5.3% over the orchard. Thus, this result indicated that the canopy temperatures (Tc) varied from the different canopy structures and would result in the variability of ET rate. The spatially distributed NDVI was generated using VHR red and NIR imagery from the UAV sensing, shown in Figure 5b. The figure presents canopy NDVI, where soil NDVI was classified and excluded using the histogram-based method for the tree segmentation process in the next step. Canopy NDVI also showed spatial variability with 4.5% of CV, which was slightly lower than the CV of the canopy temperature during the UAV campaign. The canopy NDVI in the reference sample trees was compared with the PAR-based LAI and fitted by an exponential relationship [35], resulting in a coefficient of determination (R 2 ) of 0.9 as shown in Figure 6. We assumed NDVI = 0 when LAI was reduced to zero. LAI values of the canopies over the whole field were retrieved using the exponential function of NDVI vs. LAI as shown in Figure 5. Since the UAV data acquisition was carried out to capture instantaneous features over the field at the maturity stage, NDVI values from the reference trees varied only within a narrow range of around 0.75. Ideally, a broad distribution of NDVI vs. LAI is required to accurately model the relationship. However, in this research, it was assumed that the retrieved LAI over the field follows the reference PAR-LAI range with an ignorable error, since NDVI values from all trees were in the narrow distribution and were similar to NDVI values from the reference trees [27]. Figure 7 shows the results of the tree-by-tree segmentation both in surface temperature and LAI map. The result of tree-by-tree segmentation represents the possible canopy areas of individual trees. The orchards consisted of 12 rows, and 120 trees were planted in each row. According to the field descriptions, tree and inter-row spacing are 1.0 m and 4.5 m, respectively. Using the field descriptions and geo-referenced positioning of each row and tree, each tree section was segmented. Then, the mean value of surface temperature and LAI were assigned to each section as a representative tree value. The results provided better analytic maps to feature the individual plants and interpret the variability between plants, enabling the analysis of tree-by-tree water losses.

Energy Balance Components
The energy balance components, including latent heat flux, were estimated using the HRMET model, and the hourly ET was calculated over the study field, as shown in Figure 8. The original HRMET model has been tested in a uniform surface area with a spatial resolution of 1-2.5 m of thermal and 30 m of multispectral inputs under the assumption of uniform meteorological conditions over small-sized fields (9 ha and 12 ha). Although this study site presents different site characteristics (e.g., discontinuous canopy cover, tree crops, inter-row soil), HRMET can be applied to examine high-resolution ET since only vegetated pixels were aggregated to individual trees and calculated in the energy balance model. This was possible as very-high-resolution (<8 cm) RS data and a small-sized field (<1 ha) were employed in this study, so precise canopy and soil pixels were classified on account of sufficient spatial resolution data and the removal technique of mixed pixels. The estimated ET in the reference trees was approximately 0.62 mm h −1 . A relatively higher ET was observed in the first three rows (nectarine with VL) and in the northern part of the tree rows, presenting a rate of 0.76 mm h −1 at the most. Since the study field was a small orchard and the UAV data acquisition time was less than 15 min, the meteorological variables such as incoming shortwave radiation, wind speed, and vapor pressure were regarded as consistent across the field. Thus, the different ET rates along the trees were determined, which were mainly derived by the differences of canopy temperature, crop type, and LAI (hence, fractional vegetation cover), as similar patterns were confirmed in the surface temperature map and NDVI (or LAI) map (Figure 7).
The estimated ET was compared with the leaf transpiration rate (mmol m −2 s −1 ) measured by the chamber (6 cm 2 ) of gas exchange analyzer (LI-6400) for the reference sample trees, and results with standard deviations are shown in Figure 9.
The estimated ET showed a strong relationship (R 2 = 0.89) with the average value of leaf transpiration. A coefficient of determination (R 2 ) would not be a robust indicator to show a significance of correlation with a small dataset. Hence, the magnitude of the relationship (Figure 9) could be over-estimated, as the dataset consists of a small number of measurements. However, the result indicates that the estimated canopy ET was correlated with the ground measurement of leaf transpiration. The leaf-level transpiration measurements may not be the most ideal method as ground truth to compare with treelevel ET in general cases. Since the main factor of the difference in tree-transpiration and tree-ET is related to LAI [35,36] and the experimental orchards were of the same age and maintained rather uniformly in canopy structure and size, it is acceptable to show the correlation of leaf transpiration and tree ET. Although the leaf-level transpiration is not comparable with the UAV-borne ET in absolute quantity, the method employed suits the primary focus of demonstrating the utility of UAV-borne TIR imagery to map spatially varying tree-by-tree ET in this work. Although the direct measurement of ET such as eddy covariance was not available as the validation data in this research, the crop ET (ET c ) value of the study field was compared with the estimated ET to evaluate the method. Since ET c is a single value at a field-scale, it has a limitation to present the field-scale comparison with UAV-ET without a sufficient number of repetitions. Under the circumstance of a single-day UAV campaign, the main focus of the comparison, however, was to check that the UAV-born ET is within a reasonable range derived from ET c .
As a result, the extrapolated daily (UAV) ET was obtained as 5.50 (mm d −1 ), whereas the crop ET was calculated as 6.35 (mm d −1 ). The difference between the two daily ET values was 0.85 (mm d −1 ). Considering that the extrapolated ET may include estimation errors (e.g., ranging from 0.25 to 1.17 mm d −1 ) depending on the extrapolation methods and crop types [34], the estimated daily ET value is within the range of expected ET, showing a comparable value for crop ET. Although it was challenging to evaluate the results thoroughly due to the absence of directly measured ET or multi-seasonal UAV data, the estimated results were compared with leaf transpiration and the daily crop ET to show that the method has the potential to estimate tree-by-tree ET with intra-field variability and to accommodate VHR imagery. This research method can be supported by further experiments and explored in different sites/cultivars/phenological stages.

Conclusions
This study examined ET's estimation using very-high-resolution (VHR) multispectral and thermal imagery derived from UAV sensing. The energy balance components were estimated based on the HRMET surface energy balance model. Tree-by-tree analytic maps were produced by the systematic feature segmentation method based on pure canopy extraction and statistical analysis of the distribution of surface temperatures and LAI. A strong linear relationship between the estimated ET and leaf transpiration was obtained in this work. The estimated ET presented a close value to the crop ET over the study site. As remotely sensed NDVI has a limited capability to represent the heterogenous canopy structure of trees, particularly for the vertical profile of canopies, NDVI-based LAI in the proposed method likely poses the same limitation when it is applied to the fields of trees with complex and heterogeneous vertical structures, such as a forest. Thus, the proposed method would be more suited to orchards where planted trees are of similar genetics and age and managed by an agronomically uniform method with rather homogenous canopy structure and tree size. The proposed approach can potentially provide a practical method of assessing the intra-field variability of tree-by-tree ET at a sub-field scale for precision irrigation scheduling. Further research on VHR ET estimation comparing RSEB models and crop types is required to develop robust water management information and strategies.