Spatial-Multitemporal Analysis of Heatwaves in Thailand: Discrepancies between In-Situ Air Temperature and Remote Sensing-Derived Land Surface Temperature

1. Introduction

Heatwaves, being acute climatic hazards, pose a persistent and severe threat to human well-being and the natural environment. These extreme events have extensive and harmful consequences across global, national, and local scales, impacting society, economies, and built environments. Lethal exposure to heat and heatwaves is associated with mortality and may amplify morbidities [1], resulting in particularly pronounced direct health effects among vulnerable populations, including children, the elderly, the sick, and low-income communities [1,2]. Additionally, indirect impacts encompass disturbances in ecosystem productivity [3], economic losses [4], livestock fatalities [5], increased fire risks [6], effects on water resources and agriculture [7], and power outages [8]. Between 1998 and 2013, extreme heat led to over 100,000 deaths in 164 locations across 36 countries, rising to 166,000 in 2017 [9]. The deadly heat event of 2019 claimed 5,000 lives during record-breaking temperatures in Europe, Asia, and Australia. In the summer of 2022, South Asia, North America, Europe, and China all experienced record-breaking heat, contributing to 15,000 heat-related deaths in Europe [10,11]. This trend persists into the following year, with Argentina, India, Vietnam, Laos, and Thailand experiencing record temperatures. In May 2023, twelve million people in the Pacific Northwest are under a heat advisory, and Canada faces forest fires [12]. Clearly, extreme heat conditions increasingly impact human health and healthcare systems, particularly by elevating global morbidity and mortality rates [13,14]. Studies reveal that the frequency of heatwaves changes rapidly worldwide, with trends in duration and cumulative heat accelerating since the 1950s [15]. Experts anticipate that these characteristics will become up to seven times more probable within the next 30 years [16,17]. According to the Eurostat mortality database, the number of heat-related summertime deaths may reach 68,000 by 2030, increasing to 94,000 by 2040 and 120,000 by 2050 [18]. Consequently, assessing heatwaves is a critical research area, especially in the context of climate change and the escalating global frequency and intensity of heatwaves.

In Thailand, studies demonstrate that the El Niño southern oscillation (ENSO) induces dry and warm conditions, as seen in 2015, resulting in the country’s second warmest year in 65 years and an escalation of extreme temperatures [19,20,21,22]. Consequently, heatwaves in Thailand led to 158 deaths between 2015 and 2018, particularly in the northern and central regions [23]. While extensive research on extreme warmth focuses on cities in high-income nations, such as Europe, the United States, and Australia [24,25,26,27,28], there is limited evidence and attention given to developing countries, especially Thailand [29]. As [30] shows the atmospheric scale interactions and spatiotemporal dynamics contribute to increased temperatures across the region for both urban and rural environments. In addition, comparative studies between rural and urban areas are sparse [31], despite urban areas generally experiencing more severe heatwaves, and rural areas undergoing rapid transformation into urban zones  [32,33]. To comprehend the intricate transition of urban-rural characteristics and enhance climate-informed decision-making, it is crucial to observe spatial patterns of heat events and identify hotspot areas [34,35]. Introducing a spatially focused framework and measures [16] is essential, particularly in regions like Thailand with high humidity levels that can amplify the impacts of high temperatures on human health and the environment.

The lack of a consistent definition for heatwaves persists despite continuous research efforts. Previous studies employ diverse methods, resulting in variations based on meteorological, socio-demographic, acclimatization attributes, and geographical conditions [36]. Traditional definitions, often reliant on absolute temperature thresholds and the number of affected days (typically two or more), may not suit tropical regions where temperature ranges are narrow, but humidity levels vary significantly [26,37,38,39,40]. An alternative definition considers heatwaves as extended periods of $T_{\mathrm{air}}$, measured at 1.5-2.0 m in height, exceeding a specific threshold for at least three consecutive days [15,41]. In Thailand, national agencies and researchers identify heatwaves when the daily maximum temperature exceeds the normal maximum by $5^{\circ }\mathrm{C}$ for over five consecutive days [23]. While others opt for specific values [42,43,44]. Some studies incorporate factors like $T_{\mathrm{air}}$ and relative humidity to assess heatwave impacts on human thermal comfort and stress. For example, [45] and [46] utilize the R.G. Steadman method to define the heat index (HI). However, the 90th percentile of daily maximum and minimum temperature, known as CTX90pct (daily maximum temperature) and CTN90pct (daily minimum temperature), is utilized by several researchers [41,47,48]. This framework offers consistent and universally applicable metrics for heatwave frequency, duration, and intensity, providing cutoff numbers that prove beneficial for diverse locations and impact sectors [49,50], especially in tropical areas like Thailand [51].

Satellite remote sensing technologies play a pivotal role in the comprehensive assessment of heatwaves, enabling global-scale monitoring and analysis in real-time and over extended periods. Unlike in-situ $T_{\mathrm{air}}$, satellite-derived LST serves as a crucial parameter for heatwave determination [52] and is a key factor in understanding the impacts of extreme heat [53]. It is noted that a $1^{\circ }\mathrm{C}$ difference in surface temperatures leads to a 4.5% increase in heat mortality risk [19,26,54]. Satellite-based LST data offer significant advantages, including high resolution and wide coverage, overcoming limitations posed by unavailable or poorly distributed ground station networks  [50,55]. [56] emphasize the benefits of utilizing the Moderate Resolution Imaging Spectroradiometer (MODIS) data for monitoring LST dynamics and trends, providing the longest consistent time series covering vast global regions. However, widespread data gaps in LST product retrieval due to cloud cover affecting over 60% of global MODIS LST datasets. Therefore, to overcome the restriction of missing values resulting from clouds, a range of research has focused on developing reconstruction methods [57,58,59,60,61,62,63], such as a random forest (RF) machine learning algorithm to estimate uncompleted LST data. While LST correlates highly with $T_{\mathrm{air}}$, their magnitudes and temporal behaviors exhibit substantial heterogeneity [64]. Limited exploration is dedicated to evaluating the relationship between LST thermal anomalies and heatwaves, including determining whether high LST values indicate heatwave occurrences [50]. Furthermore, it is crucial to examine how heatwaves detected at the near-surface and LST may be related, amplified, or mitigated across space and time.

Spatial and temporal heterogeneity plays a crucial role in the variations of LST between urban and non-urban areas [46,49,50,51,52]. [21], [56], [65], and [66] illustrate that the composition and configuration of land cover significantly influence LST magnitudes. The proportion of different landscape types emerges as the most influential factor affecting LST. [67] further notes that cities experience a slower warming process as impervious surfaces expand. They highlight that thermal contributions from suburban areas increase, and the rapid expansion of urban edges exacerbates local heat consequences. In Bangkok city, different land use categories derived from Landsat satellite and MODIS images reveal that high-density residential and commercial areas at the core exhibit mean UHI intensities ranging between $4^{\circ }\mathrm{C}$ and $6^{\circ }\mathrm{C}$ over the period of 2003–2016 [68,69] and up to $6.{20}^{\circ }\mathrm{C}$ in 2020 [70]. Linear regression analysis indicates that built-up land has a positive correlation with LST, with a 1% increase in built-up area resulting in a $0.{15}^{\circ }\mathrm{C}$ increase in LST in Bangkok city  [71]. [72] describe the spatial pattern of LST in Bangkok, identifying high-LST hot spots (>${40}^{\circ }\mathrm{C}$) in downtown, northern, and eastern parts. An up-to-date, localized understanding of the spatio-temporal patterns of LST as a proxy for heatwaves is lacking in this country, particularly in urban and non-urban areas. The existing literature on heatwave measurement underrepresents the implications of integrating earth observation data with other sources (e.g., ground-based measurements and model simulations).

1.1. Research Gap and Our Novel Contribution

This study addresses a significant gap arising from the alarming increase in heatwave problems and the absence of contemporary, longitudinal, area-based comparative studies in urban and non-urban areas of Thailand. To our knowledge, we are pioneers in employing a geospatial-based method, integrating remote sensing data with ground observations, to delineate spatiotemporal heatwave patterns. We systematically detect heatwave episodes and their associated metrics—number, frequency, duration, magnitude, and amplitude—across the entire year, identifying trends. In addressing data gap-filling, we crucially optimize the dataset during the preprocessing phase by systematically addressing missing values, ensuring the integrity and precision of climate-data-driven predictions. Specifically, for gap-filling in spatial and temporal LST, we employ a random forest (RF) machine learning algorithm. This technique is known for its robust predictive capacity and versatility with nonlinear data, making it effective in predicting a range of extreme climate events. The model integrates both temporal and spatial variables, offering a comprehensive evaluation of feature importance, particularly for the predicted LST data. Furthermore, we explore the amplification and attenuation effects through spatial and temporal homogeneity in correlation analysis. It has been demonstrated how this highly accurate technique can serve as a valuable tool in supporting sustainable climate practices, government policy, and decision-making processes.

Our results not only advance the utilization of satellite data in heatwave assessment, especially for residents situated far from meteorological stations, but also furnish policymakers and the general public with a nuanced comprehension of specific heatwave events and trends at the local level. Crucially, our contributions bridge a significant gap in adaptable measurement and spatial modeling of heatwave patterns and distribution across diverse dimensions in public decision-making. This work showcases substantial advancements in heatwave assessment in Thailand, enabling more effective climate-resilient spatial planning to mitigate the rise in heat-related illnesses and deaths and facilitating adaptation to the escalating magnitude of impending extreme climatic crises. Importantly, our work introduces a pioneering approach, emphasizing its novel contribution to the field in Thailand.

2. Materials and Methods

2.1. Study area description

In this study, our focus is on urban, peri-urban, and rural areas in the central region of Thailand, as illustrated in Figure 1. We categorize the types of areas based on national standards, considering administrative, demographic, prominent land use, and the majority of occupation [73,74]. These areas align with the tropical monsoon climate, dividing into three seasons: summer (March to April), rainy (May to October), and winter (November to February). The selected areas are as follows:

•

Bangkok; an urban area located at latitude ${13}^{\circ }$38’N and longitude ${100}^{\circ }$54’E. This area is flat and low-lying, with altitudes ranging from 1.50 m to 2.0 m. It occupies $1,600k{m}^2$ and consists of 50 districts with a registered population of 5,527,994, known as the most highly inhabited city in Thailand. Built-up area covers 67.36% while agricultural land shows 0.16% of total area. The annual mean temperature ranges from $28.0^{\circ }\mathrm{C}$ to ${30}^{\circ }\mathrm{C}$, the average maximum temperature in April is $33.0^{\circ }\mathrm{C}$ to $38.0^{\circ }\mathrm{C}$, and the average minimum in January is $16.0^{\circ }\mathrm{C}$ to $25.0^{\circ }\mathrm{C}$. The annual average rainfall is 1,700 mm, the average relative humidity is 73%, and the average wind speed is 12 km/h.

•

Pathum Thani; a peri-urban area is located at latitude ${14}^{\circ }$01’N and longitude ${100}^{\circ }$32’E with an average altitude of 2.30 m. It occupies $1,519k{m}^2$ and consists of 7 districts with a registered population of 1,190,060. The annual mean temperature ranges from ${28}^{\circ }\mathrm{C}$ to ${30}^{\circ }\mathrm{C}$, the average maximum temperature in April is ${32}^{\circ }\mathrm{C}$ to ${34}^{\circ }\mathrm{C}$, and the average minimum in December is $24.0^{\circ }\mathrm{C}$ to $26.0^{\circ }\mathrm{C}$. The annual average rainfall is 1,200-1,400 mm, the average relative humidity is 80%, and the average wind speed is 16-21 km/h. This setting is characterized by 29.79% of the land serving as city outskirts to accommodate dense dwellings, industrial estates, and no forest areas.

•

Saraburi; a rural area located at latitude ${14}^{\circ }$31’N and longitude ${100}^{\circ }$54’E, with altitudes ranging from 2.0 m to 10.0 m. This area occupies $3,576k{m}^2$, and 643,963 registered residents live in 13 districts of this countryside. The annual mean temperature is $28.2^{\circ }\mathrm{C}$, with the average maximum temperature in April being $31.4^{\circ }\mathrm{C}$ and the average minimum temperature in January being $23.6^{\circ }\mathrm{C}$. The annual average rainfall is 1,185 mm, the average relative humidity is 79.9%, and the average wind speed is 1.49 km/h. Approximately 78% of this province is occupied by agricultural land and forest areas.

2.2. Data Sources

2.2.1. Ground-observed air temperature data

Daily air temperature data, including maximum ($T_{\max }$) and minimum ($T_{\min }$), spanning the years 1981 to 2019, were acquired from the automatic meteorological stations operated by the Thai Meteorological Department. These observations are measured in line with the World Meteorological Organization (WMO) standard at 2.0 m above the ground [75]. In this study, 10 meteorological stations provide in-situ temperature data as shown in Figure 1a. The weather stations are categorized based on their locations into urban, peri-urban, and rural groups, as shown in Table 1. Since Sara Buri province, representing the rural area, has no observed ground station, this study appoints the meteorological stations in Ayutthaya, Lopburi, and Nakhon Ratchasima to represent the nearest neighbor weather stations for this inactive site.

2.2.2. Satellite-observed land surface temperature data

Typically, satellite-observed LST is a widely used parameter for assessing heatwaves, as it provides wider spatial coverage compared to SAT measured at weather stations. Therefore, we rely on LST data sourced from MODIS sensor. This sensor is part of the Earth Observing System (EOS) operated by the National Aeronautics and Space Administration (NASA) and is designed to provide comprehensive data that aids in understanding global dynamics and processes. Our analysis incorporates two specific MODIS-LST products, both from tile h27v07: the Terra daily land surface temperature (MOD11A1) for the period 2000-2019, and the Aqua daily land surface temperature (MYD11A) for 2002-2019. We choose these periods to give a comprehensive overview of temperature trends over nearly two decades. We utilize the application for extracting and exploring analysis ready samples, namely A$\rho$$\rho$EEARS, which is available at https://appeears.earthdatacloud.nasa.gov/, to generate the MODIS-LST as a point extraction and supplementary data, taking advantage of its ability to provide efficient access to geospatial data.

2.2.3. Elevation and land use

The Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM), publicly available on the USGS Earth Explorer, has a spatial resolution of 90 meters. The elevation of the study area is determined using the SRTM DEM with a spatial resolution of 90 m. Subsequently, the data at 90-meter resolution are averaged to achieve a 30-meter resolution, as illustrated in Figure 1b. Additionally, land use/land cover data from the Land Development Department (LDD) of Thailand are utilized. These data, in raster format for the years 2000–2019, cover specific regions including Bangkok, Pathum Thani, Samut Prakan, Ayutthaya, Lopburi, and Nakhon Ratchasima. The data are categorized according to the Level 1 LDD standard classification criteria, comprising five classes: built-up area (U), agricultural land (A), forest area (F), water body (W), and miscellaneous land use (M), as depicted in Figure 1c. This classification is instrumental in understanding regional land use patterns and their potential influence on local temperatures.

2.3. Methods

This section outlines the methodology used to examine the spatiotemporal quantification of heatwaves and to identify associated trends in different regions of Thailand. Initially, the approach encompasses data collection from diverse sources. Then, it addresses missing data through imputation due to technical malfunctions. Subsequently, the methodology involves identifying heatwave metrics. Following this, it conducts trend analysis. Finally, it visualizes the spatiotemporal pattern of heatwaves, thus presenting a comprehensive processing sequence, as illustrated in Figure 2.

2.3.1. Data gap-filling

In the preprocessing phase, addressing the challenge of missing values is vital to ensure the integrity and precision of climate-data-driven predictions. The dataset undergoes optimization through the systematic removal or statistical imputation of these missing values, rendering it more suitable for subsequent predictive modeling or cluster analysis functions.

• 2.3.1.1 Spatiotemporal ground-observed air temperature gap-filling

Due to malfunctions at the meteorological station, the $T_{\mathrm{air}}$ dataset is identified as partially incomplete. Recognizing that the missing data follows a pattern of being missing at random (MAR), [76] recommends the use of simple imputation approaches to address the gaps. Depending on the extent of the missing data—specifically, when the missing data level is below 5 percent—a mean substitution approach is employed. This involves replacing a missing value with the mean of the observed values computed across several years.

Conversely, when the percentage of missing data exceeds 5 % or spans over a month, the regression model is applied. Regression models, whether utilizing a single closest neighbor station [77] or multiple nearby stations [78,79,80], have demonstrated effectiveness in estimating daily weather observations [78,81,82,83], and the imputation processes are applied to both $T_{\max }$ and $T_{\min }$.

• 2.3.1.2 Spatiotemporal satellite-based land surface temperature gap-filling

Gaps in the daily LST data retrieved from MODIS platforms are noted. Recently, some of the related studies focused on reconstructing the LST of satellite datasets, forecasting daily LST from time series data, and fusing multisource data to estimate subpixel LST ( [84,85,86,87]. Although machine learning methods have been used in LST retrieval, they have not yet been widely used [88]. To address this issue, the RF machine learning algorithm [89] is implemented, known for its robust predictive capacity and versatility with nonlinear data. Machine learning has been employed in climate and weather forecast sciences for decades, serving various goals such as post-processing, data assimilation, and physical analysis [90]. Moreover, machine learning models based on algorithms can examine and derive insights from patterns without following explicit instructions  [91]. Performance has been shown to improve with increased correlation between target and reference variables [92,93], making it widely used in analyzing high-dimensional relationships [94] and showing promise in predicting a range of extreme climate events [91]. In our application of RF, the predicted value is the class receiving the most votes across multiple trees. This method reduces variance and controls overfitting. For simplicity, we assume that values are missing completely at random, and missing values are imputed using sequential and non-sequential imputations [93]. For data infilling, we utilize the Python-based scikit-learn module [95]. We intentionally omit the scalar function, often used to standardize and normalize numerical inputs for classification and regression, to minimize potential bias.

• 2.3.1.3 Model performance and accuracy assessment of predicting land surface temperature

A comprehensive set of statistical measures is computed on the test set to assess the accuracy and performance of the predictive model for LST. These measures include important metrics like the coefficient of determination ($R^2$), root mean square error (RMSE), minimum and maximum confidence intervals, mean absolute error (MAE), mean bias error (MBE), and the confidence level is set at 0.95.

For the assessment, the calibration phase utilizes variable data spanning from 2000 to 2017, while the subsequent validation phase involves data from 2017 to 2019, following an 80% train and 20% test split. The $R^2$ value, ranging from 0 to 1.0, provides a valuable indication of the predictive model’s accuracy in estimating outcomes. The RMSE is effectively utilized to evaluate biases in both mean and spatial variance, whereas the MAE serves as a reliable measure of error magnitude, with lower values demonstrating superior performance. Additionally, the MBE offers insightful observations regarding the direction of error bias, with a value of zero denoting an unbiased estimation by the model [96,97].

2.3.2. Heatwave definition and its metric detection

In this study, the detection and measurement of $T_{\mathrm{air}}$ heatwaves and LST heatwaves occur using established methods referenced in [47,98], and summarized in Table 2. The study employs a criterion based on the calendar-day 90th temperature percentile to identify heatwaves. The CTX90pct method is utilized for $T_{\max }$ and LST day (MOD11A1 day, MYD11A1 day), while the CTN90pct method is employed for $T_{\min }$ and LST night (MOD11A1 night, MYD11A1 night).

To accommodate the variability of extreme heat events throughout the year, the 90th percentile threshold calculates for each specific day within the study period for the corresponding data type, resulting in a total of 365 values for each day. A heatwave is considered present when temperatures exceed the 90th percentile threshold for a minimum of three consecutive days. This approach ensures the capture of sustained extreme temperature events and accounts for temporal dependencies. The finding of [99] suggests that a relative threshold is more effective in capturing the regional nuances of heatwaves, especially in Southeast Asia’s complex landscape of land-sea distribution and varying topographies. Therefore, it is advisable to employ a percentile-based method for defining climate extremes like heatwaves in this region.

Additionally, the use of this method facilitates the assessment of prolonged heatwave conditions, enabling the identification of persistent extreme temperatures over an extended period. By utilizing a calendar-day estimate of the 90th temperature percentile, heatwave days can be identified across all seasons, providing a comprehensive understanding of heatwave occurrence throughout the year. Moreover, we measure the median start/end dates of the first/last $T_{\mathrm{air}}$ and LST heatwave incidents at all stations and areas for each year, listing specified stations as representatives of each different type of area.

2.3.3. Ground-observed air temperature and satellite-based land surface temperature heatwave metric trend analysis

The nonparametric tests have been widely used in climatic time series analysis due to their simplicity and suitability for data with outliers [100]. Therefore, we employ the non-parametric Mann-Kendall (MK) test to assess the statistical significance of trends in the series of yearly time series heatwave indices for both $T_{\mathrm{air}}$ and LST datasets. This analysis enables us to identify and quantify any significant trends observed in the heatwave patterns over time. To determine the magnitude of these trends, we utilize Theil-Sen’s slope, which provides a robust estimation of the slope in the presence of outliers and non-normal distribution.

The results are statistically significant if the p-values from the MK test are less than the significance threshold of 0.05. This threshold ensures that only meaningful trends are detected in the data, considering a 95% confidence level. To perform the MK test, we employ the Python package called pyMannKendall [101], which provides efficient and reliable calculations of the test statistic. This package allows us to handle large datasets and efficiently compute the MK test for $T_{\mathrm{air}}$ and LST heatwave time series. By conducting the trend analysis using the MK test and Theil-Sen’s slope, we explore the presence and magnitude of heatwave indices trends. This analysis provides valuable insights into the long-term changes in heatwave patterns and their significance in the study areas.

2.3.4. Spatial homogeneity in correlation analysis

The point-to-pixel analysis method is employed to match the series of heatwave metrics data, calculated based on $T_{\mathrm{air}}$ , with the corresponding pixels of gridded LST data. Only pixels that have at least one available and used a ground-based gauge for the calculation are included [102,103]. Given the extensive research area and the time required for retrieving daily MODIS data, a sampling design with a buffer zone at the confluence of the regular 5 km × 5 km latitude and longitude grid is adopted for validating $T_{\mathrm{air}}$ and LST heatwaves. The use of available grids is because heatwaves may not occur at every grid, considering the restriction of three consecutive days exceeding the threshold temperature in the heatwave definition [99]. Therefore, this produces a total of 248 valid grids encompassing three provinces of study.

To establish the spatial correlation between the pixel values of $T_{\mathrm{air}}$ and the corresponding pixels in the detected heatwave from the LST dataset, the pixel-wise correlation is determined using the Pearson correlation coefficient (r) throughout the research period. Furthermore, the spatial-temporal relationship between a meteorological station grid and the four nearest LST grids is determined using the median value of the four Pearson’s r coefficients. Notably, the$T_{\max }$ is cross-validated with daytime LST (MOD11A1 (day) and MYD11A1 (day)), whereas the $T_{\min }$ is evaluated with nighttime LST (MOD11A1 (night) and MYD11A1 (night)).

3. Results

3.1. Predicting land surface temperature to fill in missing satellite data and variable importance

Our RF models to predict LST missing data include both temporal (day of the year (DOY) and year) indicators and spatial variables ($T_{\max }$, $T_{\min }$, elevation, built-up land, agricultural land, forest land, water bodies, and miscellaneous land), allowing for an evaluation of feature importance specifically for predicted LST. To operationalize this method, we grid the area at 1-km intervals near each of the 10 meteorological stations and select 4 points according to cardinal directions to generate the predicted LST daytime and nighttime models. After achieving the models for LST prediction, we grid the area at 5-km intervals, covering all seven provinces with weather stations, resulting in 11 training samples and 1,250 testing samples. Of these, only 248 points overlap with the three targeted study area (Bangkok, Pathum Thani, and Saraburi).

According to RF results, a concise overview of the model’s feature importance for LST modeling is presented in Figure 3. The factor DOY show the highest explanatory rates for daytime LST, as in MOD11A1_day (Figure 3a) and MYD11A1_day (Figure 3b), whereas the daily minimum temperature ($T_{\min }$) affects most LST night, as in MOD11A1_night (Figure 3c) and MYD11A1_night (Figure 3d). In predicting LST day, $T_{\max }$ and $T_{\min }$ rank second, followed by year, elevation, and land use. For MOD11A1 data, the water body from the land use factor is essential, while the built-up area is vital for MYD11A1 data. Later, in the case of MOD11A1 and MYD11A1, $T_{\min }$ emerges as the second most influential factor for predicting nighttime LST, following the elevation of the land and the extent of built-up areas.

A summary of LST model statistical measurement between observed LST and predicted LST is given in Table 3. The result demonstrates that MOD11A1 night has the greatest satisfied RMSE (2.09 °C) and the best linear relationship ($R^2=0.64$), followed by MYD11A1 night ($R^2=0.48$), and MOD11A1 day ($R^2=0.45$). On the other hand, MYD11A1 day is an unsatisfactory model due to the highest RMSE (5.02 °C) and lowest $R^2$ ($0.29$). These results confirm the ability of the RF machine learning algorithm; particularly the MOD11A1 night model, to predict LST, whereas the other predicted LST models indicate overestimated data.

3.2. Heatwave detection and measurement

3.2.1. Ground-observed air temperature heatwave

$T_{\max }$ heatwave episodes in the metropolitan area mostly occur in April, terminate in July in 1981-1999, move to mid-March, and continue until late November in 2000-2019. Similarly, the rural heatwave starts early in January and extends until April, with the previous occurrence happening in December and reducing to October. The first heatwave, on the other hand, is observed at the rural-urban boundary in March and lasts until November. During the years 1981-1999, the first occurrence of $T_{\min }$ heatwave series at a city site is identified in May and finishes in July, while from 2000 to 2019, it begins in February and ends in late December. The outskirts start in February with the first event and end in November. The rural region starts in late April and moves to February, with the conclusion date in June and pushed back to September. Notably, the start date of both $T_{\max }$ and $T_{\min }$ heatwaves tends to be earlier, whereas the end date is delayed in all places.

Figure 4 summarizes year-by-year and cumulative detected day heatwave results by $T_{\max }$ (Figure 4a, c, e, g, and i) and the spatial distribution of cumulative heatwave indices (Figure 4b, d, f, h, and j). As a result, the highest HWN is 24 events in PPTN (Pathum Thani) in 2019. Following this, UBKD (Don Muang Airport) has the highest total HWF, with 162 days in 1997. For the total annual longest HWD, UBKD (Don Muang Airport) and RAYT (Ayutthaya Province) exhibit the greatest number, at 28 in 1997 and 2016, respectively. The highest mean heatwave temperature (HWM) is noticed in UBKK (Khlong Toei, Bangkok) at $45.5^{\circ }\mathrm{C}$ in 1997. Lastly, the highest HWA is recorded in RBCL (Bua Chum Meteorological Station, Lop Buri) at $43.2^{\circ }\mathrm{C}$ in 2016. Notable exceeded temperatures in 1997 and 2016 contribute to the highest heatwave number and its metrics in all locations.

On the other hand, Figure 5 represents night heatwave detection by $T_{\min }$ (Figure 5a, c, e, g, and i) and the spatial distribution of cumulative heatwave indices (Figure 5b, d, f, h, and j) to indicate a heatwave reporting disparate results. The maximum number of individual heatwaves is recognized in a peri-urban area in PPTN (Pathum Thani) at 24 events in 2019. Considering the fact that urban areas are responsible for the highest values of other indices, the HWF is attained 89 times in UTMD (Bang Na, Bangkok) in 2019, the HWD reaches 59 days in UBKD (Bang Na, Bangkok) in 2013, the HWM surpasses $33.3^{\circ }\mathrm{C}$ in UBKP (Bangkok Port, Bangkok) in 2015, and the HWA exceeds $30.8^{\circ }\mathrm{C}$ in UBKD (Don Muang Airport, Bangkok) in 2013 and 2016. Here, the years 2013 and 2019 are considered to be the warmest years and induce the largest number of heatwave events.

The difference in meteorological data intervals affects the impartiality of all data comparisons, and here, we present the analysis as a percentage determined based on the detected year’s ratio to the observed year, as demonstrated in Table 4. During the observed heatwave by $T_{\max }$, PSVN (Samut Prakan) serves as a peri-urban domain, examining the highest number of heatwave cases at 100%, followed by two sites in an urban area, UBKD (Don Muang Airport, Bangkok) and UBKP (Bangkok Port, Bangkok), at 95% and 92%, respectively. The lowest phenomenon occurs in UBKK (Khlong Toei, Bangkok) and UBCL (Bua Chum, Lopburi) at 79%. Meanwhile, in the observed heatwave by $T_{\min }$, PSVN (Samut Prakan) also remains at 100%, and after that, RPCN (Pak Chong, Nakhon Ratchasima), a rural area, is 87%, and PPTN (Pathum Thani), a peri-urban plat, is 86%. UTMD (Bang Na, Bangkok) has the lowest occurrence (62%).

In addition, Table 4 also displays the annual mean number of heatwave indices for both $T_{\max }$ and $T_{\min }$ heatwaves. We find that the average number of HWN trends for daytime $T_{\mathrm{air}}$ heatwaves is 3.4-4.2 days/year in the urban area, 5.2-6.8 days/year in the peri-urban area, and 3.6-4.0 days/year in the rural area. The highest HWF occurs 15.0-16.7 days per year in urban regions, 26.4-27.1 days per year in peri-urban areas, and 16.9-21.0 days per year in rural areas. The longest HWD lasts 4.5-5.0 days per year in cities, 5.5-6.0 days in peri-urban areas, and 5.9-6.5 days in rural areas. The warmest HWM in the urban region is $35.8^{\circ }\mathrm{C}$-$36.7^{\circ }\mathrm{C}$, $35.5^{\circ }\mathrm{C}$-$36.8^{\circ }\mathrm{C}$ in the peri-urban area, and $34.3^{\circ }\mathrm{C}$-$37.6^{\circ }\mathrm{C}$ in the rural area. After all, the hottest day of HWA is $37.5^{\circ }\mathrm{C}$-$37.9^{\circ }\mathrm{C}$ in the city, $37.6^{\circ }\mathrm{C}$-$38.1^{\circ }\mathrm{C}$ in the suburbs, and $36.4^{\circ }\mathrm{C}$-$39.7^{\circ }\mathrm{C}$ in the countryside. The highest incidence of day heatwaves is recorded in the peri-urban area of Pathum Thani. Urban areas, particularly Bangkok, consistently show the highest frequency and duration of heatwaves, with these trends intensifying over time.

In terms of nighttime $T_{\mathrm{air}}$ heatwaves, the HWN in the urban region is 3.0-3.5 days/year, similar to daytime for the other two types of places. The maximum HWF occur 12.3-15.5 days per year in the urban area, 18.8-22.5 days per year in the peri-urban area, and 3.6-14.3 days per year in the rural area. The longest HWD lasts 3.3-13.5 days per year in the city, 4.6-5.4 days in the peri-urban area, and 3.8-5.0 days in the rural area. The average cumulative magnitude (HWM) in urban areas is $35.8^{\circ }\mathrm{C}$-$36.7^{\circ }\mathrm{C}$, $26.3^{\circ }\mathrm{C}$-$27.1^{\circ }\mathrm{C}$ in peri-urban areas, and $23.7^{\circ }\mathrm{C}$-$26.3^{\circ }\mathrm{C}$ in rural areas. After all, the hottest day of the entire hottest HWA is $28.5^{\circ }\mathrm{C}$-$29.3^{\circ }\mathrm{C}$ in the city, $27.8^{\circ }\mathrm{C}$-$28.8^{\circ }\mathrm{C}$ in the suburbs, and $25.6^{\circ }\mathrm{C}$-$27.4^{\circ }\mathrm{C}$ in the countryside. Nighttime $T_{\mathrm{air}}$ heatwaves also follow a similar pattern as daytime $T_{\mathrm{air}}$ heatwaves, with urban regions experiencing more frequent and intense events.

3.2.2. Satellite-based land surface temperature heatwave

The analysis of data from 2000 to 2019 indicates that daytime heatwaves in urban and peri-urban areas commonly commence in either February or April and conclude in either September or November. In rural areas, these heatwaves commence as early as January and persist until late September. Urban nighttime heatwaves typically occur from February to July, while peri-urban nighttime heatwaves span from mid-February to mid-October or November. Rural regions experience nighttime heatwaves starting in late February, lasting through April, and concluding in November.

The LST dataset, annotating 1,250 patches, promises to spatialize spatiotemporal heatwave coverage in the extended area of the meteorological station. However, only 268 valid points covering the study domains are assigned to detect heatwaves. The results of each heatwave index are slightly larger than those detected by $T_{\mathrm{air}}$. For daytime heatwaves, observed MOD11A1 (day) over the period 2000–2019 indicates that the highest total number of HWN is 23 events in 2015 (Point 262, Ayutthaya) and 2019 (Point 168, Pathum Thani), respectively. The total number of HWF is highly significant in the rural area (Point 1209, Nakhon Ratchasima) with 183 days. Similarly, the rural area (Point 1178, Nakhon Ratchasima) has the longest total length of HWD at 59 days. The largest HWM occurs at Point 531 RBCL (Lopburi) in 2002 with $41.4^{\circ }\mathrm{C}$. The hottest day of HWA resides in the rural area (Point 722, Lopburi) at $44.4^{\circ }\mathrm{C}$. In addition, MYD11A1 (day) from 2002–2019 exhibits a tremendous cumulative amount. The peri-urban areas (Point 104 and 128, Pathum Thani) manifest the largest HWN with 27 events in 2019, while the largest HWD and HWA exist in the rural areas, 67 days in 2005 (Point 908, Nakhon Ratchasima) and $47.6^{\circ }\mathrm{C}$ in 2016 (Point 678, Nakhon Ratchasima), respectively. The warmest HWM is also at Point 604 RBCL (Lopburi) in 2017 with $45.5^{\circ }\mathrm{C}$. However, HWF is lower in MYD11A1 than in MOD11A1, with 180 days in 2005 (Point 1090, Nakhon Ratchasima).

The nighttime heatwave results surpass the daytime outcomes for all variables. In the observed MOD11A1 data, peri-urban and rural areas exhibit the highest HWN, with 29 incidents in 2010 at Point 73 in Samut Prakan, and in 2019 at Point 493 in Lopburi, respectively. Point 83 in Samut Prakan recorded the maximum HWF in 2010, with 190 days. The rural area at Point 237 in Ayutthaya experienced the longest HWD of 29 days in 2005. Meanwhile, the metropolitan region of Bangkok, at Points 20 and 40, recorded the highest HWM in 2013 and the warmest HWA in 2003, with temperatures of $28.3^{\circ }\mathrm{C}$ and $29.5^{\circ }\mathrm{C}$, respectively. In fact, the nighttime heatwave corresponds to MYD11A1; the derived metrics determine a number greater than MOD-LST, resulting in MYD-LST. HWN has 32 events in 2016 (Point 1195, Nakhon Ratchasima), whereas HWF has 212 days in 2005 (Point 211, Ayutthaya and Point 1138, Nakhon Ratchasima). HWD lasts 84 days in 2005 (Point 1138, Nakhon Ratchasima). In comparison, HWM at Point 1 UBKP (Khlong Toei, Bangkok) indicates the hottest magnitude in 2016 with $27.3^{\circ }\mathrm{C}$, and HWA in urban areas (Point 7 and 62, Bangkok) and peri-urban sites are warmest at $28.4^{\circ }\mathrm{C}$ in 2016 (Point 80, Samut Prakan).

Figure 6 summarizes the yearly distribution and characteristics of daytime heatwaves. Pathum Thani and eastern Bangkok have the highest number of heatwave events (HWN), while northern Pathum Thani experiences the highest frequency (HWF). The rural regions in northern Saraburi have the longest heatwave duration (HWD), while urban areas in Don Muang, Bangkok, exhibit the highest maximum temperature (HWM) and hottest day (HWA) values. Heatwaves are more frequent in peri-urban areas, while rural areas experience longer-lasting heatwaves.

Figure 7 shows the yearly nighttime heatwave indices. Downtown Bangkok and northern Pathum Thani have the highest number of nighttime heatwave events (HWN) and frequency (HWF). The longest heatwave duration (HWD) is observed in downtown Bangkok. The nighttime maximum temperature (HWM) and hottest nighttime temperature (HWA) are also highest in downtown Bangkok. Nighttime heatwaves are more frequent and longer-lasting than daytime heatwaves, but have lower maximum temperatures. Urban areas, particularly downtown Bangkok, experience severe nighttime heatwaves.

3.3. Trend analysis for temporal heatwave metric detection

3.3.1. Ground-observed air temperature trend

To evaluate the trend over the study period, we use the nonparametric Sen’s slope estimator with statistical significance at $\alpha$ = 0.05. The results show that the peri-urban station (PPTN, Pathum Thani) exhibits the most prominent increasing trends for $T_{\max }$ on HWN at 0.45 days/year, HWF at 2.00 days/year, and HWD at 0.27 days/year, as shown in (Table 5), which presents annual climatological averaged trends of each heatwave aspect. It is clearly followed by two sites in the urban area: UBKK (Khlong Toei, Bangkok) and UTMD (Bang Na, Bangkok), with 0.20 days/year for HWN, 0.75 days/year for HWF, and 0.10 days/year for HWD, respectively. Intriguingly, none of the three regions evinces a significant trend in HWM and HWA.

The finding provides results for the heatwave nighttime and metrics obtained by $T_{\min }$ (Table 5), and it demonstrates that the most significant increasing HWN trends are indicated in the rural area at 0.39 days/year in RLBR (Lop Buri), followed by 0.25 days/year in the urban area (UTMD: Bang Na, Bangkok), and 0.24 days/year in the rural area (RPCN: Pak Chong, Nakhon Ratchasima). Also clearly validated in the rural regions is the highest HWF, which is 1.44 days/year in RLBR (Lop Buri), 1.00 days/year in RPCN (Pak Chong, Nakhon Ratchasima), and UTMD (Bang Na, Bangkok). The most extended HWD demonstrates that RPCN (Pak Chong, Nakhon Ratchasima) has the most significant rising trend with 0.15 days/year, followed by RLBR (Lop Buri) with 0.14 days/year and UTMD (Bang Na, Bangkok) with 0.08 days/year, respectively. It is interesting to note that UTMD (Bang Na, Bangkok) is proven to be the only one with a trend towards an increasing HWA with $0.{07}^{\circ }\mathrm{C}$/year. However, there is no discernible pattern in HWM across any locations.

3.3.2. Satellite-based land surface temperature trend

The trend analysis of annual average heatwave metrics, conducted using the MK test combined with Sen’s slope (Table 6), indicates the existence of both positive and negative trends at the 95% level of significance. Specifically, only 2% of LST grids exhibit a statistically significant decrease, while 20-30% of grids show a significant increase across all daytime heatwave indicators. In the case of nighttime heatwaves, approximately 41–57% of the data demonstrate a significant upward trend, with only 1-4% of the data suggesting a decreasing tendency for all heatwave parameters.

The significant increasing trends in annual daytime heatwave detection range between 0.20 and 0.83 events/year for HWN, with the maximum number being found in Pathum Thani as of the suburban area, as seen in Figure 8b. The statistically rising trend of HWF varies between 0.91 and 3.08 days per year, and the region in Khlong Toei (Bangkok) (Figure 8d). Within the range of 0.19–0.57 days per year, the outskirt location in northern Pathum Thani exhibits a mean annual dramatic rise in the most extended length of HWD (Figure 8f). The upsurge of the mean HWM is between 0.10 and $2.{05}^{\circ }\mathrm{C}$/year, and the positive trend of the mean HWA ranges between 0.21 and $2.{21}^{\circ }\mathrm{C}$/year. Both trends indicate an upward tendency in temperature. Both degrees and values are at their highest in the Chom Thong western neighborhood of Bangkok (Figure 8g and i). This demonstrates that yearly daytime heatwaves climb, and the peri-urban areas are the ones in which this shift is most prominent.

Regarding the annual trends of nighttime heatwaves, the most substantial increase in the spatial variation of HWN varies from 0.14 to 1.11 events/year, and the range for HWD is from 0.16 to 0.94 days per year. As depicted in Figure 9b and Figure 9f, Bangkok’s northeastern Pravet district exhibits the highest HWN rate trend. In the range of 0.50–3.91 days/year, an upward trend in HWF is observed, with the greatest value belonging to Saphan Sung in southeastern Bangkok as shown in Figure 9d. In addition, the HWM-increasing trend ranges from 0.07 to $1.{71}^{\circ }\mathrm{C}$/year; the grid in Khlong Toei, located south of Bangkok, reveals an enormous value, as seen in Figure 9h and Figure 9j which reveal that the annual trend for the warmest HWA ranges from 0.13 to $1.{60}^{\circ }\mathrm{C}$/year, and the peak occurs in Suan Luang, Eastern Bangkok. The peri-urban region shows a clear trend of increasing nighttime temperatures, more so than rural and urban areas. Daytime trends also indicate higher spatial or temporal variation, with the peri-urban area experiencing a notable increase in heatwave exposure day and night.

3.4. Spatial homogeneity in correlation analysis

Table 7 presents the cumulative ’median’ of Pearson’s correlation coefficient (r) for the 10 meteorological stations and 10 LST grids that match heatwave indices. Our LST modeling outperforms the existing MOD11A1 and MYD11A1 products during a daytime heatwave, notably for HWN at $r=0.55,0.62$, and HWF at $r=0.66,0.71$, respectively. In other indices, their associations are slightly moderate: HWD $r=0.48$, HWM $r=0.32$ to $0.39$ , and HWA $r=0.35$. In contrast to daytime heatwaves, nighttime heatwaves demonstrate the link between all heatwave characteristics and a weaker result.

Figure 10 illustrates the distribution of the correlation coefficient (Pearson’s r) at $p<0.05$ corresponding to heatwave indices, considering pixel-by-pixel observations of $T_{\mathrm{air}}$ and LST. In rural, peri-urban, and urban areas, the link between $T_{\max }$ and MOD11A1 and $T_{\max }$ and MYD11A1 is positively strong, with r values larger than 0.60 and up to 0.90, respectively; the highest correlated MOD11A1 result is observed in P07 (located in Ayutthaya province as a rural context). In addition, the association between $T_{\max }$ and MYD11A1, particularly in rural regions in Lopburi province, P08 and P09, remains more robust than in other places ($r=0.80$) for HWN, HWF, and HWD. At $r=0.60$, however, peri-urban and urban stations are strongly related to HWM and HWA. In a peri-urban area, the connection between nighttime heatwaves corresponding to $T_{\min }$ and MOD11A1 is highest for HWN ($r=0.85$) and HWF ($r=0.75$) in Suvarnabhumi Airport (P05). In contrast, HWD shows the strongest link in the urban region (P03 Bang Na, Bangkok) with a correlation coefficient of 0.75, while P01 (Khlong Toei, Bangkok) displays the highest HWM ($r=0.65$) and HWA ($r=0.62$). The $T_{\min }$ against the MYD11A1 product indicates that P05, as the outskirt area, has the strongest link with all heatwave indices at $r=0.60$ to $0.85$, with the exception of P02 as the urban area, which has a smaller magnitude. Nighttime heatwaves are discovered to have a negative correlation in more regions than daytime heatwaves, particularly in rural areas.

4. Discussion

4.1. Accuracy and validation of land surface temperature predictive modeling

A crucial step in heatwave assessment is the optimization of the dataset during preprocessing, which includes a systematic treatment of missing values essential for accurate heatwave detection. The quality of LST data is often limited by cloud contamination, as noted by[52], particularly in heatwave assessment studies that require daily data inputs rather than instantaneous or averaged values. In this study, LST prediction analysis was conducted using the RF model, yielding compelling results similar to those found in [57]. This success is attributed to the employment of both temporal and spatial variables as significant predictors.

The selection of appropriate input variables is crucial before deep learning model training [104,105]. Our results align with the findings of [106], who observed that LST patterns are not constant and exhibit seasonal variations. Our selected predictors, more varied than those in previous studies, highlight the critical importance of temporal factors in LST analysis. This emphasis on temporal variables represents an advancement in understanding LST, suggesting that these factors may have a more pronounced impact on LST than previously recognized.

Our methodology aligns with the [107] study, which underscores the importance of a diverse range of environmental and land-use factors to LST prediction. Moreover, our work builds upon the findings of [52,108,109], which linked LST with factors such as land cover types, terrain, vegetation, moisture conditions, solar radiation, elevation, and Julian day. Interestingly, our research has distinctively identified the DOY as the most significant feature for predicting daytime LST (MOD11A1 = 0.28 and MYD11A1 = 0.36) as in Figure 3, owing to its critical role in influencing weather patterns through its capture of seasonal changes. Similarly, the daily $T_{\min }$ is crucial for nighttime LST prediction (MOD11A1 = 0.33 and MYD11A1 = 0.36), confirmed by the study of [110] stated that the average of nighttime LST was closest to $T_{\min }$, addressing diurnal temperature variations in both urban and non-urban settings.

The superior performance of the MOD11A1 night model, exhibiting an RMSE of $2.{09}^{\circ }\mathrm{C}$ and an $R^2$ of 0.64, aligns with theories suggesting that certain algorithms are more adept at processing and analyzing complex environmental data, particularly under stable atmospheric conditions at night. Our results are consistent with those of [63], who demonstrated the good spatiotemporal continuity and relative accuracy of reconstructed nighttime MODIS-LST products covering China. However, the temporal feature extraction, which involves incorporating time-related variables into the dataset, is significant. This process is reflected in the challenges faced by the MYD11A1 day model, which showed data overestimation issues with the highest RMSE ($5.{02}^{\circ }\mathrm{C}$) and the lowest $R^2$ (0.29). Contrast with the study of [111] found that MOD11A1 and MYD11A1 products slightly underestimated daytime LST with an overall absolute bias < $0.9^{\circ }\mathrm{C}$ and RMSE < $2.9^{\circ }\mathrm{C}$. The effective use of the MOD11A1 night model and the challenges associated with the MYD11A1 day model are consistent with theories in environmental data analysis. This demonstrates a nuanced understanding of how algorithms perform under various atmospheric conditions. The issue of data overestimation in the MYD11A1 day model highlights the practical difficulties in applying theoretical models, and this can be explained by the fact that during the daytime, solar radiation affects the thermal infrared signal, and the relationship between $T_{\mathrm{air}}$ and LST becomes more complicated [109]. Furthermore, adjusting hyper-parameters in the RF model can lead to overfitting, particularly when the rules are overly complex and tailored to the training data. Employing unseen test data and cross-validation methods can help mitigate this risk and ensure improved performance on new datasets [112,113].

4.2. Heatwave detection and its characteristic measurement

4.2.1. Ground-observed air temperature assessment

The results of our study reveal significant shifts in the patterns of both $T_{\max }$ and $T_{\min }$ heatwaves (daytime and nighttime) across urban, peri-urban, and rural regions from 1981 to 2019. The extended duration and earlier onset of heatwaves in recent decades, particularly post-2000, suggest the significant impact of climatic changes. Urban areas experience heatwave periods extending into late November, unlike earlier years, underscoring the intensifying UHI effect aligns with [114]. Conversely, rural areas show heatwaves starting earlier and ending sooner. This aspect may warrant further investigation to understand the underlying causes.

Our comprehensive analysis of $T_{\mathrm{air}}$ heatwave patterns in Thailand, detailed in Figure 4 and through broader study observations, reveals significant variations across socio-economic regions. During daytime heatwaves ($T_{\max }$), characteristically, peri-urban areas, such as Pathum Thani, and urban areas like Don Muang in Bangkok, demonstrate the highest increases in HWF, HWD, HWM, and HWA. Due to the fact that urban areas typically have heat retention properties such as concrete, asphalt, and steel, and inherent morphological characteristics, they are known to experience longer durations, higher frequencies, and greater magnitudes of heatwaves than other areas [115]. The analysis further highlights the highest number of each index of heatwaves in 1997, 2016, and 2019, indicating an escalation in heatwave events, similar to the study of [99,116,117,118,119,120,121,122,123], which identified severe heatwaves coinciding with the strongest El Niño years ever recorded to date, such as 1998, 2010, and 2016.

In contrast, the result from Figure 5, focusing on heatwave detection by $T_{\min }$, reveals distinct patterns and trends that differ from those observed for $T_{\max }$. Prominently, the most significant heatwave activity, particularly in terms of frequency and duration, is noted in urban areas (Bang Na, Klong Toei, and Don Muang in Bangkok), with peri-urban areas like Pathum Thani also showing high numbers of events. The years 2013 and 2019 stand out as particularly impactful, with the highest recorded values for various night heatwave indices. Both day and night, the highest number of heatwave metrics in 1997, 2013, 2016, and 2019 occurred due to severe heatwaves identified globally during the strongest El Niño event ever recorded to date [19,20,21,22]. These findings are in line with the expansive research in this field, as exemplified by the study by [15,124], which observes significant shifts of climatic changes affecting heatwave patterns, especially in urban areas, and highlights the critical importance of understanding regional heatwave trends for biophysical, human activity, and land use changes.

Table 4 reveals the mean annual heatwave indices during the day and night in urban, peri-urban, and rural areas. Distinctively, peri-urban areas (Samut Prakan and Pathum Thani) exhibit the highest percentages of emergence and the highest index values for both daytime and nighttime heatwave incidences, suggesting more pronounced and vulnerable heat impacts in these regions. Conversely, certain urban areas, such as Bang Na in Bangkok, demonstrate the lowest occurrence rates, but other locations experience notable frequency and intensity of these events. The mean annual numbers of HWD, HWM, and HWA suggest that nighttime heatwaves have higher numbers than daytime ones, particularly in urban settings, as the study of [125] indicates that the presence of UHI in Bangkok has increased in terms of intensity. The pattern observed in urban areas is mirrored in peri-urban and rural areas, albeit with varying degrees of severity.

4.2.2. Satellite-based land surface temperature assessment

The result uncovers distinct seasonal trends in heatwaves: urban and peri-urban areas typically experience daytime heatwaves from February or April to September or November, while rural areas face them from January to September. Nighttime heatwaves follow a similar pattern, with urban areas affected from February to July, peri-urban areas from mid-February to mid-October or November, and rural areas from late February to November. Underlined, these findings underscore the significant impact of geographical and urbanization factors on the timing and duration of heatwaves, anticipating similar trends in heatwave $T_{\mathrm{air}}$ patterns but showing higher values than those detected by the $T_{\mathrm{air}}$. These localized insights are crucial for developing targeted interventions and policies. In line with this, [126] states that heatwaves have become more frequent, persistent, and intense in recent decades, and even unprecedented events that would have been very unlikely without the increasing greenhouse gas (GHG) emissions are now occurring [127].

The evaluation of 268 patches through the LST dataset covering the area of interest pronouncedly advances our understanding of heatwave distribution. The findings reveal a marked difference in heatwave characteristics between daytime and nighttime; nighttime heatwaves tend to be more intense in terms of HWN, HWF, HWD, HWM, and HWA compared to daytime. This is similar to the results of [128], which declare that nighttime and compound heatwaves exhibit stronger increases in both frequency and fraction than daytime heatwaves. Peri-urban and rural areas consistently exhibit higher metrics. Additionally, the annual pattern of heatwave characteristics in specific regions (Figure 6 (daytime) and Figure 7 (nighttime)) highlights that daytime occurrences are more frequent and intense in peri-urban zones like Pathum Thani and eastern Bangkok, with annual episodes ranging from 2 to 9 and durations spanning 10 to 39 days. Conversely, urban centers such as downtown Bangkok are more susceptible to nighttime heatwaves, exhibiting a broader occurrence range of 3 to 12 events and longer durations of 13 to 62 days annually. To understand heatwave patterns, this finding is crucial for predicting and preparing for heatwave events, which are becoming more frequent and severe with climate change.

4.3. Temporal heatwave metric trend analysis

4.3.1. Ground-observed air temperature heatwave trend

In this study, spanning nearly four decades, we observed significant trends in daytime heatwaves, aligning with the recommendations of [15,129,130], who suggest a period of at least three to four decades to robustly assess changes in heatwaves. This duration is considerably longer than that proposed for average temperature trends. Remarkably, the peri-urban station of Pathum Thani exhibited the most pronounced and significant increases in daytime HWN, with an increase of 0.45 days/year, HWF, with an increase of 2.00 days/year, and HWD, with an increase of 0.27 days/year, respectively, as shown in Table 5. This trend indicates a heightened vulnerability of peri-urban areas to heatwaves, a concern that is amplified by factors such as rapid urbanization and reduced greenery. Additionally, local climatic factors contribute to exacerbating heatwave sensitivity in these regions. This finding is particularly concerning in light of the increasing trend of urban sprawl and the transformation of rural landscapes into peri-urban areas, which may further intensify heatwave conditions.

We observed distinct variations in nighttime heatwave patterns across different settings, challenging the notion that heatwaves are primarily an urban concern. Urban areas like Khlong Toei and Bang Na in Bangkok exhibited significant increases in heatwave metrics such as HWN, HWF, and HWD, but these increases were comparatively lesser than those observed in rural areas. As confirmed by [41,99], observed spatial and temporal consistencies in HWF and HWD can be attributed to the direct impact of the total HWN. Any variation in the total days contributing to heatwaves necessitates corresponding changes in both their HWD and annual HWF. Specifically, in the rural area of Lop Buri, located near Saraburi, the increase in HWN was notable at 0.39 days/year, followed by 0.25 days/year in Bang Na, an urban area of Bangkok. Pak Chong, located near Saraburi, also showed a considerable rise in HWN at 0.24 days/year, along with the highest increase in HWD at 0.15 days/year. This pattern was mirrored in the frequency of heatwaves, with Lop Buri recording a HWF of 1.44 days/year and Pak Chong following closely at 1.00 days/year, indicating a significant shift in heatwave trends in rural areas. Furthermore, Bang Na uniquely exhibited an increasing trend in HWA at $0.{07}^{\circ }\mathrm{C}$/year, underscoring the localized impacts of urban development on heatwave characteristics. Our findings emphasize that rural areas are particularly vulnerable to heatwave occurrences, a fact also confirmed by the study of [30], that heatwaves pose significantly elevated health risks for rural communities, as the $T_{\mathrm{air}}$ differences between heatwave and non-heatwave conditions in rural regions can be substantially larger than in urban areas. This discrepancy is largely due to synoptic scale high-pressure systems, which cause increased temperatures across extensive spatial scales. Additionally, under typical conditions (i.e., non-heatwave scenarios), rural areas usually experience lower temperatures.

Of particular interest is the unique trend observed in the urban area of Bang Na (Bangkok), which exhibited an increasing trend in HWA. This could indicate a distinct microclimatic alteration within this urban environment, potentially linked to the urban heat island effect. The increasing HWA in Bang Na underscores the impact of urban architectural and infrastructural elements on local climate conditions, raising concerns about the health and well-being of urban populations, especially in densely populated areas where the effects of heightened HWA could be more pronounced.

4.3.2. Satellite-based land surface temperature trend

As detailed in Table 6, the analysis reveals both positive and negative trends in heatwave metrics at a 95% level of significance. This statistical evidence indicates climatic changes that are directly impacting heatwave patterns. Specifically, our findings show a significant increase in nighttime heatwave occurrences, ranging from 41–57%, compared to daytime heatwaves, which exhibit an increase of 20-30%. This differential trend, observed over nearly four decades, underscores the necessity for distinct approaches in both studying and managing the varying impacts of heatwaves during different times of the day. The rising nighttime trends, particularly in urban settings like Bang Na and Suan Luang district in Eastern Bangkok, also imply a potential amplification of the urban heat island effect, particularly HWA shows increases of up to $1.{60}^{\circ }\mathrm{C}$/year. This increased trend may have been influenced by the sharp elevation of new residential development projects, warehouses, and buildings, as well as proximity to industrial estates and Suvarnabhumi International Airport, as indicated in [131]’s study, which found that a 10% increase in urban built-up density led to a 0.08% to 0.95% increase in heat wave number.

Regarding city-level vulnerabilities, the study specifically underscores peri-urban areas, such as Pathum Thani. These regions function as suburban areas supporting urban development and have experienced dramatic growth in urban services, multiple transportation options, and residential and industrial areas. Consequently, this expansion has resulted in the creation of a heat sink effect, primarily due to the extensive impermeable surfaces and substantial loss in vegetation cover [132], leading to noticeable increases in daytime heatwave metrics. It pronounced increases in daytime HWN, with increases of up to 0.83 events/year, HWF by up to 3.08 days/year, and HWD by up to 0.57 days/year. Urban areas like Khlong Toei in Bangkok demonstrate significant, albeit lesser, increases in similar heatwave metrics. These variations in heatwave impact are indicative of localized susceptibility, likely driven by factors such as rapid urbanization, evolving land use patterns, increasing impervious surfaces that retain high heat, and developed public infrastructure, including sky and underground trains. The peri-urban areas, characterized by a blend of urban and rural elements, are emerging as particularly vulnerable to heatwave conditions, necessitating focused adaptation and mitigation strategies. According to the findings presented in [128] study, daytime heatwaves are linked to elevated levels of solar radiation during dry conditions and a decrease in cloud cover and humidity under clear skies. Conversely, nighttime heatwaves tend to occur in more humid conditions, with higher cloud cover, humidity, and longwave radiation during the night.

Heatwave events, days, and duration were found to be significantly increasing (p < 0.05) for many populated regions [133]. The increasing trends in heatwave intensity and frequency observed in our study, particularly in peri-urban and urban areas, raise significant public health concerns. Therefore, our findings not only enhance our understanding of heatwave dynamics but also serve as a crucial guide for policymakers and urban planners. This guidance is essential for devising strategies to combat the escalating challenge of heatwaves, especially in rapidly urbanizing regions and areas transitioning from city centers to urban fringes. Particularly vulnerable populations, including the elderly, children, and individuals with pre-existing health conditions, as well as those living in poor-quality housing—lacking amenities like air conditioning, old house, and with thin panels—are especially susceptible to the adverse health effects of urban heat islands [134], such as those in overcrowded households in the Klong Toei and Bang Na communities in Bangkok.

4.4. Spatial homogeneity in correlation analysis between detected heatwave indices from ground-observed air temperature and satellite-based land surface temperature

Our investigation, employing Pearson’s correlation analysis, reveals a significant relationship between observed heatwaves based on $T_{\mathrm{air}}$ and MODIS-derived LST across various areas, particularly in its correlation with heatwave indices, as detailed in Table 7. This dual approach is comprehensive, as it covers both the atmospheric temperature (felt by residents) and the surface temperature (which influences the local microclimate). During daytime heatwave conditions, the findings are noteworthy, especially in the context of two specific heatwave characteristics: HWN and HWF, with $r=0.55-0.71$, indicating a strong positive correlation. Other aspects, such as HWD, HWM, and HWA, demonstrate a moderate association. Consistent with the findings of [55], the observed match percentages are relatively high, especially when considering the differences in terms of HWM. Our results contrast daytime and nighttime heatwave conditions, noting weaker correlations across all heatwave characteristics during the nighttime, with $r=0.02-0.45$.

In particular, the degree that determines the spatiotemporal pairwise correlation is initially measured, as shown in Figure 11. According to the findings of our research, grids located in rural areas have the potential to form the strongest associations compared to the other grids in terms of the HWN, HWF, and HWD, as measured by $r=0.93$, $r=0.94$, and $r=0.80$, respectively. On the other hand, the correlation coefficients for HWM and HWA are found to be at their maximum in a peri-urban area (Pathum Tani), with values of $r=0.65$ and $r=0.85$, respectively. Overall, this study represents a significant step forward in our ability to model and understand heatwaves through LST data, especially during daytime heatwaves, and it is crucial as it highlights the challenges in modelling nighttime heatwaves. This agrees with the findings of [50], which showed that the proposed LST index effectively identified heatwaves in the Mediterranean region during the daytime; however, this correlation was slightly weaker during the nighttime.

The use of MODIS-retrieved LST datasets for mapping heatwaves, particularly in peri-urban and rural areas with limited meteorological data, marks a significant advancement over traditional methods, as it potentially provides the repetitive imaging that heatwave episodes require due to their frequent sampling times [50]. The strong correlations in rural and peri-urban areas, coupled with the varied correlations in urban settings, underscore the intricate interplay between environmental factors and heatwave dynamics across diverse geographical landscapes. Our approach not only allows for comprehensive mapping of areas prone to heatwaves but also emphasizes the importance of extensive geographic analysis and demographic information in understanding heatwave patterns.

4.5. Summary of comparing research findings with related existing works

Here, we meticulously present a comparative analysis, elucidating the unique contributions of our current study in the realm of spatiotemporal heatwave quantification. In this comprehensive research endeavor, we delve into the intricacies of heatwave dynamics, uncovering critical insights that distinguish our work from the existing literature.

Studies by [57,109] offer validation for our utilization of MODIS-LST data and RF models. Notably, our study places added emphasis on factors both temporal (such as Day of Year) and spatial (such as elevation), showcasing their paramount importance in predicting $T_{\mathrm{air}}$ and LST with heightened accuracy. This novel insight contributes to the existing understanding and expands the scope of predictive variables, a dimension not explicitly explored in previous works.

Our results resonate with the findings of [110], reinforcing a robust correlation between surface air temperature (SAT) and LST, particularly evident in urban settings like Bangkok. This correlation implies heightened heatwave activities linked to urbanization. Importantly, our work extends beyond this correlation by encompassing a comprehensive analysis spanning over 39 years, fulfilling the crucial criterion highlighted by [15,129,130] for effective assessment of heatwave changes. This extended temporal analysis provides a more nuanced understanding of long-term heatwave dynamics.

Aligning with [131,132], our research underscores the escalating vulnerability of urban and peri-urban areas to heatwaves due to urban density and land-use changes. Noteworthy is the corroboration of our findings by [30], highlighting substantial health risks in rural areas during heatwaves, resulting in larger temperature differences compared to urban areas. This spatial variation emphasizes the need for tailored approaches to heatwave management.

Furthermore, the effectiveness of LST data in identifying daytime heatwaves, as underscored by [50,55], aligns with our observations, reinforcing the strong correlations between $T_{\mathrm{air}}$ and MODIS-LST data. This alignment substantiates the robustness of our analytical approach and highlights the consistency of these relationships across different studies.

The integration of our findings with the insights from [99,117,118,119,120,123,135,136] on the crucial influence of the El Niño in Southeast Asia further strengthens our study. Our identification of peak heatwave metrics in years coinciding with extreme warm events on an interannual timescale aligns seamlessly with their observations, offering a comprehensive understanding of the complex interplay between ENSO and heatwave dynamics.

In synthesizing our research with existing studies, we underscore the multifaceted nature of heatwave dynamics, advocating for nuanced approaches in both urban and non-urban heatwave management. The progressive nature of heatwave patterns highlighted in our study further emphasizes the imperative for sustained research efforts in this critical domain.

4.6. Limitations and Future Research

Given the constraints posed by the limited data availability at similar temporal scales between $T_{\mathrm{air}}$ and MODIS-LST data and incomplete observation records, this paper is considered the most accessible and available dataset. However, it is important to note that the focus of this study primarily lies in understanding and analyzing heatwave patterns, and as such, comprehensive solutions were not directly addressed. Several factors contribute to this limitation:

• Scope of the Study: The primary objective of this research was to analyze and characterize heatwave patterns in different socio-economic regions of Thailand. Delving into comprehensive solutions would require an expanded scope, considering various aspects such as public health, urban planning, and community resilience, which fall beyond the intended focus of this study.
• Interdisciplinary Nature of Solutions: Addressing the impacts of heatwaves involves an interdisciplinary approach that spans fields such as climatology, public health, urban planning, and community engagement. While our study contributes valuable insights into the patterns and characteristics of heatwaves, developing holistic solutions requires collaborative efforts from experts in these diverse fields.
• Data Limitations: The availability of comprehensive and reliable data is crucial for proposing effective solutions. Incomplete observation records and the limited temporal synchronization between $T_{\mathrm{air}}$ and MODIS-LST data pose challenges in formulating robust strategies.

5. Concluding Remarks and Possible Future Works

In this research, we addressed a significant gap in understanding the challenges posed by heatwaves in Thailand, a country that required up-to-date comparative analyses between urban and non-urban areas. Our pioneering approach combined geospatial analysis with remote sensing and ground data to track heatwave patterns in three distinct socio-economic regions: urban (Bangkok), peri-urban (Pathum Thani), and rural (Saraburi). This study contributed the most comprehensive collection of $T_{\mathrm{air}}$ and satellite-based heatwave data to date, providing an updated spatiotemporal analysis of extreme heat events in Thailand. Additionally, by integrating satellite-based LST data with ground observations, our methodology offered a more accurate and reliable approach for employing LST as a proxy in heatwave assessment. This included conducting an in-depth analysis of heatwave characteristics, focusing on their frequency, duration, intensity, and how these vary seasonally across different areas.

The findings about the performance of various machine learning models in land surface temperature prediction aligned with existing theories in environmental data analysis and machine learning. The effective use of the MODIS-MOD11A1 night model and the challenges faced with the MODIS-MYD11A1 day model resonated with theoretical expectations about the varying complexities of environmental data under different conditions. It highlighted the importance of context-specific model selection and the need for continuous refinement in predictive modeling, underscoring the dynamic nature of environmental data analysis. Therefore, future studies should consider ensemble methods, reflecting the theoretical stance that combining multiple models can enhance prediction accuracy and reliability—a concept widely supported in predictive modeling research.

Heatwave patterns in Thailand revealed significant variations across urban, peri-urban, and rural areas. Maximum air temperature heatwaves predominantly occurred from April to November in urban regions, while rural areas experienced them earlier in the year, from January to April. The highest incidence of heatwaves was recorded in the peri-urban area of Pathum Thani. Urban areas, particularly Bangkok, consistently showed the highest frequency and duration of heatwaves, with these trends intensifying over time. Nighttime heatwaves also followed a similar pattern, with urban regions experiencing more frequent and intense events. Overall, the study indicated an increasing trend in both the intensity and frequency of heatwaves across all regions, with urban and peri-urban areas being the most affected. Urbanization’s impact on heatwave patterns highlights its role in guiding urban planning and public health strategies. Research is needed for effective urbanization impact mitigation, addressing unique challenges in different urban environments. A proactive approach to policy and community engagement is imperative to address increasing heatwave trends in Thailand. Collaborative efforts and a comprehensive understanding are essential for effective policy implementation.

In summary, this research demonstrates the significant benefits of integrating various methods to deepen our understanding of climate patterns. Notably, our findings represent a major advancement as the first longitudinal study to evaluate heatwaves in Thailand. Such a comprehensive analysis has provided a clearer and more detailed understanding of heatwaves in different settings across Thailand. Moreover, our approach is particularly useful in regions with sparse or irregularly distributed meteorological stations, and in areas with distinct seasonal changes. This research is crucial for identifying and quantifying heat-related risks, contributing to informed decision-making for the public welfare. Consequently, it facilitates the development of more robust strategies against climate change, a vital step in mitigating the health impacts of heat and preparing for the escalating intensity of future extreme climate phenomena.

Our study lays the groundwork for understanding heatwave patterns, emphasizing the need for future research. Ensemble models show promise for improved predictive capabilities, offering nuanced insights. Additionally, exploring emerging technologies and socio-economic factors provides a holistic view of climate patterns and human activities. For future endeavors, the continuous refinement of predictive models is crucial for accurate LST predictions, enabling a detailed analysis of heatwave characteristics. Exploring temporal-focused predictors enhances the accuracy of LST models, benefiting environmental monitoring and climate change research. Further investigation into recalibration and alternative models is essential for adapting to evolving environmental conditions. Recognizing spatial variations in heatwaves, especially in peri-urban areas, calls for comprehensive studies that prioritize community awareness and preparedness.