The study areas consist of three megacities in Asia: Dhaka, Kolkata, and Bangkok (Figure 1). The three megacities are located in tropical or subtropical regions, and the seasonality of temperature and precipitation in these cities are similar. They are all categorized as tropical savannah (Aw) climate zone in the Köppen-Geiger classification scheme (Beck et al., 2018). Dhaka is the largest city in Bangladesh and one of the major cities in South Asia with a population of over 15 million. Kolkata is the capital of West Bengal and the third largest city in India after Mumbai and New Delhi. In addition, it is the largest trading and commercial center in eastern India which is located 150 km north of the Bay of Bengal above the plains of the Ganges Delta (Nath et al., 2015). As the capital and largest city of Thailand, Bangkok is the second largest city in Southeast Asia. It is the center of politics, economy, trade, transportation, culture, science and technology, education, religion and all aspects of Thailand. Since the study focused on SUHI analysis, urban functional areas of each city were delineated as main study areas.

Landsat data of the study area covering five time periods: 2000, 2005, 2010, 2015, and 2020 were acquired in this study. It was obtained from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/). One Landsat image with good data quality in each time period was selected for land cover classification and two Landsat images were both used for data fusion. The details about Landsat images for land cover classification are listed in Table 1 and the details for data used for fusion are listed in Table 2. Radiometric calibration and atmospheric correction were performed for the preprocessing of the Landsat images. The red, near-infrared, and thermal infrared bands were extracted for data fusion.

MODIS data is provided by the United States National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/). The MODIS data we used includes surface reflectance data from MOD09GA and MOD09A1 products and daytime land surface temperature data from MOD11A1 and MOD11A2 8-day composite products. Due to the frequent cloud coverage in the study areas, cloud-free remote sensing images in winter months (December to February) were selected to predict the Landsat-like land surface temperature data (Table 2).

The original MODIS data is stored in HDF format. The bands selected from MOD09GA/MOD09A1 are B01 and B02, while Emis_31 and LST_Day_1 km bands were selected from MOD11A1/MOD11A2. Batch data processing codes were used to convert the original MODIS data to geotiff format using MODIS Reprojection Tool. The extracted data bands were converted to surface reflectance by applying the scale factor and offset using the following equation: L λ = S c a l e   F a c t o r   ∗   D N + O f f s e t (1) where L λ is the surface reflectance of band λ; Scale Factor is the reflection coefficient of band; DN is the integer value of band λ, and Offset is the Offset of band.

The blackbody radiance can be calculated according to Eq. 2: B λ ( T S ) = K 1 e K 2 / T s − 1 (2) where K 1 and K 2 are Landsat sensor parameters, and T s is the true reflectance of band LST_Day_1 km.

The radiant brightness in thermal infrared band received by MODIS satellite sensor is calculated using the radiative transfer equation: L λ = ε λ   ∗   B λ ( T S )   ∗   τ λ + L a t m _ λ ↑ + ( 1 − ε λ )   ∗   L a t m _ λ ↓   ∗   τ λ (3) where L λ is the brightness of thermal infrared radiation with wavelength λ; B λ ( T S ) is the blackbody radiation brightness when the surface temperature is T S ; τ λ is the atmospheric transmittance along the path direction from the target to the sensor; ε λ is the surface emissivity with wavelength λ; L a t m _ λ ↑ and L a t m _ λ ↓ are atmospheric upward radiation and downward radiation. Atmospheric transmittance, atmospheric upward radiation, and atmospheric downward radiation can be calculated using the atmospheric correction parameter calculator provided by NASA. The processed MODIS data were resampled to 30 m with the same resolution as Landsat data. The preprocessing process is automatically performed using IDL scripts.

Land cover types including impervious surface, bare land, green space and water were classified using Support Vector Machine (SVM) classifier in the study areas for 2000, 2005, 2010, 2015, and 2020. Training and validation samples were randomly collected on Google Earth. The definition criteria issued by the United Nations Human Settlements Programme was used to determine the urban functional boundary (Lu et al., 2022). The urban functional boundary could reflect the urbanization situation in the study area more accurately than the administrative boundary. According to the land cover classification, the pixels with above 50% built-up density was defined as the urban pixel, 25%–50% was defined as the suburban pixel, and less than 25% was defined as the rural pixel. The area within 100 m around the urban or suburban pixels was defined as the urban open space, and the urban, suburban and urban open space pixels were merged. Then, the area of each polygon was expanded by 25% to establish buffer zones. The buffer zone with the largest area is the urban functional area, and its boundary is the functional boundary of the city. The urban functional boundaries of the three study areas were created using the land cover classification results in 2020 (Figure 1).

We used ESTARFM to fuse Landsat and MODIS data in this study (Zhu et al., 2010). Landsat and MODIS images at the time of T m and T n were input. The model searched for pixels similar to the center pixel according the window size. The conversion coefficients were determined according to the weights of similar pixels. The predicted images were finally obtained by combing two predicted results.

Among the commonly used land surface retrieval algorithms, the algorithm based on radiative transfer equation has a solid physical basis and high inversion accuracy (Lu et al., 2020; Duan et al., 2021). Radiative transfer equation-based method and the fused Landsat-like data were used to retrieve the LST in our study. The red band and near-infrared band of the fused image were used to calculate the normalized difference vegetation index (NDVI) (Lu et al., 2015; Yuan et al., 2018): N D V I = ( N I R − R E D ) ( N I R + R E D ) (4) where NIR is the near-infrared band and RED is the red band of the image.

Then, the vegetation coverage P v was calculated using NDVI: P v = N D V I − N D V I S o i l N D V I V e g − N D V I S o i l (5) where N D V I S o i l is the NDVI value of the bare soil or the area without vegetation cover, and N D V I V e g is the NDVI value of the pixel completely covered by vegetation. We take the empirical value: N D V I V e g = 0.70 and N D V I S o i l = 0.05. When the NDVI of a pixel is greater than 0.70, the P v value is 1. When NDVI is less than 0.05, P v is set as 0.

The land surface emissivity was then calculated as: ε = m   ∗   P v + n (6) by using typical constant values for the emissivity of vegetation and soil, m and n were set as 0.004 and 0.986, respectively (Sobrino et al., 2008).

The blackbody radiance was calculated using the radiative transfer equation: B λ ( T S ) = L λ − L a t m λ ↑ − ( 1 − ε λ )   ∗   L a t m _ λ ↓   ∗   τ λ ε λ   ∗   τ λ (7)

The land surface temperature was calculated as: T s = K 2 ln ( K 1 / B λ ( T s ) + 1 ) (8)

The temperature, T s , was converted to degrees Celsius. The LST retrieval was performed automatically using IDL scripts.

The surface urban heat island intensity (SUHII) can be defined as the difference between the surface temperature of built-up area and green space (Estoque and Murayama, 2017; Bechtel et al., 2019). At each built-up area pixel, the SUHII was calculated using the following equation: SUHII = T u r b a n − T r u r a l ¯ (9) where T u r b a n is LST at pixel i and T r u r a l ¯ is the land surface average temperature of green space pixels. The average temperature of dense vegetation with NDVI greater than 0.55 in green patch within urban functional boundary was taken as the average temperature of green space.

The SUHII of the study area was further divided into five levels as: very low level ( S U H I I <μ-std), low level (μ-std≤ S U H I I <μ-0.5*std), moderate level (μ-0.5*std≤ S U H I I <μ+0.5*std), high level (μ+0.5*std≤ S U H I I <μ+std), and very high level ( S U H I I ≥μ+std). The average of the annual SUHII in the urban functional area (excluding water area) was calculated as: SUHII c j = ∑ i = 1 n SUHII c i j ∑ i (10) where S U H I I c j represents the average SUHII of city C in the jth year (j = 2000, 2005, 2010, 2015, and 2020), i is the pixel i, and S U H I I c i j represents the value of SUHII of city c at pixel i in the jth year.

UGS patches were manually selected using Google Earth for each study area. Due to the complex composition of green space, only green space patches with dense trees were selected in order to control the influencing variables and dense trees are recognised as having strong cooling effects. The selection of green space patches follows several rules: select green patches that are primarily covered by trees, green patches that have no internal water spots and are more than 300 m away from other green spaces or water bodies to avoid interaction between samples, and green patches smaller than 900 m² (30 m*30 m) were removed. Based on these rules, 73, 90, and 107 green space patches were delineated in Dhaka, Kolkata and Bangkok, respectively.

The intensity and range of urban cooling island (UCI) were used to quantify the cooling effect of green space patches in our study. The intensity of UCI was calculated as the temperature difference between the temperature of the green patch and the peak temperature at the first turning point when the temperature decreased (Yu et al., 2018b). The cooling range is the distance from the edge of the patch to the first turning point on the temperature curve (Yu et al., 2017). In order to study UGS cooling effect, we used Arcpy scripts to create 18 buffers for each green patch. Since the image resolution is 30 m, the buffer distance is set to 30 m, and then the average LST of each buffer was calculated. According to previous studies (Fan et al., 2019), 540 m (18*30 m) could cover far more of the cooling range that a green space patch can produce. Then, the UCI intensity and range of each green space patch were calculated. The patches which LST dropped at the beginning were treated as abnormal patches and removed in our analysis.

Yu et al. (2017) defined the threshold value of efficiency (TVoE) to find the optimal patch size of urban green space in different study areas. UCI intensity increases rapidly in the first half of the period with the increase of patch area, while the change of UCI intensity is no longer obvious after exceeding a certain value (Yu et al., 2018a; Yu et al., 2018b; Fan et al., 2019). The log function is used to fit the change in UCI intensity according to Eq. 11: y = a   ∗   ln ⁡ x + b (11) where a is the coefficient of the log function; b is the constant of the log function. The critical point where the slope of the log function is equal to one is defined as TVoE. The UCI intensity increases significantly with the increase in patch size before the TVoE, and increase insignificantly after exceeding the TVoE. The logarithmic regression was used to extract the TVoE or optimal green space patch size in each city.

Following previous studies (Sun and Chen, 2017; Fan et al., 2019; Yu et al., 2020; Shah et al., 2021; Tan et al., 2021), four dominant factors, namely normalized difference vegetation index (P_NDVI), land surface temperature (P_LST), shape index (P_SI) and area (P_Area) of UGS patches were selected to explore their influence on the cooling intensity and range of UGS. The shape index represents the complexity of patch shape. The larger the value of the shape index, the more complex the shape is. The four metrics, their definition and description are listed in Table 3.

Pearson correlation was used to examine the relationship between the UCI effect and relevant landscape metrics. In specific, The UCI effect includes cooling intensity and cooling range and landscape metrics include P_NDVI, P_LST, P_SI, and P_Area. Pearson correlation coefficients were calculated to represent the magnitude of the correlation between UCI effect and relevant landscape metrics. Landscape metrics can affect UCI effect independently as well as interactively. Due to the possible correlation between landscape metrics, partial correlation analysis was performed to exclude the interaction between the landscape metrics.

Landscape metrics were used as independent variables, and the intensity and range of UCI effect were used as dependent variables. For landscape metrics with significant correlation relationship with UCI, regression analysis was conducted to analyze their effects on cooling effects. Because the relationship between the two is not necessarily a simple linear relationship, multiple functions were used to find the optimal regression. The R ² represents the degree to which the dependent variable can be explained by the model. The coefficient reflects the expected change in the dependent variable for every unit change in the associated independent variable. The ordinary least squares (OLS) regression was performed using SPSS25.0 software.

Table 4 shows the validation results of land cover maps for these three Asian megacities in five time periods. The overall accuracies range from 85.75% to 91.78%, and the Kappa coefficients range from 0.82 to 0.88. For Dhaka and Kolkata, low overall accuracies were obtained in 2000. It might be caused by the relatively low image quality.

The three cities have experienced significant urban development from 2000 to 2020. The urban areas have continuously expanded and green space has shrunk (Figure 3). The land cover changes in the three cities experienced similar changing trajectories. From 2000 to 2020, the green space of the three cities decreased by 1138.69 ha, 13894.52 ha, and 14416.24 ha respectively. Bare land decreased by 15257.98 ha, 20315.92 ha, and 109164.43 ha, respectively. The area of water bodies changed slightly in Dhaka and Kolkata. The impervious surface of each city expanded rapidly and increased by 16062.34 ha, 34707.11 ha, and 112698.38 ha respectively. The built-up areas in the three cities are rapidly expanding outward. The green space and bare land areas were mainly converted into impervious surface areas.

According to previous studies (Weng et al., 2014; Liu et al., 2020b), the observed Landsat data can be used to evaluate the accuracy of prediction results. The LST data obtained from the Landsat image on 13 December 2003 was used to estimate the accuracy of the LST from the predicted image. At the significance level of 0.01, MAD and RMSE between the two datasets were 0.26°C and 2.02°C, respectively, and Pearson’s value was 0.99. These results indicated that there was a strong correlation between predicted and observed LST (Figure 4).

Figure 5 shows the spatial pattern changes of SUHII in the three cities from 2000 to 2020. The High Level and Very High Level areas in all three cities continue to extend outward. The high level areas are concentrated in the city center in 2000, while they spread to the whole functional area of the city in 2020. The SHUII in the three cities show different spatial patterns. Multiple high intensity areas were observed in Dhaka in 2020. Some marginal areas even have higher heat island intensity than the central areas, and most high-intensity areas are concentrated in the central and northwest of the city. The SUHI in Kolkata sprawls in all directions from the city center. In Bangkok, the SUHI gradually spreads to the east and west from 2000 to 2020. Compared with the land cover map, it is obvious that the spatial evolution of SUHI corresponds well with the spatial pattern of urban expansion.

The temporal variations of SUHII and LST differences in the three cities were analyzed from 2000 to 2020 (Figure 6). The intensity of SUHI showed a monotonic increasing trends, and experienced a fluctuation in 2010 only in Bangkok. The changing slope of SUHII indicated that all three cities has been experiencing intensified SUHI effect. The continuous rise of SUHII is related with the increase of the difference between urban temperature and background temperature. The increase of urban surface temperature is significantly faster than that of background temperature. From 2000 to 2020, the urban surface temperature in three cities increased by 4.03°C, 8.39°C, and 4.47°C respectively, and the background temperature increased by 3.05°C, 6.97°C, and 3.74°C respectively.

Figure 7 shows the logarithmic regression analysis results in each city. Based on the fitted logarithmic regression models, the TVoE values were estimated. The p-value of each city is less than 0.01, indicating that the relationship is significant. The TVoE is 0.37 ha (R ² = 0.08), 0.77 ha (R ² = 0.31) and 0.42 ha (R ² = 0.12) in Dhaka, Kolkata and Bangkok, respectively.

Table 5 shows the statistics of landscape metrics and UCI intensity in the study areas. The mean patch size of greens space is larger in Kolkata and Bangkok than Dhaka. The shape index of green patches is largest in Bangkok. The average NDVI and LST of UGS are higher in Kolkata than the other two cities. The maximum UCI intensity ranges from 4.83°C in Bangkok to 8.07°C in Kolkata. The mean UCI intensity is highest in Kolkata and lowest in Dhaka. The maximum UCI range varies from 300 m in Bangkok to 420 m in Kolkata. The mean UCI range varies from 67 m in Bangkok to 81 m in Dhaka.

The correlation between landscape metrics and cooling effects in the megacities was also analyzed (Table 6). Among the four metrics, P_NDVI has a significant positive correlation with UCI intensity in Kolkata and Bangkok, but only has a significant positive correlation with UCI range in Bangkok. P_LST has a strong influence on UCI intensity and UCI range in all the cities, and the absolute values of Pearson correlation coefficients are greater than 0.50 on UCI intensity. For P_SI index, although correlation is negative on UCI intensity, the relationship is no significant. P_Area has a strong influence on UCI intensity and UCI range in the study areas.

For the significant correlations, the effects of landscape metrics were further analyzed using OLS regression (Table 7). The increase in NDVI leads to increase of UCI intensity in Kolkata (coefficient = 8.37) and Bangkok (coefficients = 4.83) and increase of UCI range in Bangkok (coefficient = 201.99). An increase in average LST is significantly associated with the decrease of UCI intensity and range in all the megacities. The effect of LST on UCI intensity in strongest in Kolkata, and its effect on UCI range is most significant in Dhaka. In addition, the increase of patch area is significantly associated with the increase of UCI intensity and range in all the megacities.

The study areas have a tropical monsoon climate, which is cloudy and rainy almost all year round. These lead to the difficulty in obtaining continuous gap-free remote sensing images. The blending of MODIS and Landsat data produced LST data with high spatio-temporal resolution. Since MODIS data have a spatial resolution of 1 km, narrow rivers and small-sized land features cannot be identified from the LST images. Compared with Hassan et al. (2021) which used MODIS data to explore SUHI effect, our data captured more detailed information (Figure 5). The overall land cover change results are similar to previous studies (Estoque et al., 2017; Si et al., 2022). The comparison between the land cover and SUHI map indicates that the expansion of urban impervious surface areas affects the variations of heat island. Due to the different physical properties of visible band and thermal infrared band, the visible bands and thermal infrared bands were fused separately, which effectively improved the accuracy of data fusion. The data fusion methods and processing framework have been developed to improve the efficiency and accuracy of data fusion. Since data gaps often exist in MODIS data, Liu et al. (2020b) used a mean filter or linear regression to fill the missing pixels, which enhanced data availability and provided richer information for fusion. Since the moving window of ESTARFM model has a fixed size, Liu et al. (2020a) used surface heterogeneity information to develop a new spatiotemporal data fusion algorithm which could automatically adjust the size of the moving window. These methods can be tested and applied to improve the accuracy of fused LST data in future studies. Based on the multi-temporal land cover and SUHII maps, the variations of land cover and urban thermal environment from 2000 to 2020 was captured clearly in the study areas.

The maximum intensity of UCI was found to range from 4.83°C to 8.07°C, and the maximum range of UCI varies from 300 m to 420 m in the three megacities. Our result show that the largest park in Kolkata (39 ha) provided a maximum cooling island effect with an intensity of 8.07°C and range of 420 m. Surveys in different seasons, times of day, and measuring methods indicated that the maximum cooling of large-sized parks could reach 5°C–8°C (Aram et al., 2019; Yin et al., 2022). The average night-time cooling range of Kensington Gardens (111 ha) in London was found to vary between 20 and 440 m (Doick et al., 2014). The Heiwa Park (147 ha), Nagoya had a cooling distance of 200–300 m during night and 300–500 m during daytime (Aram e al., 2019). The mean cooling range and intensity of blue-green space are 150 m and 2.47°C in summer in Copenhagen (Yang et al., 2020). Different TVoE values to provide the optimal green space patch area were obtained for each megacity. Our study shows that the TVoE is 0.37 ha, 0.77 ha, and 0.42 ha in Dhaka, Kolkata and Bangkok, respectively. Geng et al. (2022) applied the TVoE method to 207 urban parks in 27 cities in East China with four different local background climates, including warm temperate sub-humid monsoon, northern subtropical sub-humid monsoon, northern subtropical humid monsoon, and middle subtropical humid monsoon climate, and found that the TVoE was significantly affected by the background climate. Fan et al. (2019) found that the TVoE was about 0.6–0.62 ha in Hong Kong, Jakarta, Mumbai, and Singapore, and 0.92–0.96 ha in Kaohsiung, Kuala Lumpur, and Taiwan. Tan et al. (2021) found that the TVoE of three-based green spaces was about 0.31 ha in Nanning, China. The TVoE was found to be 4.55 ha in Fuzhou (Yu, et al., 2017). According to Fan et al. (2019), the TVoE was influenced by the background temperature and NDVI of green patches in the seven hot humid Asian cities. With higher background temperature, the optimal size of green space will be larger. This agrees well with our study.

Our results indicate that UCI intensity is correlated negatively with P_LST and positively with P_NDVI and P_Area, while UCI range is correlated negatively with P_LST and positively with P_Area. Since dense vegetation can increase the NDVI value and UCI intensity, the planting of woody vegetation is suggested to be increased in urban green space. The results also show that patch shape index has insignificant relationship with UCI intensity and range. Thus, city planners should choose the most convenient and economic way to plan green patches in terms of green space shapes. The effective method to enhance UCI intensity and range include increasing the area and biomass of green space, and lowering its background temperature. However, different results have been reported in previous studies. Some studies found that the correlation between P_Area, P_NDVI, P_SI, and UCI was insignificant (Shashua-Bar and Hoffman, 2000; Derkzen et al., 2017). The cooling effect of UGS is a huge and complex system (Murakawa et al., 1991). It is not only related to these landscape metrics, but also to other factors such as local climate background (Oliveira et al., 2011), geographical location (Hathway and Sharples, 2012), detailed landscape components (Zhao et al., 2014), humidity and evaporation (Vaz Monteiro et al., 2016; Zhou et al., 2017). In addition, ordinary linear regression (OLS) was used in our study. The OLS assumes that a linear relationship exists between the dependent and the independent variables. It also assumes error terms are independent and normally distributed. Spatial regression models such as spatial lag model, spatial error model and geographically weighted regression can be applied in order to account for spatial dependencies in the relationship between characteristic of UGS and intensity and range of UCI (Bartesaghi-Koc., 2022; Baqa et al., 2022).

Urban ecology is a huge system, and any environmental factors are likely to affect urban cool effects (Yu et al., 2017). The proposed methods can be adopted in other tropical cities to obtain more generalizable findings about the effects of urban greening characteristics on cooling effects in tropical climates (Bartesaghi-Koc et al., 2018). The physical characteristics of trees, the patterns of planting, and arrangement can also influence the air temperature and cooling effects of green infrastructures (Hami et al., 2019). The effects of these factors can be investigated by using more precise and comprehensive approaches in future studies. Moreover, the cooling effect of city parks varies in different seasons (Yang et al., 2020). The investigation of cooling effect of urban green spaces in different seasons can provide more information for the design of UGS patches. Urban morphology/geometry parameters such as sky view factor, aspect ratio, building and tree height, vegetation structure or stratification, orientation, altitude may also influence the thermal environment in cities (Irger, 2014; Morakinyo et al., 2020). They affect air circulation, heat dissipation, and thermal energy absorption in open spaces and urban canyons (Morakinyo et al., 2017; Ahmadi Venhari et al., 2019; Bartesaghi-Koc et al., 2022). A comprehensive investigation of urban and green space characteristics on urban thermal environment will be helpful to the planning of urban spaces to mitigate heat and improve thermal environment of cities in future studies.