The Use of Land Cover Indices for Rapid Surface Urban Heat Island Detection from Multi-Temporal Landsat Imageries

The aims of this study were to determine surface urban heat island (SUHI) effects and to analyze the land use/land cover (LULC) and land surface temperature (LST) changes for 11 time periods from the years 2002 to 2020 using Landsat time series images. Bursa, which is the fourth largest metropolitan city in Turkey, was selected as the study area, and Landsat multi-temporal images of the summer season were used. Firstly, the normalized difference vegetation index (NDVI), soil-adjusted vegetation index (SAVI), modified normalized difference water index (MNDWI) and index-based built-up index (IBI) were created using the bands of Landsat images, and LULC classes were determined by applying automatic thresholding. The LST values were calculated using thermal images and SUHI effects were determined. The results show that NDVI, SAVI, MNDWI and IBI indices can be used effectively for the determination of the urban, vegetation and water LULC classes for SUHI studies, with overall classification accuracies between 89.60% and 95.90% for the used images. According to the obtained results, generally the LST values increased for almost all land cover areas between the years 2002 and 2020. The SUHI magnitudes were computed by using two methods, and it was found that there was an important increase in the 18-year time period.


Introduction
Urbanization, population growth, and industrialization lead to increases in temperature, and this situation adversely affects the world, especially in urban areas. By the year 2018, about 55% of the world's population lived in urban areas, and by 2050, this proportion is expected to increase to up to 68% [1]. Urbanization positively affects the living population with social and economic welfare but has a negative impact on the natural environment [2]. Climate change occurs on a local and global scale, particularly due to the interventions of human beings that have negative consequences for nature. The average surface temperature of the world has constantly been increasing since the late 19th century [3][4][5]. Although the temperature increases are generally seen in all types of land use/land cover (LULC) areas, the situation is even more serious in built-up areas, such as industrial, commercial, airport or dense residential areas [6,7]. Especially since the late 1980s, studies have shown that urban areas are generally warmer than rural areas around them, and this phenomenon is called the urban heat island (UHI) [8][9][10]. UHI effects occur in almost all urban areas [11], however the UHI magnitude changes according to the properties of the cities. In recent years, studies on UHI have increased considerably [12]. The UHI studies can be basically divided into two parts according to the measurement techniques [13]: (i) studies examining the effect of atmospheric UHIs by using air temperature data [14][15][16][17][18]; and (ii) studies on the effect of SUHIs using thermal remote sensing In the next step, urban, vegetation, water and other areas were extracted automatically by applying Otsu thresholding to NDVI, SAVI, MNDWI and IBI images. As mentioned previously, there are studies that have defined land cover classes using indices. However, for urban extraction, we used IBI, NDVI, MNDWI and SAVI indices together. In addition, different from the SUHI studies performed thus far, in this study, all land cover classes were generated automatically by applying Otsu thresholding. Finally, SUHI effects were determined using two techniques, and the LULC and SUHI changes were analyzed.

Materials and Methods
To investigate the SUHI effect of Bursa, Turkey, for 11 time periods from the years 2002 to 2020, multispectral and thermal bands of Landsat images were used. The flowchart of the proposed methodology is shown in Figure 1. processing was performed on the Landsat imagery. Then, NDVI, SAVI, MNDWI and IBI images were generated using the multispectral bands of the Landsat imageries, while the LST image was computed by using the thermal bands of Landsat imageries. In the next step, urban, vegetation, water and other areas were extracted automatically by applying Otsu thresholding to NDVI, SAVI, MNDWI and IBI images. As mentioned previously, there are studies that have defined land cover classes using indices. However, for urban extraction, we used IBI, NDVI, MNDWI and SAVI indices together. In addition, different from the SUHI studies performed thus far, in this study, all land cover classes were generated automatically by applying Otsu thresholding. Finally, SUHI effects were determined using two techniques, and the LULC and SUHI changes were analyzed.

Materials and Methods
To investigate the SUHI effect of Bursa, Turkey, for 11 time periods from the years 2002 to 2020, multispectral and thermal bands of Landsat images were used. The flowchart of the proposed methodology is shown in Figure 1.

Study Area
Bursa, which is located in the northwest of Turkey and southeast of the Marmara Sea, is the Turkey's fourth largest city. According to the Turkish Statistical Institute (TUIK), the 2002 and 2020 populations were 2,231,582 and 3,101,833, respectively [57]. Bursa is situated at an altitude of about 155 m and has a mild climate [58]. However, the climate also varies region by region. In the north, there is a mild and warm climate due to the Marmara Sea, and in the south, there is a harsh climate around the Uludag Mountain. The average annual precipitation is 706 mm, as of the 52-year observation period, and the average relative humidity in the province is around 69% [58]. Bursa is one of the important industrial provinces of Turkey, and there are industrial regions on three urban entrances. In addition, a natural gas combined-cycle power plant is located in the northern part of Bursa province. As a result, Bursa is developing rapidly, and its annual population growth rate was 14.85% for 2019-2020 [57]. Urbanization and industrialization are the most important triggers of the SUHI effect. Therefore, Bursa was chosen as the study area ( Figure  2).

Study Area
Bursa, which is located in the northwest of Turkey and southeast of the Marmara Sea, is the Turkey's fourth largest city. According to the Turkish Statistical Institute (TUIK), the 2002 and 2020 populations were 2,231,582 and 3,101,833, respectively [57]. Bursa is situated at an altitude of about 155 m and has a mild climate [58]. However, the climate also varies region by region. In the north, there is a mild and warm climate due to the Marmara Sea, and in the south, there is a harsh climate around the Uludag Mountain. The average annual precipitation is 706 mm, as of the 52-year observation period, and the average relative humidity in the province is around 69% [58]. Bursa is one of the important industrial provinces of Turkey, and there are industrial regions on three urban entrances. In addition, a natural gas combined-cycle power plant is located in the northern part of Bursa province. As a result, Bursa is developing rapidly, and its annual population growth rate was 14.85% for 2019-2020 [57]. Urbanization and industrialization are the most important triggers of the SUHI effect. Therefore, Bursa was chosen as the study area ( Figure 2).

Data Description and Pre-Processing
The only data used in this study were multi-temporal Landsat imageries. MODIS imageries were used for validation. The Landsat images were first obtained in 1972, and Landsat Level-1 data products are used to produce high level scientific data such as surface temperature, water area, snow-covered, and burned areas [59], and they are among the most commonly used images in remote sensing studies because they can be downloaded for free. In this study, Landsat 5 TM and Landsat 8 OLI/TIRS images were used. For the study area, the Landsat images that belonged to the summer season were investigated, and then the July images that had no clouds were selected. At the end, the images were downloaded from the USGS Earth Explorer (https://earthexplorer.usgs.gov, accessed on 1 February 2021). Landsat 5 TM images have 6 multispectral and 1 thermal bands, whereas Landsat 8 images have 8 multispectral, 1 panchromatic and 2 thermal bands. For both satellite images, the spatial resolution of the multispectral bands was 30 m. The thermal bands were 120 m and 100 m for Landsat 5 TM and Landsat 8 TIRS, respectively. However, the thermal bands were resampled to 30 m for the delivery products. The Landsat images could be rescaled to the top-of-atmosphere (TOA) reflectance, and the thermal constants (K1 and K2) needed to convert thermal band data to TOA brightness temperature are available in the metadata file [60]. Formulae for these conversions are provided in Section 2.4.
Atmospheric conditions when the images are taken may affect the quality of the LST estimate. Therefore, atmospheric correction was applied to the images. The DOS1 atmospheric correction module of QGIS open-source software was used for this process. This module is based on studies performed by Chavez [61], Sobrino et al. [62], Moran et al. [63] and Congedo [64].

Data Description and Pre-Processing
The only data used in this study were multi-temporal Landsat imageries. MODIS imageries were used for validation. The Landsat images were first obtained in 1972, and Landsat Level-1 data products are used to produce high level scientific data such as surface temperature, water area, snow-covered, and burned areas [59], and they are among the most commonly used images in remote sensing studies because they can be downloaded for free. In this study, Landsat 5 TM and Landsat 8 OLI/TIRS images were used. For the study area, the Landsat images that belonged to the summer season were investigated, and then the July images that had no clouds were selected. At the end, the Landsat 5 TM For both satellite images, the spatial resolution of the multispectral bands was 30 m. The thermal bands were 120 m and 100 m for Landsat 5 TM and Landsat 8 TIRS, respectively. However, the thermal bands were resampled to 30 m for the delivery products. The Landsat images could be rescaled to the top-of-atmosphere (TOA) reflectance, and the thermal constants (K1 and K2) needed to convert thermal band data to TOA brightness temperature are available in the metadata file [60]. Formulae for these conversions are provided in Section 2.4.
Atmospheric conditions when the images are taken may affect the quality of the LST estimate. Therefore, atmospheric correction was applied to the images. The DOS1 atmospheric correction module of QGIS open-source software was used for this process. This module is based on studies performed by Chavez [61], Sobrino et al. [62], Moran et al. [63] and Congedo [64].

Generation of Land Cover Indices
Using land cover indices together, different land cover classes can be identified [29,31,49]. In this study NDVI, SAVI, MNDWI, and IBI land cover indices were used. The NDVI is a vegetation index developed by Rouse et al. [65]. SAVI is also a vegetation index, and is used for cities which have an intense urban texture. When compared with the NDVI, the SAVI gives more accurate results for detecting less dense vegetation areas. It is useful for areas that have low vegetation to correct the NDVI for soil brightness influences [66]. The NDVI and SAVI indices are computed using Equations (1) and (2), respectively, and information about vegetation density is presented in the image. The values of these indices vary between -1 and 1. The normalized water difference index was modified (MNDWI) by Xu [67], using GREEN and MIR bands to improve open-water properties (Equation (3)). This modified version reduces built-up land noises that are frequently correlated with open water in other indices. The IBI was proposed by Xu [68] for rapidly built-up area extraction using three thematic indices, NDVI, NDBI and MNDWI. Xu stated that in the IBI image, the built-up areas are enhanced, whereas the background noise is effectively suppressed [68]. In this study, the IBI index was calculated using Equation (4), that enabled the generation of IBI immediately without having to generate the indices [68]. Figure 3 shows the obtained land cover indices for the year 2020.
(Equation (3)). This modified version reduces built-up land noises that are frequently correlated with open water in other indices. The IBI was proposed by Xu [68] for rapidly built-up area extraction using three thematic indices, NDVI, NDBI and MNDWI. Xu stated that in the IBI image, the built-up areas are enhanced, whereas the background noise is effectively suppressed [68]. In this study, the IBI index was calculated using Equation (4), that enabled the generation of IBI immediately without having to generate the indices [68]. Figure 3 shows the obtained land cover indices for the year 2020.
The formulae for land cover indices used in this study are: where GREEN is the green band, RED is the red band, NIR is the near-infrared band, MIR is the middle-infrared band, and l is the soil-adjusted correction factor (the value is between 0 and 1: 0 for very dense vegetation; 1 is very sparse vegetation). Generally, this value is accepted as 0.5.  The formulae for land cover indices used in this study are: where GREEN is the green band, RED is the red band, NIR is the near-infrared band, MIR is the middle-infrared band, and l is the soil-adjusted correction factor (the value is between ISPRS Int. J. Geo-Inf. 2021, 10, 416 6 of 20 0 and 1: 0 for very dense vegetation; 1 is very sparse vegetation). Generally, this value is accepted as 0.5.

Obtaining LST Values
Different methods have been developed to estimate LST values. Some of these methods require more than one thermal band (e.g., split window algorithm [27,69,70], multi angle algorithm [71]). However, some methods can be implemented with a single thermal band, e.g., single-channel algorithm [72,73], mono-window algorithm [74,75], and methods that are based on Planck's law [76,77]. The Landsat 5 TM has one thermal band. Uncertainties were observed in the 11 bands of the Landsat 8 images, and therefore, as the USGS has suggested, users should work with TIRS Band 10 data as a single spectral band [60,75]. Thus, in this study, a method that could be applied with a single thermal band was preferred [76]. The main reason for using this method is that it enables the calculation of LST by utilizing only emissivity values without the need for atmospheric parameters at the time of image acquisition.
Firstly, digital numbers (DN) was converted to radiance using Equation (5): where L b is the spectral radiance at the top of atmosphere in watts/m 2 *sr*µm, L min is the minimum spectral radiance for the thermal band, L max is the maximum spectral radiance for thermal band, and Q CalMax represents the maximum quantized and calibrated standard product pixel values (DN). Brightness temperature was calculated using the spectral radiance with Equation (6): where T B is the brightness temperature, and K 1 and K 2 are band-specific thermal conversion constants from the metadata; for Landsat 8 TIRS: K 1 = 774.89 and K 2 = 1321.08, and for Landsat 5 TM: K 1 = 607.76 and K 2 = 1260.56. Then, the surface temperature (T S ) was derived using Equation (7) [76,77], and it was converted to Celsius (T C ) with Equation (8).

Obtaining LST Values
Different methods have been developed to estimate LST values. Some of these methods require more than one thermal band (e.g., split window algorithm [27,69,70], multi angle algorithm [71]). However, some methods can be implemented with a single thermal band, e.g., single-channel algorithm [72,73], mono-window algorithm [74,75], and methods that are based on Planck's law [76,77]. The Landsat 5 TM has one thermal band. Uncertainties were observed in the 11 bands of the Landsat 8 images, and therefore, as the USGS has suggested, users should work with TIRS Band 10 data as a single spectral band [60,75]. Thus, in this study, a method that could be applied with a single thermal band was preferred [76]. The main reason for using this method is that it enables the calculation of LST by utilizing only emissivity values without the need for atmospheric parameters at the time of image acquisition.
Firstly, digital numbers (DN) was converted to radiance using Equation (5): where Lb is the spectral radiance at the top of atmosphere in watts/m 2 *sr*µm, Lmin is the minimum spectral radiance for the thermal band, Lmax is the maximum spectral radiance for thermal band, and QCalMax represents the maximum quantized and calibrated standard product pixel values (DN). Brightness temperature was calculated using the spectral radiance with Equation (6): where TB is the brightness temperature, and K1 and K2 are band-specific thermal conversion constants from the metadata; for Landsat 8 TIRS: K1 = 774.89 and K2 = 1321.08, and for Landsat 5 TM: K1 = 607.76 and K2 = 1260.56. Then, the surface temperature (TS) was derived using Equation (7) [76,77], and it was converted to Celsius (TC) with Equation (8).
where ε is the emissivity, ʎ is the average wavelength of the band, and α = 1.438 * 10 −2 mK.
The last step was to derive the calculation emissivity [62,78,79]. Emissivity was calculated with Equation (9): where ε is the emissivity, ε is the emissivity value of soil surfaces, ε is the emissivity value of vegetation surfaces, NDVIs is the NDVI value of soil surfaces, and NDVIv is the NDVI value of vegetation surfaces.

Obtaining LST Values
Different methods have been developed to estimate LST values. Some of these methods require more than one thermal band (e.g., split window algorithm [27,69,70], multi angle algorithm [71]). However, some methods can be implemented with a single thermal band, e.g., single-channel algorithm [72,73], mono-window algorithm [74,75], and methods that are based on Planck's law [76,77]. The Landsat 5 TM has one thermal band. Uncertainties were observed in the 11 bands of the Landsat 8 images, and therefore, as the USGS has suggested, users should work with TIRS Band 10 data as a single spectral band [60,75]. Thus, in this study, a method that could be applied with a single thermal band was preferred [76]. The main reason for using this method is that it enables the calculation of LST by utilizing only emissivity values without the need for atmospheric parameters at the time of image acquisition.
Firstly, digital numbers (DN) was converted to radiance using Equation (5): where Lb is the spectral radiance at the top of atmosphere in watts/m 2 *sr*µm, Lmin is the minimum spectral radiance for the thermal band, Lmax is the maximum spectral radiance for thermal band, and QCalMax represents the maximum quantized and calibrated standard product pixel values (DN).
where ε is the emissivity, ʎ is the average wavelength of the band, and α = 1.438 * 10 −2 mK.
The last step was to derive the calculation emissivity [62,78,79]. Emissivity was calculated with Equation (9): ε ε + (ε − ε ) * P ε NDVI < NDVI NDVI ≤ NDVI ≤ NDVI NDVI > NDVI (9) where ε is the emissivity, ε is the emissivity value of soil surfaces, ε is the emissivity value of vegetation surfaces, NDVIs is the NDVI value of soil surfaces, and NDVIv is the NDVI value of vegetation surfaces.
The last step was to derive the calculation emissivity [62,78,79]. Emissivity was calculated with Equation (9): where ε λ is the emissivity, ε sλ is the emissivity value of soil surfaces, ε vλ is the emissivity value of vegetation surfaces, NDVI s is the NDVI value of soil surfaces, and NDVI v is the NDVI value of vegetation surfaces. The P v values in Equation (9) can be calculated using Equation (10) [80]: where NDVI S = 0.2 and NDVI V = 0.5 [78]. Moreover, in Equation (9), the ε vλ is 0.99 and the ε sλ was computed with Equation (11): where RED is the red band DN values.

Extracting Land Cover Areas and Detecting SUHIs
Green plants perform photosynthesis using their chlorophyll pigments. Thus, they use radiation in the red wavelengths of the solar radiation spectrum. Therefore, they absorb a large part of the radiation in the red wavelengths. However, they reflect almost all of the radiation in the near-infrared region. As a result, NDVI shows differentiation between green areas and water surfaces, and gives different values according to the surface properties of the pixels. Hence, it is possible to mention the specific ranges of NDVI values for specific surfaces. If the surface of a pixel contains almost no vegetation in the image, the NDVI value is less than 0.2. When the surface is covered with low-density vegetation texture, the NDVI value ranges between 0.2 and 0.5, whereas when the surface is covered with high-density vegetation, the NDVI value varies between 0.5 and 1 [49,79].
In this study, the threshold values of land cover indices were determined by using the Otsu thresholding method [81]. Then, using these threshold values, urban, vegetation, water and other areas were extracted. For NDVI images, two different threshold values (lower: THR1 NDVI , and higher: THR2 NDVI ) were found by using the "multithresh" function of MATLAB R2019a programming language, which is used for determining multilevel image thresholds using Otsu's method. The vegetation areas were defined from the NDVI images using the higher threshold value (DN NDVI > THR2 NDVI ). On the other hand, for IBI, MNDWI and SAVI images, a single threshold value was computed for each image, again using Otsu's method. The IBI (DN IBI > THR IBI ) MNDWI (DN MNDWI < THR MNDWI ) , NDVI (DN NDVI < THR1 NDVI ) and SAVI (DN SAVI < THR SAVI ) images were used to determine urban areas. In addition, water areas were determined using MNDWI images (DN MNDWI > THR MNDWI ). The conditions and automatically found threshold values which were used to determine the LULC classes in this study are given in Tables 1 and 2, respectively.  After obtaining the land cover classes, morphological operations and sieve filters were used as post-processing. Dilation morphological operation was used for the detected urban areas to remove small holes, and the kernel size was determined as 3. Then, the sieve filter was applied to the whole land cover image to remove isolated class pixels.
In this study, the sizes of SUHIs for Bursa were calculated using two methods: (1) SUHI Urban-Vegetation : computing mean LST differences of urban and vegetation areas; (2) SUHI Urban-Buffer : calculating the mean LST differences of urban and surrounding buffer areas. The buffers were generated with 1 km intervals. In this method, four buffer rings were created that had 0-1 km, 1-2 km, 2-3 km, 3-4 km distances from urban areas.

Results
The results of this study were examined in three parts: (i) land cover detection; (ii) LST values and variations; and (iii) SUHI detection.

Land Cover Detection
The obtained land cover images are displayed in Figure 4. To denote the quality and usability of the produced land cover images and to assess the accuracies, test pixels were collected, and error matrices were generated ( Table 3). The determination of test pixel size is important during the classification accuracy assessment. Congalton [82] proposed using a minimum of 50 test samples for each land use category and stated that the minimum sample size should be increased to 75 or 100 samples per class if the area is large. Therefore, 500 testing pixels per class were collected in this study.
The land cover areas were successfully identified when the results were analyzed both visually and quantitatively. The overall accuracies were computed as between 89.60% and 95.90% for the 11 time periods between the years 2002 and 2020 (Table 3). These results indicate the quality and usability of the resultant images in SUHI detection studies. When the producer's and user's accuracy values were analyzed, it can be stated that the urban, vegetation and water classes which could be used for SUHI computation are quite high. As expected, the water class was determined most accurately, while the vegetation class was the second most accurately determined.   The areas of LULC classes and changes are given in Table 4. For the 18-year time period, it can be observed that the urban and vegetation areas have increased significantly. In contrast, the other areas have decreased considerably. In this time period, the water areas hardly ever changed. The urban areas were found to expand over 93.14 km 2 in 2002, whereas the extent was 210.86 km 2 in 2020. The vegetation area was computed as 383.42 km 2 and 572.44 km 2 for 2002, and 2020, respectively. On the other hand, the water area was calculated as 3.23 km 2 for the year 2002 and 3.28 km 2 for the year 2020. The increases in urban, vegetation and water areas were 117.72 km 2 , 189.02 km 2 and 0.05 km 2 , respectively, from 2002 to 2020. On the other hand, the decrease in other areas was 306.79 km 2 . These values and visual inspection indicated that, over time, the urban area has expanded into the surrounding areas.

LST Values and Variations
The produced LST images of the study area for 11 time periods from the years 2002 to 2020 are given in Figure 5; the LST graphics that show the mean LST values of urban, vegetation and water classes are given in Figure 6a; for the urban class, minimum, maximum and mean values are given in Figure 6b. When the LST values and maps were analyzed, it can be stated that in general there was an increasing trend for LST values of all classes as well as for the whole study area. However, unexpectedly, there was an extreme increase in 2003, and another unexpected slight increase in 2013.
According to the obtained results, it was observed that Bursa urban areas are warmer than vegetation and water areas in the summer season for all years. The LST values were lowest in water areas, because the daytime images were used.  The validity of the obtained LST values was assessed by using MODIS LST/Emissivity data, and the values of MODIS LST were acquired using the (0.02) rescaled factor [83].

SUHI Detection
The SUHI values were computed by using mean LST values of urban and vegetation classes and obtained to be 5.77 °C for 2002 and 7.50 °C for 2020 (Table 5). When the changes in the SUHIUrban-Vegetation results were analyzed, there was an important increase over the 18-year time period (1.73 °C). The SUHI values computed by using mean LST values of urban and buffer areas are given in Table 5. The used buffer areas are shown in Figure 7a, and the graphic that compares SUHIs computed using four buffer areas is given in Figure 7b. The obtained results indicate that, generally, there is an increasing trend in SUHI magnitudes from the years 2002 to 2020. Additionally, it can be stated that the computed SUHI magnitudes increases when the buffer ring area is further away from the urban area, and vice versa. These results can be explained by the possibility that areas farther from the city are more rural, and therefore have lower LST values. For Buffer 1 (0-1 km), Buffer 2 (1-2 km), Buffer 3 (2-

SUHI Detection
The SUHI values were computed by using mean LST values of urban and vegetation classes and obtained to be 5.77 • C for 2002 and 7.50 • C for 2020 (Table 5). When the changes in the SUHI Urban-Vegetation results were analyzed, there was an important increase over the 18-year time period (1.73 • C). The SUHI values computed by using mean LST values of urban and buffer areas are given in Table 5. The used buffer areas are shown in Figure 7a, and the graphic that compares SUHIs computed using four buffer areas is given in Figure 7b. The obtained results indicate that, generally, there is an increasing trend in SUHI magnitudes from the years 2002 to 2020. Additionally, it can be stated that the computed SUHI magnitudes increases when the buffer ring area is further away from the urban area, and vice versa. These results can be explained by the possibility that areas farther from the city are more rural, and therefore have lower LST values. For Buffer 1 (0-1 km), Buffer 2 (1-2 km), Buffer 3 (2-3 km) and Buffer 4 (3-4 km), the SUHI magnitudes were computed to be in the ranges of 0.94-3.65, 1.74-6.26, 3.07-7.90 and 3.91-9.49, respectively.  At the end of the study, instead of directly comparing images from different years, the LST was normalized and compared in order to minimize the effect of land use change and different atmospheric conditions between years. SUHI maps were created using normalized LST values using Equation (12). There are studies in the literature which have considered regions with normalized LST values greater than 1.5 °C or greater than 2 °C as SUHI areas [84,85]. SUHI = (LST − LSTmean)/std (12) where LST is the pixel LST value, LSTmean is the mean LST of the study area, and std is the standard deviation [85].
When the created SUHI maps were examined, it was seen that the normalized LST values of the urban area increased from 2002 to 2020. The SUHI sizes of some parts of the city and the impermeable surfaces near the city were higher in 2002, whereas SUHI values were higher in the urban area and its immediate surroundings after 2011. The SUHI phenomenon, which was weaker in 2002, strengthened considerably in 2020, and the thermal conditions of the city of Bursa have changed substantially with the growth of the city. It can clearly be seen that the values higher than 1.5 °C correspond to the urban heat island areas in the 2020 SUHI map (Figure 8). At the end of the study, instead of directly comparing images from different years, the LST was normalized and compared in order to minimize the effect of land use change and different atmospheric conditions between years. SUHI maps were created using normalized LST values using Equation (12). There are studies in the literature which have considered regions with normalized LST values greater than 1.5 • C or greater than 2 • C as SUHI areas [84,85]. SUHI = (LST − LSTmean)/std (12) where LST is the pixel LST value, LSTmean is the mean LST of the study area, and std is the standard deviation [85]. When the created SUHI maps were examined, it was seen that the normalized LST values of the urban area increased from 2002 to 2020. The SUHI sizes of some parts of the city and the impermeable surfaces near the city were higher in 2002, whereas SUHI values were higher in the urban area and its immediate surroundings after 2011. The SUHI phenomenon, which was weaker in 2002, strengthened considerably in 2020, and the thermal conditions of the city of Bursa have changed substantially with the growth of the city. It can clearly be seen that the values higher than 1.5 • C correspond to the urban heat island areas in the 2020 SUHI map (Figure 8). ISPRS Int. J. Geo-Inf. 2021, 10, x FOR PEER REVIEW 15 of 21 Figure 8. The UHI maps of the study area for the years: 2002,2003,2004,2007,2009,2011,2013,2015,2016, 2017 and 2020. Figure 8. The UHI maps of the study area for the years: 2002,2003,2004,2007,2009,2011,2013,2015,2016,2017 and 2020.

Discussion
When the obtained LST values for the whole study area were analyzed, it could be stated that almost all of the LST values increased from the year 2002 to 2020. On the other hand, as can be seen in the LST images ( Figure 5) and graphics (Figure 6), the 2003 LST values are extremely high. This situation may give rise to thought that the obtained results are erroneous. However, when the previous temperatures and heat waves were explored, it could be evaluated that the heat wave of 2003 in Europe may have caused these results. The 2003 summer was reported to be the hottest summer since 1500 [86][87][88]. It was said that the sweltering temperatures in summer 2003 had many negative effects, causing the extra deaths of thousands of people [88]. Therefore, according to our evaluation, it is very possible that this heat wave in Europe also affected Bursa. In addition, the 2013 LST values are also relatively high, and it could be because of a 2013 heatwave that affected Europe.
The relationships between LST and IBI and LST and NDVI images for the whole study area were analyzed, and the correlation coefficient values are presented in Table 6. The obtained correlation coefficient values indicate that LST presented a strong positive relationship with the IBI image in every situation, with correlation coefficient values higher than 0.75. Conversely, LST and NDVI images had strong negative relationships, and the correlation coefficient values were computed as between -0.69 and -0.80. Guha et al. [89] studied the relationships of NDVI, NDWI, NDBI, and normalized multiband drought index (NMDI) with LST in 11 images from the same year and four different seasons. They found that regardless of the season, LST had negative relationships with NDVI and NMDI, but positive relationships with NDWI and NDBI. These results are similar to our findings. LST generated using thermal remote sensing techniques is widely used in environmental studies. Although in this study quite high correlation values were found with the produced Landsat LST images and MODIS LST images, it should be noted that there are sources of error (such as error source due to atmospheric effects and emissivity uncertainty) when LST is retrieved from remote sensing data [90].
Although the 2003 LST values were computed to be exceptionally higher than all of the other years, for SUHI magnitudes, the situation was different. The calculated 2003 SUHI magnitudes were higher than the years 2002 and 2004, although lower than other years. This situation reveals that although sudden decreases and rises in temperatures can be observed in some years due to climatic effects, the increase in SUHI effects will continue to increase unless necessary precautions are taken, and the temperature differences between urban and rural areas are increasing.
The obtained SUHI maps show that the highest magnitudes were obtained for urban areas. In this study, residential, industrial and commercial areas were generally defined as urban areas. For detailed information and comparison of different parts of urban areas, zoning was necessary. In Figure 9, the 2020 color-infrared compositions of Landsat imagery and 2020 SUHI map are displayed. Examining these areas more closely with Google Earth images, it could be observed that the SUHI magnitudes of industrial areas had highest values (Figure 9). It was seen that especially Bursa, Demirtas, and Hasanaga organized industrial sites had the highest values. zoning was necessary. In Figure 9, the 2020 color-infrared compositions of Landsat imagery and 2020 SUHI map are displayed. Examining these areas more closely with Google Earth images, it could be observed that the SUHI magnitudes of industrial areas had highest values (Figure 9). It was seen that especially Bursa, Demirtas, and Hasanaga organized industrial sites had the highest values. Considering that SUHI is computed by differencing the mean LST values of urban and rural areas, the correct determination of urban and rural areas in these studies is a very important step that directly affects the results. The classification processes require a lot of time for this stage. The obtained results indicate that with the proposed method used in this study, urban, vegetation and water areas can be obtained quite rapidly and accurately. However, defining rural areas and determining their boundaries is more difficult than urban areas. For this reason, it can be seen that in some studies, the mean LST values of urban-vegetation, urban-water, and urban-non-urban areas are used, whereas in some other studies, urban-buffer areas around the urban area are used. There is no general consent on buffer distances to the urban area. Hence, we used four different buffer rings (0-1 km, 1-2 km, 2-3 km, 3-4 km). The obtained results indicate that the computed SUHI magnitudes were quite different from each other, and the SUHI magnitudes increased when the distance of the buffer ring to the urban area increased. Although the SUHIs were generated by using urban-buffer ring areas, it can be stated that the 3-4 km buffer to the urban area gave the most similar results to SUHIs generated by using the urban-vegetation areas. As a result, the detection of rural areas in SUHI studies is an important issue that should be addressed in future studies. For future work, we are also planning to consider seasonality, using images taken in different seasons.

Conclusions
The aims of this study were to determine SUHIs effect easily, efficiently, and rapidly, and to analyze the land cover and land surface temperature (LST) changes, solely using Landsat time series images for 11 time periods from 2002 to 2020. The proposed approach was implemented in Bursa metropolitan city, which is the fourth largest city in Turkey, using the summer season images of the years 2002,2003,2004,2007,2009,2011,2013,2015,2016,2017,2020; the LULC, LST and SUHI changes were analyzed in this 18-year time Considering that SUHI is computed by differencing the mean LST values of urban and rural areas, the correct determination of urban and rural areas in these studies is a very important step that directly affects the results. The classification processes require a lot of time for this stage. The obtained results indicate that with the proposed method used in this study, urban, vegetation and water areas can be obtained quite rapidly and accurately. However, defining rural areas and determining their boundaries is more difficult than urban areas. For this reason, it can be seen that in some studies, the mean LST values of urban-vegetation, urban-water, and urban-non-urban areas are used, whereas in some other studies, urban-buffer areas around the urban area are used. There is no general consent on buffer distances to the urban area. Hence, we used four different buffer rings (0-1 km, 1-2 km, 2-3 km, 3-4 km). The obtained results indicate that the computed SUHI magnitudes were quite different from each other, and the SUHI magnitudes increased when the distance of the buffer ring to the urban area increased. Although the SUHIs were generated by using urban-buffer ring areas, it can be stated that the 3-4 km buffer to the urban area gave the most similar results to SUHIs generated by using the urban-vegetation areas. As a result, the detection of rural areas in SUHI studies is an important issue that should be addressed in future studies. For future work, we are also planning to consider seasonality, using images taken in different seasons.

Conclusions
The aims of this study were to determine SUHIs effect easily, efficiently, and rapidly, and to analyze the land cover and land surface temperature (LST) changes, solely using Landsat time series images for 11 time periods from 2002 to 2020. The proposed approach was implemented in Bursa metropolitan city, which is the fourth largest city in Turkey, using the summer season images of the years 2002,2003,2004,2007,2009,2011,2013,2015,2016,2017,2020; the LULC, LST and SUHI changes were analyzed in this 18-year time period. For this purpose, firstly, the images for these years were obtained, pre-processed, and made suitable for analysis. Then, NDVI, SAVI, MNDWI and IBI indices were produced, and LST data were generated using only Landsat images. Subsequently, the LULC classes were extracted by applying an automatic thresholding technique. Finally, the SUHI magnitudes were calculated using: (i) urban-vegetation LST differences; and (ii) urban-buffer area LST differences.
The LULC detection results were good enough to be used in SUHI extraction purposes, with overall accuracy values between 89.60% and 95.90% for the used images. When the results of the study are examined, it can be seen that the urban area has increased significantly in the 18-year time period. Within this time period, the LST values increased for almost all land cover areas. The obtained LST values indicate that the Bursa metropolitan city is affected by urban warming due to the urbanization process, and probably also by the increase in heatwaves due to climate change. The SUHI magnitudes were computed using urban-vegetation mean LST differences and urban-buffer mean LST differences. The obtained SUHI magnitude values were computed as between 5.77 and 8.10 for the used years when urban-vegetation mean LST differences were used. On the other hand, when urban-buffer mean LST differences were used, the values were computed as between 0.94 and 10.6, and the values increased when the buffer ring area was farther away from the urban area. It was also found that there were significant changes in SUHI magnitudes from 2002 to 2020. Bursa is an industrial city and is developing very rapidly. The obtained results indicate that the urban expansion and LST increase noticeably even in just an 18-year time period, and this situation reveals the importance of considering the climate properties and taking timely precautions during urban planning processes. To mitigate the ongoing increase in the UHI effect over the years, using green/cool roofs, cool pavements, increasing the numbers of trees, and designing efficient landscape areas in urban environments could be solutions.
Funding: This research received no external funding.