Next Article in Journal
Gap-Filling of a MODIS Normalized Difference Snow Index Product Based on the Similar Pixel Selecting Algorithm: A Case Study on the Qinghai–Tibetan Plateau
Previous Article in Journal
The Delineation and Grading of Actual Crop Production Units in Modern Smallholder Areas Using RS Data and Mask R-CNN
Previous Article in Special Issue
Evaluation of Seven Atmospheric Profiles from Reanalysis and Satellite-Derived Products: Implication for Single-Channel Land Surface Temperature Retrieval
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Land Surface Temperature Retrieval for Agricultural Areas Using a Novel UAV Platform Equipped with a Thermal Infrared and Multispectral Sensor

1
Institute of Bio- and Geosciences, Plant Sciences (IBG-2), Forschungszentrum Jülich GmbH, 52428 Jülich, Germany
2
Department of Geophysics, University of Bonn, 53115 Bonn, Germany
3
German Remote Sensing Data Center (DFD), German Aerospace Center (DLR), 82234 Weßling, Germany
4
Department of Remote Sensing, University of Würzburg, 97074 Würzburg, Germany
5
Center for Remote Sensing of Land Surfaces (ZFL), University of Bonn, 53115 Bonn, Germany
6
Field Lab Campus Klein-Altendorf, University of Bonn, 53359 Rheinbach, Germany
7
Department of Geography, Ruhr-University Bochum, 44801 Bochum, Germany
8
Department of Environmental Remote Sensing & Geoinformatics, University of Trier, 54296 Trier, Germany
*
Author to whom correspondence should be addressed.
Submission received: 29 February 2020 / Revised: 21 March 2020 / Accepted: 24 March 2020 / Published: 27 March 2020

Abstract

:
Land surface temperature (LST) is a fundamental parameter within the system of the Earth’s surface and atmosphere, which can be used to describe the inherent physical processes of energy and water exchange. The need for LST has been increasingly recognised in agriculture, as it affects the growth phases of crops and crop yields. However, challenges in overcoming the large discrepancies between the retrieved LST and ground truth data still exist. Precise LST measurement depends mainly on accurately deriving the surface emissivity, which is very dynamic due to changing states of land cover and plant development. In this study, we present an LST retrieval algorithm for the combined use of multispectral optical and thermal UAV images, which has been optimised for operational applications in agriculture to map the heterogeneous and diverse agricultural crop systems of a research campus in Germany (April 2018). We constrain the emissivity using certain NDVI thresholds to distinguish different land surface types. The algorithm includes atmospheric corrections and environmental thermal emissions to minimise the uncertainties. In the analysis, we emphasise that the omission of crucial meteorological parameters and inaccurately determined emissivities can lead to a considerably underestimated LST; however, if the emissivity is underestimated, the LST can be overestimated. The retrieved LST is validated by reference temperatures from nearby ponds and weather stations. The validation of the thermal measurements indicates a mean absolute error of about 0.5 K. The novelty of the dual sensor system is that it simultaneously captures highly spatially resolved optical and thermal images, in order to construct the precise LST ortho-mosaics required to monitor plant diseases and drought stress and validate airborne and satellite data.

Graphical Abstract

1. Introduction

Land surface temperature (LST) is the primary driving force of the exchange between turbulent heat flux and long-wave infrared (LWIR) radiation (8–15 μ m) at the interface of the Earth’s surface and atmosphere. Therefore, the LST is a key parameter for the physical description of the surface energy and water balance processes at the local to global scale [1,2]. The potential of this crucial parameter has been repeatedly demonstrated in various thermal infrared-based studies and applications, such as evapotranspiration [3,4], hydrological modelling [5], vegetation monitoring [6], ‘urban heat island and urban development’ [7,8,9], climate change and weather conditions [1,10], agriculture [3,11], and the monitoring of land use changes in wetlands [12]. In agricultural applications, precision farming has increasingly required thermal remote sensing techniques to detect water-stressed crops [3,13,14], plant diseases [13,15], and for irrigation management [13,14].
However, LST significantly varies, both temporally and spatially (about 10 K), during the daytime depending on the present weather conditions and the solar irradiation that warms the land surface [1]. Surfaces with high spectral differentiated emissivities, such as bare soil surfaces and rocks (inhomogeneous heating) and sealed areas (‘urban heat islands’) have shown particularly large surface temperature amplitudes in the diurnal cycle [1,16,17]. Water-stressed vegetation also depicts this phenomenon, due to the rapidly changing transpiration rate resulting in a strongly fluctuating plant canopy temperature [18]. Therefore, high repetition rates (depending on water supply and weather conditions) or even continuous temperature measurements are required to account for the detailed and accurate thermal mapping of water-stressed crops or in irrigation management.
Satellite-based LST retrieval methods have been reviewed and compared in several case studies [1,19,20], but depending on the retrieval algorithm, large discrepancies between the retrieved LST and ground truth data can occur. Even the same sensor, such as a high precision satellite instrument (i.e., the Spinning Enhanced Visible and Infrared Imager—SEVIRI—aboard Meteosat-9) using different LST algorithms has featured discrepancies of about 6 K [21]. The reason for this is often an inaccurate or even erroneously estimated surface emissivity [1], which has a great impact on the accuracy of LST [19,20]. The lower the emissivity, the higher the impact of environmental thermal emissions on the LWIR radiation (Section 2.1), which needs to be considered to avoid errors in LST retrieval (Section 2.1.3). Moreover, the results also depend on properly performing atmospheric corrections because of changing weather conditions: A rise in relative humidity leads to decreasing transmittance, which attenuates the thermal signal’s path through the atmosphere [1]. In terms of the great land surface heterogeneity, the central issue is an incorrectly estimated emissivity. Land surface emissivity (LSE) is highly dynamic, particularly in the mid and high latitudes with seasonally changing vegetation cover of miscellaneous landscape types, such as mixed forests, grasslands, wetlands, and in agricultural areas with short-term land use and land cover (LULC) changes [1,10,12,22]. This makes the precise retrieval of LST difficult, as LSE depends on changes in both natural and man-made land surface properties. Hence, an accurate and relatively timely LSE estimation is necessary for the retrieval of correct LST results.
In addition, a scale effect, due to different pixel resolutions when using the same LST algorithm for observations with different spatial scales at different altitudes, must be taken into account. For instance, a pixel might only contain a part of the plant at the leaf level, whereas the same pixel scene might represent a mixture of multiple information of bare soil and overlapping leaves at the canopy level [1,10,17]. Resolution cells of mixed pixels increase the uncertainty if they are not taken into account. With respect to airborne and satellite sensed data, this influencing factor becomes more important, due to the lower spatial resolution at regional scales [1,17]. Therefore, high-contrast and highly spatially resolved thermal infrared mapping is urgently needed (e.g., by airborne campaigns) to investigate the mentioned scale effects and to overcome the large discrepancies observed between remotely sensed data and field measurements.
The use of thermal and multispectral optical sensors on board an unmanned aerial vehicle (UAV) has started with very simple configurations (e.g., [23]). Although some parameters in UAV-based LST retrievals are often simplified due to low flight altitudes (e.g., atmospheric transmittance). Ref. [24] performed a land use classification and assigned specific emissivity values from reviewed literature to each class in order to derive the actual temperature of different target surfaces. Sagan et al. [25] used environmental parameters measured at nearby ground stations to correct air temperature, humidity, and emissivity using linear calibration equations. Similarly, Si et al. [26] determined a single emissivity value that was applied to the entire study area. They measured the emissivity at a certain temperature and integrated it with the spectral response function of the thermal imager [26]. Malbéteau et al. [27] estimated emissivity in terms of vegetation cover, expressed by vegetation index values [17,19]. Some studies, however, do not take emissivity into account and therefore provide results that should be viewed as qualitative rather than quantitative [14,28]. To overcome the challenges mentioned, we have introduced a dual sensor system that was combined with a thermal imager that had a much higher sensor resolution than those available off-the-shelf. A multispectral optical sensor was used to provide vegetation indices to estimate emissivities, whereas, thermal images were captured almost without distortion due to the slow flight speed of the fixed-wing glider.
For LWIR applications, universal emissivity ( ϵ ) for dense and healthy vegetation surfaces of ϵ = 0.985 [29] has typically been assumed, although the emissivities of different landscape types, such as forests, grasslands, and agricultural areas, are diverse [17,30,31]. Additionally, the vegetation density has a direct impact on LSE; for instance, when there is a soil background within a resolution cell. Hence, this assumption is valid only under certain conditions. The emissivity also depends on the state of the plant, and can significantly decrease (e.g., for dry vegetation) [31]. Generally, an error in ϵ of 1% corresponds to a value of 0.75 K [32].
In the thermal spectrum, fully vegetated canopies show Lambertian behaviour, where angular dependence on emissivity is negligible [33]. Accordingly, the thermal anisotropy can be disregarded for high emissivities [33,34]. However, the prerequisite for this is that the radiative transfer model has been applied (Section 2.1). Nevertheless, thermal measurements are affected by water vapour absorption. Therefore, an atmospheric correction is reasonable to consider the water vapour. In particular, at higher relative humidities (>50%), it is expedient to determine tropospheric transmittance, even for a flight altitude (sensor-target distance) of 50 meters [35]. Thus, we performed an atmospheric correction in this study.
In summary, if an incorrect specific emissivity is assumed and environmental thermal emissions (background temperature) and tropospheric transmittance are disregarded, a potential error of about 9 K can arise (Section 3). Eventually, retrieval of the precise LST (or the kinetic temperature) is necessary; for instance, to assess evapotranspiration or to calculate agricultural indices such as the Crop Water Stress Index (CWSI). This is an indicator for the water supply status of plants, which is required in precision farming (e.g., for irrigation management) [36]. An accurate estimation of LSE to generate precise LST maps remains a key prerequisite.
Thermal remote sensing provides complete spatially averaged and highly temporally resolved information about specific radiation characteristics to describe the intricacies of the LST, rather than point values by field measurements over large areas [1]. Hence, we have investigated various algorithms using an LWIR thermal camera, as well as a multispectral sensor covering the visible and near infrared (VNIR) domain, to identify the method with the highest potential in terms of retrieving LST. Multispectral VNIR data were used to derive the Normalised Difference Vegetation Index (NDVI) [37]. We chose predefined NDVI thresholds, where we constrained the emissivity per pixel to distinguish between bare soil surfaces, dense plant canopies, and mixed surface types (Section 2.1.2).
The present research is focused on agricultural areas, which entail highly dynamic LSE values, due to the varying phenological statuses of the crops. Besides the natural ripening stages, this is mainly caused by external impacts, such as water stress, plant diseases, or frost damage. The objective of this study is to retrieve precise LST, taking into account the highly variable LSE derived by using an NDVI threshold approach, as well as the atmospheric impacts. For this purpose, we equipped a UAV with a multispectral sensor and a thermal camera to simultaneously capture VNIR and LWIR imagery of the entire study site. Based on the acquired VNIR and LWIR images, we generated corresponding ortho-mosaics to retrieve and quantify the LST values. Finally, we have developed a standardised LST retrieval algorithm which is optimised for fixed-wing UAV applications, in order to enable semi-supervised thermal mapping in agricultural areas. The study site of the UAV campaigns (April 2018) was an agricultural research field in Rheinbach (North Rhine-Westphalia), which is affiliated with the Agricultural Faculty of Bonn University in Germany.

2. Methodology of the LST Retrieval Algorithm Using UAV-Based Multispectral Data

The following sections describe the UAV-based LST retrieval methodology in detail. The various processing steps required to retrieve precise LST values are illustrated in Figure 1:
For this approach, we developed a dual sensor system to simultaneously acquire multispectral VNIR and LWIR thermal imagery from a flight altitude of 77 meters, resulting in a spatial resolution of 7 cm per pixel (px) and 10 cm/px, respectively (Figure 1i). Different reflectance panels were captured during the flights for subsequent calibrations of the VNIR and LWIR data. The background temperature was incorporated by using a panel of crumpled aluminium foil and quantified with a pyrometer [29], which was sensitive at similar wavelengths as the thermal camera (Appendix A.1). Due to its high specific heat capacity, the water temperature of a nearby pond was used as a reference for validation of the thermal observations (Figure 1ii). The advantage of a water surface reference is its quite stable and homogeneous temperature distribution during the recording time and its constant emissivity (Section 2.2.3). VNIR and LWIR digital numbers (DN) were then converted into reflectance and brightness temperature (BT), respectively. Further preprocessing steps that we performed on the VNIR data included reflectance calibration and band sensitivity normalisation. LWIR images were removed, regarding their contrast applying a 2D Fourier transform (Appendix A.2); afterwards, we conducted drift compensation and non-uniformity correction (Figure 1iii). Subsequent atmospheric correction for the thermal data was based on the radiative transfer model (Figure 1iv). Before constructing the ortho-mosaic, additional calibration of the lens parameters to rectify the LWIR imagery was necessary (Appendix A.2). For both data sets, the same tools (with adapted parameters) to generate the ortho-mosaics were chosen: aligning imagery, then dense cloud and mesh building (Appendix A.2). Then, we calculated the NDVI to constrain the LSE per pixel, based on predefined NDVI thresholds. Co-registration and resampling of the VNIR and BT maps was needed for pixelwise calculation of the stacked layer (Figure 1vi). Finally, we computed the LSE, based on the proposed methods, and retrieved the LST, considering the specific emissivity, transmittance, and background temperature (Figure 1vii).

2.1. Lst Retrieval Methodology

Based on the radiative transfer model (RTM), the radiation path through the troposphere can be written as [17,29]:
L s e n s = ( L s u r f + L r e f l ) · τ + L a t m ,
where L s e n s is the at-sensor radiance, L s u r f the surface radiance attenuated by the atmosphere, and L a t m the upwelling radiance emitted by the atmosphere. L r e f l describes the downwelling radiance emitted by the atmosphere (e.g., by water vapour re-emission), which is reflected at the surface ( ϵ 1 ) and attenuated by the atmosphere. The atmospheric transmittance ( τ ) is dimensionless (defined in the range of 0–1) and for the radiances, the unit W m 2 sr 1 is used. As L s e n s depends on both the wavelength and sensor-specific properties, the inverse Planck function must be used to convert from the raw at-sensor radiance ( L s e n s ) into the at-sensor brightness temperature ( B T s e n s ). Planck’s law (1900) describes the spectral density of electromagnetic radiation emitted by a blackbody at a given temperature [1]. However, the LWIR detector does not cover the entire thermal wavelength range. Thus, the camera manufacturer has used a modified Planck function [38], which is based on the sensor-specific spectral response curve and laboratory-derived Planck constants (Appendix A.1). B T s e n s is automatically calculated by the thermal camera software, with parameters of transmittance and emissivity set to one.
Due to the reflections of natural surfaces ( ϵ < 1 ) in the thermal domain, the radiometric temperature or brightness temperature sensed by the sensor is always lower than the equivalent measured absolute (kinetic) temperature (T). Radiant and kinetic temperatures are two different quantitative terms in thermodynamics [34,39,40], which can be equal in no way, and simply provide the best approximations to describe the temperatures [41]. The thermodynamic relationship for ϵ < 1 can be described, by the Stefan–Boltzmann law (1884), as [34,40,42]:
M = ϵ · σ · T 4 T = M ϵ · σ 4 ,
where M is the irradiance (W m 2 ), ϵ is the emissivity (ranging from 0–1; dimensionless), σ is the Stefan–Boltzmann constant ( 5.67 × 10 8 W m 2 K 4 ), and T is the absolute temperature (K).
According to [29,43], based on RTM (Equation (1)) and using the Stefan–Boltzmann law (Equation (2)), the LST retrieval (Equation (3)) can be written as (Appendix A.3):
L S T : = T s u r f ( ϵ , τ ) = B T s e n s 4 ( 1 ϵ ) · τ · T b k g 4 ( 1 τ ) · T a i r 4 ϵ · τ 4 ,
where T s u r f ( ϵ , τ ) is the retrieved surface temperature, B T s e n s the at-sensor brightness temperature, T b k g the background temperature, and T a i r the air temperature (K, respectively); and emissivity ( ϵ ) and transmittance ( τ ) range from 0–1 (dimensionless). From Equation (Equation (3)), it is obvious that an atmospheric correction (Section 2.1.1) considering τ , T a i r , T b k g , and the corresponding emissivity are needed to retrieve a precise and correct absolute temperature (T) at the surface. At altitudes above 100 meters, we recommend considering the (dry) adiabatic lapse rate to correct the air temperature to the appropriate UAV height. The background temperature was incorporated using a panel of crumpled aluminium foil (Appendix A.1). Equation (3) has been described by [43] and was used in an analogous way in [29], where we included the parameter τ to make the use of an additional commercial software for atmospheric correction unnecessary. The standard software (e.g., FLIR Thermal Studio Pro, FLIR Systems, Inc., Wilsonville, OR, USA; IRBIS 3 professional, InfraTec GmbH, Infrarotsensorik und Messtechnik, Dresden, Germay) uses Equations (4) and (5) (Section 2.1.1) to perform atmospheric correction. To validate the retrieved LST (Section 3), we used the water temperature (as measured by thermocouples) of a nearby pond as relative constant temperature reference (Section 2.2.3). In order to account for small spatial temperature variation, we measured the air temperature at two different locations. Finally, we used the estimated emissivity ( ϵ ) to retrieve the precise LST (Section 2.1.2), as proposed by [17,44].

2.1.1. Impacts of Atmospheric Conditions on Thermal Measurements

In this subsection, we address several complex influencing factors that can affect thermal measurements during data acquisition: Infrared thermography must be executed in thermal equilibrium to obtain physically correct data. However, changing the ambient conditions of the camera (FLIR Tau 2 640), can affect the thermal measurements and reduce the radiometric accuracy by up to ±5 ° C [45], because uncooled thermal sensors (without thermo-electric cooler) cannot maintain a stabilised constant temperature of the microbolometer (focal-plane array, FPA). The accuracy of the measured B T s e n s is dependent on the internal FPA (detector) temperature and housing temperature, as well as external effects caused by the UAV system itself. These parameters fluctuate throughout a flight (e.g., from shaded to unshaded conditions) and should, therefore, be taken into account [46]. The issue of fluctuating detector and housing temperatures has been solved by the ‘Advanced Radiometry’ tool (FLIR) using a built-in shutter to calibrate itself. Before and after changing temperature events, flat field correction (FFC) is performed by submitting a uniform temperature (a flat field) to the detector. During the FFC process, the internal shutter of the camera is closed and quickly reopened to estimate the shutter temperature. Therefore, each FFC event was detected during the data acquisition and tagged in the recordings. Then, the ThermoViewer 2.1.6 software (TeAx Technology GmbH, Wilnsdorf, Germany) was used to perform ‘Drift compensation’ (for temperature fluctuations) with appropriate correction coefficients. Drift of the housing temperature was also compensated by permanent cooling with an adjacently mounted fan. However, sensory effects can also lead to non-uniformity in thermal imagery, such as ‘cold corners’ (vignetting). The same software, hence, was used for a motion-based ‘Non-Uniformity Correction (NUC)’ filter, which distinguished constant errors by shifts within the frames.
Furthermore, thermal measurements are affected by water vapour absorption. Therefore, atmospheric correction in the thermal spectrum was performed to consider tropospheric conditions (Section 1). The required parameters were obtained from the in situ weather station and a microclimatic station in the field (Section 2.2.3). Initially, the present water vapour content (WVC) was computed using parameters such as the air temperature and relative humidity:
ω = ω % · e x p h 1 · T a i r 3 + h 2 · T a i r 2 + h 3 · T a i r + h 4 ,
where ω is the water vapour content (mm), ω % is the relative humidity (ranging from 0–1; dimensionless), T a i r is the air temperature ( ° C), and h 1 , h 2 , h 3 , and h 4 are specific parameters for a temperature range of −40 to +120 ° C in the LWIR domain ( h 1 = 6.8455 × 10 7 , h 2 = −2.7816 × 10 4 , h 3 = 6.939 × 10 2 , and h 4 = 1.5587) [43,47]. Subsequently, the transmittance was estimated using parameters such as the flight altitude (distance) and WVC:
τ = K a t m · e x p d α 1 + β 1 ω + 1 K a t m · e x p d α 2 + β 2 ω ,
where τ is the transmittance (0–1; dimensionless), d is the distance (m), along with specific parameters such as the scaling factor of atmospheric damping ( K a t m = 1.9), the attenuation of the atmosphere without water vapour ( α 1 = 0.0066, α 2 = 0.0126), and the attenuation of water vapour ( β 1 = −0.0023, β 2 = −0.0067) [43,47]. Atmospheric damping is caused by absorption by gaseous components and the turbidity of the troposphere, which decrease the amplitude of the electromagnetic radiation [48]. From this equation, it is obvious that the transmittance depends mainly on the distance and the WVC.

2.1.2. Lse Estimation Based on NDVI Thresholds

Remotely sensed spectral vegetation indices, such as the NDVI, allow a more representative spatial estimation of the emissivity than point measurements [30,37]:
N D V I = N I R R N I R + R ,
where N I R and R are the reflectances in the near-infrared and red spectra, respectively (−1...+1). Moreover, N D V I thresholds provide a reasonable way to differentiate the emissivities of various land surface types. This is particularly necessary for dynamic land surfaces, such as seasonally varying natural vegetation and agricultural areas.
The first predictions of LSE from NDVI values were made in [30,44]. For estimating the emissivity determined by NDVI thresholds, the vegetation cover method was initially proposed in [49,50]. Using information of the fractional vegetation cover or vegetation proportion ( P v ), along with knowledge of the emissivities of full vegetation and bare soil surfaces, this method produces effective LSE. The NDVI threshold method (NDVI-THM), using certain NDVI thresholds, was first introduced in [51], where it was applied to Advanced Very High Resolution Radiometer (AVHRR) satellite data. For applications in agriculture, we identified this method as the most suitable, as it classifies the emissivity into three different surface types: bare soil ( N D V I < N D V I s o i l ), full vegetation ( N D V I > N D V I v e g ), and mixed areas (composed of both attributes). Using information from P v and NDVI histograms to identify the convenient NDVI thresholds is acutely beneficial: No accurate atmospheric correction is essential, because corrected and uncorrected NDVI are similar [17,52]. For known emissivities of full vegetation and bare soil surfaces, P v (Equation (7)) can be obtained, according to [17,52], by:
P v = N D V I N D V I s o i l N D V I v e g N D V I s o i l 2 ,
where N D V I s o i l and N D V I v e g are the NDVI of bare soil and full vegetation, respectively. The P v values can, then, be used to compute the emissivity of mixed pixels ( ϵ ):
ϵ = ϵ v e g · P v + ϵ s o i l · ( 1 P v ) + C ,
where ϵ v e g (=0.988) and ϵ s o i l (=0.935) are the mean emissivities (in the LWIR spectrum) of dense plant canopy [30,31] and bare soil (of loamy to silty texture) [31,53], C is the term of cavity effect due to surface roughness ( C = 4 · d ϵ · P v · ( 1 P v ) ), and d ϵ is a revised parameter (d ϵ ≈ 0.01). For flat areas ( C = ! 0 ), the formula (Equation (8)) simplifies to:
ϵ = ϵ v e g · P v + ϵ s o i l · ( 1 P v ) .
Finally, we constrain the emissivity per pixel based on certain NDVI thresholds to distinguish between three different surface types, as proposed by [17]. The following NDVI thresholds were used to reclassify the LSE ( ϵ ), based on the emissivities of bare soil, full vegetation, and mixed states:
bare soil: N D V I < N D V I s o i l ⇒ NDVI < 0.157 ϵ s o i l = 0.935
dense canopy: N D V I > N D V I v e g ⇒ NDVI > 0.905 ϵ v e g = 0.988
mixed pixels: N D V I s o i l N D V I N D V I v e g ϵ (Equation (9)),
where ϵ v e g is maintained as the global emissivity of dense canopies [30,31] and ϵ s o i l is assumed for loamy to silty (moisture) soil [31,53]. For loamy sand soil, we suggest a lower emissivity of 0.914, as proposed by [30]. A N D V I s o i l value of 0.157 (standard deviation, SD = 0.227) was derived by the means of four samples from the rapeseed field studied. A N D V I v e g value of 0.905 (SD = 0.111) was established by the means of four rape plots (healthy varieties). However, for different crops (i.e., corn, wheat, barley, and rape), we recommend a global N D V I v e g value of 0.814 (as observed in summer 2018); or around 0.8, as proposed by [17]. Nevertheless, the NDVI thresholds for soil and vegetation need to be recalculated before using NDVI-THM to estimate the LSE.
For open canopies, ref. [30] have postulated a logarithmic relationship ( ϵ ) between the mean emissivity and the mean NDVI, which fits best to NDVI values between 0.157 and 0.727 (Equation (10)):
ϵ = a + b · l n ( N D V I ) ,
with a = 1.0010 and b = 0.047 (level of significance: 0.99). To analyse the results of the emissivity approach from information derived from P v (Equation (9)), we used the logarithmic relationship (Equation (10)). Finally, we performed a comparison of both approaches to determine their selectivities, as well as to study their benefits for applications in agriculture (Section 3).

2.1.3. Consequences of Omitting Essential Parameters

Subsequent assumptions were made to emphasise the importance of considering essential parameters, such as transmittance, air temperature, emissivity, and background temperature. In our first assumption, we set the emissivity to one ( ϵ = ! 1), which simplifies the formula (Equation (3)) to:
L S T ( τ ) : = T s u r f ( τ ) = B T s e n s 4 ( 1 τ ) · T a i r 4 τ 4 .
In our second assumption, which presupposes a completely transmissive atmosphere, we set the transmittance to one ( τ = ! 1). This simplifies the same initial formula to:
L S T ( ϵ ) : = T s u r f ( ϵ ) = B T s e n s 4 ( 1 ϵ ) · T b k g 4 ϵ 4 .
Based on these limiting assumptions (Equations (11) and (12)), we made the following two case distinctions: Case i: LST( τ ) includes the atmospheric transmittance, but does not take into account the emissivity or background temperature (Equation (11)); and case ii: LST( ϵ ) includes the emissivity and background temperature, but omits the transmittance (Equation (12)). In summary, we have highlighted the potential implications of omitting these parameters, which should be considered in order to minimise the uncertainties in retrieving precise LST values. The deviations between the uncorrected B T s e n s and LST( τ ) or LST( ϵ ) are shown in the results (Section 3).

2.2. Study Site and Data Sets

2.2.1. Agricultural Research Campus

The study site of the UAV campaigns was the northern part of agricultural research fields of ‘Campus Klein-Altendorf’ (CKA) in Rheinbach (North Rhine-Westphalia), which is affiliated with the Agricultural Faculty of Bonn University, Germany (Figure 2). The total agricultural area of CKA was about 181 ha, mainly cultivated with barley, wheat, corn, and sugar beet; moreover, rapeseed, potatoes and renewable resources (e.g., Miscanthus) were also cultivated. CKA had two weather stations that continuously recorded meteorological data, such as air temperature and humidity, as well as the surface temperature in the northern and southern part of the campus. In the fields, air temperature and humidity were measured continuously by data loggers (Section 2.2.3). This study was focused on the rapeseed field within the ‘Common Experiment’ (Figure 2) during the springtime in 2018. The maintenance of crop rotation and working on homogeneous fields was part of this experiment. The rapeseed field experiment consisted of 63 plots (and border plots), each with dimensions of 3 × 6 meters and planted with 10 different genotypes with a minimal repetition of 6 plots (the present genotypes were Alaska, Pirola, DH5, Markus, Expert, Olympiade, Wotan, Major, Groß Lüsewitzer, and Sensation NZ).

2.2.2. Fixed-Wing Uav System

A fixed-wing UAV (Figure 3) was used to map the heterogeneous and diverse agricultural cropping systems and to efficiently cover the large fields of the campus.
The special features of the UAV are the novel sensor platform (Figure 3) and a slow flight speed over ground (Appendix A.4), which is particularly important for thermal mapping, in order to avoid blurry images. The UAV was mounted with a dual sensor system, consisting of both a multispectral VNIR sensor (Parrot Sequoia; Parrot Drones SAS, Paris, France) and a LWIR thermal camera (FLIR Tau 2 640), attached to a ThermalCapture 2.0 module (TeAx Technology) to simultaneously capture the images (Table 1). The accuracy of the thermal camera ( d T c a m ) under stable ambient conditions was d T c a m = ±0.50 ° C, which has been laboratory-confirmed by [45]. However, in changing ambient conditions, the radiometric accuracy can be reduced by up to ±5.0 ° C [45].
An irradiance-sensing element, belonging to the VNIR sensor, was used to account for the present lighting conditions. This sensor was fixed on the UAV and oriented towards the sky. Both sensors covered the optical spectrum, separated into four bands (centre wavelength): green/G (550 nm), red/R (660 nm), red edge/RE (735 nm), and NIR (790 nm). Moreover, the sensor platform was equipped with an internal fan to cool the devices and to provide a stable ambient temperature, which is relevant for proper thermal imaging. An advanced autopilot (Pixhawk 4) connected to a GPS antenna for positioning and a pitot tube (a device measuring the flight speed with respect to wind speed) on the wing were built-in to automatically fly the programmed routes. Generally, we conducted the UAV campaigns once or twice times per month during the growing season. The flights were split into three parts (approximately 60 ha per section), due to the limited power supply of the UAV batteries, which were generally scheduled at 10:00 (focused study field), 12:00, and 14:00 (CET), to minimise the influence of shadows in the imagery. Every flight took approximately 30–45 minutes.

2.2.3. Meteorological Data and Reference Measurements

A weather station (GWU-Umwelttechnik GmbH) was situated in the northern part of the campus (near the pond), which recorded the following parameters (at 2 m above the ground): wind speed, air temperature ( T a i r ), relative humidity ( ω % ), surface temperature ( T s u r f ) at 0.05 m, and net radiation (Table 2 and Table 3, respectively). Within the rapeseed field, a microclimatic station continuously recorded (at 1 min intervals) the air temperature and relative humidity (HMP110) as shown in Table 3. A pull-up resister (DS18B20) measured the (bare soil) surface temperature with an accuracy of ±0.5 ° C. In addition to the water surface temperature of the pond, T s u r f was used to validate the retrieved LST from the thermal camera (Section 3).
As the weather conditions were quite stable (Table 2) and there were no noticeable fluctuations in T s u r f or abruptly changing plant states (such as crop water stress), no further reference measurements (aside from the water temperature) were carried out with thermocouples to validate the surface temperatures [29]. However, weather conditions have a direct impact on the surface temperature: increasing T a i r will directly increase T s u r f in analogical degree, but also change the net radiation, affecting the surface temperature. The interdependence of T a i r and T s u r f as an approach to surface temperature correction ( T s u r f _ c ) has been proposed by [29]:
T s u r f _ c = T s u r f T a i r + T a i r _ m e a n ,
where T s u r f _ c is the corrected surface temperature, T s u r f , T a i r = < T a i r >, and T a i r _ m e a n is the mean of < T a i r > ( ° C). The correction of T s u r f compensates for the noise caused by drifts in air temperature (Figure 4).
This approach (Equation (13)) was used to account for small discrepancies ( d T s u r f ) in the thermal measurements caused by the impact of < T a i r >. We calculated the differences of T s u r f and T s u r f _ c , which were defined as d T s u r f . There were no significant changes in < T a i r > and T s u r f during the campaign and, so, d T s u r f = 0.28 ± 0.1 ° C was small.
Due to stable weather conditions (Table 2), the averaged values of T b k g (from before, during, and after the flights) were used (Appendix A.1) to consider the background temperature [29]. Moreover, for cloudy conditions (diffuse scattering), the background temperature is quite similar to the measured air temperature (Table 3). Flight 3 showed a much lower background temperature ( T b k g ), due to (‘cold’) clear sky conditions. The transmittance ( τ ) was always about 0.95. Both T b k g and τ were used to perform atmospheric correction.
In addition, during the flights, data at a nearby pond were captured, which were used as reference for the temperature validation (Section 3). The Normalised Difference Water Index (NDWI) was applied to locate the surface of the pond [54,55]:
N D W I = G N I R G + N I R ,
where G and N I R are the reflectances in the green and near-infrared spectra, respectively (−1...+1). According to [56], we used an N D W I threshold of N D W I ≥ 0.3 to detect the water surface and to reclassify the corresponding emissivity ϵ w ( ϵ w = 0.985) [31,34]. Water surfaces almost totally absorb LWIR at a layer thickness of >1 mm and show a constant emissivity [40].
In terms of the water surface temperature ( T p o n d ) measurements, we constructed a floating triangle out of bamboo poles (with a side length of 3 m each) for the thermocouples, which were mounted at each corner and in the middle. The thermocouples were placed at depths of 3.0–30 mm (the dimensions of a thermocouple element is 5.1 × 33 mm). T p o n d was measured directly under the surface layer by four thermocouples, using a data logger with an accuracy of ±0.36 ° C. The mean temperature of the four thermocouples was used to validate the retrieved camera temperature ( L S T p o n d ). According to Equation (3), the captured thermal images of the pond were used to retrieve L S T p o n d . Before and after the flight over the Common Experiment (Figure 2), each set of 24 suitable thermal images of the pond was selected for validation. Both T p o n d and L S T p o n d were used to calculate the temperature deviations ( d T p o n d ) during the flight.

3. Results

The CIR composite (Figure 5a) enhanced the visibility of different varieties for the rape field studied:
Similar reddish colours showed the same varieties, depending on their biomass, and dark green plots were plants damaged by frost due to cold spells in the springtime. These plants belonged to the same genotype and turned out to be less resilient against frost. Plants affected by frost were inhibited in their growth or entirely perished. Striking greenish patches depicted bare soil areas around the plots. These patterns continued in the subsequent maps. The NDVI map (Figure 5b) was generated to find certain NDVI thresholds of bare soil ( N D V I s o i l = 0.157; SD = 0.227) and full vegetation ( N D V I v e g = 0.905; SD = 0.111), which were used to compute the vegetation proportion ( P v ) and to estimate the LSE (Section 2.1.2). The P v map (Figure 5c) was used to compute the LSE (Equation (9)), whereas LSE’ was calculated based on the logarithmic relationship (Equation (10)). High NDVI values represent a healthy and dense canopy; corresponding to P v = 1 for full vegetation. Both the NDVI and the P v maps confirm the lack of biomass for the withered varieties, resulting in low NDVI. The latter map even indicated a marginal vegetation proportion ( P v = 0) for the affected plots, similar to that of bare soil surfaces.
The LSE (Figure 6a) and LSE’ (Figure 6b) maps were computed using certain NDVI thresholds for bare soil, dense canopy, and mixed states (Section 2.1.2). For mixed pixels, LSE was estimated by information from P v (Equation (9)) and LSE’ by using the logarithmic relationship (Equation (10)). Comparison of LSE and LSE’ was performed to study the selectivity of both approaches and to evaluate the accuracy of the NDVI-THM approach (Figure 6c). The minimum emissivity was 0.935, representing loamy to silty soil, and the maximum emissivity was 0.988, representing dense canopy. The LSE (Figure 6a) map, derived by the P v approach, shows detailed information of the open canopy, such as withered plants and plot edges (white areas with an emissivity of about 0.96). The LSE’ map (Figure 6b), based on the logarithmic relationship ( ϵ ’) between the mean emissivity and the mean NDVI, fits for NDVI values between 0.157 and 0.727 better than for NDVI thresholds of up to 0.905 (i.e., for dense-canopied rapeseed).
In contrast, the P v approach could be adjusted to determine the certain NDVI values of different cereals. The comparison (Figure 6c) of both approaches shows a selectivity of dLSE = 0.03, especially for plants affected by frost and plot edges (red areas). This fact could be an indication of a discrepancy in the estimation of lower emissivities (i.e., for dry vegetation). For dense canopy and bare soil, dLSE asymptotically approached zero, which confirms the accuracy of the estimated LSE. However, we have stated that the NDVI-THM approach cannot be used to estimate the emissivity of senescent or withered plants, as the NDVI of dry vegetation (which is considerably affected by water stress) has a similar emissivity spectrum as that of bare soil [57,58]. The emissivity of dry vegetation can reach very low values, compared to healthy plants, which leads to LST values that are considerably underestimated; without additional consideration of meteorological parameters, a potential error of about 9 K can occur (Section 4).
The B T s e n s map (Figure 7a) clearly distinguishes colder areas of the dense canopy and warmer spots (up to 4.8 K) of the rape varieties affected by frost, as well as bare soil surfaces. As the measured at-sensor brightness temperature does not take into account the surface emissivity and current atmospheric conditions, B T s e n s should only be used for surfaces with high emissivities and relative temperature observations under similar weather conditions. The LST retrieval (Figure 7b) included the estimated LSE (Figure 6a), as well as the present atmospheric conditions. The minimum temperature of the dense canopy was 13.5 ° C and the maximum temperature of the bare soil and dried-up plants was 19.2 ° C. The difference between LST and B T s e n s (Figure 7c) showed a maximum deviation (dLST) of 1.0 K for bare soil surfaces and varieties that had suffered from frost. The accuracy of the LST retrieval is strongly dependent on the specific surface emissivity and differs substantially with decreasing emissivity: Without consideration of mixed-resolution cells from soil and vegetation, the LST will be underestimated for affected pixels by at least 1 K (Figure 7c). Hence, we recommend using an NDVI thresholding method (Section 2.1.2) for the emissivity estimation, in order to consider mixed pixels and to minimise the uncertainties.
The final LST retrieval (Figure 7b) depicts almost the same temperature distribution as the B T s e n s map: warmer spots (up to 5.7 K) showed frost-damaged plants and bare soil, and dense canopy occurred as colder areas. Compared to the CIR map, with specific reddish colours for same varieties, the temperature variation within the plots was larger than for different varieties. Therefore, the varieties were not distinguishable by certain LST thresholds.
Moreover, we had nearly constant weather conditions during the campaigns (Section 2.2.3) and no noticeable shifts in the plant state. Thus, the temperature variance of the varieties was generally caused by different evapotranspiration, plant water content, and soil water supply. In our case, the differing canopy temperature within the plots was mainly related to the diverse emissivities of withered and healthy plants.
Finally, we analysed the impacts of the meteorological parameters on thermal measurements. Different assumptions were made to emphasise the effects of when essential parameters, such as transmittance and environmental thermal emissions, were not taken into account in LST retrieval (Section 2.1.3). Each case was compared to the uncorrected B T s e n s (Figure 7a) and the modified LST retrievals (Equations (11) and (12)) to depict the consequences of inaccurately retrieved LSTs: Case i: Atmospheric transmittance included in LST( τ ), but without consideration of the background temperature and specific emissivity (Equation (11)); and case ii: The background temperature and specific emissivity included in LST( ϵ ), but without consideration of transmittance (Equation (12)). For B T s e n s without consideration of the atmospheric transmittance, a temperature deviance of up to 0.4 K occurred (Figure 8a). This was mainly caused by water vapour absorption in the thermal domain. In contrast, omitting the background temperature and specific emissivity led to a deviance of up to 0.6 K (Figure 8b). This result was derived during cloudy sky conditions ( T b k g = 8.8 ° C). In a modified scenario, with respect to the background temperature in clear sky conditions, the temperature deviation was up to 2.9 K (Figure 8c). In all scenarios, the largest deviation was found for bare soil areas. During clear sky conditions, T b k g had minimal temperature −40 ° C and maximal temperature 10 ° C, in the case of a diffusely scattering sky [32]. Diffuse scattering (in overcast conditions) leads to almost identical measurements between the air temperature and the observed background temperature ( T b k g ). However, the background temperature should be taken into account to minimise uncertainties and to prevent the LST from being underestimated. In conclusion, environmental thermal emission (if ϵ < 1) is the most influential meteorological parameter in thermal measurements.
The NDWI map (Figure 9a) was generated to detect the pond’s surface and to reclassify the corresponding emissivity ( ϵ w ) of the water using a determined NDWI threshold (Section 2.2.3). Based on the emissivity map (Figure 9b), the mean pond surface temperature ( L S T p o n d ) was retrieved (Figure 9c). Before and after the flight over the agricultural experiment site (flights 1 and 3, respectively), we retrieved 24 suitable LST images ( L S T p o n d ) to validate the thermal data using the measurements with the pond thermocouples ( T p o n d ). Both L S T p o n d and T p o n d were used to calculate the mean temperature deviations ( d T p o n d ) to assess the accuracy of retrieved LST.
The mean absolute error (MAE) was used to compute the errors of the temperature comparisons [59]. Flight 1 (Figure 10a) had an error of d T p o n d _ 1 = 0.53 K, with a coefficient of determination of R 2 = 0.985 (n = 24). For flight 3 (Figure 10b), the MAE was d T p o n d _ 2 = 0.54 K and R 2 = 0.964 (n = 24).
The results of the performed validation, especially for flight 1, confirmed the high significance of the performance of the retrieval algorithm (Equation (3)). During flight 3, the correlation was somewhat less significant in comparison to flight 1.
The temperature validation of the rapeseed experiment ( d T r a p e ) also confirmed the high performance of the LST retrieval (Figure 11). The MAE was d T r a p e = 0.41 K and R 2 = 0.879 (n = 27).
For the UAV campaigns using the proposed LST retrieval algorithm (Equation (3)), we found the accuracy of the thermal measurements to be around MAE = 0.5 K.

4. Discussion

In this study, we have developed an UAV-based LST retrieval algorithm (Equation (3)) suitable for conducting thermal campaigns, especially in agricultural areas. Agricultural applications increasingly demand precise and current LST maps for irrigation scheduling, evapotranspiration and drought stress monitoring, and plant disease detection. A true physical temperature, as calculated by the LST retrieval of crops, is essential for calculating indices such as the Vegetation Health Index (VHI) or CWSI [13]. The accuracy of retrieved LST depends mainly on properly estimating the emissivity of the surface type (LSE). Thus, we have used the proposed NDVI threshold method (Section 2.1.2) to distinguish the emissivities in three classes: bare soil ( N D V I < N D V I s o i l ), dense canopy ( N D V I > N D V I v e g ), and mixed states (Equation (9)). For dense canopy, we suggest a global N D V I v e g of 0.814 (as observed in summer 2018) or around 0.8 (as proposed by [17]), which was averaged for different crops (i.e., corn, wheat, barley, and rape). For bare soil, we derived N D V I s o i l to be 0.157, in agreement with [30], or around 0.2 (as proposed by [17]). In the case of higher sand content, lower NDVI values arise.
The determined emissivity of dense canopy ( ϵ v e g = 0.988) agreed with frequently cited values, as in [27,29,31,60]. However, the plants must be healthy, as the emissivity decreases significantly in the case of withered crops. In the study site, soil surfaces predominantly consisted of loamy to silty textures. Thus, we maintained an emissivity of bare soil of ϵ s o i l = 0.935 (as proposed by [31,53]). For a sandy soil texture, we suggest a lower emissivity, such as 0.914 (as proposed by [30]). Regardless of the emissivity of mixed resolution cells, the LST was underestimated for the affected pixels by at least 1 K (Figure 7c). The lower the emissivity, the higher the impact of the background temperature ( T b k g ) on the retrieved LST. Hence, we recommend atmospheric correction, especially in (‘cold’) clear sky conditions, in which observation periods commonly occur to ensure the best radiometric image quality. However, this method reaches its limits in the case of lower emissivities (e.g., ϵ = 0.88 ), such as those for dry vegetation [31,34]. Disregarding essential meteorological parameters (Section 2.1.3) and assuming an overall emissivity for dense canopy (for instance, of ϵ v e g = 0.985), a potential error of about 9 K occurs.
For data validation, we continuously measured the water temperature of a nearby pond by thermocouples, which we compared with each set of 24 suitable retrieved LST images of the pond from the flights before and after observing the experiment. Water bodies provide a convenient opportunity for reference measurements, due to their relatively homogeneous and constant surface temperature. For the pond, we assumed an emissivity of ϵ w = 0.985 (according to [31,34]). The mean emissivity (in the LWIR range of the camera) of distilled water ( ϵ = 0.986), obtained from the emissivity library of the Moderate Resolution Imaging Spectrometer (MODIS) confirms this assumption [61].
The advantages and disadvantages of conducting thermography under homogeneously cloudy or clear sky conditions remain to be discussed: Clear sky conditions provide the most contrasting results; and cloudless weather is indispensable for the validation of airborne or satellite data with UAV observations. However, the lower the emissivity, the higher the reflection of the surface. Therefore, the impact of environmental thermal emissions on bare soils and dry vegetation is the greatest, which must be taken into consideration. Cloudy conditions have a strong impact on the quality of the captured thermal imagery, since they reduce the contrast and radiometric resolution [13]. However, the impact of environmental thermal emissions is low (isotropic radiation) and the background temperature is quite similar to the air temperature. Scenarios with the biggest drawbacks are rapidly changing air and surface temperatures or cloud shadows in the recordings. In order to ensure the best lighting conditions, we recommend thermal measurements in clear sky or, at least, in homogeneously cloudy conditions.
In conclusion, the validations showed MAEs of d T p o n d _ 1 = 0.53 K ( R 2 = 0.985) and d T p o n d _ 2 = 0.54 K ( R 2 = 0.964), respectively. These results confirm the highly significant performance of the retrieval algorithm (Equation (3)). For the thermal measurements, we stated an MAE of about 0.5 K.

5. Conclusions

In this study, we present a semi-supervised algorithm that combines UAV-based thermal infrared and multispectral measurements in order to obtain precise LST maps for agricultural applications. For this purpose, we have developed an innovative dual sensor platform mounted on a fixed-wing glider that simultaneously captures VNIR and LWIR images. In the analyses, we have noted the improvements that can be achieved by performing an atmospheric correction to minimise the uncertainties. The processing chain incorporates:
  • Construction of an ortho-mosaic for NDVI maps, LST maps, and emissivity maps;
  • Atmospheric correction, including the background temperature; and
  • Emissivity estimation for three surface types (bare soil, dense canopy, and mixed states).
The high spatial resolution of the LST maps and the capacity of the UAV system to observe the entire study site (181 ha) within 3–6 hours allow for operational agricultural applications and the cost-effective validation of airborne and satellite data complying with high quality standards. Special advantages of the proposed method for agricultural applications are:
  • Low-cost UAV campaigns for the thermal mapping of large areas with high repetition rates;
  • Irrigation management, crop water stress monitoring, and plant disease detection; and
  • High spatial resolution of the LST maps (10 cm/pixel) for thermography at the leaf level.
This approach requires neither LST ground stations (except for validation procedures) nor in-situ stations to account for canopy emissivity change over time. However, the proposed NDVI threshold method needs NDVI means of the current canopy to estimate the emissivities. In order to ensure high quality LST results, we recommend performing atmospheric corrections and considering the background temperature. Further investigation is needed to find a method for the emissivity estimation of senescent and withered vegetation, in order to prevent the LST from being underestimated in the case of lower emissivities.

Author Contributions

Conceptualization, S.H., B.S. and J.M.; methodology, S.H.; software, S.H. and B.S.; validation, S.H.; formal analysis, S.H., B.S., J.S. and N.W.; investigation, S.H.; resources, U.R., T.K. and O.M.; data curation, S.H. and B.S.; writing–original draft preparation, S.H.; writing–review and editing, S.H., B.S. and F.T.; visualization, S.H.; supervision, C.J., Uwe Rascher, B.S., T.U. and A.K.; project administration, S.H.; funding acquisition, U.R. and S.H. All authors have read and agreed to the published version of the manuscript.

Funding

This project was financed by a scholarship of the ’Evangelisches Studienwerk e.V. Villigst’, supported by the German Federal Ministry of Education and Research (BMBF). Part of this study was supported by the German Plant Phenotyping Network (DPPN), which is funded by the BMBF (project identification number: 031A053).

Acknowledgments

Additionally, we would like to thank Kai Wittneben for conducting the flights and Simon König, Lars Grünhagen, Matthias Langen, and Ines Munoz-Fernandez for their assistance in the field during data acquisition and data pre-processing, as well as Valerie Graw, Stefanie Steinbach, and Vivien Czepok for the English revision.

Conflicts of Interest

The contents of this article are the sole responsibility of the authors. The authors declare no conflict of interest.

Appendix A

Appendix A.1. Pre-Processing of VNIR and LWIR Data

Immediately before take-off and after landing, we captured three images (in all VNIR bands) of a calibrated reflectance panel (MicaSense; Inc., Seattle, WA, USA). Four additional greyscale tarps (sailcloth, calibrated by a field spectrometer) were captured during the flights. All images contained GPS, roll, pitch, and yaw angle information as metadata (exchangeable image file format, Exif) obtained from the VNIR sensor. The stored GPS and solid angle information were used by the ortho-mosaic software to correct the position of the imagery. The total captured VNIR imagery (including all four bands) comprised about 7300 recordings (resulting from a frame rate of 2/3 Hz) with a volume of about 17 GB, stored as a 16-bit (binary raw data, raw) tagged image file format (TIFF). The total captured LWIR imagery comprised about 24,000 recordings (with a frame rate of 8.33 Hz) with a volume of about 5 GB, stored as 14-bit raw radiometric data stream (TMC). The ThermoViewer 2.1.6 software was used to convert the TMC files into the 16-bit TIFF format, including metadata with GPS and heading information from the autopilot as Exif data. The amount of captured imagery was necessary to generate ortho-mosaics and to compensate for the loss of discarded images (i.e., those taken during the launch/landing maneuvers and blurred recordings). The VNIR images of the calibrated reflectance panels were incorporated for radiometric calibration. The PhotoScan 1.4.2 software (Agisoft PhotoScan, St. Petersburg, Russia) was used for reflectance calibration and for normalisation of the band sensitivity.
For thermal sensors, radiometric calibration is essential to correctly convert raw DN into (spectral) radiance. The DN of the thermal camera (FLIR Tau 2 640) consists of specific instrument factors, such as the self-emission of the detector and its spectral response curve (FLIR’s Tau 2 spectral response curve [62]). Calibrated blackbody sources at exactly known temperatures allow us to derive the spectral response curve [63]. This camera is the most sensitive at around 9.75 μ m. Blackbody calibration, however, was done by the manufacturer FLIR. The ThermoViewer 2.1.6 software was used for NUC and drift compensation of the LWIR imagery (Section 2.1.1). Thermal reflections at the land surface by emitting background objects (e.g., buildings) and re-emission—for instance, by lower clouds—were determined to correct the LWIR data (Equation (3)). A pyrometer (optris MSPro; Optris GmbH, Berlin, Germany), with emissivity set to ϵ = 1, was used to quantify the background temperature ( T b k g ) on a panel of crumpled aluminium foil [29,32,64], which had a very low emissivity of ϵ = 0.03 (0.04 at 10 μ m [53]) and acted as an isotropically reflecting Lambert radiator (see Lambert’s cosine law, 1760). The pyrometer indirectly captured the background temperature in the range of 8–14 μ m (accuracy: ±1 ° C). As a pre-check, we directly measured the overhead sky (not including the sun) to compare with the T b k g value obtained from the aluminium panel [65]. Averaged temperatures of T b k g from before, during, and after the flights were used for atmospheric correction.

Appendix A.2. Ortho-Mosaic Construction and Used Tools and Parameters

Overlaps of the captured imagery were standardised for the smallest field of view (FOV), therefore based on the LWIR camera. Generally, overlaps result from the FOV and flight altitude, which was about 77 meters. In addition, forward overlaps depend on the frame rate and the UAV’s ground speed. The resulting distance between each flight line was about 25.5 meters, which was relevant for the side overlaps. Based on the ground speed of 5–8 ms 1 , the FOV of the LWIR camera lens (see Table 1), and the frame rate (LWIR: 8.33 Hz), the forward overlap was about 98.8% and side overlap was about 60%. More detailed information is listed in Table A1.
Table A1. Details on the forward and side overlap in the captured imagery at the corresponding altitude of 77 meters.
Table A1. Details on the forward and side overlap in the captured imagery at the corresponding altitude of 77 meters.
SensorGround Image SizeFrame RateGround SpeedFrontlapSidelap
LWIR/Tau 263.79 × 51.53 m8.33 Hz5 ms 1 98.8%60.0%
VNIR/Sequoia92.35 × 69.37 m2/3 Hz5 ms 1 89.2%72.3%
The Calibrate Reflectance tool (PhotoScan) with selected parameters of reflectance panels and sun sensor was executed, taking into account the captured VNIR images of the calibrated reflectance panels. The required reflectance values of each band were either obtained from the company of the used panel or from field spectrometer measurements of the four greyscale tarps described above (Section 2.2.3). The Normalize band sensitivity tool was utilised to include the spectral response (for each band) of the VNIR sensor (contained as metadata). For the VNIR sensor, the following pre-calibrated lens parameters were used: focal length (f, pixels), co-ordinates of the main point (cx and cy), and radial distortion coefficients (k1–k4). For the LWIR camera, the following most fitting geometric lens parameters, which were generated during self-calibration in the alignment process (of the Adaptive camera model fitting tool), were used as fixed pre-calibrated values: f, cx and cy, k1–k4, biased transformation coefficients (b1 and b2), and tangential distortion coefficients (p1–p4). After the pre-processing steps, the ortho-mosaic was constructed, generating both the VNIR and LWIR maps. Ortho-mosaic construction is a common method in photogrammetry. Various structure from motion (SFM) software, such as PhotoScan, Pix4D (Pix4D, Lausanne, Switzerland), and VisualSFM (a freeware GUI application) provide tools to generate digital surface models (DSM) and ortho-images while compensating for undesirable lens distortion [46]. Even when using a high-performance computer with two graphics cards (for a total of 24 GB of GDDR5 memory), mosaic construction is a time-consuming process due to the quantity of high-resolution imagery and the required processing operations (the appropriate parameters are listed in Table A2). Ortho-mosaic construction of the entire study site, using the software PhotoScan 1.4.2, took us about three days. The remaining calculations (Figure 1) were conducted using the R language. In this context, before generating the BT ortho-mosaic, a method was applied to exclude blurry images from further analysis. As the frame rate of the thermal image data acquisition was quite high, the sharpest image from each sequence of five consecutive recorded images was identified. The identification of the sharpest image was done by calculating an image quality measure from the representation of every single image in the frequency domain after applying a 2D Fourier transform. A detailed description of this method can be found in [66]. The image with the highest image quality measure within each sequence was selected. Thus, only the sharpest images were used in the subsequent process of image mosaic construction.
Table A2. Tools and parameters selected for ortho-mosaic construction.
Table A2. Tools and parameters selected for ortho-mosaic construction.
Align PhotosVNIR ImageryLWIR Imagery
AccuracyHighestHighest
Generic preselectionEnabledEnabled
Reference preselectionEnabledEnabled
Key point limit50,00050,000
Tie point limit40000
Adaptive camera model fittingEnabledEnabled
Build Dense Cloud
QualityUltra HighHigh
Depth filteringModerateMild
Build Mesh
Surface typeHeight field (2.5D)Height field
Source dataDense cloudDense cloud
Face countHighHigh
InterpolationDisabledDisabled
Build Orthomosaic
SurfaceMeshMesh
Blending modeMosaicMosaic
Hole fillingDisabledEnabled
Back-face cullingDisabledDisabled
Export Orthomosaic
Raster transformIndex valueIndex value
TIFF compressionNoneNone
Write BigTIFF fileEnabledEnabled
Generate TIFF overviewEnabledEnabled

Appendix A.3. Derivation of the LST Retrieval Algorithm (Equation (3))

From the radiative transfer model (RTM) to the at-sensor temperature ( T s e n s ) [29,43]:
R T M : = L s e n s = ϵ · τ · L s u r f + ( 1 ϵ ) · τ · L r e f l + ( 1 τ ) · L a t m ,
where term 1 is the emitted radiance from the surface ( ϵ · L s u r f ), term 2 the downwelling radiance ( L r e f l ) reflected at the surface (1 − ϵ ), and term 3 the emitted radiance from the atmosphere ( L a t m ); the radiances (W m 2 sr 1 ) are attenuated by the atmosphere ( τ ). The emissivity ( ϵ ) and atmospheric transmittance ( τ ) are dimensionless (0–1). Simplification of Equation (A1) by setting ϵ = 1 (term 1), (1 − ϵ ) = 1 (term 2), and (1 − τ ) = 1 (term 3):
L s e n s = ( L s u r f + L r e f l ) · τ + L a t m .
The Stefan–Boltzmann law:
L = ϵ · σ · T 4 T = L ϵ · σ 4 ,
where L is the radiance (W m 2 sr 1 ), ϵ is the emissivity (0–1), σ is the Stefan–Boltzmann constant ( 5.67 × 10 8 W m 2 K 4 ), and T is the absolute temperature (K). By setting the steradian (sr) to one, Equation (A3) is solved for irradiance (M) Equation (2).
From the at-sensor temperature ( T s e n s ) to the surface temperature ( T s u r f ):
T s e n s 4 · ϵ · σ = ϵ · σ · τ · T s u r f 4 + ( 1 ϵ ) · σ · τ · T b k g 4 + ( 1 τ ) · σ · T a i r 4 ,
where term 1 consists of the surface temperature ( T s u r f ), term 2 the background temperature ( T b k g ) reflected at the surface ((1 − ϵ σ · τ ), and term 3 the air temperature ( T a i r ). At altitudes above 100 meters, the adiabatic lapse rate should be taken into account to derive the air temperature of the appropriate height. The background temperature can be indirectly measured from the thermal reflection at a panel of bottled aluminium foil using a pyrometer (Appendix A.1). With ϵ = 1, T s e n s 4 · ϵ · σ B T s e n s 4 · σ , solving the formula (Equation (A4)) for T s u r f :
T s u r f 4 · ϵ · σ · τ = σ · B T s e n s 4 ( 1 ϵ ) · σ · τ · T b k g 4 ( 1 τ ) · σ · T a i r 4
T s u r f 4 = σ · B T s e n s 4 ϵ · σ · τ ( 1 ϵ ) · σ · τ · T b k g 4 ϵ · σ · τ ( 1 τ ) · σ · T a i r 4 ϵ · σ · τ
T s u r f ( ϵ , τ ) = B T s e n s 4 ( 1 ϵ ) · τ · T b k g 4 ( 1 τ ) · T a i r 4 ϵ · τ 4 ,
where T s u r f ( ϵ , τ ) is the retrieved surface temperature, B T s e n s the at-sensor brightness temperature, T b k g the background temperature, and T a i r the air temperature (K, respectively); and emissivity ( ϵ ) and transmittance ( τ ) range from 0–1 (dimensionless).

Appendix A.4. Detailed Technical Information About the UAV (Table A3)

Table A3. Technical details of the fixed-wing UAV.
Table A3. Technical details of the fixed-wing UAV.
ItemUnit
span4 m
length2 m
take-off weight5 kg
payload700 g
ground speed5–8 ms 1
slowest airspeeddown to 4.5 ms 1
general airspeed7–10 ms 1
flight timeup to 60 min
flight distanceup to 20–25 km
battery (LiPo 3S)11.1 V; 10,000 mAh (UAV)
battery (LiPo 2S)7.4 V; 3000 mAh (sensor system)
current in steady flight10 A
current at full throttle70 A

References

  1. Li, Z.L.; Tang, B.H.; Wu, H.; Ren, H.; Yan, G.; Wan, Z.; Trigo, I.F.; Sobrino, J.A. Satellite-derived land surface temperature: Current status and perspectives. Remote Sens. Environ. 2013, 131, 14–37. [Google Scholar] [CrossRef] [Green Version]
  2. Sobrino, J.A.; Del Frate, F.; Drusch, M.; Jiménez-Muñoz, J.C.; Manunta, P.; Regan, A. Review of thermal infrared applications and requirements for future high-resolution sensors. IEEE Trans. Geosci. Remote Sens. 2016, 54, 2963–2972. [Google Scholar] [CrossRef]
  3. Gerhards, M.; Schlerf, M.; Mallick, K.; Udelhoven, T. Challenges and Future Perspectives of Multi-/Hyperspectral Thermal Infrared Remote Sensing for Crop Water-Stress Detection: A Review. Remote Sens. 2019, 11, 1240. [Google Scholar] [CrossRef] [Green Version]
  4. Raoufi, R.; Beighley, E. Estimating Daily Global Evapotranspiration Using Penman–Monteith Equation and Remotely Sensed Land Surface Temperature. Remote Sens. 2017, 9, 1138. [Google Scholar] [CrossRef] [Green Version]
  5. Su, Z. The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes. Hydrol. Earth Syst. Sci. 2002, 6, 85–100. [Google Scholar] [CrossRef]
  6. Berni, J.A.; Zarco-Tejada, P.J.; Suárez, L.; Fereres, E. Thermal and narrowband multispectral remote sensing for vegetation monitoring from an unmanned aerial vehicle. IEEE Trans. Geosci. Remote Sens. 2009, 47, 722–738. [Google Scholar] [CrossRef] [Green Version]
  7. Deilami, K.; Kamruzzaman, M.; Liu, Y. Urban heat island effect: A systematic review of spatio-temporal factors, data, methods, and mitigation measures. Int. J. Appl. Earth Obs. Geoinf. 2018, 67, 30–42. [Google Scholar] [CrossRef]
  8. He, B.J.; Zhao, Z.Q.; Shen, L.D.; Wang, H.B.; Li, L.G. An approach to examining performances of cool/hot sources in mitigating/enhancing land surface temperature under different temperature backgrounds based on landsat 8 image. Sustain. Cities Soc. 2019, 44, 416–427. [Google Scholar] [CrossRef]
  9. Li, H.; Zhou, Y.; Li, X.; Meng, L.; Wang, X.; Wu, S.; Sodoudi, S. A new method to quantify surface urban heat island intensity. Sci. Total Environ. 2018, 624, 262–272. [Google Scholar] [CrossRef]
  10. Merchant, C.J.; Matthiesen, S.; Rayner, N.A.; Remedios, J.J.; Jones, P.D.; Olesen, F.; Trewin, B.; Thorne, P.; Auchmann, R.; Corlett, G.K.; et al. The surface temperatures of Earth: Steps towards integrated understanding of variability and change. Geosci. Instrum. Methods Data Syst. 2013, 2, 305–321. [Google Scholar] [CrossRef] [Green Version]
  11. Sousa, D.; Small, C. Mapping and Monitoring Rice Agriculture with Multisensor Temporal Mixture Models. Remote Sens. 2019, 11, 181. [Google Scholar] [CrossRef] [Green Version]
  12. Muro, J.; Strauch, A.; Heinemann, S.; Steinbach, S.; Thonfeld, F.; Waske, B.; Diekkrüger, B. Land surface temperature trends as indicator of land use changes in wetlands. Int. J. Appl. Earth Obs. Geoinf. 2018, 70, 62–71. [Google Scholar] [CrossRef]
  13. Khanal, S.; Fulton, J.; Shearer, S. An overview of current and potential applications of thermal remote sensing in precision agriculture. Comput. Electron. Agric. 2017, 139, 22–32. [Google Scholar] [CrossRef]
  14. Tucci, G.; Parisi, E.I.; Castelli, G.; Errico, A.; Corongiu, M.; Sona, G.; Viviani, E.; Bresci, E.; Preti, F. Multi-sensor UAV application for thermal analysis on a dry-stone terraced vineyard in rural Tuscany landscape. ISPRS Int. J. Geo-Inf. 2019, 8, 87. [Google Scholar] [CrossRef] [Green Version]
  15. Mogili, U.; Deepak, B. Review on Application of Drone Systems in Precision Agriculture. Procedia Comput. Sci. 2018, 133, 502–509. [Google Scholar] [CrossRef]
  16. Gillespie, A.; Rokugawa, S.; Matsunaga, T.; Cothern, J.; Hook, S.; Kahle, A.A.B.; Steven Cothern, J.; Hook, S.; Kahle, A.A.B. A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images. IEEE Trans. Geosci. Remote Sens. 1998, 36, 1113–1126. [Google Scholar] [CrossRef]
  17. Sobrino, J.A.; Jimenez-Muoz, J.C.; Soria, G.; Romaguera, M.; Guanter, L.; Moreno, J.; Plaza, A.; Martinez, P. Land Surface Emissivity Retrieval From Different VNIR and TIR Sensors. IEEE Trans. Geosci. Remote Sens. 2008, 2, 316–327. [Google Scholar] [CrossRef]
  18. Leuzinger, S.; Körner, C. Tree species diversity affects canopy leaf temperatures in a mature temperate forest. Agric. For. Meteorol. 2007, 146, 29–37. [Google Scholar] [CrossRef]
  19. Li, Z.L.; Wu, H.; Wang, N.; Qiu, S.; Sobrino, J.A.; Wan, Z.; Tang, B.H.; Yan, G. Land surface emissivity retrieval from satellite data. Int. J. Remote Sens. 2013, 34, 3084–3127. [Google Scholar] [CrossRef]
  20. Sekertekin, A.; Bonafoni, S. Land Surface Temperature Retrieval from Landsat 5, 7, and 8 over Rural Areas: Assessment of Different Retrieval Algorithms and Emissivity Models and Toolbox Implementation. Remote Sens. 2020, 12, 294. [Google Scholar] [CrossRef] [Green Version]
  21. Gao, C.; Jiang, X.; Li, Z.L.; Nerry, F. Comparison of the Thermal Sensors of SEVIRI and MODIS for LST Mapping. In Thermal Infrared Remote Sensing; Remote Sensing and Digital Image Processing; Springer: Dordrecht, The Netherlands, 2013; pp. 233–252. [Google Scholar] [CrossRef]
  22. Leuzinger, S.; Vogt, R.; Körner, C. Tree surface temperature in an urban environment. Agric. For. Meteorol. 2010, 150, 56–62. [Google Scholar] [CrossRef]
  23. Bendig, J.; Bolten, A.; Bareth, G. Introducing a low-cost mini-uav for thermal- and multispectral-imaging. ISPRS Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2012, XXXIX-B1, 345–349. [Google Scholar] [CrossRef] [Green Version]
  24. Naughton, J.; McDonald, W. Evaluating the Variability of Urban Land Surface Temperatures Using Drone Observations. Remote Sens. 2019, 11, 1722. [Google Scholar] [CrossRef] [Green Version]
  25. Sagan, V.; Maimaitijiang, M.; Sidike, P.; Eblimit, K.; Peterson, K.T.; Hartling, S.; Esposito, F.; Khanal, K.; Newcomb, M.; Pauli, D.; et al. UAV-based high resolution thermal imaging for vegetation monitoring, and plant phenotyping using ICI 8640 P, FLIR Vue Pro R 640, and thermomap cameras. Remote Sens. 2019, 11, 330. [Google Scholar] [CrossRef] [Green Version]
  26. Si, M.; Tang, B.H.; Li, Z.L. Estimation of land surface temperature from unmanned aerial vehicle loaded thermal imager data. Int. Geosci. Remote Sens. Symp. IGARSS 2018, 2018, 1210–1213. [Google Scholar] [CrossRef]
  27. Malbéteau, Y.; Parkes, S.; Aragon, B.; Rosas, J.; McCabe, M.F. Capturing the diurnal cycle of land surface temperature using an unmanned aerial vehicle. Remote Sens. 2018, 10, 1407. [Google Scholar] [CrossRef] [Green Version]
  28. Raeva, P.L.; Šedina, J.; Dlesk, A. Monitoring of crop fields using multispectral and thermal imagery from UAV. Eur. J. Remote Sens. 2019, 52, 192–201. [Google Scholar] [CrossRef] [Green Version]
  29. Maes, W.H.; Huete, A.R.; Steppe, K. Optimizing the Processing of UAV-Based Thermal Imagery. Remote Sens. 2017, 9, 476. [Google Scholar] [CrossRef] [Green Version]
  30. VandeGriend, A.A.; Owe, M. On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces. Int. J. Remote Sens. 1993, 14, 1119–1131. [Google Scholar] [CrossRef]
  31. Lillesand, T.; Kiefer, R.W.; Chipman, J. Remote Sensing and Image Interpretation; John Wiley & Sons: Hoboken, NJ, USA, 2014. [Google Scholar]
  32. Jones, H.; Archer, N.; Rotenberg, E.; Casa, R. Radiation measurement for plant ecophysiology. J. Exp. Bot. 2003, 54, 879–889. [Google Scholar] [CrossRef]
  33. Sobrino, J.; Jimenezmunoz, J.; Verhoef, W. Canopy directional emissivity: Comparison between models. Remote Sens. Environ. 2005, 99, 304–314. [Google Scholar] [CrossRef]
  34. Kuenzer, C.; Dech, S. Theoretical Background of Thermal Infrared Remote Sensing. In Thermal Infrared Remote Sensing; Remote Sensing and Digital Image Processing; Springer: Dordrecht, The Netherlands, 2013; pp. 1–26. [Google Scholar] [CrossRef]
  35. Minkina, W.; Klecha, D. Atmospheric transmission coefficient modelling in the infrared for thermovision measurements. J. Sens. Sens. Syst. 2016, 5, 17–23. [Google Scholar] [CrossRef] [Green Version]
  36. Sepaskhah, A.R.; Kashefipour, S.M. Relationships between leaf water potential, CWSI, yield and fruit quality of sweet lime under drip irrigation. Agric. Water Manag. 1994, 25, 13–21. [Google Scholar] [CrossRef]
  37. Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 1979, 8, 127–150. [Google Scholar] [CrossRef] [Green Version]
  38. Pour, T.; Miřijovský, J.; Purket, T. Airborne thermal remote sensing: The case of the city of Olomouc, Czech Republic. Eur. J. Remote Sens. 2019, 52, 209–218. [Google Scholar] [CrossRef]
  39. Bott, A. Synoptische Meteorologie: Methoden der Wetteranalyse und-Prognose; Springer: Berlin/Heidelberg, Germany, 2012; p. 486. [Google Scholar]
  40. Demtröder, W. Experimentalphysik 2: Elektrizität und Optik, 7th ed.; Springer: Berlin/Heidelberg, Germany, 2017. [Google Scholar]
  41. Norman, J.M.; Becker, F. Terminology in thermal infrared remote sensing of natural surfaces. Agric. For. Meteorol. 1995, 77, 153–166. [Google Scholar] [CrossRef]
  42. Wang, K.; Wan, Z.; Wang, P.; Sparrow, M.; Liu, J.; Zhou, X.; Haginoya, S. Estimation of surface long wave radiation and broadband emissivity using Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature/emissivity products. J. Geophys. Res. Atmos. 2005, 110, D11109. [Google Scholar] [CrossRef]
  43. Tran, Q.H.; Han, D.; Kang, C.; Haldar, A.; Huh, J. Effects of Ambient Temperature and Relative Humidity on Subsurface Defect Detection in Concrete Structures by Active Thermal Imaging. Sensors 2017, 17, 1718. [Google Scholar] [CrossRef]
  44. Valor, E.; Caselles, V. Mapping land surface emissivity from NDVI: Application to European, African, and South American areas. Remote Sens. Environ. 1996, 57, 167–184. [Google Scholar] [CrossRef]
  45. Kelly, J.; Kljun, N.; Olsson, P.O.; Mihai, L.; Liljeblad, B.; Weslien, P.; Klemedtsson, L.; Eklundh, L. Challenges and best practices for deriving temperature data from an uncalibrated UAV thermal infrared camera. Remote Sens. 2019, 11, 567. [Google Scholar] [CrossRef] [Green Version]
  46. Ribeiro-Gomes, K.; Hernández-López, D.; Ortega, J.F.; Ballesteros, R.; Poblete, T.; Moreno, M.A. Uncooled Thermal Camera Calibration and Optimization of the Photogrammetry Process for UAV Applications in Agriculture. Sensors 2017, 17, 2173. [Google Scholar] [CrossRef] [PubMed]
  47. Minkina, W.; Klecha, D. 1.4-Modeling of Atmospheric Transmission Coefficient in Infrared for Thermovision Measurements. In Proceedings of the IRS2 2015, Nürnberg, Germany, 19 May–21 May 2015; pp. 903–907. [Google Scholar] [CrossRef]
  48. Hirschel, E.H.; Prem, H.; Madelung, G. Aeronautical Research in Germany: From Lilienthal Until Today; Springer Science & Business Media: Berlin, Germany, 2004. [Google Scholar]
  49. Caselles, V.; Sobrino, J. Determination of frosts in orange groves from NOAA-9 AVHRR data. Remote Sens. Environ. 1989, 29, 135–146. [Google Scholar] [CrossRef]
  50. Peres, L.F.; DaCamara, C.C. Emissivity maps to retrieve land-surface temperature from MSG/SEVIRI. IEEE Trans. Geosci. Remote Sens. 2005, 43, 1834–1844. [Google Scholar] [CrossRef]
  51. Sobrino, J.A.; Raissouni, N. Toward remote sensing methods for land cover dynamic monitoring: Application to Morocco. Int. J. Remote Sens. 2000, 21, 353–366. [Google Scholar] [CrossRef]
  52. Carlson, T.N.; Ripley, D.A. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 1997, 62, 241–252. [Google Scholar] [CrossRef]
  53. Minkina, W.; Dudzik, S. Appendix B: Normal Emissivities of Various Materials (IR-Book 2000, Minkina 2004). In Infrared Thermography: Errors and Uncertainties; Wiley: Oxford, UK, 2009. [Google Scholar] [CrossRef] [Green Version]
  54. Gao, B.C. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 1996, 58, 257–266. [Google Scholar] [CrossRef]
  55. McFeeters, S.K. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 1996, 17, 1425–1432. [Google Scholar] [CrossRef]
  56. McFeeters, S.K. Using the Normalized Difference Water Index (NDWI) within a Geographic Information System to Detect Swimming Pools for Mosquito Abatement: A Practical Approach. Remote Sens. 2013, 5, 3544–3561. [Google Scholar] [CrossRef] [Green Version]
  57. Jiménez-Muñoz, J.C.; Sobrino, J.A.; Gillespie, A.; Sabol, D.; Gustafson, W.T. Improved land surface emissivities over agricultural areas using ASTER NDVI. Remote Sens. Environ. 2006, 103, 474–487. [Google Scholar] [CrossRef]
  58. Olioso, A.; Soria, G.; Sobrino, J.; Duchemin, B. Evidence of Low Land Surface Thermal Infrared Emissivity in the Presence of Dry Vegetation. IEEE Geosci. Remote Sens. Lett. 2007, 4, 112–116. [Google Scholar] [CrossRef]
  59. Willmott, C.J.; Matsuura, K. Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Clim. Res. 2005, 30, 79–82. [Google Scholar] [CrossRef]
  60. Torres-Rua, A.; Aboutalebi, M.; Wright, T.; Nassar, A.; Guillevic, P.; Hipps, L.; Gao, F.; Jim, K.; Alsina, M.M.; Coopmans, C.; et al. Estimation of surface thermal emissivity in a vineyard for UAV microbolometer thermal cameras using NASA HyTES hyperspectral thermal, Landsat and AggieAir optical data. In Proceedings of the SPIE—The International Society for Optical Engineering, Baltimore, MD, USA, 14–18 April 2019; Volume 11008. [Google Scholar] [CrossRef] [Green Version]
  61. Zhang, Y. MODIS UCSB Emissivity Library. 1999. Available online: https://icess.eri.ucsb.edu/modis/EMIS/html/em.html (accessed on 10 September 2019).
  62. FLIR Systems, Inc. FLIR’s Tau 2 Spectral Response Curve; FLIR Systems, Inc.: Wilsonville, OR, USA, 2018. [Google Scholar]
  63. Montanaro, M.; Lunsford, A.; Tesfaye, Z.; Wenny, B.; Reuter, D.; Montanaro, M.; Lunsford, A.; Tesfaye, Z.; Wenny, B.; Reuter, D. Radiometric Calibration Methodology of the Landsat 8 Thermal Infrared Sensor. Remote Sens. 2014, 6, 8803–8821. [Google Scholar] [CrossRef] [Green Version]
  64. Svensson, T.; Hallberg, T. Infrared absorption bands measured with an uncooled interferometric LWIR hyperspectral camera. Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XXIV. Int. Soc. Opt. Photonics 2018, 10644, 106440Y. [Google Scholar] [CrossRef]
  65. Maes, W.H.; Steppe, K. Estimating evapotranspiration and drought stress with ground-based thermal remote sensing in agriculture: A review. J. Exp. Bot. 2012, 63, 4671–4712. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Kanjar, D.; Masilamani, V. Image Sharpness Measure for Blurred Images in Frequency Domain. Procedia Eng. 2013, 64, 149–158. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Flowchart of the proposed Unmanned Aerial Vehicle (UAV)-based method to retrieve precise Land Surface Temperature (LST): (i) Visible and near infrared (VNIR) and long-wave infrared (LWIR) data acquisition; (ii) Reference measurements for VNIR and LWIR data during the flights; (iii) Preprocessing of the acquired data; (iv) Atmospheric correction for the thermal domain; (v) Ortho-mosaic generation; (vi) Co-registration, resampling and layer stacking of the VNIR and brightness temperature (BT) maps; (vii) Computation of the Land Surface Emissivity (LSE) and LST maps.
Figure 1. Flowchart of the proposed Unmanned Aerial Vehicle (UAV)-based method to retrieve precise Land Surface Temperature (LST): (i) Visible and near infrared (VNIR) and long-wave infrared (LWIR) data acquisition; (ii) Reference measurements for VNIR and LWIR data during the flights; (iii) Preprocessing of the acquired data; (iv) Atmospheric correction for the thermal domain; (v) Ortho-mosaic generation; (vi) Co-registration, resampling and layer stacking of the VNIR and brightness temperature (BT) maps; (vii) Computation of the Land Surface Emissivity (LSE) and LST maps.
Remotesensing 12 01075 g001
Figure 2. VNIR-generated colour infrared (CIR) map (NIR-red-green composite) of the rapeseed experiment (white border) as part of the Common Experiment (yellow border) with bare soil patches, at recording time (11 April 2018, 13:42 CET), prepared for maize and summer wheat experiments. The Common Experiment follows crop rotation with homogeneously planted crops; in this image, summer barley (dotted yellow margin) is planted every two years. Furthermore, south of the Common Experiment, sugar beet and, to the north, barley were grown, according to standard agricultural practice. A data logger was placed within the rapeseed field (yellow x; 50.6237N, 6.9873E). The spatial resolution of the CIR map is 7 cm per pixel.
Figure 2. VNIR-generated colour infrared (CIR) map (NIR-red-green composite) of the rapeseed experiment (white border) as part of the Common Experiment (yellow border) with bare soil patches, at recording time (11 April 2018, 13:42 CET), prepared for maize and summer wheat experiments. The Common Experiment follows crop rotation with homogeneously planted crops; in this image, summer barley (dotted yellow margin) is planted every two years. Furthermore, south of the Common Experiment, sugar beet and, to the north, barley were grown, according to standard agricultural practice. A data logger was placed within the rapeseed field (yellow x; 50.6237N, 6.9873E). The spatial resolution of the CIR map is 7 cm per pixel.
Remotesensing 12 01075 g002
Figure 3. The fixed-wing UAV used during our thermal campaigns (left). The sensor platform of the UAV consists of a dual sensor system with a multispectral VNIR sensor, including an irradiance sensing element and a thermal camera (right). The LWIR camera is connected to the autopilot of the UAV to record GPS data and heading information, stored as metadata.
Figure 3. The fixed-wing UAV used during our thermal campaigns (left). The sensor platform of the UAV consists of a dual sensor system with a multispectral VNIR sensor, including an irradiance sensing element and a thermal camera (right). The LWIR camera is connected to the autopilot of the UAV to record GPS data and heading information, stored as metadata.
Remotesensing 12 01075 g003
Figure 4. From air temperature (< T a i r >) and surface temperature ( T s u r f ), the resultant T s u r f _ c from the campaign on 11 April 2018, 13:30–14:20 CET. At 13:55, the surface temperature is higher than the corrected surface temperature due to increasing air temperature is compensated by the mean air temperature ( T a i r _ m e a n ).
Figure 4. From air temperature (< T a i r >) and surface temperature ( T s u r f ), the resultant T s u r f _ c from the campaign on 11 April 2018, 13:30–14:20 CET. At 13:55, the surface temperature is higher than the corrected surface temperature due to increasing air temperature is compensated by the mean air temperature ( T a i r _ m e a n ).
Remotesensing 12 01075 g004
Figure 5. Rapeseed experiment with 63 plots of 3 × 6 meters planted with 10 different genotypes with a minimal repetition of 6 plots. Differences between the varieties are highlighted by the CIR map (a); the NDVI (b) and P v (c) maps were used to estimate the emissivity (11 April 2018, 13:42 CET).
Figure 5. Rapeseed experiment with 63 plots of 3 × 6 meters planted with 10 different genotypes with a minimal repetition of 6 plots. Differences between the varieties are highlighted by the CIR map (a); the NDVI (b) and P v (c) maps were used to estimate the emissivity (11 April 2018, 13:42 CET).
Remotesensing 12 01075 g005
Figure 6. The LSE map (a) represents the emissivity considering P v ; whereas the LSE’ map (b) depicts the emissivity derived by the logarithmic relationship. The dLSE map (c) accentuates the selectivity of both approaches.
Figure 6. The LSE map (a) represents the emissivity considering P v ; whereas the LSE’ map (b) depicts the emissivity derived by the logarithmic relationship. The dLSE map (c) accentuates the selectivity of both approaches.
Remotesensing 12 01075 g006
Figure 7. The temperature distribution between rape plots and surrounding soil areas: B T s e n s map (a) generated from the LWIR imagery (11 April 2018, 13:42 CET). The LST map (b) based on estimated LSE after the corrections were made. The dLST map (c) was computed to quantify the temperature deviation between B T s e n s and the retrieved LST.
Figure 7. The temperature distribution between rape plots and surrounding soil areas: B T s e n s map (a) generated from the LWIR imagery (11 April 2018, 13:42 CET). The LST map (b) based on estimated LSE after the corrections were made. The dLST map (c) was computed to quantify the temperature deviation between B T s e n s and the retrieved LST.
Remotesensing 12 01075 g007
Figure 8. Aberrations of B T s e n s and LST for different assumptions: the d L S T t a u map (a), only including the atmospheric transmittance (case i); the dLS T b k g map (b) only considering the background temperature (case ii) for overcast sky conditions; and the modified dLS T b k g ’ map (c) for clear sky conditions, to highlight the effects of leaving out essential meteorological parameters.
Figure 8. Aberrations of B T s e n s and LST for different assumptions: the d L S T t a u map (a), only including the atmospheric transmittance (case i); the dLS T b k g map (b) only considering the background temperature (case ii) for overcast sky conditions; and the modified dLS T b k g ’ map (c) for clear sky conditions, to highlight the effects of leaving out essential meteorological parameters.
Remotesensing 12 01075 g008
Figure 9. The NDWI map (a) for the border of the water surface, the LSE map of the water surface (b), and the retrieved water temperature L S T p o n d (c). The yellow triangle marks the area in which we measured reference temperatures with thermocouples to validate the thermal data 11 April 2018, 13:32 CET).
Figure 9. The NDWI map (a) for the border of the water surface, the LSE map of the water surface (b), and the retrieved water temperature L S T p o n d (c). The yellow triangle marks the area in which we measured reference temperatures with thermocouples to validate the thermal data 11 April 2018, 13:32 CET).
Remotesensing 12 01075 g009
Figure 10. Correlation of the retrieved pond temperature ( L S T p o n d ) and measured pond temperature ( T p o n d ), with flight 1 (a) on the left and flight 3 (b) on the right 11 April 2018, 13:32 and 16:32 CET, respectively).
Figure 10. Correlation of the retrieved pond temperature ( L S T p o n d ) and measured pond temperature ( T p o n d ), with flight 1 (a) on the left and flight 3 (b) on the right 11 April 2018, 13:32 and 16:32 CET, respectively).
Remotesensing 12 01075 g010
Figure 11. Correlation of the retrieved surface temperature ( L S T r a p e ) and measured surface temperature ( T r a p e ) during flight 1 over the rapeseed experiment (11 April 2018, 13:34–13:49 CET).
Figure 11. Correlation of the retrieved surface temperature ( L S T r a p e ) and measured surface temperature ( T r a p e ) during flight 1 over the rapeseed experiment (11 April 2018, 13:34–13:49 CET).
Remotesensing 12 01075 g011
Table 1. Details of the thermal camera (uncooled microbolometer, thermal resolution: 0.04 K) and the multispectral sensor used. According to the spatial resolution of the sensors, the pixel size (px) for VNIR at a height of about 80 meters is about 7 cm/px and for LWIR 10 cm/px.
Table 1. Details of the thermal camera (uncooled microbolometer, thermal resolution: 0.04 K) and the multispectral sensor used. According to the spatial resolution of the sensors, the pixel size (px) for VNIR at a height of about 80 meters is about 7 cm/px and for LWIR 10 cm/px.
SensorSpectral BandField of ViewResolutionGreyscale
LWIR/TauLWIR: 7.5–13.5 μ m45.0 ° × 37.0 ° 640 × 51214-bit
VNIR/SequoiaG: 550 nm, R: 660 nm,61.9 ° × 48.5 ° 1280 × 96010-bit
(centre wavelengths)RE: 735 nm, NIR: 790 nm
Table 2. Meteorological data (as means) during the flight on 11 April 2018, 13:30–14:20 CET (at 2 m above the ground): wind speed, air temperature, surface temperature ( T s u r f ) at 0.05 m and net radiation. Taking into account small spatial temperature variations, T a i r of both stations (area of flown routes) was averaged (< T a i r >).
Table 2. Meteorological data (as means) during the flight on 11 April 2018, 13:30–14:20 CET (at 2 m above the ground): wind speed, air temperature, surface temperature ( T s u r f ) at 0.05 m and net radiation. Taking into account small spatial temperature variations, T a i r of both stations (area of flown routes) was averaged (< T a i r >).
TimeWind Speed (m s 1 )< T air > ( ° C) T surf ( ° C)Net Radiation (W m 2 )
13:301.1 ± 1.511.2 ± 1.216.6368 ± 24
13:401.0 ± 1.111.2 ± 1.016.4360 ± 16
13:500.9 ± 1.011.5 ± 1.117.0423 ± 35
14:000.9 ± 0.911.7 ± 1.017.6494 ± 89
14:100.5 ± 0.911.9 ± 0.917.7406 ± 126
14:200.6 ± 0.911.9 ± 0.817.4309 ± 17
Means:0.8 ± 1.111.6 ± 1.017.1393 ± 51
Table 3. Meteorological data as means from the weather station (a) and microclimatic station (b) during the flights on 11 April 2018 at 13:20–13:50 CET (flight 1) and 16:20–16:40 CET (flight 3). The background temperature ( T b k g ) and transmittance ( τ ) were used for atmospheric correction.
Table 3. Meteorological data as means from the weather station (a) and microclimatic station (b) during the flights on 11 April 2018 at 13:20–13:50 CET (flight 1) and 16:20–16:40 CET (flight 3). The background temperature ( T b k g ) and transmittance ( τ ) were used for atmospheric correction.
Time T air ( ° C) ω % T bkg ( ° C) τ
Flight 1 (a): pond (time 13:30–13:34)13:2012.4 ± 0.378.7 ± 1.7
13:3012.2 ± 0.277.4 ± 1.3
13:4012.6 ± 0.376.1 ± 1.8
Means:12.4 ± 0.377.4 ± 1.68.8 ± 0.20.95
Flight 1 (b): experiment (time 13:34–13:49)13:3010.0 ± 0.197.1 ± 1.7
13:4010.3 ± 0.191.0 ± 1.0
13:5010.4 ± 0.187.8 ± 1.1
Means:10.2 ± 0.192.0 ± 1.38.8 ± 0.20.95
Flight 3 (a): pond (time 16:30–16:34)16:2013.5 ± 0.272.3 ± 1.9
16:3013.7 ± 0.473.7 ± 1.7
16:4013.6 ± 0.472.5 ± 2.2
Means:13.6 ± 0.372.8 ± 1.9−25.2 ± 0.40.95

Share and Cite

MDPI and ACS Style

Heinemann, S.; Siegmann, B.; Thonfeld, F.; Muro, J.; Jedmowski, C.; Kemna, A.; Kraska, T.; Muller, O.; Schultz, J.; Udelhoven, T.; et al. Land Surface Temperature Retrieval for Agricultural Areas Using a Novel UAV Platform Equipped with a Thermal Infrared and Multispectral Sensor. Remote Sens. 2020, 12, 1075. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12071075

AMA Style

Heinemann S, Siegmann B, Thonfeld F, Muro J, Jedmowski C, Kemna A, Kraska T, Muller O, Schultz J, Udelhoven T, et al. Land Surface Temperature Retrieval for Agricultural Areas Using a Novel UAV Platform Equipped with a Thermal Infrared and Multispectral Sensor. Remote Sensing. 2020; 12(7):1075. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12071075

Chicago/Turabian Style

Heinemann, Sascha, Bastian Siegmann, Frank Thonfeld, Javier Muro, Christoph Jedmowski, Andreas Kemna, Thorsten Kraska, Onno Muller, Johannes Schultz, Thomas Udelhoven, and et al. 2020. "Land Surface Temperature Retrieval for Agricultural Areas Using a Novel UAV Platform Equipped with a Thermal Infrared and Multispectral Sensor" Remote Sensing 12, no. 7: 1075. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12071075

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop