Diversification of Land Surface Temperature Change under Urban Landscape Renewal: A Case Study in the Main City of Shenzhen, China

Abstract

Unprecedented rapid urbanization in China during the past several decades has been accompanied by extensive urban landscape renewal, which has increased the urban thermal environmental risk. However, landscape change is a sufficient but not necessary condition for land surface temperature (LST) variation. Many studies have merely highlighted the correlation between landscape pattern and LST, while neglecting to comprehensively present the spatiotemporal diversification of LST change under urban landscape renewal. Taking the main city of Shenzhen as a case study area, this study tracked the landscape renewal and LST variation for the period 1987–2015 using 49 Landsat images. A decision tree algorithm suitable for fast landscape type interpretation was developed to map the landscape renewal. Analytical tools that identified hot-cold spots, the gravity center, and transect of LST movement were adopted to identify LST changes. The results showed that the spatial variation of LST was not completely consistent with landscape change. The transformation from Green landscape to Grey landscape usually increased the LST within a median of 0.2 °C, while the reverse transformation did not obviously decrease the LST (the median was nearly 0 °C). The median of LST change from Blue landscape to Grey landscape was 1.0 °C, corresponding to 0.5 °C in the reverse transformation. The imbalance of LST change between the loss and gain of Green or Blue landscape indicates the importance of protecting natural space, where the benefits in terms of temperature mitigation cannot be completely substituted by reverse transformation.

1. Introduction

The unprecedented rapid urbanization in developing countries, e.g., China, has been accompanied not only by extensive urban landscape renewal, but has also been linked to severe environmental and health problems [1]. The urban heat island (UHI), the best example of anthropogenic climate modification resulting from urbanization, is one of the most severe environmental problems affecting human health and ecosystem functioning [2]. Although there is no doubt that the form and expansion of urban landscape directly increases the threats of UHIs, the quantitative spatiotemporal variations of land surface temperature (LST) under urban landscape transformation are diverse [3,4,5,6]. This diversification could be caused by differences in seasons, hours, climate zones, and urbanization types [7,8,9,10]. Thus, the relationship between landscape pattern and LST is not consistent among various research conditions. How to quantify the effect of landscape transformation on LST change is still a question that must be answered by making multiple observations on various cities suffering high thermal risks.

Optimizing urban landscape planning, and designing the distribution of green infrastructure depending on the relationship between urban landscape and LST, are practical ways of relieving urban thermal environmental risk [11]. Within this practical objective, the first task is to clarify the correlation of landscape pattern on LST [12]. The spatial distribution of landscape and LST in Shanghai showed that residential land made the biggest contribution to the UHI, followed by industrial land [13]. A case study in Beijing using remotely sensed imagery observed that a 10% increase in the coverage of green space resulted in approximately a 0.86 °C decrease in LST [14]. Besides the area proportion of landscape type, the expansion of urban size was reported to be a significant factor affecting Beijing’s UHI [15]. A coincident study in Shanghai revealed the significant cooling effect of park size on LST [16]. Moreover, with the agreement about landscape composition influences, the cases in Beijing and Shanghai showed that landscape configuration could make an additional contribution to LST variations [17,18,19]. Knowledge of these relationships detected at the city scale can undoubtedly support the landscape management in China’s metropolitan areas.

Nevertheless, the correlation between landscape pattern and LST distribution cannot provide unique or perfect evidence to support landscape management. Notably, the spatial correlation of a static landscape pattern with LST cannot act as a complete substitute for the actual correspondence between landscape transformation and temporal LST change. In particular, the spatial difference of LST due to different landscape types cannot refer to the fact that the amplitude of LST change would be the same in the inter-transformation of two landscape types. Therefore, besides the evidence of a relationship between urban landscape pattern and LST distribution, the observed change of LST during actual landscape transformation should be closely examined. The most interesting observation from the previous case study of Beijing was that the amplitude of LST change might be very different as a result of green space gain and loss [20]. Whether Beijing’s LST change is an occasional effect of landscape transformation, or whether the evidence of this change can be confirmed in other metropolitan areas, are important questions that remain unanswered. If the evidence of a widespread effect could be established by multiple observations, knowledge about sustainable landscape management would be improved.

Among China’s three urban megaregions, Shenzhen is an economic center in south China, ranked third by economic size after Beijing (north China) and Shanghai (east China) [21]. Compared with Beijing and Shanghai, the subtropical climate in Shenzhen causes the city to suffer a higher risk of thermal environment. This study took the main city of Shenzhen as a case, and interpreted landscape types and retrieved LST from 49 Landsat images. The study objective was to determine the landscape transformation impact on LST change in the view of median amplitude, either spatially or temporally. The analyses were divided into two parts: (1) identifying the change trajectory of both urban landscape renewal and LST; and (2) determining the median amplitude of LST change due to landscape transformation.

2. Materials and Methods

2.1. Study Area

Shenzhen City is located on the south coast of China, east of the Pearl River and north of Hong Kong (113°46′–114°37′E, 22°27′–22°52′N). The territory area of the city is approximately 1997 km². After an adjustment of administrative divisions, Shenzhen has eight administrative districts (Bao’an, Nanshan, Longhua, Futian, Luohu, Yantian, Longgang, and Pingshan) and two functional districts (Guangming and Dapeng, which formerly belonged to Bao’an and Longgang, respectively) (Figure 1). In less than 40 years since China’s reform and opening-up policy, Shenzhen has evolved from a small town in the southern frontier to one of the economic centers of China. The gross domestic product of Shenzhen is almost equivalent to that of Hong Kong. Accompanying the rapid urbanization, Shenzhen has experienced drastic changes in landscape patterns [22], and the urban thermal environmental risk has continuously increased [23]. Because of the deterioration of ecosystem health during landscape transformation [24], Shenzhen has become the priority area for sustainable landscape planning [25]. In this study, the main city of Shenzhen, which is consistent with the width of Landsat images, was selected as the study area.

2.2. Datasets

Landsat satellites have been in use for more than 40 years. The thermal infrared radiation band of Landsat images has the advantages of a long period of record and good spatial resolution, making the images the most suitable for long time-series urban ecological remote sensing studies. The imaging of Landsat satellites in the mid-latitudes of the northern hemisphere is acquired at moderate sun angles (25–30°) during satellite travel in the same direction at the same time each morning, ensuring consistent observation conditions. Except for the panchromatic band and the thermal infrared band, the spatial resolution of the Thematic Mapper-5, ETM+, and OLI bands of Landsat 5–8 satellites is 30 m, and this is the basic spatial resolution of landscape pattern analysis adopted in this study. However, the spatial resolution of TM5 and ETM+ thermal infrared bands is 120 m and 60 m, respectively. The Landsat-8 thermal infrared band is derived from a TIRS Sensor, and the original data have a spatial resolution of 100 m, while the geometrically corrected data have a spatial resolution of 30 m with the effect of resampling in the image stack.

Because cloud and haze seriously affect the retrieval of LST data, it is necessary to select the cloudless clear-sky days for LST analysis. In addition, the full administrative area of Shenzhen City is captured on two Landsat images, and splicing the two unsynchronized images would result in uncertain LST values, especially for the overlapping areas in the two images. Therefore, only the Landsat images with the tile number 122/44 in the main urban area of Shenzhen City were selected (http://landsat.usgs.gov) in this study, of which 49 images were acquired during full clear-sky conditions with 0% cloud in the study area in the period from 1987 to 2015 (

Table 1. Acquisition dates of remote sensing data.

Time	Sensor	Time	Sensor	Time	Sensor	Time	Sensor
12-08-1987	TM5	11-15-1999	ETM+	03-09-2004	TM5	01-02-2009	TM5
11-24-1988	TM5	01-02-2000	ETM+	06-13-2004	TM5	02-03-2009	TM5
07-06-1989	TM5	09-14-2000	ETM+	10-19-2004	TM5	10-17-2009	TM5
09-14-1991	TM5	11-01-2000	ETM+	11-04-2004	TM5	11-02-2009	TM5
01-20-1992	TM5	03-01-2001	TM5	11-20-2004	TM5	03-26-2010	TM5
12-05-1992	TM5	11-20-2001	ETM+	09-16-2005	TM5	06-01-2011	TM5
12-24-1993	TM5	12-30-2001	TM5	11-23-2005	TM5	11-29-2013	OLI
01-25-1994	TM5	01-7-2002	ETM+	11-10-2006	TM5	10-15-2014	OLI
10-24-1994	TM5	11-7-2002	ETM+	12-28-2006	TM5	11-16-2014	OLI
03-03-1996	TM5	01-10-2003	ETM+	01-13-2007	TM5	10-08-2015	OLI
11-30-1996	TM5	01-18-2003	TM5	01-29-2007	TM5
08-29-1997	TM5	12-04-2003	TM5	03-04-2008	TM5
11-04-1998	TM5	01-21-2004	TM5	12-17-2008	TM5

). Due to the rainy weather associated with the subtropical climate, the number of suitable summer images was relatively limited. Summer was defined as May–September, together with winter as December–February, and spring and autumn (the transition season) as March–April and October–November, respectively. Accordingly, seven images were acquired in the summer, with 16 images in the winter, and 26 images in the transition season. Because nearly all of the plants in Shenzhen are evergreen and the winter season is warm enough to sustain plant growth, the temperature difference between urban and rural areas, namely the heat island phenomenon, exists in the winter, as well as in other seasons.

2.3. Landscape Types Interpretion

Taking advantage of the multi-spectral bands of Landsat images, hundreds of biophysical parameters have been derived to identify vegetation, soil, water, construction, rock, and other land covers from the imagery. In recent years, many mathematical statistical algorithms such as linear spectral separation and principal component analysis have been introduced to analyze remotely sensed data. To ensure comparability of the results with those from related studies, this study used the Normalized Difference Vegetation Index (NDVI) and Normalized Difference Built-up Index (NDBI) to describe vegetation information and non-vegetation information, respectively [26,27], and the Normalized Difference Moisture Index (NDMI) to describe water information. In addition, two relatively new remote sensing biophysical parameters were extracted for landscape type interpretation [28,29,30], i.e., the Automated Water Extraction Index (AWEI) and Biophysical Composition Index (BCI).

The NDVI is one of the earliest proposed biophysical parameters and has been widely used in remote sensing, geography, and ecology for more than 30 years. This study used NDVI as a biophysical parameter to identify the Green landscape, as shown in Equation (1).

$$NDVI=\frac{NIR-Red}{NIR+Red}$$

(1)

In Equation (1), NIR is the near infrared reflectance and Red is the red reflectance.

The identification of building construction and impervious surfaces has always been a challenge when applying biophysical parameters. Although the NDBI index does not meet the accuracy of NDVI for vegetation information, the NDBI is still the most widely used architectural biophysical parameter for resources and environmental studies, and is described by Equation (2).

$$NDBI=\frac{SWI{R}_1-NIR}{SWI{R}_1+NIR}$$

(2)

In Equation (2), NIR is the near infrared reflectance and SWIR₁ is the first shortwave infrared reflectance.

Similar to NDBI, NDMI is the most widely used water biophysical parameter even though it is not very accurate. Feyisa et al. [29] proposed a spectral operation rule (named the Automated Water Extraction Index, AWEI) that enhanced water information by separating background shadows, which had proved to be a practical parameter for extracting coastal zones [31]. To prevent the use of a single index leading to some uncertainties, both the AWEI (Equation (3)) and NDMI (Equation (4)) were extracted as cross-checks for water identification.

$$AWEI=4\times \left(Green-SWI{R}_1\right)-\left(0.25\times NIR+2.75\times SWI{R}_2\right)$$

(3)

$$NDMI=\frac{NIR-SWI{R}_1}{NIR+SWI{R}_1}$$

(4)

In Equations (3) and (4), Green is the green reflectance, NIR is the near infrared reflectance, SWIR₁ is the first shortwave infrared reflectance, and SWIR₂ is the second shortwave infrared reflectance.

In addition to the arithmetic band operations, statistical methods were used to extract multi-band information in a certain dimension as a means to enhance the effectiveness of biophysical parameters in presenting landscape types. Based on the brightness, greenness, and humidity information from a tasseled cap transformation, Deng and Wu [30] extracted three dimensions of the band combination layers to identify bright impervious surface, dark soil and mixed objects, and bright vegetation. The BCI composite parameter was built as shown in Equation (5). When the value of BCI is positive, higher values indicate a higher degree of impervious surface. When the value of BCI is negative, lower values indicate a higher degree of green vegetation. When the BCI value is near 0, the pixel is considered to have a mixture of both bare soil and vegetation. The BCI effectively compensates for the uncertainties of NDVI and NDBI, thus forming a mutual check.

$$BCI=\frac{\left(H_{cap}+L_{cap}\right)/2-V_{cap}}{\left(H_{cap}+L_{cap}\right)/2+V_{cap}}$$

(5)

In Equation (5), H_cap is the normalized result of the High Albedo layer, L_cap is the normalized result of the Low Albedo layer, and V_cap is the normalized result of the Green layer.

Spectral-based supervised classification and object-oriented classification are widely used for landscape type interpretation, and the object-oriented classification method has unique advantages in the interpretation of patch boundaries. However, taking into account the differences among the spectral features of the multi-modal images from different years and seasons, considerable work is required to extract the signature of training features from every image. If an interpretation is unsatisfactory and the signature must be re-adjusted, the amount of work is multiplied. In contrast, decision tree classification based on expert knowledge can be conducted without the need to determine the signature of an interpretation image. Although the accuracy of patch boundary determination is weaker than that using object-oriented classification, the efficiency of decision tree classification is significantly improved. If the results do not meet the requirements, the initial interpretation steps can be quickly revised, which significantly reduces the manual workload. Because the decision tree algorithm does not need prior spectral information, the flexibility of the algorithm is more suitable than other classification methods for analyzing remote sensing data with high temporal frequency.

However, the difficulty of using a decision tree algorithm lies in constructing and debugging expert knowledge. To ensure classification accuracy, this study divided a landscape into three types: “Blue”, “Green”, and “Grey”. The Blue landscape included the sea, rivers, reservoirs, and other water bodies. The Green landscape encompassed vegetation cover including woodland, grassland, garden, and farmland. The Grey landscape included non-vegetation covers except water bodies, including construction land, built land, sand, and bare bedrock. In this classification, the five selected remote sensing biophysical parameters can effectively identify the spectral differentiation characteristics of each landscape type, thus effectively reducing the possibility of misclassification. Furthermore, to reduce the uncertainty due to using a single parameter, two different biophysical parameters were set simultaneously as the classification criteria in the decision tree. In addition, to weaken the effect of seasonal fluctuation on the value of biophysical parameters, the double standard of absolute value and proportional value were combined to determine a landscape type. The thresholds were delineated by manual sampling of the test results (Figure 2).

Using the BCI, the integration of impervious surface information and vegetation information served as a starting point for judgment. The rules for interpretation were as follows:

• When the BCI < 0 or less than 40% of the order, a pixel may have vegetation information; if this pixel also has NDVI > 0.2 or greater than 50% of the order, the pixel is classified as Green. If the pixel fails the NDVI criterion, then go to the next step to judge the water.
• If AWEI > 0 or greater than 95% of the order, and then NDMI > 0 or greater than 50% of the order, the pixel is classified as Blue. If the pixel fails to pass the AWEI or NDMI criterion, then go to the next step to judge impervious surface and bare soil.
• When BCI > 0.1 or greater than 50% of the order, the pixel may have non-vegetation information. If the pixel also has NDBI > 0 or greater than 50% of the order, the pixel is classified as Grey; however, if the pixel fails one of the two criteria, it is classified as Unknown.
• For pixels classified as Unknown, the two nearest images in temporal dimension are used to correct the Unknown classification. If the pixel classified as Unknown is instead classified as the one of the three types (Green, Blue, or Grey) in the previous or subsequent images, the specific type is assigned to the Unknown pixel. However, if the pixel is classified as Unknown in all three images, the classification is assigned based on what landscape type comprises the highest proportion of the landscape, which is Green landscape in the application.

To eliminate the interference of a solitary cell, filter windows were set to be 19 × 19, 15 × 15, 11 × 11, 7 × 7, 5 × 5, and 3 × 3 pixels. These filters were applied three times. The filter window was the largest in the first filtering. Classification mistakes of the Blue landscape often appear in the shaded area of buildings, and thus the filtering radius was relatively large to avoid error caused by small building shadows. If a solitary pixel was less than 18/361, 14/225, or 10/121 of the Blue landscape window, the pixel was merged into the surrounding landscape type to ensure the contiguous distribution of water and removal of building shadows. The blend of Green landscapes or Grey landscapes is mainly influenced by mixed pixels of both landscape types. Compared to the concentrated urban green space, the Grey landscape distribution was slightly more dispersed due to the presence of roads, sandy land, bare soil, and bedrock. Thus the Green landscape window filter threshold was set as 8/225 or 6/121, and the Grey landscape window filter threshold was 7/225 or 4/121. In the second filtering, the Blue landscape window filter threshold was 7/49 or 4/25, the Green landscape window filter threshold was 6/49 or 3/25, and the Grey landscape window filter threshold was 5/49 or 3/25. In the final filtering, a 3 × 3 window was used to merge any solitary pixels. The decision tree algorithm was programmed in Matlab (Version 2014a, Mathworks, Inc., Natick, MA, USA).

2.4. Land Surface Temperature Retrieval

The change of surface temperature is the direct consequence of landscape pattern change affecting the surface energy balance, especially in urban areas with high landscape heterogeneity. In view of the limited number of climate stations, LST retrieval is the main observation to present the spatial distribution of urban temperature. The thermal infrared bands including TM Band 6, ETM+ Band 61, and TIRS Band 10 were used in this study. The LST retrieval process consisted of three steps, namely, radiation calibration and atmospheric correction, emissivity extraction, and LST conversion, which were described as follows.

(1)

Radiation calibration and atmospheric correction. The digitized value of the thermal infrared band can be converted to the top atmospheric radiation using Equation (6). The value of the gain and offset is already available in the image information, and the top atmospheric radiation is corrected to the blackbody radiation using Equation (7). The required parameters were determined using the United States Government’s space agency calculator based on the image information (http://atmcorr.gsfc.nasa.gov/).

$$L_{\lambda }=DN\times Multiple+Add$$

(6)

$$L_T=\frac{L_{\lambda }-L_{\mu }-\tau \left(1-\epsilon \right)L_d}{\tau \epsilon }$$

(7)

In Equations (6) and (7), DN is the digitized value of the band, Multiple is the gain, Add is the offset, L_λ is the top atmospheric radiation, L_T is the blackbody radiation, L_μ is the upward atmospheric radiation, and L_d is the downward atmospheric radiation. τ is the atmospheric transmittance, and ε is the surface emissivity. The measurement unit for all L series parameters is W/(m² sr μm).

(2)

Emissivity extraction. Surface emissivity is characterized by the proximity of surface thermal radiation and blackbody thermal radiation, which is mostly in excess of 0.9. As the interpretation of real objects often introduces interpretation error, the general NDVI is used as a surrogate surface. The conversion ratios from NDVI to emissivity were determined using Equation (8) [32]. The water emissivity was set as 0.995.

$$\epsilon =\{\begin{array}{c}0.979-0.035\rho red NDVI<0.2\\ 0.986+0.004Pveg 0.2\le NDVI\le 0.5\\ 0.99 NDVI>0.5\\ 0.995 Water\end{array}$$

(8)

In Equation (8), ε is the surface emissivity, NDVI is the normalized difference vegetation index, ρred is the red band reflectivity, and Pveg is the vegetation coverage based on NDVI normalization. The upper limit for NDVI is 0.99 and the lower limit is 0.

(3)

LST conversion. Depending on the blackbody radiation, LST can be converted according to the sensor scaling constant [33], as shown in Equation (9). In Landsat 5, 7, and 8, the constant K₁ is 607.76, 666.09, and 774.89 W/(m² sr μm), respectively, and the constant K₂ is 1260.56, 1282.71, and 1321.08 K, respectively.

$$LST=\frac{K_2}{ln\left(K_1/L_T+1\right)}$$

(9)

In Equation (9), K₁ and K₂ are the calibration constants and L_T is the blackbody radiation.

In addition, it should be mentioned that the above calculator cannot support the images before 2000 at present. Therefore, an alternative method has been introduced for 14/49 of the images in this study, which has the following three steps: converting to at-sensor spectral radiance, converting to at-sensor brightness temperature, and converting to LST [13,34,35,36]. Please see the references for detailed equations.

2.5. Spatial Statistics

The pixel changing time was identified for urban renewal mapping. Three spatial statistical indicators, namely the hot-cold spot, gravity center, and transect line, were applied to detect changes in the LST distribution.

(1)

Urban renewal mapping. In previous studies, the time track of landscape transformation has not received enough attention, and results have only reflected the number, proportion, and location of changes, while failing to show the gradual change process from the inner city to the outskirts over a long time. Based on the sequence of landscape types obtained through the interpretation of images in this study, a simple criterion (Equation (10)) for the time track of landscape transformation was developed. The premise of Equation (10) is that when checking one pixel in each image, if a certain pixel is the same (i.e., has the same digital characteristics) in the previous two images, and the landscape type of this pixel is the same in the next two subsequent images (but different from that in the previous image), the pixel is identified as the transformation point. The corresponding year of the transformation point is assigned. Under this criterion, the requirement that a pixel has the same characteristics in two consecutive images is used to weaken the seasonal disturbance. Nevertheless, if the sequence number of consecutive images is too long, the sensitivity of the transformation point identification will be weakened.

If (Num − 1 == Num − 2) & (Num == Num + 1== Num + 2) & (Num ≠ Num − 1)
Year = Numyear

(10)

In Equation (10), Num is the serial number of each image, the conditional statement is for the landscape type of each pixel, Numyear is the year that satisfies the judgment condition, and Year is the output of the transforming year for each pixel.

(2)

Hot-cold spot. Although LST is different in different times, this difference does not affect the spatial distribution of relatively high and low values. Because the spatial resolutions are inconsistent in different Landsat images, the sampling scale was set as a 360 m “fishnet” to smooth noise. Based on the ArcGIS (ESRI Inc., Redmonds, CA, USA) Spatial Statistics Toolbox, the Getis-Ord Gi* local statistics were used to extract the hot and cold spots of the LST in the vector fishnet [37].

(3)

Gravity center. Based on the ArcGIS Spatial Statistics Toolbox, the spatial geometric center of gravity of the LST in the fishnet was extracted using the Mean Center tool [38], and the LST was the weight value. Combined with all the gravity centers, the spatial movement of the LST gravity center can be tracked.

(4)

Transect line. To obtain an intuitive observation of the LST response to the landscape transformation, four digital survey lines were established in the directions of east-west, south-north, southeast-northwest, and northeast-southwest, intersecting at the Shenzhen Civic Center (Figure 3). The sampling interval was 120 m. Six images were sampled, which were acquired in December 1987, March 1996, January 2000, September 2005, March 2010, and October 2015.

3. Results

3.1. Presicion of Landscape Types Interpration and Urban Renewval Identification

Based on 500 randomly generated verification points, the interpretation of landscape type was validated using Google Earth^TM (Google, Inc., Mountain View, CA, USA) (Figure 3). In view of the time limit for historical high resolution images, 2010 and 2015 were selected as the verification years. The validation showed that the Kappa coefficient was 0.82 in 2015, and the total accuracy was 90.60%. Among the landscape types, the Green landscape interpretation had the highest accuracy (up to 96.08%), and the Blue landscape interpretation had the lowest accuracy (72.73%). For the 2010 validation, the Kappa coefficient was 0.88, and the overall accuracy was 94.00%. The highest classification accuracy in 2010 was also for the Green landscape (up to 97.50%). The error types were mainly rivers, roads, and road green belts, which belonged to mixed pixels or were solitary points in the medium-resolution remote sensing images.

Because the landscape type in each pixel may exhibit multiple changes in a long time sequence, only the starting change time point and the ending change time point were extracted (Figure 4). For the starting change time, there were obvious changes in 1989, 1992, 1994, and 1996, and the number of changed pixels in these years far exceeded the number in any year after 2000. Correspondingly, the end change time points were mainly found after 2005, and a large number of pixels changed in the years of 2009, 2013, and 2014. The rapid transformation of urban landscape not only occurred at the beginning of urbanization, but also after the promulgation of the Basic Ecological Line policy in 2005; moreover, the urban renewal of Shenzhen City was continuous. Judging landscape changes by the differences between the end and start years, approximately 62.77% of the pixels exhibited no change; these pixels were mostly undisturbed mountain vegetation, some stable water bodies, and small built-up areas. For the changed pixels, the mean time span between the start and end years was 11.2 years, but the quartile distribution ranged from five to 16 years. The difference between the start and end years was probably due to the various phases of construction, including vegetation removal (arable land was converted to bare land), land storage (bare land was covered by natural grass), the construction of impervious surface (grass-covered area was converted to hard surface), and built area greening (hard surface was converted to artificial grassland). In addition, long-term landscape transformation occurred during the pre-construction and post-renovation of city villages, resulting in a change sequence in which cultivated land was transformed to construction land, then to bare land, and finally to artificial grassland. This process corresponded to more than 10 years between the start and end years of landscape change. In short, under long-term continuous urban renewal, the transformation of the landscape occurred virtually uninterrupted.

3.2. Spatial-Temporal Variation of Land Surface Temperature

Although the values of LST were different at a daily time step, the spatial distribution of clusters of high-low differences should have been regular. The distributions of hot and cold spots in most of the years were relatively stable, except for the northern hot spots, which occasionally changed to be insignificant (Figure 5). Before 2006, the distribution of LST showed the regularity of hot temperatures in the east and south, and cold temperatures in the west and north. In the last 10 years, the spatial differentiation of hot-cold-temperature between the east and west has strengthened. However, the difference between hot and cold temperatures in the south and north were gradually replaced by the opposite pattern of hot in the north and cold in the south, reflecting the rapid urbanization in the north of Shenzhen in recent years since the cancellation of exterior administrative boundaries. At the same time, with the implementation of land reclamation projects, the western waters changed into hot spots from cold spots. Thus, the total number of hot spots was not significantly reduced, although some hot spots became insignificant in recent years. Before 2000, the hot spots were mainly connected to each other to form several contiguous zones for most of the time. However, with the rapid development of the western coastal area, the hot spots gradually expanded to the western and northern regions, leading to a fully integrated hot spot pattern after 2000. In contrast, depending on the construction of greening projects, hot spots in the early-developed parts of the core city became insignificant. Nevertheless, the widespread formation and expansion of hot spots indicated that the thermal environmental risk of Shenzhen is becoming more serious.

The geometric gravity center of LST moved slightly westward (Figure 6). Before 2005, the development of the core city was mainly due to the expansion to the surrounding area, and there was no obvious direction of movement. Although the LST hot spots were expanding, the distribution of relative difference values was stable. The gravity center of LST did not show a significant shift, and was mainly confined in an area of 35 m wide in the east-west direction, and 30 m wide in the north-south direction. The westward moving trend of the gravity center for LST was obviously exacerbated after 2005, corresponding to the rapid rise of hot spots in the west and the decline of hot spots in the central area (Figure 5). Although the gravity center of LST shifted westward, the overall distance was short. However, even though the western LST was slightly higher than before, the LST rise of other areas should not be neglected. The mitigation of thermal environmental risk is a concern not only for the western part; a landscape optimization plan for thermal environmental risk adaptation in the whole city is needed.

Samples of LST along the east-west transect showed the LST change caused by landscape transformation in the periods of 1987–1996, 1996–2000, and 2005–2010. Blue landscape transformed to Grey landscape, resulting in an LST rise (Figure 7, window a), and Grey landscape transformed to Green landscape, which led to an LST decrease (Figure 7, windows b and c). Samples along the south-north transect showed the LST changes in the periods of 1987–1996, 1996–2000, and 2010–2015, and reflected the transformation of Grey landscape to Green landscape leading to an LST decrease. The northeast-southwest transect also showed a consistent phenomenon in the periods of 1987–1996, 2000–2005, and 2010–2015, in which Blue landscape transformed to Grey landscape, leading to an LST rise (Figure 7, windows a and c), and Green landscape transformed to Grey landscape, also causing an LST rise (Figure 7, window b). In general, the transect analysis further confirmed the realistic LST changes caused by various kinds of landscape transformation, and fully demonstrated the possibility of thermal environmental risk adaption through landscape optimization. Moreover, the landscape transformation did not show a one-to-one correspondence with LST changes. That is, the landscape transformation directly caused the LST changes, but the LST changes were also largely affected by other factors besides landscape transformation. For example, heat generated by anthropogenic sources would affect the LST, and might result in the heat effect of the landscape transformation being offset. Thus in some situations, the landscape transformation would not be accompanied with corresponding LST change.

3.3. Land Surface Temperature Change Resulting from Landscape Transformation

Based on the aforementioned analysis, the LST fluctuation caused by seasonal differences cannot affect the overall consistency between landscape transformation and LST change. Accordingly, a logical question is “what are the differences observed among the effects of various landscape transformations on LST change?” As the seasonal fluctuation of LST was not easy to eliminate, the observations were based on a number of images. The average LST of the changed and stable landscape type in every two consecutive images was extracted. The average LST of the stable landscape type was subtracted from the average LST of the changed landscape type to determine an average change of LST (∆LST) due to landscape transformation. For example, for two consecutive images, the average LST of the water body transformed into green land was extracted, and the average LST of the stable water body in the two images was also obtained. Thus, the ∆LST of water body transforming to green land could be quantified. This analysis resulted in 48 groups of ∆LST for six types of landscape transformation (Figure 8).

The median of the LST changes showed that the increase in LST when the landscape transformed from Blue to Grey was as much as approximately 1 °C. In contrast, the cooling effect of the transformation from Grey landscape to Blue landscape was approximately 0.5 °C, which was much less than the warming effect resulting from the reverse landscape change. The Blue landscape transformation into Green landscape led to warming, and the reversed transformation led to cooling with nearly the same magnitude. This observation indicated that even though a large amount of vegetation was planted during the reclamation of coastal areas, the warming effect of this transformation was nearly inevitable. Moreover, it was noteworthy that, although the Green landscape transformed into Grey landscape led to an LST increase, the transformation of Grey landscape into Green landscape did not achieve a significant cooling effect (the median ∆LST was near 0). The warming effect resulting from green space reduction was greater than the cooling effect resulting from a green space increase, consistent with the conclusions of Sun and Chen [20]. Generally speaking, the temperature change resulting from different landscape transformations is not linearly reversible. The transformation of Green landscape to Grey landscape often happened in highly vegetated areas, which weakened the role of plant transpiration and caused a warming effect. In comparison, the transformation of Grey landscape to Green landscape often happened in flat unbuilt areas or in redeveloped brownfield sites, where the vegetation coverage was not extensive. Therefore, although the LST in Green landscape was largely lower than that in Grey landscape, the transformation from Grey landscape into Green landscape did not necessarily achieve a significant cooling effect. The vegetation cover density and impact threshold should be further considered in analyzing the effects of landscape changes on LST. In summary, the imbalance of LST changes between the loss and gain of Green and Blue landscapes indicates the importance of natural space protection, where the lost benefits of a converted natural landscape cannot be completely ameliorated by re-greening or re-bluing, at least in terms of temperature mitigation.

4. Discussion

4.1. Spatial Imbalance of Land Surface Temperature Change and the Planning Inferences

The LST deviations were grouped by the elevation and slope to show the spatial imbalance of LST change (Figure 9). From 1897 to 2015, the absolute values of LST deviations increased, either positively or negatively. At the elevations above 100 m and slopes greater than 15°, the LST changes were generally coincident, with the values mostly ranging from −1 °C to −3 °C. In contrast to the rapid urbanization that occurred in the plain, the negative LST deviations were strengthened after 2009. In the groups of 50–100 m elevation and 5–15° slope, the LST deviations were nearly coincident. The value of the deviations was near 0 until 2009, especially for the elevation group. The two curves in the 10–50 m group and 2–5° group of LST deviations oscillated prior to 2003, but tended to stabilize near 0.5 °C after 2004. Unlike the other three groups, the LST deviations in the 0–10 m and 0–2° groups showed an increasing trend that exceeded 0.5 °C after 2000 and more than 1 °C after 2011. It was also worth noting that the landscape types in the 0–10 m and 0–2° ranges contained a large amount of marine water, which had much lower LST; thus, the positive LST deviations reflected the huge impact of urbanization. In addition, the change of the landscape area and the change of LST did not fully correspond to each other in time, confirming that there were other factors that affected the spatial imbalance of LST.

In fact, after delineation of the Basic Ecological Line in Shenzhen, the high elevation and steep slope areas were located in a large proportion of the Basic Ecological Line. The LST deviations in different elevation and slope groups showed the obvious cooling effect that occurred within the ecological control line. However, the increase of the positive LST deviations in the flat urban area indicated that the risk to the thermal environment in the core area of the city might be increased. Therefore, although the Basic Ecological Line effectively protected the city’s green infrastructure, the mitigation of thermal environmental risk was far from sufficient. It can be concluded that the number of green spaces is not the only criterion that should be assessed for temperature mitigation. This study corroborated the opinion of Sun and Chen [20] who pointed out that the “ecosystem services of temperature mitigation were not equal between the loss and gain of green spaces even within the same area”. Besides controlling the city’s total number of green spaces through urban planning, other indicators should be introduced. Numerous studies in different countries have shown that the association of vegetation coverage and impervious surface parameters might be more closely related to thermal indicators than landscape area [39,40,41,42,43,44]. This study showed that the conversion between Green and Grey landscapes might not occur in the same intensity per unit area, resulting in a mismatch of LST increase and decrease. Faced with urban thermal environmental risk mitigation, landscape planning is necessary to emphasize the integration of quantity, quality, and space in ecosystem-based assessment and management.

4.2. Uncertainties of Landscape Types Interpretation and Land Surface Temperature Retrieval

In terms of the resolution of Landsat satellite images, the accuracy of this study met the requirements. Unfortunately, the spectral distinction between Blue landscape and urban architectural shadows was relatively poor (

Table 2. Confusion matrix and total accuracy in interpreting landscape types.

	Validation	Interpretation			Total	Accuracy
		Blue	Green	Grey
2015	Blue	16	4	2	22	72.73%
	Green	1	196	7	204	96.08%
	Grey	0	33	241	274	87.96%
	Total	17	233	250	500	90.60%
2010	Blue	17	2	5	24	70.83%
	Green	0	180	16	196	91.84%
	Grey	0	7	273	280	97.50%
	Total	17	189	294	500	94.00%

). To effectively distinguish Blue landscape from architectural shadows, this study set a relatively large range of pixel neighborhoods to remove “noise” caused by shadows, based on the fact that the architectural shadows were less prone to a concentrated distribution than the water surface. The problem faced by supervised classification was that a historical signature for interpretation was not easy to establish, and the interpretation of the recent points could not provide training samples for the interpretation of long-term images. This was the main reason for choosing decision trees for the interpretation of dozens of images. However, the drawback of the decision tree approach was that it was difficult to outline the detailed landscape types, such as distinguishing woodland, grassland, and farmland, as well as separating urban settlements, rural settlements, and facility buildings based on a 30-m resolution image. Under the trade-off between classification accuracy and classification type, this study adopted the most intuitive classification of “Blue”, “Green”, and “Grey", and the LST difference caused by vegetation or building types was left for analysis when higher resolution data become available.

Validating LST has been the core difficulty in the field of thermal infrared remote sensing. Whether observing LST from the ground surface or atmospheric profile, it is necessary to rely on ground measurements at verification sites for long-term observation. However, the historical data are difficult to obtain, and the precision of observations (especially historical ones) cannot be guaranteed. A weak but easy to implement method is the use of cross-validation. This study used the latest MODIS/Terra sensor V006 version MOD11A1 LST products (https://lpdaac.usgs.gov/), which are the world’s most widely used daily LST data. As Shenzhen is a subtropical coastal city with much water vapor in the atmosphere, MODIS daily LST products have a large number of missing values after quality control. Finally, eight MODIS LST images were chosen to correlate with the Landsat LST data (

Table 3. Correlation between land surface temperatures retrieved using Landsat and Terra (MODIS LST) satellites.

Date	Determination Coefficient	Samples	Date	Determination Coefficient	Samples
09-14-2000	0.832 **	1491	12-04-2003	0.758 **	215
03-01-2001	0.837 **	69	01-21-2004	0.787 **	183
11-07-2002	0.845 **	1505	10-19-2004	0.741 **	101
01-10-2003	0.823 **	1493	10-08-2015	0.822 **	326

** significant at p = 0.05.

). The numbers of MODIS samples for some dates were limited by the number of valid pixels. As the transit time of the different satellites is not exactly the same, the LST values obtained from different satellites cannot be exactly matched. In this study, the significance of correlations ranged from 0.74 to 0.85, indicating that the overall LST distributions were consistent, and the detailed local distributions were different. The retrieval algorithm and spatial resolution both induced differences between the LST values obtained from the two sources. In the absence of an LST on-site observation, it cannot be considered that Landsat LST was more distorted than MODIS LST. The consistency in values from the two sources may be used as an important prerequisite for assuring data effectiveness.

Additionally, it should be noted that the error of LST inversion is often no more than ±1 °C. In this study, 49 images were adopted and the median LST value had been extracted to remove the abnormal values. Although the same conclusion was reached in a former study [20], the research design in this study might be more strict due to the use of mass data. However, the result can still be improved as the error has not been completely excluded. The spatiotemporal spectral fusion of remote sensing images in cross validation may be a solution to reducing LST inversion error [45,46]. Meanwhile, this study has not considered the seasonal influence of the LST value, largely because Shenzhen is an ever-blooming city. However, the mean value of LST in winter is still obviously lower than that in summer in this city. It is considered that the variation of the mean value of LST may be the inherent result of seasonal influence. Thus, a consequential scientific question exists as to whether the relationship between landscape pattern and LST depends on the mean value of LST. A quantitative response on this question should be an important extension of this study.

5. Conclusions

In recent years, innovative remote sensing interpretation algorithms have been constantly developed. It is now necessary to explore how to reduce artificial participation in the interpretation of data to enhance the efficiency of interpretation. Especially for the data mining needs of massive geographic data sets, the artificial time cost should be further reduced. In this study, we developed a new decision tree algorithm, combined with the newly proposed BCI and AWEI indexes. Although the accuracy was not optimal, more than 90% of the interpretation accuracy was sufficient to meet the research objectives at the 30-m resolution. This decision tree algorithm has few indexes and no requirement for training signatures, making it especially suitable for analyzing high dimensional remote sensing images. In the study of landscape change, spatially rasterizing the time information is not a common process. This study constructed a judgment criterion to map the landscape transformation type, so as to identify the urban renewal process. The criterion provided a new way to display landscape change results. At the same time, based on hot-cold spots, the gravity center, and the transect line of LST, the multi-period LST change was characterized to produce multidimensional confirmation of the assumption that landscape change was a sufficient but not necessary condition for LST variation. In summary, this study confirmed the importance of landscape changes to LST in the main city area of Shenzhen, and deduced the important warming effect of urbanization. The results can provide strong support for the analysis of relationships between landscape patterns and LST.