Evaluating the Spectral Indices E ﬃ ciency to Quantify Daytime Surface Anthropogenic Heat Island Intensity: An Intercontinental Methodology

: The surface anthropogenic heat island (SAHI) phenomenon is one of the most important environmental concerns in urban areas. SAHIs play a signiﬁcant role in quality of urban life. Hence, the quantiﬁcation of SAHI intensity (SAHII) is of great importance. The impervious surface cover (ISC) can well reﬂect the degree and extent of anthropogenic activities in an area. Various actual ISC (AISC) datasets are available for di ﬀ erent regions of the world. However, the temporal and spatial coverage of available and accessible AISC datasets is limited. This study was aimed to evaluate the spectral indices e ﬃ ciency to daytime SAHII (DSAHII) quantiﬁcation. Consequently, 14 cities including Budapest, Bucharest, Ciechanow, Hamburg, Lyon, Madrid, Porto, and Rome in Europe and Dallas, Seattle, Minneapolis, Los Angeles, Chicago, and Phoenix in the USA, were selected. A set of 91 Landsat 8 images, the Landsat provisional surface temperature product, the High Resolution Imperviousness Layer (HRIL), and the National Land Cover Database (NLCD) imperviousness data were used as the AISC datasets for the selected cities. The spectral index-based ISC (SIISC) and land surface temperature (LST) were modelled from the Landsat 8 images. Then, a linear least square model (LLSM) obtained from the LST-AISC feature space was applied to quantify the actual SAHII of the selected cities. Finally, the SAHII of the selected cities was modelled based on the LST-SIISC feature space-derived LLSM. Finally, the values of the coe ﬃ cient of determination (R 2 ) and the root mean square error (RMSE) between the actual and modelled SAHII were calculated to evaluate and compare the performance of di ﬀ erent spectral indices in SAHII quantiﬁcation. The performance of the spectral indices used in the built LST-SIISC feature space for SAHII quantiﬁcation di ﬀ ered. The index-based built-up index (IBI) (R 2 = 0.98, RMSE = 0.34 ◦ C) and albedo (0.76, 1.39 ◦ C) performed the best and worst performance in SAHII quantiﬁcation, respectively. Our results indicate that the LST-SIISC feature space is very useful and e ﬀ ective for SAHII quantiﬁcation. The advantages of the spectral indices used in SAHII quantiﬁcation include (1) synchronization with the recording of thermal data, (2) simplicity, (3) low cost, (4) accessibility under di ﬀ erent spatial and temporal conditions, and (5) scalability.

Remote Sens. 2020, 12,2854 3 of 23 Firozjaei, et al. [48] developed a physical approach based on a triple-source surface energy balance (triple-SEB) to model LST due to AHF and SAHI intensity (SAHII). They showed that LST due to AHF in Beijing, Tehran, Istanbul, Athens, Atlanta, and Los Angeles over the past three decades ranged from 0.72, 0.58, 0.64, 0.61, 0.55, and 2.02 to 2.76, 2.32, 1.19, 1.66, 1.73, and 2.99 • C, respectively. Additionally, the SAHII value for these cities increased by 1.32, 0.95, 0.98, 0.95, 0.92, and 0.73 • C, respectively. They showed a high spatial correlation between ISC and LST due to AHF. Single date Landsat 8 images in each year were used to model LST due to AHF and SAHII variations over the past three decades.
Various studies have shown that the ISC can well reflect the degree and extent of human activity in an area. However, the accurate extraction of ISC from satellite imagery is a major challenge. Different actual ISC (AISC) datasets are available for different parts of the world. For example, the National Land Cover Database (NLCD) dataset represents surface imperviousness information for the United States of America (USA) for 1992, 2001, 2006, 2011, and 2016. The High Resolution Imperviousness Layer (HRIL) database also contains information on European impervious surfaces for 2006, 2009, 2012, and 2015. However, the temporal and spatial coverage of available and accessible AISC datasets are limited. Therefore, it is necessary to use remote sensing (RS)-based indices and methods to extract ISC information for different environmental applications.
The objective of this study was to evaluate the spectral indices efficiency to daytime SAHII (DSAHII) quantification. The innovations and distinguishing features of the present study are (1) SAHII modelling based on spectral indices and (2) evaluation of the DSAHII of some European and American cities.

Study Area
To evaluate and compare the performance of different spectral indices through SAHI modelling, 14 test sites with different conditions were selected. The test sites were Budapest, Bucharest, Ciechanow, Hamburg, Lyon, Madrid, Porto, and Rome in Europe and Minneapolis, Dallas, Phoenix, Los Angeles, Chicago, and Seattle in the USA. The geographical locations of these cities are shown in Figure 1.

Data
A set of Landsat 8 satellite image data, MODIS products, and AISC datasets were used. Details on the data used are shown in Table 2.
Landsat 8 images were used to model surface properties such as LST and various built-up indices. According to previous studies, Landsat images are suitable data for modelling and monitoring environmental conditions due to their spatial, temporal, and radiometric resolution [24,31]. The characteristics of the Landsat 8 bands are given in Table 3.   The Landsat Provisional Surface Temperature product with 30 m spatial resolution was used for USA cities. This product is generated from the Landsat Collection 1 Level-1 thermal infrared bands, Top of Atmosphere (TOA) Reflectance, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Emissivity Database (GED) data, ASTER Normalized Difference Vegetation Index (NDVI) data, and atmospheric profiles of geopotential height, specific humidity, and air temperature were extracted from reanalysis data (https://www.usgs.gov/land-resources/nli/landsat/ landsat-surface-temperature).
MOD11A1 and MOD07 products were also used to calculate and evaluate LST based on Landsat 8 images for European cities. The HRLI and NCLD datasets were used as the AISC for European and American cities, respectively.

Methods
In this study, a conceptual model with four main sections was designed ( Figure 2). First, the Landsat 8 images were preprocessed. Second, the spectral index-based ISC (SIISC) and LST were modelled based on different built-up indices (as described in Section 3.2.2), tasseled cap transformation (TCT), the biophysical composition index (BCI), and a Single-channel algorithm, from the Landsat 8 images. Additionally, a linear least squares model (LLSM) was obtained from the LST-AISC feature space was applied to quantify the actual DSAHII of the selected cities. In the third step, the DSAHII of selected cities was modelled based on the LST-SIISC feature space-derived LLSM. Finally, the value of the coefficient of determination (R2) and the root mean square error (RMSE) between the actual and modelled DSAHII were calculated to evaluate and compare the performance of the different spectral indices in DSAHII quantification.

Preprocessing
To model surface characteristics using Landsat imagery, the digital numbers of the reflective and thermal bands must be converted to top-of-atmosphere radiance and top-of-atmosphere brightness temperature (BT) based on the calibration data provided via metadata [63,64]. Then, the Fast Lineof-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model was applied to perform atmospheric correction on the Landsat reflective bands [65].

Preprocessing
To model surface characteristics using Landsat imagery, the digital numbers of the reflective and thermal bands must be converted to top-of-atmosphere radiance and top-of-atmosphere brightness temperature (BT) based on the calibration data provided via metadata [63,64]. Then, the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model was applied to perform atmospheric correction on the Landsat reflective bands [65].

Modelling LST and SIISC
Single-channel algorithm, TCT, BCI, and built-up indices were used to model LST and SIISC. The single-channel algorithm method presented by [66] was used to calculate LST. The Landsat 8 band 11 has a bias that causes an error in calculating LST [67,68]. Hence, in this method, the top-ofatmosphere BTs obtained from band 10 of Landsat 8 were used to LST calculation based on singlechannel algorithm. This algorithm can be presented as: where L is the quantity of recorded spectral radiance in the sensor for the thermal band, LSE is the amount of land surface emissivity coefficient related to the wavelength of the thermal band used, γ and δ are parameters related to the Planck function and ψ , ψ , and ψ are atmospheric functions.
MOD07 was used to calculate the amount of water vapor in the atmosphere. The NDVI threshold method was also used to calculate the pixel-scale LSE.
MOD11A1 was used to evaluate the accuracy of the LST obtained from the Landsat 8 images. First, the spatial resolution of the Landsat 8 image-derived LST was up-sampled to 1000 m. Then, the R2 and RMSE were calculated between the up-sampled LST values of the Landsat 8 images and the LST values obtained from the MOD11A1 product for each city. Additionally, the Landsat Provisional Surface Temperature product with 30 m spatial resolution was used for USA cities.
TCT is a method based on the linear combination of different spectral bands to extract information about the main surface characteristics. Equations (2)-(4) were, respectively, used to extract surface brightness, greenness, and wetness information based on Landsat 8 image bands [69].

Modelling LST and SIISC
Single-channel algorithm, TCT, BCI, and built-up indices were used to model LST and SIISC. The single-channel algorithm method presented by [66] was used to calculate LST. The Landsat 8 band 11 has a bias that causes an error in calculating LST [67,68]. Hence, in this method, the topof-atmosphere BTs obtained from band 10 of Landsat 8 were used to LST calculation based on single-channel algorithm. This algorithm can be presented as: where L sen is the quantity of recorded spectral radiance in the sensor for the thermal band, LSE is the amount of land surface emissivity coefficient related to the wavelength of the thermal band used, γ and δ are parameters related to the Planck function and ψ 1 , ψ 1 , and ψ 2 are atmospheric functions. MOD07 was used to calculate the amount of water vapor in the atmosphere. The NDVI threshold method was also used to calculate the pixel-scale LSE.
MOD11A1 was used to evaluate the accuracy of the LST obtained from the Landsat 8 images. First, the spatial resolution of the Landsat 8 image-derived LST was up-sampled to 1000 m. Then, the R2 and RMSE were calculated between the up-sampled LST values of the Landsat 8 images and the LST values obtained from the MOD11A1 product for each city. Additionally, the Landsat Provisional Surface Temperature product with 30 m spatial resolution was used for USA cities.
TCT is a method based on the linear combination of different spectral bands to extract information about the main surface characteristics. Equations (2)-(4) were, respectively, used to extract surface brightness, greenness, and wetness information based on Landsat 8 image bands [69].
where Bi indicates the surface reflectance in the i band of the Operational Land Imager (OLI) sensor. Deng and Wu [70] showed that the combination of brightness, greenness, and wetness information obtained from the TCT method based on BCI indicates the ISC. For this purpose, in the first step, the standardized brightness, greenness, and wetness maps were calculated using Equations (5)-(7).
In previous studies, various spectral indices have been developed for the extraction of built-up lands. A number of these indices were used in this study (Table 4). Information from two or more spectral bands and different spectral indices was combined to calculate these indices. Table 4. Spectral indices used in this study.

Spectral Index Equation
Equations (9) and (10) were used to calculate the normalized difference vegetation index (NDVI) and the modified normalized difference water index (MNDWI).
In this study, the mean and standard deviation (SD) values for the surface characteristics obtained from the different spectral indices were calculated for the different cities and compared with each other.

Quantifying DSAHII
Human activities such as the conversion of natural surfaces to urban surfaces are the most important factor affecting the change in SAHI. The conversion of natural surfaces into impervious urban lands increases the value of LST. ISC datasets such as HRLI and NLCD can be used to represent urban lands and human settlement regions [71,72]. In this study, HRLI and NLCD datasets were used to build the LST-AISC feature space (Figure 3). The fitted linear regression function slope, i.e., the increment of LST versus AISC, was used to quantify the DSAHII. The value of the slope indicates how much the LST value increases with increasing AISC. A higher slope value indicates a higher value of DSAHII. The process for DSAHII quantification is composed of the following four steps: (a) rescale the AISC values to between 0 and 1; (b) classify pixels based on the standardized AISC values per 100 classes with a class length of 0.01; (c) calculate the mean values of LST and rescale the AISC in each group of pixels to reduce the uncertainty caused by the heterogeneity of urban surfaces in modelling; (d) adapt an LLSM between the mean values of the LST and the rescaled AISC, in which the slope value of the fitted function indicates the value of the DSAHII. Additionally, the R2 value indicates the accuracy of the LLSM in DSAHII modelling. urban lands and human settlement regions [71,72]. In this study, HRLI and NLCD datasets were used to build the LST-AISC feature space (Figure 3). The fitted linear regression function slope, i.e., the increment of LST versus AISC, was used to quantify the DSAHII. The value of the slope indicates how much the LST value increases with increasing AISC. A higher slope value indicates a higher value of DSAHII. The process for DSAHII quantification is composed of the following four steps: (a) rescale the AISC values to between 0 and 1; (b) classify pixels based on the standardized AISC values per 100 classes with a class length of 0.01; (c) calculate the mean values of LST and rescale the AISC in each group of pixels to reduce the uncertainty caused by the heterogeneity of urban surfaces in modelling; (d) adapt an LLSM between the mean values of the LST and the rescaled AISC, in which the slope value of the fitted function indicates the value of the DSAHII. Additionally, the R2 value indicates the accuracy of the LLSM in DSAHII modelling. In this study, the DSAHII values of different cities were calculated and compared based on the LLSM obtained from the LST-AISC feature space.

Evaluating the Efficiency of SIISC for DSAHII Quantification
To evaluate and compare the performance of the SIISC in DSAHII quantification, the SIISC was used instead of AISC in the conceptual model presented in Figure 3. The SIISC parameters include UI, BI, BAEM, BU, NBBSI, SI, IBI, albedo, NDBI, brightness, ABEI, and BCI. A DSAHII value was modelled for each city based on each spectral index. To evaluate the performance of the spectral indices in DSAHII quantification, the R2 and RMSE between the modelled DSAHII based on SIISC and the actual DSAHII obtained from the AISC were calculated. In this study, the DSAHII values of different cities were calculated and compared based on the LLSM obtained from the LST-AISC feature space.

Evaluating the Efficiency of SIISC for DSAHII Quantification
To evaluate and compare the performance of the SIISC in DSAHII quantification, the SIISC was used instead of AISC in the conceptual model presented in Figure 3. The SIISC parameters include UI, BI, BAEM, BU, NBBSI, SI, IBI, albedo, NDBI, brightness, ABEI, and BCI. A DSAHII value was modelled for each city based on each spectral index. To evaluate the performance of the spectral indices in DSAHII quantification, the R2 and RMSE between the modelled DSAHII based on SIISC and the actual DSAHII obtained from the AISC were calculated.

Spatial Distribution of Spectral Index Values
The mean values of R 2 and RMSE between the LST values obtained from the Landsat 8 images and MOD11A1 for the selected cities were obtained to be 0.91 and 1.58 • C, respectively. These values indicate a reasonable accuracy of the Landsat 8-derived LST for these cities [2,73]. The spectral index values of selected cities were spatially heterogeneous (Figures 4 and 5). The values of built-up land indices, BCI-derived characteristics, and LST in the central areas of the cities were higher than those in the suburbs.
The mean values of R 2 and RMSE between the LST values obtained from the Landsat 8 images and MOD11A1 for the selected cities were obtained to be 0.91 and 1.58 °C, respectively. These values indicate a reasonable accuracy of the Landsat 8-derived LST for these cities [2,73]. The spectral index values of selected cities were spatially heterogeneous (Figures 4 and 5). The values of built-up land indices, BCI-derived characteristics, and LST in the central areas of the cities were higher than those in the suburbs.    cities, respectively. In general, the spatial variation of AISC was higher in USA cities than in European cities. Among the various indices, AISC and ABEI had the highest and lowest CV, respectively.

Quantifying DSAHII
The LST-AISC feature space formed for the different cities is shown in Figure 6.

Quantifying DSAHII
The LST-AISC feature space formed for the different cities is shown in Figure 6.      The mean DSAHII values for selected European cities were 4.5, 6.6, 4.3, 3.5, 3.0, 3.0, 2.3, and 1.9 • C, respectively, and those for USA cities were 3.0, 6.7, 5.6, 4.4, 5.4, and 1.9 • C, respectively. The impact of human activity on LST varied among the selected cities. Among those, Rome and Seattle had the highest and lowest negative impacts of human activities on LST, respectively. The mean DSAHII values for the selected cities in Europe and USA were 3.6 and 4.5 • C, respectively. Generally, in green cities (with large fraction of vegetation coverage) including Hamburg, Budapest, Porto, Bucharest, Minneapolis, Seattle, Chicago, and Dallas due to high surface wetness and vegetation cover, and low heat and dryness, DSAHII is more intense. While in desert cities including Ciechanow, Madrid, Lyon, Rome, Los Angeles, and Phoenix, the DSAHII is lower.

Evaluating the Effectiveness of SIISC for DSAHII Quantification
The performance of the SIISC parameters in DSAHII quantification differed (Figure 7).

Discussion
SAHIs are one of the important negative effects of human activity in the natural environment [37,48]. Increasing human activity increases the percentage of impermeable surfaces and increases the LST of these areas compared to that in natural areas (Figures 4 and 5).
Marando, et al. [74] investigated the effect of green infrastructure elements such as urban and peri-urban forests, street trees, as well as the effect of vegetation cover and tree diversity in the reduction of the SUHI effect in Rome, Italy. The results of this study show that the green infrastructure significantly reduces the SUHI phenomenon in a Mediterranean city. Grigoras , and Urit , escu [75] conducted an analysis based on multi-time remote sensing data to investigate the impact of land use change in Bucharest's SUHI. The results suggest that the increase in built-up lands and the decrease in vegetation cover due to anthropogenic activities caused an increase in surface temperature and expansion of the area affected by SUHI. Arnds, et al. [76] analyzed the spatio-temporal variance of the SUHI of Hamburg. In summary, the SUHI showed a radial gradient in the center, which is mostly corresponding to the urban densities. Dian, et al. [77] studied the relationship between SUHII and local climate zones (LCZ) classes for Budapest. The results of this investigation indicate that as the density of the building decreases, the intensity of SUHI also decreases. The highest SUHII is in the city center and the lowest intensity of SUHI with negative values can be found in vegetation-covered LCZ classes.
Due to the negative consequences of the SAHI effect on various aspects of human quality of life, its quantification is of great importance. Firozjaei, Weng, Zhao, Kiavarz, Lu and Alavipanah [48] used a triple-SEB to model SAHII. The results showed that the triple-SEB could be highly effective for SAHII modelling. However, triple-SEBs are highly complex and require many calculations. Additionally, the implementation of this model requires many input datasets, including land cover parameters, surface digital models, climatic conditions, and so on.
Various studies have shown that ISC information is a good index for the degree of urban-related human activity in an area [38,41,43,44]. Zhang and Cheng [72] and Li, Zhou, Li, Meng, Wang, Wu and Sodoudi [71] used the LST-ISC feature space to model SUHII. The most important challenge of this method is using appropriate ISC information. Existing ISC databases have serious drawbacks, including spatial and temporal coverage constraints. However, satellite imagery can be used to address these challenges. In previous studies, various spectral indices and methods have been proposed for ISC modelling and built-up land extraction [51,53,61].
The results of this study showed that the IBI, BU, and NBBSI indices show good performance in DSAHII modelling (Figure 7). The TCT-derived brightness did not perform well in DSAHII modelling. Combining brightness with greenness and wetness information in the BCI increases the accuracy of DSAHII modelling. Some studies have shown that BCI can be effective in demonstrating spatial changes in the ISC in urban environments [24,70]. Firozjaei, Sedighi, Kiavarz, Qureshi, Haase and Alavipanah [61] showed that the ABEI is more effective than other indices for separating built-up lands from other land covers, especially bare lands. However, in this study, the ABEI accuracy for DSAHII quantification was lower than that of other indices. Therefore, this study showed that the ABEI is not suitable for heterogeneous modelling within built-up lands.
In general, SIISC has advantages for quantifying DSAHII, such as concurrency with thermal data recording, simplicity, low cost, accessibility under different spatial and temporal conditions, and scalability. The results showed that the use of the LST-SIISC feature space was highly effective for DSAHII modelling. However, one of the limitations of this method is that it is unable to model DSAHII changes in different geographical locations within a city. Therefore, to increase the spatial resolution of the modelled DSAHII, the LST-SIISC feature space must be implemented locally, such as for different urban regions.

Conclusions
SAHI modelling and quantification are very important to the quality of urban life. In this study, to evaluate and compare spectral indices used for DSAHII modelling, 14 cities in Europe and the USA with different conditions were selected. The DSAHII was quantified using the LST-AISC feature space and the LST-SIISC feature space. The results showed that the DSAHII in the selected cities in Europe and the USA was different. The DSAHII in cities with humid climates was higher than that in cities with dry climates. The performance of the spectral indices in DSAHII quantification varied. The results showed that IBI had the best performance for DSAHII quantification. In general, regarding the advantages of SIISC, it can be useful in identifying and characterizing the effects of human activity on the urban environment. It is suggested that in future studies, based on the approach presented in this study, the DSAHII in cities worldwide should be examined multi-temporally. Providing an appropriate model for future DSAHII prediction is also an important area for future studies.