Surface Urban Heat Island Analysis of Shanghai ( China ) Based on the Change of Land Use and Land Cover

In this paper, we present surface urban heat island (SUHI) analysis of Shanghai (China) based on the change in land use and land cover using satellite Landsat images from 2002 to 2013. With the rapid development of urbanization, urban ecological and environmental issues have aroused widespread concern. The urban heat island (UHI) effect is a crucial problem, as its generation and evolution are closely related to social and economic activities. Land-use and land-cover change (LUCC) is the key in analyzing the UHI effect. Shanghai, one of China’s major economic, financial and commercial centers, has experienced high development density for several decades. A tremendous amount of farmland and vegetation coverage has been replaced by an urban impervious surface, leading to an intensive SUHI effect, especially in the city’s center. Luckily, the SUHI trend has slowed due to reasonable urban planning and relevant green policies since the 2010 Expo. Data analyses demonstrate that an impervious surface (IS) has a positive correlation with land surface temperature (LST) but a negative correlation with vegetation and water. Among the three factors, impervious surface is the most relevant. Therefore, the policy implications of land use and control of impervious surfaces should pay attention to the relief of the current SUHI effect in Shanghai.


Introduction
Since the 20th century, urbanization has become the most significant human activity [1].The most intuitive expression of the rapid development of urbanization is the transformation of land cover types.Transformations in land use change the physical characteristics of the Earth's surface, affect the energy exchange between the ground surface and the atmosphere, impact the cycle of biogeochemistry, and have a profound influence on the structure and function of the regional or even global ecosystem [2].
Due to rapid urbanization, urban ecological and environmental problems have evoked widespread concern from the public, government and scientists.The urban heat island (UHI) effect is a most crucial issue, as its generation and evolution are closely related to social and economic activities.Studies on the distribution of UHI and its evolutionary mechanism have become a hot topic in multi-disciplines [3][4][5].Furthermore, UHI also leads to an urban rain island effect, which concentrates heavy rain during the flood season and causes regional water logging in megacities such as Shanghai [6].
UHI is a phenomenon in which the urban surface and atmospheric temperature are warmer than the surrounding non-urban environment [7].Usually, the temperature of the urban suburbs subtracted from that of the urban center acts as a measure of the intensity of the heat island.
Over the years, there have been many studies researching the causes [8,9], shape and structure [10], process and change [11], mechanism and simulation [12,13] of UHI formation.Urbanization has been shown to change the dynamic characteristics of the atmosphere and the heat exchange properties of the underlying surface, resulting in the rapid change of surface cover and land use and promoting the UHI [14].The larger the city and its population, the stronger the intensity of the UHI.The formation and development of an UHI is closely related to the geographical location and geometric shape of a certain city.In urban areas, factories, mines, enterprises, institutions and human activities releasing living heat promote the formation of UHI as well [15].
Many studies have shown that the formation of UHI and weather conditions have a strong correlation [16][17][18].UHI is closely related to the wind speed and varies with changes in the amount of clouds [19].Weather conditions such as sunny skies, quiet winds and low-pressure gradients can further intensify the UHI effect [17].
Except for the above air temperature components of UHIs, surface temperature components also matter considerably [20].Previous studies have also illustrated that SUHI is directly related to land surface types and surface modifications [21][22][23].Each type of land has its own thermal characteristics, radiation features and anthropologic heat, significantly affecting the interchange of surface energy, and then affecting the urban climate [24,25].For example, cement and tile structured buildings, squares, residential areas, bridges, roads and other urban land use types release more heat and cause higher temperatures, while bare land mostly consisting of soil, vegetation and water lead to lower temperatures [23].Therefore, with the expansion of the city and the changes in land use type, the SUHI effect will also produce corresponding changes.Vegetation can regulate energy exchange by transpiration.It is suggested that LST is associated with NDVI, while the results of relevancy vary considerably [26,27].In addition, most studies demonstrated a positive correlation between LST and IS [28][29][30].The effect of IS was inferior to that of NDVI [30], while some studies argued that IS had a strong correlation with SUHI, even in an exponential relationship [28].Therefore, how IS, vegetation and water affect energy absorption of the land surface and the extents of impacts need to be studied.
Remote sensing technology has been widely applied and has contributed much to assess SUHI with LST patterns from advanced very high resolution radiometer (AVHRR) [31], moderate resolution imaging spectroradiometer (MODIS) [32] to advanced spaceborne thermal emission and reflection radiometer (ASTER) [20].
Although many previous studies evaluated the LST and SUHI effect [23,33], there are few reports of SUHI in Shanghai [34] using satellite images in recent years.Shanghai, one of China's major economic, financial and commercial centers, has the highest level of urbanization in the country, which increased from 73.84% in 1999 to 89.8% in 2012, far more than the national average [35].In the first ten years of the 21st century, Shanghai has experienced a great change in land use.In addition, extreme hot weather has occurred in Shanghai with increasing frequency since that time.
In this study, we concentrate on analyzing the SUHI effect based on land-use and land-cover (LUCC) analysis in Shanghai using Landsat images.Therefore, this study will analyze the relationship between the impervious surface, land use and SUHI of Shanghai and draw general rules from its findings.
The LUCC is the key in analyzing the SUHI effect [22].The relationship between land surface temperature in Figure 1 and indicators like vegetation coverage and impervious surface area (ISA) will be the main research contents in this study.

Study Area
Shanghai is a municipality and one of the first open coastal cities in China.It is located at the confluence of the Yangtze River and Huangpu River and in the center of China's north and south coast.Shanghai is a part of the alluvial plain of the Yangtze River Delta.The Yangtze River Delta city group, comprising Shanghai, Jiangsu, Zhejiang and Anhui provinces, has become one of the six major world-class city groups.
Shanghai administers 16 municipal districts covering a total area of 6340 square kilometers.Shanghai, which has a subtropical humid monsoon climate, exhibits the characteristics of four distinct seasons, full sunshine and abundant rainfall.Table 1 gives detailed climate information for Shanghai which is presented on the website of EnergyPlus [36].All the climate data for Shanghai were obtained based on the decades of statistical climate data (1983 to 2010).Shanghai is China's economic, transportation, technological, industrial, finance, trade, exhibition and shipping center.The throughput of cargo and containers in Shanghai ranks first in the world.As the world's leading financial center, Shanghai's GDP ranked first among China's cities and second among Asian cities in 2015, second only to that of Tokyo, Japan.Shanghai is an immigrant city.Its unique conditions have attracted a large number of people over the long history of its development process as this sea town has become a world metropolis.By the end of 2016, the population of Shanghai was 24.197 million.Figure 2 shows the location of Shanghai in China.

Image Data
In this study, all raw image data of Shanghai were downloaded from the United States Geological Survey (USGS) website.Three phases of image data in 2002, 2007 and 2013 all have good quality with no cloud cover found in the selected area.

Image Data
In this study, all raw image data of Shanghai were downloaded from the United States Geological Survey (USGS) website.Three phases of image data in 2002, 2007 and 2013 all have good quality with no cloud cover found in the selected area.
Table 2 shows the raw data for Shanghai.The images are named after their imaging time by year and day of the year.Every image includes the information of its imaging sensor, date, resolution and wave bands.The OLI image consists of visible bands, the near infrared (NIR) band, thermal infrared (TIR) band and short wave infrared (SWIR) band, which are present in TM images, and also the coastal band, panchromatic (Pan) band and cirrus band.TIRS bands are also thermal infrared bands with a higher resolution compared with TIR bands.
In addition, LST retrieval requires radiometric calibration for thermal band.For TM images, thermal band is band 6.For OLI images, thermal band is band 10 and 11.The specific implementation can be seen in Section 4.4.

Radiometric Calibration
Radiometric calibration is used to determine the exact radiation brightness value at the sensor entrance and to further convert the radiance value to the outer surface reflectivity [38].The formula is presented as follows.
is the radiation value at the sensor entrance for band i.
is the brightness value of band i output by the sensor.
is the absolute calibration gain coefficient.is the absolute calibration bias value.The values of gain and bias are available in the header file of the remote sensing image.

Atmospheric correction
Atmospheric correction is used to convert the radiation brightness value or the outer surface reflectivity to the actual reflectivity of land surface, and the purpose is to eliminate the error caused by the atmospheric scattering, absorption and reflection.The method used in this project is based on the radiative transfer models.* ( , , , ) = ( , is the reflectivity formed by the path of molecular scattering and aerosol scattering.( , ) is the reflectivity formed by atmospheric absorption.S is the atmospheric spherical albedo.
is the reflectivity of the land target object.( ) is the scattering transmittance from the sun to the ground.and are solar altitude and sensor altitude, respectively.

Impervious Surface
The impervious surface is extracted from the vegetation coverage based on the normalized difference vegetation index (NDVI), combined with the application of the water mask extracted by the modified normalized difference water index (MNDWI).

Radiometric Calibration
Radiometric calibration is used to determine the exact radiation brightness value at the sensor entrance and to further convert the radiance value to the outer surface reflectivity [38].The formula is presented as follows.
L i is the radiation value at the sensor entrance for band i. DN i is the brightness value of band i output by the sensor.A i is the absolute calibration gain coefficient.B i is the absolute calibration bias value.The values of gain and bias are available in the header file of the remote sensing image.

Atmospheric correction
Atmospheric correction is used to convert the radiation brightness value or the outer surface reflectivity to the actual reflectivity of land surface, and the purpose is to eliminate the error caused by the atmospheric scattering, absorption and reflection.The method used in this project is based on the radiative transfer models.
ρ r+a is the reflectivity formed by the path of molecular scattering and aerosol scattering.T g (θ s , θ V ) is the reflectivity formed by atmospheric absorption.S is the atmospheric spherical albedo.ρ s is the reflectivity of the land target object.T(θ s ) is the scattering transmittance from the sun to the ground.θ s and θ V are solar altitude and sensor altitude, respectively.

Impervious Surface
The impervious surface is extracted from the vegetation coverage based on the normalized difference vegetation index (NDVI), combined with the application of the water mask extracted by the modified normalized difference water index (MNDWI).

NDVI
NDVI is an approach to assess whether the target on the ground surface has live green vegetation.The calculation of NDVI is relevant to the red band and near-infrared band because of spectral signatures of vegetation, which can be written as: The value of NDVI ranges from −1 to +1.The higher the result, the higher the density of green leaves.

MNDWI
In this study, the ratio method is applied to extract water information.The ratio method can make use of the difference of the object in different bands, and then highlight the information by the ratio calculation.The most common water index is the modified difference water index (NDWI) [39].NDWI is calculated using the difference of the green band and near-infrared band, which effectively eliminates the vegetation information to highlight the water information.However, NDWI neglects the influence caused by construction areas like commercial buildings and housing estates.Considering that urban area covers a lot of the study area, the modified normalized difference water index is applied in this study [40], whose calculation is represented as follows: The calculation of MNDWI helps to obviously display water regions.The pixel values of water tend to be higher than those of other ground objects.

Extraction of Impervious Surface
The dimidiate pixel model is a commonly used remote sensing estimation model to classify the spectral range of the sample as the dividing line to determine whether the terminal pixel falls into the spectral range of the classification sample [41].In the estimation of the vegetation coverage ratio, assuming that the pixels completely covered by the vegetation and soil as the dividing line, the vegetation coverage (V c ) expression can be obtained according to the NDVI, which reflects the information of vegetation growth status on the ground surface.
NDV I veg can be approximately equal to the maximum value of NDVI, and NDV I soil can be approximately equal to the minimum value of NDVI.Therefore, Equation ( 5) can be expressed as:

Land Use Classification
A decision tree based on the CRUISE algorithm is used to classify land use types.The full name of CURISE is classification rule with unbiased interaction selection estimation [42].It is a statistical decision tree algorithm used for data classification and data mining, in which there are four main features.Each node is divided into multiple child nodes, and the number of child nodes is the total number of response variable classes.The bias when the variable is selected is negligible.The algorithm has a variety of ways to handle missing values.The algorithm can detect the local interaction between the predicted variables.

Land Surface Temperature
The theoretical method of land surface temperature retrieval is to solve radiative transfer equations, eliminate atmospheric impact, and then obtain the land surface temperature.There are three kinds of commonly used methods: single band algorithms [43][44][45], split window algorithms [43,46,47] and multi-angle algorithms [48,49].Single band algorithms are effective and need fewer parameters and are thus are used in this study.
Additionally, taking into account the limited scope of the study area and considering that the remote sensing images were taken under clear and cloudless weather conditions, the degree of atmospheric impact in space that is consistent and the relative temperature of the ground temperature distribution would not be affected.Therefore, an image-based inversion algorithm using thermal infrared band is chosen in the study [50].Specific steps are as follows; the first is to convert the DN value to the radiance value.
L is the radiation value at the sensor entrance of the thermal band.DN is the data value of image pixels.gain and bias are the gain coefficient and bias value of thermal band, respectively, which are available from the header file of the remote sensing image.
Second is to convert the radiance value to brightness temperature.
T b is the brightness temperature, K 1 and K 2 are preset constants before launch.Third is to calculate the land surface temperature.
λ is the wavelength of emission radiation and ε is emissivity.Due to the complexity of the underlying surface type, for the thermal infrared band with resolution of 60 meters, a pixel corresponding to the underlying surface often contains a variety of materials.Various materials have different emissivity values.Correspondingly, the calculation of emissivity is quite complicated.The following method is used to calculate the emissivity.Firstly, the remote sensing image is divided into three types: the construction land, the water body and the natural ground surface.The emissivity of water pixels is assigned to 0.995.The emissivity of construction land and ground surface is calculated according to formulas (11), ( 12) and ( 13) [50].For a simple calculation, NDV I soil is defined as 0, and NDV I veg is defined as 0.7.12)

Experiment
In the study, the ENVI 5.1 software helps to process and analyze the geospatial remote sensing images and includes both a new interface and classic tools.

Preprocessing
Due to the limited geographical coverage of a remote sensing image, the Shanghai area should be made up of two remote sensing images.To ensure the comparability of images from different sensors or from the same sensor on different dates, and to eliminate the radiation error caused by atmospheric scattering, the images required radiometric calibration and atmospheric correction.In ENVI, both of these procedures have corresponding operation modules in the toolbox.Then, the projected boundary of Shanghai was imported for the selected regions of interests.Figure 4 shows the result of image preprocessing and is displayed in true color.

Impervious Surface Calculation
The calculation of NDVI is a ratio operation between the red band and near infrared band.There is no need to input the computational formula in ENVI due to the available tool in the toolbox.For TM images, band 3 is the red band and band 4 is the NIR band, and for OLI images, band 4 is the red band and band 5 is the NIR band. Figure 5 shows the result of NDVI.The lighter or white pixels have a higher probability of live green leaves.Different from NDVI, there is no ready function for MNDWI.Therefore, it has to be done in band math.The input formula is (b1 − b2)/(b1 + b2).Both for TM and OLI images, b1 is assigned as the green band and b2 is assigned as the middle infrared band.Figure 6 shows the result of MNDVI.The darker or black pixels are regions of water.

Impervious Surface Calculation
The calculation of NDVI is a ratio operation between the red band and near infrared band.There is no need to input the computational formula in ENVI due to the available tool in the toolbox.For TM images, band 3 is the red band and band 4 is the NIR band, and for OLI images, band 4 is the red band and band 5 is the NIR band. Figure 5 shows the result of NDVI.The lighter or white pixels have a higher probability of live green leaves.

Impervious Surface Calculation
The calculation of NDVI is a ratio operation between the red band and near infrared band.There is no need to input the computational formula in ENVI due to the available tool in the toolbox.For TM images, band 3 is the red band and band 4 is the NIR band, and for OLI images, band 4 is the red band and band 5 is the NIR band. Figure 5 shows the result of NDVI.The lighter or white pixels have a higher probability of live green leaves.Different from NDVI, there is no ready function for MNDWI.Therefore, it has to be done in band math.The input formula is (b1 − b2)/(b1 + b2).Both for TM and OLI images, b1 is assigned as the green band and b2 is assigned as the middle infrared band.Figure 6 shows the result of MNDVI.The darker or black pixels are regions of water.Different from NDVI, there is no ready function for MNDWI.Therefore, it has to be done in band math.The input formula is (b1 − b2)/(b1 + b2).Both for TM and OLI images, b1 is assigned as the green band and b2 is assigned as the middle infrared band.Figure 6 shows the result of MNDVI.The darker or black pixels are regions of water.According to Section 3.2.3, the formula of impervious surface calculation in band math is 1 − (b1 − NDVImin)/(NDVImax + NDVImin), where b1 is the NDVI result in Section 4.2.1.Because this result has not eliminated the influence of water, masking of water should be applied.The binariation formula is (b1 ge threshold) × 0 + (b1 lt threshold) × 1.Here, b1 is the result of MNDWI calculation.The threshold is 0.92, 0.91 and 0.93 for 2002, 2007 and 2013, respectively.After selection of water, the water mask can be applied on the calculation result of reversed vegetation coverage.Based on the ground features of the study area and combining visual interpretation with high resolution of QuickBird images in Shanghai, the ground objects are divided into three categories: impervious surface, water area and vegetation cover area; then, the training samples are determined.Using QuickBird images as reference data, randomly generate 200 sampling points in the three phases of images and check their categories.The random sampling points are compared with the corresponding points in the QuickBird images at the same latitude and longitude to figure out the similarities and differences.Table 4 shows the accuracy.According to Section 3.2.3, the formula of impervious surface calculation in band math is 1 − (b1 − NDVI min )/(NDVI max + NDVI min ), where b1 is the NDVI result in Section 4.2.1.Because this result has not eliminated the influence of water, masking of water should be applied.The binariation formula is (b1 ge threshold) × 0 + (b1 lt threshold) × 1.Here, b1 is the result of MNDWI calculation.The threshold is 0.92, 0.91 and 0.93 for 2002, 2007 and 2013, respectively.After selection of water, the water mask can be applied on the calculation result of reversed vegetation coverage.Figure 7 shows the result of IS extraction, which demonstrates the expansion of urban impervious surface in Shanghai.The percentages of imperious surface are 19.47%,36.23% and 37.09% in 2002, 2007 and 2013, respectively as Table 3 displays.Based on the ground features of the study area and combining visual interpretation with high resolution of QuickBird images in Shanghai, the ground objects are divided into three categories: impervious surface, water area and vegetation cover area; then, the training samples are determined.Using QuickBird images as reference data, randomly generate 200 sampling points in the three phases of images and check their categories.The random sampling points are compared with the corresponding points in the QuickBird images at the same latitude and longitude to figure out the similarities and differences.Table 4 shows the accuracy.

Accuracy Assessment
Based on the ground features of the study area and combining visual interpretation with high resolution of QuickBird images in Shanghai, the ground objects are divided into three categories: impervious surface, water area and vegetation cover area; then, the training samples are determined.
Using QuickBird images as reference data, randomly generate 200 sampling points in the three phases of images and check their categories.The random sampling points are compared with the corresponding points in the QuickBird images at the same latitude and longitude to figure out the similarities and differences.Table 4 shows the accuracy.In this experiment, a decision tree classification based on the CRUISE algorithm is used.Before building a new decision tree, classification training samples should be specified.There are four classes, which are vegetation, water, urban and others.In these four classes, urban is defined as an impervious surface while other three classes are pervious land.
Then, the decision trees can be created with the application of the plug-in RuleGen.The trees created by RuleGen tend to be too complex and need to be trimmed.
Figure 8   In this experiment, a decision tree classification based on the CRUISE algorithm is used.Before building a new decision tree, classification training samples should be specified.There are four classes, which are vegetation, water, urban and others.In these four classes, urban is defined as an impervious surface while other three classes are pervious land.
Then, the decision trees can be created with the application of the plug-in RuleGen.The trees created by RuleGen tend to be too complex and need to be trimmed.

Confusion Matrix
In this study, the high-resolution QuickBird images are also used to examine the land use classification result of 2002, 2007 and 2013.In ENVI5.1, this procedure can be realized by the module of the confusion matrix using ground truth images.Selecting ground truth images and adding the combination with the four classes lead to the results.Table 5 displays the classification accuracy.Figure 9 shows that land surface temperature retrieval is a complicated process.

Confusion Matrix
In this study, the high-resolution QuickBird images are also used to examine the land use classification result of 2002, 2007 and 2013.In ENVI5.1, this procedure can be realized by the module of the confusion matrix using ground truth images.Selecting ground truth images and adding the combination with the four classes lead to the results.Table 5 displays the classification accuracy.To calculate vegetation coverage is a little bit different from Section 4.2.1.The formula input in band math is (b1 gt 0.7) × 1 + (b1 lt 0) × 0 + (b1 ge 0 and b1 le 0.7) × ((b1 − 0)/(0.7 − 0)), where b1 is the NDVI.Then, the calculation result of vegetation coverage is obtained.
The radiance value is relevant to emissivity and the radiometric calibration result for the thermal band.The formula is represented as below, (b2 − 0.86 − 0.87 × (1 − b1) × 1.42)/(0.87× b1), and b1 is emissivity when b2 is radiometric calibration result for thermal band.Then, the calculation result of radiance value is obtained.

Impervious Surface Distribution in Shanghai
Impervious surface is one of the main types of land cover in urban areas.It is also an important component of urban ecosystems.
By a series of processes involving the impervious surface, binarization diagrams in different years were obtained (Figure 7).These three pictures vividly revealed that there is a rapid expansion of impervious surface in Shanghai, especially in the city's center.
In general, the statistical data present an ascending trend from 2002 to 2013.In detail, there is a booming from 2002 to 2007 in terms of proportion, after which the figure for impervious surface remains stable with minute growth from 2007 to 2013, as shown in Table 3.
The expansion of the built-up area depends on the balance between increasing urban construction land scale and the limited environmental capacity of proper conditions, and its developing direction concentrates in a special region.Therefore, the study area is divided into 10 sub-regions according to the administrative division: downtown (including Huangpu, Xuhui, Changning, JingAn, Putuo, Hongkou and Yangpu seven districts), Pudong (include old Pudong and Nanhui), Minhang, Baoshan, Jiading, Jinshan, Songjiang, Qingpu, Fengxian and Chongming districts.The following data analysis will be deeper based on these 10 sub-regions.
After clipping the result for the impervious surface, the entire study area was divided into 10 sub-regions.By binarization processing with the suitable threshold value, the three impervious surface statistics information data were obtained for 2002, 2007 and 2013 in Table 7. Statistical data of the generated cylindrical statistical diagram can directly reflect the diverse regional impervious surface growth conditions.The development condition of the impervious surface in each period within the study area reflects the significant difference.The transformation of Pudong new district, Jiading, Jinshan, Songjiang, Qingpu and Fengxian districts is most distinct among all sub-regions, which showed an almost twofold increase from 2002 to 2007.Generally, the tendency between 2007 and 2013 is a very small increase in all areas.There is a mild and regular rise in all the non-downtown areas with a great increase in Pudong new district.Specifically, the change in downtown is with a tiny decline partly due to more vegetation covers after the 2010 Expo and its environmental improvement plan.
Table 8 shows the impervious surface proportion table of 10 sub-regions.In order to proceed with further research of the transformation law based on differences in the regional impervious surface information over time, the three phase of images are divided into two periods (namely 2002-2007, 2007-2013) so we can understand the urban growth over time more specifically.Table 9 lists the impervious surface area growth in two periods.Apparently, there is only a minute decline in downtown areas; conversely, five figures (>100) reflect a dramatic increase in Pudong, Minghang, Jiading, Songjiang and Fengxian districts.Based on the statistical data above, the growth percentage of the impervious surface in the study area within each period can be calculated, which could represent the overall growth rate of the urban area within two periods.According to Table 10, both downtown and other districts display an obvious increasing rate from 2002 to 2007.In particular, the figures of Pudong new district, Jiading, Jinshan, Songjiang, Qingpu and Fengxian even exceed 100%.From 2007 to 2013, apart from the slight decline for the downtown, the other nine districts all gain slight increase in the coverage of impervious surface.Generally speaking, with the rapid development of urbanization, the density of buildings in urban built-up areas in Shanghai is quite high, while the greening rate is low.Statistics indicate that the average impervious surface rate in 2013 inside Shanghai's outer ring road reached 70%, much higher than that of most foreign mega-cities (around 40%).
From the spatial distribution map of impervious surface, it is estimated that the overall impervious surface rate of Shanghai is generally high, especially within the inner ring of the city center, and its impervious surface rate has reached over 90%.In addition, Zhabei and Wusong, industrial zones along Huangpu River, and some of the new industrial areas in Pudong new district have also climbed to more than 85%.There is a significant difference between the urban land outside the inner ring and the downtown area, and its impervious surface rate of ranges from 50% to 85%.Besides, parks and green space in the urban area and the remaining farmland in Pudong new district have an impervious surface rate below 50%.

Urban Development Density Analysis
Impervious surface is an important indicator of urban spatial distribution and development density.A previous study reported the characteristics of thermal environment and city development condition in Tampa Bay and Las Vegas, the United States [51], in which the selected regions were the impervious surface rate that included more than 10% as urban land (e.g., 10% to 40% as low, 40% to 60% as medium and 60% to 100% as a high development density area).
Shanghai as one of the largest cities in China has been under high-density development for a long time.Its impervious surface rate is higher than other cities at home and abroad.Therefore, the threshold is different from the above rules.In this study, urban development density is divided into four degrees.Regions whose impervious surface rate is less than 40% are considered as non-urban land, mainly including bare land, farmland and water.Forty to 60% is considered a low development density area.Sixty to 80% is a medium development density area.Eighty to 100% is a high development density area.Selecting the downtown area as regions of interests, Figure 11 shows the distribution of urban development density.the average impervious surface rate in 2013 inside Shanghai's outer ring road reached 70%, much higher than that of most foreign mega-cities (around 40%).
From the spatial distribution map of impervious surface, it is estimated that the overall impervious surface rate of Shanghai is generally high, especially within the inner ring of the city center, and its impervious surface rate has reached over 90%.In addition, Zhabei and Wusong, industrial zones along Huangpu River, and some of the new industrial areas in Pudong new district have also climbed to more than 85%.There is a significant difference between the urban land outside the inner ring and the downtown area, and its impervious surface rate of ranges from 50% to 85%.Besides, parks and green space in the urban area and the remaining farmland in Pudong new district have an impervious surface rate below 50%.

Urban Development Density Analysis
Impervious surface is an important indicator of urban spatial distribution and development density.A previous study reported the characteristics of thermal environment and city development condition in Tampa Bay and Las Vegas, the United States [51], in which the selected regions were the impervious surface rate that included more than 10% as urban land (e.g., 10% to 40% as low, 40% to 60% as medium and 60% to 100% as a high development density area).
Shanghai as one of the largest cities in China has been under high-density development for a long time.Its impervious surface rate is higher than other cities at home and abroad.Therefore, the threshold is different from the above rules.In this study, urban development density is divided into four degrees.Regions whose impervious surface rate is less than 40% are considered as non-urban land, mainly including bare land, farmland and water.Forty to 60% is considered a low development density area.Sixty to 80% is a medium development density area.Eighty to 100% is a high development density area.Selecting the downtown area as regions of interests, Figure 11 shows the distribution of urban development density.From Figure 11, it is obvious that most area within the inner ring belongs to the high development density areas.Tongji University, People Square, Century Park and the old residential area near Xujiahui are formed into some rarely seen medium development density areas and low development density areas within the inner ring.Medium development density areas are mainly located in the new towns outside the inner ring, except for Wusong Industrial Zone and Zhabei Industrial Zone.Low development density areas are mainly located in Pudong new district and Dachang Town in Baoshan district.In general, the density of urban development is significantly different between the inner ring and outer ring.

Land-Use and Land-Cover Change
As illustrated in Section 2.2, the images ownloaded in this study were all taken in winter, from December, January and February.Therefore, the accuracy of land use classification depends on the actual situation on that day.The result of land classification (Figure 8) is obtained by the means of a decision tree, which vividly showed a change in urban, water, vegetation and other land (mainly soil), especially in the urban area.
In detail, the data from the Table 11 reflect an increase in terms of urban and water areas from to 2007, while the figure for urban area declines from 2007 to 2013.The trend also complies with the result of impervious surface extraction.From the proportion diagram (Figure 12), the trend of the greater decrease in vegetation area is quite obvious, as it was from 2002 to 2007.The change in water stays in acceptable ranges.

Land-Use and Land-Cover Change
As illustrated in Section 2.2, the images ownloaded in this study were all taken in winter, from December, January and February.Therefore, the accuracy of land use classification depends on the actual situation on that day.The result of land classification (Figure 8) is obtained by the means of a decision tree, which vividly showed a change in urban, water, vegetation and other land (mainly soil), especially in the urban area.
In detail, the data from the Table 11    Through LUCC maps complying with the impervious surface maps in 2002, 2007 and 2013, the regions inside the outer ring were found to be basically urbanized, and only a small amount of farmland existed in Pudong new district.Urban areas gradually replace the original vegetation and soil land.The built-up area of the city center has extremely high coverage and has been under Through LUCC maps complying with the impervious surface maps in 2002, 2007 and 2013, the regions inside the outer ring were found to be basically urbanized, and only a small amount of farmland existed in Pudong new district.Urban areas gradually replace the original vegetation and soil land.The built-up area of the city center has extremely high coverage and has been under intensive development for a long time.
Furthermore, the average LSTs of four types of land use in 2002, 2007 and 2013 were marked in Figure 12.The difference between the average temperature of urban and vegetation areas in 2007 exhibited a maximum value among the three phases of images, which was exactly in accordance with the differences of their areas.Thus, a general result could be drawn that the missing vegetation coverage that turned into ISA led to the intensifying of SUHI, which also verified the impact of LUCC on LST.

Surface Urban Heat Island Intensity Distribution
According to the land surface temperature retrieval result, the temperatures are divided into seven intervals for fair analysis considering that different days in divergent months have different land surface temperatures.The distribution maps are shown in Figure 10.Their main characteristic is that the high temperature zone expands from the central urban area to the outskirts.Through the comparison of SUHI intensity maps, areas of heat islands are shown to increase, mainly in the city's center and the centers of districts from 2002 to 2007.However, the heat island effect shows signs of abating from 2007 to 2013, and the expansion of the medium temperature zone and sub-high temperature zone generally focuses on Pudong, Nanhui and Fengxian.The area of extreme high temperature zones and high temperature zones is decreasing, and their distribution becomes scattered, not solely focusing on the city's center and centers of districts.Low temperature zones are located in the outskirts of Shanghai, mainly in Fengxian, Jinshan and parts of Pudong and Nanhui.
To test the phenomenon that the heat island intensity first increased and then decreased, the heat island intensity of the city center and nine other districts was extracted in this study.Here is the calculation formula of the SUHI intensity: T Avg(Downtown) is the average temperature of the downtown area.T Avg(Outskirts) is the average temperature of outskirts, and outskirts are defined as regions outside the outer ring.S d is the standard deviation of temperatures.
Table 12 shows the average LSTs and SUHI intensities of the downtown area and the remaining districts in Shanghai.Obviously, at the imaging time, the average LSTs and SUHI intensities both declined from downtown to outskirts.It is known that LSTs vary in different seasons or even at different times on the same day, and winter nights exhibit the most intensive SUHI in Shanghai due to the specific location and structure [52].Therefore, the result of winter is chosen to distinctly analyze the relationship of LST and LUCC. Figure 13 demonstrates that the heat island intensity gradually decreases from the city center to the suburbs, and reaches the maximum positive value in the city center, while the heat island intensity in the outer suburbs is negative.From the trend of the entire curve, it is easy to see that the heat island effect is the strongest in 2007, while in 2013, it dropped slightly.Both the urban construction and increasing population and human activities contribute to the intensity of the heat island effect.Therefore, the SUHI effect was particularly evident in 2009.Nevertheless, with the coming Expo in 2010, Shanghai not only carried out more rational planning for urban construction but also made a breakthrough in the ecological environmental construction.This progress was maintained after the Expo in 2010.Therefore, the SUHI effect was not that intensive in 2013.

Correlation Analysis of SUHI and IS, LUCC
SUHI is a local climate characteristic affected by global climate change and is determined by human factors and local geographical conditions.The land surface temperature is closely related to the land cover which is an important factor leading to the SUHI.In urban built-up areas, land cover types are mainly impervious surface, vegetation and water.These three ecological factors, meteorological factors, artificial heat have a comprehensive impact on the SUHI effect.Therefore, this section will demonstrate the correlation between LST and IS, NDVI, MNDWI in 2013 with the help of MATLAB software.
A total of 500 sample points were randomly generated by means of MATLAB.The impervious surface rates, NDVI values, MNDWI values and corresponding land surface temperature were used to show the relationships (Figures 14-16).[52].Both the urban construction and increasing population and human activities contribute to the intensity of the heat island effect.Therefore, the SUHI effect was particularly evident in 2009.Nevertheless, with the coming Expo in 2010, Shanghai not only carried out more rational planning for urban construction but also made a breakthrough in the ecological environmental construction.This progress was maintained after the Expo in 2010.Therefore, the SUHI effect was not that intensive in 2013.

Correlation Analysis of SUHI and IS, LUCC
SUHI is a local climate characteristic affected by global climate change and is determined by human factors and local geographical conditions.The land surface temperature is closely related to the land cover which is an important factor leading to the SUHI.In urban built-up areas, land cover types are mainly impervious surface, vegetation and water.These three ecological factors, meteorological factors, artificial heat have a comprehensive impact on the SUHI effect.Therefore, this section will demonstrate the correlation between LST and IS, NDVI, MNDWI in 2013 with the help of MATLAB software.
A total of 500 sample points were randomly generated by means of MATLAB.The impervious surface rates, NDVI values, MNDWI values and corresponding land surface temperature were used to show the relationships (Figures 14-16).meteorological factors, artificial heat have a comprehensive impact on the SUHI effect.Therefore, this section will demonstrate the correlation between LST and IS, NDVI, MNDWI in 2013 with the help of MATLAB software.
A total of 500 sample points were randomly generated by means of MATLAB.The impervious surface rates, NDVI values, MNDWI values and corresponding land surface temperature were used to show the relationships (Figures 14-16).In the correlation analysis of IS and LST, both linear fitting and quadratic fitting gave a high goodness of fit.The linear fitting model passed the F test and the t test under the 95% confidence interval.The regression coefficient is positive, which means that the impervious rate has a linear positive correlation with land surface temperature.That is, land surface temperature increases with the addition of the impervious surface rate.
As for the correlation analysis of NDVI and LST, an apparent linear or exponential correlation does not exist, while a negative trend does.In addition, where rocks and brick exist (NDVI = 0) has higher land surface temperature, while vegetation-covered land (NDVI > 0.7) and water (NDVI around 0.2) have lower land surface temperatures, which become the cold island in Shanghai.In the correlation analysis of IS and LST, both linear fitting and quadratic fitting gave a high goodness of fit.The linear fitting model passed the F test and the t test under the 95% confidence interval.The regression coefficient is positive, which means that the impervious rate has a linear positive correlation with land surface temperature.That is, land surface temperature increases with the addition of the impervious surface rate.
As for the correlation analysis of NDVI and LST, an apparent linear or exponential correlation does not exist, while a negative trend does.In addition, where rocks and brick exist (NDVI = 0) has higher land surface temperature, while vegetation-covered land (NDVI > 0.7) and water (NDVI around 0.2) have lower land surface temperatures, which become the cold island in Shanghai.As for the correlation analysis of MNDWI and LST, the value of LST generally drops in accordance with the increasing value of MNDWI.When the MNDWI value is greater than 0.6, it has lower land surface temperature.Therefore, MNDWI and LST also have a negative correlation, which is not that obvious.The temperatures are distinct depending on regional conditions.
Generally speaking, impervious surface has a positive correlation with land surface temperature and is the strongest between the three factors.Vegetation and water have negative correlations with land surface temperature, and vegetation is more relevant than water.The linear fitting model passed the F test and the t test under the 95% confidence interval.The regression coefficient is positive, which means that the impervious rate has a linear positive correlation with land surface temperature.That is, land surface temperature increases with the addition of the impervious surface rate.
As for the correlation analysis of NDVI and LST, an apparent linear or exponential correlation does not exist, while a negative trend does.In addition, where rocks and brick exist (NDVI = 0) has higher land surface temperature, while vegetation-covered land (NDVI > 0.7) and water (NDVI around 0.2) have lower land surface temperatures, which become the cold island in Shanghai.
As for the correlation analysis of MNDWI and LST, the value of LST generally drops in accordance with the increasing value of MNDWI.When the MNDWI value is greater than 0.6, it has lower land surface temperature.Therefore, MNDWI and LST also have a negative correlation, which is not that obvious.The temperatures are distinct depending on regional conditions.
Generally speaking, impervious surface has a positive correlation with land surface temperature and is the strongest between the three factors.Vegetation and water have negative correlations with land surface temperature, and vegetation is more relevant than water.

Discussion
With the development of urbanization especially in developing countries, buildings and structures have been continuously built over decades.The resulting increase of impervious surface area and decline of water and vegetation have gradually affected land surface temperature and contributed to the SUHI effect.
Studies demonstrate that IS has a warming effect which is a key indicator in urban planning [28][29][30], consistent with the conclusion of this study.Xu stated that LST and IS have exponential relationships [28], while Cao illustrated their linear positive correlation [30].Based on the conclusion drawn from Figure 14, both linear fitting and quadratic fitting methods showed a high goodness of fit.This means that an exponential correlation exists due to variations of imaging time and actual on-site conditions.
In comparison, vegetation and water are able to cool the built-up area.Yuan and Zhou both reported that NDVI has a negative linear correlation with LST [26,27], while Yuan argued that it was a weak correlation, similar to our result in this study.Cao also demonstrated that MNDWI has a strong negative linear correlation with LST [30], but our study shows a weaker dependency than IS.
Nonetheless, slight differences among similar studies are possible.As mentioned above, the images used in this paper were from winter nights.Since the underlying surface of Shanghai transforms from season to season, the change of surface structure led to the variation of LSTs.Thus, the specific exponential relationships or their power might be a little different although they share the same trends between LUCC and LST.

Conclusions
In this study, we analyzed the general trend of impervious surface growth from 2002 to 2007 and its minor increase from 2007 to 2013.The greatest increase in the impervious surface was from farmland, causing a loss of vegetation coverage and leading, to a certain degree, to local SUHI effects.Simulations also demonstrate that the relationship between LST and IS is a positive correlation, while the relationships with NDVI and MNDVI are negative.This result means that the construction of impervious surface strengthens the SUHI effect, while vegetation and water bring about some degree of relief.
Conversely, the driving force of land-use change varies, including the influence of policy, economy, population distribution, and the effect of climate or topography.In short, with the rapid development of urbanization, the policy implications for land use and the control of impervious surface should receive more attention from the urban planning and design used in Shanghai in the future.
The remote sensing data used in this study are images from winter for determining the influence of vegetation shades, which leads to the seasonal influence on the SUHI effect being ignored.Therefore, analysis of different seasons is taken into consideration in the future study.This study focuses on the macroscale problem, which is the speculated relationship between SUHI and IS, water and vegetation.In the near future, it is still necessary to include multiple factors such as ventilation and microclimate for further study in local areas of Shanghai.Moreover, multi-date meteorological data should be introduced to determine how synoptic conditions synergistically affect SUHI.

Figure 1 .
Figure 1.Flow chart of data processing.

Figure 1 .
Figure 1.Flow chart of data processing.

Figure 2 .
Figure 2. The location of Shanghai, China.

Figure 2 .
Figure 2. The location of Shanghai, China.

Figure 7 .
Figure 7. (a) IS extraction result of Shanghai in 2002; (b) IS extraction result of Shanghai in 2007; (c) IS extraction result of Shanghai in 2013.

Figure 7 .
Figure 7. (a) IS extraction result of Shanghai in 2002; (b) IS extraction result of Shanghai in 2007; (c) IS extraction result of Shanghai in 2013.

Figure 7 .
Figure 7. (a) IS extraction result of Shanghai in 2002; (b) IS extraction result of Shanghai in 2007; (c) IS extraction result of Shanghai in 2013.
shows the land use classification results of 2002, 2007 and 2013.The red pixels represent urban areas.The green pixels represent vegetation.The blue pixels represent water.The light-yellow pixels represent other types, mainly bare land.Sustainability 2017, 9, 1538 10 of 22

Figure 8 .
Figure 8. Land use of Shanghai in 2002; (b) Land use of Shanghai in 2007; (c) Land use of Shanghai in 2013.

Figure 8 .
Figure 8. Land use of Shanghai in 2002; (b) Land use of Shanghai in 2007; (c) Land use of Shanghai in 2013.

4. 4 .
Figure9shows that land surface temperature retrieval is a complicated process.

Figure 9 .Figure 10 .
Figure 9. Flow chart of the procedures of land surface temperature retrieval.

Figure 9 .
Figure 9. Flow chart of the procedures of land surface temperature retrieval.
14)/alog (480.89/b1+ 1) − 273.15.The first formula is applicative for images of 2002 and 2007.The later ones are used for 2013.Then, the land surface temperature is obtained.4.4.2.Land Surface Temperature Intervals Since the images used in the study are from different years and different months, normalization of land surface temperature is applied for fair comparison.In addition, the temperatures are divided into seven intervals to display the difference of temperature in different regions of Shanghai in 2002, 2007 and 2013.They are low temperature, sub-low temperature, sub-medium temperature, medium temperature, sub-high temperature, high temperature and extreme high temperature.The division rules are shown in Table 6. Figure 10 is the surface urban heat island effect distribution of 2002, 2007 and 2013.

Figure 10 .
Figure 10.Distribution of urban heat island effect of 2002; (b) Distribution of urban heat island effect of 2007; (c) Distribution of urban heat island effect of 2013.Figure 10.Distribution of urban heat island effect of 2002; (b) Distribution of urban heat island effect of 2007; (c) Distribution of urban heat island effect of 2013.

Figure 10 .
Figure 10.Distribution of urban heat island effect of 2002; (b) Distribution of urban heat island effect of 2007; (c) Distribution of urban heat island effect of 2013.Figure 10.Distribution of urban heat island effect of 2002; (b) Distribution of urban heat island effect of 2007; (c) Distribution of urban heat island effect of 2013.

Figure 11 .
Figure 11.Distribution of urban development density.

Figure 11 .
Figure 11.Distribution of urban development density.
reflect an increase in terms of urban and water areas from 2002 to 2007, while the figure for urban area declines from 2007 to 2013.The trend also complies with the result of impervious surface extraction.From the proportion diagram (Figure 12), the trend of the greater decrease in vegetation area is quite obvious, as it was from 2002 to 2007.The change in water stays in acceptable ranges.

Figure 13 .
Figure 13.Changing of heat island intensity in different districts.

Table 2 .
Experimental data collection.

Table 3 .
Dynamic change of the impervious surface of Shanghai.

Table 3 .
Dynamic change of the impervious surface of Shanghai.

Table 3 .
Dynamic change of the impervious surface of Shanghai.

Table 4 .
Accuracy assessment of impervious surface extraction of Shanghai.

Table 4 .
Accuracy assessment of impervious surface extraction of Shanghai.

Table 5 .
Land use classification accuracy.

Table 5 .
Land use classification accuracy.

Table 6 .
The range of 7 land surface temperature intervals.> T a + 2S d High temperature T a + S d T S ≤ T a + 2S d Sub-high temperature T a + S d /2T S ≤ T a + S d Medium temperature T a − S d /2T S ≤ T a + S d /2 Sub-medium temperature T a − S d T S ≤ T a − S d /2 Sub-low temperature T a − 2S d T S ≤ T a − S d Low temperature T S < T a + 2S d T S represents land surface temperature.T a is the average land surface temperature.S d is standard deviation.

Table 8 .
Regional proportion of impervious surface.

Table 10 .
Growth rate of impervious surface area.

Table 12 .
Urban heat island intensive distribution in different dstricts ( • C).