Evaluation of spatiotemporal dynamics of land cover and land surface temperature using spectral indices and supervised classification: A case study of Jobai Beel Area, Bangladesh

This study aims to evaluate the spatiotemporal change of land cover (LC) and surface temperature of the Jobai Beel area, an exclusive agriculture zone, situated in the far-flung area of northwest Bangladesh using satellite data. Multi-temporal Landsat series of data from 1989 to 2020 and geospatial techniques have been employed to evaluate the LC change and land surface temperature (LST) variation. Different spectral indices such as NDVI, MNDWI, NDBal have been used to retrieve individual LC. Corresponding LST has also been extracted using the thermal bands. Supervised Classification and the post-classification change detection technique were employed to determine the temporal changes and validate the individual LC. The results were employed to assess the LST variation associated with LC changes. The results reveal that the area had undergone a drastic and inconsistent heterogeneous LC transformation during the study period. Water and vegetation areas have expanded at a rate of 0.24km2/year and 0.45km2/year respectively, while bare lands have shrunk at a rate of 0.70km2/year. In general, Bare land exhibits a significant positive correlation, when Vegetation areas show a significant negative correlation with LST. However, the correlation between water areas and LST was found statistically insignificant. Agriculture in the form of vegetation has been found the most dominating land cover character throughout the study period, which has been regulating the LST variation across the area.


Introduction
*Bio-geophysical changes due to land-cover and land-use alteration are considered as one of the main controlling factors of micro-climate on local, regional, and global scales. Because of its different effects of the degree of solar radiation absorption, albedo, surface temperature, evaporation rates, the transmission of heat to the soil, storage of heat, etc., surface water bodies and vegetation have an important role in cooling the land surface temperature (LST); conversely, bare lands contribute to rise in LST (Kleerekoper et al., 2012). Recent decades have been experiencing significant heterogeneous landscape changes due to huge anthropogenic pressure. Multiple previous studies suggest that conversion of natural lands into impervious surfaces or transformation of vegetated and wetlands to agriculture lands or bare waste lands are the major causes of this change of LST.
Exploring the means of LST responses to LC change (LCc) has become a popular research topic among the researchers, which mainly concentrated in different urban landscapes of the globe (e.g. Indianapolis, USA (Liu and Weng, 2008), Melbourne, Australia (Harmay et al., 2021), Tehran in Iran (Shafizadeh-Moghadam et al., 2020), Guangzhou in China (Sun et al., 2012), Cameron Highlands in Malaysia (Aik et al., 2020), Makassar, Indonesia (Suriana et al., 2020), Lahore, Pakistan (Khan et al., 2020) West Bengal in Indian (Das and Angadi, 2020), Chattogram in Bangladesh (Roy et al., 2020), Bangkok in Thailand (Khamchiangta and Dhakal, 2020), etc. However, irrigated agricultural land is more widespread than urban land, and has significant potential for altering climate (Kueppers et al., 2007). Unlike urban areas, rural areas have been encountering different types of conversions of lands into agricultural lands and others. However, research initiatives on the relation between LCc and LST in rural areas, and especially in extremely agricultural zones at different spatial scales are relatively scarce compared to urban area studies. Hence, it is very important to evaluate the effect of LCc on LST in rural agricultural zones or areas.
Multi-temporal and multi-spectral optical remote sensing data nowadays have been used successfully and efficiently to reveal the changing characteristics of spatiotemporal dynamics of LC. In addition, the integration of thermal infrared bands with multispectral bands in optical remote sensing data has enabled to extract the LST and LC simultaneously from the same dataset (Roy et al., 2020), which helps to calculate the response of LST for LCc from corresponding temporal viewpoints. However, the limitation of these satellite data is the different dates and times of acquisitions, which could potentially affect the result of a spatiotemporal change of the study. Landsat multispectral remote sensing optical data have been used widely by scientists to assess and evaluate the spatial relation between LCc and its response by LST due to its high resolution, consistent and repetitive coverage, long historical archives, and its free availability (da Cunha et al., 2020;Rathnayake et al., 2020;Zhang et al., 2013).
Different spectral indices are very significant to reflect more and more information on water, vegetation, and bare lands. Modified Normalized Difference Water Index (MNDWI) proposed by Xu (2006) has been found very effective to delineate open water. Normally the positive value of MNDWI shows water bodies, while the negative value of MNDWI indicates soil and vegetation. Normalized Difference Vegetation Index (NDVI) is generally used to determine the vegetation abundance (Martín et al., 2020), and in LST-related studies (Julien et al., 2006). In the case of water, NDVI shows negative values; a positive value of NDVI up to 0.2 shows bare lands and built-up areas. NDVI values with more than 0.2 indicate vegetation; higher values of NDVI indicate healthy vegetation (Guha et al., 2018;Guha and Govil, 2021;Bhandari et al., 2012).
To retrieve bare land from Landsat imageries, Zhao and Chen (2005) have proposed an index namely: Normalized Difference Bareness Index (NDBaI). Value of NDBaI>0 refers to primary bare lands. Using the different threshold values, different bare lands e.g., primary, secondary bare lands could be identified (Chen et al., 2006).
Jobai Beel area, a remote rural agriculture zone, situated in the northwestern part of Bangladesh under Sapahar Upazila, now comes across huge unplanned landscape transformation due to extensive agricultural activities which ultimately lead to environmental degradation, especially in terms of climate and surface water availability. Therefore, it is very imperative to evaluate the LCc and its effect on the LST of the Jobai Beel area for future LC, and environmental planning.
This study aims to quantify and evaluate the spatiotemporal change of LC and the impact of individual LCc effects on the LST of the Jobai Beel area. Landsat series of available imageries (TM5, ETM+, and L8 OLI) have been used to reveal the characteristics of spatiotemporal dynamics of LC. Thermal bands of these corresponding imageries have been used to extract the LST of the individual land cover area. It is expected that this study would be helpful for the future environmental planning and development of sustainable rural-based agroeconomy zones like the Jobai Beel area.

Study area
Jobai Beel area is located at the transboundary point of Bangladesh and India which lies in a remote rural area of Sapahar Upazila (Sub-district) under Naogaon districts of Rajshahi division, Bangladesh (Fig. 1). The area is well-known as one of the exclusively intensified agriculture zones of the country. It consists of two beels (lake-like wetland, typically formed by the inundation of low-lying lands) namely "Domroil" and "Mahil" along with vast plain fertile lands. The area includes the parts of Shiranti, Goala, Aihai, and Pathari unions (Union consists of several villages) of the Upazila. It extends between 25 0 05ʹ15.09ʺN to 25 0 12ʹ19.22ʺN, and 88 0 27ʹ2.16ʺE to 88 0 33ʹ54.02ʺE with an average elevation of 20m (Reza and Mazumder, 2005). It enjoys the humid subtropical and dry climate. It is marked as monsoons, high temperature, considerable humidity, and moderate rainfall. The mean annual temperature of this area is around 24.88 0 C; when January is the coldest month having an average temperature of 17.5 0 C, and June with 28.4 0 C is the hottest month. The annual average rainfall is around 154.83mm; when December is the driest month having only 6mm average rainfall and June is the wettest month with an average of 369mm rainfall. The pre-monsoon month of February has been selected for this study as this is the dry and cloud-free month with most agricultural activities especially by means of irrigation. Around 100000 people live in the nearby areas. 79.03% of income comes from agriculture. The main crops are paddy, wheat, vegetable, and mustard. Horticulture like mango, jackfruit, banana, papaya, and watermelon are seen in this area. Concentrating the Beel area, fisheries, dairies, and poultries (especially duck culture) farming is very prominent here (Islam, 2003). To minimize the noise of rural human settlement and other non-agricultural activities on the image, the area has been subsetted by 13.08Km X 8.76Km size centering the wetland area to keep the maximum agriculture zone within this subset. Fig. 1 represents the location map of the study area.

Data and method
Integration of thermal infrared bands with multispectral bands in Landsat (TM, ETM+, and OLI) imageries brings about the enhanced aptitudes to derive LST and LC information simultaneously from the same image (Roy et al., 2020). Moreover, its broader swath (185km), higher spatial and temporal (30m, and 16 days respectively) resolution, and free availability of long-term historical data have made Landsat data popular among the researchers (Pain, 1985;Tang and Di, 2019).

Fig. 1: Study area map of Jobai Beel area
Therefore, considering all these through literature review, Landsat series of different sensor data (TM, ETM+, and OLI) have been used in this study. All these images have been obtained from the United States Geological Survey (USGS)'s official website https://glovis.usgs.gov. The limitation of these image data is that these images were not acquired on the same date. However, near date data of February month for each year has been used considering the minimum cloud cover and maximum agriculture activities, and almost similar climatic conditions. This different data could affect the result of the study. Hence, for better insight, online meteorological air temperature (2 meters above from surface) collected from the NASA Langley Research Centre (LaRC) "POWER" project has also been employed. The detailed specification of the employed images has been described in Table 1, and the whole methodologies have been illustrated in

Data pre-processing and processing
Selected data were reprojected to Universal Traverse Mercator Coordinate system zone 43N, using GCS_WGS_1984 datum. To improve the interpretability and quality of the selected satellite data, radiometric calibration and atmospheric correction using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction module has been performed using Envi 5.3 software to convert digital numbers (DN) to reflectance and to filter the interference from the path radiation, such as aerosol, dust particles, and water vapor.

Image classification
Selected Landsat images of six years: 1989, 2000, 2005, 2010, 2015, and 2020 have been employed to produce the land cover maps (Fig. 3). Considering its vast plain agriculture land characteristics, three broad land cover types namely: Waterbody, Vegetation, and Bare land have been selected to classify. In these classifications, the water body includes all types of surface water e.g., wetlands, rivers, streams, and canals, natural or artificial ponds, and irrigated water, etc. Vegetation includes agricultural croplands, forest and plantation patches, and other types of green vegetation. The rest of the areas from these two have been classified as Bare land which includes fallow lands, sandy bar, and a very few areas of impervious surface, and built-up areas (which is not expedient to classify separately), etc. Table 3 describes the details of selected land cover classes. Supervised classification with maximum likelihood classifier was employed using all bands (except band 6 of TM and ETM+, and band 10 and 11 of OLI/TIRS) in each case. Different 112/113 random training signatures for each class and each year were selected with the personal previous field visit experience of the author. Finally, an accuracy assessment of these classified LC maps was performed using Google Earth Pro. Error matrix table with Kappa value has been evaluated to ensure the level of satisfaction. The results reveal high acceptance of the classification for each year. Detail result has been summarized in Table 2. The results of these supervised classified images have been used to compare the corresponding land cover areas that have been extracted using different indices. Table 6 and Table 7 illustrate the comparative results of areas of distinct three LC types. Also, Table 4 shows band characteristics of used bands of Landsat different sensors.

Spectral indices
Three remote sensing-based spectral indices namely: NDVI, MNDWI, and NDBaI have been employed using Landsat imageries in an Envi 5.3 and ArcGIS 10.3 environment to characterize the threeland cover types of viz vegetation, water, and bare lands respectively to explore its relationship with LST in the study area.

NDVI
Normalized Difference Vegetation Index (NDVI) is a dimensionless index that describes the difference between visible and near-infrared reflectance of vegetation cover and can be used to estimate the density of green on an area of land (Weier and Herring, 2000). It is computed as the ratio between the TOA reflectance of a near-infrared band and a red (RED) band. The resulting values range from -1 to +1 (Guha and Govil, 2021). Vegetated areas show positive values, and densely vegetated areas are more close to value 1; when water and built-up areas show zero or negative values (Viana et al., 2019). The optimum threshold value has been selected after a precise trial and error procedure to delineate the vegetated area.

MNDWI
Substituting the middle wave infrared (MIR)/ short wave infrared (SWIR1) band for near-infrared (NIR) band, Modified Normalized Difference Water Index (MNDWI) was proposed by Xu (2006) to discriminate built-up areas from water bodies. This index is used to enhance and delineate the open water features as it can efficiently suppress and even remove the built-up land noise and as well as vegetation and soil noise. MNDWI can be expressed as follows: The resulting values range from -1 to +1; among which the positive values represent the water bodies because of their higher reflectance in the green band compared to those of the MIR band (Nath et al., 2021). A seemly threshold value has been set for the detection of water bodies after a prudent long trial and error practice.

NDBaI
To retrieve bare land from the Landsat imageries, Zhao and Chen (2005) introduced an index known as Normalized Difference Bare land Index (NDBaI). Using the band ratio between Short Wave Infrared (SWIR1) and the Thermal Infrared bands (TIR), NDBaI characterizes bare land from other LC classes. TIR can differentiate high and low levels of reflectance from built-up objects. It also reveals a high level of contrast for vegetation. The normalized difference between SWIR and TIR return near-zero value for water, negative values for vegetation, and positive values for built-up and bare land (Viana et al., 2019). NDBal can be expressed as follows:

Land surface temperature retrieval
The thermal bands of different Landsat images (band 6 of Landsat TM5 and Landsat ETM+ 7, and band 10 of Landsat OLI8) were used to extract the LST of the study area. The following steps have been performed using the software Envi 5.3 with a mathematically derived formula to retrieve the LST of the domain pixel:  Step 1: Radiometric Correction: In this step, the remotely sensed digital number (DN) is converted to spectral radiance value and data comparable. Eq. 1 is the expression of this conversion: (1) where, =Spectral radiance; = = 1; = = 255 = Eq. 1 converts DN values into spectral radiance ( ).


Step 2: Conversion of at sensor spectral radiance: At sensor, spectral radiance conversion is performed by converting the pixel values (Q) in the raw data and unprocessed image data into absolute radiance values in the line of radiometric calibration. Eq. 2 elaborates the procedure of this conversion and the process for scaling the satellite data into 8 bits (Qcalmax=255).
Step 3: Conversion of TOA reflectance ( -to-ρp): Top of Atmosphere Reflectance (TOA) of Landsat TM requires correction and processing due to its variation in solar zenith angles due to the time difference between data acquisitions. This TOA reflectance is computed following Eq. 3: where, Step 4: Sensor Radiance to effective at-sensor brightness temperature ( to T): The thermal band of TM and ETM+ (Band 6) are required to convert from at-sensor spectral radiance to effective at-sensor brightness temperature (The prelaunch calibration constant used for at-sensor temperature has been given in Table 5). The conversion formula for the at sensor's spectral radiance to at-sensor brightness temperature is as the Eq. 4.
Three spectral indices e.g., MNDWI, NDVI, and NDBaI have been employed individually to extract the three-land cover types of namely water, vegetation, and bare land, respectively Previous studies suggest that to delineate water, vegetation, and bare land exclusively MNDWI, NDVI, and NDBaI has been found competent, respectively (Chen et al., 2006;Gautam et al., 2015;Han-Qiu, 2005). After several iterations through a trial-and-error method, individual optimal thresholds were chosen for each of the indexes for each year (as the sensors are different for each year) to delineate the targeted land cover types. To verify the outputs of these indices induced LC classes, all individual results of each year's land cover types have been compared and rationalized using previously supervised classified images of corresponding classes and years. It may be mentioned that all these maximum likelihood supervised classified images have been validated earlier for each of their accuracy using Google Earth Pro (Guha and Govil, 2021). Thus, the land cover change detection over these 31 years from 1989 to 2020 has been performed. LST for the corresponding year has also been extracted using the thermal infrared band of the corresponding image. Finally, to extract the value of individual indices and their corresponding LST within each LC, Zonal Statistic tools in ArcGIS 10.3 were deployed. This LC-specific LST describes the characteristics of surface reflectance of the individual land cover responses. Spatiotemporal dynamics of the Land cover and corresponding LST have been evaluated applying different statistical analyses. The interval of the study is 5 years from 2000 to upwards. However, due to the unavailability of the archive of Landsat imageries, the interval of 1989 and 2000 could not be maintained for 5 years (though it is suggested to maintain a consistent interval for the whole study period). 1989 has been considered as the baseline year since agriculture expansion, as well as the irrigation activities, were not so high during this time. Table 6 and Table 7 are the summarized results of three land cover extents derived from supervised classifications, and three different indices, respectively. The results articulate very close similitude between the same correspondence classes in the same spatiotemporal domain. Here, vegetation has been found the main LC change driver due to the rising demand for food for the growing population. Fig. 3 depicts the land cover maps derived from supervised classification, while Figs. 4, 5, and 6 portray the water, vegetation, and bare lands areas, respectively derived from three respective indices described earlier.  Significant changes have been observed in the areal extent in LC during the time of 31 years from 1989 to 2020, though these changes were not in a consistent way throughout the time. One of the reasons for these inconsistencies could be the difference among the dates of image acquisition. Normally, December to early February is the season of paddy cultivation. As the month proceeds, the harvesting season arrives (Shelley et al., 2016). So, the difference among the LC will appear if the images are on different dates, especially if the gap is big among the dates. Another probable cause may be the diversification of agriculture. Now, due to high demand, different types of agriculture in different seasons have been introduced in this area (Assefa et al., 2021;Hoque, 2001;Mahmud et al., 1994). Therefore, plantation and harvesting time differs with the times which could bring changes in reflectance of the different date images. However, it is observed that an increase of areal coverage of surface water has occurred by 29.53% which is 7.77km 2 in the area from 1989 to 2020 (Fig. 4). The extension of agricultural irrigation and the reexcavation of wetland areas could be the cause of this change. The highest areal coverage of the water body was 41.54km 2 in 2015 which was 36.25% of the total area at that time. In 2015, the government took up a project to re-excavate the Jobai Beel wetland area, and regulate the water used for irrigation. During this time, the stakeholders of this wetland have been given the authority to monitor the water area of this wetland. All these initiatives could bring about the highest expansion of water areas. However, vegetation has been found to be the most dominant type of land cover in the study area throughout the whole 31 years of this study period. The highest area of vegetation was in 2005 which was 70.55km 2 , and 61.57% of the total area. In 2005, indiscriminate agriculture expansion was observed on the image. Even in this year, a big portion of the wetland, and bare land areas were seen to be converted into agricultural land. The overall increase of vegetation was found 35.82% which is 14.52km 2 within this 31 years' time (Fig. 5).

Fig. 3: LC maps from supervised classification
Bare land has been found as the most affected LC type and predisposed to change in the study area. It has been decreasing significantly throughout the study period from 1989 to 2020, except in 2010. The cause of the increase of bare lands in 2010 could be due to the decrease of vegetation (either it was the early period of cultivation or the lands have been kept as fallow) from the previous year. Bare lands have decreased by 22.29km 2 which is 46.70% during the study period from 1989 to 2020 (Fig. 6). The highest area of bare land was 47.73km 2 in 1989, which was 41.66% of the total area. In the initial state in 1989, due to comparatively less demand of food for less population, and lack of widespread of irrigated agriculture, a big portion of the study area were kept remaining as fallow lands or bare lands. However, in course of time, the growing population and its growing demand for food have converted more and more bare lands/fallow lands, even some parts of wetlands into agricultural land across the study area. Therefore, the result reveals that during these 31 years of this study period, water areas, and vegetation has been increasing at a rate of 0.24Km 2 /year and 0.45Km 2 /year correspondingly, when bare land has been decreasing at a rate of 0.70km 2 /year due to the emergent anthropogenic activities like irrigated agriculture and population growth. Expansion of agriculture activities reclaiming the bare lands converted to vegetated areas. At the same time as this agriculture depends on irrigation, it has increased the water area especially during the growing period of agriculture. Khan et al. (2019) has reported that in Bangladesh, the fluctuations of surface water take place due to flooded land and highly flooded irrigation lands. As the acquisition date of the imageries in this study is different, the growth of crops, amount of irrigated land, and water differ. Therefore, it is challenging to find out the exact expansion of each type of land cover as this happens arbitrarily and different agriculture starts at different times. Fig. 7 illustrates the land cover changes, and Table 8 describes the details of these change characteristics of LC.  Table 9 presents the distribution of indices' values for the study period from 1989 to 2020. It elaborates the minimum, maximum, and mean values for the individual indices. It is observed that the mean MNDWI has been changed from -0.27 to -0.14 during the study period of 31 years which indicates the increase of water bodies in the study area. The increase of mean NDVI also indicates the increase of greenness i.e., the vegetation across the study area. From 1989 to 2020 the mean NDVI increased from 0.24 to 0.37. The images that have been used in this study were acquired during the high time of irrigated agriculture. Moreover, the introduction of horticulture (mango, banana, etc.), other agriculture like corn and wheat could be the cause of this change of MNDWI and NDVI. However, the mean NDBaI has been changed from -0.18 to -0.45, which indicates the decrease of the bare lands during this time of the study. Fig. 8 portrays the LC trends of the study area. For a better understanding, it may be mentioned that for MNDWI, values>0 indicates water bodies area, and as the value proceeds to higher positive, means the higher water extent. In the case of NDVI, value <0.5 indicates the fully vegetated area, and in the case of bare land, Value > 0 is fully bare land. Therefore, all these changes of these indices support the swapping of the LC and its underlying cause. The trend of three types of land cover change has been depicted in Fig. 8.

Spatiotemporal evaluation of LST
Thermal bands of the different sensors of different Landsat images were employed to extract LST for the study area for the corresponding years of LC delineation. The LST retrieved maps for the whole study periods have presented in Fig. 9. The average aggregated different LST values for the years 1989, 2000, 2005, 2010, 2015, and 2020 have been calculated (Table 10) in 1989, 15.93 to 24.78 in 2000, 18.45 to 27.93 in 2005, 16.56 to 25.40 in 2010, 22.39 to 31.45 in 2015, and 19.43 to 26.52 0 C in 2020. The average minimum LST declined by 2.07 0 C in 31 years which is around 0.06 0 C/year. The highest average minimum LST observed in 2015 was 22.39 0 C. When the lowest average minimum LST was 15.93 0 C in 2000. The average maximum LST has also declined by 6.34 0 C within this stipulated time of 31 years at a rate of around 0.20 0 C/year. The highest LST in this group was 32.86 0 C in 1989, and the lowest LST of this group was 24.78 0 C in 2000.
The mean LST has decreased by 5.49 0 C during the last 31 years from 1989 and 2020 at a rate of around 0.17 0 C/year. This change of LST was found in the line of literature reviewed in Khan et al.'s (2019) study (Murakawa et al., 1991;Wong et al., 2012). This mean LST was highest with 27.58 0 C in 1989, and lowest was 20.44 0 C in 2010. Alike LC change, LST shows no consistent trend of change during this study period.
It is observed that the changes of all classes of temperature are too high. It may be due to the variation among the image acquisition dates. Normally, from the first week to the last week of February, the temperature varies from 5 0 C to 10 0 C. Moreover, sometimes there are cold waves, which prevail for 7 to 14 days and bring down the temperature to an abnormal scale. However, meteorological air temperature collected from the NASA Langley Research Center (LaRC) POWER Project funded through the NASA Earth Science/Applied Science Program also reveals inconsistent temperature changing pattern of air temperature at 2meter above the surface among the years of the whole study period, and big variation between the dates in line with the date of satellite image acquisition dates (28 February 1989 and 02 February 2020). It may be mentioned that satellite extracted LST provides a definite day and time's land surface temperature of a certain area when meteorological daily air temperature provides the whole day temperature for the specific day/s of the month of a certain area at a defined resolution. Therefore, these two temperatures are not the same and may vary. However, for rationalizing the condition, these two temperatures could be compared for more insight. Table 14 represents the air temperature at 2 meters above the surface of "Prediction of Worldwide Energy Resources" (Power) data of Jobai Beel area at 0.5 0 x0.5 0 latitude and longitude global grid.

Spatiotemporal evaluation of LC and LST relation
Land cover-specific extracted LST reveals that the water area exhibits the lowest land surface temperature across the area during the entire period of this study. The mean LST of water bodies were 24.00 0 C, 18.71 0 C, 20.41 0 C, 19.00 0 C, 24.36 0 C, and 21.44 0 C in 1989, 2000, 2005, 2010, 2015, and 2020 correspondingly. All the highest mean temperatures were observed on the Bare land areas throughout the study period. Bare lands showing the mean LST of 29.46 0 C, 22.20 0 C, 24.82 0 C, 22.05 0 C, 26.67 0 C, and 22.37 0 C in 1989, 2000, 2005, 2010, 2015, and 2020, respectively. In case of vegetation, the mean LSTs were 27.70 0 C, 20.54 0 C, 21.81 0 C, 20.53 0 C, 25.87 0 C, and 22.37 0 C correspondingly in 1989, 2000, 2005, 2010, 2015, and 2020. These vegetation-specific LSTs always remain in between the mean LST of water area and bare lands. However, the trends of all land cover specific LST show declining trends during these study periods (Fig. 10). The difference between the mean LST of water area and Bare lands ranges from 0.93 0 C in 2020 to 5.45 0 C in 1989. The mean difference of LST between the Water and Vegetation area ranges from 0.93 0 C in 2020 to 3.69 0 C in 1989. Mean LST difference between Bare land and vegetation area was found highest in 2005 with a reading of 3.01 0 C, though there was no difference in 2020. The lowest reading of minimum and maximum LST also predominates within the water area during the entire study period. These findings of land cover specific LST indicate that water bodies and vegetation have the cooling effect when bare lands contribute to a rise in LST, which is consistent with the previous different studies on LST and LC relation. Table 11 and Table 12 describes the details:   Table 13 illustrates the spatial distribution of LST values for the area above the mean LST and below the mean LST. It reveals that the mean LST of water areas was below the man LST of the whole area which ranges from 0.65 0 C in 2020 to 3.58 in 1989. Vegetation areas show LST below the mean LST in 2000 and 2005. However, for the rest of these two years, the mean LST of the vegetation area was mild to moderately higher than the mean LST of the whole area. The mean LST of Bare land has been always higher than the mean LST of the whole area which ranges from 0.28 0 C to 2.54 0 C in 2020 and 2005, respectively. Landscape heterogeneity and its unpredictable change would be the cause of these inconsistencies of LST among the land covers. Fig. 10 depicts the spatiotemporal distribution of LST values for the area above and below the mean LST value of the whole area. It is observed that in 1989, LST above the mean LST was dispersed almost all over the areas except the area covered by water. Water areas were remaining under the mean LST. This water area was a big wetland. Therefore, this higher temperature could be due to the abundant existence of bare lands and the absence of enough agricultural areas. Another reason for this would be due to seasonal variation of the acquisition of satellite images. Normally from the first week to the last week of February, the temperature varies from 5 0 C to 10 0 C. The image of 1989 was captured on 28 February (Table 1), which is almost the end of winter, and when the temperature normally starts to rise. In 2000, most of the areas exhibited near and below the mean LST, except the southeast part, where LST was above the mean temperature. Most of the areas show temperature below mean LST and near to mean LST in 2005. This could be due to the extensive expansion of irrigated agriculture this year. Though a big part of the southeast portion shows LST above the mean, as this area was remaining as bare lands during the time of this year 2005. A major part of the area shows LST below or near the mean LST in 2010, though LST above the mean LST prevails scattered across the area. It is observed that in 2015 and 2020, most of the areas remained below the mean LST when very few areas exhibit LST above the mean. In these two years, agriculture grabbed a significant area along with water bodies replacing the bare land areas. However, the time of acquisition of satellite image would be another reason as it was the early period of February (Table 1). Fig. 11 shows land cover-based LST trends.    Table 15 illustrates the details of correlations between LC and LST in the study area. Bare lands show a significant and positive correlation with LST which is 0.72. Vegetation also reveals a significant negative correlation with LST which is -0.65. But the Water area shows statistically a very insignificant positive correlation with LST which is 0.29. Among the Land covers, almost the same negative correlation is observed. Between the water area and vegetation area, the correlation is -0.63, and between Bare land and vegetation area, the correlation is 0.62. It reveals that vegetation has been encroaching water and bare land areas. Water area and bare land area have a very insignificant correlation. Bare land and vegetation have been found to be the main controlling factor of LST. Again, the vegetation shows the most dominant land cover type. The result also reveals that vegetation is the main driver of change which has been encroaching the bare land and water area and controlling the LST throughout the study period across the area.

Conclusion
This study was designed to investigate individual LC influence on LST especially in a rural agro-based area in a developing country like Bangladesh employing remote sensing and GIS techniques, and statistical methods. It is observed that indices viz, MNDWI, NDVI, and NDBaI can effectively quantify the land covers like water, vegetation, and bare land in that order. The study has also explored the areal change of land covers from 1989 to 2020. In general, the existing study exhibits a similar relation between LC and LST in line with the relationship that cooling source can mitigate, while hot source can aggravate the LST effect. Vegetation was found to be the most dominating land cover type. This vegetated area has been expanding at a rate of 0.45km 2 /year. As an agro-based area, this vegetation mostly indicates agriculture activities. So, agriculture has been apprehending and changing the LC pattern in this area. By the invasion process of agriculture, bare land areas have been found most affected, hence it has been declining at a rate of 0.70km 2 /year as agriculture expands. Moreover, as this agriculture is fully dependent on irrigation, therefore, apparently though it reveals the increase of water area, parament water bodies like wetlands, ponds, tanks, etc. have been shrinking. Water from these Land cover specific LST Trends (1989 Water Area LST Vegetation Area LST Bare Land Area LST permanent water bodies has also been exploited in an increasing manner for irrigation purposes and has been converted to agriculture or other land covers. Therefore, the natural setting of this area like wetlands, and parament reservoirs like ponds, tanks, etc. are under threat. Though it reveals a negative change in LST during this time of the study, several studies, however, suggest that at some other stage of the year, especially when there are no irrigation activities, the LST could be changed to increase as irrigated surface water is not available and its cooling effect is deficient at that time (Basak et al., 2013;Rahman and Lateh, 2017;Shahid et al., 2012). This study also reveals that as the increase of cooling factors like water bodies and vegetation, the LST shows downward trends, though the magnitude of these relationships between LC types and LST were not always consistent, as mixed land cover types exhibit complex results. LST also varies among each LC type along with its geographical location, seasonal variation, depth of water, vegetation density, and characteristics of bare lands. However, the study reveals the overall decrease of mean LST by 5.49 0 C from 1989 to 2020 though the image dates were not identical. The overall scenario indicates inconsistent and unplanned land use activities in the Jobai Beel area. The big natural setting of wetland in this area along with the small tanks and ponds is also perceived as vulnerable to human despoil. Therefore, proper land use planning is proximately entailed to protect the natural environment of this area. It is very challenging to explore the exact relation of different LC changes and LST. However, a more accurate result of this perceptive issue needs the same dated high-resolution satellite images with continuous field data. Meteorological variables like precipitation, potential evapotranspiration, etc. also could be employed for better insight into these relations. Therefore, it is expected that the results of this study combined with more field data could have significant implications for the policymakers and the communities to provide sustainable land use planning and land management of this area.