Impact of spatiotemporal land-use and land-cover changes on surface urban heat islands in a semiarid region using Landsat data

ABSTRACT

Many factors are involved in urban heat island development, such as lack of green spaces, improper choice of building materials, densification, and other human activities. The aim of this research was to quantify the effects of land-use/land-cover (LU/LC) changes on urban land surface temperature (LST) during a 25-year period (1993–2018) for the semiarid Shiraz City in southern Iran using Landsat data (TM, ETM⁺, and OLI/TIRS) and machine learning algorithms. Five main LU/LC classes, such as orchard, vegetation, bare surface, asphalt cover, and built-up areas, were identified using a support vector machine algorithm. Landsat images were used to retrieve normalized difference vegetation index (NDVI) and normalized difference built-up index (NDBI). The results showed that the mean LST over the entire study domain increased considerably between 1993 and 2018, due to urbanization, decrease of green areas, and increasing industrial areas. Built-up areas increased considerably by 25.8% from 80 to 100.6 km² between 1993 and 2018, while vegetation cover decreased dramatically by 69.3%. Mean LST increased from 38.4 to 40.2°C during the 25-year period with a significant increase of 3.9°C between 2013 and 2018. In addition, the Urban heat island Ratio Index (URI) displayed a substantial upward trend during the 25-year period.

KEYWORDS:

• Urbanization
• urban heat island
• land use
• land cover
• semiarid
• population density

1. Introduction

The global urban population represented 54.5% in 2016 (Rousta et al. 2018). It is expected that 2/3 of the global population will be concentrated on urban areas by the mid-twentieth century (Jaber 2018; Hosseini, Pourahmad, and Pajoohan 2016). This means that 2/3 of the global population will not only experience increased temperature from global warming but also increased temperature due to Urban Heat Island (UHI) effects. In this regard, land-use and land-cover (LU/LC) changes are two main factors for urban planning, environmental monitoring, and land management (Chaudhuri and Mishra 2016).

One of the main environmental consequences of urbanization is the UHI effect, a phenomenon that significantly increases the temperature in urban areas as compared to surrounding rural areas (Ranagalage, Estoque, Handayani et al. 2018; Wang, Derdouri, and Murayama 2018). Some of the negative impacts of the UHI include increased energy consumption, air pollutant emissions, and greenhouse gases causing negative environmental effects for human health and comfortability (Simwanda et al. 2019).

The effects of UHI are probably more significant for arid and semiarid regions. These regions are already experiencing extreme summer temperatures combined with the absence of rainfall that could clean the air and reduce the temperature. Iran is a typical country in this respect with fast urbanization. The largest city in southern Iran is Shiraz with an increasing population and densification of the urban area as well as the growing industry. The growing urban areas and UHI effects are deemed to lead to negative consequences for the inhabitants.

In general, UHI can be divided into atmospheric and surface UHIs. The atmospheric UHI is calculated based on air temperature while the surface UHI is calculated based on land surface temperature (LST) (Mathew, Khandelwal, and Kaul 2016; Ranagalage, Estoque, and Murayama 2017; Ranagalage, Estoque, Handayani et al. 2018). The main contributor to surface UHI in metropolitan cities is land surface warming due to heat exchange effects from building materials involved in the urbanization (Hansen 2010). The surface UHIs are as well affected by the geographical location and local weather conditions as influenced by changes in the thermal properties of the urban infrastructure (Karakuş 2019).

A key element in challenges associated with surface UHI is the land-use/land-cover (LU/LC) changes over a period of time (Chaudhuri and Mishra 2016). The LU/LC change is the interface between human activities and surrounding environment that can impact the urban climate (Silva, da Silva, and Santos 2018). Observations of LU/LC change at various temporal and spatial scales are important for predicting changes in the urban environment (Gashu and Gebre-Egziabher 2018). Generally, LU/LC change is a function of urbanization defining changes in surface use (Amiri et al. 2009). For instance, orchard and vegetation cover can be replaced by impervious surfaces and buildings. This may lead to increased surface temperature (Steeneveld et al. 2014). Rise in surface temperature has the potential to expose inhabitants to heat-related stress (Guha et al. 2018; Guha 2017; Mushore et al. 2017).

To quantify surface UHI and LU/LC changes, the use of geographical information systems (GIS) and remote sensing (RS) techniques is becoming increasingly crucial. These techniques are economical to monitor magnitude and, spatial and temporal variation of UHI (Karakuş 2019). The most common indices to estimate the spatial and temporal variation of LST are the normalized difference vegetation index (NDVI) and the normalized difference built-up index (NDBI), which are derived from remote sensing data (Kumar and Shekhar 2015). Both of these are used as indicators of LU/LC changes (Yuan and Bauer 2007).

Many studies have investigated relationships between LU/LC and LST using RS imagery and GIS at regional and global scales (Bokaie et al. 2016; Chen et al. 2017; Dissanayake et al. 2019; Sultana and Satyanarayana 2020). However, very few studies have focused on arid and semiarid urban regions where extreme surface UHIs may be expected. Shiraz city, which is one of the biggest cities of Iran, is located in the arid and semiarid area and under the pressure of rapid urbanization and population increase. As Shiraz has developed rapidly, impermeable surfaces have continuously replaced the green area. Impermeable surfaces for living, industrial, and commercial uses, produce a high amount of heat especially in the newly developed districts in the city center. To efficiently reduce UHI in Shiraz, it is necessary to investigate the LU/LC changes and their impacts on LST. Hence, the main objective of the current study was to investigate the relationships between LU/LC changes with surface UHI using RS and GIS techniques. To enhance our analysis, the Gram-Schmidt method was used to increase the resolution of Landsat data (from 30*30 cell size to 15*15) that, in turn, enhances the accuracy of the LU/LC classification. Mean LST of each LU/LC class was extracted to assess the distributed LST changes. The urban heat island ratio index (URI) was used to compare the UHIs intensity for the study period (1993–2018). The analysis provides a framework to investigate the impact of urbanization on surface UHI. The results are useful for urban planners and decision makers to manage the negative temperature effects of urbanization in semiarid regions.

2. Materials and methods

2.1. Study area

Shiraz is the fifth most populous city in Iran and the capital of Fars Province with over 1.5 million inhabitants. It is located in southern Iran with an area of 145 km² between longitude 52° 25′ 45″ and 52° 37′ 60″ E and latitude 29° 32′ 16″ and 29° 42′ 42″ N (Figure 1). It is known as one of the Iranian cities with the fastest-growing urban population (Sabet Sarvestani, Ibrahim, and Kanaroglou 2011). The average surface elevation of the city is about 1500 m above mean sea level. Shiraz is divided into eight districts with different management. The city includes large industrial areas as well as tourist and historical attractions. The annual mean temperature is 23.3°C and mean annual precipitation is about 300 mm. The climate characteristic of the study is shown in Figure 2 (1993–2018). The figure shows that the dry period spans from April to August and most of the rain (255.8 mm; 81%) falls from October to February.

Figure 1. Location of Shiraz City in Fars Province, Iran.

Figure 2. Ambrothermic graph of Shiraz City for the 25-year study period.

2.2. Data availability and pre-processing

To estimate the surface UHI influenced by urbanization, four satellite images of Landsat TM/ETM+ (thermal band 6) and Landsat 8 TIRS (thermal band 10) for 1993, 2003, 2013, and 2018 (May–June period) were acquired from the United States Geological Survey (USGS) (Table 1). We focused on daytime images since these are the only available data. The time steps were selected to investigate the LU/LC changes and their impact on UHI over a 25-year period. The QUick Atmospheric Correction module (QUAC) was used to calculate atmospheric correction using the ENVI 5.3 software as a pre-processing step. Multispectral band data for all images were converted into surface reflectance (0–1), while thermal band data were converted into at-sensor brightness temperature, expressed in degrees Celsius.

Table 1. Details of the Landsat images used in this study.

Sensor	Path/raw	Acquisition date	Time (GMT)	Thermal conversion constant		Sun-elevation	Sun-azimuth
				K1	K2
Landsat 5 TM	163/39	29/05/1993	6:26:03	607.56	1260.56	61.037	100.059
Landsat 7 ETM^+	162/40	26/05/2003	6:46:25	666.09	1282.71	66.500	103.760
Landsat 8 OLI/TIRS	162/40	29/05/2013	6:59:47	774.88	1321.07	69.520	105.790
Landsat 8 OLI/TIRS	163/39	03/06/2018	7:02:17	774.88	1321.07	68.860	106.970

2.3. Land-use/land-cover classification and accuracy assessment

To increase the resolution of the multispectral bands in the Landsat images (ETM+/OLI) from 30 to 15 m resolution, the Gram-Schmidt method was used for pan-sharpening in ENVI 5.3. This was done to increase the lower spatial resolution of the multispectral Landsat imagery based on a higher resolution of panchromatic imagery. To detect and classify the LU/LC changes during the study periods (1993, 2003, 2013, and 2018), multispectral bands of TM/ETM+/OLI images and a support vector machine algorithm (supervised classification) were used (Dissanayake et al. 2019). Shiraz City was categorized into five major LU/LC classes encompassing orchard area, vegetation, bare surface, asphalt cover, and built-up area (Table 2). Orchard area is defined as tree cover such as pomegranate, pine, etc. while vegetation-covered area is defined as farming land and grass-covered areas. The accuracy of LU/LC classification was evaluated by selecting more than 300 checkpoints (random sampling) on the ground for each LU/LC class compared with the corresponding grid in the satellite images using the overall accuracy, kappa coefficient, user accuracy, and producer accuracy indices. These indices are commonly used for quantitative assessment of the accuracy of LU/LC classification (Sultana and Satyanarayana 2020; Dissanayake et al. 2019).

Table 2. Definition of land-use/land-cover (LU/LC) classes in Shiraz.

Class of LU/LC	Description LU/LC class
Orchard area	Tree cover such as pomegranate, pine, etc.
Vegetation	Farming land, grass covered areas
Bare land	Bare soil surface without any cover
Asphalt cover	Asphalt cover such as streets, highways, and parking lots
Built-up area	Buildings such as residential, commercial, and industrial

2.4. Retrieval of land surface temperature

To estimate LST, raw thermal band data (band 6 in Landsat 5 (TM) and Landsat 7 (ETM⁺), and band 10 in Landsat 8 (TIRS)) were converted into spectral radiance using (Shirani-bidabadi et al. 2019; Karakuş 2019):(1) $L_{\varphi }=M_{\mathrm{L}}\cdot Q_{\mathrm{C}\mathrm{a}\mathrm{l}}+A_{\mathrm{L}}$(1) where L_ϕ is the top of atmospheric spectral radiance W/(m²·sr·μm), M_L is the band-specific multiplicative rescaling factor from the metadata, Q_Cal is quantized and calibrated standard product pixel values (DN), and A_L is the band-specific additive rescaling factor from the metadata.

In the next step, the at-sensor (brightness) temperature was evaluated by Bokaie et al. (2016, 2019), El-Hattab, Amany, and Lamia (2018) and Nurwanda and Honjo (2020):(2) $T_{\mathrm{B}}=\left[\frac{K_2}{\mathrm{L}\mathrm{n}\left({\displaystyle \frac{K_1}{L_{\varphi }}+1}\right)}\right]$(2) where T_B is the brightness temperature, and K₁ and K₂ are the thermal conversion constants from the metadata presented in Table 2.

To calculate the LST, we firstly evaluated the land surface emissivity (ε) according to Dissanayake et al. (2019) and Simwanda et al. (2019):(3) $\epsilon =m{P}_{\mathrm{V}}+n$(3) where ε is the land surface emissivity, n = 0.004, m = 0.986, and P_V is the proportion of vegetation (Sultana and Satyanarayana 2020; Shirani-bidabadi et al. 2019):(4) $P_{\mathrm{V}}={\left(\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}+\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}\right)}^2$(4)

NDVI is one of the major indicators to quantify vegetation cover. The range of NDVI is between −1 and +1, where small positive values indicate built-up or bare land, large positive values indicate vegetation (Yuan and Bauer 2007), and negative values denote non-vegetated surfaces and water (Silleos et al. 2006). NDVI was calculated from the red and near-infrared bands:(5) $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\left(\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}\right)$(5) where NIR refers to surface reflectance of band 4 (Landsat 5 TM) and band 5 (Landsat 8), and Red refers to surface reflectance of band 3 (Landsat 5 TM) and band 4 (Landsat 8).

Finally, LST in °C was estimated by Sultana and Satyanarayana (2020):(6) $\mathrm{L}\mathrm{S}\mathrm{T}=\left(\frac{T_{\mathrm{B}}}{1+\left(\lambda \times {\displaystyle \frac{T_{\mathrm{B}}}{\rho }}\right)\ln \epsilon }\right)-273.15$(6) where λ is the wave length of emitted radiance (λ = 11.5 μm for Landsat TM and ETM⁺ band 6, and λ = 10.8 μm for Landsat OLI band 10), ρ = h × c/σ (1.438 × 10⁻² m K), h = Planck’s constant (6.626 × 10⁻³⁴ Js), c = velocity of light (2.998 × 108 m/s), and σ = Boltzmann constant (1.38 × 10⁻²³ J/K).

2.5. Retrieval of normalized difference built-up index (NDBI)

The NDBI is an index for quantifying built-up areas and impervious surfaces, which were derived from Landsat satellite images using reflectance in the NIR and mid-infrared (MIR) bands (Faqe Ibrahim 2017; Ranagalage, Estoque, and Murayama 2017). The NDBI was calculated using:(7) $\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{I}=\left(\frac{\mathrm{M}\mathrm{I}\mathrm{R}-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{M}\mathrm{I}\mathrm{R}+\mathrm{N}\mathrm{I}\mathrm{R}}\right)$(7) where MIR = band 5 (for Landsat 5 TM and 7 ETM⁺) and band 6 (for Landsat 8 OLI), and NIR = band 4 (for Landsat 5 TM and7 ETM⁺) and band 5 (for Landsat 8 OLI).

2.6. Change detection of LST and land-use/land-cover

The primary goal of this study was to determine changes of LU/LC and their impacts on the UHI. To retrieve LU/LC and LST changes, the phase shift of a given pixel over the time was compared using cross-tabulation through gain and loss analysis (Pal and Ziaul 2017).

2.7. Retrieval of urban heat island ratio index

Spatiotemporal variation of UHI was evaluated by calculating LST during the 25-year study period. For this, the LST data at different time steps (1993, 2003, 2013, and 2018) were normalized (Shirani-bidabadi et al. 2019; Xu and Chen 2004):(8) $\mathrm{N}\mathrm{L}\mathrm{S}{\mathrm{T}}_i=\frac{\left(\mathrm{L}\mathrm{S}{\mathrm{T}}_i-\mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}}\right)}{\left(\mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}}-\mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}}\right)}$(8) where NLST_i is the normalized pixel value, LST_i is the LST value of pixel i, and LST_minimum and LST_maximum are the minimum and maximum LST value of each time step, respectively. The normalized LSTs were then classified using mean standard deviation into five thermodynamic levels including low, sub-medium, medium, sub-high, and high (Table 3).

Table 3. Normalized LST classification using mean standard deviation.

Interval	Classification
T ≤ Tmean − 1.5 std	Low
Tmean − 1.5 std < T < Tmean − std	Sub-medium
Tmean − std < T ≤ Tmean + std	Medium
Tmean + std < T ≤ Tmean + 1.5 std	Sub-high
T > Tmean + 1.5 std	High

Thus, the distribution of LST classes for Shiraz City was established and the area within each class was determined according to Table 3. Finally, the URI was calculated to assess the contribution of urbanization to UHI by Qiao et al. (2014), Shirani-bidabadi et al. (2019) and Xu and Chen (2004):(9) $\mathrm{U}\mathrm{R}\mathrm{I}=\frac{1}{100m}\sum _{i=1}^n{w}_i{p}_i$(9) where m is the number of normalized LST classes, i is the class of temperature higher than medium, n is the number of LST levels higher than the medium class, w_i is the weighted value using the sub-high and high classes, and p is the per cent of i.

3. Results and discussion

This section focuses on the effect of spatiotemporal land-use and land-cover changes on UHI following the accuracy assessment of LU/LC classification. To achieve this aim, the results are given in four main sub-sections including (і) classification of LU/LC and their spatiotemporal changes; (іі) relationship between LU/LC changes and LST; (ііі) relationship between LST and, NDVI, NDBI, and population; and (іv) calculation of UHI’s intensity over time using URI.

3.1. Accuracy assessment of LU/LC classification

As mentioned above, Shiraz City was classified into five different LU/LC categories including orchard, vegetation, bare land, asphalt surface, and built-up areas at four different time periods using the support vector machine algorithm (1993, 2003, 2013, and 2018). The results of the accuracy assessment are presented in Table 4. The table shows that user and producer accuracy indices yielded more than about 90% modeling accuracy for all time periods. However, the overall accuracy index of the 4-year classification using the support vector machine algorithm was about 98%. The kappa coefficient for the four-year classification was 97, 98, 97, and 97, respectively. According to previous research, kappa coefficients greater than 0.75 indicate that classification and reference data are compatible (Nigatu, Dick, and Tveite 2014). According to Table 4, the overall accuracy of LU/LC classification is well preserved.

Table 4. Accuracy assessment of LU/LC classification for the studied period.

Indices	LU/LC	1993	2003	2013	2018
User accuracy (%)	Orchard	89	96	95	98
	Vegetation	99	93	90	96
	Bare land	95	99	99	99
	Asphalt	97	99	99	99
	Built-up	95	99	99	97
Producer accuracy (%)	Orchard	82	98	97	99
	Vegetation	99	86	84	94
	Bare land	85	100	100	99
	Asphalt	99	99	99	99
	Built-up	97	99	99	98
Overall accuracy (%)		98	98	98	98
Kappa coefficient (%)		97	98	97	98

3.2. LU/LC classification and change in patterns between 1993 and 2018

Figure 3 shows the spatial distribution of LU/LC classification for the study area in 1993, 2003, 2013, and 2018, respectively, produced by the support vector machine algorithm. As mentioned before, the entire domain was categorized into five typical land uses. The calculated area and percentage of each land-use class are presented in Table 5, and a summary of LU/LC changes for the studied period is shown in Table 6. In Table 6, a minus (−) sign indicates that a particular land-use class decreased and a plus (+) sign indicates an increase compared to the previous time step. The results show that the built-up areas increased considerably by 24% from 80 km² to 99.4 km² between 1993 and 2003, while a more gradual increase was noted from 2003 to 2018 (Tables 5 and 6). However, the total area of bare land decreased dramatically by 64% from 35.9 to 12.9 km² between 1993 and 2018 (Table 5), with a main change in the first time period by 44% reduction (Table 6). Moreover, the majority of bare land was converted into built-up areas and asphalt cover by 57% and 10% during the whole study period, respectively. Furthermore, the total orchard area declined slightly by 13.7% during the entire 25-year period (Table 6 and Figure 3). The vegetation cover decreased three-fold from 1993 to 2018, with a rapid change between 2013 and 2018 by nearly 54% reduction (Table 6). However, the results also show that 47% of the vegetation cover was converted into built-up areas from 1993 to 2018. Furthermore, Table 5 shows that the asphalt cover steadily increased during the period of study. The total area of asphalt cover increased dramatically from 6.7 km² in 1993 to 17.3 km² in 2018. Generally, the outcomes of LU/LC change detection reveal that the spatial extent of the urban area of the city increased from 1993 to 2018. However, impervious surfaces such as built-up and asphalt areas expanded to the east, south, and west parts of the city. Green areas, such as orchards and vegetation cover, and bare land were converted into impervious surfaces. Hence, rapid urbanization had detrimental effects on the temperature-stabilizing areas of the city. Therefore, similar results are reported by the researchers in previous studies (Balew and Korme 2020; Fatemi and Narangifard 2019; El-Hattab, Amany, and Lamia 2018; Rousta et al. 2018).

Figure 3. Land-Use/Land-Cover classification pattern in 1993, 2003, 2013, and 2018.

Table 5. Land-use/land-cover areas in 1993, 2003, 2013, and 2018.

LU/LC types	1993		2003		2013		2018
	km^2	%	km^2	%	km^2	%	km^2	%
Orchard	13.6	9.3	9.1	6.2	11.8	8.1	11.7	8
Vegetation	9.2	6.3	8.2	5.7	6.1	4.2	2.8	1.9
Bare land	35.9	24.7	20.0	13.7	15.5	10.6	12.9	8.9
Asphalt cover	6.7	4.6	8.6	5.9	11.7	8.0	17.3	11.9
Built-up areas	80.0	55.0	99.4	68.4	100.2	69.0	100.6	69.2

Table 6. Land-use/land-cover change statistics of Shiraz City in 1993, 2003, 2013, and 2018.

LU/LC types	1993–2003		2003–2013		2013–2018		1993–2013		1993–2018
	km^2	%	km^2	%	km^2	%	km^2	%	km^2	%
Orchard	−4.5	−33.2	2.7	30.2	−0.1	−0.9	−1.8	−13.0	−1.9	−13.7
Vegetation	−1.0	−10.6	−2.1	−25.3	−3.3	−54.1	−3.1	−33.2	−6.4	−69.3
Bare land	−15.9	−44.3	−4.5	−22.7	−2.6	−16.5	−20.4	−56.9	−23.0	−64.0
Asphalt cover	1.9	27.8	3.1	35.6	5.6	48.0	5.0	74.9	10.6	158.3
Built-up areas	19.5	24.3	0.8	0.8	0.4	0.4	20.2	25.3	20.6	25.8

3.3. Relationship between LST and LU/LC

Much research has been devoted on quantifying general relationships between LU/LC changes and surface thermal signature. The LST distribution of Shiraz for the studied period is depicted in Figure 4, and the statistical measures of LST are presented in Table 7. These findings reveal that the LST changed from 26 to 51°C, 16–57°C, 27–51°C, and 31–55°C during 1993, 2003, 2013, and 2018, respectively. Besides, it can be seen that the mean LST increased from 38.4 to 40.2°C between 1993 and 2018, while in the first decade it decreased by 1.8°C (Table 7). Furthermore, the mean LST did not change significantly from 2003 to 2013. However, a significant change with 3.9°C increase occurred between 2013 and 2018. Mean LST for each LU/LC class is summarized in Table 8. The results show that the highest and lowest LST are related to the bare land and orchard classes, respectively. The surface temperature for bare land increased by 3.5°C between 1993 and 2018. Moreover, the maximum change in LST was recorded for the orchard class by an increase of 4.8°C. The LST for orchard and vegetation cover in the first decade decreased by 3.4°C and 2.5°C, respectively, while it increased dramatically by 8.2 and 7.0°C from 2003 to 2018, respectively. The LST for built-up areas increased by 4.2°C from 2013 to 2018 but a decrease by 3.4°C was observed between 1993 and 2013. In addition, surface temperature for asphalt cover increased as well especially between 2013 and 2018 by 3.6°C. In principal, the results show that the largest increase occurred during the last five years from 2013 to 2018.

Figure 4. Distribution of LST over Shiraz City in 1993, 2003, 2013, and 2018.

Table 7. Summary of LST (°C) for Shiraz City during May–June (1993–2018).

Year	Minimum (°C)	Maximum (°C)	Mean (°C)	Standard deviation
1993	26.5	51.5	38.4	2.8
2003	16.9	57.7	36.6	3.3
2013	27.3	51.5	36.3	3.0
2018	31.3	55.3	40.2	2.6

Table 8. Mean LST for each LU/LC class for Shiraz City (1993–2018).

LU/LC type	Mean LST (°C)				Mean difference (°C)
	1993	2003	2013	2018	1993–2003	2003–2013	2013–2018	1993–2013	1993–2018
Orchard	32.2	28.8	33.7	37.0	−3.4	4.9	3.3	1.5	4.8
Vegetation	35.1	32.6	36.0	39.5	−2.5	3.5	3.5	0.9	4.5
Bare land	39.4	39.4	40.3	42.9	−0.0	0.9	2.6	0.9	3.5
Asphalt cover	38.2	37.4	37.6	41.2	−0.9	0.3	3.6	−0.6	3.0
Built-up areas	39.3	37.0	35.9	40.1	−2.3	−1.1	4.2	−3.4	0.8

Generally, the results reveal that the largest increase in LST is related to buildings and urban infrastructure development. The lowest increase in LST is related to orchard and vegetation cover. This is primarily due to high evapotranspiration from vegetation that decreases temperature (Simwanda et al. 2019; Priyankara et al. 2019; Balew and Korme 2020). Research has shown that green cover tends to reduce the temperature by absorption and reflection of radiation and adjustment of latent and sensible heat exchange (Yue et al. 2007; Lin et al. 2017). As well, impervious surfaces and human-made constructions play a striking role in increasing the temperature (concrete, asphalt, and metal) (Frazer 2005; Nakayama and Fujita 2010). The LST spatiotemporal variation is associated with urban expansion and the LU/LC variation, and this affects the UHI of the city (Silva, da Silva, and Santos 2018; Correia Filho et al. 2019; Balew and Korme 2020).

3.4. Relationship between NDVI and LST

The NDVI has been extensively used in research to measure vegetation cover. The index has also been used in many UHI studies as the cooling effect of vegetation related to latent heat vaporization and the convective cooling of tree canopies facilitate the air cooling process (Chun and Guldmann 2014). Accordingly, the distribution of NDVI for Shiraz City during the 25-year period is shown in Figure 5, and the related statistical summary is given in Table 9. An NDVI between 0.1 and 0.75, in general, denotes vegetation cover, while NDVI more than 0.75 indicates a dense canopy cover. The NDVI for bare land and soil is close to zero, and negative values depict water surface. The calculated mean NDVI was 0.2, 0.17, 0.18, and 0.13 in 1993, 2003, 2013, and 2018, respectively. The lowest values were observed in areas that are composed of bare land and buildings. Hence, the results illustrate that most of the city is covered by impervious surfaces and buildings. Typically, it is difficult to distinguish between the NDVI of bare land relative to buildings (Varshney 2013; Ma et al. 2019; Khanal et al. 2019). On the other hand, the area covered by vegetation will have a lower surface temperature (LST) and higher NDVI as shown in Figures 4 and 5. Moreover, the highest NDVI was found in the northwest and southeast parts of the city, which is covered by orchard and vegetation, respectively. The results of the NDVI analysis show that this index, in general, decreased during the 25-year period that is a consequence of the LU/LC changes from green to built-up areas (e.g. Karakuş 2019; Bokaie et al. 2016; Sultana and Satyanarayana 2020).

Figure 5. Distribution of NDVI over Shiraz City in 1993, 2003, 2013, and 2018.

Table 9. Summary statistics for NDVI for Shiraz during the 25-year study period.

Year	Minimum	Maximum	Mean	Standard deviation
1993	−1.00	1.00	0.20	0.22
2003	−1.00	1.00	0.17	0.18
2013	−0.45	0.92	0.18	0.18
2018	−0.60	0.88	0.13	0.16

The correlation between LST and NDVI is shown in Figure 6. The figure shows that there is a negative correlation between LST and NDVI for all time steps with R² = 0.84, 0.66, 0.07, and 0.52 in 1993, 2003, 2013, and 2018, respectively. Earlier studies for cities in different countries indicate comparable. Figure 7 depicts the average air temperature anomaly in May for the entire study period. According to the figure, the low R² between LST and NDVI in 2013 could be due to an unusually low temperature in this year (May) as compared to the 25-year average. However, due to the inverse correlation, it can be stated that decreasing vegetation cover results in an increased surface temperature (Alemu 2019; Ning, Gao, and Chen 2017; Deng et al. 2018; Balew and Korme 2020).

Figure 6. Correlation between LST and NDVI in 1993, 2003, 2013, and 2018.

Figure 7. Temperature anomaly for the 25-year period.

3.5. Relationship between NDBI and LST

Spatiotemporal variation of NDBI for the four time steps (1993, 2003, 2013, and 2018) is shown in Figure 8, and a statistical summary is given in Table 10. One of the main indicators for understanding the urban climate is NDBI (Ranagalage, Estoque, and Murayama 2017; Rousta et al. 2018; Grover and Singh 2016). NDBI ranges between −1 and +1, where small positive values denote bare land, large positive values built-up areas, and negative values denote vegetation cover (Ranagalage, Estoque, Zhang et al. 2018). Ranagalage, Estoque, Handayani et al. (2018) suggested that an NDBI larger than −0.22 represents built-up areas (Ranagalage, Estoque, Handayani et al. 2018). Based on our calculations, mean NDBI decreased from −0.22 in 1993 to −0.26 in 2003, while it increased to −0.12 in 2018. The significant increase in NDBI is due to the reduction of bare land and lower vegetation, as well as increase in the built-up areas. The similar results are reported by the researchers working on cities of different parts of the world (Guha 2017; Chen and Zhang 2017; Correia Filho et al. 2019; Balew and Korme 2020).

Figure 8. Distribution of NDBI over Shiraz in 1993, 2003, 2013, and 2018.

Table 10. Summary statistics of NDBI for Shiraz during the 25-year study period.

Year	Minimum	Maximum	Mean	Standard deviation
1993	−1.00	1.00	−0.22	0.21
2003	−1.00	0.34	−0.26	0.16
2013	−1.00	0.43	−0.21	0.16
2018	−1.00	0.58	−0.12	0.16

The correlation between LST and NDBI is shown in Figure 9. According to the figure, there is a positive correlation between LST and NDVI for all time steps, where R² was 0.82, 0.62, 0.05, and 0.44 in 1993, 2003, 2013, and 2018, respectively (for similar results see e.g. Balew and Korme 2020; Ranagalage, Estoque, and Murayama 2017; Faqe Ibrahim 2017; Rousta et al. 2018). The insignificant correlation between LST and NDBI in 2013 could be due to considerable lower temperature in that year relative to the long-term average. According to the positive correlation, it can be expected that the increase in built-up areas results in increased land surface temperature (Simwanda et al. 2019; Dissanayake et al. 2019; Wang, Derdouri, and Murayama 2018; Rousta et al. 2018; Ranagalage, Estoque, Zhang et al. 2018). It is noted that there was no positive correlation between NDVI and NDVI indices during the whole study period.

Figure 9. Correlation between LST and NDBI in 1993, 2003, 2013, and 2018.

3.6. Relationship between population increase and LST

The population of Shiraz City has increased significantly during the four last decades (Sabet Sarvestani, Ibrahim, and Kanaroglou 2011). This growth has expanded built-up areas on the expense of bare land and green areas (Li, Zhang, and Kainz 2012; Sultana and Satyanarayana 2020; Zhou et al. 2014; Bokaie et al. 2016; Patra et al. 2018). However, since the population density for the entire study period was inaccessible, we used the population density of each district for Shiraz in 2016, as shown in Figure 10. The figure shows that the density of population in district 7 is higher than other districts as it is the commercial center of the city. However, due to the occurrence of orchards and vegetation cover in this district, it had a lower LST than other districts for all time steps (Figure 11). As can be seen from Figure 10, district 3 is low-populated while the mean LST was the same as for districts 4, 5, and 6, which are located in the central parts of Shiraz. The LST increased in these districts due to extensive building activities combined with the absence of green areas (Bokaie et al. 2019). Population increase has caused the growth of residential structures, at the same time, more energy is trapped in building structures because of high heat capacity. This leads to the need for cooling systems that again increase the anthropogenic heat.

Figure 10. Spatial distribution of population density in Shiraz in 2016.

Figure 11. Mean LST for different districts of Shiraz in 1993, 2003, 2013, and 2018.

3.7. Urban heat island ratio index

The results of normalized LST classification are shown in Figure 12. The figure illustrates that the normalized LST classes higher than medium were mostly located in urban areas. This illustrates that the higher classes expanded in the urban areas through the entire study period, especially, in the southern part of Shiraz City, where the LU/LC changed from bare land to built-up areas. Furthermore, the result of urban heat island ratio index (URI) calculations is shown in Figure 13. The URI graph shows that the intensity of UHI increased continuously during the entire period of study. The intensity rose dramatically from 0.045 in 1993–0.17 in 2018. This can be compared to the URI calculated by (Shirani-bidabadi et al. 2019) for Isfahan City in central Iran during a 17-year period. Although the trend of URI in both studies is the same, the magnitude differs. Both studies indicate that the decision makers simply ignored enhancing green spaces and other nature-based cooling measures for these two major cities in Iran.

Figure 12. Normalized LST classification for 1993, 2003, 2013, and 2018.

Figure 13. Urban heat island Ratio Index (URI) for 1993, 2003, 2013, and 2018.

4. Conclusions

This study analyzed the effects of LU/LC change on UHI based on LST during a 25-year period (1993–2018) for Shiraz City using Landsat data (TM, ETM+, and OLI/TIRS). The spatiotemporal pattern of LU/LC changes was quantified for four different time steps and five different LU/LC classes including orchard, vegetation, bare surface, asphalt cover, and built-up areas. To achieve the objectives, the spatiotemporal pattern of LST, NDVI, and NDBI, variation of LST following LULU types, and variation of LST within the NDVI and NDBI were examined. Moreover, URI was used to depict the variation of the intensity of UHI in the study area.

The results of LU/LC change illustrate that the built-up areas increased substantially from 80 km² to 100.6 km² between 1993 and 2018. The total bare land area decreased dramatically by 64% between 1993 and 2018. Thus, the major part of bare land was converted to built-up areas. Besides, the vegetation cover decreased considerably and 47% of the vegetation cover was converted to built-up areas.

According to the LST analysis, the bare land and orchards were identified as the hottest and the coolest places for all time steps, respectively. The mean LST in 2003 was 36.6°C that increased to 40.2°C in 2018. The results showed that the LST pattern changed with the LU/LC. The maximum change in LST was recorded for the orchard cover by an increase of 4.8°C between 1993 and 2018. The LST for vegetation cover increased by 4.5°C from 1993 to 2018. In addition, the LST for built-up areas increased by 4.2°C from 2013 to 2018. In principal, the results illustrate that the largest increase occurred during the last five years from 2013 to 2018.

Additionally, the results showed a significant negative correlation between mean LST and NDVI, while a positive correlation was found between mean LST and NDBI for all time steps. The results of normalized LST classification revealed that higher classes expanded in the urban areas where the LU/LC changed from bare land to built-up areas. The URI results showed that the intensity of UHI increased dramatically from 0.045 in 1993 to 0.17 in 2018. In addition, areas of high temperature increased for all time steps due to a considerable growth in population and subsequent urbanization.

In general, the results of this study indicate that different types of LU/LC have a substantial impact on temperature. Urbanization played a key role in creating an intensified UHI during the 25-year study period. The emergence and expansion of this phenomenon depend on many factors such as human activities, air pollution, LU/LC types, and weather conditions. For this, measures are needed to reduce their adverse effects on the urban environment. The most efficient measure to reduce these negative effects is green space expansion in the urban and industrial context. Promoting green roof technology with the aim of maximizing the absorption of radiation and establishing urban gardens are other ways to reverse the negative trend of UHIs.

Acknowledgements

The authors are grateful to the United States Geological Survey (USGS) team for providing the Landsat images via https://earthexplorer.usgs.gov/. We also would like to thank the Centre for Middle Eastern Studies at Lund University and Ministry of Science, Research and Technology of the Islamic Republic of Iran. The authors also thank three anonymous reviewers and the journal editor for constructive comments on the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).