Scale effect analysis of the relationships between urban heat island and impact factors: case study in Chongqing

Abstract Several indices, including the normalized difference vegetation index (NDVI), normalized difference build-up index (NDBI), and normalized difference water index (NDWI), were retrieved from the Landsat 5 Thematic Mapper (TM) images, and the land surface temperature (LST) was reversed by the thermal radiation transfer model. Next, with the help of the thermal field variance index from LST, characteristics of spatial distribution of urban heat island were analyzed. Then, the regression models between LST and the three indices were formed by the least-squares method after they were aggregated from 30 to 120, 240, 480, 960, and 1200 m, and the relationships between the accuracies of the regression models and spatial scales were analyzed quantitatively. Finally, results from the experiment exemplified by Chongqing show that at all these spatial scales, both NDVI and NDWI are negatively correlated with LST and NDBI is positively correlated with LST. At the same spatial scale, cooling effect by NDWI acting upon LST is superior to that by NDVI, while inferior to that by NDBI to enhance LST. The fitting determination coefficients ( R 2 ) of the regression models of LST and NDVI, NDBI, and NDWI increase from 0.23 (30 m) to 0.57 (1200 m), from 0.34 (30 m) to 0.70 (1200 m), and from 0.31 (30 m) to 0.68 (1200 m), respectively. Further, the R 2 of all of the regression models and spatial scales have positive logarithmic function relations, while the standard deviation and spatial scale have negative logarithmic function relations correspondingly, and the R 2 of all of the logarithmic models are > 0.95 . The heat islands of Chongqing are roughly along northeast–southwest directions, while the heat islands of the urban core area, such as Yuzhong district, are not obvious due to the influence of the Yangtze River and the Jialing River.


Introduction
The monitoring and evaluation of the distribution and influence of urban heat island (UHI) have great significance for urban microclimate and ecological environment assessment, weather forecasting, urban planning, and so on. UHI is mainly divided into three categories: canopy layer heat island, boundary layer heat island, and surface heat island (SHI). 1 The technology of remote sensing has become one important and efficient method in monitoring and evaluating UHI, especially for SHI.
Most researches on UHI are focused on spatial distributions, driving factors, changing trends, and so on. Among the above, the relationships between driving factors and UHI had been studied deeply. By now, the analysis on the causes of UHI is mainly studied at some spatial scale.
Price and Gallo used the Advanced Very High Resolution Radiometer (AVHRR), of which the resolution is 1100 m, to investigate the relationships between surface temperature, air temperature and vegetation cover, and normalized difference vegetation index (NDVI). 2, 3 Streutker exploited a two-dimensional Gaussian model to estimate the relationships between UHI magnitude and driving factors based on AVHRR data. 4,5 Kloka et al. found that the daytime SHI intensity, namely the average surface temperature difference between urban and rural, of Rotterdam could be as large as 10°C, and the SHI is largest for neighborhoods with scarce vegetation that have a high fraction of ISA and a low albedo, with the help of AVHRR data. 6 The spatial pattern of the downscaled AVHRR land surface temperature (LST) was found to be improved when compared to that of the original AVHRR LST and to resemble more than that of Thematic Mapper (TM) LST. 7 Sua and Foodya used geographically weighted regression to test the spatial stationarity of the relationships between land cover types and the surface temperature, and found that the strength of these relationships was markedly higher than from a conventional regression analysis based on AVHRR image. 8 Schwarz et al. compared 11 different indicators for quantifying surface UHIs using MODIS data whose infrared band' resolution is 1000 m. 9 Marc et al. used ISA and LST by spatial analysis to assess the amplitude of UHI (urban-rural temperature difference), and its relationship to land development intensity, size, and eco-environment for the 38 most populated cities in the United States by MODIS data. 10 Carlson and Arthur took advantage of AVHRR and Landsat TM data, whose thermal infrared band' resolution is 120 m, to analyze the urbanization in Costa Rica region and its impact on the tropical rain forest. 11 Lo and Dale found that both temperature and NDVI had high relevance with volatile organic compounds and nitrogen oxides (NO x ) emissions but very low with the rates of cardiovascular and chronic respiratory diseases using Landsat multi spectral scanner and TM images. 12 Yuan and Bauer compared NDVI and ISA as indicators of UHI effects by investigating the relationships between LST, ISA, and NDVI by TM and Enhanced Thematic Mapper Plus (ETMþ) data. 13 Three dates of Landsat TM∕ETMþ images were utilized to document the historical morphological changes in impervious surface and vegetation coverage and to analyze the relationship between these changes and those occurred in LST. 14 Li et al. used ETMþ images to investigate how landscape composition and configuration would affect UHI in the Shanghai metropolitan region of China, based on the analysis of LST in relation to NDVI, vegetation fraction (Fv), and percent ISA. 15 Scholars like Chen, Igor, and Xiong used TM∕ETMþ images analyzed the influence to LST by NDVI, normalized difference build-up index (NDBI) and NDWI at some spatial scale. 16,17,18 The correlation between LST and NDVI, NDBI is analyzed to explore the impacts of the green land and the build-up land on the UHI with the help of TM data and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data. 19 In virtue of the frontal area index, Wong et al. had access to urban wind ventilation mapping and analyzed the relationship between the frontal area index and land use types by ASTER data and TM data. 20 Jusuf, Buyantuyev, and Callejas probed and verified that land use types had rather more influence on the increase of LST with the help of TM, ASTER, respectively. 21-23 Tiangcco used ASTER data to examine the relationship between the amount of vegetation and temperature in Metro Manila. 24 Furthermore, as far as we know, Weng et al. had initiated the study of the relationships between LST and NDVI, vegetation cover, and ISA at different spatial scales. Weng et al. used ETMþ image to study the relationships between LST and NDVI at different spatial scales, 25 but the correlations between the spatial scale and the fitting accuracy of the regression model need to be investigated further. Weng et al. used ASTER image to investigate the UHI intensity in Indianapolis at different spatial scales and the relationships between the average UHI intensity and greenness and imperviousness. 26 To conclude, comprehensive studies have been underway on the causes of UHI, as well as the relationships between UHI and land cover types, NDVI, NDBI, NDWI, and so on. However, most studies above were mainly focused on the causes of UHI at some specific spatial scale. As far as we know, currently only Weng et al. had studied the relationships between LST and NDVI at different spatial scales, 25,26 but a little pity is that they neglected to make a quantitative description on the relationships between spatial scale and accuracy of the LST-NDVI regression models. Therefore, the problems, namely the relationships between LST and impact factors at different spatial scales and the changing trends of the relationships with spatial scale, have been unresolved until now.
To solve the above problems, we retrieved NDBI, NDVI, and NDWI from Landsat 5 TM images, and reversed LST by means of the thermal radiation transfer model. Then, the regression models were constructed between LST with NDVI, NDBI, and NDWI at six different spatial scales from 30 to 1200 m, and the relationships between the spatial scale and the accuracies of the regression models were analyzed quantitatively.

Study Area
With many low mountains and hills, Chongqing is known as a "mountain city," and for its hot weather in summer, Chongqing is also nicknamed as one of the three "stove" cities in China. We chose the area within the inner ring of Chongqing, namely the core urban area and its surrounding rapid expansion area, as the area of interest. The specific location of the study area is shown in Fig. 1.
The study area is located in north latitude 29.46 to 29.65 deg and in east longitude 106.43 to 106.62 deg. The distance along south to north is about 21 km, the width along east to west is about 16 km, and its area is about 325 km 2 . The maximum elevation is 529 m, the lowest elevation is 109 m, and the average elevation is 241 m over the study area. The Yangtze River and the Jialing River embrace the Yuzhong district from south and north directions, respectively, intersecting at the ChaoTianMen pier and forming the Yuzhong peninsula. Separated by the terrains and the rivers naturally, the study area forms a multicenter structure. Being an old urban city, with a high population and building density, its heat island effects caused by rapid urbanization are very prominent.

Image Preprocessing
Landsat 5 TM images of May 23, 2010 were used in this research. Image bands through 1 to 5 and 7 have a spatial resolution of 30 m, and the thermal infrared band has a spatial resolution of 120 m. We also employed some auxiliary data such as land use map and meteorological data. Data preprocessing was mainly carried out in ENVI 4.3. First, the raw images were geo-referenced to a common universal transverse mercator projection coordinate system based on the SPOT reference images and the 1:25,000 scale topographic maps, and were resampled using the nearest neighbor algorithm with a pixel size of 30 × 30 m for all bands, including the thermal band. The root mean square error of rectification was <0.5 pixels in this study. Then, the surface reflectivities of six optical bands of the raw images were obtained after radiometric calibration and atmospheric correction by the ENVI FLAASH tool.

Derivation of NDVI, NDBI, and NDWI and Land Use Classification
The normalized indices NDVI, NDBI, and NDWI were used to characterize the land cover types and quantify the relationships between land cover type and UHI effect. NDVI is widely used because of its advantages, such as lower influence of atmospheric variations, greater sensitivity to chlorophyll, and reduction of noise by normalization. NDVI is expressed as follows: 27 (1) The NDBI index was introduced by Zha to analyze increments of reflectance on TM bands 4 and 5 for images of urbanized and barren land areas, 28 which can reflect the building information more accurately. The NDBI index is described as NDWI is a normalized difference water index, also called leaf area water-absent index, which implies the water content within the vegetation. The NDWI index is described as 29,30 Among the above three equations, ρ represents the reflectivity, DN represents the digital number. Subscript R, NIR, and MIR represents red band, near-infrared band, and middle infrared band, respectively. For Landsat TM image, R, NIR, and MIR represents band 3, band 4, and band 5, respectively.
According to the SPOT images and background information of the study area, there are eight land use types in this area, namely farmland, woodland, grassland, residential land, commercial land, traffic land, bare land, and water. From a pixel scale, the earth surface is generally regarded to be composed of three types, namely natural surface, urban, and water. 31 As the spectral features of woodland, grassland, and farmland are similar, these types belong to the natural surface; the spectral features of residential land, commercial land, traffic, and bare land are similar, these four types belong to urban.
We established the classification model in the study area using the decision tree classification tool of ENVI. The classification model is described as if ðNDVI > 0.4Þ class ¼ natural surface else if ½ðNDVI < 0.0Þ and ðNDBI < 0.0Þ class ¼ water else class ¼ urban. (4)

Derivation of Land Surface Temperature
The TM thermal infrared band (10.4 to 12.5 μm) data were utilized to derive the LST. First, the DN is converted to top-of-atmospheric (TOA) radiance measured by the instrument. Next, the TOA radiance is converted to surface-leaving radiance by removing the effects of the atmosphere. 32 where ε is the surface emissivity, τ is the atmospheric transmission, L TOA represents the TOA radiance, L upper is the upwelling or atmospheric path radiance, L down is the downwelling radiance. Three atmospheric parameters, namely τ, L upper , and L down , are inversed by the atmospheric correction tool developed by Barsi et al.,32 which is available at the NASA website (http://atmcorr.gsfc.nasa.gov).
Moreover, the ε in Eq. (5), namely the surface emissivity, was estimated according to the classification result and NDVI. The surface emissivity from satellite pixel scale is calculated by following equation: 33 where P V is the vegetation coverage, R V is the temperature ratio of vegetation, ε V is the emissivity of vegetation, R X is the temperature ratio of bare land or building, and ε X is the emissivity of bare land or building. Among them, the vegetation coverage P V is estimated as 33 In Eq. (7), NDVI V and NDVI S represent the NDVI value of vegetation and bare land, respectively.
Finally, the surface radiance is converted to LST using the Landsat specific estimate of the Planck curve, expressed as In Eq. (8), LST is the surface temperature in Kelvin (K), for Landsat 5 TM images, K1 ¼ 607.76 (Wm −2 sr −1 μm −1 ) and K2 ¼ 1260.56 K.

Relationships Between LST and Driving Factors at Different Spatial Scales
This section mainly analyzes the relationships between LST and NDVI, NDBI, NDWI, and the changing trends of these relationships at different spatial scales. The main ideas are as follows: First, the LST, NDVI, and NDBI at different spatial scales are obtained by aggregation processing, and the expressions are shown below: Xðk; lÞ; X¼ LST; NDVI; NDBI; NDWI: In Eq. (9), Xðk; lÞ represents the original pixel value, Xði; jÞ represents the aggregated pixel value, X represents LST, NDVI, NDBI, S is the aggregated area, and M is the number of pixels of aggregated area.
Then, the regression models of LST and NDVI, NDBI, NDWI are built by the least-squares method at different spatial scales, and the relationships between spatial scale and the fitting accuracy of the regression models are further analyzed.

Experimental Results
The NDVI, NDBI, NDWI, and the land use classification result of the study area are shown in Fig. 2.
Then the atmospheric upwelling radiance, atmospheric downwelling radiance, and the atmospheric transmittance of the study area of May 23, 2010, were obtained from NSAS website, which were 2.35; 3.73 Wm −2 sr −1 μm −1 , and 0.7, respectively. Also, the LST of the study area was inversed as shown in Fig. 3.
The surface temperatures of urban areas are mainly distributed in 30 to 49°C, the average temperature is 35.6°C; the surface temperature of rural areas are mainly distributed in 22 to 32°C, the average temperature is 26.8°C; and the water temperature is distributed in 17 to 23°C, the average temperature is 18.5°C.

Spatial Distribution of Surface UHIs
Heat island intensity was defined as the average temperature difference between urban and rural. The average temperature of urban area and rural area was 37.6°C and 27.9°C, respectively. So, the heat island intensity of the study area was 9.7°C, and the heat island effects were obvious.
Simultaneously, we used the thermal field variance index to analyze the heat island effects' spatial distribution of the study area quantitatively. Also, the thermal field variance index (HI) is defined as follows: 34 HI ¼ ðT s − T mean Þ∕T mean : where T s and T mean denote the surface temperature and the mean surface temperature, respectively. The thermal field variance index is a relative variable based on the surface temperature, and the higher the value is, the more obvious the heat effect is, so it can better describe UHI effect than the surface temperature. The thermal field variance index was graded into four degrees in this paper: HI ≤ 0.0 represented no heat island effect; 0.0 < HI ≤ 0.1 represented weak heat island effect; 0.1 < HI ≤ 0.3 represented stronger heat island effect; and HI > 0.3 represented strongest heat island effect. The spatial distribution of heat island of the study area was shown in Fig. 4. Fig. 2 The NDVI, NDBI, NDWI, and classification result of the study area. Fig. 3 The land surface temperature of the study area (°C).
From Fig. 4, we can find the spatial distribution patterns of UHI of the study area are as follows: (1) UHIs are distributed approximately along northeast and southwest directions.
(2) The UHI centers are located in the industrial park in the Dadukou district and the railway station in the JiangBei district, and other high energy consumed and densely populated areas, where UHI intensity is about between 5 and 10°C. (3) The heat island effects of the densely built areas such as the Yuzhong district are not obvious, mainly because of the effect of the Yangtze River and the Jialing River.

Relationships Between LST and Impact Factors at Different Spatial Scales
The surface temperature is mainly affected by the land use types and is also controlled by the various land cover parameters. The indexes NDVI, NDBI, and NDWI are the three main influencing factors. In this paper, we selected the three index factors to analyze how they affect the surface temperature and the spatial scale effects.
First, the 30-m LST, NDVI, NDBI, and NDWI images were aggregated to the resolutions of 60, 120, 240, 480, 960, and 1200 m. Due to the limited space, we only gave the aggregated results of LST, as shown in Fig. 5. With the increase of spatial scale, images became fuzzy. Especially after spatial scale was increased to 960 m, the overall contour of LST almost vanished. Figures 6-8 are the scatter plots of LST and NDVI, NDBI, NDWI at different spatial scales, respectively. Then, the regression models between LST and NDVI, NDBI, NDWI at different spatial scales were constructed. As the variation of water temperature is relatively small and the relationships between LST and NDVI, NDBI and NDWI are special, the water pixels were masked from the samples. All of the regression models are significant at 0.05 level.
It is apparent from Fig. 6 that with the increase of spatial scale, the distribution range between LST and NDVI continuously decreases. The distribution range of LST changes from 16 to 51°C (30 m) to 25 to 35°C (1200 m) while the distribution range of NDVI changes from −0.3 to 0.9 (30 m) to about 0.2 to 0.7 (1200 m). At all these spatial scales, LST tends to negatively correlate with NDVI, which proves consistent with many previous research results. [13][14][15][16][17]25 Due to NDVI, standing for normalized differential vegetation index, the more the numerical value is, the higher the vegetation coverage degree is, the lower the LST is in corresponding regions. The regression coefficient (RC) ranges from the values of   It can be found from Fig. 7 that the LST is positively correlated with NDBI at all spatial scales, which conforms to many previous research results. 16,17,18 It results from that NDBI as normalized differential building index, the bigger the NDBI value is, the higher the building coverage degree is, the higher the LST is. The RC increases from the values of 17.96 (30 m) to 31.68 (1200 m). The result shows that at the spatial scale of 30 m, as NDBI increases 0.1, LST will heighten 1.80°C; when the spatial scale rises to 1200 m, as NDBI increases 0.1, LST will heighten 3.17°C. With the increase of spatial scale, the effect by NDBI to enhance LST grows. Under the circumstance of same spatial scale, the effect by NDBI to enhance LST is greater than that by NDVI to decrease LST.  The above results clearly show that the relevance between LST and NDBI increases with the increase of spatial aggregation scale. Figure 8 throws light on the result that with the incessant increase of aggregation scale, the distribution range of NDWI is on the incessant decrease from −0.4 to 0.6 (30 m) to 0.0 to 0.5 (1200 m). Meanwhile, very like NDVI, NDWI is negatively correlated with LST at all spatial scales with its RC falling from −13.2 (30 m) to −25.8 (1200 m) at intervals. This illustrates that at the spatial scale of 30 m, as NDWI increases 0.1, LST will fall 1.32°C; at the spatial scale of 1200 m, as NDWI increases 0.1, LST will fall 2.58°C. Another discovery is that at the same spatial scale, the effect by NDWI to decrease LST is greater than that by NDVI, but is less than that by NDBI to enhance LST. Moreover, the regression model of NDWI and LST is similar with that of NDBI and LST in fitting accuracy, with R 2 increasing from 0.31 (30 m) to 0.68 (1200 m), SD decreasing from 2.34°C (30 m) to 0.94°C (1200 m). The result makes clear that the relevance between LST and NDWI is mounting up with the increase of space aggregation scale.
In order to further study the quantitative relation between fitting accuracy (including R 2 and SD) and spatial scale in the respective regression models made up of LST with NDVI, NDBI, and NDWI, we constructed regression models and two scatter charts, one of which is R 2 and spatial scale, and the other is SD and spatial scale, as is seen from Fig. 9 and the following equations: 8 < : In Eqs. (11) and (12), R 2 T−NDVI , R 2 T−NDBI , and R 2 T−NDWI represent the R 2 of the regression models of LST and NDVI, NDBI, NDWI, respectively; SD T−NDVI , SD T−NDBI , and SD T−NDWI represent the SD of the regression models of LST and NDVI, NDBI, NDWI, respectively; and x represents the spatial scale.
From Fig. 9 and Eqs. (11) and (12), R 2 in the regression models made up of LST with NDVI, NDBI, and NDWI keep increasing with the increase of spatial scale, which reveals that there is apparent positive logarithmic function relation between R 2 and spatial scale, and fitting degrees of all these three logarithmic models are relatively high, with R 2 reaching 0.975, 0.995, and 0.997, respectively. On the contrary, the SD values of the regression models of LST with NDVI, NDBI, and NDWI keep decreasing with the increase of spatial scale, which represents distinct negative logarithm function relation between SD and spatial scale, with R 2 of three logarithmic models 0.987, 0.984, and 0.951, respectively. The reason why we made such change is maybe that we adopt the method to acquire a mean value resulting in the resolution change from Fig. 9 The relationships between spatial scale and the fitting accuracies of the regression models. high to low, and the fitting accuracies of the regression models of LST with NDVI, NDBI, and NDWI continue to rise.

Conclusions
In this paper, the NDVI, NDBI, and NDWI were extracted from Landsat 5 TM images, and the LST was inversed by thermal radiation transfer model to analyze the distribution feature of UHI. Then, the regression models of the LST and NDVI, NDBI, NDWI were built at different spatial scales, and the relationships between the fitting accuracies (including R 2 and SD) of the regression models and spatial scales were discussed quantitatively. Chongqing was selected to experiment for a case study, and the results are as follows: 1. The SHIs of Chongqing are approximately along northeast and southwest directions distribution, and the heat island centers are mainly distributed nearby the high energy consumed and densely populated areas; but the heat island effects of the Yuzhong district, as the urban core area, are not obvious due to the influence of the Yangtze River and the Jialing River. 2. At all spatial scales, both NDVI and NDWI were negatively correlated with LST, and NDBI was positively correlated with LST. At the same spatial scale, cooling effect by NDWI acting upon LST is superior to that by NDVI, while inferior to that by NDBI to enhance LST. 3. The R 2 of all of the regression models and spatial scales had positive logarithmic function relations while the SD and spatial scales had negative logarithmic function relations correspondingly, and the R 2 of all of the logarithmic models are >0.95.
In future research, several works are needed to be further studied. First, the retrieval method of surface temperature needs to be improved to reduce the influence of the atmospheric impacts. Second, other aggregation methods, compared to the simple averaging aggregation, need to be further studied. Finally, because of the limitations of the complexities of UHI impact factors and the data, the conclusions of this paper just reflect the relationships between UHI and several surface indexes; the comprehensive effects of the underlying surface properties, terrains, human activities, and other impact factors on the UHI of mountainous city will be the emphasis of future works.