Correcting an Off-Nadir to a Nadir Land Surface Temperature Using a Multitemporal Thermal Infrared Kernel-Driven Model during Daytime

Abstract

Land surface temperature (LST) is a fundamental parameter in global climate, environmental, and geophysical studies. Remote sensing is an essential approach for obtaining large-scale and frequently updated LST data. However, due to the wide field of view of remote sensing sensors, the observed LST with diverse view geometries suffers from inconsistency caused by the thermal radiation directionality (TRD) effect, which results in LST products being incomparable, especially during daytime. To address this issue and correct current off-nadir LSTs to nadir LSTs, a semi-physical time-evolved kernel-driven model (TEKDM) is proposed, which depicts multitemporal TRD patterns during the daytime. In addition, we employ a Bayesian optimization method to calibrate seven unknown parameters in the TEKDM. Validation results using the U.S. Climate Reference Network (USCRN) sites show that the RMSE (MBE) for GOES-16 and MODIS off-nadir LST products is reduced from 3.29 K (−2.0 K) to 2.34 K (−0.02 K), with an RMSE reduction of 0.95 K (29%) and a significant reduction in systematic bias. Moreover, the proposed method successfully eliminates the angular and temporal dependence of the LST difference between the satellite off-nadir LST and in situ nadir LST. In summary, this study presents a feasible approach for estimating the high-accuracy nadir LST, which can enhance the applicability of LST products in various domains.

1. Introduction

The land surface temperature (LST) can reveal the thermodynamic properties of the Earth’s surface and plays an essential role in surface energy balance [1,2], global land warming [3], global water cycle [4], and other geophysical and environmental studies. Thermal infrared remote sensing is the only observation method to generate large-scale, long-time-series, and pixel-averaged LST products [5,6]. However, due to the wide field of view (FOV) of satellite sensors, the viewing geometries are quite different for the same scan line in one image or for the same pixel within one revisit period. Taking the Moderate Resolution Imaging Spectroradiometer (MODIS) as an example, the view zenith angle (VZA) can differ by up to 65° for the same pixel within one revisit period. However, the LST significantly varies with viewing geometries when observing the same scene at the same time, i.e., the widely applicable thermal radiation directionality (TRD) effect [7]. Recent satellite-scale studies reported that the influence of the TRD effect can reach 5 K in urban regions during summer daytime [8,9] and an average of 4 K for sparse vegetation canopies in summer [10]. This severe discrepancy will induce large uncertainty and limit subsequent studies, such as urban heat island and urban planning studies [11,12]. Therefore, it is necessary to correct the current LST products to a reference direction, which is typically the nadir direction [13,14].

To account for the TRD effect in LST products, semi-physical kernel-driven models (KDMs) have been proposed, which consider two anisotropy kernels representing the gap fraction and hotspot effect [13,15,16,17,18]. However, due to the temporal variation in the TRD effect, KDMs with three or four unknown parameters require three or four simultaneous multi-angle LST observations as inputs, which cannot be fulfilled by current satellite sensors [7]. Consequently, rigorous hypotheses have been proposed by reducing the spatial and temporal heterogeneity to obtain more angular information or removing one anisotropy kernel to weaken the complexity of KDMs. They include: (1) the TRD effect is consistent within a single land cover type or spatial cluster [17,19,20]. (2) The TRD effect remains constant for the same pixel throughout the year [18,21]. (3) The TRD effect observed in sparse in situ sites can represent the global-scale TRD effect [14]. (4) The gap fraction effect can be ignored, i.e., discard the temperature difference between soil and vegetation surface components [22,23]. (5) The hotspot effect can also be ignored, i.e., discard the temperature difference between sunlit and sunshade surface components [24]. These hypotheses induce large uncertainty as they discard the basic TRD features: the TRD effect is induced by a heterogeneous temperature distribution within one pixel, including both a gap fraction and hotspot feature, which varies with time and surface structure.

To depict the temporal variation in the TRD effect and achieve TRD correction within a single pixel, Qin et al. [23] introduced the diurnal temperature cycle (DTC) model to the unknown kernel coefficients in the KDM and developed a time-evolved kernel-driven model (TEKDM). However, this model neglected the canopy gap fraction kernel as it was initially designed for a single geostationary satellite with a constant VZA. This condition is not satisfied for polar-orbiting satellite LST products with diverse VZAs. Furthermore, Qin et al. [25] introduced a gap fraction kernel for the TRD correction in an overlapping area of two geostationary satellites and then generated a high-accuracy hemispherical integrated LST. In this study, a TEKDM was employed to estimate the nadir LST within one pixel during the daytime, considering widespread spatial and temporal heterogeneity. Moreover, this multitemporal TEKDM was utilized to correct Geostationary Operational Environmental Satellite-16 (GOES-16) and MODIS off-nadir LSTs to nadir LSTs based on a Bayesian optimization method, which is widely applied in the thermal infrared remote sensing community [26].

The structure of this paper is as follows: Section 2 introduces the satellite LST products and the in situ LST measurements. Section 3 presents the TEKDM used for TRD correction. Section 4 displays the corrected nadir LST results. Section 5 offers a comparison with previous methods and discusses the TRD effect at nighttime and the applicability of the TEKDM. Finally, Section 6 provides the conclusions of this study.

2. Materials

2.1. Remote Sensing LST Products

This study utilized one geostationary and two polar-orbiting LST products to obtain multitemporal observations with different observation geometries. The geostationary satellite provides high-frequency temporal observations with varying solar directions. However, the VZA remains fixed for each pixel. To address this limitation, polar-orbiting satellites were incorporated, which offer varying VZAs ranging from 0° to 65° but their temporal sampling is sporadic. The joint utilization of geostationary and polar-orbiting satellites enhances both temporal and angular information, enabling a robust estimation of the nadir LST.

For this study, the NOAA’s GOES-16 Advanced Baseline Imager (ABI) LST product covering the conterminous United States (CONUS) was employed. This product was generated using a split-window algorithm with the input of top-of-atmosphere thermal infrared radiance measurements [27]. The spatial resolution of this product is 2 km at nadir, and the imaging frequency is 60 min.

Additionally, two LST products (MOD11A1 and MYD11A1) derived from the MODIS sensor onboard the TERRA and AQUA satellites were selected due to their ability to provide diverse viewing angles within a revisit period. The MODIS LST products were generated using a general split-window algorithm, incorporating a classification map-based emissivity as the input [28]. The spatial resolution is 1 km at nadir and the VZA varies from 0° to 65° for a pixel within one revisit period. The combination of GOES-16 and MODIS products provides both temporal and angular information in one day for correcting the off-nadir LST.

2.2. In Situ LST Data

In situ radiometer measurements from the U.S. Climate Reference Network (USCRN) in the western part of CONUS from 2019/12/01 to 2020/11/30 were employed to evaluate the proposed TEKDM. The USCRN radiometers are placed vertically downward and situated approximately 1.3 m above the ground, providing land surface brightness temperature measurements every 5 min [29]. After correcting for the impact of broadband emissivity and reflected atmospheric surface downward radiation, the in situ measured nadir LST is obtained, which is recommended for the evaluation of LST products [29,30,31].

Considering that the observed LST at in situ sites does not usually represent the surrounding areas within the 2 km satellite sensor footprint, the evaluation results suffer large uncertainty. To test the thermal homogeneity within one pixel, the 100 m Landsat-8 LST product was first downscaled to 10 m using the Sentinel-2 visible bands, NDVI, and NDWI based on a linear regression method [32,33]. Then, the standard deviation (STD) of LST values within a 200 × 200 window centered at each site was estimated to assess the 2 km scale spatial representativeness. The STD was calculated only when all the 200 × 200 LST values were under clear-sky conditions after removing cloud and cloud shadow pixels using the Quality Assurance (QA) band.

Those sites with a median STD of less than 1.5 K and a maximum LST of no more than 2 K were used for LST validation. These two thresholds were selected as the balance between the number of valid sites and spatial homogeneity. Finally, six sites in the western region of CONUS were selected. Their locations and photographs are shown in Figure 1a–g. They are situated in sparsely vegetated areas with relatively homogeneous land cover. The boxplots in Figure 1h illustrate the STD values for these sites, with the median and maximum STD values represented by the blue lines and upper whiskers. More detailed information is shown in

Table 1. The detailed information about the six selected USCRN sites.

Site ID	Site Name	Longitude	Latitude	Elevation	Land Cover
S1	AZ Williams 35 NNW	112.34° W	35.76° W	1809 m	GRA
S2	CO La Junta 17 WSW	103.82° W	37.86° W	1342 m	GRA
S3	CO Nunn 7 NNE	104.76° W	40.81° W	1633 m	GRA
S4	NM Socorro 20 N	106.89° W	34.36° W	1479 m	BSV
S5	NM Las Cruces 20 N	106.74° W	32.61° W	1325 m	OSH
S6	OR Riley 10 WSW	119.69° W	43.47° W	1407 m	GRA

, including the elevation and land cover type extracted from the GTOP30 DEM and MCD12Q1 IGBP classification products. There are four grassland (GRA) sites, one barren (BSV) site, and one open shrubland (OSH) site in total.

3. Methods

3.1. Time-Evolved Kernel-Driven Model

The kernel-driven model simulates the directional radiometric temperature using the linear combination of one isotropy kernel and two anisotropy kernels depicting the gap fraction and hotspot effect, which is shown as follows [15]:

$$T({\theta }_s,{\theta }_v,\Delta \phi )=T_N+f_{Gap}\cdot K_{Gap}({\theta }_v)+f_{Hotspot}\cdot K_{Hotspot}({\theta }_s,{\theta }_v,\Delta \phi ,width)$$

(1)

where T(θ_s,θ_v,Δφ) is the observed LST with the sun zenith angle (SZA) of θ_s, the VZA of θ_v, and the relative azimuth angle (RAA) of Δφ. T_N is the nadir LST at time t when the gap fraction kernel (K_Gap) and hotspot kernel (K_Hotspot) are normalized to 0 at nadir. f_Gap is the unknown parameter related to the temperature difference between soil and a leaf. This parameter controls the extent of the gap fraction effect. f_Hotspot and width determine the hotspot signature. f_Hotspot is related to the temperature difference between sunlit and sunshaded surface parts. The width parameter controls the influence region of the hotspot, which is determined by scene structure (e.g., leaf area index, clumping index, leaf inclination angle, canopy height, etc.). The kernel K_Gap represents the angular variation in the gap fraction, which is related to a_soil and is only determined by the VZA [34]. The kernel K_Hotspot depicts the variation in sunlit surface parts, which will be the maximum when the view direction is coincident with the solar direction. In summary, the kernels record the angular variation trend of LST and the scene structure information, while the kernel coefficients retain the component temperature information. Based on this clear separation of kernels and kernel coefficients, Qin et al. [23,25] depicted the temporal variation in kernel coefficients using a DTC model and proposed a general TEKDM modeling framework shown as follows:

$$\begin{array}{c}T({\theta }_s,{\theta }_v,\Delta \phi ,t)=T_N(t)+f_{Gap}(t)\cdot K_{Gap}({\theta }_v)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill +f_{Hotspot}(t)\cdot K_{Hotspot}({\theta }_s,{\theta }_v,\Delta \phi ,width)\end{array}$$

(2)

The time-varying kernel coefficients are parametrized as follows:

$$T_N(t)=T_0+T_a\cdot \cos \left(\frac{\pi }{\omega }(t-t_m)\right)$$

(3)

$$f_{Gap}(t)=A\cdot T_N(t)$$

(4)

$$f_{Hotspot}(t)=B\cdot T_N(t)\cdot \cos {\theta }_s$$

(5)

where T_N(t) is expressed by a DTC model, T₀ is the nadir LST around sunrise time, T_a is the temperature amplitude, ω is the daytime hours, and t_m is the moment when the maximum LST occurs. There are 4 unknown temporal parameters to be estimated. A and B are unknown parameters related to the component temperature difference. By substituting Equations (3)–(5) into Equation (2), we can obtain the TEKDM used in this study:

$$\begin{array}{c}T({\theta }_s,{\theta }_v,\Delta \phi ,t)=\left[T_0+T_a\cdot \cos \left(\frac{\pi }{\omega }(t-t_m)\right)\right]\cdot \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \left[1+A\cdot K_{Vinnikov}({\theta }_v)+B\cdot K_{RL}({\theta }_s,{\theta }_v,\Delta \phi ,width)\right]\end{array}$$

(6)

where T₀, T_a, ω, t_m, A, B, and width are 7 unknown parameters to be estimated with the input of no less than 7 non-real-time LST observations within one day. In addition, the Vinnikov kernel (i.e., K_Vinnikov) and RL kernel (i.e., K_RL) were selected to fill the roles of gap fraction and hotspot kernel, respectively [18,22], which are expressed as follows:

$$K_{Vinnikov}({\theta }_v)=1-\cos {\theta }_v$$

(7)

$$K_{RL}({\theta }_s,{\theta }_v,\Delta \phi ,k)=\frac{e^{-kf}-e^{-k{f}_N}}{1-e^{-k{f}_N}}$$

(8)

$$f=\sqrt{{\tan }^2{\theta }_s+{\tan }^2{\theta }_v-2\tan {\theta }_s\cdot \tan {\theta }_v\cos \Delta \phi }$$

(9)

$$f_N=\tan {\theta }_s$$

(10)

where k is the hotspot width parameter.

The shapes of K_Vinnikov and K_RL in the solar principal plane (SPP) with the SZA of 40° are shown in Figure 2. The K_Vinnikov is a bowl-shaped function symmetrical about the VZA (Figure 2a), which can capture the pattern of the gap fraction function well [34]. In Figure 2b, the adjustable hotspot width parameter k in K_RL improves the fitting accuracy in the hotspot region to better depict the hotspot signatures [15]. Both K_Vinnikov and K_RL were normalized to 0 at the nadir.

3.2. Solving the TEKDM

To solve the 7-parameter model with a relatively high degree of freedom, a Bayesian optimization algorithm was adopted in this study. Bayesian optimization can synthesize the prior knowledge for unknown parameters and the information from observation data. To achieve this, a Markov chain Monte Carlo (MCMC) method was used to solve the non-linear optimization problem. This method can reduce the influence of observation errors and produce a robust estimation of unknown parameters, meaning it is widely employed in thermal infrared remote sensing studies [26]. The initial values and boundaries of the unknown parameters are shown in

Table 2. The initial values and boundaries of 7 unknown parameters in the TEKDM.

Parameter	Initial Value	Lower Bound	Upper Bound
T0	$T_0^{\prime }$	$T_0^{\prime }$ − 5	$T_0^{\prime }$ + 5
Ta	$T_a^{\prime }$	$T_a^{\prime }$ − 5	$T_a^{\prime }$ + 5
tm	$t_m^{\prime }$	$t_m^{\prime }$ − 1	$t_m^{\prime }$ + 1
ω	${\omega }^{\prime }$	${\omega }^{\prime }$ − 1	${\omega }^{\prime }$ + 1
A	0.015	0	0.03
B	−0.015	−0.03	0
k	0.5	0.0001	1

.

$T_0^{\prime }$, $T_a^{\prime }$, $t_m^{\prime }$, and ${\omega }^{\prime }$ (i.e., the initial values of T₀, T_a, t_m, and ω) in Table 2 are the parameters for a daytime DTC model directly driven by LSTs before TRD correction according to Hong et al. [35]. The parameters A and B control the amplitude of the TRD effect, and we set them to a broad region. When T_N equals 300 K, the maximum influence of the gap fraction and hotspot effect reaches 9 K, which is relatively large according to Cao et al. [7]. The boundaries of k cover a broad region according to Ermida et al. [17]. After solving 7 unknown parameters in the TEKDM, the f_Gap and f_Hotspot can be acquired, and the nadir LST was then estimated using Equation (1) with the input of directional LST.

3.3. Evaluation of the Nadir LST

An in situ radiometer records the radiance from the land surface at nadir. The recorded brightness temperature should be corrected to the LST using the Stefan–Boltzmann law with the inputs of broadband surface emissivity and downward longwave radiation. Therefore, the LST is calculated by the following equation:

$$T_s={\left(\frac{{\sigma }_s{T}_b^4-\left(1-\epsilon \right)\cdot R^{\downarrow }}{\epsilon {\sigma }_s}\right)}^{1/4}$$

(11)

where T_s is the in situ measured LST, T_b is the recorded brightness temperature, ε is the broadband emissivity, σ_s is the Stefan–Boltzmann constant (i.e., $5.67\times {10}^{-8}$), and ${\mathrm{R}}^{\downarrow }$ is the downward longwave radiation. In this study, ε is estimated using a vegetation cover method [36,37] with the inputs of the ASTER GED soil broadband emissivity and the land-cover-based vegetation broadband emissivity. $R^{\downarrow }$ is extracted from the ERA5-Land 0.1° product.

Before TRD correction using the TEKDM, outliers were first removed. The outliers widely exist in satellite LST products because operational cloud mask products cannot remove all the clouds, especially small clouds within one pixel. Therefore, a widely employed “3σ-Hampel identifier” method was adopted to reduce this uncertainty [38,39]. This method assumes that the error of LSTs follows a normal distribution, and those points that deviate more than 3σ from the median value are identified as outliers. σ was calculated by the following equation:

$$\sigma =1.4826\times median\left\{\left|x_i-x_m\right|\right\}$$

(12)

where x_m is the median value of ∆LST (i.e., satellite LST minus in situ LST) and x_i is the i-th data point. This method can mitigate the influence of outliers when calculating standard deviation and has been widely employed in LST product evaluation studies [39,40].

To evaluate the off-nadir LST and the corrected nadir LST, the root mean square error (RMSE) and mean bias error (MBE) were selected in this study. The RMSE quantifies the average distance between satellite LST values and reference in situ LST values. The MBE depicts the bias of satellite LST values.

The general process for correcting the TRD effect is shown in Figure 3. Firstly, the GOES-16 and MODIS LST data were extracted using a nearest neighbor method to match the in situ LST; this was acceptable as the impact of resolution difference was weak, which we knew because the 2 km scale thermal homogeneity of in situ sites was carefully evaluated. Moreover, the GOES-16 LST was corrected to the MODIS LST using a linear regression method (i.e., T_MODIS = a * T_GOES + b). The parameters a and b were calibrated using nighttime LST matchups when the absolute difference in VZA was less than 15° and the time difference was no more than 30 min. This relationship was further applied to daytime LSTs. Then, the corrected off-nadir LSTs, viewing geometries (i.e., VZA, SZA, and RAA), and local solar times were jointly applied to drive the multitemporal kernel-driven model using a Bayesian optimization method. After calibrating the 7 unknown parameters in Equation (6), f_Gap and f_Hotspot can be obtained to estimate the nadir LST with the input of the off-nadir LST.

4. Results

The accuracies of the GOES-16 and MODIS daytime LST products when compared against in situ nadir LSTs before and after TRD correction are presented in this section. Additionally, the dependence on the VZA when calculating the LST difference between MODIS off-nadir and corrected nadir LSTs is analyzed in detail. Due to the temporal variation in TRD amplitude, the corresponding accuracy variation is also analyzed. Finally, the spatial distribution of off-nadir and nadir LSTs is shown, to illustrate TRD correction principles.

4.1. The Accuracy of Off-Nadir and Nadir LST

A direct comparison between the in situ LST and satellite LST before (off-nadir LST) and after (nadir LST) TRD correction is shown in Figure 4. Notably, a significant underestimation can be observed in both the GOES-16 and MODIS LST products, particularly for LST values exceeding 300 K (Figure 4a–c). The RMSE and MBE for off-nadir LST are 3.29 K and −2.0 K for all LST values (Figure 4a). Following TRD correction, the corrected nadir LST exhibits an RMSE reduction of 0.95 K (29%) and a significant reduction in absolute MBE (Figure 4d), with an RMSE (MBE) of 2.34 K (−0.02 K). The histogram in Figure 4a,d demonstrates the effective reduction in systematic bias and variance in ∆LST (i.e., the satellite LST minus in situ LST). This correction trend is consistent for both GOES-16 and MODIS off-nadir LST products. Specifically, the RMSE of GOES-16 LST is reduced from 3.44 K to 2.33 K with a reduction of 1.11 K (32%), while the MBE is reduced from −2.22 K to −0.09 K after TRD correction. Similarly, the RMSE (MBE) for MODIS LST decreases from 2.66 K (−1.18 K) to 2.35 K (0.21 K) with an RMSE reduction of 0.31 K (12%) and a notable elimination of systematic bias. Overall, the proposed correction method removed the severe underestimation in off-nadir LST products, achieving a high accuracy after TRD correction.

4.2. The Angular Dependence of MODIS Off-Nadir and Nadir LST

MODIS provides LST values with a VZA ranging from 0° to 65°, resulting in a significant gap fraction effect. Figure 5 displays the distribution of ∆LST (i.e., satellite LST minus in situ LST) under different VZAs for MODIS LST products before and after TRD correction. The off-nadir ∆LST demonstrates a declining trend with the increase in VZA (Figure 5a). In the case of a pixel mixed with vegetation and soil components, a large VZA leads to more vegetated surface parts with lower temperatures observed. Consequently, the off-nadir LST is lower than the nadir LST, with the maximum average ∆LST approaching −3 K (Figure 5a). After TRD correction, the mean ∆LST consistently varies around 0 K, successfully eliminating the dependence on VZA (Figure 5b).

4.3. The Temporal Variation in Off-Nadir and Nadir LST

The TRD effect arises from the inhomogeneous temperature distribution within one pixel. Temporal variation in component temperatures contributes to varying TRD amplitudes. In addition, this variation happens on various temporal scales, which cover hourly, daily, and monthly scales.

Figure 6 shows the in situ LST, satellite off-nadir LST, and corrected nadir LST for a specific open shrubland site (i.e., the NM Las Cruces 20 N site with the maximum number of clear-sky observations) during four selected days in the middle of each season. This hourly scale result shows that the off-nadir LST generally exhibits lower values, while the corrected nadir LST effectively mitigates this underestimation. In addition, the temporal variation trend of nadir LST aligns well with that of in situ LST as the uncertainty caused by the angular effect is removed. However, a discrepancy of approximately 2 K still remains on 2020/01/13 at around noon (Figure 6a). Generally, the corrected nadir LST is more consistent with the in situ LST, indicating the validity of the proposed TEKDM.

Figure 7a shows detailed information about the daily variation in ∆LST. The results show that the off-nadir ∆LST exhibits a significant seasonal variation trend. The ∆LST is usually within −2 K in autumn and winter but exceeds −2 K during spring and summer, with the maximum absolute value approaching −6 K in summer. In contrast, the corrected ∆LST generally varies within the range of ±2 K over one year, indicating a significant reduction in temporal dependence.

Figure 7b shows the variation in ∆LST across different months. Generally, the off-nadir LST exhibits a declining ∆LST from winter (−0.85 K ~ −1.15 K) to spring (−1.85 K ~ −3.50 K) and summer (−3.11 K ~ −3.75 K), followed by an increasing trend towards autumn (−0.88 K ~ −2.26 K). After TRD correction, the corrected nadir LST aligns well with the in situ measured nadir LST. The corrected ∆LST varies around 1 K and remains relatively constant across different seasons. Furthermore, the proposed method exhibits larger correction in summer and less correction in winter, indicating its ability to capture the seasonal variation pattern of the TRD effect.

4.4. The Spatial Distribution of Off-Nadir and Nadir LST

To exhibit the corrected nadir LST more clearly, Figure 8 shows the spatial distribution of the LST before and after TRD correction and some relevant auxiliary variables. A scene of the MOD11A1 LST product (h08v05) on 2020/06/09 with lots of clear-sky pixels was selected. To improve computational efficiency and reduce uncertainties caused by the resolution difference between MODIS and GOES-16 products, the LST images were resampled to 0.05°.

Figure 8a shows the MOD11A1 off-nadir LST product. The VZA and ξ (the angle between solar and view direction) are large in the lower right part, as shown in Figure 8d,e. After TRD correction, the nadir LST (Figure 8b) is much higher in this region than the directional LST (Figure 8a). This temperature difference reaches 7 K, as shown in Figure 8c. On the one hand, the large VZA induces a severe gap fraction effect as more vegetation surface parts with lower temperatures were viewed, resulting in an underestimation for the off-nadir LST compared to the nadir LST. On the other hand, the ξ became smaller after TRD correction, as shown in Figure 8e,f (note that the SZA is the ξ between the solar and nadir directions). With the decrease in ξ, more sunlit surface parts with higher temperatures were viewed, resulting in higher LST values due to the enhancement of the hotspot effect. Therefore, the nadir LST is much higher than the off-nadir LST in the lower right part with the combined influence of the gap fraction and hotspot effect. In the upper left part of Figure 8a, the VZA is around 0° to 20° (Figure 8d), and the ∆LST is small, as shown in Figure 8c. This weak correction is attributed to the observed LST being close to nadir, thus requiring a weak correction. In general, the TRD correction is obvious when the difference in VZA or ξ is large between off-nadir and nadir directions. The spatial distribution of the corrected nadir LST is consistent with the cause of the TRD effect.

5. Discussion

5.1. Comparison with Previous TRD Correction Method

A detailed comparison between the proposed TEKDM and a TRD correction method proposed by Hu et al. [19] is shown in

Table 3. The accuracies of the original LST product and the two TRD correction methods.

Product	Metric	Off-Nadir LST	TEKDM	Hu’s Method
All Data Points	RMSE (K)	3.29	2.34	2.37
	MBE (K)	−2.0	−0.02	−0.68
GOES-16	RMSE (K)	3.44	2.33	2.41
	MBE (K)	−2.22	−0.09	−0.82
MODIS	RMSE (K)	2.66	2.35	2.22
	MBE (K)	−1.18	0.21	−0.17

. Hu’s method employed the calibrated KDM parameters from Ermida et al. [13,17], which were estimated with the inputs of SEVIRI and MODIS simultaneous LST observations over 15 spatial clusters in 2011. This method was then extended to correct the TRD effect for a single LST product [19]. Table 3 shows that the two TRD correction methods mitigated the angular effect successfully. The RMSE after correction almost remains the same. However, there is still a large underestimation for Hu’s method (i.e., an MBE value of −0.68 K) compared to the TEKDM. Generally, the RMSE is similar for the two methods, but the TEKDM can better reduce the MBE value.

5.2. The TRD of the Nighttime LST Product

In this study, the proposed method only employed the daytime LST to solve the TEKDM. However, the nighttime LST may contain valuable angular information, which is ignored in this study. Figure 9 shows the accuracy of nighttime LST products and their angular dependence for the MODIS LST. The RMSE of the nighttime LST is 1.74 K (Figure 9a), which is consistent with the reported product accuracy [41]. However, the nighttime LST products exhibit an underestimation with the MBE of −1.02 K, which may be attributed to uncertainties in emissivity, undetected clouds, and the TRD effect. After sunset, the hotspot effect tends to diminish, and the TRD effect is dominated by the gap fraction effect. Figure 9b shows the influence of the gap fraction effect in MODIS LST products. With the increase in VZA, the ∆LST shows a slight downward trend, and the average ∆LST varies by approximately 1 K. Numerous studies have employed this nighttime TRD information, assuming that the gap fraction kernel coefficient remains constant both during the day and night [17,18]. However, the reliability of this assumption requires further comprehensive investigations. Adopting reasonable hypotheses and developing high-accuracy models that integrate both daytime and nighttime angular information is a promising approach to generating nadir LST products with enhanced temporal consistency on a daily scale.

5.3. The Applicability of the TEKDM

The proposed TEKDM requires at least seven clear-sky LST observations to calibrate unknown parameters, which is hard to satisfy as approximately 68% of the Earth’s surface is covered by clouds each day [42]. Therefore, more practical methods are desired to make the model more applicable. Jointly utilizing multi-source LST products from various geostationary and polar-orbiting satellites is helpful to access more LST observations. In addition, developing a TEKDM that can employ nighttime LST values can also mitigate the lack of daytime LST values. Moreover, assessing and adopting various hypotheses to acquire more angular information without much loss of accuracy is a feasible way to solve the TEKDM. With the development of machine learning algorithms, combining TEKDM and data-driven methods has the potential to achieve TRD correction in a single moment, which deserves further investigation. Generally, incorporating more information and developing advanced models are helpful to achieve TRD correction at a global scale and further promote subsequent applications.

In addition to the limitation concerning clear-sky observation numbers, the TEKDM also suffers in terms of spatial applicability over diverse land surfaces. The employed Vinnikov–Chen model in the TEKDM is derived under the gap fraction and hotspot framework, which was originally designed for vegetation and soil structure with a random distribution [15]. This framework may not accurately capture the TRD pattern over cropland covers with a regular spatial distribution or over non-vegetated areas such as bare soil, sand, permanent snow and ice covers, urban areas, water bodies, and so on. Currently, many TRD models have already been developed in these areas, which can be integrated into the TEKDM [43,44]. Additionally, the performance of the TEKDM over complex land surfaces such as heterogeneous scenes and rugged terrains deserves more comprehensive studies.

6. Conclusions

In this study, we employed a multitemporal TEKDM and a Bayesian optimization method to correct GOES-16 and MODIS off-nadir LSTs to nadir LSTs within a single pixel during the daytime. The main conclusions can be summarized as follows:

• The evaluation using the USCRN in situ nadir LST showed that the proposed TEKDM effectively reduced the RMSE (MBE) of off-nadir LST products from 3.29 K (−2.0 K) to 2.34 K (−0.02 K), resulting in an RMSE reduction of 0.95 K (29%) and a significant reduction in systematic bias.
• The TEKDM successfully eliminates the angular and temporal dependence when calculating the LST difference between the original off-nadir LST and the in situ nadir LST, which indicates that the TEKDM can depict the TRD patterns well and adapt to their temporal variations.

In summary, the proposed TEKDM can generate nadir LST values with high accuracy and consistency. However, this model is driven by at least seven daytime LST observations, which severely hinders its widespread application. Therefore, developing a TRD correction model suitable for both daytime and nighttime observations with fewer unknown parameters and investigating its spatial applicability will help to generate time-continuous nadir LST products, which can find relevant applications where high-accuracy LST products are needed.