Dynamic Changes of Local Climate Zones in the Guangdong–Hong Kong–Macao Greater Bay Area and Their Spatio-Temporal Impacts on the Surface Urban Heat Island Effect between 2005 and 2015

Local climate zones (LCZs) emphasize the influence of representative geometric properties and surface cover characteristics on the local climate. In this paper, we propose a multi-temporal LCZ mapping method, which was used to obtain LCZ maps for 2005 and 2015 in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA), and we analyze the effects of LCZ changes in the GBA on land surface temperature (LST) changes. The results reveal that: (1) The accuracy of the LCZ mapping of the GBA for 2005 and 2015 is 85.03% and 85.28%, respectively. (2) The built type category showing the largest increase in area from 2005 to 2015 is LCZ8 (large low-rise), with a 1.01% increase. The changes of the LCZs also vary among the cities due to the different factors, such as the economic development level and local policies. (3) The area showing a warming trend is larger than the area showing a cooling trend in all the cities in the GBA study area. The main reasons for the warming are the increase of built types, the enhancement of human activities, and the heat radiation from surrounding high-temperature areas. (4) The spatial morphology changes of the built type categories are positively correlated with the LST changes, and the morphological changes of the LCZ4 (open high-rise) and LCZ5 (open midrise) built types exert the most significant influence. These findings will provide important insights for urban heat mitigation via rational landscape design in urban planning management.


Introduction
As one of the most significant changes in the development of human society, urbanization is the inevitable result of socio-economic development [1]. The variation of land use/land cover (LULC) and the enhancement of human activities in rapid urbanization [2,3] have caused changes to the urban thermal fields, which have led to the emergence of the urban heat island (UHI) effect [4,5]. The UHI effect refers to the higher temperature that occurs in urban areas compared to that in suburban or rural areas. Therefore, the UHI intensity is defined as the temperature difference between the urban and rural areas. Most current calculations of the UHI intensity are based on the urban-rural temperature differences. Urban theorists claim that the relation between city and country is more accurately described as a continuum, or a dynamic, rather than as a dichotomy [6]. The complexity of urban-rural systems and the difficulty of defining urban-rural boundaries [7,8] have brought great uncertainty to the quantification of the UHI intensity. In addition, only using urban and rural areas to describe the major differences in the urban thermal environment could result in the loss of numerous surface details (e.g., land cover, urban morphology and function) that determine the local climate differences. Hence, Stewart and Oke proposed a detailed surface classification framework for urban thermal field research named the local climatic zone (LCZ) classification scheme. LCZs are regions with uniform surface cover, structure, materials, and human activity, in a horizontal scale of hundreds of meters to several kilometers, including 10 built types and seven land-cover types. The LCZ classification scheme does not include climate factors, but is based on the attributes that affect climates, namely the surface structure (height and spacing of buildings and trees) and surface cover (pervious or impervious). The surface structure affects the local climate through its modification of airflow, atmospheric heat transport, and shortwave and longwave radiation balances, while surface cover modifies the albedo, moisture availability, and heating/cooling potential of the ground [8]. Thus, although the description of the LCZ classification is similar to the land cover or land use classification, the LCZ characteristics are more focused on climate-related factors. Being highly correlated with urban structure and the land-cover types, the LCZ classification scheme emphasizes the influence of the representative geometric properties and surface cover characteristics on the local climate [9]. As consequence, LCZs have been used to explore how urbanization affects the thermal fields of human life, mainly referring to the urban heat island (UHI) effect.
The land surface temperature (LST) regulates the lower-layer air temperature of the urban atmosphere [10]. When obtained from thermal infrared remote sensing imagery, LST can represent the continuous daytime or nighttime temperatures over a larger geographic area than air temperature [11]. Hence, LST has become a key variable in the study of the surface UHI (SUHI) effect during the urbanization process [12][13][14][15][16][17].
To date, the studies of the relationship between LCZs and LST have mainly focused on the correlation analysis of a single time phase [18][19][20][21], including the descriptive analysis of LCZs and LST [22], the comparative analysis of LST in the LCZs of different regions [7], and the correlation analysis of surface features and LST [23]. However, there is no doubt that the single-time analyses ignore the changes in urban structure and urban thermal fields during the urbanization process. Multi-temporal LCZ analysis has the potential to be used to investigate the causes and trends of urban climate changes by studying the long-term dynamics of LCZs and LST, which can compensate for the shortcomings of a single-year analysis. Despite the apparent importance of multi-temporal LCZ analysis for SUHI research, there have been few studies to date of multi-temporal LCZ and LST analysis. Moreover, these studies have only described LCZ changes from the macroscopic point of view [24][25][26][27]; for example, by simply studying the variation of the spatial composition of the LCZs in a single city over a long period of time, or the LST statistics for multiple LCZ categories at different times. To sum up, due to the lack of microscopic and refined quantitative analysis, the studies conducted to date cannot clearly reflect the changes of the LCZs and LST within a city. Consequently, a more fine-grained study of the relationship between LCZs and LST over a long time series is of great importance for an in-depth understanding of the impacts of urbanization on the SUHI effect.
On the other hand, multi-temporal LCZ mapping is the basis of dynamic LCZ analysis. To date, the method of single-temporal mapping has commonly been adopted to obtain multi-temporal LCZ maps. Specifically, the LCZ training samples with different time phases have been selected separately, and then the World Urban Database and Access Portal Tools (WUDAPT) method [27,28] or a supervised classification method [25,29] has been applied to obtain multi-temporal LCZ maps. However, because of the relatively large number of LCZ categories and the large amount of training samples required for the mapping, selecting samples for each separate classification process is a time-consuming and labor-intensive task. Therefore, this kind of sample selection method is difficult to apply for multi-temporal or large-geographic-scale LCZ mapping, especially for time phases for which samples are difficult to obtain. Moreover, this approach it is not conducive to the accurate and consistent classification of time-series images in a more automated manner. In this context, a time-series sample migration strategy is proposed in this paper to reduce the workload of sample collection and facilitate faster and more efficient multi-temporal LCZ classification.
As one of the most important bay areas in the world, and one of the most open and economically dynamic regions in China, the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) stands at an important strategic position in the overall development of China. According to the Outline Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area (https://www.bayarea.gov.hk/filemanager/en/share/pdf/Outline_ Development_Plan.pdf, accessed on 28 January 2021), the aim is for the GBA to be not only a vibrant world-class city cluster and an important part of the "One Belt One Road" initiative, but also a high-quality living area that is pleasant to live in, work in, and travel to, which will eventually become a model for high-quality urban development. However, the urban agglomeration in the GBA, as the main spatial carrier of urbanization, has brought about a series of ecological and environmental problems while promoting rapid socio-economic and political development. One of these problems is the UHI effect, which is caused by the dramatic increase of artificial surfaces and the enhancement of human activities. Accordingly, how to mitigate the UHI effect and promote the sustainable development of the GBA region has become an important goal of regional development planning.
Recently, several studies related to the impacts of urbanization on ecological quality have been conducted in the GBA [30][31][32][33], for which the SUHI effect has been one of the important concerns. Some of these studies have analyzed urban outdoor comfort from the perspective of resident perception [34,35], while other studies have investigated the effects of building layout or urban morphology on LST [36][37][38][39]. Some studies have compared the multiple-city correlations [22,40] based on surface coverage, whereas other studies have been limited to single-temporal LCZ mapping [41][42][43]. Only a few studies have evaluated the SUHI effect in a long time series, such as the study of Wang [27], in which the dynamic distribution of LST was explored.
In summary, the previous studies on the SUHI effect carried out in the GBA have mostly focused on the relationships between mono-temporal LULC characteristics and the SUHI effect, which is an approach that cannot reveal the change patterns of thermal fields in a time series during the urbanization process. Furthermore, the studies that have considered long time series have failed to address microscopic descriptive analysis, and they have ignored the impacts of fine surface structure changes on the SUHI effect in the urbanization process. In addition, it is worth noting that most of the previous studies have only covered parts of the cities in the GBA, and there is a lack of research and analysis covering the whole of the GBA. Given these circumstances, it is clear that carrying out research on the dynamic changes of the LCZs in the GBA and their spatio-temporal effects on the SUHI effect is of great significance.
The main purposes of this study were (1) to propose a multi-temporal LCZ mapping method that can be used to achieve the detection of multi-temporal LCZ changes in the GBA region, and (2) to analyze the relationship between the multi-temporal LCZ changes and LST changes.
The rest of this paper is structured as follows: Sections 2 and 3 introduce the overview of the study area, the data, and the methods; Section 4 provides the results and a discussion of the results; Section 5 summarizes the findings of this study.

Study Area
The GBA (21 • 32 -24 • 26 N, 111 • 20 -115 • 24 E) is located in southern China and includes nine cities in the Pearl River Delta city cluster (i.e., Guangzhou, Shenzhen, Foshan, Dongguan, Zhongshan, Zhuhai, Huizhou, Jiangmen and Zhaoqing), as well as two Special Administrative Regions (i.e., Hong Kong and Macao) ( Figure 1). With multiple financial services, educational resources, cultural industries, and high-tech industries, the GBA is representative of China's efforts to build a world-class city cluster and plays an important strategic role in the overall development of China. However, the rapid urbanization has led to prominent human-nature conflicts, especially environmental problems such as the UHI effect [31,44,45]. During the development of ecological civilization in China, establishing the concept of "people-oriented" development and handling the relationships between "development" and "protection" have become the key to the high-quality development of the GBA. Therefore, how to relieve the SUHI effect in order to make the GBA region more conducive to living and sustainable development is an urgent issue that needs to be addressed.

LST Retrieval
The UHI is strongest in summer due to the high temperature. Thus, the summer temperature is usually analyzed in the UHI studies. The GBA is rainy in summer. In order to reduce the interference of other factors (clouds, etc.), we used the average surface temperature obtained based on multiple remote sensing images. The average daily air temperatures of the GBA in 2005 and 2015 were obtained from the China Meteorological Data Network ( Figure 2). It can be seen that, in both 2005 and 2015, the temperature difference between July and August was the smallest and the fluctuation range was minimal, which means that the climate in July and August was the most stable. However, in 2005, the temperature in June was close to that in May, and the difference was about 3 • C from that in July. Therefore, in order to reduce the influence of the difference in June on the average value, we analyzed the spatio-temporal impacts of LCZ changes on LST changes in the process of urbanization by utilizing the average LST in July and August in 2005 and 2015. The LST was acquired from the EOS-Aqua-MODIS eight-day (MYD11A2) product for daytime (13:30), which has a spatial resolution of 1000 m. The accuracy of this product has been reported to be better than 1 • C [46]. The advantages of the high precision and wide area coverage have resulted in MODIS LST data being widely employed in SUHI research [47][48][49].

Time-Series LCZ Mapping
The common LCZ types in the GBA are shown in Figure 3. LCZ7 (lightweight low-rise) is usually used as informal settlements, shantytowns, smallholder lots, and so on. The building surface fraction and impervious surface fraction of LCZ9 (sparsely built) are both less than 20%. Due to the high level of economic development in the GBA, there are very few lightweight low-rise buildings, and the buildings are relatively dense as well as the impervious surfaces in the built-up area. Consequently, the LCZ7 (lightweight low-rise) and LCZ9 (sparsely built) are almost non-existent in the GBA. LCZC (bush, scrub) refers to the open bushes, shrubs, and short, woody trees on bare soil or sand, which are also rare in the GBA. As a result, we selected 14 common LCZ types. As suggested in existing studies [27,42], LCZ maps of the GBA were generated with a spatial resolution of 100 m in our research. To reduce the interference of factors such as clouds and improve the accuracy of the classification, we obtained the atmospherically corrected (and almost cloud-free) Landsat 8 surface reflectance (SR) products for the GBA for the whole of 2015 from the United States Geological Survey (USGS), and then performed image mosaicking and clipping. After careful visual interpretation, a total of 2987 samples were selected based on the high-resolution images of the same period from Google Earth, of which 25% were randomly used as test samples and the others were used as the training samples. The preparation of the input features consisted of two parts: (1) The multispectral and thermal infrared bands of the Landsat 8 SR products were bilinearly resampled to a 10-m resolution, and then resampled to a 100-m resolution by a zonal mean function based on the LCZ grid area. To consider the contextual characteristics of the features, the mean, minimum, maximum, median, 25th, and 75th quantile values of the nine pixels in a 3 × 3 window were calculated ( Figure 4b) [42]. (2) We also produced input features based on the China Land Use/Cover Dataset (CLUD, http://www.igsnrr.ac.cn/, accessed on 28 January 2021) for the same year. The CLUD dataset was mainly generated from Landsat TM/ETM+/OLI and HJ-1A/1B imagery by a human-computer interactive method and has been reported as having an overall accuracy of more than 90% [50,51]. Specifically, as with the Landsat data, we also calculated the feature values, including the fractions of impervious surfaces, trees, and water, at a spatial scale of 100 m. The random forest classifier was finally employed to complete the LCZ mapping in the GBA for 2015. (1) Unchanged: the change vector analysis (CVA) method was applied to determine whether a pixel had changed using the multispectral and thermal infrared bands of the two time phases [52,53]. The gray values of a pixel for 2005 and 2015 were, respectively, described by G = (g 1 , g 2 , . . . , g n ) T and H = (h 1 , h 2 , . . . , h n ) T , where n represents the number of bands. The change vector was described as: The change magnitude was calculated by: A greater change magnitude corresponds to a higher possibility of image change. When the change intensity exceeds a specific threshold, it can be treated as a changed pixel.
(2) Reliable: in a region where the images had not changed, only the positions with higher reliability in the time 1 image classification could be considered as samples for time 2. The reliability of a pixel was measured by the classification probability of the image [54], and the classification certainty w(x) of pixel x was defined as: where p k (x), p k+1 (x) are descending probabilistic outputs; k refers to the number of probabilistic outputs. (3) Representative: the candidate samples were further refined according to their representativeness in order to reduce the redundancy. Specifically, on the basis of the invariability and reliability of the samples, the k-nearest neighbor algorithm was adopted to further select representative samples [55].
Due to the insufficient samples for the high-rise buildings (LCZ1 and LCZ4), heavy industry (LCZ10), and bare rock (LCZE) categories, some samples needed to be added manually after the sample migration ( Table 1). The classification input features for 2005 were obtained in the same way as for 2015, based on the multispectral and thermal infrared bands of the Landsat 5 SR products and the CLUD dataset, with random forest utilized for the classification. A total of 735 test samples were randomly collected in the region, outside of the training samples, to ensure that the test and training samples were spatially independent of each other.

Analysis of the LCZ and LST Changes
The landscape metrics are quantitative indicators that can effectively reflect the spatial composition and spatial configuration of the landscape, so they could be used to assist with the quantification of the spatial morphology of the LCZs in the GBA. Based on previous studies [40,56,57], four widely used landscape indices were employed (Table 2): the percentage of landscape area (PLAND), the edge density (ED), the largest patch index (LPI), and the aggregation index (AI). PLAND is a spatial composition index, which measures the number of different patch types in the landscape, and the others are spatial configuration indices, which estimate the spatial structure of the patches. FRAGSTATS software was adopted to calculate the landscape indices (http://www.umass.edu/landeco/ research/fragstats/fragstats.html, accessed on 28 January 2021). In order to identify the influence of the changes in the spatial morphology of the LCZs on the changes of LST, the most commonly used multiple linear regression models were applied [58][59][60]. The correlation between each landscape metric of LCZ changes and the LST changes from 2005 to 2015 was analyzed separately. Table 2. The landscape metrics selected in this study.

Metric (Abbreviation) Calculation Description
Percentage of landscape (PLAND) ∑ n j=1 a ij /A × 100 Measures the proportional abundance of each class in the landscape (unit: %).
Largest patch index (LPI) max a ij /A × 100 Measures the dominance of each class in the landscape (unit: %).
Aggregation index (AI) Measures the degree of spatial aggregation (unit: %).
a ij = area of patch ij; n i = number of patches in the landscape for class i; e ik = total length of edge in landscape involving class i, including landscape boundary and background segments involving class i; A = total landscape area; g ii = number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method; max → g ii = maximum number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method. There are 13 categories of LCZ types detected in the GBA. Theoretically, the 13 classes involve 156 (13 × 12) change transition types but, in fact, only 25 change types occur. Some change categories are rare, e.g., buildings (built types) to LCZG (water) or vegetation (LCZB/D), LCZG (water) to LCZA (dense trees), and compact high-rise buildings (LCZ1) to other types. Due to the large number of change categories and the distinct areas of some change types, approximately 30 samples were randomly selected for each change transition type to ensure that the test samples could cover all the change types, with a total of 738 samples. The overall accuracy of the change detection is 76.83% and the Kappa coefficient is 0.759 (Figure 7). Figure 8 shows the percentage of the LCZ area changes in the total area from 2005 to 2015, indicating the degree of the LCZ area changes. As shown in Figure 8, the area of the built types (LCZ1-10) in the GBA shows an increase of 2.38% from 2005 to 2005, among which LCZ8 (large low-rise) shows the largest increase in area of 1.01%, followed by open buildings (LCZ4, LCZ5 and LCZ6), while LCZ3 (compact low-rise) is the only built category that shows a decrease, of 0.51%. Among the land-cover types (LCZA-G), the largest area change is in LCZD (low plants), which shows a significant reduction, and accounts for about 2.17% of the total area of the GBA.

LCZ Classification Results and Analysis of the LCZ Changes
In the following, the results are analyzed from the perspective of each city, considering the changes in LCZ area and gross domestic product (GDP) and gross industry product (GIP) (Figure 9) As inland first-tier and new first-tier cities, Guangzhou, Shenzhen, and Dongguan show the largest increase in the proportion of LCZ4 area (more than 1%), which can be attributed to the high incremental GDP, the strong economic dynamics, and the high verticalization rate of the urban buildings. Among the four second-tier cities, Foshan shows the highest GIP growth, which is close to that of the first-tier cities (i.e., Guangzhou and Shenzhen). More importantly, Foshan's GIP accounts for 58% of the GDP, ranking first in the Pearl River Delta city cluster in 2015. The rapid industrial development in Foshan is partly on account of its geographical location, i.e., the proximity of Guangzhou, making it vulnerable to the economic radiation of Guangzhou. Under these circumstances, the percentage increase (3.46%) in the area of LCZ8 (large low-rise buildings), representing industrial buildings, in Foshan is much higher than in the other cities. In addition, Jiangmen and Zhaoqing, which are two third-tier cities, are relatively distant from Guangzhou and Shenzhen, with weak radiation effects from the big cities. As a result, their GDP and GIP are relatively low, and the amount of increased area in the built types is the smallest among the inland cities. With respect to Macau and Hong Kong, which are two Special Administrative Regions of China, because of the high density of population and the scarcity of available land, the variations in the area change of the built types are the smallest and the urban buildings are extremely compact. Consequently, the LCZ1 (compact high-rise) areas show an increase of more than 0.5%, which is larger than the increase seen in the nine inland cities. Furthermore, the area proportion of LCZG (water) in Macao shows a significant reduction (11.36%), as a result of the construction of the Zhuhai-Macao artificial island.     It should be noted here that economic activity and LCZ changes differ between urban centers and suburban areas. However, LCZ is proposed to eliminate the impact of the inaccurate division between urban and rural areas on the study of the UHI. In addition, we have focused on the variation between cities, and discussed the possible influence of the economy and the policy on the change of the LCZs in the GBA from the macro level, so we did not indicate the urban centers with active economic activities. We will study the urban center delineation method and the relationship between the local changes in the LCZ and the economy in depth in the future.
It is important to pay attention to the fact that the growth rate of the compact building types in 11 cities in the GBA was less than 1% in the process of urbanization. The compact building types are LCZ1 (compact high-rise), LCZ2 (compact midrise) and LCZ3 (compact low-rise), all of which have the building surface fraction greater than 40% and the impervious surface fraction of less than 20% with few or no trees. Therefore, the compact building categories increased in size by less than 1%, indicating that new buildings are rarely laid out in a compact manner, taking into account the feelings of residents and emphasizing the concept of livability, workability and people-orientation. In addition, the forest coverage areas of all the cities in the GBA do not show significant decreases, with the LCZA (dense trees) area of Zhuhai showing the greatest reduction, at only 3.1%. That may be related to the fact that the nine inland cities in the GBA have successfully achieved the titles of National Forest City or National Garden City, adopted and implemented policies and regulations related to the protection of forestland, and explicitly proposed to build a forest city cluster in the GBA. The results signify that, in the process of urban development, the GBA has not overly expanded in a disorderly manner, and there has been a commitment to the protection of forest vegetation and the improvement of the ecological environment.
In the following, we further analyzed the characteristics and reasons for the LST changes based on the analysis of the LCZ changes. Since the main purpose of the LCZ concept is to avoid the influence of the urban-rural division on the study of the UHI, the whole urban area is used as the study object, and no distinction is made between the central urban area and the surrounding rural areas in the LST change analysis.

LST Variation Analysis
The average LST of July to August in 2005 and 2015 in the GBA was obtained based on the MYD11A2 product ( Figure 10). The LST differences between the two time phases were then divided into six levels, and the LST change categories for each city and the composition of the LST change categories in the GBA were analyzed (Figure 11).
It can be observed from Figure 11a that the area with raised temperature is larger than that with reduced temperature in all the cities in the GBA. Among these cities, Zhongshan shows the largest proportion of warming area, accounting for about 94.38% of the total area of Zhongshan, and the level 5 (2 to 4 • C) and level 6 (more than 4 • C) categories have the largest area proportions, at 42.44% and 6.76%, respectively. According to the Overall Urban Planning of Zhongshan City (2005 to 2020), Zhongshan has taken high-tech and modern manufacturing industry development as the central task and has formed a clustered development structure. The warming areas of Zhongshan are mainly found in the Torch High-tech Industrial Development Zone, as a representative of the eastern district, and Henglan Town, as a representative of the northwest district. The eastern district is the concentrated development area for industrial innovation, high-tech industries, heavy industries, and port industries, and the northwest district has enlarged the industrial agglomeration effects and actively developed comprehensive service industries. Both areas have been developed as sub-centers of the city. Therefore, during the urban development process, both areas have mostly transformed from LCZG (water) and LCZD (low plants), which both have cooling effects, to LCZ8 (large low-rise), which produces a large amount of high-temperature wastewater and exhaust gases, which have caused the LST to rise. The areas of rising LST have also had a thermal radiation effect on the surrounding areas. At the same time, the process of creating sub-centers of the city has led to an increased concentration of population and activities, resulting in more anthropogenic heat and an increase in the area of the habitable open built types (LCZ4, LCZ5 and LCZ6). Accordingly, the LST has increased in many regions, and the thermal environment has been strongly influenced by urbanization in Zhongshan (Figure 12).  The city with the largest proportion of cooling area is Guangzhou, where it accounts for 49.05% of the total urban area. The cooling areas are mainly concentrated in the dense tree areas such as Nankunshan and Shimen. By paying attention to the construction and protection of forest areas, Guangzhou became the first city in the GBA to obtain the title of National Forest City. During the gradual expansion of LCZA (dense trees) and the increase of the vegetation canopy, the shading and evapotranspiration effects of vegetation have become pronounced, which has contributed to the significant cooling of the ground surface. However, as a whole, the warming area that accounts for 68.44% of the total area is much larger than the cooling area in the GBA, where the level 4 (0 to 2 • C) and level 5 categories (2 to 4 • C) of LST change account for 55.39% and 12.22%, respectively. In the warming area, 5.97% of the region has been converted from land-cover types to built types, resulting in an increase in impervious surfaces and a decrease in specific heat capacity, as well as an enhancement of the anthropogenic heat generated by human activities, which has led to the rise in LST. In addition, 3.11% of the warming area consists of new large low-rise buildings (LCZ8), which are mainly found in industrial areas, and 4.9% consists of new built categories, with high population concentration and strong human activities (LCZ1, LCZ2, LCZ4 and LCZ5). It should be noted that 75% of the warming area does not show any LCZ change, so the reason for the increase in LST might be enhanced human activities, increased heat generation, or the influence of thermal radiation from the surrounding high-temperature areas. Of the cooling area, 95% is covered by land-cover types, 91.59% of which is made up of LCZA (dense trees), LCZD (low plants), and LCZG (water), indicating that vegetation and water have the most dominant role in mitigating the SUHI effect.
The composition of the LST change levels in the GBA is plotted in Figure 11b. From level 1 to level 6, the area ratio of each city was gradually similar. The level 1 category is mainly found in Jiangmen and Zhaoqing, and the area ratios of Jiangmen and Zhaoqing in the cooling area (levels 1-3) are the largest, accounting for about 31.56% and 32.09%, respectively. In accordance with the Overall Urban Planning of Jiangmen City approved in 1996 and 2011, Jiangmen was positioned as a light industrial city with high-tech industries and a waterfront city with modern manufacturing industries, commercial logistics industries, and cultural tourism. In addition, Zhaoqing was planned as a national historical and cultural city, in terms of the Overall Urban Planning of Zhaoqing City approved in 2010. Under the influence of these policies, the pace of the economic development of these two third-tier cities has been relatively slow, and the development of heavy industries has lagged behind the other cities, leading to slower urban expansion and a lesser area of new buildings. As a result, the new impervious surface contributing to surface warming [61,62] and the anthropogenic heat promoting surface thermal radiation have remained at a relatively low level. Simultaneously, the two cities have paid attention to the protection of ecological resources, and LCZA (dense trees) accounts for 33.4% of the GBA, which benefits the development of cultural tourism. As time has passed by, the density of the vegetation canopy has increased, so that the shading and evapotranspiration effects of vegetation have become more prominent, which can significantly reduce surface temperatures [10,40,56]. Therefore, the cooling areas of Jiangmen and Zhaoqing are the largest, and they are playing an important role in mitigating the SUHI effect of the entire GBA. Among the warming regions (levels 4 to 6), the proportions in Zhaoqing, Jiangmen, and Huizhou are higher, because of the larger total area (15,000 km 2 , 9505 km 2 , and 11,300 km 2 , respectively) and the larger area of buildings with warming effects, even though these three cities have the largest area of LCZA (dense trees) (12,337 km 2 , 6197 km 2 , and 8478 km 2 , respectively), which has substantially impacted the cooling effect. Therefore, it can be stated that when establishing measures to mitigate the SUHI effect, it is important to focus on not only the regions with rapid industrial development and high-intensity human activities, but also the regions with lagging economies but large built-up areas. In addition, we compared the proportions of the warming areas and cooling areas in each city. The warming area is 11% larger than the cooling area in Zhuhai and Zhongshan, and 6% larger in Shenzhen, Dongguan, and Foshan. In Hong Kong and Macao, the warming area is about the same as the cooling area. Nevertheless, the cooling area of Guangzhou, Jiangmen, and Zhaoqing is 10 to 20% larger than the warming area. Therefore, Zhuhai and Zhongshan are the cities in the GBA where the SUHI effect still needs to be addressed most urgently.

Influence of LCZ Changes on LST Changes
The PLAND, ED, LPI, and AI metrics were applied to characterize the spatial morphology of the LCZ distribution on the data of 1822 samples at a grid size of 5000 m × 5000 m. Table 3 shows the results of the multiple linear regression analysis for the LCZ spatial morphology changes and LST changes, indicating the effect of LCZ spatial composition and configuration changes on LST changes. The variance inflation factors of all the variables are between 1 and 2, demonstrating that multicollinearity does not exist. Table 3. Multiple linear regression analysis of LCZ spatial morphology changes and LST changes. Each column represents a multiple linear regression model for the landscape metric changes of different LCZ categories and LST changes. The explanatory variables are the values of LCZ landscape metric changes, i.e., the differences in each landscape metric of different LCZs in Table 2. from 2005 to 2015. The objective variable is LST change from 2005 to 2015 in Figure 10. The values in the color region refer to the standard regression coefficients, which reflect the importance of the corresponding explanatory variable. The color represents the magnitude of the regression coefficients, where the redder the color, the stronger the positive effect of the variable on LST, and the bluer the color, the stronger the negative effect of the variable on LST. "**" refers to the significance less than 0.01, indicating that the data is extremely significant, meanwhile, "*" refers to the significance less than 0.05, indicating that the data is significant. The results reveal that, for the built types (LCZ1-10), the degree of LST increases as PLAND, ED, LPI, and AI increase. The main reasons for this are as follows: (1) the larger the increased area of buildings, the greater the range of human activities and the larger the area of anthropogenic heat generation. At the same time, the increase of impervious surface area with a low specific heat capacity causes LST to rise after absorbing solar radiation heat and artificial heat. (2) ED measures the complexity of the shape and degree of the isolation between the patches. An increase in ED represents an increase in the fragmentation of the patches, resulting in an increase in the thermal radiation range of the built types and an increase in LST. (3) LPI can recognize the dominant category in the spatial range. An increase in LPI means that the dominant position of the built types is becoming more and more prominent, resulting in a greater degree of surface temperature increase. (4) An increase in AI indicates intensification of the population concentration and human activity in this region, which promotes an increase in anthropogenic heat emissions. It should be pointed out that, although the ecological significances of ED and AI are mutually exclusive to a certain extent, the R 2 of the AI model is much smaller than that of the ED model, indicating that the ED model has a stronger explanatory power for the LST changes. For this reason, we should try to avoid a fragmented distribution of built types and reduce the effective contact area of buildings with the surrounding environment in the urban planning and design.

PLAND
Among the built types (LCZ1-10), the morphology changes of LCZ4 (open high-rise) and LCZ5 (open midrise) are both strongly correlated with the changes of LST, which might be due to the fact that LCZ4 (open high-rise) and LCZ5 (open midrise), as the main built types for urban residential areas, are the categories with the newest area, and less area converted to other types. The increased impervious surfaces and anthropogenic heat has led to a significant increase in LST in these two categories.
With regard to the land-cover types (LCZA-G), LCZA (dense trees), LCZD (low plants), and LCZG (water) show opposite trends. The results reveal that the greater the variation in area, patch fragmentation, dominance, and aggregation of these categories, the more favorable the LST reduction, which is correlated with the evapotranspiration effect of vegetation and water and the shading effect of trees. What is more, it should be noted that the LPI changes in LCZA (dense trees) correlate most strongly with the LST changes, which indicates that the extent of the aggregation of dense trees has a significant cooling effect on LST. Therefore, attention needs to be paid to the aggregation effect of forest in mitigating the SUHI effect.
When comparing the difference between natural and built-up landscapes with regard to their effects on LST, it is found that the standardized regression coefficients of LCZ4 (open high-rise) and LCZ5 (open midrise), with the strongest positive correlation with LST changes, are basically larger than those of LCZD (low plants) and LCZG (water), which have a cooling effect on LST. The conclusion we can draw is that the influence of the increase of LCZ4 (open high-rise) and LCZ5 (open midrise) on LST is much greater than that of LCZD (low plants) and LCZG (water). As analyzed above, to mitigate the urban warming effect effectively, rational layouts should be based on the correlation between built-up landscapes and LST. More specifically, LCZD (low plants) and LCZG (water) should be expanded within a certain range, and a large patch area for LCZ4 (open high-rise) and LCZ5 (open midrise) should be avoided, to weaken the dominant position of the built types. Meanwhile, the degree of patch fragmentation in the built types should be reduced to mitigate the effect of thermal radiation on the surrounding area.
In addition, it is commonly recognized that the grid size can have a significant influence on the relationship between temperature and urban surface structure, so a suitable analysis cell is crucial for the study of the spatial characteristics of the SUHI effect [40,63,64]. Since the pixel size of the LCZ classification images produced for this study was 100 m × 100 m and the acquired MODIS LST products were 1000 m × 1000 m, the side length of the grid should be an integer multiple of 1000 m. As a consequence, the strength of the linear regression (R 2 values) was used to examine the effect of grid sizes with side lengths ranging from 2000 m to 5000 m in increments of 1000 m. As presented in Figure 13, the linear regression strength (R 2 values) of the model shows a slowly increasing trend with the increase of the grid size, and the peak value appears at the grid side length of 5000 m. Considering the slow increase of R 2 and the sudden decrease of the sample size, it was not advisable to increase the grid size anymore. Thus, 5000 m × 5000 m was determined as the appropriate grid size in this study. Figure 13. Sensitivity of the correlation between LCZ spatial morphology changes and LST changes to the grid size.

Conclusions
In this study, with the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) as the study area, we derived LCZ classification maps of the GBA by a time-series sample migration method and obtained the LST based on MODIS data from 2005 and 2015. We then systematically investigated the dynamic change characteristics of the LCZs and LST in the urbanization process of the GBA, and quantitatively analyzed the influence of LCZ changes on LST variation through calculating the changes of LCZ spatial composition and configuration by the use of landscape indices.
Our findings suggested that the time-series sample migration method can be used to achieve rapid and efficient mapping of LCZs with multiple time series in a large area. The accuracy of the LCZ mapping for the GBA in 2005 and 2015 was 85.03% and 85.28%, respectively. The built type category with the largest increase in area in the GBA from 2005 to 2015 was LCZ8 (large low-rise), with a 1.01% increase, while LCZ3 (compact low-rise) was the only built type with a decrease in area, of 0.51%. Among the land-cover types (LCZA-G), the largest area change was in LCZD (low plants), which decreased by 2.17%. The changes in the LCZ categories varied among the cities due to the different factors, such as the economic development level and local policies.
During the urbanization process, we found that the area of warming was larger than that of cooling in all the cities of the GBA. The main reasons for the increase of LST were the conversion of land-cover types into built types, the enhancement of human activities, and the heat radiation from surrounding high-temperature areas. The characteristics of the LST changes varied from city to city because of the influence of local planning, the area of the region, and the economic development level.
Another contribution of this study is that we have revealed the different effects of LCZ spatial morphology changes on LST variation. The spatial morphology changes of the built categories were found to be positively correlated with LST changes, and the morphological changes of LCZ4 (open high-rise) and LCZ5 (open midrise) exerted the most significant influence, while the change of the LPI of LCZA (dense trees) had the strongest negative correlation with LST variation. Consequently, there is a need to focus on the rationalization of the layout of the built types, and attention also needs to be paid to the agglomeration effect of forest in urban planning.  Institutional Review Board Statement: The study did not involve humans or animals.

Informed Consent Statement:
The study did not involve humans.

Data Availability Statement:
The study did not report any data.