Spatiotemporal Analysis of Land Use / Cover Patterns and Their Relationship with Land Surface Temperature in Nanjing, China

: Rapid urbanization is one of the most concerning issues in the 21st century because of its signiﬁcant impacts on various ﬁelds, including agriculture, forestry, ecology, and climate. The urban heat island (UHI) phenomenon, highly related to the rapid urbanization, has attracted considerable attention from both academic scholars and governmental policymakers because of its direct inﬂuence on citizens’ daily life. Land surface temperature (LST) is a widely used indicator to assess the intensity of UHI signiﬁcantly a ﬀ ected by the local land use / cover (LULC). In this study, we used the Landsat time-series data to derive the LULC composition and LST distribution maps of Nanjing in 2000, 2014, and 2018. A correlation analysis was carried out to check the relationship between LST and the density of each class of LULC. We found out that cropland and forest in Nanjing are helping to cool the city with di ﬀ erent degrees of cooling e ﬀ ects depending on the location and LULC composition. Then, a Cellar Automata (CA)-Markov model was applied to predict the LULC conditions of Nanjing in 2030 and 2050. Based on the simulated LULC maps and the relationship between LST and LULC, we delineated high- and moderate-LST related risk areas in the city of Nanjing. Our ﬁndings are valuable for the local government to reorganize the future development zones in a way to control the urban climate environment and to keep a healthy social life within the city.


Introduction
Following the rapid development of urbanization and industrialization, the quantity, structure, and degree of land use/cover (LULC) have been significantly changed over the past years, especially in developing countries [1,2]. Meanwhile, the urban heat island (UHI) phenomenon has become severely causing excessive consumption of energy, an increase of air pollutants, and deterioration of the environment [3]. Therefore, by determining the relationship between land surface temperature (LST) and LULC distribution, we can evaluate the impact of urbanization on LST, which is crucial for sustainable urban development.
The key physical drivers of the UHI are fivefold as summarized by Coffel et al. [4]: (1) the increase of flat and cemented surface at the expense of frictional surfaces in urban areas, which leads to the a result of the Chinese economic reform. According to the National Bureau of Statistics, the total population of Nanjing increased from 6.13 million in 2000 to 8.34 million in 2018, with an annual increase rate of 12.28% [21]. In 2016, the area of the built-up in Nanjing expanded to 773.79 km 2 , pushing the city to rank as the ninth-largest among all Chinese cities [21].
Nanjing has a subtropical monsoon climate with an average annual temperature of 16.1 • C and annual precipitation of 1106.8 mm in 2017 [22]. The Yangtse River is crossing the area from southwest to east. In this study, a range of 31 • 34 -32 • 33 N and 118 • 15 -119 • 20 E, which is the central part of Nanjing city is decided as the study area ( Figure 1). The topography of the study area ranges between an elevation of -92 m and 441 m with hilly and mountainous areas surrounding the city from almost all directions.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 17 change process. From 1978, Nanjing, like other Chinese cities, has experienced breakneck urbanization as a result of the Chinese economic reform. According to the National Bureau of Statistics, the total population of Nanjing increased from 6.13 million in 2000 to 8.34 million in 2018, with an annual increase rate of 12.28% [21]. In 2016, the area of the built-up in Nanjing expanded to 773.79 km 2 , pushing the city to rank as the ninth-largest among all Chinese cities [21]. Nanjing has a subtropical monsoon climate with an average annual temperature of 16.1°C and annual precipitation of 1106.8 mm in 2017 [22]. The Yangtse River is crossing the area from southwest to east. In this study, a range of 31°34'-32°33'N and 118°15'-119°20'E, which is the central part of Nanjing city is decided as the study area ( Figure 1). The topography of the study area ranges between an elevation of -92 m and 441 m with hilly and mountainous areas surrounding the city from almost all directions.

Workflow
The workflow of this study is shown in Figure 2. The final objectives are to simulate the future LULC development in 2030 and 2050 and to assess the risk of high LST distribution in order to support urban planning initiatives and sustainable development in the study area. To achieve these goals, several steps have been conducted. First, we applied the maximum likelihood classification to extract LULC maps from Landsat images. The principle of maximum likelihood supervised classification is based on training samples and is used to classify pixels into suitable LULC categories [23,24]. The selection of training samples was based on our previous knowledge of the geographical settings of the study area. Second, we conducted an urban-gradient analysis to monitor the urbanization trends in Nanjing. Third, we retrieved LST maps and conducted the regression analysis with LULC maps to examine the existing relationship. Lastly, we simulated the LULC in 2030 and 2050 using a hybrid model of Markov and CA models with the aim of assessing the LST related risk areas as a final step [25,26].

Workflow
The workflow of this study is shown in Figure 2. The final objectives are to simulate the future LULC development in 2030 and 2050 and to assess the risk of high LST distribution in order to support urban planning initiatives and sustainable development in the study area. To achieve these goals, several steps have been conducted. First, we applied the maximum likelihood classification to extract LULC maps from Landsat images. The principle of maximum likelihood supervised classification is based on training samples and is used to classify pixels into suitable LULC categories [23,24]. The selection of training samples was based on our previous knowledge of the geographical settings of the study area. Second, we conducted an urban-gradient analysis to monitor the urbanization trends in Nanjing. Third, we retrieved LST maps and conducted the regression analysis with LULC maps to examine the existing relationship. Lastly, we simulated the LULC in 2030 and 2050 using a hybrid model of Markov and CA models with the aim of assessing the LST related risk areas as a final step [25,26].

Data Processing
LULC and LST maps were extracted from three Landsat images captured on September 16, 2000 (Landsat 7 ETM+), November 18, 2014 (Landsat 8 OLI/TIRS) and October 12, 2018 (Landsat 8 OLI/TIRS) ( Table 1). All the images were obtained from the United States Geological Survey (USGS) [27]. To avoid cloud-coverage issues, images were acquired during the autumn season with a cloudcoverage of less than 10%. Additionally, it is worth mentioning that seasonal variations also influence NDVI and LST values. Usually, from March to November, it is easy to distinguish the differences in LST intensities between different LULC classes [28]. Consequently, in this study, all the Landsat images were collected between these two months to exclude any possible seasonal variability effects. Other auxiliary data, mainly the modeling variables for driving LULC changes, were collected as well, including the central business district (CBD), main roads, water areas, protected areas, Digital Elevation Model (DEM) and slope (See Figure A1 in Appendix A). The CBD is located in the central region of the study area, which has potential effects on economic development and urbanization. The main roads are always recognized as one of the essential factors influencing urban construction. The permanent water bodies and protected areas were selected for their quasi-stable status during the LULC change process as these areas are difficult to be converted into other LULC types. The DEM and slope data were considered due to the varied terrain topography of the study area, characterized by the disparity in elevation. The DEM was obtained from the USGS which was used to extract the slope map.

Data Processing
LULC and LST maps were extracted from three Landsat images captured on September 16, 2000 (Landsat 7 ETM+), November 18, 2014 (Landsat 8 OLI/TIRS) and October 12, 2018 (Landsat 8 OLI/TIRS) ( Table 1). All the images were obtained from the United States Geological Survey (USGS) [27]. To avoid cloud-coverage issues, images were acquired during the autumn season with a cloud-coverage of less than 10%. Additionally, it is worth mentioning that seasonal variations also influence NDVI and LST values. Usually, from March to November, it is easy to distinguish the differences in LST intensities between different LULC classes [28]. Consequently, in this study, all the Landsat images were collected between these two months to exclude any possible seasonal variability effects. Other auxiliary data, mainly the modeling variables for driving LULC changes, were collected as well, including the central business district (CBD), main roads, water areas, protected areas, Digital Elevation Model (DEM) and slope (See Figure A1 in Appendix A). The CBD is located in the central region of the study area, which has potential effects on economic development and urbanization. The main roads are always recognized as one of the essential factors influencing urban construction. The permanent water bodies and protected areas were selected for their quasi-stable status during the LULC change process as these areas are difficult to be converted into other LULC types. The DEM and slope data were considered due to the varied terrain topography of the study area, characterized by the disparity in elevation. The DEM was obtained from the USGS which was used to extract the slope map.

LST Retrieval
LST is a crucial factor that can reflect the surface temperature among different LULC types. For retrieving LST maps from Landsat images, several steps need to be conducted. Initially, the Digital Numbers (DNs) of thermal bands (band 6 in Landsat TM and ETM+, and band 10 in Landsat OLI/TIRS) should be converted into radiance values [29][30][31]. These values should be collected for calculating the at-satellite brightness temperature [3,32]. All images have been undergone atmospheric and emissivity correction procedures prior to the analysis. In this study, the preprocessed band 6 of Landsat 7 in 2000 and band 10 of Landsat 8 in 2014 and 2018 were used for expressing at-satellite brightness temperature in Kelvin. The mathematic equation of land surface emissivity (ε) can be expressed as follows [33]: where P v is the proportion of vegetation, which can be calculated using the following equation: where NDVI refers to the Normalized Difference Vegetation Index, which reflects the vegetation density of land. It can be derived from the red and near-infrared bands of Landsat images [34] using the following formula: Lastly, LST values can be calculated using the following equation [35]: where TB is band 6 in Landsat ETM+, and band 10 in Landsat OLI/TIRS; λ is the wavelength of emitted radiance (11.5 µm for Landsat ETM+ band 6 and 10.8 µm for Landsat OLI/TIRS band 10) [36];

LULC Classification
LULC classification is a significant challenge when it comes to obtaining complex LULC distribution information from Landsat data [18,37]. Several classification methods have been employed in different studies depending on the target area and the spatial resolution of the satellite images [38][39][40]. From a performance perspective, the maximum likelihood supervised classification method provides an accurate outcome [23]. Therefore, we selected it as the primary classification method in this study. Based on the statistic features of remote sensing data, the maximum likelihood supervised classification method can calculate the mean value and variance to establish a function. According to this function, each pixel can verify their belonging category by plugging it into this equation.
The maximum likelihood supervised classification method was employed to classify LULC maps using ArcGIS 10.4 software. Based on the purpose of this study, four LULC categories were considered, namely built-up, cropland, green, and water. The built-up category includes impervious surfaces such as buildings, roads, airports, parking lots, and sidewalks. The cropland category consists of agricultural lands without tree canopy. The category of green comprises forest, woodland, shrubland, grassland, and green spaces with the tree canopy. The water category contains all water bodies, such as rivers, lakes, ponds, reservoirs, and wetlands. The number of training samples selected on Landsat images for classification was 368, 395, and 386 for 2000, 2014, and 2018, respectively. The range of pixels' size is between 80 and 140 for all the training samples. Finally, an accuracy assessment was conducted to confirm the exactitude of the LULC maps.

Urban-Rural Gradient Analysis
A prefecture-level city always consists of several spaces. Each space serves certain or multiple urban purposes. Since these sectors have different characteristics, their LULC composition differs. Keeping in mind that LULC composition strongly influences the future development potential in terms of sustainability and economy, it is essential to analyze the LULC composition for each zone to perceive the future trends of LULC. Urban-rural gradient analysis is a method to analyze the spatial distribution of different LULC classes by predefined distance ranges from the city center all the way to the nearby rural areas [1,39]. The method can help to understand the urbanization structure and discover further opportunities for potential development in the city. In this study, 50 ring buffers were created from the city center to the rural areas in Nanjing with a distance interval of 1 km.

Markov Model
Markov model is a stochastic model that can trace the changing process from one state to another. The future state is only influenced by the current state, and each state should keep the same time step [41]. Markov model can be applied to reflect the change of LULC process [25,42]: (1) all the LULC types can change from each other under the historical influence; (2) LULC change is a complicated process which is influenced by several parameters and cannot be described by simple functions; and (3) the transformation of LULC is stable during the time period. Based on its characteristics, the Markov model is an essential method for simulating future LULC changes. Therefore, it was considered for the LULC simulation. The Markov process can be written as follows: where P (n) is the transition state probability in time period n, and P (n−1) is the current state probability in time period n − 1.
During the LULC transition process in the Markov model, the condition of LULC in the next time period is only influenced by the current state and the transition rule subordinated to the transition probability matrix [43]. The original transition probability matrices can be calculated from LULC maps, and the LULC category in a one-time period should be defined before employing the Markov model. This can be expressed as follows: where P ij refers to the transformation probability matrix of LULC changes occurred from ith category into jth within one time period, and n is the land use category in the study area. P should meet the following condition: By considering that LULC maps do not keep the same time period, the "Hazard rate approach" was applied during the Markov process. The per unit time transition matrix should be calculated using the following equation: where H ij is the hazard rate in one year. In this study, the transition probability matrix from 2014 to 2018 was calculated using the Markov model at first. Next, the hazard rate approach was combined with a Markov model to extract a one-year transition probability matrix from 2014 to 2018. Finally, the employed model was set to calculate the one-year transition probability matrix 12 and 32 times to simulate the LULC map in 2030 and 2050, respectively.

Cellular Automata-Markov Model
The CA-Markov model is a hybrid model combining Markov and cellular automata models. Because of its ability to allocate spatial information of landscape, it has been widely used in land cover modeling and landscape prediction [1,44]. It has already gained an excellent reputation and has a prospective future in several fields, for instance: computer science, mathematics, and geographical science [45,46]. The CA-Markov model can be applied to spatiotemporally simulate future LULC distribution. The CA model is established of spatial cells, which are all independent and interactive with each other. The transition process of the CA model also can reflect the interaction among spatial cells. Because of these characteristics, it has also been applied for simulating the LULC process [47,48]. In the CA model, all cells should obey the specific transition rules. The transition rules employed a 3 × 3 neighborhood to affect the central cell condition in the future. Using the CA model to simulate the future LULC, the property of LULC types should be considered, for instance, the built-up area cannot transfer into the water in the near future [25].
LULC change is a complex process affected by historic trends, environmental effects, transition rules, and land use policy. To comprehensively simulate the LULC process, the CA model and the Markov model were combined to provide spatiotemporal changes. The following equation illustrates the cell changing process from time t to time t + 1 [49]: where and NI (i,j)t is the interaction by surrounding cells under the transition rules.  (Table 2). Cropland areas are mainly located in rural areas. As a result of the expansion of urban areas to meet the growing urbanization demand, cropland areas decreased from 713,951.55 ha in 2000 to 556,336.68 ha in 2018 (Table 2). On the other hand, no changes have been observed for green areas during the study period, and this is due to the location of these areas in the mountains. Similarly, no huge changes were witnessed in the total area of water bodies consisted mainly of the Yangtse River.

Cellular Automata-Markov Model
The CA-Markov model is a hybrid model combining Markov and cellular automata models. Because of its ability to allocate spatial information of landscape, it has been widely used in land cover modeling and landscape prediction [1,44]. It has already gained an excellent reputation and has a prospective future in several fields, for instance: computer science, mathematics, and geographical science [45,46]. The CA-Markov model can be applied to spatiotemporally simulate future LULC distribution. The CA model is established of spatial cells, which are all independent and interactive with each other. The transition process of the CA model also can reflect the interaction among spatial cells. Because of these characteristics, it has also been applied for simulating the LULC process [47,48]. In the CA model, all cells should obey the specific transition rules. The transition rules employed a 3  3 neighborhood to affect the central cell condition in the future. Using the CA model to simulate the future LULC, the property of LULC types should be considered, for instance, the built-up area cannot transfer into the water in the near future [25].
LULC change is a complex process affected by historic trends, environmental effects, transition rules, and land use policy. To comprehensively simulate the LULC process, the CA model and the Markov model were combined to provide spatiotemporal changes. The following equation illustrates the cell changing process from time t to time t + 1 [49]: ( , ) is the transition probability of cell (i, j) at time t, and ( , ) is the interaction by surrounding cells under the transition rules.  Table 2). Cropland areas are mainly located in rural areas. As a result of the expansion of urban areas to meet the growing urbanization demand, cropland areas decreased from 713,951.55 ha in 2000 to 556,336.68 ha in 2018 (Table 2). On the other hand, no changes have been observed for green areas during the study period, and this is due to the location of these areas in the mountains. Similarly, no huge changes were witnessed in the total area of water bodies consisted mainly of the Yangtse River.   The results of the urban-rural gradient analysis are shown in Figure 4. In general, the percentage of built-up areas decreased from the city center to rural areas.  The results of the urban-rural gradient analysis are shown in Figure 4. In general, the percentage of built-up areas decreased from the city center to rural areas.

Relationship between LST and LULC
The spatial pattern of NDVI of Nanjing for 2000, 2014, and 2018 are presented in Figure 5

Relationship between LST and LULC
The spatial pattern of NDVI of Nanjing for 2000, 2014, and 2018 are presented in Figure 5.       Figure 7, the locations of water dots are easy to distinguish from other LULC types as they are scattered with negative values of NDVI and LST values inferior to 0.5. The built-up dots are concentrated with high LST values. Contrarily, green dots have obviously high NDVI values with low LST values. To clearly understand the correlation between LULC and LST, we summarized the statistical values of the sample points in Table 3. We observed a continuous trend that built-up has the highest mean LST, followed by cropland, green, and water in the three years. In addition, the standard deviation ranged between 1.81 • C (green) and 2.  Table 3. We observed a continuous trend that built-up has the highest mean LST, followed by cropland, green, and water in the three years. In addition, the standard deviation ranged between 1.81°C (green) and 2.

Future LULC and LST in 2030 and 2050
The simulation results of LULC change in 2030 and 2050 are illustrated in Figure 8. Comparing the LULC maps of 2000, 2014, and 2018 (Figure 3), it can be observed that built-up areas are expected to continually expand and even occupy the nearby croplands surrounding them. The area

Future LULC and LST in 2030 and 2050
The simulation results of LULC change in 2030 and 2050 are illustrated in Figure 8. Comparing the LULC maps of 2000, 2014, and 2018 (Figure 3), it can be observed that built-up areas are expected to continually expand and even occupy the nearby croplands surrounding them. The area percentages of each of the four LULC classes from 2000 to 2050 are shown in Figure 9. Same with other developing countries, built-up areas are expanding, to meet the demand of urbanization, at the expense of cropland that is consistently decreasing. Most of the green areas are expected to remain intact from 2018 to 2050. During the LULC change process, water is expected to preserve its share in the total area. Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 17 percentages of each of the four LULC classes from 2000 to 2050 are shown in Figure 9. Same with other developing countries, built-up areas are expanding, to meet the demand of urbanization, at the expense of cropland that is consistently decreasing. Most of the green areas are expected to remain intact from 2018 to 2050. During the LULC change process, water is expected to preserve its share in the total area.  Based on the relationship between LULC and LST, the future distribution of LST has been estimated ( Figure 10). It can be observed that areas with high LST risk are expected to be distributed mainly in the central area in 2030 resulted from the forecasted built-up distribution. Whereas, moderate LST risk areas are anticipated to be located in several scattered settlements. The overall LST risk areas are expected to increase from 2030 to 2050. It is projected that the high LST risk areas will be expanded by 2050, and some new moderate LST risk areas are expected to be generated. The LST differences between the central area and the surrounding areas are forecasted to become higher, which may cause environmental problems that can have bad effects on human health. These simulation results may provide insights to the local government and city planners to take corrective measures for wiser sustainable urban development policies.  Figure 9. Same with other developing countries, built-up areas are expanding, to meet the demand of urbanization, at the expense of cropland that is consistently decreasing. Most of the green areas are expected to remain intact from 2018 to 2050. During the LULC change process, water is expected to preserve its share in the total area.  Based on the relationship between LULC and LST, the future distribution of LST has been estimated ( Figure 10). It can be observed that areas with high LST risk are expected to be distributed mainly in the central area in 2030 resulted from the forecasted built-up distribution. Whereas, moderate LST risk areas are anticipated to be located in several scattered settlements. The overall LST risk areas are expected to increase from 2030 to 2050. It is projected that the high LST risk areas will be expanded by 2050, and some new moderate LST risk areas are expected to be generated. The LST differences between the central area and the surrounding areas are forecasted to become higher, which may cause environmental problems that can have bad effects on human health. These simulation results may provide insights to the local government and city planners to take corrective measures for wiser sustainable urban development policies. Based on the relationship between LULC and LST, the future distribution of LST has been estimated ( Figure 10). It can be observed that areas with high LST risk are expected to be distributed mainly in the central area in 2030 resulted from the forecasted built-up distribution. Whereas, moderate LST risk areas are anticipated to be located in several scattered settlements. The overall LST risk areas are expected to increase from 2030 to 2050. It is projected that the high LST risk areas will be expanded by 2050, and some new moderate LST risk areas are expected to be generated. The LST differences between the central area and the surrounding areas are forecasted to become higher, which may cause environmental problems that can have bad effects on human health. These simulation results may provide insights to the local government and city planners to take corrective measures for wiser sustainable urban development policies.

Discussion
Existing studies have pointed out that the increase of LST is related to urban growth, especially in developing countries [50]. The agglomeration and expansion of urban areas have a significant effect on the increase of impervious surfaces [51,52]. Moreover, previous studies have found that green space can maintain soil moisture to reduce LST [53]. In this study, we analyzed the relationship between LULC and LST by using remote sensing data. Our results, based on quantitative analysis and spatial distribution, showed that built-up, cropland, green, and water have different effects on LST. In addition, an interesting result has been found, which is the fact that the central LST of builtup has been reduced considerably from 2000 to 2018 (Figure 7). This might be attributed to the existence of natural obstacles in the form of water bodies and green spaces (high altitude) located in the central area. Meanwhile, built-up areas have expanded beyond these natural obstacles in almost all directions, creating what can be called a mixed LULC area ( Figure 3).
In recent years, the relationship between LST and LULC has received more and more attention from researchers [14,54]. However, few of them have reached similar results because of regional differences and geographical settings of the targeted areas. In this research, Nanjing was selected as a case study reasoned by the fact that it is one of three furnace cities in China. Local people are suffering from high temperatures during summertime. Therefore, reducing the LST to moderate levels has become a critical target for improving the living standards and supporting future sustainable development. In this work, by combining remote sensing techniques and geographic information system capacities, we examined the impacts of multiple LULC categories on LST quantitatively and spatially. Furthermore, a simulation of future land use/cover changes in 2030 and 2050 has been conducted. ArcGIS and TerrSet software were used for LULC classification and examining LST changes in each LULC category at different spatial scales. The Landsat data of three different years during autumn season (2000, 2014 and 2018) of the Chinese city of Nanjing were acquired to analyze the LULC changes and their relationship with LST.
The monitoring result shows that during the urbanization and industrialization process, the built-up area had been expanding from the central area of the city to the suburban areas from 2000 to 2018. At the same time, the UHI phenomenon started to be observed in Nanjing. Many studies have shown that green space can help reducing LST, especially in the central area [15,55]. Additionally, mixed land use not only can enrich the living environment but also can decrease LST to support sustainable development [56,57]. Since it is prohibited to change protected green areas into other

Discussion
Existing studies have pointed out that the increase of LST is related to urban growth, especially in developing countries [50]. The agglomeration and expansion of urban areas have a significant effect on the increase of impervious surfaces [51,52]. Moreover, previous studies have found that green space can maintain soil moisture to reduce LST [53]. In this study, we analyzed the relationship between LULC and LST by using remote sensing data. Our results, based on quantitative analysis and spatial distribution, showed that built-up, cropland, green, and water have different effects on LST. In addition, an interesting result has been found, which is the fact that the central LST of built-up has been reduced considerably from 2000 to 2018 (Figure 7). This might be attributed to the existence of natural obstacles in the form of water bodies and green spaces (high altitude) located in the central area. Meanwhile, built-up areas have expanded beyond these natural obstacles in almost all directions, creating what can be called a mixed LULC area ( Figure 3).
In recent years, the relationship between LST and LULC has received more and more attention from researchers [14,54]. However, few of them have reached similar results because of regional differences and geographical settings of the targeted areas. In this research, Nanjing was selected as a case study reasoned by the fact that it is one of three furnace cities in China. Local people are suffering from high temperatures during summertime. Therefore, reducing the LST to moderate levels has become a critical target for improving the living standards and supporting future sustainable development. In this work, by combining remote sensing techniques and geographic information system capacities, we examined the impacts of multiple LULC categories on LST quantitatively and spatially. Furthermore, a simulation of future land use/cover changes in 2030 and 2050 has been conducted. ArcGIS and TerrSet software were used for LULC classification and examining LST changes in each LULC category at different spatial scales. The Landsat data of three different years during autumn season (2000, 2014 and 2018) of the Chinese city of Nanjing were acquired to analyze the LULC changes and their relationship with LST.
The monitoring result shows that during the urbanization and industrialization process, the built-up area had been expanding from the central area of the city to the suburban areas from 2000 to 2018. At the same time, the UHI phenomenon started to be observed in Nanjing. Many studies have shown that green space can help reducing LST, especially in the central area [15,55]. Additionally, mixed land use not only can enrich the living environment but also can decrease LST to support sustainable development [56,57]. Since it is prohibited to change protected green areas into other LULC categories, no significant changes have occurred in these areas during the study period. The LST values in these areas are lower than the ones observed in the central area.
Based on the extracted LST and LULC maps, the characteristics of UHI in Nanjing has been discussed. According to the simulated LULC maps of 2030 and 2050, similarly to other cities in developing countries, urban growth is expected to continue in the city (Figure 8). Additionally, maps of the LST risk area have been extracted and have been analyzed based on the thermal radiation characteristics of each LULC type (Figure 9). The result shows that if the urbanization sustained from 2018 to 2050, the high LST risk areas are expected to expand.
Based on the results of simulating LST risks in Nanjing, we propose the following three suggestions to alleviate the LST increase: (1) to promote the use of renewable energy resources and environmental protection materials: For example, using high reflectance materials on buildings' walls or roofs to reduce absorption from solar radiation energy, and using water permeability materials to construct the urban roads in order to improve transpiration; (2) to build ecological cities that support sustainable development, where urban green areas can significantly decrease LST and improve the living standards; and (3) to properly redistribute the transportation network and to encourage the utilization of public transportation.
The Chinese government has always governed and regulated urbanization. In recent years, it has started to become aware of the risks posed by UHI on its sustainable development strategy. The sponge city, forest city, and climate-sensitive city are examples of the proposed cities of the future. All of them can decrease UHI by establishing built-up mixed green areas [58,59]. In this study, the government regulations related to environment protection has been considered to reduce the LST in the future simulation in 2030 and 2050. Compared with the past period 2000-2018, the impact of the population on urban planning is expected to getting smaller. In contrast, the limited supply of land and economical regulations have an increasing effect on LULC change. Therefore, to examine the current development situation of LULC and simulate the future LULC change under the green area protection regulations, it is possible to verify the feasibility of the current city plan to reduce the potential pressure on the environment in the future. According to the LULC maps in 2030 and 2050 (Figure 8), the built-up is continually increasing by occupying other LULC types. Compared with the period from 2000 to 2018, the urbanization speed has slowed down during the period from 2018 to 2050. This trend is conforming to the existing city planning forms in China.
Establishing a sustainable development in the future is always recognized as the ultimate objective in many LULC studies [60,61]. There are many factors that influence LST to engender the UHI phenomenon, such as psychrometrics, temperature, and intensity of illumination [62]. Compared with other factors, LULC is a relatively constant factor that predominantly influences the causation of UHI. By simulating future land use/cover changes and analyzing the relationship with LST, UHI effects can be effectively reduced. Land use/cover simulation not only can provide several scenarios of what can be occurred in the future but also avoid disparate development between social-economic and environmental protection.
Although increasing green space can relieve LST, there are still some limitations that should be pointed out. There are many extraneous and uncontrollable factors influencing LST such as precipitation, illumination, season, etc. This is the reason why LST is different day by day, even hour by hour. In this study, we only focused on the spatiotemporal correlation land use/cover types with LST. Therefore, in future research, we will assess the impacts of the density of green space and built-up areas, which impact the living environment to accommodate government policy and social problems.

Conclusions
In this study, we monitored the spatio-temporal dynamics of land use/cover and LST of Nanjing from 2000 to 2018. The results demonstrated that during the first 18 years of the 21st century, Nanjing had experienced rapid urbanization with the cost of a tremendous decrease in the agricultural lands.
The expansion of impervious surfaces resulted from the process of rapid urbanization, has caused an increase in urban LST indicating that the UHI phenomenon has become more severe in Nanjing. Through a correlation analysis between LST and landscape composition, we found that forest and cropland had the highest potential in cooling the LST in Nanjing. A significant impact of water bodies on the temperature of the surrounding areas was not recognized. The cooling effects of cropland and forest differed when the surrounding landscape composition differed. Considering the abovementioned findings, it would be more efficient to find a suitable combination of landscape composition that would be more efficient in cooling the LST instead of merely increasing the cropland or forest.
The CA-Markov model was employed to simulate the land use/cover conditions of Nanjing in 2030 and 2050. From the empirical outcomes, we observed a continuous and fast urban expansion accompanied with the loss of cropland and forest. Based on the simulation, the central areas in Nanjing along the Yangtze River would be supposed to face severe UHI impacts in the coming future if no effective policies or regulations are established. The results of the carried analyses are an alarming information for the local government that the landscape of the central areas of Nanjing should be reorganized to cope up with the expected upcoming severe effects of UHI on the residents. Our study proposed a valuable framework to delineate and predict the risk of UHI in the future based on the changing land use/cover conditions. The method developed in this study would be applicable in other megacities to assess the areas with high LST risk.  (c) distance to CBD; (d) distance to green; (e) distance to road; (f) distance to water; (g) slope; (h) protected area.