Automatic Detection and Dynamic Analysis of Urban Heat Islands Based on Landsat Images

Abstract

Given rapid global urban development, increases to impervious surfaces, urban population growth, building construction, and energy consumption result in the urban heat island (UHI) phenomenon. However, the spatial extent of UHIs is not clearly mapped in many UHI studies based on a remote sensing approach. Therefore, we developed a method to extract the spatial extent of the UHI during the period from 2000 to 2021 in Nanjing, China, and explored the impact of urban two- and three-dimensional expansion on UHI spatial extent and UHI intensity. After cropland effects (i.e., bare soil) were eliminated, our proposed method combines the Getis-Ord-Gi* and the standard deviation of the normalized difference vegetation index (NDVI STD) to extract the UHI area from Landsat 5 and Landsat 8 images using land surface temperature (LST) spatial autocorrelation characteristics and the seasonal variation of vegetation. Our results show the following: (1) Bare farmland has a large influence on the extraction results of UHI—combined with the seasonal variation characteristics of NDVI STD, the impact of bare soil on UHI extraction was highly reduced, strongly improving the accuracy of UHI extraction. (2) The dynamics of the UHI area are consistent with the changes in the built-up area in Nanjing at both spatial and temporal scales, but with the increase of the urban green ratio, the UHI area of mature urban areas trends to decrease due to the cooling effect of green space. (3) The accumulation of population and GDP promote the vertical expansion of urban buildings. When the two-dimensional expansion of the city reaches saturation, the UHI intensity is primarily affected by three-dimensional urban expansion.

1. Introduction

The phenomenon of the urban heat island (UHI) effect describes when the air and land surface temperatures (LST) are higher in an urban area than that in the surrounding rural area. This is presumed to be caused by the increasing in impervious surfaces, population growth, buildings, energy consumption, and decreasing green space in urban areas [1,2,3]. With accelerating urbanization, a large number of natural landscapes have been replaced by impervious surfaces, changing the distribution of solar radiation and aggravating the urban thermal environment [4]. UHI effects, enhanced by global warming and urban growth, have a large impact on the urban ecological environment, the health of residents [5], and the sustainable development of the economy [6,7,8,9,10].

Quantitative and qualitative monitoring of UHI is fundamental for UHI research [11]. Quantitative UHI measurements are usually characterized by atmospheric temperature and LST [12]. The traditional assessment for UHI is to identify the difference in air temperature measured from ground meteorological stations between urban and suburban regions [11,13]. This method has high temporal resolution, but lacks the potential for large-scale monitoring with spatial variability due to low densities of meteorological stations [12,14]. The expansion of urban built-up areas replacing urban green spaces gradually results in land surface energy accumulation and eventually leads to changes in LST [6,15,16]. Therefore, LST is the most direct manifestation of UHI [17], also making it an effective indicator of land surface energy balance, a key parameter of land surface physical processes at regional or global scales [12], and an important index for UHI evaluation [18,19]. Remotely sensed thermal infrared imagery has high potential to monitor LST at multiple temporal and spatial resolutions [11,13,20]. Compared to the conventional monitoring method using air temperature, remote sensing with much higher spatial resolution has great potential to effectively monitor long-term spatial variation in UHI, making it suitable for research at a large temporal and spatial scale [6,21]. Identifying the spatial extent of UHI is fundamental for both quantitative analysis and qualitative monitoring of UHI structure and change mechanisms. Generally, the method for mapping UHI is to divide the LST into multiple temperature groups (i.e., high, sub-high, medium, sub-low, low) first based on a threshold method and then identify the spatial extent of the group with the highest temperature as the UHI region [4,22,23,24]. However, the results from this method include many commission errors, as there are other land features with high LST including farmland without surface crops or other exposed soil [25,26]. More importantly, the recognition results of UHI varies from different threshold algorithms. Other studies identify the spatial extent of UHI by calculating the threshold of urban–rural temperature differences on the basis of defining urban–rural extents [3,27,28,29,30]. Peng et al. [31] used the radius method to determine the threshold of LST between urban and rural areas (i.e., the break point of temperature change between urban and rural areas) using the center of mass of the core urban area as the circle, thus improving the accuracy of UHI range identification. But with rapid and persistent urbanization, it is often difficult to clearly delineate the urban–rural boundary. Currently, some scholars use the spatial autocorrelation method to show the distribution of hot and cold spots of LST based on the LST so as to explore the spatial distribution characteristics of a UHI [25,32,33,34]. However, it is still a challenge to accurately map the spatial range of a UHI instead of describing the spatial distribution of LST. Therefore, it is important to develop a method to delineate UHI spatial extent accurately from thermal infrared imagery.

The spatial distribution of a UHI corresponds to that of built-up areas, including the urban central business districts, industrial areas, roads, and high-density residential areas [22,25]. The spatial distribution of a UHI is related to the morphological characteristics of the city, including two-dimensional morphological characteristics (horizontal direction) and three-dimensional morphological characteristics (vertical direction) [35]. Urban expansion is manifested as two-dimensional expansion and three-dimensional expansion of the city [17]. The two-dimensional expansion of the city is mainly reflected in the change in land cover types, which could be characterized by changes in indicators such as impervious surface ratio, soil ratio, normalized vegetation index, and building coverage ratio [17,36,37]. Among the studies of the driving factors of UHIs, most have focused on the relationship between land use and land cover change [20,38], urban impervious surface change [39,40], urban traffic development [41], and LST. The UHI intensity is positively correlated with the area of urban agglomeration and the proportion of built-up areas [15]. Those driving factors are mainly related to the expansion of UHI regions instead of UHI intensity. The three-dimensional expansion of the city depends on the change in building composition, which can be described by characteristic factors such as floor area ratio, sky view factor, and building height [17,36,42]. The study suggests that UHI intensity is positively correlated with building density, urban agglomeration area, and built-up area ratio [15]. The increase in urban building density hinders ventilation and enhances UHI effects [17]. When the urban population continues to grow, to alleviate the shortage of land resources, buildings expand vertically, and high-rise and super-high-rise buildings increase in number, promoting three-dimensional urban expansion [12,42]. Therefore, it is very important to clarify the influencing mechanisms of urban expansion (two- and three-dimensional) on the spatial range and intensity of UHI.

Therefore, our aim was to develop an automated method to map the spatial extent of UHIs and explore the impact of two- and three-dimensional urban expansion on the spatial extent and intensity of UHI. To achieve this goal, we take Nanjing city as a case-study research area and develop a method to extract the UHI range based on the spatial autocorrelation of LST. The specific research objectives are as follows: (1) combine the spatial autocorrelation of LST and the standard deviation of normalized difference vegetation index (NDVI STD) to map the spatial extent of UHI range, (2) analyze the spatial and temporal variation of UHI spatial extent and UHI intensity in Nanjing from 2000 to 2021, and (3) explore the impact of two- and three-dimensional urban expansion on the spatial extent and intensity of the UHI.

2. Materials and Methods

2.1. Study Area

Nanjing (31°14′–32°37′N, 118°22′–119°14′E; Figure 1), the capital of Jiangsu Province, China, is an important node city located in the strategic intersection of the Eastern Coastal and the Yangtze River economic belts [43]. Nanjing falls in a subtropical humid climate region with abundant rain, short spring and autumn, and long summer and winter [38]. There are 11 administrative districts in Nanjing, including Pukou, Qixia, Yuhuatai, Xuanwu, Gaochun, Liuhe, Lishui, Jianye, Qinhuai, Gulou, and Jiangning (Figure 1a). Nanjing is one of the four major “furnace” cities in China [38], describing the clear UHI effects (Figure 1b), making Nanjing a good case for UHI studies. The resident population of Nanjing in 2020 was 9.31 million, of which the urban population accounts for 8.08 million, and it has an urbanization rate of 86.80%. Compared with its population in 2010, the urban population in 2020 has increased by 8.86% (Nanjing Bureau of Statistics) [44] as the UHI effect has intensified. Based on the degree of urban expansion in each administrative district, Nanjing is divided into mature urban areas (Qixia, Xuanwu, Gulou, Qinhuai, Jianye, and Yuhuatai; Figure 1) and developing urban areas (Jiangning, Pukou, Lishui, Gaochun, and Liuhe; Figure 1). The main farmlands which may influence the extraction of UHI in Nanjing are rice and winter wheat croplands. The winter wheat is sown in the mid-to-late September and harvested at the end of May of the following year, while rice is mainly sown in mid-June and harvested in early November.

2.2. Datasets

2.2.1. Landsat Images

A total of 39 cloudless images from Landsat 5 and Landsat 8 collection-1 level-1 in Nanjing from 2000 to 2021 were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 5 August 2022; See the Supplementary Table S1 for detailed information about Landsat images, including the acquisition dates and the sensors). This includes 17 Thematic Imager (TM) images (2000–2011) and 22 Operational Land Imager (OLI) images (2013–2021). All collection-1 level-1 images obtained from the USGS website underwent both geometric and radiometric correction; the FLAASH module in ENVI 5.3 software was used for atmospheric correction.

2.2.2. Air Temperature Data

Daily air temperature data for Nanjing from 2000–2021 (April to October only) were obtained from the National Oceanic and Atmospheric Administration (NOAA; https://www.ncei.noaa.gov/maps/daily/, accessed on 14 September 2022).

2.2.3. Additional Data

The permanent population data, gross domestic product (GDP) data, greening rate data, and completed floor area of building of Nanjing were obtained from Nanjing Statistical Yearbook (http://tjj.nanjing.gov.cn/bmfw/njtjnj/, accessed on 5 July 2023). Nanjing floor data is obtained from the housing information website (https://nanjing.anjuke.com/, accessed on 5 July 2023). The surface coverage types in Nanjing are divided into five categories based on the world’s first global 30 m surface coverage fine classification product in 2020 [45] (the dataset is free to download at http://doi.org/10.5281/zenodo.4280923), including farmland, vegetation, water body, impervious surface, and wetlands.

2.3. Methodology

After preprocessing, including LST inversion, NDVI calculation, clipping to the main study area, and converting raster to points, the spatial autocorrelation analysis of LST was conducted using hot-spot analysis (Getis-Ord-Gi*). Hot spots with a high confidence interval (i.e., greater than 95%) were identified as UHI areas. However, the UHI areas extracted this way include bare soil areas in farmland. Therefore, the NDVI STD value was calculated using vegetation phenological characteristics to eliminate the effects of the exposed bare soil. Getis-Ord-Gi* was used to analyze the spatial autocorrelation of NDVI STD as well, and hot spots with confidence greater than 90% (i.e., the bare soil regions) were removed. Finally, the post-processing—including point densities calculation, pixel selection with a density greater than 0.0003, converting selected pixels to polygons, and removing small polygons—was performed to extract the UHI area. Factors influencing UHI spatial extent and UHI intensity were also explored using time series analysis (Figure 2).

2.3.1. Land Surface Temperature Inversion

LST were calculated using the inverse function of the Planck formula (Equation (1)) [17,36,46]. The atmospheric profile parameters were obtained using the unified calculation from the NASA website (http://atmcorr.gsfc.nasa.gov/, accessed on 5 August 2022). The calibration of band 11 from Landsat 8 was unstable, so band 10 was selected for the inversion of LST [4].

$$T_S=\frac{K_2}{\ln (K_1/B\left(T_S\right)+1)}-273.15$$

(1)

where T_s refers to LST (unit: °C), calculated from the radiance of blackbody in the thermal infrared band (band 6 of Landsat 5 or band 10 of Landsat 8) and calibration constants K₁ and K₂ (K₁ = 607.76 W·m⁻²·sr⁻²·μm⁻¹, K₂ = 1260.56 K for Landsat 5, and K₁ = 774.89 W·m⁻²·sr⁻²·μm⁻¹, K₂ = 1321.08 K for Landsat 8, which were obtained from the image metadata). $B\left(T_S\right)$ (unit: W·m⁻²·sr⁻²·μm⁻¹) is the land surface radiance of the blackbody having the real LST T_s (Equation (2)) [36].

$$B\left(T_S\right)=\left[L_{\lambda }-L_{\uparrow }-\tau \left(1-\epsilon \right)L_{\downarrow }\right]/\epsilon \tau$$

(2)

where L_λ (unit: W·m⁻²·sr⁻²·μm⁻¹) is the brightness value of thermal infrared radiation; τ, L_↑, L_↓ were obtained from the USGS website (https://atmcorr.gsfc.nasa.gov/, accessed on 5 August 2022) and are atmospheric profile parameters. τ is the atmospheric transmission of the thermal infrared band, L_↑ (unit: W·m⁻²·sr⁻²·μm⁻¹) is the atmospheric upward radiation brightness, and L_↓ (unit: W·m⁻²·sr⁻²·μm⁻¹) is the atmospheric downward radiation brightness [47]. ε is the surface emissivity calculated by vegetation coverage (Equation (3)) [46,48].

$$\mathsf{\epsilon }=0.004P_V+0.986$$

(3)

where P_V is the vegetation coverage (Equation (4)).

$$P_V=\left[\right(NDVI-{NDVI}_{Soil})/({NDVI}_{Veg}-{NDVI}_{Soil}\left)\right]$$

(4)

where NDVI is the normalized vegetation index (Equation (5)) [49], NDVI_Soil is the NDVI of bare soil area without vegetation coverage with a value of 0.05, and NDVI_Veg is the NDVI of vegetation with a value of 0.70 [46].

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$$

(5)

where the ρ_NIR is the surface reflectance of the near-infrared band (band 4 of Landsat 5 and band 5 of Landsat 8), and ρ_RED is the surface reflectance of the red band (band 3 of Landsat 5 and band 4 of Landsat 8).

2.3.2. Spatial Autocorrelation Analysis

UHI effects cause high LST clustered in urban areas, resulting in spatial clusters with relatively homogeneous characteristics in the LST images different from other land cover types [32]. Spatial autocorrelation analysis is a method to identify spatial clusters with relatively homogeneous phenomena [50] (e.g., LST and air temperature [21]). Getis-Ord-Gi* (also called hot-spot analysis), one of the most commonly used spatial autocorrelation analytical tools [6], uses the distance-defined space method to generate a spatial weight matrix [51]. From this matrix, significant spatial clusters with high and low values were extracted as hot and cold spots, respectively. In this study, the initial results of UHI spatial extent (i.e., the results before bare soil regions were removed) were extracted by Getis-Ord-Gi* (Equation (6)). UHI effects are magnificent in Nanjing during the period from May to October, so we selected images during this time to extract the UHI spatial extent. The hot spots (p ≤ 0.05; Gi_Bin ≥ 2, points with confidence interval greater than 95%) were extracted as the initial results of the UHI spatial distribution.

$$G_i=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^{N}W(i,j)x_i{x}_j}{{\sum }_{i=1}^N{\sum }_{j=1}^N{x}_i{x}_j}$$

(6)

where N is the total number of points, W(i, j) is the spatial weight matrix of the inverse distance method, and X_i is the LST of point i. X_j is the spatial neighbor of X_i, and G_i is the spatial correlation index (the larger the G_i value, the higher the degree of aggregation) [32].

2.3.3. The Standard Deviation of NDVI

After crops in farmland regions are harvested, exposed bare soil area become spatial clusters with high LST [25,52]. The phenology effects of vegetation results in clear seasonal variation in NDVI [49], the most commonly used vegetation index to evaluate the growth status and coverage of vegetation [53]. Therefore, the standard deviation of annual NDVI (Equation (7)) is useful to differentiate the bare soil regions of farmland from urban areas [51]. Due to Landsat images being limited within a year by cloud effects, NDVI STD was calculated from images acquired in two consecutive years to capture seasonal variation characteristics. Finally, hot spots of high NDVI STD (p ≤ 0.05; Gi_Bin ≥ 1, points with confidence interval greater than 90%) were calculate to be eliminated from the initial results of the UHI spatial distribution.

$$\mathsf{\sigma }=\sqrt{\frac{{\sum }_{i=1}^n{(x_i-\overline{x})}^2}{(n-1)}}$$

(7)

where n is the number of images used to calculate annual NDVI, $\overline{x}$ is the average NDVI value over n images, and x_i is the NDVI value of image i.

2.3.4. Post-Processing for Urban Heat Island Extraction

Post-processing included calculating point density, selecting the pixels with a density greater than 0.0003 (i.e., the number of midpoints per square meter), and converting raster to polygon. The small polygons (with an area less than 9,000,000 m²) were removed. The remaining polygons were considered UHI areas.

2.3.5. The Changes of UHI Area and Intensity

After extracting the spatial distribution of UHI, we explored the dynamics of UHI intensity and UHI spatial extents during 2000 and 2021 and their driving factors. The change in UHI area in each district of Nanjing was counted, and the ‘segmented’ piecewise linear relationship in R software was used to analyze the inter-annual dynamic change characteristics of UHI areas. This piecewise linear relationship is widely used in ecological research [50].

The dynamics of UHI intensity was analyzed by time series analysis of LST data inversed from Landsat thermal bands. First, LST of each month from May to October, 2000 and 2021 were calculated (if multiple images are obtained in the same month, the average value was used). For the missing values of LST for some months, it requires interpolation processing to further process data using time series analysis. The Kalman filter method in the “zoo objects” package in the R software package was used to fill in missing values and form a continuous series of LST data from April to October during 2000 and 2021. The seasonal Kalman filter can interpolate missing values from a dataset with repeated seasonal variability [54]. Time series analysis is then performed using the “ts” function in the “rts” package in R software to decompose the continuous LST series into trend, seasonal, and random noise components. With time series analysis, seasonal variation is minimized to create a consistent temporal trend for long-term study. Finally, we extracted the average value of the trend component for each year for further temporal analysis of UHI intensity.

The temporal changes of UHI intensity are not only affected by the local climate characteristics but also by two- and three-dimensional urban expansion. To estimate the effects of urban expansion on UHI intensity, air temperature was analyzed using time series analysis. The monthly average temperature from May to October in 2000 and 2021 in Nanjing was obtained by averaging the daily temperature within the same month in each year. The linear model (lm ()) in R software was used to analyze the trend of air temperature and LST (i.e., average value of the trend component from time series analysis) by simple linear regression.

3. Results

3.1. UHI Spatial Distribution Mapping

The initial results of UHI spatial extent determination (Figure 3b,e) suggest reasonable performance of spatial autocorrelation for LST (Figure 3a,d). However, most of the bare soil areas in Liuhe and Gaochun districts were mistakenly extracted as UHI areas (Figure 3b,e). Combined with the seasonal variation characteristics of farmland vegetation, the spatial autocorrelation analysis of NDVI STD significantly reduced the effects of farmland on the extraction of UHI spatial distribution, which significantly improves the accuracy of UHI extraction (Figure 3c,f).

3.2. Temporal and Spatial Variation of UHI Spatial Extent

The spatial and temporal variation of UHIs were analyzed from the UHI image extractions obtained on 10 October 2000, 12 July 2002, 20 May 2006, 18 May 2011, 9 October 2011, 18 May 2017, 9 October 2017, 6 June 2018, and 4 October 2021 (Figure 4a–f). In 2000, the UHI area of Nanjing was concentrated in the mature urban area (Figure 4a). By 2021, the UHI area of the mature urban area of Nanjing gradually expanded to the south and east in Qixia district, to the west in Jianye district, and to the south of Yuhuatai district (Figure 4i). With rapid urbanization, the spatial distribution of UHI in developing areas is expands correspondingly (Figure 5). The spatial extent of UHI in Jiangning district gradually expanded from the north to the south and west (Figure 4a–i). In Lishui district, the UHI distribution has primarily, but gradually, expanded from the central to the north and south sides (Figure 4a–i). The spatial extent of UHI expanded mainly from the east to the west along the river in Pukou district (Figure 4a–i). The UHI area of Liuhe and Gaochun districts primarily expanded from the southwest to the north, and from the center to the north and east, respectively (Figure 4a–i).

The segmented linear regression demonstrated that the UHI area of the mature urban area increased by 190 km² in 2021 from 2000, but with a downward trend after 2008 (Figure 5a). The UHI area of the developing urban area has been continuously rising from 2000 to 2021 (Figure 5b). Among all developing urban districts, the UHI area of Jiangning district has increased the most from 2000 to 2021 (381 km², Figure 6a). The UHI area of Pukou, Lishui, Liuhe, and Gaochun districts has increased by 153 km² (Figure 6b), 183 km² (Figure 6c), 140 km² (Figure 6e), and 91 km² (Figure 6d), respectively. Among these, however, the increasing rate of UHI area in Lishui district did not change significantly after 2011 (Figure 6c).

3.3. Dymaics of UHI Intensity

UHI intensity is characterized by LST, and the temporal trend of LST in mature urban areas and developing urban areas of Nanjing is generally consistent with the increasing air temperature trend in Nanjing (Figure 7). However, the increasing rate of LST in both mature and developing urban areas is significantly faster than that of air temperature (0.34, 0.35, and 0.03 for LST in mature urban areas, LST in developing urban areas, and air temperature, respectively). UHI intensity in mature urban areas is slightly higher (about 2 °C) than that of developing urban areas (Figure 7). And the increasing rate of UHI intensity in developing urban areas is slightly greater than that in mature urban areas (Figure 7).

4. Discussion

4.1. Extraction of Urban Heat Island by Getis-Ord-Gi*

LST in the UHI region is significantly spatially autocorrelated [21]. The initial results of mapping UHI spatial extent using LST hot-spot analysis show that a large amount of cropland in the north of Liuhe district and Baguazhou area were misidentified as UHI areas (Figure 3; see all the extracted results in Supplementary Table S2). Urban building materials and roads have low reflectivity and evapotranspiration rate, generating a large amount of anthropogenic heat and resulting in high LST in urban built-up areas [6,17]. Therefore, urban areas are shown as UHI in initial extraction results, as are farmlands as spatial clusters of LST when they are bare soil before crop planting and after harvesting [25,52,55]. Thus, many farmlands distributed in Liuhe district and Baguazhou were extracted incorrectly as UHI regions in the initial explorations (Figure 8).

Farmlands without crops have similar spectral characteristics as urban built-up areas, so they are also easily mapped as urban using classification methods of multispectral images [52]. The main crops planted in Nanjing are winter wheat and rice. Winter wheat is sown in mid-to-late September and harvested by the end of May. Rice is sown in mid-June and harvested in early November. Therefore, lands sown in winter wheat affect the initial extraction of UHI spatial distribution over the study period (May to October) (Figure 9a–d). There is a large amount of winter wheat planted in the Baguazhou area, which is routinely extracted as a UHI area, while rice croplands influence the preliminary results from May to June (Figure 9a,c,d); the northern part of the northern Liuhe district is also extracted as a UHI area due to the large amount of rice cultivation.

The UHI effect is significant in summer and autumn, so we select Landsat imagery from May to October to extract the UHI range. Due to the different growing period of crops (winter wheat and rice), some farmland areas are exposed in summer and autumn (i.e., the winter wheat has been harvested when rice is sown), which impacts the accuracy of UHI range extraction. Therefore, our study considers seasonal changes of vegetation NDVI because the bare soil area of farmland may have a high NDVI STD value. The NDVI STD is extracted as LST hot spots influencing the outcomes of the spatial autocorrelation analysis of NDVI STD. The significant hot-spot areas are the bare soil fields in farmland—the real UHI area can be determined by removing it. Most of the farmland in the north of Liuhe district and the bare farmland in Baguazhou can be accounted for by eliminating the NDVI STD hot spots (Figure 3c,f; Supplementary Table S2).

There are multiple error sources in the methods used to automatically extract heat islands that can reduce the accuracy of UHI extraction from Landsat images. Although the UHI effect is more obvious in summer and autumn, and the effect of UHI extraction based on summer and autumn images is significant, Nanjing has more precipitation in summer resulting in fewer clear summer images. This increases the uncertainty of time series analyses. Other error sources originate from variation in image quality and sensors (i.e., TM and OLI)—the methods used to extract UHI (e.g., spatial autocorrelation distance and small polygon elimination) [50]. Removing the small polygons with areas less than 9,000,000 m² after raster conversion to polygons during post-processing may result in the loss of some fragmented small UHI areas. The original spatial resolutions of the thermal infrared bands of Landsat 5 and Landsat 8 are 120 × 120 m and 100 × 100 m, respectively, and the difference in spatial resolution affects the spatial autocorrelation [56], which will influence the results of UHI extraction. To correspond to the spatial resolution of multispectral bands (i.e., for NDVI images), the thermal bands of Landsat 5 and Landsat 8 were resampled to a final image pixel size of 30 × 30 m using a triple convolution difference method. In this process, redundant image information is generated, which biases the observed correlation [26,43].

4.2. The Impact of Two-Dimensional Urban Expansion on UHI Area

The expansion of UHI distribution is generally consistent with two-dimensional urban expansion caused by population growth and economic development [40,55]. In 2000, UHIs were primarily found in mature urban areas, because commercial districts (e.g., Xinjiekou and Taipingmen) that were located in mature urban areas (Figure 4a) accumulated heat by attracting a large number of people and having densely distributed buildings [57]. From 2000 to 2021, the distribution of UHIs in mature urban areas expanded from the center to the surrounding areas (Figure 4), mainly related to the construction and development in the Hexi area (Jianye and Gulou districts) and the Xianlin University Park of the Qixia district. Intense construction in the Hexi area began in 2002, and the northern and central parts of Hexi area were fully completed by 2015 [58].

The increasing of UHI area in mature urban areas was dramatic before 2008, which is consistent with the population growth (Figure 10). Population in mature urban areas is quite stable during the period from 2008 to 2015 (Figure 10a). The UHI expansion in mature urban areas is primarily reflected in the eastern part of Qixia district, together with development of the Qixia district beginning with the Xianlin University Park in 2002 (Figure 4b). The UHI distribution in mature urban areas has a downward trend after 2008 (Figure 5a), which corresponds to the construction of urban green space and parks (the urban green ratio in Nanjing for 2005, 2010, 2015, 2018, 2019, and 2020 is 11.42%, 12.88%, 14.71%, 15.33%, 15.38%, 15.53%, and 15.62%, respectively). The increase in urban green space area has a cooling effect on UHI [44], potentially decreasing the trend and spatial expansion of UHI in mature urban areas. In addition, the original resolution of the thermal infrared bands of Landsat 5 and Landsat 8 is different (i.e., 120 m for Landsat 5 and 100 m for Landsat 8), which might lead to a slight decrease in the UHI range measured after 2012.

From 2000 to 2021, the UHI area in developing urban areas in Nanjing showed an overall upward trend (Figure 5b). After 2011, and the rate of expansion of UHI area also slowed (Figure 5b), which is consistent with the slower urbanization rate in the developing area of Nanjing [44]. The population of the developing urban areas became stable between 2010 to 2015 (Figure 10b). The Pukou district in developing urban areas even shows a downward trend in UHI area after 2014 (Figure 6b). When the developing urban regions become mature gradually (i.e., population of Jiangning, Pukou, and Gaochun districts became stable during 2010 to 2014; Figure 10c–e), the increasing urban green ratio would reduce the UHI area. Nanjing adjusted some administrative boundaries of the developing urban area in 2013, establishing the Lishui and Gaochun districts from what was Lishui and Gaochun counties [43], which promoted an increase in urban construction density and UHI expansion in the two districts (i.e., reflected in the population increasing rate after 2014 in the two districts; Figure 10e,f). The unevenly distributed Nanjing road network is primarily concentrated to the south of the Yangtze River, gradually reaching the urban periphery [59], which affects UHI effects in the developing urban areas of Nanjing. The accessibility and proximity of transportation infrastructure in Jiangning district (Lukou Airport and highways to neighboring provinces) also promotes land development and use [60]. Therefore, the decline in cultivated land and the expansion of impervious area in Jiangning district has resulted in the southward expansion of the measured UHI region (Figure 4) [60]. Jiangbei New District (new constructed district since 2015) [58], located in the northwest part of the Yangtze River in Nanjing, spans the urban area of Pukou, Liuhe, and Baguazhou districts, which form a UHI ‘strip’ along the Yangtze River (Figure 4f–i).

4.3. Influence of Three-Dimensional Expansion on UHI Intensity

The UHI intensity in Nanjing, described by LST, has been on the rise (Figure 7). The changing UHI intensity trend in the mature urban and developing urban areas is generally consistent with the changing air temperature trend. However, the increasing rate of UHI intensity is much higher than that of air temperature because air temperature is not the root cause of the UHI intensity [61,62]. Therefore, there must be other factors affecting UHI intensity change.

The increasing of UHI area slowed after 2011, and the urban expansion of mature urban areas was almost saturated between 2016 and 2021 [44]. However, UHI intensity in mature urban areas is still on the rise (Figure 7). The mature urban areas is districted with a major business circle, and the concentration of population has attracted more construction investment, which promotes the expansion of buildings in the vertical direction. Population and GDP are the influencing factors of urban expansion [63,64]. GDP is increasing year by year, and the GDP of mature urban areas is higher than that of developing urban areas, which is the same as the trend of UHI intensity (Figure 11a). Three-dimensional urban expansion is reflected in a change in building density and height, parameters that are positively correlated with LST [17]. The current studies indicate that the increase in floor area ratio, building volume, building height–width ratio, and building height promote the increase in LST [17,36,37,57,65]. There are differences in building height and density among administrative districts in the Nanjing mature urban areas. The area of completed housing construction in Nanjing has been on an increasing trend over the period 2000–2021 (Figure 11b), but urban expansion in the region tends to stabilize from 2016 to 2021. According to the information of urban buildings in Nanjing, the height of buildings in mature urban areas is increasing. The average number of floors of new buildings in Qinhuai district in 2000 was less than 10, and it reached 12 in 2008 and 21 in 2018. In addition, the building height in some areas (mainly along the Qinhuai River) of Qinhuai district is limited in order to protect the history and culture. The number of floors of new buildings in Qixia, Jianye, and Gulou districts has already reached 33 in 2000. Although the number of new floors in those three districts was mainly 33 in 2021, the number of such high rises increased. In 2000, the average number of new floors in Yuhuatai district and Xuanwu district was mostly less than 10, while it was greater than 20 in 2021. Meanwhile, the height of a single floor gradually increased as well. There is a positive correlation between building height and LST, and the increase in building height in mature urban areas promotes the increase of UHI intensity [57]. From 2000 to 2021, the UHI area in developing urban areas has been increasing and is greater than that in mature urban areas (Figure 5). The UHI intensity of mature urban areas is greater than that in developing cities, while the increasing rate of developing urban areas is slightly higher than that of the mature areas (Figure 7). The influence of three-dimensional expansion on UHI intensity is stronger than that of two-dimensional expansion, corresponding to the findings of Yang et al. [42]. While their increase is among the factors that lead to an LST increase, the shade cast by high-rise buildings also has a cooling effect; there is a given threshold when exceeded, reduces the cooling effect [12]. The proportion of super-high-rise buildings in Gulou district is very high (e.g., The height of Nanjing Zifeng Tower is 450 m), and the proportion of high-density buildings is only relatively high, while Yuhuatai and Qinhuai districts are home to mainly low-density with some medium-high density—but mostly low-rise—buildings [57]. With the increase in building height, the height difference of new buildings in mature urban areas is small. The uniform height of buildings will trap heat in a compact space, which can inhibit ventilation and increase UHI intensity [12,17].

5. Conclusions

The spatial autocorrelation of LST measured by combining Getis-Ord-Gi* and NDVI STD has high capacity to accurately extract UHI spatial distribution, laying a foundation for the follow-up studies of UHI.

Our research reveals two phenomena. First, urban population promotes urban expansion, and the changes in UHI distribution and urban built-up area in Nanjing are spatially and temporally consistent. With the increase of the urban green ratio, the UHI area of mature urban areas trends to decrease due to the cooling effect of green space. Second, the accumulation of population and GDP promote the vertical expansion of urban buildings. And UHI intensity is primarily affected by urban three-dimensional expansion of the city when the two-dimensional expansion of the city becomes saturated.