Understanding spatiotemporal evolution of the surface urban heat island in the Bangkok metropolitan region from 2000 to 2020 using enhanced land surface temperature

Abstract

The urbanization process has significantly intensified surface urban heat island (SUHI) effects in the Bangkok Metropolitan Region (BMR). However, understanding the evolution of the urban thermal environment is challenging due to the difficulty in obtaining consistent remote sensing data of the cloud-prone landscape in the BMR. In this study, the data fusion algorithm was utilized to fill cloud-induced data gap and create high spatiotemporal-resolution data by blending Landsat and MODIS remote sensing images. The fused data was used to retrieve land surface temperature (LST) for winter months from 2000 to 2020. The spatiotemporal variations in SUHI were then captured using spatial cluster analysis. Finally, gradient analysis and geographically weighted regression (GWR) were employed to analyse the effects of land cover composition on LST. The SUHI intensity in winter increased from 4.40 °C in 2000 to 5.76 °C in 2020. The areal percentage of SUHI hot spots increased from 24.86% to 29.13%. The gradient analysis results indicated that vegetation with a density higher than 0.3 had a significant effect on LST compared to low-density areas. The woody lands were more effective in lowering LST than cultivated lands. These results provide useful information for developing heat mitigation strategies in the metropolitan regions.

Keywords:

• Spatiotemporal data fusion
• land surface temperature
• surface urban heat island
• sustainable development goals
• megacity

1. Introduction

As a complex socio-economic process, urbanization shifts the population distribution from dispersed rural areas to dense urban areas. The urbanization process converts natural land to impervious surfaces, which changes the thermal properties and energy budget of land surface. These changes further lead to changes in the thermal environment of urban areas and the urban heat island(UHI) effect (Zhao et al. 2014; Manoli et al. 2019). The UHI effect can exacerbate the intensity and frequency of extreme heat events (Marcotullio et al. 2021; Mentaschi et al. 2022), cause higher human health risks (Estoque et al. 2020), enhance environmental pollution (Henao et al. 2020), and increase energy consumption (Lu et al. 2019a). The investigation of UHI is of great importance for developing thermal mitigation strategies and enhancing thermal resilience of urban settlements, thereby contributing to the achievement of United Nation’s sustainable development goals in urban areas (United Nations 2015).

The UHI effect can be represented using the boundary layer heat island, canopy layer heat island, and surface urban heat island (SUHI) (Stewart and Oke 2012; Venter et al. 2021). The former two are also known as atmospheric heat islands, which are measured based on the temperature acquired from weather stations. SUHI is described using land surface temperature (LST) differences between urban and neighbouring areas based on remote sensing data (Zhou et al. 2019). With the development of thermal infrared remote sensing technologies, LST obtained from satellite sensors has been widely used in SUHI studies (Weng et al. 2018). In the 1990s, Landsat TM images were applied to reveal the LST differences among different land cover types (Xian et al. 2022). In the 2000s, the MODIS Terra sensor was used to retrieve LST at a ∼1 km resolution, but imagery with such a low spatial resolution may fail to capture urban structural details (Sobrino et al. 2012; Cheval et al. 2022). Remote sensing data with a spatial resolution lower than 50 m tend to insufficiently distinguish the thermal differences within cities (Stewart and Oke 2012; Ezimand et al. 2021). To obtain remote sensing data with both high spatial and temporal resolution, spatiotemporal data fusion models have been developed (Weng et al. 2014; Ghamisi et al. 2019). These models blend the frequent temporal information and high-resolution spatial information from different sensors. The spatial and temporal adaptive reflectance fusion model (STARFM) proposed by Feng et al. (2006) has been a widely adopted approach due to its high efficiency. The spatiotemporal adaptive data fusion algorithm for temperature mapping (SADFAT) was developed by introducing a conversion coefficient (Weng et al. 2014). Compared with these two methods, the enhanced spatiotemporal adaptive reflectance fusion model (ESTARFM) can further improve the fusion accuracy for complex and heterogeneous surfaces (Zhu et al. 2010). The ESTARFM model has been used to fuse Landsat and MODIS satellite imagery and described the spatial details of high-resolution LST data clearly (Liu et al. 2020). Using Landsat-like images produced by blending satellite data from different sensors, fine-scale spatiotemporal variations of surface urban heat islands and their relationship with geophysical parameters were revealed in Rasht city, Iran (Ezimand et al. 2021).

Since land cover categories such as impervious surfaces and vegetation can affect land surface temperature, the spatial distribution of LST shows high variability in urbanized regions (Lu et al. 2019b; Kedia et al. 2021). The urban-rural gradient approach has been used to explore urban growth patterns and its linkage with LST in Jakarta, Bangkok, Manila (Estoque et al. 2017), and the Greater Bay Area, China (Xie et al. 2021). The influence of land cover features on LST and UHI intensity has been explored in numerous studies (Lu et al. 2020; Dilawar et al. 2021; Singh et al. 2022). LST was found to be closely related to the density of impervious surface (positive) and vegetation (negative) areas (Tariq et al. 2020; Mondal et al. 2021). Despite the importance of understanding urban climate in the tropical zones, the knowledge on the impact of urban characteristics on urban heat island lags behind that of temperate climate (Giridharan and Emmanuel 2018). Although the association between land cover and LST has been reported in multiple tropical cities, such a relationship varies with different city sizes and different climatic characteristics (Chakraborty and Lee 2019; Cao et al. 2022). Comparing 150 major cities in India, Mohammad and Goswami (2021) reported that different variables such as vegetation, evapotranspiration and thermal inertia were responsible for the daytime SUHI intensity variations during summer/winter season in cities in different climatic zones. Climatic factors such as rainfall and solar radiation affect the shading and evapotranspiration intensity of urban vegetation, and their cooling effects on thermal environment are different in different regions (Cheung et al. 2021).

As the capital of Thailand, Bangkok has experienced severe SUHI effects due to extensive urban development in recent decades (Xu et al. 2019; Pimonsree et al. 2022). Air temperature records from meteorological stations revealed that the maximum SUHI intensity could reach up to 6-7 °C during the dry season in Bangkok (Arifwidodo and Tanaka 2015). The enhanced urban temperature, urban heat island, and global warming exacerbated the heat stress effects and brought negative impacts on the health and well-being of urban residents in Bangkok in recent years (Arifwidodo and Chandrasiri 2020a). The majority of UHI studies for the Bangkok Metropolitan Region (BMR) have used in situ meteorological data, since frequent cloud coverage makes it difficult to obtain consistent and continuous remote sensing data (Pakarnseree et al. 2018). Due to the uneven and sparsely distribution of weather stations, the meteorological data cannot reflect the spatiotemporal variation of the UHI phenomenon in this region accurately (Pakarnseree et al. 2018; Kim and Brown 2021). MODIS data has been used to explore SUHI trends in six major cities in Thailand (Keeratikasikorn and Bonafoni 2018a). However, temporal and spatial changes in SUHI were not clearly captured due to the low spatial resolution(1 km) of the MODIS satellite sensors. Therefore, a comprehensive investigation of SUHI using remotely sensed land surface temperature data with high spatial and temporal resolution is crucial to fill the knowledge gap for the development of heat mitigation and adaptation strategies in the BMR.

This study aimed to reveal the detailed spatiotemporal evolution of SUHI using an enhanced LST dataset in the BMR in last two decades. The second objective of this study is to examine the relationship between LST and land cover composition. The Landsat and MODIS remote sensing data were blended using the ESTARFM algorithm to create an enhanced LST dataset of the study area for the 2000–2020 period with a 5-year interval. Maps of SUHI were then obtained, and spatial cluster analysis was used to examine the spatiotemporal variation in SUHI intensity (SUHII). The relationships between LST and the densities of the built-up and vegetation areas were explored using gradient analysis and geographically weighted regression.

2. Study area and data

2.1. Study area

Bangkok is the largest metropolis on the Indo-China Peninsula and Southeast Asia’s second-largest city (Figure 1). It is located in the Chao Phraya River delta in Thailand’s middle plain with an average elevation of less than 2 m (Keeratikasikorn and Bonafoni 2018b). The built-up areas are concentrated in Bangkok’s central and eastern areas along the Chao Phraya River. There are six industrial and commercial districts in the BMR. Among them, the most prosperous is Bangkok Grand Kyoto, which is the political centre of Thailand. Bangkok has a tropical monsoon climate which is hot and rainy all year round. The annual average temperature is 30 °C, and the annual average precipitation is 1,500 mm. The year can be divided into three distinct seasons: summer (March to April), the rainy season (May to October), and winter (November to February).

Figure 1. Geographical location of the Bangkok Metropolitan Region. The land use map of BMR(c) is obtained from Department of Public Works and Town & Country Planning, Bangkok (https://dpt.go.th/en/).

2.2. Data source and pre-processing

In the study, the analysis was set in winter season, since more data were available and cloud-free remote sensing images cannot be obtained in the rainy season. Cloud-free Landsat 5/8, MOD11 A2, and MOD09 A1 images were used for LST retrieval in winter during the 2000–2020 period with a 5-year time interval. Landsat 8 Collection 2 level-2 surface temperature products were used to evaluate accuracy of LST prediction results. Satellite imagery from 2008 was used to supplement the imagery from 2010 and from 2014 supplemented 2015 due to the high cloud coverage in 2010 and 2015. Landsat 5 TM and Landsat 8 OLI/TIRS Collection 1 level-1 images were downloaded from the USGS Earth Explorer (http://earthexplorer.usgs.gov/). The MOD11 A2 and MOD09 A1 data products were obtained from NASA's Earth Observation Data and Information Systems (https://earthdata.nasa.gov/).

The remote sensing images were radiometrically calibrated and atmospheric corrected, and then was projected to UTM-WGS84 coordinate system and resampled to 30 m spatial resolution. For Landsat 5/8 satellite imagery, the surface reflectance data in the red and near-infrared bands corresponded to bands 1 and 2 in MOD09A1 imagery. Using the radiative transfer equation (RTE) based method (Chatterjee et al. 2017), spectral radiance was derived from the LST and emissivity data of MOD11A2 products. It corresponds to the spectral radiance in the Landsat 5/8 thermal infrared band. Using the ESTARFM model, Landsat-like LST data were generated from November to January in each time period (Table 1).

Table 1. Remote sensing datasets used in this study.

Landsat 5/8		MODIS 11A2/09 A1		Prediction date (8-day daytime)
First pair	Second pair	First pair	Second pair	November	December	January
28/02/2000	05/12/2003	05/03/2000-12/03/2000	03/12/2003-10/12/2003	24/11/2000-01/12/2000	18/12/2000-25/12/2000	09/01/2001-16/01/2001
07/12/2004	29/12/2006	10/12/2004-17/12/2004	27/12/2006-03/01/2007	17/11/2005-24/11/2005	11/12/2005-18/12/2005	09/01/2006-16/01/2006
30/01/2007	25/04/2009	25/01/2007-01/01/2008	14/01/2009-21/01/2009	08/11/2008-15/11/2008	02/12/2008-09/12/2008	17/01/2009-24/01/2009
17/01/2014	05/02/2015	17/01/2014-24/01/2014	02/02/2015-09/02/2015	07/11/2014-14/11/2014	11/12/2014-18/12/2014	01/01/2015-08/01/2015
17/12/2019	21/02/2021	19/12/2019-26/12/2019	18/02/2021-25/02/2021	16/11/2020-23/11/2020	02/12/2020-09/12/2020	09/01/2021-16/01/2021

3. Methods

Figure 2 shows the data processing workflow used in this study. The enhanced LSTs were first retrieved through the fusion of Landsat and MODIS data from the winter months of 2000, 2005, 2008, 2015, and 2020. The study then analyzed the temporal and spatial changes in SUHII in the BMR from 2000 to 2020. Finally, the relationship between SUHI and land cover was explored using gradient analysis and regression methods.

Figure 2. Data processing workflow.

3.1. Land cover classification

A support vector machine (SVM) classifier was employed to map the study area into five categories, including built-up areas, cultivated lands, woodlands, water bodies, and bare lands (Deliry et al. 2021). A stratified random sampling strategy was adopted to select reference samples. More than 150 samples were selected for land cover classification and accuracy evaluation for each land cover type. Among the samples collected for each land cover, 60% were training samples, and 40% were used for validation. The quantitative metrics such as users’ accuracy, producers’ accuracy, overall accuracy, and the kappa coefficient were used to evaluate classification accuracy. Finally, the area of each land cover type was extracted from the classified imagery.

3.2. Land surface temperature retrieval

The ESTARFM algorithm is specifically tailored for data fusion of complex terrain regions and includes four main steps (Figure 2). First, two Landsat images at $t_w$ and $t_q$ were used to find pixels within a sliding window that are spectrally similar to the centre pixel. ESTARFM uses the threshold approach to find similar neighbourhood pixels based on Equation (1)(1) $$\left|F\left(x_i,y_i,t_p,B\right)-F\left(x_{\omega /2},y_{\omega /2},t_p,B\right)\right|\le 2\times \sigma \left(B\right)$$(1) (Zhu et al. 2016). (1) $$\left|F\left(x_i,y_i,t_p,B\right)-F\left(x_{\omega /2},y_{\omega /2},t_p,B\right)\right|\le 2\times \sigma \left(B\right)$$(1) where $x_i,y_i$ defines the i^th pixel; $x_{\omega /2},y_{\omega /2}$ is the centre pixel of the sliding window; $\omega$ is the sliding window size; $F\left(x_i,y_i,t_p,B\right)$ and $F\left(x_{\omega /2},y_{\omega /2},t_p,B\right)$ respectively represent the i^th pixel and centre pixel at $t_p$ high spatial resolution images on band B; $\sigma \left(B\right)$ represents the standard deviation of band B at $t_p$ time; and $m$ refers to the estimated feature category.

The conversion coefficient was obtained using a linear regression between Landsat and MODIS images at $t_w$ and $t_q.$ The weight of a spectrally similar pixel was computed using the location, reflectance, and thermal similarity between Landsat and MODIS data. Subsequently, two approximated Landsat-like images were created at the prediction time $t_p.$ Finally, the 30 m resolution image was estimated using a weighted combination of the two Landsat-like images.

This study used fused Landsat-like data and the radiative transfer equation (RTE) method to retrieve land surface temperatures. The RTE method is a widely used LST retrieval algorithm with high accuracy (Malakar et al. 2018). It is described using Equation (2)(2) $$L_{\lambda }=B_{\lambda }\left(T_s\right){\epsilon }_{\lambda }{\tau }_{0\lambda }+L_{0\lambda }^{\uparrow }+\left(1-{\epsilon }_{\lambda }\right)L_{0\lambda }^{\downarrow }{\tau }_{0\lambda }$$(2) (Sekertekin and Bonafoni 2020): (2) $$L_{\lambda }=B_{\lambda }\left(T_s\right){\epsilon }_{\lambda }{\tau }_{0\lambda }+L_{0\lambda }^{\uparrow }+\left(1-{\epsilon }_{\lambda }\right)L_{0\lambda }^{\downarrow }{\tau }_{0\lambda }$$(2) where $L_{\lambda }$ represents the radiation received by the sensor; $B_{\lambda }\left(T_s\right)$ is the blackbody radiance given by Plank’s law; $T_s$ is LST; ${\epsilon }_{\lambda }$ and ${\tau }_{0\lambda }$ represents land surface emissivity and atmospheric transmissivity, respectively; and $L_{0\lambda }^{\uparrow }$ and $L_{0\lambda }^{\downarrow }$ represent upwelling and downwelling atmospheric radiance, respectively. The Normalized differential vegetation index (NDVI) threshold method was used to calculate land surface emissivity. The land surface temperature can be calculated by Equation (3)(3) $$T_{\mathit{LST}\left(K\right)}=K_2/Ln\left[1+K_1/L_{\lambda }\right]$$(3) . (3) $$T_{\mathit{LST}\left(K\right)}=K_2/Ln\left[1+K_1/L_{\lambda }\right]$$(3) where $T_{\mathit{LST}\left(K\right)}$ represent LST in Kelvin. $K_1$ and $K_2$ are the calibration constants for the sensors.

By comparing the predicted Landsat-like LST with observed Landsat data on the same date, the accuracy of the predicted LSTs was evaluated. The root mean square error (RMSE), mean absolute difference (MAD), and correlation coefficient (CC) between estimated and observed LST were calculated to quantitatively evaluate the estimation error of the fused LST data.

3.3. Surface urban heat island intensity quantification

Following the local climate zone framework (Stewart and Oke 2012), the difference in LST between urban and dense trees was used to calculate the surface urban heat island intensity (SUHII). This LST difference was measured using Equation (4)(4) $$\mathit{SUHII}=T_i-\mathrm{ }{\overline{T}}_{gs}$$(4) (Schwaab et al. 2021). (4) $$\mathit{SUHII}=T_i-\mathrm{ }{\overline{T}}_{gs}$$(4) where $T_i$ represents the land surface temperature at pixel i and ${\overline{T}}_{gs}$ represents the average temperature of green space with dense trees. Based on land cover classification results, woodland pixels with NDVI > 0.5 were regarded as densely-vegetated green space. SUHII was further classified using combinations of the mean value and standard deviations of μ (mean), μ ± 0.5 std (standard deviation), and μ ± 1 std, yielding five grades (Table 2).

Table 2. Classification of SUHII.

SUHII grade	Range
Very low level	SUHII < μ - std
Low level	μ - std ≤ SUHII ≤ μ − 0.5 std
Moderate level	μ − 0.5 std ≤ SUHII ≤ μ + 0.5 std
High level	μ + 0.5 std ≤ SUHII ≤ μ+ std
Very high level	μ + std ≤ SUHII

3.4. Spatial cluster analysis

Multi-scale segmentation was used to create image objects to analyze the spatial pattern of LST in the study area. The segmentation of LST images was performed using eCognition Developer software (http://www.ecognition.com/). A scale parameter of 20 was used since it produced meaningful objects with 'natural’ geographic borders (Guo et al. 2015). The study area was divided into more than 7,500 objects using Landsat data. Getis-ord $G_i^{\mathrm{*}}$ is a spatial cluster analysis method. Its basic principle is to calculate each element in the data set with Equation (5)(5) $$G_i^{\mathrm{*}}=\frac{{\sum }_{j=1}^n{\omega }_{i,j}\mathrm{*}{x}_i}{{\sum }_{j=1}^n{x}_i}$$(5) (Lemus-Canovas et al. 2020) and obtain the Z-score and p-value. (5) $$G_i^{\mathrm{*}}=\frac{{\sum }_{j=1}^n{\omega }_{i,j}\mathrm{*}{x}_i}{{\sum }_{j=1}^n{x}_i}$$(5) where $G_i^{\mathrm{*}}$ is the spatial correlation index; ${\omega }_{i,j}$ is the spatial weights matrix; and $x_i$ is the observation $x.$ A feature with a high Z-score and modest p-value implies a substantial hot spot. A substantial cold spot is indicated by a low negative Z-score and a small p-value. This method can be used to evaluate the spatial characteristics of LST at fine spatial scales.

3.5. Relationship analysis of land surface temperature and land cover

Based on the SVM classification results, the built-up land and vegetation areas (woodland and cultivated land) were selected to explore the relationship between LST and land cover on the urban-rural gradient. Following the previous urban-rural gradient analysis studies (Xu et al. 2019; Kaminski et al. 2021) and considering the spatial resolution of Landsat thermal images, this study created a total of 45 concentric buffer zones with 1 km intervals around the city center of Bangkok (13°45'7"N, 100°31'48"E) (Bonafoni and Keeratikasikorn 2018). The areal proportions of the built-up and vegetation areas to the total area within each buffer were calculated as the built-up and vegetation density, respectively.

The geographically weighted regression (GWR) method, which considers the spatial heterogeneity of the correlation relationship, was used to explore the spatial association between LST and land cover (Fotheringham et al. 2017). Firstly, the study area was divided into 6,737 grid cells with a size of 1 km × 1 km. Built-up area, cultivated land and woodland accounted for a large proportion in the study area. The density of three land cover types in each grid cell were calculated and selected as explanatory variables of LST. The Equation (6)(6) $$y_i={\beta }_0\left(u_i,v_i\right)+{\sum }_k{\beta }_k\left(u_i,v_i\right)x_{k,i}+{\epsilon }_i$$(6) was used for GWR analysis (Fotheringham et al. 2017). (6) $$y_i={\beta }_0\left(u_i,v_i\right)+{\sum }_k{\beta }_k\left(u_i,v_i\right)x_{k,i}+{\epsilon }_i$$(6) where ${\mathrm{ }y}_i$ is the response variable, ${\beta }_0$ and ${\beta }_k$ are the fitting parameters, $x_{k,i}$ is the kth explanatory variable in region $i,$ $\left(u_i,v_i\right)$ is the geospatial location, and ${\mathrm{ }\epsilon }_i$ is the random error in region $i.$

4. Results

4.1. Land cover changes

Figure 3 shows the land cover classification results for the BMR in 2000, 2005, 2008, 2014, and 2020. Most of the built-up area was located in the central part of the metropolitan area, while cultivated lands and woodlands were mainly distributed in the surrounding areas. Water bodies were scattered in the southern coastal area. Most of the built-up area increase was converted from vegetated areas.

Figure 3. Land cover classification results in Bangkok Metropolitan Region from 2000 to 2020.

Table 3 shows the accuracy assessment results (producer’s accuracy (PA) and user’s accuracy (UA)) for the land cover classification maps. The overall accuracy of the land cover maps ranges from 90.45% to 93.90%. The kappa values range from 0.87 to 0.92. The overall classification accuracy was above 90% for each year. The land cover change analysis indicates that the built-up and bare land areas increased, while cultivated lands, woodlands, and water bodies decreased from 2000 to 2020. Among the different land cover types, the built-up areas exhibited the maximum areal increase, with 71.97% of the area converted from cultivated land.

Table 3. Accuracy assessment results for the land cover classification maps.

Land cover type	2000		2005		2008		2014		2020
	PA. (%)	UA. (%)	PA. (%)	UA. (%)	PA. (%)	UA. (%)	PA. (%)	UA. (%)	PA. (%)	UA. (%)
Built-up area	97.62	91.62	98.31	90.67	99.44	99.44	95.40	92.74	94.25	98.20
Bare land	79.75	94.03	100.00	88.52	99.01	98.96	61.11	84.62	99.90	93.02
Cultivated land	85.35	97.58	93.39	90.75	91.67	86.27	83.45	80.00	92.63	79.28
Woodland	97.12	87.83	68.52	85.38	75.00	94.44	75.36	92.86	76.19	82.76
Water body	96.19	80.37	98.01	98.34	99.01	95.22	98.69	94.97	96.30	99.52
Overall accuracy (%)	90.57		92.09		93.90		90.88		93.20
Kappa coefficient	0.87		0.89		0.92		0.87		0.91

4.2. Land surface temperature retrieval results

The fused Landsat-MODIS images were used for LST retrieval in the study area from November to January in each time period. Figure 4 shows the predicted Landsat-like LST in the study area. The fused LST more clearly reflected the spatial distribution of rivers, roads, and urban gatherings in comparison to MODIS. The texture features of the spatiotemporal fusion data were clearly observable.

Figure 4. Observed Landsat on 19/12/2020 (a), predicted Landsat-like LST on 09/01/2021-16/01/2021(b), and MODIS LST on 09/01/2021-16/01/2021(c).

Figure 5 shows the comparison between two Landsat 8 LST products and the concurrent fusion results. Through visual comparison, the spatial distribution of predicted LST and Landsat 8 LST products are highly consistent. The RMSE between the two datasets was 1.73 °C and 1.64 °C, MAD was 0.86 °C and 0.79 °C, and CC was 0.82 and 0.89, respectively. It indicates that there was a strong correlation between the predicted LST and the Landsat LST product in most regions. The predicted LST on both dates are slightly higher than that of Landsat 8 products.

Figure 5. Landsat 8 LST and predicted Landsat-like LST on 01/04/2015 (a, b) and 17/11/2020 (d, e), and scatter plots (c, f).

4.3. Spatiotemporal variation of the surface urban heat island

Figure 6 shows the spatial pattern of the average winter SUHII from 2000 to 2020. SUHII was divided into five grades based on the mean standard deviation. The SUHII generally increased across time as built-up areas expanded except in 2008. In 2008, Bangkok showed the lowest value of SUHII. The enhanced precipitation may result in the low overall SUHII and large uncertainty in LST retrieval results in 2008 (Limsakul and Singhruck 2016; Pimonsree et al. 2022). The very high-level areas gradually expanded and the most noticeable SUHII increases occurred in the eastern part of the study area. In 2000, very low-level SUHII areas were in Bangkok’s urban core. In 2020, the urban core exhibited moderate-level and high-level SUHII values, while large areas with high and very high-level SUHII values were observed in the eastern BMR.

Figure 6. Mean SUHII in the Bangkok Metropolitan Region in winter season from 2000 to 2020.

The areal proportion of SUHII in different grades during the study period is shown in Figure 7. Aside from 2008, the proportion of high-level and very high-level SUHII areas increased from 24.89% in 2000 to 46.78% in 2020, with an annual average growth rate of 1.10%. The proportion of low-level and very low-level SUHII areas decreased from 56.36% to 19.94%. The moderate-level SUHII areas accounted for the largest proportion of the study area with little temporal change. The proportion of very low-level SUHII is the highest in year 2008, while the areal percentage of high-level SUHII is the lowest.

Figure 7. Areal percentage of different SUHII grades in the Bangkok Metropolitan Region from 2000 to 2020.

Figure 8 shows the LST spatial cluster results for the study area. Hot spots were mostly concentrated in the urban core, and cold spots were scattered in the urban fringe and southern coastal areas. The hot spot areas were observed to spread from the urban centre to the surrounding areas over time. In addition, the hot spots in the city centre have decreased slightly. One of these areas is Bang Kachao park which area has increased during the study period. This indicates that the government has made some efforts to improve public parks and green spaces in the central area of the city and the SUHI intensity was reduced in these areas. The spatial variation of cold and hot spots captured the changes in the thermal environment of the study area over time effectively.

Figure 8. LST spatial clusters in the Bangkok Metropolitan Region from 2000 to 2020.

The proportion of hot and cold spot cluster areas showed an increasing trend over time (Table 4). Overall, the area of hot spot clusters was larger than cold spot clusters in the BMR. The areal percentage of LST hot spots increased from 24.86% to 29.13% during the study period. The average LST of cold spot clusters increased from 28.38 °C in 2000 to 31.56 °C in 2020. The average LST of hot spot clusters increased from 34.61 °C to 38.79 °C.

Table 4. Cold and hot spot clusters identified in the study area from 2000 to 2020.

Year	Cold Spot		Hot Spot
	Areal proportion (%)	LSTmean (°C)	Areal proportion (%)	LSTmean (°C)
2000	17.54	28.38	24.86	34.61
2005	18.33	29.71	23.29	35.39
2008	23.64	27.24	27.42	31.29
2014	23.39	30.39	27.58	34.78
2020	22.79	31.56	29.13	38.79

4.4. Relationship between land cover and land surface temperature

The relationship between LST and land cover types was analysed in five time periods (Figure 9). Generally, the LST in built-up areas was much higher than in other land cover types, and LST was lowest in water bodies. The LST of the four land cover types showed an increasing trend from 2000 to 2020 except for 2008. Among them, the LST in built-up areas showed the greatest change, increasing by about 6.20 °C from 2000 to 2020.

Figure 9. Box plots of land surface temperatures for four land-cover types: (a) built-up areas; (b) cultivated lands; (c) woodlands; and (d) water bodies.

Gradient analysis shows that the density of built-up areas decreased while the vegetation density increased from the city centre to the surrounding suburbs and rural areas (Figure 10a and b). High LST values with little variation were observed in the core region with high built-up density, which is within 15 km from city center. The built-up area density increased over time. Intensive built-up area sprawl in the transitional areas between urban and outskirt areas was observed from 2000 to 2020 (Figure 10a). The mitigation effect on LST was stronger in the area with vegetation density ranging from 0.30 to 0.80 than in low-density areas (Figure 10b). The relationship between built-up density, vegetation density, and LST was analysed within each buffer (Figure 10c and d). There were significant correlations between land cover densities and LST, with R² values above 0.87. The density of built-up areas had a significant positive relationship with LST, while vegetation density had a negative relationship with LST. The slope value of the vegetation density-LST relationship was greater than in the built-up density-LST relationship. This indicates that the effect of vegetation density on LST was greater than that of built-up land. The results are in line with those obtained in previous studies (Adulkongkaew et al. 2020).

Figure 10. The mean LST vs. built-up density with the distance from the city center (a), the mean LST vs. vegetation density with the distance from the city center (b), relationship between built-up density and LST from 2000 to 2020 (c), and relationship between vegetation density and LST from 2000 to 2020 (d).

In 2020, the GWR results spatially varied for different land cover types (Figure 11). In the core regions, the built-up density was close to 1, and the correlation between LST and urban built-up density was low. In the urban-rural fringe, the density of built-up areas greatly changed across space and featured a high correlation with LST. Cultivated land was scattered in patches in the core urban area and was contiguously distributed in the suburbs. Woodlands had a wider distribution of low regression coefficients compared to cultivated land. The land cover composition had an insignificant relationship with LST in regions with a high proportion of mixed land use.

Figure 11. GWR analysis results for 2020.

5. Discussion

5.1. Implications for urban heat island mitigation

Our analysis showed that 71.97% of the newly increased built-up areas in the BMR were converted from agriculture land between 2000 and 2020. The land in the urban core was dense and compact; however, there were many fragmented built-up areas in the outskirts and hinterland (Xu et al. 2019; Adulkongkaew et al. 2020). The built-up area became more aggregated across time. During the study period, urban expansion grew radially from the urban core to the surrounding areas along the Chao Phraya River. The gradient analysis results indicate that urban land density showed an inverse S-shaped function trend based on the distance from the urban center.

The temporal changes of LST showed an increasing trend in the four land cover categories. The average LST of built-up land increased most significantly from 31.03 °C in 2000 to 37.23 °C in 2020. The intensified SUHI have been revealed in a number of cities in developing countries (Mathew et al. 2018; Sultana and Satyanarayana 2020). The intensities of UHI were observed to be increasing during winter season and the hotspots with high LSTs were identified in high-density constructed region or dry barren lands in Indian cities such as Kolkata, Mumbai and Hyderabad (Sultana and Satyanarayana 2020). Bankgkok has experienced wetter conditions and a concomitant increase in the frequency and magnitude of heavy precipitation events from 1955 to 2014 (Limsakul and Singhruck 2016; Pimonsree et al. 2022). Located at the mouth of the Chao Phraya River, serious floods (such as the 2011 flood) usually take place in October and last for one or two months when the upstream flood and local rainstorm occur at the same time (Worawiwat et al. 2021). The low SUHII in Bangkok in 2008 might be caused by the heavy precipitation and flooding during the study period (Arifwidodo and Tanaka 2015).

The spatial clustering analysis results of Bangkok show that the high LST clusters are mainly distributed in the central urban construction land intensive areas, while the low LST clusters are mainly scattered in the surrounding areas of the city. The area of hot spots increased from 24.86% in 2000 to 29.13% in 2020, showing an obvious growth trend. Previous studies found that the built-up density greatly impacted SUHI in Bangkok, and the building coverage and floor area ratio were significantly linked with SUHI (Pakarnseree et al. 2018). Despite the continuous expansion of the high SUHI intensity area in the mostly urbanized central area, few SUHI hotspots were identified in the newly developed areas with proper planning and green space in Kolkata, India (Halder et al. 2021). With the improvement of road transportation facilities such as subway and light rail in Bangkok in recent years, residents tend to settle in Nonthaburi and Pathum Thani areas, which intensifies the LST in northern Bangkok. Public green spaces such as Bang Kachao Park should be built to alleviate heat stress in the newly developed settlements in the north (Arifwidodo and Chandragiri 2020b).

The effects of land cover composition on LST have been studied extensively (Mondal et al. 2021). Keeratikasikorn and Bonafoni (2018b) found that the average SUHI was lowest in rural and agricultural areas and highest in commercial areas in the BMR. This was related to the concentration of commercial land use in the city center. The association between LST and vegetation density indicates that dense and aggregated vegetation areas were very effective at reducing the temperatures in the BMR. An analysis of 44 major cities in India found that the SUHII was intense during daytime in cities in tropical climatic zones such as Chennai and Kolkata and the impact of vegetation on the determination of SUHII is obvious (Raj et al. 2020). A recent study in Delhi, India also reported a high degree of positive correlation between LST and Normalized difference built-up index(NDBI) (R² = 0.8369) and negative correlation between LST and NDVI (R² = −0.8872) (Singh et al. 2022). This is consistent with the results that the mitigation effect of vegetation on SUHI is stronger than the warming effect of buildings in this study. The GRW analysis further reveals that woody vegetation is more effective in lowing LST than agricultural lands. Therefore, urban green space and forests should be planned and provided to mitigate SUHI in BMR. According to Adulkongkaew et al. (2020), the LST in the BMR can be greatly reduced by expanding the tree canopy to around 20%, water bodies to around 30%, or green yards to about 40% in urban areas. As reported by Estoque et al. (2017), urban planning and management should be strengthened in the areas 2–12 km away from the city center. The mitigation strategies of UHI effect can be classified into three types: roof strategies (high reflectivity roofs, vegetation roofs), non-roof strategies (shading structures with energy generation/architectural devices, high reflectance paving, plantation, and water bodies), and covered parking strategies (Khare et al. 2021). The UHI effect can be considerably minimized by employing these measures.

5.2. Limitations and future perspectives

The spatiotemporal characteristics of SUHI provide valuable information for land use planning and management to mitigate SUHI in the study area. However, this study has several limitations. First, the study only investigated the daytime SUHII in winter due to the lack of remote sensing data. SUHI effects showed significant diurnal and seasonal variability (Mohammad and Goswami 2021). Multi-source remote sensing data can be applied to explore the diurnal and nocturnal LST changes in different seasons. Secondly, detailed land use types can be classified to examine the impact of land use changes on SUHI using remote sensing imagery with high spatial resolution (Keeratikasikorn and Bonafoni 2018b). Finally, SUHI patterns can be influenced by different elements, such as urban morphology, vegetation cooling, and background climate. Fine-scale datasets such as three dimensional morphology data of buildings and field measurements can facilitate the exploration of SUHI influencing factors and the implementation of practical heat mitigation measures (Pakarnseree et al. 2018; Arifwidodo and Chandrasiri 2021).

6. Conclusions

In this study, high spatiotemporal resolution remote sensing data were used to explore the variations in the thermal environment and their interactions with land cover features in the Bangkok Metropolitan Region. Monthly LST data were estimated using Landsat-like imagery from 2000 to 2020 with a five-year interval in the winter months. The rise in mean SUHII from 4.40 °C in 2000 to 5.76 °C in 2020 demonstrated an overall surface warming trend in the study area. The areas with high-level and very high-level SUHII values were mainly distributed in the urban core and eastern BMR. Low-level and very low-level SUHII areas were located on the urban fringe and coastal regions. The areal percentage of LST hot spots increased from 24.86% in 2000 to 29.13% in 2020. The gradient analysis results indicate that vegetation had a stronger effect on LST than the built-up area, while the effects of built-up land increased from 2000 to 2020. Woodlands had a more widely distributed cooling effect compared to agricultural lands. The increase in medium and high-density woody vegetation and green space would effectively mitigate the negative effects of SUHI in the BMR, especially in suburban and newly developed areas.

Although the present study provides insight for SUHI mitigation in the study area, the seasonal and nocturnal variations of urban thermal environment can be investigated in future studies. The application of more detailed data such as high-resolution remote sensing data, building morphology and field measurements can enable better understanding of the drivers of SUHI and the development of practical thermal mitigation strategies.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The datasets used in this study is available from the corresponding author (Linlin Lu) upon reasonable request.

Additional information

Funding

The authors are grateful to U.S. Geological Survey for providing Landsat dataset and NASA's Earth Science Data Systems (ESDS) program for providing MODIS products. The authors also would like to express their appreciation to the anonymous reviewers and editors for their helpful comments and suggestions. This research was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant number XDA19030101); the Director Fund of the International Research Center of Big Data for Sustainable Development Goals (grant number CBAS2022DF016); and the National Natural Science Foundation of China (grant number 42071321).