Examining the Relationship Between Spatial Configurations of Urban Impervious Surfaces and Land Surface Temperature

The urban heat island (UHI) effect has significant effects on the quality of life and public health. Numerous studies have addressed the relationship between UHI and the increase in urban impervious surface area (ISA), but few of them have considered the impact of the spatial configuration of ISA on UHI. Land surface temperature (LST) may be affected not only by urban land cover, but also by neighboring land cover. The aim of this research was to investigate the effects of the abundance and spatial association of ISAs on LST. Taking Harbin City, China as an example, the impact of ISA spatial association on LST measurements was examined. The abundance of ISAs and the LST measurements were derived from Landsat Thematic Mapper (TM) imagery of 2000 and 2010, and the spatial association patterns of ISAs were calculated using the local Moran’s I index. The impacts of ISA abundance and spatial association on LST were examined using correlation analysis. The results suggested that LST has significant positive associations with both ISA abundance and the Moran’s I index of ISAs, indicating that both the abundance and spatial clustering of ISAs contribute to elevated values of LST. It was also found that LST is positively associated with clustering of high-ISA-percentage areas (i.e.,>50%) and negatively associated with clustering of low-ISA-percentage areas (i.e.,<25%). The results suggest that, in addition to the abundance of ISAs, their spatial association has a significant effect on UHIs.


Introduction
Owing to rapid urbanization and economic growth in China during recent years, a large amount of natural land has been transformed into developed land (e.g., for residential, industrial, or transportation purposes). One of the direct environmental consequences of this transformation is modification of underlying land cover composition, with a significant increase in urban impervious surface area (ISA). ISA represents land surface area that is correlated with anthropogenic features, such as asphalt, concrete, metal, or glass. There has been a corresponding significant increase in urban land surface temperature (LST) (Shao and Xiao, 2018). An analysis of urban warming characteristics associated with impervious surfaces has been reported for Wuhan in China (Fan et al, 2019), with the results showing that the increase in ISA in this region has led to elevations of the LST. Likewise, Weng and Lu (2008) examined the interplay between ISA fraction, vegetation fraction, and LST at the subpixel level. They concluded that LST increased concurrently with increasing ISA. Tang and Xu's study (2013) showed a strong positive correlation between LST and ISA in six cities. Furthermore, they observed that a higher percentage of ISA corresponded to a greater acceleration in the increase of LST. It has also been found that different ISA materials, such as low-albedo and high-albedo materials have different thermal characteristics (Wu and Murray, 2003;Deng and Wu, 2013). Klok et al. (2012) explored the urban surface heat island (UHI) of Rotterdam and found that the LST was higher within clusters of low-albedo urban surface characteristics compared with other areas. Models based on spectral unmixing thermal mixing (SUTM) have provided quantitative estimates of LST through consideration of the different thermal properties of dark and bright impervious surfaces (Deng and Wu, 2013). In addition to ISA abundance, the spatial configurations of ISAs may also have an impact on LST. Landscape configuration and structure, spatial arrangement of ISAs, the normalized difference vegetation index (NDVI), and population distribution are all related to urban warming (Yang et al., 2016). Feng et al. (2017) also found that UHI was affected by a combination of various driving factors and spatial scales. Zheng et al. (2014) assessed the effect of the spatial configuration of anthropogenic land cover on the urban warming environment, and found that large patches of paved surface (>50%) caused a significant UHI.
Although the aforementioned studies were performed to examine both the relationship between LST and ISA abundance and the relationship between LST and ISA spatial configuration, the latter relationship remains unclear. That is, the ways in which the patterns of clustering or dispersal of areas with low, medium, and high percentage of ISA affect LST are yet to be adequately addressed. Therefore, the objective of this study is to examine the impact of spatial configurations (especially spatial association) of ISAs on LST. Taking Harbin, China as an example, we examined the impact of ISA spatial association on LST measurements. In particular, ISA abundance and LST measurements were derived from Landsat Thematic Mapper (TM) imagery from the years 2000 and 2010. Further, the spatial association patterns of ISAs were calculated using the local Moran's I index. Finally, the impacts of ISA abundance and spatial association on LST were examined using correlation analysis.

Study area and data
Harbin is the capital of Heilongjiang Province, China, and is located by the Songhua River in Northeast China. It has a semi-moist continental monsoon climate, with an annual mean temperature of 3.4℃ and an annual precipitation of 500 mm, of which 85% falls between June and September (Huang and Wu, 2009;Wang et al, 2009). The study area included the central city region and some regions located in southern Harbin, specifically Pingfang, Nangang, Daoli, Xiangfang, Songbei, Daowai, and Hulan. These regions represent diverse categories of urban land use and land cover types, including commercial, industrial, and residential buildings, transportation, agricultural, water, and other open lands, such as wetland, barren land, and grassland.
Two cloud-free Landsat TM images were acquired on 26th September, 2000 and 22nd September, 2010. These were reprojected on a Universal Transverse Mercator (UTM) projection with WGS84 datum and UTM zone 52. The digital numbers (DNs) of multispectral bands were converted to at-satellite reflectance according to the Landsat Calibration model from ENVI 5.1 software. In addition, thermal band DNs were converted to brightness temperature to retrieve LST. The original thermal band of the TM images was resampled to a 30-m resolution for consistency with the spatial resolution of the other bands. An image acquired in 2010 from Google Earth was obtained for impervious surface validation.

Land surface temperature retrieval and standardized
To obtain accurate LST values, we adopted the monowindow algorithm (MWA) developed by Qin et al. (2001). The LST (T s ) can be expressed as where ε is the emissivity for the Landsat TM thermal band 6, τ is the atmospheric transmittance on the image acquisition date, and T is the brightness temperature for band 6. T a is the effective mean atmospheric temperature. The LST (℃) is given by 273.15 We adopted the NDVI thresholds method used by Sobrino et al. (2001). This method considers different impacts of distinct land cover types (e.g., water, vegetation, bare soil, and impervious surface) to determine emissivity. Prior to the calculation using this method, an iterative self-organizing data analysis (ISODATA) unsupervised classification was performed to obtain water pixels, which were then assigned an emissivity of 0.990 (Snyder et al., 1998;Okwen et al., 2011). Fully vegetated pixels were assigned an emissivity of 0.985. Similarly, non-water pixels were considered as non-vegetated land cover and assigned an emissivity of 0.972 (Snyder et al.,1998;Sobrino et al., 2004). Finally, mixed pixels included various degrees of vegetation and non-vegetation (e.g., bare soil and impervious surfaces). The emissivity of these mixed pixels was calculated as follows (Sobrino et al., 1990;Carlson and Ripley, 1997;Sobrino and Raissouni, 2000): where NDVI s and NDVI v are the NDVI values for non-vegetation and full vegetation, respectively, p v is the scaled NDVI, v  and  n are the emissivities of vegetation and non-vegetation, respectively, and F = 0.55, is a shape factor taking account of the geometrical distribution (Roberts et al., 1998). For the atmospheric transmittance, we used information from the atmospheric correction model (http://atmcorr. gsfc.nasa.gov/) on image acquisition date, latitude, and longitude. Using the Landsat Calibration model from the ENVI 5.1 software, thermal band DNs were converted to brightness temperature T to determine the effective mean atmospheric temperature T a . As the Landsat image was acquired in September and Harbin is located in a mid-latitude region, we adopted the regression model developed by Qin et al. (2001) for the summer season of mid-latitude regions: where T 0 is the near-surface air temperature (K), which was acquired from China field station climate data sets (http://www.gscloud.cn/). The intensity of urban heat island was obtained by standardizing the surface temperature retrieved from the images. min max min where UHII is urban heat island intensity, T d is the surface temperature, T max is the maximum surface temperature, T min is the minimum surface temperature.

Estimation of ISA
To estimate the abundance of ISA, we applied a spectral unmixing model (Roberts et al., 1998;Wu and Murray, 2003;Wu, 2004). The linear spectral mixture analysis (SMA) model assumes that the reflectance of each pixel is a linear weighted sum of the reflectances of various homogeneous land cover types (called 'endmembers') (Smith et al., 1990;Adamset al., 1995;Lu and Weng, 2004). The specific model of linear SMA with full constraints can be formulated as follows: where R i is the reflectance for each band i (1,…,7) of the remote sensing image, N is the total number of endmembers, f k is the resultant fraction of endmember k (1,2,…, N), R ki is the reflectance of endmember k in band i, and e i is the model residual.
The tasseled cap (TC) transformation was originally developed by Kauth and Thomas (1976) to identify agricultural crop developments. In our study, we attempted to obtain each endmember between TC components (Crist, 1985;Huang et al., 2002). Following the approach of Ridd's conceptual V-I-S triangle model, four spectral endmembers, namely, vegetation, impervious surface, soil, and water, were first selected from the feature space scatter plot of the TC component analysis. The first three TC components were all normalized to the range from 0 to 1. Feature space scatter plots were then generated using these first three TC components. The Landsat TM images of the study area were used to carry out the TC transformation. Following this, each derived TC component was linearly normalized within the range from 0 to 1 (Deng and Wu, 2012). This methodology can be expressed as min max min where TC i (i=1, 2, 3) are the first three TC components, and TC imin and TC imax are the minimum and maximum values of these components. Samples of each endmember were selected and plotted in three feature space images.

Spatial pattern of ISA
The spatial pattern of ISA was characterized using the local Moran's I (Anselin, 1995), which has been found to be an effective measure of spatial patterns of land cover . A binary map was generated for ISA, with ISA being assigned a value of 1 and other land cover classes a value of 0. The local Moran's I index (I i (d)) was calculated according to where Z i is the value of the attribute of interest at location i, Z mean is the mean of the attribute z values, w ij (d) is a spatial weight matrix, and d represents the size of the neighborhood. If a pixel is located within a distance of d from the center pixel i, it is considered in the calculation. In this context, d is equal to 30 m, with the assumption that all pixels observed within a particular spatial extent are spatially dependent. The local Moran's I values of impervious surfaces were normalized in the range of −1 to 1. Local Moran's I values of −1 represent highly dispersed patterns, values of 0 random patterns, and values of 1 highly clustered patterns. Fig. 1 shows hypothetical spatial patterns of impervious surface (value 1 in Fig. 1) and their corresponding local Moran's I values.

Accuracy assessment
To assess the accuracy of the calculated ISA, 200 randomly distributed pixels were selected. A sampling unit of 3 × 3 pixels was employed to reduce the effect of misregistration of Landsat TM and Google Earth images. Reference abundances of ISAs were visually digitized on the Google Earth aerial photographs. Three measures of precision were applied to quantify the relative estimation error at the pixel level: root mean square error RMSE, system error SE, and mean absolute error MAE. These metrics can be expressed as follows: where x ei is the estimated ISA for pixel i, x ri is the actual ISA of pixel i from the Google Earth image, and N is the total number of pixels of the samples.

Intensity and magnitude of UHI in 2000-2010
The absolute LST differed from year to year, which made it difficult to compare LSTs for different periods. Therefore, to explore the variance of the UHI in Harbin, all LST values had to be normalized to the range from 0 to 1. One potential problem was that the results might be affected significantly by abnormal values (e.g., for cloud, metal, or water). To reduce the impact of extreme values on the normalized LST results, we visually identified pixels with abnormal values and removed these pixels from the calculation. During the 10 yr from 2000 to 2010, the mean UHI intensity within the study area increased by 18.54% (Table 1 and Fig. 2). Pingfang District had the smallest increase in UHI intensity (6.77%), whereas the central city region (the traffic ring) had the largest increase (24.30%). Hulan, Daoli, and Daowai had similar increases in UHI intensity. By 2010, UHI areas in Harbin had clearly expanded, with a magnitude of UHI greater than 0.7 in the larger areas.

Intensity and magnitude of ISAs in 2000-2010
Through applying the spectral unmixing method described in Section 2.2.2, we derived the percent ISA maps for 2000 and 2010 (Fig. 3). To evaluate the performance of spectral unmixing for mapping ISA abundance and distribution, a reference ISA was derived by manual interpretation of 200 randomly selected samples from the Google Earth imagery. A 90 m × 90 m sampling size on the original Landsat TM image was adopted to calculate the fraction of ISA (Powell et al., 2007;Roberts et al., 2012  impervious surface area was greater than 0 within the study areas, and was shown to vary within different parcel sizes. Spatial autocorrelation was lower with larger grain size, and spatial autocorrelation indices were in the range from 0.51 to 0.82, although lagging by one or two pixels (Fig. 4). Local Moran's I values of ISAs located in the Songhua River were less than 0 in both 2000 and 2010, but in general the local Moran's I index of ISA increased from 2000 to 2010 and represented dispersed to relatively clustered patterns, with high-ISApercentage regions being clustered while other ISAs were dispersed owing to grassland, water, and soil situated around the central study areas.

Correlation analysis
To examine the relation between LST and ISA abunda-nce, as well as between LST and ISA spatial association, we applied Pearson's correlation analysis (Table 3 and  . This indicates that ISA spatial configuration also significantly affects LST, although the effect is smaller than that of ISA abundance. Such a relationship, however, cannot be observed in the scatter plots (Fig. 5), since, for an individual pixel, the complex variation in LST in our study cannot be explained solely by the presence of impervious surface fractions or local Moran's I index of ISA. To address this issue, we grouped the pixels into different groups based on ISA fractional intervals. We controlled for the ISA fraction   Fig. 4 Moran's I of impervious surface area (ISA) for Harbin City by grouping observations in 10% intervals: 0-10%, 10%-20%, 20%-30%, 30%-40%, 40%-50%, 50%-60%, 60%-70%, 70%-80%, 80%-90%, 90%-100%. ISA percentages were additionally grouped in the intervals 0-5%, 5%-25%, 25%-50%, 50%-75%, 75%-100% (Table 4). Interestingly, the mean local Moran's I of the 75%-100% ISA fraction was found to be higher than that of the other fractions, indicating that high-ISA-percentage regions were clustered (Table 4). However, in both the 0-5% and 5%-25% ISA fractions, the local Moran's I was observed to be lower than in all other regions (local Moran's I = -0.06 and -0.024 respectively), indicating that the spatial pattern was random. For 0-5% ISA, the mean LST was 17.670℃. The mean LST increased substantially at ISA fractions greater than 50% (e.g., LST = 20℃ for ISA fractions >75%). For instance, for every 10% increase in the ISA fraction, LST would increase by 0.25℃ to 0.37℃, and every 0.1 increase in the local Moran's I index of ISA yielded 0.25℃ and 0.51℃ increase in LST from 2000 to 2010. A greater coverage of impervious surface will lead to an increase in LST, as will a clustered pattern of impervious surfaces, and there was a greater effect on LST in 2010, a time when the ISAs exhibited spatial heterogeneity. These results permitted visual comparisons of spatial patterns of ISAs among various types of urban environment with different ISA fractions. We found that higher LST measures were associated with clusters of higher-ISA-fraction categories (20% and higher), and lower LST measures were associated with lower-ISA-fraction categories (0-20%), whereas the spatial configuration did not play a major role (Table 4).

Effects of fraction and spatial pattern of ISA on LST
Because higher percentage cover of ISA tended to result in a more clustered pattern, the relatively strong correlations between local Moran's I index of ISA and LST were also attributed to the impervious surface fractions, when ISA was not controlled ( Table 5). Because of the weak correlation between ISA fraction, local Moran's I index of ISA, and LST in 2000 ( Fig. 6 and Table 5), we examined relationships between local Moran's I index of ISA and LST for different groups at the pixel level in 2010. All observations were separated into 15 groups using the 5%, 10%, 20%, or 25% fraction rule. Some groups had larger numbers of observations than others. Groups with more than 4000 pixels were considered for analysis because a small sample size could potentially lead to unreliable or unrepresentative results. This allowed us to accurately quantify the effects of the spatial arrangement of ISA features on LST. The Pearson correlation coefficient r between local Moran's I index of impervious surface area and LST ranged from −0.009 to 0.266 (Table 5). A weak correlation (r = −0.009) between local Moran's I index of impervious surface area (0-5% ISA fraction) and LST was found, but strong correlations (r > 0.25) were shown to be statistically significant for the 40%-50% and 50%-75% groups. We considered that the spatial arrangement of the ISAs had a substantial impact upon LST when r > 0.25. These results highlighted that clustered ISA patterns increased the LST more dramatically. The impact of ISA spatial arrangement on LST was minimal for areas with a low fraction of ISA (0-10%). However, for ISA >80%, the Pearson correlation coefficients were lower (−0.009) than for other ISA fractions. For 0-50% ISA, mean local Moran's I index of impervious surface area values ranged from −0.219 to 0.034, representing relatively dispersed to random patterns of ISAs (Table 5). The spatial pattern ranged from random to somewhat clustered (local Moran's I index of ISA with 0.004 and 0.407, respectively) in land-cover regions with 20%-80% ISA. However, in land-cover regions with > 80% ISA, we observed random and highly clustered spatial patterns (local Moran's I index of ISA with −0.027 and 0.535, respectively; Table 5). The spatial arrangement of the ISA had a significant impact on LST even when the ISA fraction was low (20%-30%). As we observed stronger correlations between local Moran's I index of ISA and LST (e.g., r = 0.266 and 0.253 for 50%-75% and 50%-60% ISA, respectively), we consistently found that warming effects were strongest when the ISA fraction exceeded 50%. Interestingly, ISA fraction had a negative correlation with LST when the landscape was dominated by 80%-100% ISA.

Discussion
This study applied the mono-window algorithm to retrieve LST in Harbin in the period from 2000 to 2010.
LST was normalized to analyze the UHI intensity within this timeframe. We found that urban warming environment levels significantly changed within the central city. ISA was extracted using the spectral unmixing method, and then a spatial autocorrelation index, the local Moran's I, was used to quantify dispersed or clustered patterns of ISAs. The spatial configuration of ISAs exhibited obvious differences at various scales, with high-percentage ISAs having a relatively clustered pattern. The magnitude of the impact of the spatial association of ISAs on LST depended on the ISA fraction. In 2010, the ISA distribution was spatially heterogeneous, and the local spatial association had a more significant relationship with LST. ISA fractions reached 50%, indicating that the spatial pattern of ISAs was very clustered, and it had a strong effect on LST. When the ISA fraction was greater than 50%, both this fraction and the spatial patterns of ISAs were significantly correlated with LST, so the redline value of ISA with regard to ecological impact should be 50% in urban areas of Harbin. Regions of greater than 80% ISA showed more clustered patterns; however, impervious surfaces had a minimal or even negative influence on LST in our study area. This was likely because some of the ISAs were composed of bright materials such as light-colored rooftops, and also because ISAs may provide shade, resulting in a cooling effect. The majority of ISAs observed included lightcolored rooftops, most of which were on new buildings, whereas older buildings were mostly constructed with dark-colored rooftops.
The results of this study suggest that city managers should target areas with large patches of ISA (> 50%) and control the ISA fraction to minimize the effects of land cover fraction on LST. They should also optimize the spatial pattern and configuration of ISAs, and replace old buildings with new ones to facilitate the cooling effect produced by shade and light-colored rooftops.
Numerous studies have focused on the relationship between the UHI effect and the increase in urban ISAs (Fan et al, 2019;Zhou et al, 2018). Based on this previous research, the impact on LST of both the fraction and spatial pattern of ISAs were examined in this paper. Furthermore, the impact of the spatial configurations of ISAs on UHI were explored. It was shown that LST is positively correlated with ISA abundance. The increasing abundance of ISAs contributes to increasing LST, which is consistent with the results of previous studies (Lei et al, 2019). Furthermore, it was shown that LST is significantly correlated with the spatial configuration of ISAs. When the fraction of ISAs reached 50%, their spatial pattern became clustered and had a strong effect on LST. The results of this study suggest that, in addi-tion to the abundance of ISAs, their spatial association has a significant effect on UHI. This extends the results of previous studies in an important way. It also provides scientific support for urban managers in their efforts to reduce the UHI effect by controlling ISA fractions and optimizing the spatial pattern and configuration of impervious surfaces.

Conclusions
Our results showed that LST has a significant positive association with both ISA abundance and the Moran's I index of ISAs, demonstrating that both the abundance and spatial clustering of ISAs contribute to elevated LST. Furthermore, LST is positively associated with clustering of high-percentage ISAs, while it is negatively associated with clustering of low-percentage ISAs, suggesting that UHIs are significantly affected by the spatial association of ISAs as well as by their abundance. This study considered the influence of the land component on the LST, but neglected the influences of ground elevation. The elevation as well as other factors need to be considered in future work to provide more precise estimates of LST. Furthermore, the appropriate spatial distributions of impervious surfaces, vegetation, water, and soil will need to be determined in order to optimize the urban thermal environment.