key: cord-0021309-tvaaynb3
authors: Abir, Farhan Asaf; Ahmmed, Sabbir; Sarker, Soykot Hossain; Fahim, Ashraf Uddin
title: Thermal and ecological assessment based on land surface temperature and quantifying multivariate controlling factors in Bogura, Bangladesh
date: 2021-09-17
journal: Heliyon
DOI: 10.1016/j.heliyon.2021.e08012
sha: e5f3bb570f5557ada1c12ef8afbe5dde66de5444
doc_id: 21309
cord_uid: tvaaynb3

In recent years, the world has shown considerable concerns about environmental degradation accompanied by urban expansion. In terms of size, Bogura is equivalent to most of the major cities in Bangladesh, yet no thermal and ecological assessment has ever been conducted here. This study uses multitemporal Landsat satellite images between 2001 and 2020 to investigate the thermal and ecological conditions of Bogura Sadar (sub-district). Land surface temperature (LST) is obtained from Landsat images using the widely used radiative transfer equation. The thermal and ecological conditions are evaluated by computing urban heat island (UHI) and urban thermal field variance index (UTFVI) from LST data. The influence of vegetation, built-area, water-body, and bare soil on LST are examined using land cover indices through pixel-level multivariate linear regression analysis. According to the findings of this sub-district-scale (urban and rural areas) study, the mean LST has increased by 0.62 °C in the last 20 years. As per local administrative-wise findings, LST has increased in most areas, regardless of their urban or rural function. The difference between the urban area and the rest of the surroundings was significant (1.74 °C) in 2020. In 2001, UHI affected area was 5.65 km(2), which expanded to 8.84 km(2) in 2020. Thermal and ecological conditions are worse in urban areas than its surrounding areas. The regression models of the LST and land cover indices could explain more than half (R(2): 0.66 to 0.73) of LST variation over the years. Land cover could explain the LST in 2020 to the least extent implying that anthropogenic activities have greater influence than earlier. Land cover could explain less than half of the LST variation in the urban area.

Rapid urbanization facilitates economic and social progress, but this urban growth is accountable for numerous environmental problems ranging from local to global scale (Cui and Shi, 2012; Kim and Baik, 2005; Nguyen et al., 2019) . Serious complications such as social and environmental degradation, biodiversity loss, destruction of ecosystems, and health disasters emerge from uncontrolled urbanization (Anbalagan, 1993; Dewan and Yamaguchi, 2009; Gazi et al., 2020) . Bangladesh is rapidly urbanizing compared to other South Asian countries, despite being predominately rural (Muzzini and Aparicio, 2013; Rahman et al., 2019) . Many city-scale studies suggest that major cities in Bangladesh are going through uncontrolled urban growth with negative environmental consequences (Ahmed et al., 2013; Fattah et al., 2021a, b; Gazi et al., 2020; Kafy et al., 2020 Kafy et al., , 2021 Naim and Kafy, 2021) . (Rahman et al., 2019) classified 331 cities of Bangladesh using night light intensity and placed Bogura in the same cluster with seven major cities. Therefore, the impact of urban development on the thermal and ecological environment of Bogura needs to be assessed to help the authorities for better environmentally-conscious urban planning.

Rapid urbanization replaces natural surfaces with impervious surfaces, which raises land surface temperature (LST) than surrounding rural areas due to their high heat absorption and retaining capacity (Alexander, 2021; Fan and Sailor, 2005; Oke, 1976; Santiago and Gomes, 2016; Voogt and Oke, 2003; Ziter et al., 2019) . The phenomenon of excess heating of urban and suburban areas is widely known as urban heat island (UHI), which has adverse consequences on the environment (Wemegah et al., 2020) . The formation of UHI depends on the nature of the urban landscape, population density, anthropogenic activity, and meteorological conditions (Grimmond, 2007; Mohajerani et al., 2017) . The UHI effect intensifies air pollutants and greenhouse gases and deteriorates water quality (Phelan et al., 2015; Watkins et al., 2007; Wemegah et al., 2020) . UHI alters the pollination functions of urban ecosystems (Turrini and Knop, 2015) , which has a detrimental impact on urban ecology and biodiversity (Venter et al., 2020) . Urban hotspots are developed by intense UHI, which experience heat stress, excessive energy use, and lower comfortability (Guha et al., 2017; Wemegah et al., 2020) .UHI is also responsible for the deterioration of human health, rise of heat-related illness (heat exhaustion, heat syncope, and heat stroke), fatalities, and spread of vector-borne diseases in tropical cities (Alfraihat et al., 2016; Brousse et al., 2019; Kovats and Hajat, 2008; Moda et al., 2019) .

Urban thermal field variance index (UTFVI) is another widely used indicator derived from LST to evaluate the urban environment's ecoenvironmental quality or thermal well-being (Ahmed, 2018; Alcantara et al., 2019; Alfraihat et al., 2016; Kafy et al., 2021; Maithani et al., 2020; Nguyen et al., 2019; Portela et al., 2020; Sharma et al., 2021; Sobrino and Irakulis, 2020) . Changes in local wind patterns and humidity, declining air quality and comfort, rising mortality rates, rain, and thunderstorm activity are adverse effects caused by high UTFVI Sejati et al., 2019; Singh et al., 2017) .

Satellite provides decent data at local and global scales, which can be easily processed and analyzed utilizing GIS and remote sensing techniques for monitoring (Faridatul and Wu, 2018) . Satellite thermal infrared data are commonly utilized to derive LST with relatively high accuracy and high spatial resolution (Mao et al., 2021) . LST obtained from satellite imagery has been used in numerous studies of UHI and UTFVI due to the advantages of high geographic and temporal resolution, availability, and free and easy access (Ahmed, 2018; Guha et al., 2017; Kafy et al., 2021; Naim and Kafy, 2021; Nguyen et al., 2019; Sharma et al., 2021; Wemegah et al., 2020) . These studies examined the relationship between LST and land surface, which is considered one of the major factors in LST change. In Bangladesh, a number of studies used simple linear regression (Abir and Saha, 2021; Naim and Kafy, 2021) and Pearson's correlation (Ahmed et al., 2013; Kafy et al., 2020; Trotter et al., 2017) to understand the effects of land cover on LST. However, no study has yet to be found explaining the effects of land cover on LST using multivariate linear regression (MLR). MLR model could best fit the distribution, while urban LST fitting based on a single variable could underestimate the distribution Yang et al., 2019) . MLR is often used in recent studies to explain spatial parameters of LST (Unger et al., 2010) .

The primary purpose of this study is to investigate the thermal and ecological conditions of Bogura Sadar, a sub-district of Bogura district, which has not been explored before. The study is vital for better planning and design strategies because bad thermal and ecological conditions negatively impact public health and quality of life. This is the first study in Bangladesh that measures the effect of various land covers on LST using multivariate linear regression. To conduct this sub-district-scale study, 30m spatial resolution image Landsat satellite images are used. Four land cover indices are used to evaluate the spatial and temporal variations in LST induced by land cover. The following are the specific objectives of this study: (1) LST estimation from thermal infrared Landsat data, (2) identifying UHI and Non-UHI areas of Bogura, (3) evaluating the ecological condition of Bogura using UTFVI, (4) finding the correlation of LST with the different land cover indices and (5) exploring the influence of different land cover indices on LST variation.

The present study is focused on the Bogura Sadar, which encompasses the city of Bogura ( Figure 1 ). Bogura thana was founded in 1821 and was elevated to the sub-district in 1983. It is an administrative unit that is part of the Rajshahi division in Bangladesh. Bogura is an important transportation hub since the highway through Bogura connects Dhaka, Rajshahi, and Rangpur. Geographically, the area covers an area of 176.28 square kilometers and is located between latitudes 24 49 0 00 00 N and 24 59 0 00 00 N and longitudes 89 15 0 30 00 E and 88 26 0 30 00 E. The study area is divided into twelve smaller local administrative units (one municipality and eleven unions/rural areas). The information about the area, population, and density of local administrative units are given in Table 1 . The Karatoya River flows through this sub-district to the south, which is a part of the Brahmaputra valley. Bogura falls in the tropical monsoon zone marked by high temperatures, significant rains, and often high humidity (Banglapedia, 2014) . According to the Bangladesh Meteorological Department (BMD) report, the average temperature in the district is 30.6 C, with a peak temperature of 34.3 C in April. The annual total rainfall is 1843 mm, and the average relative humidity is 77 %.

Freely available Landsat 5 (Thematic Mapper) and Landsat 8 (Operational Land Imager and Thermal Infrared Sensor) images are obtained from the USGS archive. The obtained images belong to path code 138 and row code 043. Images of March-April are used in this study which represents the early summer season. Since April has the highest average temperature (BMD, 2021), images close to April are used to get a decent sample of the summer season. Although some scenes have a higher percentage of scene cloud cover, the area of interest is not affected by the cloud. The data values of the two Landsat images are used in this study to avoid inconsistencies. Wide variation of LST exists across different days of the same season in a year hence estimating LST from multiple Landsat images for each year improves the credibility of findings (Abir and Saha, 2021) . Detailed information of the data (year, acquired date, satellite sensor, spatial resolution, and cloud cover) is shown in Table 2 .

The UHI and UTFVI are evaluated using estimated LST from the thermal bands to answer research objectives. Land cover indices are generated to comprehend the change in vegetation, built-area, waterbody, and bare soil. The multivariate linear regression is applied to evaluate the influence of each index on a different part of the area. Details methods are shown in Figure 2 and describing below:

Landsat images should be preprocessed before using for any quantitative analysis to increase the accuracy. Landsat satellite images will always have some distortion due to the sensor, solar, atmospheric, and topography factors (Young et al., 2017) . The USGS Landsat Level 1 satellite data is geometrically and radiometrically corrected, and therefore no additional analysis is performed. Atmospheric correction of the Landsat image is performed using ERDAS IMAGINE before the main operation begins.

Digital numbers (DNs) are used in Landsat sensors to store thermal data (Aik et al., 2020) . DNs are converted to radiance using provided equation in the respective sensor's handbook. Landsat 8 have two thermal band, band 10 (10.6 μm-11.19 μm) and band 11 (11.5 μm-12.51 μm).

Stray light from far out-of-field has affected Landsat 8, and USGS has recommended not to use band 11 (USGS, 2017) . For this reason, the only band 10 TIRS has been used in the present study for LST estimation. There is only one thermal band (band 6) in Landsat 5 used in LST estimation. The underlying model is referred to as the radiative transfer equation (RTE), and it is one of the widely utilized processes for LST retrieval.

DNs are converted to spectral radiance using Eq. (1) for Landsat 5 (USGS, 2019).

Here, L λ ¼ Spectral radiance at the sensor's aperture; LMAX λ ¼ Spectral radiance scaled to QCALMAX; LMIN λ ¼ Spectral radiance scaled to QCALMIN; QCAL ¼ Quantized calibrated pixel value in DN; QCALMAX ¼ Maximum quantized calibrated pixel value; QCALMIN ¼ Minimum quantized calibrated pixel value. DNs are converted to spectral radiance using Eq. (2) for Landsat 8 (USGS, 2016).

Here L λ ¼ Spectral radiance at the sensor's aperture; M L ¼ Radiance multiplicative scaling factor for the band; A .L. ¼ Radiance additive scaling factor for the band; Q cal ¼ Pixel value in DN. The obtained L λ values are used to calculate the brightness temperature. The given equation is subtracted by 273.15 to achieve results in Celsius

Here T ¼ Brightness temperature in Celsius, K 1 and K 2 ¼ calibration constant L λ ¼ spectral radiance at the sensor's aperture. The calibration constants K 1 and K 2 are obtained from satellite data metadata are tabulated below in Table 3 . Depending on the land cover, spectral emissivity (ε) adjustment needs to be made, which is influenced by several elements such as texture, roughness, and structure (Sharma et al., 2021) .

The surface emissivity is an inherent property of land that converts surface heat energy into radiant energy (Sobrino et al., 2004) . The present study used NDVI-based emissivity correction using Eq. (4).

Here, P v is the proportion of vegetation and estimated by Eq. (5)

The ratio of red and near-infrared spectral bands is used to calculate NDVI.

Finally, the LST is calculated using corrected land surface emissivity values using Eq. (7).

Here λ ¼ wavelength of the band; ρ ¼ h* c σ ¼ 6:626* 10 À34 *2:998* 10 8 1:38* 10 À23 ¼ 14380*10 À2 m K; Average LST value of each year retrieved using cell statistics tool in ArcGIS 10.8.

Urban hotspots in urban heat islands are small areas with high temperatures that are uncomfortable for human activity (Sharma et al., 2021) . LST values are used to delineate urban hotspots in the study area.

Here, μ ¼ mean LST and σ ¼ standard deviation of LST. This formula identifies urban hotspots as areas with LST values greater than the mean at the 95 percent confidence interval.

To better assess the presence of UHI in the study area, UTFVI is quantitatively analyzed following an equation used by many researchers (Guha et al., 2017; Naim and Kafy, 2021; Zhang et al., 2006) . According to the ecological evaluation index (EEI), the UTFVI values can be classified into six categories range from excellent to worst, shown in Table 4 ( Naim and Kafy, 2021; Renard et al., 2019; Sharma et al., 2021; Singh et al., 2017) .

Here, T s ¼ LST value of pixels and T m ¼ the Mean LST of the area.

BSI is an improved bare soil index that uses SWIR2 instead of SWIR1 proposed by (Diek et al., 2017) used in the present study for bare land extraction. MNDWI is the modified normalized difference water index proposed by (Xu, 2006) which improved open water features identification while effectively suppressing and even eliminating built-up land, vegetation, and soil noise. SAVI is the soil-adjusted vegetation index proposed by (Huete, 1988) to eliminate the soil background and noise. NDVI can work in areas where plant cover is greater than 30%, while SAVI can work with only 15% plant cover, making it better suited for urban areas (Ray, 2006) . NDBI is widely used for built-up area extraction improved by (Zha et al., 2003) . However (Xu, 2008) , proposed an index-based built-up index (IBI) which some edge over NDBI to extract information of the built-up area. This index extracts complex built-up land cover using remote sensing data derived using index-derived bands such as NDBI, SAVI, and MNDWI. IBI enhances the contrast between vegetation and water bodies and built-ups, thus providing accurate built land information . Table 5 provides detailed information on the bands which are used in the study. 

IBI can be alternatively written as,

3.2.6. Multivariate linear regression Multivariate linear regression models are used in this study to understand the relationships between the land cover variables and LST variation. Instead of using a single variable like simple linear regression, the MLR uses multiple variables to achieve the best fit. MLR can provide a thorough explanation because multiple variables often affect simultaneously. An MLR equation follows as

Here, y is the dependent variable; X 1, X 2 , X 3 , and X 4 are explanatory variables; a 1 , a 2 , a 3 , and a 4 are absolute coefficients that describe the influence factor of the corresponding explanatory variable; and a 0 is constant.

The LST of the Bogura sub-district is estimated by combining two Landsat images for each year. Low and high LST is represented by green to red gradual color scheme in Figure 3 . The visual illustration shows red areas were widely spread in 2001, but they are getting concentrated over time. On the other hand, the magnitude of the green areas reduced in the east and increased in the west. The LST finding showed that the mean was 25.06 C, 25.29 C, and 25.68 C for the years 2001, 2009, and 2020, respectively ( Table 6 ). The lowest LST value (21.72 C) was observed in 2001, whereas the highest LST value (33.19 C) was in 2020. The mean LST of the study area is gradually increasing. Over the span of 20 years, mean LST increased by 0.62 C, which means compared to 2001, the LST increased by 2%. The mean LST increases by around 0.3 C every ten years. LST levels have also risen in Rajshahi, which is located in the same division (Kafy et al., 2020) . Similar findings have also been found in other cities such as Dhaka (Ahmed et al., 2013) , Chattogram (Gazi et al., 2020; Naim and Kafy, 2021) , Khulna (Fattah et al., 2021a, b) .

A suitable validation is essential for any satellite retrieved LST result by in-situ measured data or other satellite-driven measurements (Guha et al., 2017) . Widely available MODIS, AVHRR, ASTER, and Sentinel imageries can be freely and conveniently used to validate LST data. This study applied 1000m spatial resolution MODIS Terra data (MOD11A2) to validate Landsat LST values. The MOD11A2 product is generated by computing the average LST values from the MOD11A1 product of eight days under a clear sky situation. The MODIS satellite sensor has been collecting data since 2000 and is still in operation, making it a suitable instrument for data validation. For both Landsat and MODIS average LST value of two satellite images is considered as the final value in Table 7 .

The study area was categorized into four LST zones to explain further LST distribution (24 C, 24-28 C, 28-32 C, and >32 C). LST distribution pattern over the study area is demonstrated in Table 8 . During 2001, 33.85% of the land was under 24 C, 59.61% was between 24-28 C and 6.53% was in 28-32 C, and a negligible amount of land was in 32 C. During 2020, 13.60% of the total land was under 24 C, 75.82% was between 24-28 C and 10.56% was in 28-32 C, and 0.02% was in 32 C zone. In 2020 larger land areas were fall under in the hottest LST zones than the previous years. Compared to 2001, the amount of land in the 28-32 C LST zone increased by more than 6 km 2 . This pattern also demonstrates that LST in the study area is rising over time.

Spatial change in LST distribution is determined using raster calculation. The visual illustration (Figure 4) shows a larger portion of the area has witnessed an LST increase. The LST was unchanged or decreased in the blue symbol area near the middle region of the study area between 2001 and 2009. Between 2009 and 2020, the blue symbol area was scattered mostly in the northern part of the area. Compared to 2001 in 2020, the northern part experienced an increase of 1 C LST, the centre region witnessed a decline, and the southern part experienced an increase of more than 1 C.

Compare to 2001 in 2009, 28.01% of total land experienced a reduction of LST, 57.26% experienced an increase of 1 C, and 14.25% experienced more than 1 C increase (Table 9 ). Compare to 2009 in 2020, 30.58% of total land experienced a reduction of LST, 51.21% experienced an increase of 1 C, and 17.95% experienced more than 1 C increase. Compare to 2001 in 2020, 25.53% of total land experienced a reduction of LST, 38.58% experienced an increase of 1 C, and 35.63% experienced more than 1 C increase.

The distribution of urban hotspots for the years 2001, 2009, and 2020 is expressed in Figure 5 . Urban hotspots were delineated using the values tabulated in Table 10 . The threshold values are 28.62 C, 28. 16 C, and 28.93 C for 2001, 2009, and 2020 . LST values greater than the threshold values are termed as UHI-affected areas and highlighted in red. LST values smaller than the threshold values are termed as Non-UHI area and highlighted in green. The UHI affected area was 5.65 km 2 , 8.80 km 2, and 8.84 km 2 for 2001, 2009, and 2020. Despite the slight change between 2009 and 2020, results suggest that UHI-affected areas grew steadily between 2001 and 2020. An increase of the UHI effect was also observed in Chattogram, Rajshahi, Cumilla, and Rajshahi (Kafy et al, 2020 Naim and Kafy, 2021) . The urban heat island phenomenon is known as the urban heat archipelago, when the affected regions are sparsely scattered (Bottya'n and Unger, 2003) . Therefore, the hotspot distribution of 2001 can be best described as an urban heat archipelago. In 2001, hotspot pockets were small and distributed over a larger area while pockets were gradually getting fused and concentered in a place 2020.

Urban thermal field variance index (UTFVI) explains the UHI effect in the study area quantitatively. It evaluates the environment's well-being with respect to ecological thermal comfort while taking into account the UHI effect. The visual illustration ( Figure 6 ) express excellent to worst ecological condition represented by green to red gradual color scheme. In 2020, the worst ecological condition concentered in a small area, whereas it was dispersed in 2001. Between 2001 and 2020, the increase in the thermal comfort or improvement of eco-environmental quality is the major characteristics of the study area. In this study, the excellent, good and normal categories are termed "favorable ecological conditions" while the bad, worse, and worst categories are termed "adverse ecological conditions". The area under favorable ecological conditions was 164.76 km 2 in 2001, and it increased to 165 km 2 by 2020 (Table 11) . Over 20 years, 0.24 km 2 of land has been transformed from adverse to favorable conditions. Subsequently, adverse ecological conditions had 11.53 km 2 land coverage in 2001, and it decreased to 11.25 km 2 in 2020. Between 2001 and 2009, a large amount of land converted into favorable to adverse conditions. In 20 years, around 5 percent of total normal ecological land transformed into excellent ecological land. Conversion of fallow and agricultural land into the residential area may contribute to this change. The findings of UTFVI are different from other studies done in Bangladesh, where the worst ecological condition was the dominant attribute Naim and Kafy, 2021) . However, changes occurring in suburban and rural areas were beyond the scope of these city-level studies. On the other hand, the absolute land cover under the worst ecological state increased in this study, which is consistent with earlier studies.

The sub-district of Bogura is divided into twelve local administrative units (one municipality and eleven unions), each with its own governing body. Bogura Pourosova is a major urban area located in the midsection of the region's south. Shabgram, Nishindra, and Fapore are three adjacent unions of Bogura Pourosova situated from east to west shown in Figure 1 . Other unions located opposite sides of the sub-district predominately function as rural areas. Administrative unit-wise mean LST distribution and their change from 2001 to 2020 is shown in Table 12 . During 2001, the highest mean LST was observed in Bogura Pourosova (27.89 C), and the lowest LST was in Namuja (24.18 C). Compared to the secondhighest area (Nishindara), LST in Bogura Pourosova was higher by 1.12 C. During 2020, the highest mean LST was observed in Bogura Pourosova (29.12 C), and the lowest LST was in Noongola (24.57 C). Compared to the second-highest area (Nishindara), LST in Bogura Pourosova was higher by 1.74 C. In 2001, 2009, and 2020, the mean LST in Bogura Pourosova was 2.80 C, 2.96 C, and 3.52 C higher than the rest of the area. Over the 20 years, the highest increase of LST was observed in Fapore (1.63 C), followed by Shabgram (1.26 C) and Bogura Pourosova (1.23 C). Shakharia is the only area that observed a reduction of LST (-0.52 C).

Administrative unit-wise distribution and changes of UHI affected areas and adverse ecological condition areas from 2001 to 2020 are shown in Table 13 . During 2001, the highest amount of UHI affected area was observed in Nishindara (1.62 km 2 ), followed by Bogura Pauroshava (1.40 km 2 ) and Shakharia (1.02 km 2 ). During 2020, the highest amount of UHI affected area was observed in Bogura Pourosova (3.14 km 2 ), followed by Nishindara (2.34 km 2 ) and Fapore (1.51 km 2 ). Increase of UHI of observed in Bogura Pourosova > Fapore > Shabgram > Nishindara > Erulia > Gokul. Decrease of UHI of observed in Shakharia > Rajapur > Sekherkola > Lahiri Para > Noongola. About 51% of Bogura Pourosova was under adverse ecological conditions in 2001, and it increased to 68% in 2020. More than 20% land of Bogura Pourosova, Nishindara, and Shakharia were under adverse ecological conditions in 2001, while in 2020, Bogura Pourosova and Nishindara met this criterion. Ecological conditions in most locations have remained unchanged or improved during the last 20 years. While the ecological situation in Shakharia improved the most, it deteriorated in Bogura Pourosova, Erulia, Fapore, and Shabgram.

The administrative units are classified into five groups using mean LST and adverse ecological land coverage in 2020 (Figure 7) . Mean LST shows Bogura Pourosova and Nishindara had the highest LST intensity. Adjacent unions Fapore, Rajapur, Shabgram had lesser LST intensity. LST intensity reduced gradually from south to north. Ecological land coverage shows Bogura Pourosova and its adjacent unions are affected by adverse ecological conditions. All the other areas are under good ecological condition. The most densely populated Bogura Pourosova had the highest intensity of LST and worst ecological condition. The second most densely populated is also associated with intense LST and worse ecological conditions than low, dense areas. 

Land cover indices SAVI, IBI, MNDWI, and BSI, are generated using spectral bands to quantify vegetation, built-up, water, and bare-soil land cover for the study area. SAVI, IBI, MNDWI, and BSI are illustrated in green, red, blue, and yellow for the years 2001, 2009, and 2020 in Figure 8 .

The pattern of the correlation between LST and land cover indices in the whole study area is tabulated in Table 14 . Over time, the strength of the correlation between the land surface and LST has changed slightly. SAVI, IBI, and BSI strongly correlated with LST, while MNDWI had a moderate correlation. IBI and BSI had a positive correlation with LST, while SAVI and MNDWI had a negative correlation. Similar kinds of relationships were found in numerous studies. Built-in areas and bare lands raise LST and contribute to UHI, whereas vegetation and water lower LST (Amiri et al., 2009; Guha et al., 2017; Song et al., 2014) .

The multivariate linear regression model of Eq. (16) provides a thorough explanation of LST variation as follows-

Here, y is the dependent variable (LST); X 1, X 2 , X 3 , and X 4 are numerical values of explanatory variables (SAVI, IBI, BSI, and MNDWI); and a 1 , a 2 , a 3 , and a 4 are absolute coefficients of the corresponding explanatory variable. The positive and negative coefficients of the explanatory variable are responsible for the increase or decrease of LST. The coefficient indicates the magnitude of LST changes when the dependent variable is changed by one unit. The coefficient of determination (R 2 ) demonstrates the model's explanatory power. All pixel values (more than 190000) of the study area are used in the statistical analysis using SPSS software. For simplicity, SAVI, IBI, BSI, and MNDWI are represented in the regression equations as V, B, S, and W. During different stages of development, different LST variables may vary. A stepwise MLR analysis is performed to understand influential land covers for the years 2001, 2009, and 2020 in Table 15 . Stepwise regression is an automatic computational process that helps to generate the "best" multivariate regression model by considering only statistically significant variables from a larger number of potential predictive variables. In 2001, the built-up land cover could explain 71% of the LST variation alone. Adding vegetation cover to the model increases the model's ability to explain the variation by 2%. In 2009, vegetation cover alone explained 70%, Whereas the model could explain 73% of LST variation when all four land covers were used. In 2020, the bare-soil cover alone explained 63%, Whereas the model could explain 66% of LST variation when used all four land covers. The explanatory power of land cover for LST variation decreased in 2020, indicating that other anthropogenic activities play a larger role in 2020 than in previous years. Table 16 shows LST regression equations for different years using four land covers in the model to understand each land cover's influence better. The standardized coefficient (β) was considered significant when the confidence level was p at level 0.01. In 2001, the most dominating heating and cooling factors were built-up cover and vegetation cover, respectively. The standardized coefficient (β) was 0.40, which means for one unit increase in built-up, LST will increase by 0.40 units. Similarly, for one unit of vegetation cover, LST will decrease by 0. 45 units. The built-up was the most dominant heating land cover in 2009, increasing LST by 0.35 units for one unit. Whereas vegetation was the most dominant cooling land cover, decreasing LST by 1.06 units for one unit. The coefficient values are different for the year 2020, where all the variables are assigned with negative values. However, the coefficient can be compared with the coefficient of other variables. The most dominating heating and cooling factor was bare-soil and vegetation cover, respectively. Despite the fact that the LST reduced by 0.33 units for each unit of bare-soil yet represents the peak LST value. While vegetation was the most dominant cooling land cover, decreasing LST by 3.02 units for one unit.

Each local administrative unit's main characteristics and functions differ; therefore, their land cover type also differs. Administrative unitwise MLR analysis is performed for 2020 using four land covers to understand the impact of land covers on each administrative unit in Table 17 . The regression model could explain 80% of the LST in Fapore, which is the best among the twelve administrative units. Land cover pattern in Nishindara, Lahiri Para, Gokul, Shekherkola, Shabgram, Namuja, and Nungola could explain more than half of the variation in LST. Only 41% of the LST variation could be explained by land cover in Bogura Pourosova, implying other factors have greater influence in urban areas. Urban geometry (urban canyon, built-up height, and street width) increases LST by trapping heat inside the compact structure and blocking wind circulation (Mohamed et al., 2021; Wang and Xu, 2021) . Anthropogenic activities (population density, economic activity, industrial infrastructure, trip frequency, and travel distance) cause changes in LST in urban and suburban areas (Qiao et al., 2014) . The bare-soil cover was the dominant factor in the increase of LST for most of the administrative area. In Noongola, the built-up cover has a more significant influence on LST increment than bare-soil, increasing LST by 1.15 units for one unit. The built-up cover had an insignificant influence on three unions, namely Erulia, Lahiri Para, and Namuja located at the edge of the study area. In most of the administrative units, vegetation had a greater cooling effect. The cooling effect of vegetation was less than the built-up area in Nishindara, which is not desirable. The cooling effect of water had an insignificant influence in Erulia, Lahiri Para, and Sekherkola. However, the cooling effect of water was less than built-up in Bogura Pourosova, Gokul, Nishindara, and Shabgram, which is not desirable. This information will be helpful for policymakers, planners, and governing authorities for making better development and management plans.

The findings of the study can help in the decision-making of relevant agencies for the improvement of ecological environments. However, this study does not account for the influence of urban geometry and anthropogenic activity in LST variation, which is one of the study's major limitations. Urban geometry and anthropogenic activities such as city centre, high-rise built-ups, urban density, street pattern, street connectivity, street centrality, and land use significantly influence LST (Erdem et al., 2021; Parvez et al., 2021) . Furthermore, the UHI effect is more intense at night (Basara et al., 2010; Oke, 1988) , but Landsat collects data in the morning, which is a downside. Additional research is required in the future to strengthen the understanding. LST obtained from another satellite image or in-situ data with better spatial and temporal resolution can improve the accuracy of the results. The effects of urban geometry and anthropogenic activities on LST should be investigated to understand the urban thermal environment.

Multitemporal Landsat satellite images are used in this study to investigate the thermal and ecological condition of the Bogura subdistrict for 2001, 2009, and 2020. The UHI effect and ecological condition are quantified through urban hotspots and UTFVI value based on LST. Quantitative models describing the land cover and LST relationship are constructed using multivariate linear regression analysis considering each pixel value. The results demonstrate that LST in most regions has increased regardless of their urban and rural function. In 2001, the distribution of urban heat hotspots was similar to the archipelago, which in 2020 resembled islands. The total area covered by UHI remained almost the same in 2009 and 2020, indicating that restrictive factors limit higher-degree of urban transformation. High UHI intensity and adverse ecological conditions are concentrated in the southern side, where a higher degree of urban agglomeration exists. Ecological condition is undergoing an improvement in the northern part of the region designated as rural. Land-use change caused by urbanization and industrialization worsens the ecological condition at a relatively high pace, which greatly affects ecosystem service quality. The regression models of the LST and land cover indices (SAVI, IBI, BSI, and MNDWI) could explain more than half of LST variation over the years. In 2001 and 2009, the built-up land cover had a more significant influence than vacant land in LST rise, while bare-soil had a more significant influence in 2020. On the other hand, vegetation consistently remained the dominant land cover with the most cooling effect on LST. Administrative unit-wise regression models of 2020 revealed that land cover could explain the LST variation in Bogura Pourosova to the least extent, implying that urban geometry and anthropogenic activities have greater influence in urban areas. Land cover could explain the LST variation to a greater extent in other areas, implying that land cover has greater control of LST in rural areas. The study's knowledge of the variables influencing LST will assist authorities in translating these aspects into policies to reduce the UHI effect. Green strategies need to be adopted to mitigate UHI through the tree canopy system, urban forests, street trees, green roofing, and vertical greenery. Built-up and construction materials of artificial structures need to be changed with low heat capacity and high albedo materials. Water-body preserve policies need to be developed depending on type and geometries.

Farhan Asaf Abir: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Sabbir Ahmmed, Soykot Hossain Sarker: Analyzed and interpreted the data; Wrote the paper.

Ashraf Uddin Fahim: Contributed reagents, materials, analysis tools or data.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. 

Data will be made available on request.

The authors declare no conflict of interest.

No additional information is available for this paper.

Assessment of land surface temperature and land cover variability during winter: a spatio-temporal analysis of Pabna municipality in Bangladesh

Simulating land cover changes and their impacts on land surface temperature in Dhaka

Assessment of urban heat islands and impact of climate change on socioeconomic over Suez Governorate using remote sensing and GIS techniques

Land use/land cover changes and the relationship with land surface temperature using Landsat and MODIS imageries in cameron highlands

Geospatial assessment and modeling of urban heat islands in Quezon City, Philippines using ols and geographically weighted regression

Influence of the proportion, height and proximity of vegetation and built-ups on urban land surface temperature

Ecological evaluation of urban heat island in chicago city

Spatial-temporal dynamics of land surface temperature in relation to fractional vegetation cover and land use/cover in the Tabriz urban area

Environmental hazards of unplanned urbanization of mountainous terrains: a case study of a Himalayan town

Climate -Banglapedia

The impact of the urban heat island during an intense heat wave in Oklahoma city

Monthly Maximum Temperature | Bangladesh Meteorological Department

A multiple linear statistical model for estimating the mean maximum urban heat island

Using local climate zones in subsaharan africa to tackle urban health issues

Urbanization and its environmental effects in Shanghai

Land use and land cover change in Greater Dhaka, Bangladesh: using remote sensing to promote sustainable urbanization

Barest pixel composite for agricultural areas using Landsat time series

An analysis of urban form factors driving Urban Heat Island: the case of Izmir

Modeling the impacts of anthropogenic heating on the urban climate of Philadelphia: a comparison of implementations in two PBL schemes

Automatic classification of major urban land covers based on novel spectral indices

Multi-layer perceptron-Markov chain-based artificial neural network for modelling future land-specific carbon emission pattern and its influences on surface temperature

Spatio-temporal dynamic land cover changes and their impacts on the urban thermal environment in the Chittagong metropolitan area

Urbanization and global environmental change: local effects of urban warming

Dynamic analysis and ecological evaluation of urban heat islands in Raipur city

A soil-adjusted vegetation index (SAVI)

Prediction of seasonal urban thermal field variance index using machine learning algorithms in Cumilla

Modelling future land use land cover changes and their impacts on land surface temperatures in Rajshahi

Spatial and temporal structure of the urban heat island in seoul

Heat stress and public health: a critical review

Investigating the effect of lockdown during COVID-19 on land surface temperature: study of dehradun city

Resolution enhancement of remotely sensed land surface temperature: current status and perspectives

Impacts of climate change on outdoor workers and their safety: some research priorities

The urban heat island effect, its causes, and mitigation, with reference to the thermal properties of asphalt concrete

Urban heat island effects on megacities in desert environments using spatial network analysis and remote sensing data: a case study from western Saudi arabia

Bangladesh, Directions in Development -General

Assessment of urban thermal field variance index and defining the relationship between land cover and surface temperature in Chattogram city: a remote sensing and statistical approach

The urban energy balance

The distinction between canopy and boundary-layer urban heat islands

The influence of urban form on the spatiotemporal variations in land surface temperature in an arid coastal city

Urban heat island: mechanisms, implications, and possible remedies

Impact of urban and industrial features on land surface temperature: evidences from satellite thermal indices

Influences of urban expansion on urban heat island in beijing during 1989-2010

Classification of cities in Bangladesh based on remote sensing derived spatial characteristics

Vegetation in Remote Sensing FAQs in ER Mapper Applications. ER Mapper

Evaluation of the effect of urban redevelopment on surface urban heat islands

Heat islands in the city of macei o/AL using orbital data from Landsat 5

The spatio-temporal trends of urban growth and surface urban heat islands over two decades in the Semarang Metropolitan Region

Assessing urban heat islands and thermal comfort in Noida City using geospatial technology. Urban Climate 35

Impact of Land Use Change and Urbanization on Urban Heat Island in Lucknow City, Central India. A Remote Sensing Based Estimate

A methodology for comparing the surface urban heat island in selected urban agglomerations around the world from sentinel-3 SLSTR data

Land surface temperature retrieval from LANDSAT TM 5

The relationships between landscape compositions and land surface temperature: quantifying their resolution sensitivity with spatial regression models

Effects of rapid urbanization on the urban thermal environment between 1990 and 2011 in Dhaka Megacity

A landscape ecology approach identifies important drivers of urban biodiversity

Modeling of the Urban Heat Island Pattern Based on the Relationship Between Surface and Air Temperatures undefined

Landsat 7 (L7) Data Users Handbook

Landsat 8 OLI and TIRS Calibration Notices [WWW Document

Landsat 8 (L8) Data Users Handbook

Hyperlocal mapping of urban air temperature using remote sensing and crowdsourced weather data

Thermal remote sensing of urban climates

The impact of built-up height on urban thermal environment in summer: a case study of Chinese megacities

Spatial distribution and influencing factors on urban land surface temperature of twelve megacities in China from

Increased Temperature and Intensification of the Urban Heat Island: Implications for Human Comfort and Urban Design

Assessment of urban heat island warming in the greater accra region

A new index for delineating built-up land features in satellite imagery

Modification of normalized difference water index (NDWI) to enhance open water features in remotely sensed imagery

Local climate zone ventilation and urban land surface temperatures: towards a performance-based and wind-sensitive planning proposal in megacities

A survival guide to Landsat preprocessing

Use of normalized difference built-up index in automatically mapping urban areas from TM imagery

Land surface temperature retrieval from CBERS-02 IRMSS thermal infrared data and its applications in quantitative analysis of urban heat island effect

Scale-dependent interactions between tree canopy cover and impervious surfaces reduce daytime urban heat during summer

The authors wish to express their gratitude to the "Institute of Market and Geospatial Research" for supporting the GIS analysis. The authors are grateful to the USGS for making Landsat and MODIS data publicly accessible. The authors appreciate the reviewers for their insightful suggestions that enhanced the quality of the manuscript.