A Time Series Analysis of Forest Cover and Land Surface Temperature Change Over Dudpukuria-Dhopachari Wildlife Sanctuary Using Landsat Imagery

Forest cover change is an important criterion as it affects the environmental balance whereas land surface temperature is a significant parameter within the earth climate system. Spatio-temporal change of forest cover can be detected and land surface temperature can be retrieved by applying remote sensing technology. The present study aimed to capture the impact of forest cover change on land surface temperature in Dudpukuria-Dhopachari Wildlife Sanctuary (DDWS), Bangladesh, using multi-spectral and multi-temporal satellite data. To avoid the biasness in the calculation, leaf flash time was targeted for collecting Landsat images from United States Geological Survey (USGS) Earth Explorer and, based on availability, images were collected purposively which ones had closer time period:1990 (March 5, 1990), 2000 (February 5, 2000), 2010 (February 24, 2010) and 2020 (March 23, 2020). Unsupervised classification was applied over the images Landsat 4–5 Thematic Mapper (TM), 7 Enhanced Thematic Mapper Plus (ETM+), and 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) data for detecting forest cover change. To retrieve the land surface temperature, Mono Window Algorithm (MWA) method was applied over similar images. Maximum forest degradation was observed in 2010 and the change found was 17% as compared to 1990. After 2010, the forest started to flourish. Land surface temperature dramatically changes over the time period. The highest land surface temperature in the forested area was observed in 2020 (32.2°C) and it was changed 7.7°C from that of the 1990 (24.5°C). In every 10 years, almost 2.3°C–3.0°C temperature change was detected. In the first three decades, a reverse relationship was observed between land surface temperature and forest cover; however, in the last decade, land surface temperature was found to increase with the increase of forest cover. Thus, the results of the study revealed that land surface temperature may not be relevant with the local forest cover change directly. It can be estimated from the results that local forest cover change may have limited impact on local temperature rather than global forest cover change, whereas global warming could play a vital role in changing land surface temperature locally as well as globally.


INTRODUCTION
Remote sensing is an important tool to analyse the spatial and temporal changes of land coverage (Petropoulos et al., 2015;Chew et al., 2016;Sekertekin and Bonafoni, 2020). This is a potential technology for inventory and monitoring of world's natural resources (Goetz and Dubayah, 2011;Sharma et al., 2021). Forest is an important natural component. It maintains the environmental and ecological balance (Liu et al., 2020;Yuan et al., 2020) by preventing unwanted climate change and keeping ecosystems balance between human beings, plants, animals, and other abiotic components. Satellite imagery is a useful unbiased space borne photography as it detects the changes of land coverage in a periodic repetitive manner. Thus, forest cover changes and damages are eventually assessed by applying remote sensing technology (Zeng et al., 2020).
Remote sensing technology is an important new source of data that offers accurate geo-referencing information for harvest and analysis with comparatively low-cost but faster precise results (Lein, 2011). Among other branches, thermal remote sensing is the one that deals with Thermal Infrared (TIR) data of the Electromagnetic (EM) spectrum (Al-doski et al., 2016;Messina and Modica, 2020). The ground emits radiation and the thermal remote sensing captures it. The land surface temperature represents the temperature of the earth's surface and is estimated from radiance measurements (Sekertekin and Bonafoni, 2020;Wan Mohd Jaafar et al., 2020). Remotely sensed thermal infrared data analyse the large scale temporal and spatial land surface temperature (LST) from local as well as global levels.
Forests in developing countries are under greater risk as the per capita land is limited here (Fagan et al., 2020;Huang et al., 2020). Agricultural expansion in the forest land and illegal harvesting from the natural forest area are the major drivers of deforestation (Islam and Sato, 2012;Kissinger et al., 2012;Hishe et al., 2021). Forest degradation in an alarming rate is a global issue now (Halme et al., 2013;Morales-Hidalgo et al., 2015;Hite and Seitz, 2021;Wagner et al., 2021).
Bangladesh is a small developing country with massive population that puts the forest coverage under greater dangers Van Schendel, 2020). Forest coverage here is as an open treasure to collect resources. To protect these resources, the government of Bangladesh has declared different forest areas as protected area at different times (Chowdhury and Koike, 2010;Chowdhury, 2014;Mukul et al., 2016;Rahman and Islam, 2021). In this special consideration, there are 49 types of protected areas under the broader classification of national park, wildlife sanctuary, special biodiversity conservation area, marine protected area, vulture safe zone, botanical garden, safari park, eco-park, and aviary park (Rahman et al., 2016;BFD, 2021).
Protected area is a clearly defined marine and terrestrial conservation place, recognised, dedicated, and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values (Dudley et al., 2010;IUCN and UNEP-WCMC, 2014;Mitchell et al., 2018). Protected areas are very important as they help to protect species from extinction, boost carbon storage, maintain regular weather condition, and provide natural habitats for plants and animals (Mackey et al., 2008;Stolton et al., 2015;Xu et al., 2017). Thus, it is important to analyze Landsat data with land cover classification methods and LST retrieval methods over protected areas.
Several studies have been done by different researchers on land cover classification (Rodriguez-Galiano et al., 2012;Cetin and Sevik, 2020;Chowdhury et al., 2020;Kaya and Görgün, 2020;MohanRajan et al., 2020;Munthali et al., 2020;Olorunfemi et al., 2020) and LST retrieval (Kayet et al., 2016;Sekertekin and Bonafoni, 2020;Thakur et al., 2021;Rashid et al., 2021). Land cover changes and resulting land surface temperature in Bangladesh have been studied sporadically by various researchers (Dewan and Corner, 2012;Ahmed et al., 2013;Chaudhuri and Mishra, 2016;Ullah et al., 2019;Gazi et al., 2020;Kafy et al., 2020). In this continuation, the GIS and RS techniques were applied in the present study to detect the forest damages and land surface temperature changes, and to invent their interrelationship over Dudpukuria-Dhopachari Wildlife Sanctuary (DDWS) since 1990 to 2020 as no research have been done on land cover classification and LST retrieval over this protected area.

Study Location
Dudpukuria-Dhopachari Wildlife Sanctuary (DDWS) is situated along the borderline of Chittagong, Rangamati, and Bandarban districts, Bangladesh. The area extends from N 22 • 09 to N 22 • 22 and E 92 • 05 to E 92 • 10 (Shahadat et al., 2014;Hossain et al., 2017). This evergreen and semi-evergreen tropical natural forest was declared as a protected area of Bangladesh in 2010 and covers 4,717 ha of land area (Shahadat et al., 2014;BFD, 2021). This vital hilly forest ecosystem also serves as a watershed of two significant rivers of Bangladesh: Karnaphuli and Sangu (Shahadat et al., 2014). Dudpukuria-Dhopachari Wildlife Sanctuary also preserves 386 tons of CO 2 per hectare and is the harbor of 600 plant species and 400 vertebrates including 200 birds and 45 mammals. This tropical hilly forest is generally dominated by two tree species Dipterocarpus turbinatus (Garjan) and Artocarpus chama (Chapalish) (Hassan et al., 2011). The forest is also the herd of globally endangered Elephas maximus (Asian Elephant), Hoolock hoolock (Western Hoolock Gibbon), and globally vulnerable Arctictis binturong (Binturong) (Karim and Ahsan, 2016). People around this wildlife sanctuary are living in 28 villages, five of which are ethnic villages. These people work with Bangladesh Forest Department (BFD) to protect the forest from destruction in a co-managed manner. They help to expand the forest cover area by doing afforestation and assisting natural regeneration (ANR). These people also keep away intruders, illegal trespassers, and hunters by regular patrolling in the forest.

Data Collection
Landsat images were collected from the website of United States Geological Survey (USGS) Earth Explorer. In the study, forest cover changes were detected from 1990 to 2020. Landsat 4-5 Thematic Mapper (TM) data sets were downloaded and used for analysing the forest coverage of 1990 and 2010. For the years 2000 and 2020, Landsat 7 Enhanced Thematic Mapper Plus (ETM+) image and Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) images were analysed correspondingly. To avoid bias, the imageries were downloaded in the similar available time period. Mostly the immediate time around new leaf flushing was preferred for image downloading. Hence, the downloaded images were deliberately chosen and acquired in February and March based on image availability. In order to get the best results, image data were collected on March 5, 1990, February 5, 2000, February 24, 2010, and March 23, 2020 Classification After download, the images were pre-processed and then analysed using ArcGIS 10.3 software. In order to process and classify smoothly, a mask of the forested area was obtained from the image. To understand the severity of forest degradation, the image was broadly classified into only two land coverage, namely, natural forest coverage, and cultivated and bare land. The homesteads were too sporadic and scattered on the agricultural land or under the natural forest coverage. Thus, all homesteads did not cover at least three pixels of the image. That is why the supervised classification was avoided, as it needed at least three pixels for count. A canal-type river was found in the natural forest area and some sporadic ponds were detected throughout the villages. All the waterbodies get dry during February and March each year. Calculation of the area of waterbodies could be biased in this time period. Moreover, the present study tends to find out the forest cover change and the river situated within the forested area. In order to avoid misinterpretation and for proper detection of natural forest area changes, minor classified groups are intentionally evaded. Thus, the settlements were excluded from counting and waterbodies were also deliberately included in the broader natural forest class. In the unsupervised classification, the base map was used for proper analysis and verification. The classification process has been outlined in Figure 1.

Retrieving Land Surface Temperature
Several methods are applicable in retrieving LST. Among all methods, four are mentionable: Mono Window Algorithm (MWA) (Qin et al., 2001;Kumari et al., 2020), Single Channel Algorithm (SCA) (Dong et al., 2017), Radiative Transfer Remote Sensing 2020, 12, 294 5 of 32 Equation (RTE) method (Saratoon et al., 2013), and Split Window Algorithm (SWA) (Rozenstein et al., 2014;Wang et al., 2019) using Landsat 5 Thematic Mapper (TM), 7 Enhanced Thematic Mapper Plus (ETM+), and 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) data (García-Santos et al., 2018;Hao et al., 2019). In this FIGURE 1 | Flow chart of land use and land cover classification process in Landsat 5, 7, and 8. study, the MWA method (Qin et al., 2001;Guha et al., 2019;Guha and Govil, 2020;Kumari et al., 2020) was chosen as it is easily applicable to Landsat 5 TM, 7 ETM+ and 8 OLI/TIRS data. Mono-window algorithm method is different from SWA method and more suitable to apply as it retrieves LST from only one thermal band data. Moreover, compared to other methods, where several atmospheric parameters are required for LST retrieval, the mono-window algorithm requires only two: transmittance and mean atmospheric temperature parameters (Qin et al., 2001). This method has the following steps to follow. After retrieving land surface temperature from each images it was classified into five classes.

Conversion of the Digital Number to Top of Atmosphere Spectral Radiance
Top of Atmosphere (TOA) Reflectance is a unitless measurement that provides the radiation ratio reflected from the incident solar radiation on a given surface. The Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) sensors capture reflected solar energy. These data are then converted to radiance and finally, rescale this data into an 8-bit digital number (DN) with a range between 0 and 255 (Qin et al., 2001;USGS, 2021). Top of Atmosphere can be calculated by using the equations (1a) and (1b): For Landsat 8 OLI/TIRS where Lλ = Top of atmosphere spectral radiance [Watts/ (m 2 × sr × µm)], ML = Radiance multiplicative Band (No.), AL = Radiance Add Band (No.), Qcal = Quantized and calibrated standard product pixel values (DN), Oi = correction value for Band 10 is 0.29.
For Landsat 5 TM and 7 ETM+

Conversion of Spectral Radiance to At-Satellite Brightness Temperatures
Thermal band data can be converted from spectral radiance to top of atmosphere brightness temperature using the thermal constants K1 and K2 (USGS, 2021). Equation (2) can simplify the whole process: where BT = Top of atmosphere brightness temperature (

Normalized Different Variance Index
Normalized Different Variance Index (NDVI) measures surface reflectance and widely used to estimate vegetation growth and biomass (Bhandari et al., 2012;Wu et al., 2016). It is calculated as the ratio between TOA reflectance of a red band and a nearinfrared (NIR) band (Viana et al., 2019). NDVI can easily be clarified in the equation (3): where RED = DN values from the Red band, NIR = DN values from Near-Infrared band.

Land Surface Emissivity
Land surface emissivity (LSE) is an important surface parameter (Li et al., 2013). It is the average emissivity of an earth surface element and is calculated from the emitted radiance measured from space and land surface temperature (Norman and Becker, 1995;Pal and Ziaul, 2017). The algorithms can be explained in detail in equations (4) and (5): where: P υ = Proportion of Vegetation, NDVI = DN values from NDVI Image, NDVImin = Minimum DN values from NDVI Image, NDVImax = Maximum DN values from NDVI Image.
where ε = Land Surface Emissivity, P υ = Proportion of Vegetation, 0.004 and 0.986 correspond to correction values of the equation.

Land Surface Temperature
Land surface temperature (LST) is an important reference index used to measure the radiative skin temperature of the land derived from infrared radiation. It has impacts on energy cycles, ecological system balance, and human life (Yan et al., 2020). Procedure of retrieving LST can be elaborated in Equation (6): where BT = Top of atmosphere brightness temperature

Accuracy Assessment
Accuracy assessment is one of the best ways to determine the quality of information derived from remotely sensed data, and to remove confusion about the predicted classification with the actual one (Yusuf et al., 2014;Zhang et al., 2016). In this process, land cover classification results are compared with the geospatial data. In the present study, resulted land cover classes were linked with Google Earth maps to test the mistakes. The pixels derived from the classified map image were compared with the same location in the geospatial map. Thus, accuracy assessment provides an overall accuracy of map as well as accuracy for each class in the map. The overall accuracy percentage was calculated by using following formula: Overall Accuracy = Total Number of Correctly Classified Pixels (Diagonal) Total Number of Reference Pixels

×100
A total of 118 sample sites from Google Earth were selected for accuracy assessment. In order to calculate the class-based Hasnat Analysis of LULC and LST Over DDWS accuracy besides overall accuracy, accuracy of individual classes was calculated for each class in each map. The user's accuracy and producer's accuracy were two approaches applied. Producer's accuracy is the ratio of the number of correctly classified pixels in each category to the total number of reference pixels in that category. It takes into account the error of omission that refers the real land cover type that were erroneously omitted from the classified map. User's accuracy is calculated by dividing the number of correctly classified pixels in each category by the total number of classified pixels in that category (Pal and Ziaul, 2017). Producer's accuracy and user's accuracy are calculated by using the following formula: Producer Accuracy

×100
Besides overall accuracy, Kappa coefficient (K) is also used as an alternative measurement proposed by Verma et al. (2020). Kappa coefficient (K) was calculated by multiplying the total correctly classified pixels in the verification classes by the total pixels, from which subtracting the sum of the ground verification pixels in a class (column total) times the sum of the classified pixels in that class (row total) summed over all classes (sum of row total × column total), and dividing by the sum of the ground verification pixels in a class (column total) times the sum of the classified pixels in that class (row total) summed over all classes (sum of row total × column total) subtracting from the square of the total pixels (Pal and Ziaul, 2017). The value of Kappa generally lies between 0 and 1, where 0 denotes agreement due to chance only and 1 denotes complete agreement between the two data sets. If negative values occur, they are spurious (Congalton, 1991). Kappa coefficient (K) is formally expressed as a percentage (%). A true Kappa value is always less than or equal to 1. The more the value means, the more the accuracy.
Value equal to 1 represents the perfect agreement. More precisely, poor or very poor agreement between the images determined when Kappa value is less than 0.4, fair agreement determined when values between 0.4 and 0.55, values from 0.55 to 0.70 represent good agreement, values from 0.70 to 0.85 represent very good agreement, and values higher than 0.85 represent excellent agreement (Kikon et al., 2016;Zhang et al., 2016;Pal and Ziaul, 2017;Hua and Ping, 2018). The calculation of Kappa statistic is as follows: Where, TP = Total Pixels, and TCP = Total Corrected Pixels.

RESULTS
The forest coverage damage was analysed and detected from the Landsat images at 10 years interval since 1990 to 2020 (Figure 2). Images were classified into two broad classes: natural forest coverage, and cultivated and bare land. In 1990, the natural forest coverage was 85% which declined to 78% in 2000, and 68% in 2010 (Figures 2, 3). The cultivated and bare land increases with time. It is observed that maximum forest damages occurred during 2000-2010 which was lately covered by natural regeneration and with the help of some plantation by Bangladesh Forest Department (BFD) and United States Agency for International Development (USAID) funded projects, i.e., Integrated Protected Area Co-management (IPAC) and Climate Resilience Ecosystems and Livelihoods (CREL) projects by participating local people from 2010 to 2020.

Hasnat
Analysis of LULC and LST Over DDWS    The classification of the Landsat images for the years 1990, 2000, 2010, and 2020 revealed that the wildlife sanctuary area has been experiencing declining trends from 1990 to 2010 and dramatically changed after 2010. Since 2010 to 2020, the forest coverage had been increased and the cultivated and bare land had been decreased in the area (Figures 2, 3).
The analysis of forest coverage revealed a dynamic change in the forested area. Figure 2 shows the significant variation of forest cover changes over the study period. The natural forest coverage indicates a declining rate at first and then increased, while the cultivated and bare land increased previously with the decline of forested area and decreased when the forest coverage improved (Figure 3).
In 1990, the highest forest coverage (20938.72 acre) and lowest uncovered land (3722.90 acre) observed (Figure 3). In 2010, the lowest forest coverage (16883.95 acre) and highest forest degraded and converted land (7773.62 acre) found after analysis. From 1990 to 2000 an amount of 1787.05 acre forest land has been converted to other land uses. During 2000-2010, total converted forest land was 2267.72 acre. Since 2010 to 2020, the scenario is completely reversed. It is observed by analysing Figure 3 that a total of 1514.50 acres of bare land has been converted to forest land and increased the total forest coverage area from 16883.95 acre (2010) to 18398.45 acre (2020).
The overall accuracy and Kappa coefficients for land use and land cover maps for 1990, 2000, 2010, and 2020 were compared  (Table 1). Concomitantly, the producer's accuracy and user's accuracy for all the classes from 1990 to 2020 exceeded 90% which ensures the highest accuracy at the time when classification was done ( Table 1).
With the time period, the land surface temperature has changed significantly. It may be due to two reasons: change in forest coverage area and global warming. Global warming was considered as the most important factor in this study as after 2010 though the forest coverage increased, the temperature also increased at the same time rather than decreasing. For avoiding the biasness of analysing land surface temperature, the images were collected between the time period February 5 and March 23 based on availability of required images in the website. It was found that the temperature increased from 2.3 to 3.0 • C in each decade (Figures 4-6) which is a great threat to any country.
In 1990(March 5, 1990, the highest temperature found was 24.5 • C and the lowest was 17.9 • C. In 2000 (February 5, 2000), maximum temperature was 26.9 • C and minimum was 18.7 • C. During 2010 (February 24, 2010(March 23, 2020, the uppermost value of temperature was analysed 29.2 • C and 32.2 • C, whereas the lowermost was 21.1 • C and 21.6 • C, respectively (Table 2 and Figure 6). Mean value of temperature show a gradual upward trend (Table 2 and Figure 5). Overall temperature raised 7.7 • C from 1990 to 2020 and changes found a rise of 2.3-3 • C in each decade (Figure 7). The lowest temperature in 1990 (March 5, 1990) was 17.9 • C and in 2020 (March 23, 2020), it was 21.6 • C. In 30 years, a drastic change in temperature occurred. The highest value of temperature has changed 7.7 • C valuing from 24.5 • C in 1990 to 32.2 • C in 2020. The lowest value change was 3.7 • C esteeming from 17.9 • C in 1990 to 21.6 • C in 2020 (Table 2 and Figures 6, 7). From the analysis it can be concluded that only local forest coverage change is not responsible for dramatic temperature change in the study area.

DISCUSSION
The relation between land use land cover change and land surface temperature change is complex and multidirectional as land use change affects the local, regional, and global climate. Forest cover change in the Dudpukuria-Dhopachari Wildlife Sanctuary was observed throughout the study time. Decrease in forest cover before 2010 and increase after 2010 due to the declaration of DDWS as protected area by the government changed the whole scenario of the studied forest. An improvement of deforestation scenario in DDWS after declaration was also reported by Rahman and Islam (2021). Forest cover seemed as strongest driver of climate change and variability Winckler et al., 2019). The present research showed an increase in temperature by roughly 2.3-3.0 • C in 10 years in DDWS. This findings support the research of Aik et al. (2021) on Cameron Highlands, where a rise in temperature found 2.0-3.5 • C in each decade and the increase was introduced by conversion of forest cover in to urban and urban in to agriculture. Choudhury et al. (2019) also reported similar research findings, where the land surface temperature increased from 0.6 to 4.3 • C per decade. The studies by Gohain et al. (2021); Rashid et al. (2021); and Thakur et al. (2021) also found a negative relationship between LST and vegetation coverage. But, surprisingly, in the present study, during the first two decades, though the land surface temperature of DDWS increased with the decrease of forest cover, in the last decade, the land surface temperature increased, however, the forest cover increased.
Forest conversion has a very sensitive relationship with climate change. Commonly, land surface temperature decrease when forest cover increased. Hence, the present study provides a new turn to the existing thinking of the researchers. The present study concluded that local land cover change is not only responsible for climatic change or change of its variability. The impact of local land use and land cover change in the forest area is minor on the land surface temperature change; rather, it is the global warming which is suspected as the main reason behind the LST changes in this study. The results of the study coincide with the findings by Hansen et al. (2010) who researched on urban effects on global temperature change and found a small effect on global temperature. The similar kind of results observed by Parker (2010) and Hansen et al. (2001), who reported a modest effect of urban warming on near surface temperature. Global temperature change also includes the ocean data and polar region data. This means, the local level urbanization or forest cover change has a small impact on LST. Instead, globally, all type of local level land use and forest cover change along with ocean and polar region transition as a whole increase the local, regional, as well as global temperature. The global temperature then affects the LST. Global mean surface temperature (GMST) is poorly correlated with net radiation at the top of the atmosphere (Xie and Kosaka, 2017). It is estimated from the study that only sporadic increase of forest coverage will not be able to reduce local temperature but a combined timely increase of local forest coverage by each country may be able to halt the gradual increase of global temperature as well as local land surface temperature.

CONCLUSION
The present study detected the forest damages and retrieved land surface temperature change in Dudpukuria-Dhopachari Wildlife Sanctuary by applying GIS & RS for 30 years since 1990 to 2020. The forest cover change revealed a significant damages of forest from 1990 to 2010 and a dramatic forest restoration since 2010 to 2020. The land surface temperature showed the increasing pattern in the studied time period. Since 1990 to 2020, a highest temperature increases found was 7.7 • C and in each decade it was 2.3-3.0 • C. The present research revealed that local forest cover change is not the only reason behind the land surface temperature change, rather, global warming seems to have major influence for this high and rapid rise in temperature.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.