High Resolution Geospatial Evapotranspiration Mapping of Irrigated Field Crops Using Multispectral and Thermal Infrared Imagery with METRIC Energy Balance Model

: Geospatial crop water use mapping is critical for field-scale site-specific irrigation management. Landsat 7 / 8 satellite imagery with a widely adopted METRIC (Mapping Evapotranspiration at high Resolution with Internalized Calibration) energy balance model (LM approach) estimates accurate evapotranspiration (ET) but limits ﬁeld-scale spatiotemporal (30 m pixel − 1 , ~16 days) mapping. A study was therefore conducted to map actual ET of commercially grown irrigated-ﬁeld crops (spearmint, potato, and alfalfa) at very high-resolution (7 cm pixel − 1 ). Six small unmanned aerial system (UAS)-based multispectral and thermal infrared imagery campaigns were conducted (two for each crop) at the same time as the Landsat 7 / 8 overpass. Three variants of METRIC model were used to process the UAS imagery; UAS-METRIC-1, -2, and -3 (UASM-1, -2, and -3) and outputs were compared with the standard LM approach. ET root mean square di ﬀ erences (RMSD) between LM-UASM-1, LM-UASM-2, and LM-UASM-3 were in the ranges of 0.2–2.9, 0.5–0.9, and 0.5–2.7 mm day − 1 , respectively. Internal calibrations and sensible heat ﬂuxes majorly resulted in such di ﬀ erences. UASM-2 had the highest similarity with the LM approach (RMSD: 0.5–0.9, ET dep,abs (daily ET departures): 2–14%, r (Pearson correlation coe ﬃ cient) = 0.91). Strong ET correlations between UASM and LM approaches (0.7–0.8, 0.7–0.8, and 0.8–0.9 for spearmint, potato, and alfalfa crops) suggest equal suitability of UASM approaches as LM to map ET for a range of similar crops. UASM approaches (Coe ﬃ cient of variation, CV: 6.7–24.3%) however outperformed the LM approach (CV: 2.1–11.2%) in mapping spatial ET variations due to large number of pixels. On-demand UAS imagery may thus help in deriving high resolution site-speciﬁc ET maps, for growers to aid in timely crop water management.


Introduction
Satellite-based remote sensing (RS) has been extensively used with energy balance models for regional scale evapotranspiration (ET) mapping [1][2][3][4][5][6][7][8]. One such widely adopted model is Mapping ET at high-resolution with Internalized Calibration (METRIC) [9,10]. The METRIC model uses Landsat 5/7/8, or Terra and Aqua satellite based multispectral and thermal infrared imagery data to compute ET as an energy balance residue. METRIC is advantageous for its independence from surface
UAS flights were configured using a ground control software (MissionPlanner, version 1.3.49, Ardupilot, USA) for altitude of 100 m above ground level (AGL) to acquire multispectral and thermal images at ground sampling resolutions of 7 and 13 cm pixel −1 , respectively. The software aided in configuring the multispectral sensor to acquire images at 85% front and 75% side overlaps. The thermal imaging sensor was configured independently to acquire frames at 3 Hz to ensure image overlaps of 90-95%. Imagery data was stored on-board in respective memory cards for the two sensors. A calibrated reflectance panel (CRP, Micasense Inc., Seattle, WA, USA) was imaged before and after Drones 2020, 4, 52 4 of 18 each flight and these reference images were used for radiometric calibration of the multispectral imagery ( Figure 2).
Total six UAS flights were conducted in summers of 2018 (two missions/site × 3 sites) during the dates and approximate times of the Landsat 7/8 overpass. UAS imagery for spearmint crop was collected 10 days before the first harvest (dataset-1) and 37 days before the second harvest (dataset-2). Imagery for the potato crop was acquired 72 days (dataset-3) and 48 days before harvest (DBH, dataset-4). Data for the alfalfa crop was collected 2 days before the first harvest (dataset-5) and 7 days before the second harvest (dataset-6).
UAS flights were configured using a ground control software (MissionPlanner, version 1.3.49, Ardupilot, USA) for altitude of 100 m above ground level (AGL) to acquire multispectral and thermal images at ground sampling resolutions of 7 and 13 cm pixel −1 , respectively. The software aided in configuring the multispectral sensor to acquire images at 85% front and 75% side overlaps. The thermal imaging sensor was configured independently to acquire frames at 3 Hz to ensure image overlaps of 90-95%. Imagery data was stored on-board in respective memory cards for the two sensors. A calibrated reflectance panel (CRP, Micasense, Inc., WA, USA) was imaged before and after each flight and these reference images were used for radiometric calibration of the multispectral imagery ( Figure 2).
Total six UAS flights were conducted in summers of 2018 (two missions/site × 3 sites) during the dates and approximate times of the Landsat 7/8 overpass. UAS imagery for spearmint crop was collected 10 days before the first harvest (dataset-1) and 37 days before the second harvest (dataset-2). Imagery for the potato crop was acquired 72 days (dataset-3) and 48 days before harvest (DBH, dataset-4). Data for the alfalfa crop was collected 2 days before the first harvest (dataset-5) and 7 days before the second harvest (dataset-6).

Satellite-Based Imagery and Weather Data
Landsat 7/8 imagery datasets and 1-arc SRTM (Shuttle Radar Topography Mission) digital elevation models (DEM) were downloaded for the data collection days (Table 2). Weather data logged at every 15 min was downloaded from the nearest (1-5 km) stations (AgWeatherNet, Washington State University, Washington, WA, USA). All the data collected on a day were grouped to have a total of six datasets (see Table 2) ready for analysis.

Satellite-Based Imagery and Weather Data
Landsat 7/8 imagery datasets and 1-arc SRTM (Shuttle Radar Topography Mission) digital elevation models (DEM) were downloaded for the data collection days ( Table 2). Weather data logged at every 15 min was downloaded from the nearest (1-5 km) stations (AgWeatherNet, Washington State University, Washington, WA, USA). All the data collected on a day were grouped to have a total of six datasets (see Table 2) ready for analysis.

Preprocessing
UAS imagery was initially processed to obtain orthomosaics of the study site(s) through a series of image stitching operations ( Figure 3a) in a photogrammetry and mapping software (Pix4D Mapper, Pix4D, Inc., Lausanne, Switzerland). Independent surface temperature orthomosaic was then georeferenced and resampled to the multispectral orthomosaic (7 cm pixel −1 ) in Quantum Geographic Information System (QGIS) platform (ver.2.18.16, Open Source). The "Georeferencer" tool with "Thin Plate Spline" as the transformation type was used for georeferencing and the "Nearest Neighborhood" method was used for resampling. Five surfaces reflectance (B, G, R, RE, and NIR of image stitching operations ( Figure 3a) in a photogrammetry and mapping software (Pix4D Mapper, Pix4D, Inc., Lausanne, Switzerland). Independent surface temperature orthomosaic was then georeferenced and resampled to the multispectral orthomosaic (7 cm pixel −1 ) in Quantum Geographic Information System (QGIS) platform (ver.2.18.16, Open Source). The "Georeferencer" tool with "Thin Plate Spline" as the transformation type was used for georeferencing and the "Nearest Neighborhood" method was used for resampling. Five surfaces reflectance (B, G, R, RE, and NIR [ Figure 3b]), one surface temperature ( Figure 3c) and a DEM orthomosaic were obtained for each UAS mission.

METRIC Model Implementation
The METRIC model is detailed in Allen et al. [9]. Briefly, model computes Rn, G, H and residual LE within the blending height (100-200 m AGL) which can be converted to actual water loss to the atmosphere (ET). METRIC eliminates the need for accurate surface temperature (Ts) measurements. It uses ground weather station-based reference ET (ETr, alfalfa crop based) for internal calibration and reduction of biases related to satellite imagery. The internal calibration indexes near-surface temperature gradients to the radiometric Ts using extreme hot and cold anchor pixels. The model was primarily developed for agricultural fields and does not need crop specific inputs. This enhances METRIC's applicability to high-resolution RS data. Small UAS based RS data captures distinct features of soil, vegetation, or mixed regions that can improve internal calibration of METRIC model with respect to the study environment [26][27][28][29]33].
Conventional METRIC maps regional ET for 185×185 km scene size and all its assumptions might not necessarily apply at the field-scale. Therefore, METRIC model was implemented with three variant approaches. The resulting approaches are named as UAS-METRIC-1, -2, and -3

METRIC Model Implementation
The METRIC model is detailed in Allen et al. [9]. Briefly, model computes R n , G, H and residual LE within the blending height (100-200 m AGL) which can be converted to actual water loss to the atmosphere (ET). METRIC eliminates the need for accurate surface temperature (T s ) measurements. It uses ground weather station-based reference ET (ET r , alfalfa crop based) for internal calibration and reduction of biases related to satellite imagery. The internal calibration indexes near-surface temperature gradients to the radiometric T s using extreme hot and cold anchor pixels. The model was primarily developed for agricultural fields and does not need crop specific inputs. This enhances METRIC's applicability to high-resolution RS data. Small UAS based RS data captures distinct features of soil, vegetation, or mixed regions that can improve internal calibration of METRIC model with respect to the study environment [26][27][28][29]33].
Conventional METRIC maps regional ET for 185×185 km scene size and all its assumptions might not necessarily apply at the field-scale. Therefore, METRIC model was implemented with three variant approaches. The resulting approaches are named as UAS-METRIC-1, -2, and -3 (UASM-1, -2, and -3) with pertinent modification highlights summarized in Table 3. UASM-1 ( Figure 4) was implemented to observe relative differences only due to UAS imagery specific inputs that include (i) scene metadata from flight missions, (ii) surface albedo from on-board multispectral imager, and (iii) high-resolution DEM. All other equations were identical as the conventional METRIC. Similar to LM, leaf area index (LAI) in UASM-1 was calculated using soil adjusted vegetation index (SAVI) with background adjustment factor (L) of 0.1 for high crop cover [34,35]. UASM-2 was employed for two reasons: (i) possibility to reduce model input data needs from the grower perspective and (ii) to assess pertinent differences caused by the DEM typical to flat agricultural fields located in non-terrain regions. Implementation of UASM-2 was identical to UASM-1 except that UASM-2 used a flat DEM for local field conditions. This scenario forced surface slopes and aspects to zero and elevation as the mean of all pixels in the DEM. include (i) scene metadata from flight missions, (ii) surface albedo from on-board multispectral imager, and (iii) high-resolution DEM. All other equations were identical as the conventional METRIC. Similar to LM, leaf area index (LAI) in UASM-1 was calculated using soil adjusted vegetation index (SAVI) with background adjustment factor (L) of 0.1 for high crop cover [34,35]. UASM-2 was employed for two reasons: (i) possibility to reduce model input data needs from the grower perspective and (ii) to assess pertinent differences caused by the DEM typical to flat agricultural fields located in non-terrain regions. Implementation of UASM-2 was identical to UASM-1 except that UASM-2 used a flat DEM for local field conditions. This scenario forced surface slopes and aspects to zero and elevation as the mean of all pixels in the DEM. UASM-3 considered site-specific energy balance inputs typical to agricultural field plots imaged in this study. Pertinent modifications include (i) DEM as per field conditions (as in UASM-2), (ii) LAI calculation from spatial fraction canopy cover (FCC), (iii) incoming shortwave radiation (ISWR) from the nearest open field weather station, (iv) incoming longwave radiation (ILWR) calculation using air temperature (Ta), and (v) momentum roughness length (MRL) calculated without surface slope (S) adjustments. The LAI calculation in UASM-3 were based on the fact that surface features (i.e., soil, vegetation, and mixed) could be distinctively captured in every pixel [28][29][30] of high-resolution imagery. FCC was first derived from the normalized difference vegetation index (NDVI) [28,36,37] and soil was segmented using NDVI vegetation threshold mask (>0.3) to ensure LAI for the vegetation only. ISWR measurements from the nearest weather station may be assumed as actual shortwave radiation on the surface for the fact that these are open-field stations installed at 2 m AGL and located within 5 km radius from the study sites. Moreover, conventional ISWR calculation might get affected by cloud cover and other atmospheric dynamics [38,39]. Ta from nearest weather station was used to calculate ILWR and can be assumed constant for a field (Stefan Boltzmann's law, [40]).
For comparisons, Landsat imagery was processed through the conventional METRIC model as described in a "R package" (water, [41,42]), hereafter abbreviated as LM (Landsat-METRIC) approach. Level-2 surface reflectance products were downloaded for Landsat 8 while such products were calculated within the package for Landsat 7 datasets. Similar to the LM approach, all UASM UASM-3 considered site-specific energy balance inputs typical to agricultural field plots imaged in this study. Pertinent modifications include (i) DEM as per field conditions (as in UASM-2), (ii) LAI calculation from spatial fraction canopy cover (FCC), (iii) incoming shortwave radiation (ISWR) from the nearest open field weather station, (iv) incoming longwave radiation (ILWR) calculation using air temperature (T a ), and (v) momentum roughness length (MRL) calculated without surface slope (S) adjustments. The LAI calculation in UASM-3 were based on the fact that surface features (i.e., soil, vegetation, and mixed) could be distinctively captured in every pixel [28][29][30] of high-resolution imagery. FCC was first derived from the normalized difference vegetation index (NDVI) [28,36,37] and soil was segmented using NDVI vegetation threshold mask (>0.3) to ensure LAI for the vegetation only. ISWR measurements from the nearest weather station may be assumed as actual shortwave radiation on the surface for the fact that these are open-field stations installed at 2 m AGL and located within 5 km radius from the study sites. Moreover, conventional ISWR calculation might get affected by cloud cover and other atmospheric dynamics [38,39]. T a from nearest weather station was used to calculate ILWR and can be assumed constant for a field (Stefan Boltzmann's law, [40]).
For comparisons, Landsat imagery was processed through the conventional METRIC model as described in a "R package" (water, [41,42]), hereafter abbreviated as LM (Landsat-METRIC) approach. Level-2 surface reflectance products were downloaded for Landsat 8 while such products were calculated within the package for Landsat 7 datasets. Similar to the LM approach, all UASM approaches performed an energy balance within the blending height and used the conventional automated approach of selecting anchor pixels.

Output Comparisons
This study aimed at evaluating spatial ET maps from UAS based remote sensing (RS) data relative to satellite-based RS data for crops grown in commercial farming operations. As the conventional METRIC model is well-established for spatial ET mapping for its corroboration in over 25 countries and in almost all types of agroclimatic zones [22], the additional field validation was deemed not necessary, i.e., outside the scope of the study. Thus, we did not conduct ground-reference ET measurements.  Since the energy balance computations with LM approach were performed for images of 180 km × 180 km dimension, intermediate and final output maps were clipped to the study site areas covered by the UAS flights. The mean, standard deviation (SD), and coefficients of variation (CV, %) to assess spatial variability capturing potentials of ET pixels from LM (ET LM ) and UASM-1, -2, and -3 approaches (ET UASM ) were calculated. Percentage absolute departures of mean ET (ET dep,abs , Equation (1)) from all UAS-based approaches with respect to that from the LM approach were also calculated. Next, ET maps from all UAS-based approaches were resampled to the resolution equal to ET maps from LM approach (30 m pixel −1 ) using the "Nearest Neighborhood" method and region of interest (ROI) samples were randomly extracted. Sample ET from these ROIs were compared using root mean square difference (RMSD) and Pearson's linear correlation (r) at 5% significance (RStudio, Inc., Boston, MA, USA).

Results
Major energy balance computations from the UASM and LM approaches are presented in the following subsections. Results pertinent to spearmint (

Crop Vigor
High spatial variability in crop vigor was captured by the UAS imagery. Soil background adjustment factor (L) of 0.1 as in the LM calculated similar SAVI with the UAS imagery for selected crops [34,35]. Higher NDVI compared to SAVI for some datasets may be due to the saturation effects ignored in SAVI [43]. Soil-segmented FCC mapped similar LAI from UASM-3 as from the LM approach. MRL estimates from UAS and Landsat imagery were also similar. Slight differences in MRL from UASM-1 and -2 approaches were mainly due to surface slope factor in the case of UASM-1 approach and that between UASM-1 and -3 approaches were due to LAI and DEM inputs.

Net Radiation and Soil Heat Flux
The ISWR from LM, UASM-1, and -2 approaches were similar and slight differences were due to the DEM inputs. ISWR from UASM-3 approach was different due to its direct measurement source (nearest weather station). Such ISWR measurements may be assumed actual compared to the standard computations for the weather stations installed near study sites (within 5 km) approximately at 2 m height [38,39]. Higher ILWR estimates from the LM approach compared to UASM-1 and -2, were primarily due to higher T s . Minor ILWR differences between LM and UASM-1 or -2 approaches would be due to slightly different DEMs used to calculate atmospheric transmissivities. ILWR computation in UASM-3 approach used T a as the input (Stefan Boltzmann's law) instead of the T s [40] and no variation for use of a single value may be assumed realistic at field level. OLWR differences between LM and UASM approaches were primarily due to differences in T s from pertinent imaging systems. Similar OLWR from UASM-1, -2, and -3 approaches showed no effect of LAI differences that contribute only 1% to the surface emissivity calculations. R n and G estimate differences between UASM and LM approaches can be attributed to differences between ISWR, ILWR, OLWR and surface albedo computations. While no R n and G differences between UASM-1 and UASM-2 approaches were due to negligible differences between ISWR, ILWR, and OLWR. Such differences between UASM-1 and UASM-3 approaches were due to the ISWR and ILWR differences. Since high LAI was observed in all imaging campaigns, slightly high variations in OLWR and R n may be expected during initial days of plantation (Low LAI).

Sensible Heat Flux
Differences between H estimates from LM and UASM approaches were primarily due to T s differences, and secondarily due to different internal calibrations. As the Landsat 7/8 imagery has coarse resolution (~30 m pixel −1 ), it might not be readily possible to find uniform hot and cold anchor pixels. This demands for additional flexibility in anchor pixel identification also suggested by Jaafar et al. [42]. Contrarily, it was easier with the UAS-based approaches (~7 m pixel −1 ) to identify anchor pixels within the imaged field area (Table 7). Irrigated crops would have also facilitated identification of highly transpiring cold anchor pixels alike the reference alfalfa. H differences between UASM-1, -2, and -3 could be due to LAIs, MRLs, R n , and G differences in the different cold pixels that resulted in different internal calibrations.

Daily Evapotranspiration
Similar ET was mapped from all UASM approaches compared to the LM approach (Figures 5-9). However, UASM approaches showed larger potential for assessing spatial variations in field ET maps ( Figure 10) due to large number of pixels [21]. The ET differences between LM, UASM-1, and -3 approaches were large for dataset-5 (alfalfa field) and dataset-3 (potato field) primarily due to the internal calibration differences, H and R n differences. The ET r F for alfalfa crop in dataset-6 was high compared to dataset-5 showing more matured crop at 2 DBH than 7 DBH. Notable ET r F differences between LM and UASM approaches could also be due to accumulated differences between the intermediate outputs. ET r Fs for selected crops did not exceed 1.1, complying to the cold anchor pixel assumption (1.05 × ET r , [9]). UASM-3 approach (METRIC model modified to local conditions) can also provide reliable geospatial ET maps as LM. As also reported by Paul et al. [33], G and H differences between LM and UAS approaches did not influence daily ET maps for irrigated field crops due to small magnitudes. UASM approaches indicate their versatility for geospatial ET mapping for a range of similar crops as spearmint, potato and alfalfa (r: 0.8-0.9). However, additional parametrization may be needed for different crop physiologies [21,27]. An average 10% uncertainty of UASM approaches for point scale ET might be acceptable as discussed in prior research studies that used similar UAS imagery with METRIC model [29,33]. Flight missions within solar noon ± 2 h could have produced comparable and reasonable ET r F used for 24-h [44].  Figure 10. Spatial daily ET variation potential assessments from LM, UASM-1, -2, and -3 approaches.

Conclusions
Agricultural growers are highly dependent on generalized crop coefficients or point measurements (soil moisture, sap flow, etc.) for irrigation scheduling and are restrictive towards satellite-based remote sensing approaches due to low spatiotemporal resolution and cloud cover issues. Considering the critical need for geospatial ET mapping towards precision irrigation management, ET of irrigated field crops was mapped using very high resolution multispectral and thermal infrared imagery. Three variants of standard METRIC energy balance model modified as per UAS and local conditions were used to process UAS imagery and their performance was evaluated relative to standard LM approach.

Conclusions
Agricultural growers are highly dependent on generalized crop coefficients or point measurements (soil moisture, sap flow, etc.) for irrigation scheduling and are restrictive towards satellite-based remote sensing approaches due to low spatiotemporal resolution and cloud cover issues. Considering the critical need for geospatial ET mapping towards precision irrigation management, ET of irrigated field crops was mapped using very high resolution multispectral and thermal infrared imagery. Three variants of standard METRIC energy balance model modified as per UAS and local conditions were used to process UAS imagery and their performance was evaluated relative to standard LM approach.
Increased spatial resolution did not influence the mean field-scale ET (low RMSD: 0.2-1.0 mm day −1 , ET dep,abs : 0.4-21.2%, and r: 0.7-0.9). However, high spatial ET variation potential was assessed in UAS derived maps (6.7-24.3%) compared to the Landsat satellite derived maps (2.1-11.2%). Variants of UASM models performed very similar to the LM approach for irrigated field crops; spearmint (ET dep,abs : 0.4-6.2% and r: 0.7-0.8), potato (ET dep,abs : 0.9-21.2% and r: 0.7-0.8), and alfalfa (ET dep,abs : 7.2-9.9% and r: 0.8-0.9). UASM-2 showed highest similarity with the LM approach (RMSD: 0.5-0.9 mm day −1 and r: 0.7-0.9) for very similar energy balance outputs. ET differences between LM and UAS based approaches were majorly due to the internal calibrations differences that could have been affected slightly by the spatial resolution. UASM-1 could be advantageous for agricultural fields with elevation variability and resulting temperature differences due to lapse rates. UASM-1 uses high-resolution UAS imagery derived DEM inputs that can be accurately compared to the satellite based DEMs. UASM-2 could be advantageous for the crop fields where spatial variability in ground elevation is negligible, as it also reduces additional DEM data requirement for ET mapping. UASM-3 uses actual meteorological parameters of ISWR and T a (for ILWR) measured by the automated weather station near the field site. UASM-3 could also be advantageous for ensuring LAI calculation for the areas with vegetation only, as supported by high spatial resolution to distinctly capture soil, vegetation and mixed pixels. UASM-1 and UASM-3 could also be merged for field-relevant energy balance inputs such as meteorological parameters and ground elevation variabilities to obtain more realistic crop water use maps.
ET mapping at high spatiotemporal resolution provides more control over timely monitoring of spatial crop water requirements where restrictions of cloud cover interference and imagery acquisition at 16-day interval could be avoided unlike LM approach. LM based approaches may assist in reviewal of seasonal water distribution rights while UASM approaches may assist in irrigation water savings through grower friendly and on-demand site-specific irrigation prescription maps/tools.