Long‐Term Trend and Seasonal Cycles of Gap‐Free Downscaled Diurnal/Nocturnal LST and the Interaction to Functional Plant Trait Under Tropical Monsoon Climate

Abstract

Land surface temperature (LST) monitoring via Earth observation constellation will become optimized and consistent with spatiotemporal‐explicit characteristics. Besides, scientific evidence for the interaction between LST and vegetation biophysical variables remains limited through spatial large‐scale assessment and seamless long‐term tracking. This study addresses this gap by utilizing gap‐filled fine spatial resolution LST products in understanding the dynamic over the period 2000–2023 and the spatiotemporal relationship with leaf area index (LAI). Firstly, Moderate Resolution Imaging Spectroradiometer (MODIS) LST 1,000 m of both daytime and nighttime were downscaled to a finer resolution of 250 m using the Random Forest algorithm. The Whittaker algorithm was then applied to obtain gap‐free LST products due to the typical cloud cover under tropical monsoon climate. Time series decomposition of gap‐filled fine resolution LST revealed slight warming trends in daytime (0.005°C year⁻¹), nighttime (0.036°C year⁻¹), and mean of all‐day time (0.02°C year⁻¹) over recent 24 years, while seasonal amplitude in daytime (−3.7°C–4.8°C) is more fluctuated than in nighttime (−2.5°C–1.9°C). Spatial correlations of monthly LSTs and LAI indicated a consistent negative correlation (R ranging from −0.717 to −0.45). These findings shed light on the quantitative relationship between vegetation LAI and LST, contributing to a more unified theoretical framework for understanding functional vegetation responses under diverse climatic conditions.

Introduction

Climate extremes nowadays are no longer hypothetical scenarios but also are actually happening worldwide. The summer of 2023 was determined as the hottest season since the global mean temperature was recorded in 1880 by NASA (Karen Fox et al., 2023), and the warming signal seems to be continuing (Tollefson, 2023). Belonging one of the 55 Essential Climate Variables (ECVs) of the Global Climate Observation System (GCOS, 2018), land surface temperature (LST) plays an indispensable role in the studies concerning global warming, as well as vegetation stand traits response and dynamics (W. Kustas & Anderson, 2009; Z. L. Li et al., 2023), improve comprehensive understanding and move toward empowering us to treat Earth ecosystems in the changing climate context.

To date, satellite-borne thermal infrared (TIR) instruments provided an unprecedented data stream with diverse temporal and spatial resolution through long-term observations. Notably, a major number of LST-related applications focused on diurnal products, while nocturnal products have received less attention, primarily due to the rare availability of TIR sensors captured during nighttime (Qi et al., 2020; Yoo et al., 2022). Among widely available and freely accessible TIR sensors, Terra/Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) has continuous global tracking compared to the remaining higher resolution TIR sensors such as Landsat generations, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), or recent Sentinel 3 Sea and LST Radiometer (SLSTR) (Z. L. Li et al., 2023; Pu & Bonafoni, 2023; Trigo et al., 2008).

Broadly speaking, retrieval of LST through Earth observation TIR sensors can be done by two primary approaches: (a) assumption of land surface emissivity (LSE) is known (e.g., single-channel, split-window, dual-window), and (b) simultaneously retrieving LST and LSE associated with available atmospheric information (e.g., temperature and emissivity separation, physics-based day/night) (Cao et al., 2019; Z. L. Li et al., 2023). For most of the available satellite-based LST products presently, coarse spatial resolutions are likely acceptable for interpretations of climate phenomena by global or continent scales, while limited for regional or local applications (Pu & Bonafoni, 2023; Zhan et al., 2013). Additionally, high pixel sizes may represent inadequately the effect of LST across land cover types especially in heterogeneous landscapes (Ebrahimy & Azadbakht, 2019). Consequently, to address fine-scale applications, image sharpening offers a suitable solution to overcome these remarkable challenges.

Focus on spatial sharpening methods, either downscaling or fusion can be addressed effectively, yet the downscaling based on kernel-driven regression was applied more frequently (Pu & Bonafoni, 2023). Pioneer linear regression method is the disaggregation of radiometric surface temperature (DisTrad) (W. P. Kustas et al., 2003), which utilizes a second-order polynomial relationship between Normalized Different Vegetation Index (NDVI) and LST. This approach was further refined by incorporating a stronger relationship between fractional vegetation cover and LST, leading to the development of thermal imagery sharpening (TsHARP) (Agam et al., 2007). Hybrid approaches such as integrating DisTrad with an artificial neural network algorithm (Bindhu et al., 2013), or a combination of geographically and temporally neural network weighted autoregression model (Wu et al., 2022) were introduced with the motivation of keeping advantages while tackling disadvantages in each isolated method.

The era of artificial intelligence witnessed a plethora of machine learning algorithms applied to LST-related applications, particularly demonstrated for improvement in downscaling performance compared to long-established methods (Ebrahimy & Azadbakht, 2019; Guo et al., 2022; Hutengs & Vohland, 2016). In general, nonlinear have outperformed linear methods (Guo et al., 2022). Comparisons between several machine learning algorithms resulted in different concluding remarks depending on feature selection and environmental conditions (Ebrahimy & Azadbakht, 2019). Among potential candidates, Random Forest (RF) (Breiman, 2001) was recognized as a straightforward and robust algorithm for LST image sharpening (Ebrahimy & Azadbakht, 2019; Ebrahimy et al., 2021; Hutengs & Vohland, 2016; Yoo et al., 2022), even remote sensing applications generally (Belgiu & Drăgu, 2016).

While image sharpening procedures can tackle coarse spatial resolution, the downscaled products still lack spatial integrity due to cloud covers. Here, the phenomenon is more serious over tropical climate regions with clouds presenting year-round, and tremendous frequency in rainy seasons. Hence, missing data reconstruction is a definite post-processing technique to implement in order to obtain gap-filled products (Mo et al., 2021). Several well-established data reconstruction methods have been proposed for generating gap-free LST data, including harmonic analysis of time series (Wei et al., 2021), the Savitzky-Golay (Y. Li et al., 2018), and the Whittaker algorithm (Xie & Fan, 2021). Notably, the Whittaker algorithm has emerged as a versatile choice for both gap-filling and smoothing time series data. This technique utilizes a controlled threshold to achieve different outcomes tailored to specific objectives, such as noise reduction or smoothing of long-term time series at regional or planetary scales (Kovács et al., 2023).

Achieving LST products with fine spatial resolution through a downscaling procedure, accompanied by the reconstruction of missing data, enables advanced experiments to understand how LST interacts with vegetation biophysical variables, which is currently still limited through spatial large-scale assessment and long-term tracking. Offering insights into these relationships enables the capabilities toward different vegetation-based applications such as crop calendar optimization, and forest monitoring come along with urban spatial planning (Pu & Bonafoni, 2023; Zhan et al., 2013). Particularly in the satellite-based applications of terrestrial ecosystems, the relevance between two crucial climate and vegetation biophysical variable both belonging to ECVs, namely LST and leaf area index (LAI) have been explored with different fundamentally scientific backgrounds following specific conditions and retrieved methods. Understanding vegetation-related biophysical climate effects is crucial for accurately modeling atmospheric and terrestrial ecosystems (Y. Li et al., 2023). For instance, focused on mixed temperate forests depicting heterogeneous landscapes, LST has a weak relationship with LAI, using retrieved products from Unmanned Aerial System (Stobbelaar et al., 2022) and Landsat-8 (Neinavaz et al., 2019). Based on annual MODIS products of 2003–2018, weak to strong links between LST and LAI were found depending on individual continents (Rasul et al., 2020). In contrast, winter wheat LAI (Sentinel-2) showed a positive correlation with LST during nighttime (ASTER) while no relationship during daytime (Landsat-8) (Abdullah et al., 2020).

Within the mentioned perspectives, the objectives imposed in the present study are twofold. First, we acquire gap-filled fine spatial resolution LSTs by downscaling coarse-resolution LST data (1,000 m) to a finer spatial scale (250 m) and implementing the Whittaker algorithm to fill gaps caused by cloud cover. Second, we investigate the relationship between daytime and nighttime LST and functional trait LAI under tropical monsoon climate, analyzing spatial and temporal variations over the recent 24 years (2000–2023). Our analysis is implemented over a highland region characterized by diverse land cover types and varying altitudes. The comprehensive quantification deepens our understanding of ecosystem responses in tropical climate conditions across a range of environmental settings.

Materials and Methods

An illustration of all processing and analysis steps is depicted in Figure 1, with three major workloads of data collecting, processing, and analysis. Details are described further in the next sections.

[IMAGE OMITTED. SEE PDF]

Study Area

Located in the center of Vietnamese Central Highlands, Dak Lak is a well-known province for coffee production with geographical coordinates between 12°9'45"N and 13°25'06"N latitude and 107°28'57"E to 108°59'37"E longitude. The study site covers an area of around 13,125 km², with elevation varying from 112 to 2,428 m. Characterized by a tropical monsoon climate background, the region experiences a pronounced seasonal transition with a significant increase in precipitation. In statistics, the wet (or rainy) seasons (May to October) account for 90% of the annual rainfall compared to dry seasons (November to April). The area also represents different land cover types (i.e., forest, bare soil, agriculture, savannas, and urban) exhibited over complex terrain surfaces (Figure 2).

[IMAGE OMITTED. SEE PDF]

Downscaling LST Data Set

Both diurnal and nocturnal LSTs were obtained from MODIS daily 1 km LST V6.1 products (MOD11A1) onboard Terra constellation, considered as target variables for downscaling from 1 km to 250 m spatial resolution. To fully reflect the relationship between output variable LST and predictor variables, we involved seven explanatory factors to downscaling models including (a) elevation, (b) slope, (c) aspect, (d) land cover, surface reflectance in (e) red and (f) near-infrared (NIR) bands, and (g) NDVI, belonging to groups of topography, land cover, and spectral characteristics.

It should be noted that so far there are no standard regulations for the appropriate determination of explanatory variables as input kernels for downscaling models. All in all, the effect of topographical factors (i.e., elevation, slope, aspect) and surface attributes (i.e., land cover, NDVI) were considered in our downscaling procedure since the strong relationship with both diurnal LSTs (Ebrahimy & Azadbakht, 2019; Zakšek & Oštir, 2012) and nocturnal LSTs (Qi et al., 2020; Yoo et al., 2022). Besides, visible-near infrared bands (i.e., NIR and red) can also be used due to the stability of spatial variations in surface characteristics from day to night (Qi et al., 2020). Moreover (Qi et al., 2020; Yoo et al., 2022), have reasonably proved the use of bands acquired in the daytime as input kernels for downscaling nocturnal LSTs. To facilitate direct access and analysis of multi-sensor remotely sensed data, we leveraged the cloud-based computing capabilities of Google Earth Engine (GEE). Topography-related maps encompass elevation, slope, and aspect were created based on the Advanced Land Observing Satellite (ALOS) digital surface model (30 m) (Tadono et al., 2014). The annual land cover was provided by MCD12Q1 MODIS products (500 m) (Strahler et al., 1999). Daily surface reflectance in red (620–670 nm) and NIR (841–876 nm) wavelength were obtained by Terra surface reflectance daily products at 250 m spatial resolution (MOD09GQ.061) (Vermote & Vermeulen, 1999). Subsequently, NDVI was calculated from these two mentioned bands. The nearest neighbor resampling technique was employed to achieve spatial alignment between all independent variables (i.e., 250 and 1,000 m). This ensured the compatibility of independent variables when co-located with the 1,000 m resolution LST data to establish the regression model, and predicted the target 250 m resolution LST.

LAI Data Set

LAI data set, developed and processed by the Copernicus Global Land Service (CGLS), was utilized to characterize vegetation stand traits (Fuster et al., 2020). The CGLS LAI products were derived from imagery acquired by the Project for On-Board Autonomy—Vegetation (PROBA-V) and Sentinel-3 Ocean and Land Color Instrument (OLCI) satellites using neural network algorithm (Baret et al., 2016; Verger & Descals, 2022). We downloaded all available 300 m resolution products covering January 2014 to July 2023 from the CGLS website (). Focus on our study area, monthly maps were processed to deploy time series and explore the relationship with corresponding LST products.

Downscaling Algorithm

The machine learning regression algorithm Random Forest (RF) (Breiman, 2001) was selected to establish the relationship between target variables LST and seven explanatory input kernels. The well-performed algorithm enables the capacity to aggregate multiple decision trees accompanied by random factors to form a nonlinear fitted function (f_RF). 1 $$\text{LST}={\mathrm{f}}_{\text{RF}}(\text{Elevation},\text{Slope},\text{Aspect},\text{Land}\,\text{cover},\text{Red},\text{NIR},\text{NDVI})$$ Initially, RF regression between the original LST ($${\text{LST}}_{1000}^{\text{ori}}$$) and explanatory data sets with the same 1,000 m resolution was performed to predict new LST ($${\text{LST}}_{1000}^{\text{RF}}$$). Residual is then calculated using the equation: 2 $${{\increment}\text{LST}}_{1000}={\text{LST}}_{1000}^{\text{ori}\,}-{\text{LST}}_{1000}^{\text{RF}}$$ The model is continuously applied to explanatory data sets (250 m) to deliver downscaled LST products and the $${{\increment}\text{LST}}_{1000}$$ is subsequently compensated as a residual correction to deliver the final gridded-based LST.

Daily downscaling was accomplished by sampling all available pixels within the data sets. To ensure adequate model diversity and prevent overfitting, each RF regression model was constructed with 500 decision trees (Belgiu & Drăgu, 2016).

Intercomparisons

Independent Satellite -Based LSTs

In the absence of a well-developed regional network for ground-based LST measurements, intercomparison with independent remotely sensed LST data sets offers a viable alternative (Hutengs & Vohland, 2016; Pu & Bonafoni, 2023; Yoo et al., 2022). Here, the newest mission Landsat 9 (LS9) and ASTER were employed to verify against diurnal and nocturnal downscaled LSTs, respectively. For each of the referenced satellite-based LSTs, five image pairs, ensuring representation of both dry and rainy seasons, were randomly selected within the last 5 years (i.e., 2019–2023) to perform the comprehensive intercomparison across different periods. Pixel quality bands of LS9 (QA_PIXEL) and ASTER (QA_DataPlane) were used for masking, as only clear-sky pixels were used for further evaluation.

Moreover, to diminish the considerable dissimilarity due to sensor characteristics and LST retrieval algorithms between products (Hutengs & Vohland, 2016; Yoo et al., 2022), just before the direct comparison, linear relationships between original MODIS and ASTER/LS9 LSTs with the same 1,000 m were performed, then the obtained slope and intercept of regression equation were applied to the ASTER/LS9 LSTs 250 m resolution.

Several widely known statistical metrics, namely Pearson correlation coefficient (R), Root Mean Square Error (RMSE), Normalized RMSE (NRMSE), and Mean Absolute Percentage Error (MAPE) are used to evaluate qualitatively the performance (formulations in Table S1 of Supporting Information S1).

Gauged Station Air Temperature Data

We then evaluated coarse and fine spatial resolution products against near-surface air temperature (T_air) collected in 8 different stations with available daily observation from 2000 to 2019. T_air is conventionally measured at approximately 1.5–2 m above the ground by standard meteorological stations. Whereas T_air and LST have distinct physical meanings and exhibit different responses to land cover characteristics and complex terrain (Lian et al., 2017), they maintain a strong correlation, primarily due to their shared dependence on energy exchange, surface attributes, and atmospheric conditions (Cristóbal et al., 2008; Kilibarda et al., 2014; Mutiibwa et al., 2015). Here, we investigate the linear relationship between these two variables to assess the impact of downscaling processes on different terrains and landscapes. Since T_air was provided by the aggregation over 24hr, mean values of diurnal and nocturnal LSTs of the same day were calculated to perform match-up exercises.

Reconstruction of Missing Pixels for Gap-Free Mapping Purpose

As located in a tropical monsoon climate region, cloud contamination is normally represented in almost time, especially in rainy seasons. Consequently, missing data will certainly appear in almost all spatial products due to the removal of cloud-affected pixels under quality control. Tackling the mentioned issues and obtaining gap-free data sets for further analysis, the Whittaker algorithm (Whittaker, 1922) was applied to all monthly downscaled LST data sets. Based on penalized least squares, the idea of the algorithm is to find the optimal smoothed time series while balancing the fidelity of the original time series and the roughness of the fitted time series. A full description of methods can be found in (Atzberger & Eilers, 2011; Eilers et al., 2017), with general equations as follows: 3 $$S=\sum\limits _{\mathrm{t}}{\left({y}_{\mathrm{t}}-{z}_{\mathrm{t}}\right)}^{2}$$ 4 $$R={\sum }_{\mathrm{t}}{\left(\left({z}_{\mathrm{t}}-{z}_{\mathrm{t}\mbox{-}1}\right)-\left({z}_{\mathrm{t}-1}-{z}_{\mathrm{t}-2}\right)\right)}^{2}={\sum }_{\mathrm{t}}{\left({z}_{\mathrm{t}}-2{z}_{\mathrm{t}-1}+{z}_{\mathrm{t}-2}\right)}^{2}$$ 5 $$Q=S+\lambda R$$ where y is the original time series, z is the smoothed time series, S is the fidelity measured by the sum of squares of the difference, R is the roughness of the smoothed data, Q is the different matrix, and λ is the smoothing parameter as determined by the user. Here, the objective is identifying the optimal value for λ that minimizes Q. The selection of λ often depends on the desired outcome: preserving the original trend in the raw data or achieving a smoother time series with potentially reduced high-frequency variations (Atzberger & Eilers, 2011; Eilers et al., 2017). In this present study, the appropriate λ is set to 0.1 in order to keep the original trends of the time series.

Time Series Analysis

Taking into account the alteration of dry and rainy seasons in LST/vegetation data sets, time series decomposition is used to separate quantitatively the observed values into precise trends, seasonal components, and residuals. Here, additive decomposition using locally estimated scatterplot smoothing (LOESS) function, is approached with the assumption of negligible variations in seasonal-cycle over examined periods (Cleveland et al., 1990). 6 $${y}_{\mathrm{t}}={S}_{\mathrm{t}}+{T}_{\mathrm{t}}+{R}_{\mathrm{t}}$$ Where $${\mathit{y}}_{\mathrm{t}}$$ represents the mean values of monthly products at time t, $${\mathit{S}}_{\mathrm{t}}$$ is the seasonal component (i.e., magnitude of dry and rainy seasons), $${\mathit{T}}_{\mathrm{t}}$$ is the trend of considered time series, and $${\mathit{R}}_{\mathrm{t}}$$ is the residual related to the anomalies in LST data sets due to the extreme events. In the present study, we applied the function to monthly time series of both LST in daytime, nighttime, all-day time (average of daytime and nighttime), and vegetation traits LAI data sets.

Results

Performance of Downscaling Models

Evaluation against referenced LS9 and ASTER LSTs strongly indicates that the daytime and nighttime downscaling model achieved high correlations through corresponding R of 0.944 and 0.857 aggregated for five selected days (individual day comparisons for both original 1,000 m and downscaled 250 m resolution products are provided in Tables S2 and S3 of Supporting Information S1). On average for both diurnal and nocturnal downscaling models, RMSE reaches 1.9°C, along with NRMSE and MAPE are 8.3% and 6.7%, respectively. The daytime model surpasses the nighttime model, as slightly better statistical metrics are presented in Figure 3.

[IMAGE OMITTED. SEE PDF]

Regarding the relationship with gauged station T_air, fine resolution LSTs show higher correlations compared to those of coarse resolution. This suggests that the downscaling process has improved the ability of LSTs to capture the spatial variability associated with land cover and terrain conditions. Separated by dry and rainy seasons, achieved R values demonstrate stronger correlations during dry seasons compared to rainy seasons. Moreover, the number of match-up observations (N) in rainy seasons is noticeably lower than in dry seasons (nearly 90%), implying the missing data due to cloud contamination, leading to the challenges in verifying the relationship between T_air and LST normally from May to October over the study regions (Figure 4).

[IMAGE OMITTED. SEE PDF]

Gap-Filling Process

All daily LST products from 2000 to 2023 were applied to quantify the percentage of missing data separating daytime versus nighttime, dry versus rainy seasons (Figure 5). Statistically, the missing pixel proportion during dry seasons calculated average for both diurnal and nocturnal data sets is approximately 66%, while the value for rainy data sets is nearly 90%, highlighting the dramatic expansion of clouds during rainy seasons. Results of the Whittaker algorithm showed anomalous peaks and troughs seem to be smoother going associated with time series while maintaining absolutely the original trend and seasonal dynamics in both data sets (Figure 6).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Mapping and Time Series Analysis

To generate products suitable for mapping and subsequent analyses, we calculated the mean monthly of daily downscaled LST products. These products facilitated the representation of typical monthly LST patterns, identification of the hottest and coldest months, and time series decomposition spanning a 24-year period. Maps of monthly composite 2000–2023 show the range of LSTs in daytime normally spanning from 15–45°C, while in nighttime varying from lowest 10 to highest 30°C depending on particular periods and locations (Figure 7). The consistent pattern of LST distribution in both diurnal and nocturnal products, as well as across monthly maps hint obviously the region of lowest LST located in the Southeast with major forest cover (see Figure 2), while highest LST belong to Northwestern region with dominated savannas.

[IMAGE OMITTED. SEE PDF]

Distributed histograms of LST during day and night for individual months depicted in Figure 8 denote obviously the hottest month is April, while the coldest month is recognized as December. The periodic transitions of mean diurnal LST increase from lowest (25°C) in December to highest (34.3°C) in April, while nocturnal LST is reasonably lower (17.4–22.1°C). The offset of LST between day and night in each month is dramatically different, with the largest change recorded in March, April, and May (on average of over 12°C), while decreases mainly in the wet season (just above 7°C).

[IMAGE OMITTED. SEE PDF]

Variations from daytime and nighttime, the hottest and coldest months, and the period 2000–2023 show different spatial LST patterns since different rates of heat absorbed based on various surface types (Figure 9). The Northwestern region shows significant changes, while the area in the Southeast gives inconsequential dissimilarity in both three considered aspects. The major density of histogram-based pixels tends to be larger than mean values when comparing between day–night and April–December, while there seems to be a normal distribution for the annual difference 2000–2023.

[IMAGE OMITTED. SEE PDF]

The outcomes of decomposition functions clearly show that the time series witnessed four breakpoints of trends over the observed period 2000–2023 (Figure 10). Through the estimated trend slope, four specific subperiod (i.e., 2000–2004, 2005–2009, 2010–2017, and 2018–2023) are determined with the confirmation for the same upward/downward trend in both daytime, nighttime, and all–day time of LST data sets. Overall, the trends grew 0.11°C in daytime, 0.87°C in nighttime, and 0.489°C for all-day time within 24 years, while seasonal amplitude in daytime (−3.7°C–4.8°C) is more fluctuated than in nighttime (−2.5°C−1.9°C).

[IMAGE OMITTED. SEE PDF]

Relationship Between LAI and LST

Monthly composite LAI maps obtained from PROBA-V/OLCI data sets (2014–2023) were performed in Figure 11. Similar to land cover (Figure 2), high LAI represented forest cover in the Southeast, while lower LAI spanning from the center to the Northwestern. Figure 12 illustrates the temporal patterns of downscaled LST (all-day time) compared to LAI by extracting the originally observed values into trend and seasonal amplitudes. Consistent time series authenticate solidly the response of LAI following LST variations. Through 10 years of observation, the lowest LAI is found in March or April (on average of 2.192 m².m⁻²) belonging to the hottest period annually, while the highest period LAI is generally October (on average of 4.102 m².m⁻²) during the middle of rainy seasons. The seasonal magnitudes also indicate LAI drops by approximately 1.038 m².m⁻² when LST rises by 1.976°C in the hottest month of April, in contrast to the increase of LAI by 0.925 m².m⁻² with LST decrease around 0.864°C during October.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Sampling all pixel-based land cover attributes of monthly composite products, Figure 13 shows the same negative correlation with R ranging from −0.717 to −0.45 between LAI and diurnal/nocturnal/all-day time downscaled LSTs. Correlation in nighttime is slightly lower than in daytime approximately 15%, while all-day time shows the closest relationship in all months. The highest correlation belongs to April (the hottest month in dry seasons) while this relationship reduces in the middle of rainy seasons, that is, July, August, and September. This suggests that the interplay between vegetation function and LSTs undergoes changes during these months, potentially due to the complex dynamics of vegetation in rainy periods.

[IMAGE OMITTED. SEE PDF]

Discussion

Obtaining the Gap-Filled Fine Spatial Resolution LST Data Sets

Our selection of MODIS data was driven by its daily revisit frequency and provision of both daytime and nighttime LSTs, while its coarse spatial resolution was resolved through a downscaling procedure. Intercomparisons of fine resolution LSTs against referenced LS9/ASTER source, along with a slightly higher correlation with T_air after downscaling, implying the impact of topography and land cover in coarse LSTs. These findings underscore the necessity of a downscaling framework, especially in complex landscape areas. It is also important to acknowledge the scaling effect, a common concern in downscaling processes (Zhou et al., 2016). This effect may be less significant in our case (downscaling from 1,000 to 250 m), while becoming more pronounced in high and very high resolution downscaling scenarios, such as those aiming for 20–30 m resolution (Pu, 2021).

The absence of a standardized approach for selecting explanatory variables in downscaling LST models necessitates area-specific evaluation of their importance for understanding the influence of these variables. In this study, the justifications for the selection of seven explanatory variables (Section 2.2) and their importance were quantified (Figure S1 in Supporting Information S1) as supporting our choice. While the performance of input kernels varies with employed downscaling models, it is revealed that the use of RF resulted in the highest accuracy while testing with different explanatory factors (Dong et al., 2020). Additionally, it is important to note that pre-processing resampling methods can be impacted by minor georeferencing inaccuracies, such as those arising from using products with different projections (e.g., MODIS and ALOS) or applying nearest neighbor resampling to the input kernels (Hutengs & Vohland, 2016; Tan et al., 2006).

Limitations inherent to downscaling algorithms RF prevent extrapolation of target values beyond the training data range, particularly when considering sampling techniques. To address this limitation and capture a wide range of LST values, we comprehensively sampled all available pixels within the downscaling data sets. Reliable prediction requires sufficient sample inputs during training. Consequently, model performance may be hindered in case of large proportions of missing pixels, as evidenced by our accuracy assessment based on reference data from LS9/ASTER and Tair, which revealed a noticeable decrease in performance during rainy seasons. This highlights the sensitivity of RF models to missing data. Additionally, exploring improved RF versions like local linear forest (Yoo et al., 2022) or adaptive RF (Ebrahimy et al., 2021) may be a potential alternative to overcome this limitation.

We noticed the major causative factor of cloud cover over tropical regions. Given the estimated proportion of missing pixels (Figure 5), we highlight the inevitable role of data reconstruction in obtaining finalized maps. While spatial data reconstruction has already been demonstrated as a valuable tool for satellite-based products with cloud contamination (Xie & Fan, 2021), the synergy of multi-sensors could also be a promising solution (Desai et al., 2021), yet requires tackling the asynchronous time of satellite overpass (Crosson et al., 2012), or integrating meteorological station data for spatiotemporal-explicit products (Zhao et al., 2020).

Understanding Mean of LST Over Spatial Pattern and Temporal Profile

By tracking the spatiotemporal LST data sets, we recognized the repeat of the hottest month (April) and the coldest month (December) both found in dry seasons. While April was detected as the middle dry period showcasing the highest LST, the reason for the lowest recorded LST in December is the effect of the winter monsoon formed by the meteorological phenomenon “Siberian High” (Mori et al., 2014; Nguyen et al., 2014).

The present study is bound by examining the mean LST for all land cover classes, without insight into individual land cover attributes. We aimed here to approach a general overview of LST dynamics over complex terrains and landscapes, accompanied by the effect of seasonal variations contributing to trends using long-term time series decomposition. From our empirical results, a preliminary statement is the mainland with major forest cover (Figure 2) plays an indispensable role in controlling the harmonic of LSTs between day and night, the hottest and coldest months. Future research should focus on the influence of specific land cover (e.g., forests, urban areas, and agriculture) to gain deeper insights into how different landscape types interact with LST patterns and seasonal trends.

We emphasized the effect of tropical monsoon climate over the examined area, resulting in significant LST changes during dry and rainy seasons. As a result, time series decomposition was applied to eliminate seasonal contribution and deliver precise trends (Quan et al., 2016). Hereby, results of decomposition function showed a warming trend for 2000–2023 with an average of 0.02°C year⁻¹ in all-day time, and the faster trend in nighttime (0.036°C year⁻¹) compared to daytime (0.005°C year⁻¹).

Interaction and Relationship of Vegetation Stand Traits With LST Over Spatiotemporal Dimensions

Among various biophysical vegetation variables, LAI was selected due to its critical role in representing plant health and vegetation cover (Fang et al., 2019). The spatiotemporal analysis allows us to establish a robust observational constraint for the inverse relationship between LAI and LST under tropical monsoon climates. As both LAI and LST belong to the domain of land-atmosphere in the ECV framework, understanding the underlying mechanisms of vegetation-related biophysical climate effects provides valuable insights for effective climate change mitigation and adaptation strategies (GCOS, 2018). However, our results of moderate to strong correlation coefficients (Figure 13) indicated that LST partially influences LAI and vice versa. Under tropical monsoon climates, other environmental factors such as atmospheric and soil conditions also play an indispensable role in the bidirectional interaction of vegetation cover and temperature.

Previous studies have mostly focused on spectral indexes to analyze the relationship between LST and vegetation (Akomolafe & Rosazlina, 2022; Ambika & Mishra, 2019; Sun & Kafatos, 2007). With the emergence of reliable methods to retrieve biophysical and biochemical variables using remote sensing data (Berger et al., 2022; Verrelst et al., 2019), we suggest future analysis to improve our comprehensive understanding of the interaction between leaf/canopy functional plant characteristics (i.e., chlorophyll, water, dry matter, carotenoid, anthocyanin) and LSTs. Importantly, the analysis should consider a fine-scale application within a complex area, highlighting the capability of the downscaled LST data to capture intricate spatial variations of vegetation cover.

Despite LST fluctuates daily, it is challenging to detect the LAI response through daily observation due to the slower dynamics of vegetation. In addition, monthly composite LSTs, which average out short-term anomalies, offer a more stable representation of temperature variations (Liu et al., 2023). Particularly, monthly MODIS LST product was recommended for climate studies (Chen et al., 2017) as well as elucidating systematic behavior of vegetation dynamics under LST variations interannually (Sun & Kafatos, 2007; Yu et al., 2023). In the present study, monthly calculations also mitigated the issue of extensive cloud-induced missing pixels in daily LSTs, especially during the rainy season, enabling a robust Whitter gap-filling process. We therefore adopted a monthly monitoring approach to facilitate seasonal analysis under tropical monsoon climate using a time series decomposition technique. Our findings indicated that appropriate LST promotes vegetation development in the middle of rainy seasons, generally in October. Conversely, higher LST leads to vegetation stress as LAI reduces significantly mostly in April during dry seasons.

Regarding the result of moderate to strong inverse correlations between LST and LAI, we noticed some discrepancies in previous studies, for example, a positive relationship while considering in continental scale (Rasul et al., 2020), non-correlation in heterogeneous forest landscapes (Neinavaz et al., 2019; Stobbelaar et al., 2022), or different interaction between diurnal and nocturnal LSTs (Abdullah et al., 2020). Notwithstanding, to our knowledge, these dissimilar relationships imply the variety between land cover attributes, terrain, and climate conditions, as well as distinct retrieval methods upon separated TIR sensors. Our findings stated here taking into account: (a) sampling all exhibited land cover classes, (b) using two independent data sets MODIS-based LST (250 m) and PROBA/OLCI-based LAI (300 m), and (c) confirmation through time series 2014–2023 within the tropical monsoon climate background. Challenges remain to be elucidated partially due to the limited number of studies, suggesting more future attempts toward trustful and unified authentications of the influence of climate factors on vegetation characteristics (Neinavaz et al., 2021).

Conclusion

Our attempt in the present study utilized different processing steps to obtain high-resolution and gap-free MODIS LSTs acquired in both daytime and nighttime. This enables meaningful and reasonable analysis to understand LST dynamics over the long-term period 2000–2023, as well as the relationship with vegetation stand traits over a tropical plateau. Messages achieved include: (a) Original images captured in TIR domain over tropical regions remain limitations in not only coarse spatial resolution but also deficiency data affected mostly by cloud contamination. Therefore, techniques including image sharpening and missing pixel reconstruction should be essentially implemented. Our experiments indicated that high accuracy downscaled products, as well as reliable gap-filling procedure, resulted in spatiotemporal-explicitly LST data sets for further analysis, (b) Through different experiments, we detected the hottest and coldest month with the spatial LST differences of day and night, dry and rainy season, and the period 2000–2023. Besides, time series decomposition indicates conspicuously the seasonal cycle and warming trend through the recent 24 years, and (c) Fine-scale application was demonstrated by spatial correlation analysis and time series monitoring between high-quality diurnal, nocturnal, and all-day time LSTs and the most crucial biophysical vegetation variable LAI, leading to the confirmation of negative relationship and interaction between LST and LAI dynamics. Moving toward deep insights into LST-related applications, we encourage the amalgamation of different TIR sensors to retrieve consistent LSTs, independent investigation of LST dynamics under individual land cover attributes, as well as empirical verification the quantitative relationship between LST and diverse biophysical and biochemical vegetation variables in both leaf and canopy level.

Acknowledgments

This study was funded by the project “Nghiên cứu tác động của yếu tố khí hậu tới các loại hình sử dụng đất nông nghiêp tại tỉnh Đắk Lắk” under the program of Dak Lak Department of Science and Technology (2022–2024). The authors also would like to thank for the supporting of project VAST01.06/22–23, and the project CSCL22.02/25–26 as funded by Vietnam Academy of Science and Technology (VAST).

Data Availability Statement

Data used for this study can be freely downloaded from public websites. MOD11A1.061 Terra Land Surface Temperature and Emissivity Daily Global 1 km from . MOD09GQ.061 Terra Surface Reflectance Daily Global 250 m from . ALOS DSM: Global 30 m v3.2 from . MCD12Q1.061 MODIS Land Cover Type Yearly Global 500 m from . Leaf Area Index from . USGS Landsat 9 Level 2, Collection 2, Tier 1 from . ASTER from . Data and code are available in .

References

Abdullah, H., Omar, D. k., Polat, N., Bilgili, A. V., & Sharef, S. H. (2020). A comparison between day and night land surface temperatures using acquired satellite thermal infrared data in a winter wheat field. Remote Sensing Applications: Society and Environment, 19, [eLocator: 100368]. [DOI: https://dx.doi.org/10.1016/j.rsase.2020.100368]