Spatial Estimation of the Latent Heat Flux in a Tropical Dry Forest by Using Unmanned Aerial Vehicles

Abstract

In this paper, we address the retrieval of spatially distributed latent heat flux ($\lambda$E) over a tropical dry forest using multi-spectral and thermal unmanned aerial vehicle (UAV) imagery. The study was carried out in the Santa Rosa National Park Environmental Monitoring Super-Site, Costa Rica, in June 2016. The triangle method was used to derive $\lambda$E from the UAV imagery and the results were compared to $\lambda$E measurements of an eddy covariance system within the coincident eddy flux tower footprint. The tower footprint was derived using a two-dimensional parameterization model for flux footprint prediction. The comparisons with the flux tower measurements showed a mean relative difference of 10.98% with a slight overestimation of the UAV-based flux retrievals by nearly 7.7 Wm${}^{-2}$. The results are in good agreement with satellite-based retrievals, as provided by the literature, for which the triangle method was initially developed and mostly used so far. This study proved to be a promising approach for transferring the triangle method to UAV imagery in ecosystems such as tropical dry forests. With the presented approach, new details in spatially distributed latent heat flux estimates at ultra-high resolution are now possible, thereby potentially closing the gap in spatial resolution between satellites and flux towers. Even more, it allows tracing the latent heat flux from single trees at leaf level. Besides, this approach also opens new perspectives for the monitoring of latent heat fluxes in tropical dry forests.

1. Introduction

In recent years remote sensing has become an important element to the provisioning of spatially distributed input data for Earth system models. With various global initiatives, such as the Essential Climate Variables (ECV) defined by the World Meteorological Organization (WMO) and the European Space Agency’s Climate Change Initiative (ESA-CCI), global monitoring programs through satellite data are envisaged. Space-borne remote sensing (RS) platforms contribute to global environmental monitoring. However, limitations of satellite systems due to technical capabilities set the boundaries for the temporal and spatial resolution of many data sets [1]. This is a limitation for the validation of emerging remote sensing products with wide disparities among ground truth data [2], especially where there are significant inconsistencies between sensor and ground-based information. This is evident for ECVs such as evapotranspiration; i.e., the latent heat flux ($\lambda$E). $\lambda$E is a key element for improving our understanding of climate-mediated changes of energy and water cycling at global to local scales. Knowledge about the spatial distribution of $\lambda$E is essential for many applications, such as water resource management, agriculture, and forestry. Abnormal changes in $\lambda$E can be caused by plant water stress and reduced plant vitality. Therefore, spatially explicit estimates of $\lambda$E are of high interest for decision makers, especially in tropical dry forests, where climate change might have a negative impact on tree vitality [3].

Spatially distributed $\lambda$E estimates are usually derived from satellite-based (e.g., MODIS, MERIS, or Landsat) data products that have pixel resolutions ranging from 90 to 1000 m [2,4,5]. These data products are often validated with measurements of $\lambda$E from eddy covariance (EC) systems. However, measurements of $\lambda$E using the common EC approach might be derived for a much smaller eddy flux footprint than the sizes of one of these pixels. This difference in spatial resolution could possibly introduce large uncertainties due to spatial scaling effects [6,7]. In addition, measurements of EC-based $\lambda$E represent an average value for the total predicted flux footprint rather than plant-specific information, while the flux footprint prediction in turn is a function of wind speed, wind direction, and sensor height above a given surface [8]. Therefore, there is a lack of spatially explicit $\lambda$E information of high spatial resolution that would enable the monitoring of latent heat fluxes of forest canopies, plant individuals, or plant compartments over larger areas than those measured in-situ, which could also be used to validate satellite-based $\lambda$E estimates.

One emerging source of remote sensing information with significant potential for the calibration and validation of ECVs is near-surface observation from, e.g., multi-spectral sensors attached to autonomous unmanned aerial vehicles (UAVs). According to [9], UAVs can be used for closing the gap in scale between field and satellite-based observations by offering more flexible spatial, temporal, and spectral resolutions, potentially allowing the retrieval of information on single individual plants or land-surfaces at leaf scale. In addition, UAVs provide the unique opportunity to integrate thermal infrared (TIR) and visible and near-infrared (VIS-NIR) sensors into one single platform at low cost. This enables larger scale mapping (<0.1 km${}^2$) of thermal properties of sites that are aimed to support calibration and validation efforts of emerging satellites. Recently, studies have addressed the derivation of latent heat flux estimates from UAVs by using different approaches. For example, reference [10] estimated $\lambda$E based on a two-source energy balance model by integrating a thermal infrared sensor with a fixed-wing UAV. The study shows the feasibility of using high-resolution UAV products instead of satellite-based products as model input. The study of [11] used UAV-based multi-spectral and thermal infrared data to develop a high resolution spatial evapotranspiration system. Within their study, they retrieved actual crop evapotranspiration from five different evapotranspiration models. The best results were achieved when applying input variables based on vegetation indices such as the normalized difference vegetation index (NDVI). In contrast, refs. [10,12] used only thermal UAV imagery to derive the latent heat flux.

Usually, the models applied to retrieve $\lambda$E from remote sensing data can be grouped into contextual/scene-based models or pixel based models [13]. While the models of the first group require the reflectance and/or land surface temperature (LST) information from the whole image scene to derive a so-called hot/dry edge and cold/wet edge to constrain the retrievals, models of the latter group do not consider such information and thus need more input parameters [13]. However, contextual models need a broad variety of different reflectance values to define those dry and wet edges, which is usually a result of surface heterogeneity [14]. The triangle method, first introduced by [15], is a well-established method representing the contextual models, while the “two source energy balance” (TSEB) [16,17] and the “surface energy balance algorithm for land” (SEBAL) [18] represent the pixel-based models. According to [19], contextual models have been extensively tested over large areas with sufficient surface heterogeneity, which restricted their applicability to coarser resolution satellite data. This in turn impedes the validation due to the mismatch of in-situ footprints and the sensor’s spatial resolution, as elaborated by [2]. Using high-resolution UAV data would allow for a better understanding of the spatial distribution of $\lambda$E.

Recently, ref. [13] tested three different $\lambda$E flux retrieval approaches from thermal and red-green-blue (RGB) UAV imagery. Although they did not include the near-infrared (NIR) domain, they concluded that simpler models such as the one-source or the triangle method can produce similar results compared to more complex inversion schemes such as a two-source energy balance model [13]. In addition, reference [20] used a UAV-based triangle method approach to retrieve root zone soil moisture information. In [21], the authors used the same UAV and triangle method approach to retrieve land surface energy, water, and CO${}_2$ fluxes at very high resolution. However, development and validation of such models applied to UAV imagery taken in different ecosystems such as tropical dry forests are still ongoing research and pending.

To better understand spatially explicit $\lambda$E values of tropical dry forests at ultra-high resolution, we modeled $\lambda$E using UAV-based visible and thermal reflectance information and the triangle method approach. Modeled $\lambda$E results were compared to EC measurements of the $\lambda$E flux. The objectives of the presented study are

• To derive spatially distributed estimates of $\lambda$E;
• To show the performance of a contextual method (triangle method [22]) in high-resolution imagery in tropical dry forests;
• To trace the different contributions of $\lambda$E at tree and branch level in a tropical dry forest environment.

2. Methods

2.1. Study Site

This study was carried out at the Santa Rosa National Park Environmental Monitoring Super-Site (SRNP-EMSS, Figure 1b). The SRNP-EMSS is located in the province of Guanacaste in the northwest of Costa Rica. The area has a mean annual temperature of 25 °C and receives a mean annual precipitation of 1750 mm. The dry season (December to May) remains rain free and marks the leaf-off period for about 90% of all trees. The tropical dry forest of Costa Rica, one of the most endangered ecosystems, consists of a mosaic of new-growth to old-growth forests depending on the year of agricultural abandonment, and as such has various forest structural properties and species [23]. The trees in the study area are on average about 18 m tall and have a relative age since abandonment of 30–35 years. The SRNP-EMSS hosts two EC systems plus four optical phenology stations, which have been in operation since 2013 [3]. The two EC systems have a distance of 1200 m. In addition, about 50 sensors are distributed within the footprints of both EC systems, measuring photosynthetic active radiation, temperature, relative humidity, and soil moisture. The UAV imagery that was acquired for this study is only coincident to one of these EC systems, due to difficulties of accessing the site and launching the UAV.

2.2. UAV and Camera Payload

2.2.1. UAV System

For this study, we used a Rotorkonzept-RKM 8X octocopter. The RKM 8X was designed by the German company Rotorkonzept and is a flexible lightweight UAV which can carry payloads up to 1100 g. The UAV is controlled by a Pixhawk flight controller which is an open hardware project, making it very flexible and extendable and allows to fly according to way-points. To fit the needs of simultaneously carrying two camera systems, a special two-axis stabilizing gimbal was built. The gimbal corrects for yaw and roll effects by using brush-less engines that allow an always nadir looking imaging geometry for both cameras. The mean flight duration of the RKM-8X is around 12 min with a payload of approximately 1000 g, depending on wind and humidity conditions. The nominal flight speed in standard flight mode is 5 m/s. The multi-spectral and thermal infrared cameras, described below, were mounted on the gimbal and complete the UAV system for the field experiments, as shown in Figure 1a.

UAV acquisitions were carried out at the 29th of June 2016 within four consecutive flights in order to cover a larger area of roughly 0.1 km${}^2$. One of the four flights was directly coincident to the EC flux footprint and was carried out from 4:00 pm to 4:10 pm. The UAV images were acquired in autopilot mode in an altitude of 100 m above the canopy with a forward overlap of 80%, a sidelap of 70% between single images, and in nadir looking geometry. Image acquisition was carried out during constant overcast conditions.

2.2.2. Multi-Spectral Camera System

The Micasense RedEdge camera is a multi-rig camera allowing for the simultaneous acquisition of 5 channels. These 5 spectral bands are set up as an array of 2-1-2 containing the sensors from upper left to lower right: blue–near-infrared–rededge (middle)–green–red, covering the electromagnetic wavelength spectrum from visible to near-infrared. The band sensitivities differ from one to another, which needs to be addressed in the image pre-processing. Specifications of the Micasense RedEdge camera are shown in

Table 1. Specification of the Micasense RedEdge and FLIR TAU2 640 camera systems. GSD: ground sampling distance; HFOV: horizontal field Of view.

Parameter	Micasense RedEdge	FLIR TAU2 640
Dimensions [m]	0.121 × 0.066 × 0.046	0.045 × 0.045 × 0.060
Weight	150 g	95 g
Spectral bands
(Center wavelength and bandwidth per channel [nm])	Blue 475 (20), green 560 (20), red 668 (10), red edge 717 (10), near IR 840 (40), narrow band	7500–13,500, broad band
Resolution	1280 × 980 pixels	640 × 512 pixels
Focal length	0.0055 m	0.013 m
GSD at 100 m above ground	0.07 m (per band)	0.13
Field of View	47.2° HFOV	45° × 37° 1.308 mr

. Images are acquired with 12-bit depth and stored and saved as 16-bit tiff images for each channel separately in single files for further processing. To obtain calibrated images from the RedEdge, Micasense’s original calibration panel with highly Lambertian properties and known spectral albedo values was imaged before each flight and used for calibration inputs in the calibration procedure, as described in Section 2.4.1.

2.2.3. Thermal Camera System

The thermal infrared imagery was acquired using a FLIR TAU2 640. The TAU2 cores are manufactured to capture thermal video streams. In order to retrieve still imagery, the FLIR core was extended with a backend, designed by the German TEAX GmbH. The backend grabs the single frames from the video stream and converts it to a still image as a raw binary file of 14-bit.

2.3. Data

During the experiments, the following data

Table 2. Acquired and used data within the Santa Rosa National Park Environmental Monitoring Super-Site.

Parameter	Sensor	Remarks
LST	UAV-FLIR TAU2 640	Images were acquired at an altitude of 100 m above ground (canopy layer)
VIS-NIR	UAV-Micasense	Images were acquired at an altitude of 100 m above ground (canopy layer)
$R_n$	Kipp & Zonen CNR4 Net Radiometer	Net-radiation was measured every 10 min at the EC system
G	HukseFlux HFP01SC	Soil heat flux was measured every 10 min within the footprint of the EC system
$T_{air}$	Decagon VP-3	Air temperature was measured every 10 min at the EC system and within the EC footprint
RH	Decagon VP-3	Relative humidity was measured every 10 min at the EC system and within the EC footprint
$\lambda$E	Licor LI-7500DS & Gill Windmaster	Fluxes were measured at 10 Hz and integrated over a time span of 30 min and analyzed using Licor’s EddyPro software

were acquired by the UAV, the EC system, and the wireless sensor network, as described in Section 2.1.

2.4. Camera Calibration

To convert the raw digital numbers (DNs) of the acquired imagery into real measurements, a calibration of both cameras was carried out. In the following, the approaches for the two camera systems are outlined in detail.

2.4.1. Multi-Spectral Camera System

The Micasense comes pre-calibrated and calibration values and constants were provided in the Exif-Header section of each raw tiff-file per band. In order to convert the raw 16-bit DNs into reflectance values, a calibration procedure was carried out according to Micasense recommendations. The calibration procedure included correcting for vignetting effects of the single camera lenses, black current compensation, and the conversion from radiance to reflectance by using a Lambertian calibration target with known spectral albedo values. For each flight, the calibration target was imaged with the Micasense camera assuming constant conditions of the atmosphere during the flight [24].

2.4.2. Thermal Camera System

The calibration of the FLIR TAU2 640 is relatively straightforward as it comes already calibrated off the shelf. A linear equation was used that allows the conversion of the DNs registered by the camera into degrees Celsius (°C) depending on the expected land surface temperature range. Equation (1) describes the conversion of DNs to °C for a so-called low gain case. This calibration was applied to all UAV digital terrain models (DTM) after orthomosaic generation. It is notable that the original FLIR camera stores the data in 14-bit DN, while the TEAX backend stores the imagery in 16-bit tiff format. Therefore, a conversion factor $\frac{2^{16}-1}{2^{14}-1}$ was implemented in Equation (1):

$$LST=\frac{DN}{\frac{2^{16}-1}{2^{14}-1}}·0.04-273.15$$

(1)

Under cloudy or overcast conditions the effect of emissivity on the thermal signal does not show any significant contribution and can be neglected according to [25]. Thus, no corrections for emissivity were carried out in this study. On average, the atmospheric column between the forest canopy and the UAV was 100 m in height above the canopy layer. The narrow atmospheric column minimizes the effects of scattering and absorption of radiation by atmospheric gases and aerosols on potential differences in the radiance that is reflected by the land surface and the radiance which is measured by the sensor [26,27]. Therefore, we decided to not atmospherically correct the thermal infrared data in this study. This decision was supported by the short flight time of less than 5 min, during which no changes in air temperature were observed, as indicated by simultaneous field measurements.

2.5. Orthomosaic Generation

The UAV orthomosaics were generated by using the individually calibrated images in the Agisoft Photoscan Professional software package (v1.3.2) (Agisoft LLC, St. Petersburg, Russia). Note that we used already radiometrically-calibrated images to have full control over the calibration procedure. Consequently, for the stitching of the single imagery, we used Agisoft without any further radiometric altering/filtering. To enhance the positional accuracy of the final orthomosaic, five full ground control points (GCPs) were measured in the field using a differential GPS (Trimble GeoXT6000, (Trimble Inc., Sunnyvale, CA, USA)) with an average horizontal error of 0.50 m and an average vertical accuracy of 0.54 m. These marked GCPs were included in the Agisoft workflow and enhanced the accuracy of the aligned image positioning significantly. To obtain a common resolution between the VIS-NIR and TIR orthomosaics, the derived TIR orthomosaic was oversampled by a factor of roughly two.

2.6. Derivation of Latent Heat Flux Using the Semi-Empirical Triangle Method

The triangle method introduced by [15] and further developed by [22,28] is a widely used approach [14,29,30,31,32] based on the Priestley–Taylor Equation (2) for the estimation of the $\lambda$E flux [Wm${}^{-2}$] with the energy balance equation of the following form:

$$\lambda E=\varphi \left[(R_n-G)\frac{\mathsf{\Delta }}{\mathsf{\Delta }+\gamma }\right]$$

(2)

where $\varphi$ is a dynamical equivalent to the dimensionless Priestley–Taylor coefficient [-] [22,29], replacing the resistance terms of the Penman–Monteith Equation. $\varphi$ is derived in the following by Equation (6). $R_n$ is the net radiation [Wm${}^{-2}$], G is the soil heat flux [Wm${}^{-2}$], $\mathsf{\Delta }$ is the slope of the saturated pressure curve derived after [33] [kPa K${}^{-1}$], and $\gamma$ is the psychromatic constant derived after [34] in [kPa K${}^{-1}$]. The theory implies that a high $\lambda$E flux provides cooling of the surface area resulting in a cooler sunlit vegetation surface compared to a surface consisting of bare soil under similar conditions [14]. In comparison to the well-known SEBAL approach [16,18], fewer input variables are required and thus fewer sources of uncertainties are obtained. Although the triangle method approach was performed using predominantly satellite remote sensing products as model input variables [22,29,30,31,32], this study uses the approach with input variables derived from the aforementioned UAV-based VIS-NIR and TIR data sets. The VIS-NIR orthomosaic is hereby used to obtain the vegetation index, while the TIR orthomosaic provides the LST information. To derive $\lambda$E according to the approach proposed by [22], the NDVI is typically used. However, the near-infrared sensor of the Micasense camera showed high saturation effects across the acquired UAV imagery due to the vitality of a majority of the trees. For this study, the NDVI can not represent the fine spatial differences of vegetation greenness and results in a homogeneous land surface. A heterogeneous observation area is, however, one main requirement for the application of the triangle method [14,22,29]. Therefore, the NDVI was replaced by the rededge NDVI (RENDVI, [35,36]) which is defined as:

$$RENDVI=\frac{(rededge-red)}{(rededge+red)}$$

(3)

where the rededge band is used instead of the NIR band due to the sensitivity of the rededge band to small changes in canopy foliage content, gap fraction, and vegetation senescence. According to [36], the allowed ranges of wavelengths within the RENDVI are 670–740 nm for the rededge band and 663–673 nm for the red band, respectively. With 712–722 nm and 600–699 nm spectral bandwidth, the Micasense’s rededge and red band perfectly match the allowed ranges. With the spatially distributed RENDVI and LST maps, the Priestley–Taylor coefficient $\varphi$ is derived. For the calibration of the semi-empirical model, this parameter can be scaled between ${\varphi }_{min}$ and ${\varphi }_{max}$ by linearly interpolating $RENDV{I}_{min}$ and $RENDV{I}_{max}$ (adjusted after [33]). The so obtained range of $\varphi$ represents the maximum and minimum values for $\lambda$E such that the minimal value of $\lambda$E implies a minimum value of $\varphi$ [29]. Potential evapotranspiration (ET${}_p$) conditions occur when no plant-water stress is present, and therefore $\lambda$E can reach maximal values [22]. This is modeled by the wet-edge representing a maximum value of $\varphi$ and thus of $\lambda$E. The true dry-edge is obtained when $\lambda$E is absent and thus zero, whereas an observed dry-edge represents limited $\lambda$E values with limited ${\varphi }_{min}$. In order to model increased $\lambda$E over higher vegetated areas in comparison to areas of sparse vegetation cover, a non-linear scaling of $\varphi$ is performed with the following equation:

$${\varphi }_{(i,min)}={\varphi }_{max}{[\frac{{RENDVI}_i-{RENDVI}_{min}}{{RENDVI}_{max}-{RENDVI}_{min}}]}^2$$

(4)

where ${\varphi }_{i,min}$ represents ${\varphi }_{min}$ for each pixel i and RENDVI value, and ${\varphi }_{max}$ is according to [22] given by:

$${\varphi }_{max}=\frac{\mathsf{\Delta }+\gamma }{\mathsf{\Delta }}$$

(5)

Based on this equation, the next step spreads the RENDVI-LST triangular space to the minimum and maximum extent of evaporative cooling, represented by ${\varphi }_{i,min}$, $LS{T}_{i,max}$, and ${\varphi }_{i,max}$, $LS{T}_{i,min}$ for each RENDVI value i respectively [29,30]. The specific ${\varphi }_i$ value is obtained by:

$${\varphi }_{\left(i\right)}=\frac{LS{T}_{(i,max)}-T_i}{LS{T}_{(i,max)}-LS{T}_{(i,min)}}·({\varphi }_{max}-{\varphi }_{(i,min)})+{\varphi }_{(i,min)}$$

(6)

where ${\varphi }_i$ represents the Priestley–Taylor coefficient for a specific RENDVI value i; LST${}_{(i,max)}$ and LST${}_{(i,min)}$ are the corresponding minimum and maximum LSTs; and $T_i$ is the observed surface temperature at pixel i.

The triangle method was conducted within R [37] using the R-package “Triangle Method”—which was freely provided by the author of the package [38] on demand. The additional inputs of air temperature, soil heat flux, and net radiation were provided by the EC system and the wireless-sensor network [3]. Although the EC system was covered by one flight line only, $\lambda$E was derived for the one full scene consisting of all four flight lines in order to capture a higher land-cover heterogeneity.

2.7. $\lambda$E-Flux Measurements and Footprint Estimation

Flux equipment was installed on a horizontal boom at a height of 35 m above the forest floor and in the footprint of the tower; e.g., soil heat flux. Different corrections for, e.g., time-lags, high frequency losses, and outlier detection were applied according to [3], and fluxes (CO${}_2$, H${}_2$O, momentum, and energy) were calculated as 30-min block averages of the 10 Hz single measurement intervals. According to [3] and after additional testing, a 30-min averaging of the 10 Hz intervals allowed an adequate capture of low-frequency fluxes. For this study, only the H${}_2$O, momentum, and energy fluxes were used. $\lambda$E was calculated as the product of the latent heat of vaporization and the water vapor flux, expressed in energy units (Wm${}^{-2}$) [3].

The EC footprint was derived using a two-dimensional parameterization of flux footprint prediction (FFP), developed by [39] as a one-dimensional approach and elaborated by [8] as a two-dimensional parameterization. The FFP models the spatial extents, widths, and shapes of weighted eddy fluxes per pixel in 1 m × 1 m increments using the measurements of the installed EC system within the test-site as inputs. Figure 2a,b) shows the output of the FFP. The FFP output is a grid of the spatio-temporal EC flux footprint, where each 1 m${}^2$ pixel is weighted based on its probability of contributing to the total flux measurement, the so-called probability density function (PDF). The extent and location of an area that contributes with a certain percentage to the total measured flux can be obtained [7,8]. The x and y coordinates in the model intern footprint coordinate system describe the distance to the source of measurement (EC system) in the specific direction. The value of the pixel in m${}^2$ describes the area cumulated under the integral in the specific class between 0 and 0.99 with 0.1 increments based on the PDF. The model was applied by using the open-source FFP algorithm in R, which is provided by [8].

Figure 2 illustrates the EC footprint area divided into its PDF contributions, derived with a 30 min average sampling rate temporally coincident to the UAV flight overpass. The upper two images of Figure 2 visualize the distribution in a two-dimensional and three-dimensional directions. The color scheme and the extent in z direction respectively, indicate the contribution to the total measured flux within the specific weighted class of probability, precisely a part of the area [m${}^2$] below the integral of the PDF. Due to the limited battery capacity of the UAV, the observation area covers the EC footprint only for a footprint area where 63% of the total measured flux originates (red and yellow ellipses in Figure 2a). Thus, it represents more than half of the complete measured flux (Figure 2c,d). As the rest of the flux contributing footprint area lies outside the UAV coverage, only this area is used for validation purposes. Therefore, the modeled $\lambda$E raster file was cropped to the size of the specific footprint of the time and day.

2.8. Validation of $\lambda$E Retrievals

For the validation of the UAV based $\lambda$E retrievals, the measurements of the EC systems were used. In particular, as mentioned in Section 2.7, the mean over all the pixels within the considered footprint was benchmarked against the measurements of the EC system. For the accuracy assessment, the mean absolute difference and the mean relative difference were calculated.

3. Results

Figure 3 shows the VIS-NIR orthomosaic in true-color representation for the complete study area and for the region of interest (Figure 4).

The spatial resolution of the VIS-NIR orthomosaic is 0.08 m and consists in total of 437 single images that were aligned successfully. While 10% of the total 437 TIR images could not be aligned due to insufficient contrast in the images, we were able to close specific alignment gaps due to the applied forward overlap of 80% during the UAV image acquisition. Consequently, VIS-NIR and thermal orthomosaics have been produced without any spatial gaps. With a positional accuracy of below one pixel (less than two times the spatial resolution), the VIS-NIR orthomosaic is of high qualtiy [40,41]. The co-registration of the VIS-NIR and TIR orthomosaics through cross-correlation resulted in 93 tie-points and after applying a polynomial shift in a RMSE of 0.0014 m.

In general, different landscape features, such as roads, clearings, single trees, and even single branches, are visible. For further analysis, the generated orthomosaics have been cropped to the estimated footprint of the EC system, as shown in Figure 5 and Figure 6.

The simultaneously measured $\lambda$E flux of the EC-station is 70.2 Wm${}^{-2}$. With a mean absolute difference of 7.7 Wm${}^{-2}$ of the modeled UAV-based $\lambda$E, the triangle method approach slightly overestimates the measured flux with a mean relative difference of 10.98%.

Figure 7 shows the corresponding RENDVI-LST triangle plot for the retrieved $\lambda$E estimates. Note that the plot was calculated from the total covered flight area to ensure a maximum heterogeneity of the LST and RENDVI pixel values [22,29]. In this plot, the clustering of the pixel values for RENDVI and LST around the right corner of the not ideally developed triangle is clearly visible. This indicates near-potential evapotranspiration conditions along the wet-edge with low LST and high RENDVI values and leads to a saturation of the spatially distributed $\lambda$E estimates.

The modeled spatially explicit $\lambda$E estimates are illustrated in Figure 8 and are cropped to the coincident EC footprint area (Figure 8). Due to the high spatial resolution of the UAV data sets, even single trees and individual branches can be detected. In general, lower values of $\lambda$E correspond to the open patches of soil and leafless tree branches while higher values correspond to fully developed and vital vegetation that show nearly potential evapotranspiration. This is indicated by the RENDVI-LST-plot from Figure 7, and in particular, by the clustering at low LST and high RENDVI values. However, shadowed forest areas, forest occlusion, and understory vegetation also indicate high values in $\lambda$E close to potential evaporation. This can be observed especially in the center of the $\lambda$E footprint. A zoomed-in analysis reveals the spatial differences due to the different trees and leaf conditions. Even the more productive parts of single trees can be distinguished from dry or dead branches, which reduce the latent heat flux significantly. In addition, branches of liana infested trees, which are very common in that area [42], can be detected by reduced values of $\lambda$E. As lianas develop more woody biomass [27] and thus consequently reduce the fraction of $\lambda$E of a canopy surface. This is especially visible in the eastern part of the footprint, where a larger liana infested tree shows reduced $\lambda$E values.

4. Discussion

Within this paper, UAV-based estimates of the latent heat flux ($\lambda$E) have been carried out using a VIS-NIR and a TIR orthomosaic in combination with a conceptual model—the triangle method-over a tropical dry forest in Costa Rica. Comparisons to actual in-situ measurements using an EC system showed a good accuracy of the approach with an absolute error of 7.7 Wm${}^{-2}$ and a relative error of 10.98%. With this approach, $\lambda$E estimates are traceable to individual tree, allowing one to distinguish between leafless branches of trees, dead trees, and lianas. Although further ground-truthing needs to be conducted, the proposed approach seems to be promising for advancing the analysis of liana distribution and development, a topic that is highly relevant with expected climate-mediated increases in liana abundance in tropical dry forests [23,42].

However, the retrieval of $\lambda$E by using the triangle method was challenging due to the relatively homogeneous surface area. As described in [22,29], the triangle method provides the best results when the imaged surface is heterogeneous. Consequently, the LST-RENDVI triangle did not develop properly in comparison to other studies that surveyed a more heterogeneous surface area. In addition, our results show that the forest mostly shaded the soil surface due to its dense canopy layer. Therefore, near-potential evaporation was derived for most of the area besides some gravel roads, trails, and dead trees or woody biomass. However, this might be an overestimation because inherently shaded areas are cooler, which potentially could result in the false indication of higher $\lambda$E values that lead to potential evapotranspiration. Nevertheless, the retrieved results of this study are in good accordance with other studies that use UAV or satellite-based estimates of $\lambda$E [4,13,22,29,32,33,43]. The proposed approach provides, therefore, a very flexible method for the monitoring of spatially-distributed $\lambda$E fluxes complementary to EC measurements with a clear focus on the characterization of spatially explicit sub-footprint fluxes.

A further critical step of the proposed study was the co-registration of the VIS-NIR and TIR imagery. For the alignment, a large number of GCPs were drived in order to provide exact matching of the two orthomosaics. GCPs thereby consisted predominantly of intersecting branches that showed different thermal signals compared to the background information. Since placing the GCPs in the orthomosaics is time-consuming, cost-intensive, and in the case of this study challenging in a forest canopy dominated area, reference [41] developed an automatic co-registration procedure for multi-temporal and high-resolution UAV image data sets, without the use of GCPs. By doing so, the results are comparable to those using the manual strategy, resulting in reliable results [41]. The development of such a fully automatic processing chain would have been an advantage. However, such an implementation for the integration of a VIS-NIR with a TIR camera is not yet available. Besides, future implementations of such an algorithm should also consider the x, y and z-shifts of the two camera systems in order to overcome the co-registration problem.

Other sources of uncertainty are the flight duration and therewith the changing environmental conditions and the saturation of the NDVI. To address the changing environmental conditions, a reflectance calibration of the VIS-NIR data is crucial. Reference [44] outline that a change in illumination can be neglected for flights less than 15 min. However, continuous reflectance measurements in a one-minute interval revealed changing conditions even within such a short flight period for our study. Therefore, refs. [24,44] concluded that best reflectance calibration results are achieved using the continuous panel method by either recording spectral responses of a Lambertian surface with a spectrometer throughout the flight time and across the area of interest, or by mounting an incoming irradiance measurement device, e.g., Micasense’s day-light-sensor (DLS), on top of the UAV. Reference [9] also emphasizes the importance of a continuous reflectance calibration method. So far, first attempts of integrating the DLS with our UAV system only showed small changes in reflectance responses. However, further field experiments need to be conducted to explore the possibilities of the DLS sensor. In order to address the saturation effects of the NDVI above fully developed forest canopies [43], which predominantly originate from the saturation of the NIR band of the Micasense camera, modifications of the triangle method were conducted by exchanging the NDVI with the RENDVI. This modification was successful for the described UAV system setup and is recommended for similar studies. Furthermore, the triangle method approach requires only a small amount of measured input parameters in addition to the UAV input data sets, which makes this approach very suitable for $\lambda$E retrievals and facilitates applications also in data scarce areas [31,32].

Looking forward, the applied and tested method can pave the way for the integration of UAV systems in a multi-data analysis approach. Hereby, the proposed approach could function as a gap closure between different spatial and temporal scales which could also be used for the validation of remote sensing based $\lambda$E products across different scales [9]. As shown in Figure 2, the estimated $\lambda$E flux of a coarse resolution remote sensing product (e.g., MODIS 500 m × 500 m resolution) can be highly biased if compared to a time series of flux footprint measurements. Assuming that such a MODIS pixel would integrate all pixels of one scene to a single value, the UAV-based approach reveals much more detail, exceeding spectral mixing analysis as well, which can be seen in Figure 8.

Furthermore, thanks to the high spatial resolution of the UAV imagery, the proposed approach can be regarded as a complementary source of information compared to traditional EC measurements. Because EC measurements perfectly match the field scale but only provide rough estimates of the flux contributing areas, new insights are provided in our study. Here we were able to trace the single individual flux contributions (e.g., trees, branches, leaves, gaps) to the total flux as measured by the EC system. Our study is therefore of great value for understanding spatially explicit flux contributions in ultra-high spatial resolution as shown in Figure 2 and Figure 8. The added value lies in particular in the characterization of the different forest constitutions, such as tree species or fraction of liana infested or dead trees.

5. Conclusions

In this paper, a study is introduced which used UAV-based VIS-NIR and TIR data in combination with a contextual semi-empirical model for the retrieval of spatially explicit latent heat flux information. It was demonstrated that $\lambda$E retrievals with a mean relative difference of 10.98% can be achieved using such a setup over a tropical dry forest. This is of great value, especially as the tropics are often prone to cloud coverage, and thus satellite-based retrievals of $\lambda$E are limited. Furthermore, UAV-based, spatially-distributed latent heat estimates derived with the semi-empirical triangle method over a tropical dry forest are so far unique and allow the provision of new insights into spatially-distributed fluxes and their corresponding footprints, especially when using advanced camera systems such as a multi-spectral and a thermal camera system flown simultaneously on a UAV. By doing so, this paper demonstrates the derivation of reliable $\lambda$E estimates in an ultra-high spatial resolution. The applied method closes the gap in scale between space-borne remote sensing products and ground-based flux measurements. Consequently, it allows tracing the spatially-distributed, individual flux contributions from individual trees, branches, and leaves to the total flux in much greater detail. This is evident, especially when compared to the currently available flux footprint estimation tools or satellite remote sensing products, and thus adds complementary details of great value to current flux studies. Finally, the proposed approach provides a future method for monitoring the changing tropical dry forests of the SRNP-EMSS in Costa Rica in terms of tracking changes in evapotranspiration due to an increase in the fraction of liana infested trees.