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

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. Remote Sens. 2020, 12, 1075; doi:10.3390/rs12071075 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 1075 2 of 27

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.

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.

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

Lst Retrieval Methodology
Based on the radiative transfer model (RTM), the radiation path through the troposphere can be written as [17,29]: where L sens is the at-sensor radiance, L sur f the surface radiance attenuated by the atmosphere, and L atm the upwelling radiance emitted by the atmosphere. L re 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 sens 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 sens ) into the at-sensor brightness temperature (BT sens ). 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). BT sens 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]: where M is the irradiance (W m −2 ), is the emissivity (ranging from 0-1; dimensionless), σ is the , 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): where T sur f ( , τ) is the retrieved surface temperature, BT sens the at-sensor brightness temperature, T bkg the background temperature, and T air 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 air , T bkg , 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].

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 BT sens 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: where ω is the water vapour content (mm), ω% is the relative humidity (ranging from 0-1; dimensionless), T air 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: 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 atm = 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.

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]: where N IR and R are the reflectances in the near-infrared and red spectra, respectively (−1...+1). Moreover, NDV 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 (NDV I < NDV I soil ), full vegetation (NDV I > NDV I veg ), 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: where NDV I soil and NDV I veg 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 ( ): where veg (=0.988) and soil (=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: 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: where veg is maintained as the global emissivity of dense canopies [30,31] and soil 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 NDV I soil value of 0.157 (standard deviation, SD = 0.227) was derived by the means of four samples from the rapeseed field studied. A NDV I veg 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 NDV I veg 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)): 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).

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: In our second assumption, which presupposes a completely transmissive atmosphere, we set the transmittance to one (τ ! = 1). This simplifies the same initial formula to: 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 BT sens and LST(τ) or LST( ) are shown in the results (Section 3).

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

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. 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.
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 (dT cam ) under stable ambient conditions was dT cam = ±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]. 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. 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.

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 air ), relative humidity (ω%), surface temperature (T sur 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 sur 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 sur 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 air will directly increase T sur f in analogical degree, but also change the net radiation, affecting the surface temperature. The interdependence of T air and T sur f as an approach to surface temperature correction (T sur f _c ) has been proposed by [29]: where T sur f _c is the corrected surface temperature, T sur f , T air = <T air >, and T air_mean is the mean of <T air > ( • C). The correction of T sur f compensates for the noise caused by drifts in air temperature ( Figure 4). This approach (Equation (13)) was used to account for small discrepancies (dT sur f ) in the thermal measurements caused by the impact of <T air >. We calculated the differences of T sur f and T sur f _c , which were defined as dT sur f . There were no significant changes in <T air > and T sur f during the campaign and, so, dT sur f = 0.28 ± 0.1 • C was small. Due to stable weather conditions (Table 2), the averaged values of T bkg (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 bkg ), due to ('cold') clear sky conditions. The transmittance (τ) was always about 0.95. Both T bkg 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]: where G and N IR are the reflectances in the green and near-infrared spectra, respectively (−1...+1). According to [56], we used an NDW I threshold of NDW 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 pond ) 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 pond 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 (LST pond ). According to Equation (3), the captured thermal images of the pond were used to retrieve LST pond . 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 pond and LST pond were used to calculate the temperature deviations (dT pond ) during the flight.

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 (NDV I soil = 0.157; SD = 0.227) and full vegetation (NDV I veg = 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 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 BT sens 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, BT sens 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 BT sens (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 BT sens 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 BT sens (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 BT sens 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 bkg = 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 bkg 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 bkg ). 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 (LST pond ) 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 (LST pond ) to validate the thermal data using the measurements with the pond thermocouples (T pond ). Both LST pond and T pond were used to calculate the mean temperature deviations (dT pond ) 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 dT pond_1 = 0.53 K, with a coefficient of determination of R 2 = 0.985 (n = 24). For flight 3 (Figure 10b), the MAE was dT pond_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 (dT rape ) also confirmed the high performance of the LST retrieval ( Figure 11). The MAE was dT rape = 0.41 K and R 2 = 0.879 (n = 27). Figure 11. Correlation of the retrieved surface temperature (LST rape ) and measured surface temperature (T rape ) during flight 1 over the rapeseed experiment (11 April 2018, 13:34-13:49 CET).
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.

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 (NDV I < NDV I soil ), dense canopy (NDV I > NDV I veg ), and mixed states (Equation (9)). For dense canopy, we suggest a global NDV I veg 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 NDV I soil 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 ( veg = 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 soil = 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 bkg ) 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 veg = 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 dT pond_1 = 0.53 K (R 2 = 0.985) and dT pond_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.

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:

1.
Construction of an ortho-mosaic for NDVI maps, LST maps, and emissivity maps; 2.
Atmospheric correction, including the background temperature; and 3.
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: 1.
Low-cost UAV campaigns for the thermal mapping of large areas with high repetition rates; 2.
Irrigation management, crop water stress monitoring, and plant disease detection; and 3.
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.
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 bkg ) 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 bkg value obtained from the aluminium panel [65]. Averaged temperatures of T bkg 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. 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.  (3)) From the radiative transfer model (RTM) to the at-sensor temperature (T sens ) [29,43]: where term 1 is the emitted radiance from the surface ( · L sur f ), term 2 the downwelling radiance (L re f l ↓) reflected at the surface (1 -), and term 3 the emitted radiance from the atmosphere (L atm ↑); 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 sens = (L sur f + L re f l ↓) · τ + L atm ↑ .
The STEFAN-BOLTZMANN law: