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

: In this paper, we address the retrieval of spatially distributed latent heat ﬂux ( λ 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 λ E from the UAV imagery and the results were compared to λ E measurements of an eddy covariance system within the coincident eddy ﬂux tower footprint. The tower footprint was derived using a two-dimensional parameterization model for ﬂux footprint prediction. The comparisons with the ﬂux tower measurements showed a mean relative difference of 10.98% with a slight overestimation of the UAV-based ﬂux 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 ﬂux estimates at ultra-high resolution are now possible, thereby potentially closing the gap in spatial resolution between satellites and ﬂux towers. Even more, it allows tracing the latent heat ﬂux from single trees at leaf level. Besides, this approach also opens new perspectives for the monitoring of latent heat ﬂuxes in tropical dry forests.


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 (λE). λE is a key element for improving our understanding of climate-mediated changes of energy and water cycling at global to Recently, [13] tested three different λ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 λE values of tropical dry forests at ultra-high resolution, we modeled λE using UAV-based visible and thermal reflectance information and the triangle method approach. Modeled λE results were compared to EC measurements of the λE flux. The objectives of the presented study are 1.
To derive spatially distributed estimates of λE; 2.
To show the performance of a contextual method (triangle method [22]) in high-resolution imagery in tropical dry forests; 3.
To trace the different contributions of λE at tree and branch level in a tropical dry forest environment.

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 1200m. 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.

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.

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. 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. 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.

Data
During the experiments, the following data ( Table 2) were acquired by the UAV, the EC system, and the wireless sensor network, as described in Section 2.1.

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.

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].

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 2 16 −1 2 14 −1 was implemented in Equation (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.

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.

Derivation of Latent Heat Flux Using the Semi-Empirical Triangle Method
The triangle method introduced by [15] and further developed by [28] and [22] is a widely used approach [14,[29][30][31][32] based on the Priestley-Taylor Equation (2) for the estimation of the λE flux [Wm −2 ] with the energy balance equation of the following form: where φ is a dynamical equivalent to the dimensionless Priestley-Taylor coefficient [-] [22,29], replacing the resistance terms of the Penman-Monteith Equation. φ is derived in the following by Equation (6). R n is the net radiation [Wm −2 ], G is the soil heat flux [Wm −2 ], ∆ is the slope of the saturated pressure curve derived after [33] [kPa K −1 ], and γ is the psychromatic constant derived after [34] in [kPa K −1 ]. The theory implies that a high λ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 λ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: 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 φ is derived. For the calibration of the semi-empirical model, this parameter can be scaled between φ min and φ max by linearly interpolating RENDV I min and RENDV I max (adjusted after [33]). The so obtained range of φ represents the maximum and minimum values for λE such that the minimal value of λE implies a minimum value of φ [29]. Potential evapotranspiration (ET p ) conditions occur when no plant-water stress is present, and therefore λE can reach maximal values [22]. This is modeled by the wet-edge representing a maximum value of φ and thus of λE. The true dry-edge is obtained when λE is absent and thus zero, whereas an observed dry-edge represents limited λE values with limited φ min . In order to model increased λE over higher vegetated areas in comparison to areas of sparse vegetation cover, a non-linear scaling of φ is performed with the following equation: where φ i,min represents φ min for each pixel i and RENDVI value, and φ max is according to [22] given by: Based on this equation, the next step spreads the RENDVI-LST triangular space to the minimum and maximum extent of evaporative cooling, represented by φ i,min , LST i,max , and φ i,max , LST i,min for each RENDVI value i respectively [29,30]. The specific φ i value is obtained by: where φ 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, λE was derived for the one full scene consisting of all four flight lines in order to capture a higher land-cover heterogeneity.

λ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-minute block averages of the 10 Hz single measurement intervals. According to [3] and after additional testing, a 30-minute 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. λ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 λE raster file was cropped to the size of the specific footprint of the time and day.

Validation of λE retrievals
For the validation of the UAV based λ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. Figure 3 shows the VIS-NIR orthomosaic in true-color representation for the complete study area and for the region of interest (Figure 4).

Results
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 Figures 5 and 6.   The simultaneously measured λ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 λ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 λ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 λE estimates. The modeled spatially explicit λ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 λ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 λE close to potential evaporation. This can be observed especially in the center of the λ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 λE. As lianas develop more woody biomass [27] and thus consequently reduce the fraction of λE of a canopy surface. This is especially visible in the eastern part of the footprint, where a larger liana infested tree shows reduced λE values.

Discussion
Within this paper, UAV-based estimates of the latent heat flux (λ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, λ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 λ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 λ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 λE [4,13,22,29,32,33,43]. The proposed approach provides, therefore, a very flexible method for the monitoring of spatially-distributed λ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, [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 λ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 λE products across different scales [9]. As shown in Figure 2, the estimated λ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 Figures 2 and 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.

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 λ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 λ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 λ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.