Temporal and Spatial Variation of Anthropogenic Heat in the Central Urban Area: A Case Study of Guangzhou, China

: The urban heat island effect caused by the rapid increase in urban anthropogenic heat has gradually become an important factor affecting the living environment of urban residents. Studying the temporal and spatial variation characteristics of urban anthropogenic heat is of great signiﬁcance for urban planning and urban ecological service systems. In this study, the urban anthropogenic heat ﬂux (AHF) in 2004, 2009, 2014, and 2020 in the central urban area of Guangzhou was retrieved based on Landsat data and the surface energy balance equation, and the temporal and spatial characteristics of different types of anthropogenic heat were explored by combining the transfer matrix and the migration of the gravity center. The results showed that: (1) The overall change trend of anthropogenic heat in the central urban area of Guangzhou was enhanced, and the degree of enhancement was related to the type of urban functional land. (2) Different types of anthropogenic heat had different characteristics in terms of area expansion and spatial changes. Low-value anthropogenic heat (zero-AHF zone, low-AHF zone, medium-AHF zone) changed drastically in terms of area expansion. High-value anthropogenic heat (medium-AHF zone, high-AHF zone) changed more drastically in space. The increase in urban population, rapid economic development, and increased industrial production activities have stimulated the emission of anthropogenic heat, which has a positive impact on the intensity of anthropogenic heat.


Introduction
Urbanization, modifying the natural environment with an artificial environment, has several impacts on urban environment and climate, which are global concerns [1][2][3]. The heat emission of human production activities is increasing rapidly, making the city a huge heating element and exacerbating the temperature difference between the city and its surrounding suburbs, which is an important factor causing the increase of the urban heat island effect [4][5][6]. Anthropogenic heat is an additional source of energy in urban areas compared to the other surfaces of the earth and results mainly from human metabolism and the consumption of energy by human activity, such as heating of residential areas, commercial centers, industry, and the discharge from vehicles [7][8][9][10]. Urban anthropogenic heat is an accompanying product of urbanization, changing the energy flow of the urban ecosystem and affecting the urban ecological processes such as the urban regional climate and atmospheric environment [11]. Therefore, quantification of anthropogenic heat is of great significance to mitigation and control of the impact of urbanization on the environment and human society.
With the growing concern over anthropogenic heat flux (AHF), many researchers have conducted quantitative studies in this field. The methods of estimating anthropogenic heat can be grouped into three major categories: the energy-consumption inventory method, impact on the thermal environment. In this study, taking the urban area of Guangzhou as an example, based on Landsat time-series image data and the surface energy balance equation, the urban anthropogenic heat for winter from 2004 to 2020 was estimated. The objective of this study was to analyze the spatial distribution and the characteristics of the evolution and explore the evolutionary laws and driving factors of the urban anthropogenic thermal environment.

Study Area
Guangzhou is the economic and political center of South China, located in the northern edge of the Pearl River Delta. The study area is the central urban area of Guangzhou (as shown in Figure 1), including 6 districts of Liwan, Yuexiu, Haizhu, Tianhe, Baiyun, and Huangpu, with a range of 113 • 8 to 113 • 36 E and 23 • 2 to 23 • 25 N. The topography of the study area is high in the northeast and low in the southwest, and it has a subtropical maritime monsoon climate with an average annual temperature of about 21 • C. It covers an area of approximately 1480 km 2 and has a permanent population of approximately 9.67 million in 2019. As a core city in South China, Guangzhou is in a period of rapid economic development, significantly promoting the development of urbanization. The massive increase in urban population and the enhancement of human activities have a huge impact on urban anthropogenic heat, such as the heat island effect, ecological damage, and eco-environmental effects and incipient crises in energy availability.

Data and Preprocessing
The main data used in this study include remote sensing images and meteorological data. The ancillary data include statistical data of population and gross domestic product (GDP) in 2004,2009,2014, and 2020 from the statistical yearbooks, and the distribution map of the urban functional land from Guangzhou Municipal Planning and Natural Resources Bureau.
Remote sensing images include a Landsat 5 TM image (acquired on 21 January 2004 and 2 January 2009) and a Landsat 8 OLI/TRS image (acquired on 16 January 2014 and 18 February 2020), provided by United States Geological Survey (https://earthexplorer. usgs.gov/, accessed on 13 January 2021) and Geospatial Data Cloud Platform of Computer Network Information Center of Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 13 January 2021), the path/row of which is 122/044. The spatial resolution for the sensor of Landsat5 TM and Landsat8 OLI/TRS is 30 m in the short-wave band and 120 m (TM image),100 m (TRS image) in the thermal band. The detailed data are shown in Table 1. The projection of all the images are UTM WGS84. They were preprocessed with radiometric calibration and atmospheric correction by using FLAASH atmospheric correction method based on ENVI 5.3 software. The land cover classification has been taken by the supervised classification method, which was divided into six types: construction land, forest land, cultivated land, grassland, water body, and bare land. Accuracy assessment is performed by the true land cover types, which is digitized using ArcGIS based on Google Earth image. The overall accuracy of each year is above 90%, and the Kappa coefficient is above 0.8.
The meteorological data came from China Meteorological Data Network (http://data. cma.cn/site/index.html, accessed on 13 January 2021). The selected weather station is Guangzhou Station (23 • 10 N, 113 • 20 E) and its district station number is 59278. The meteorological data have undergone strict quality control and inspection, which includes atmospheric temperature, vapor pressure, saturated vapor pressure, wind speed, air density, and total solar radiation. The detailed data are shown in Table 2. The radiative transfer equation method [35] is used to estimate the surface temperature. The algorithm is specifically expressed as: where T s is the surface temperature in K, and L λ is the radiance value of the infrared band. ε represents the surface emissivity. τ, L ↑ , L ↓ are the atmospheric transmissivity, the upwelling radiation and downwelling radiation respectively, which can be obtained from NASA (http://atmcorr.gsfc.nasa.gov/, accessed on 13 January 2021). K 1 , K 2 represent the calibration parameter of the infrared band: K 1 =607.76 W·m −2 ·µm −1 ·sr-1, and K 2 = 1260.56 K for the band 6 of Landsat TM; K 1 = 774.89W·m −2 ·µm −1 ·sr-1 and K 2 = 1321.08 K for the band 8 of Landsat OLI.
The surface emissivity ε is estimated based on the following equation [36] in this study: where NDVI is the normalized difference vegetation index, and NDV I min and NDV I max are the minimum and maximum value, respectively. B N IR and B R represent the near-infrared band and red band of the Landsat images, respectively.

Estimating Anthropogenic Heat Flux Based on the Surface Energy Balance Model
Anthropogenic heat flux (AHF) was generally used to represent the intensity of anthropogenic heat. Here, the surface energy balance model [18] was used to calculate AHF, which can be expressed as: where R n is the net radiation (W/m 2 ), A is the AHF (W/m 2 ), G is the ground heat flux (W/m 2 ), H is the sensible heat (W/m 2 ), and E is the latent heat flux (W/m 2 ). In the energy balance equation, the net radiation (R n ) is the sum of solar shortwave radiation and longwave radiation reaching the land surface. It is the net solar radiation energy actually received by the earth's surface after being reflected by the surface. R n is determined by parameters such as surface albedo, specific emissivity, and surface temperature. It is calculated as: where R S ↓ is the downwelling shortwave radiation (W/m 2 ), which can be described as the total solar radiation retrieved from the meteorological data. α and ε represent the surface albedo and the surface emissivity, respectively, and are derived from the remote sensing data [37]. ε a is the atmospheric emissivity, which is estimated using an empirical equation from Kato and Yamaguchi [38]. σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·K 4 /m 2 ), and T a represents the atmospheric temperature in K, which can be obtained from meteorological data. T s is the land surface temperature in K, which is estimated based on Equation (2). The ground heat flux (G) represents the heat exchange energy between the surface and deep layers of the soil. In order to accurately estimate soil heat flux, it is necessary to understand the thermal conductivity of the surface material and the vertical temperature distribution. However, it is difficult to obtain the parameters throughout the study area. There is a linear relationship between soil heat flux and net surface radiation in the existing studies [19,22]. Therefore, it can be estimated from the net surface radiation flux as the following equation: where C g is the coefficient, which changes according to factors such as land cover type and season. C g is related to the heat capacity of the surface material and its thermal conductivity. The higher the heat capacity, the lower the C g value. The higher the thermal conductivity, the higher the C g value [19,22,39]. The value of C g in this study is shown in Table 3. The sensible heat (H) is the energy exchanged between the surface and the atmosphere, characterizing the energy of heating the air. It is estimated based on the following equation: where ρ is the air density in kg/m 3 . C P is the specific heat of air at constant pressure (J/kg·K). r a is the aerodynamic resistance in s/m which is calculated using the following equation: where z m and z h represent the height of wind and humidity measurement in meters, respectively. d is the zero-plane displacement height in meters. z 0m is the roughness length governing momentum transfer in meters, and z 0h is the roughness length governing transfer of heat and vapor in meters. k is the von Karman constant (k = 0.41), and u z is the wind speed in m/s. The values of z 0m , z 0h , and d in this study refer to the existing research results, as shown in Table 3 [19,39]. Latent heat flux (L) refers to the heat flux exchanged between the underlying surface and the atmosphere, including surface water evaporation and vegetation transpiration, which is an important part of the water balance. Latent heat flux is always related to parameters such as air density, saturated vapor pressure of the underlying surface, and aerodynamic resistance. Here, we used the formula proposed by Kato [19] to calculate L: where e s is the saturation vapor pressure; e a is air vapor pressure; and γ is the psychrometric constant in hPa/K, which is related to the vapor pressure in the air and the temperature. r a is aerodynamic resistance in s/m, and r s is the stomatal resistance in s/m, which is estimated using the calculation method proposed by Nishida [38,40,41]. ρ and C P have the same meanings as in Equation (8).
The workflow diagram for the overall methodology is shown in Figure 2.

Transition Matrix
In this study, we used the mathematical differences to quantify the changes in anthropogenic heat, which could be demonstrate as: where S A is anthropogenic heat change, and  In order to quantitatively explore the magnitude and law of conversion in the anthropogenic heat from 2004 to 2020, the spatial analysis tool was used to establish the transfer matrix of different levels of anthropogenic heat. According to its values, AHF was divided into four types: (1) low-AHF zone, where the value of AHF is 0-100 W/m 2 ; (2) medium-AHF zone, where the value of AHF is 100-200 W/m 2 ; (3) high-AHF zone, where the value of AHF is greater than 200 W/m 2 ; and (4) zero-AHF zone, where the value of AHF is below 0 W/m 2 .
The general form of the transition matrix is as follows: where S is the land area of anthropogenic heat; n represents the number of types of anthropogenic heat; i is the anthropogenic heat type at the beginning of the time period; and j is the anthropogenic heat type at the end of the time period. When i = j, S ij represents the unchanged area of the anthropogenic heat type; when i = j, S ij represents the area transferred from anthropogenic heat type i to anthropogenic heat type j.2.3.4. Migration of Gravity Center. The gravity center refers to the point where the forces in various directions can maintain balance in the regional space, which is an important spatial index describing the spatial transformation of geographic elements. The gravity center can reflect the overall heterogeneity of the spatial orientation, and its dynamic migration can reflect the overall migration trajectory of the spatial distribution of geographic elements. We can clearly understand the spatiotemporal evolution pattern in terms of space through researching the shifting of the gravity center [42]. Different anthropogenic heat level types are introduced into the concept of gravity center in this research. The gravity center of anthropogenic heat is the geometric center of the same anthropogenic heat level type in the area, and the gravity center of anthropogenic heat deviates from the geometric center for an uneven distribution of anthropogenic heat in the area, for which the offset direction points to the "high density" location of anthropogenic heat. For multiple geographic objects, the coordinates of the gravity center can be expressed as: where X and Y represent the latitude and longitude coordinates of the gravity center of different types of anthropogenic heat. x i and y i indicate the latitude and longitude coordinates, respectively, of each patch of anthropogenic heat, and n is the number of patches. The shift distance of the gravity center can be used to express the spatial change characteristics of the feature type, and the expression is as follows: where (X t+1 , Y t+1 ) and (X t , Y t ) are the gravity center of elements at different times.

Temporal and Spatial Changes of AHF
The spatial distribution characteristics of AHF in the central urban area of Guangzhou were obtained and are shown in Figure 3. It is clear from the results that the AHF distribution in the research area changed significantly from 2004 to 2020. The anthropogenic heat expanded rapidly in the space and increased rapidly in terms of value.  The intensity of anthropogenic heat in 2004 was lower than that in 2020, and the types of anthropogenic heat change are shown in Figure 4a. From the statistics in Table 5, it can be seen that the largest area is the AHI type, with an area of 323.44 km 2 , accounting for 21.99% of the research area. Next are the AHIS and AHS types, with 231.7 and 77.54 km 2 , accounting for 15.75% and 5.27%, respectively. The smallest area is the AHD type, with 4.58 km 2 , accounting for 0.31%.  A comparison was made between the urban functional land and the anthropogenic heat change. The distribution map of the urban functional land in the central city of Guangzhou is shown in Figure 4b. It can be seen that the type of AHI was mainly distributed in the southwest of the center of the study area, and the main functional land types here are residential land and commercial land. The anthropogenic heat in residential and commercial areas mainly comes from electrical energy consumption and human metabolism. There are densely populated and densely built areas in the southwest of the center of the study area, which have been the main reason for a slight increase in anthropogenic heat. The AHIS type was mainly distributed in the outskirts of the central city where the main types of functional land are industrial land and airport land. The anthropogenic heat of industrial land mainly comes from the use of electricity and coal, and frequent air transportation in the airport area causes a lot of waste heat emissions. The new Baiyun International Airport and other transportation facilities have been constantly improving, and their convenient transportation attracted a large number of industrial migration and construction to the mainland [43]. Therefore, the release of a large amount of waste heat because of industrial production and large-scale logistics and transportation made anthropogenic heat rise sharply [39].

Analysis of the Transition Matrix
The distribution of different levels of AHF in 2004 and 2020 is shown in Figure 5. It is clear to see that the area of medium-AHF zone expanded significantly, while the area of low-AHF zone shrunk. The area changes of different levels of AHF are shown in Table 6. During the period 2004-2020, the zero-AHF zone and low-AHF zone accounted for the largest area, exceeding 80%. The next is the medium-AHF zone, and the smallest area is the high-AHF zone, with less than 0.3% of the total research area. From the overall dynamic process, it can be seen that the area of the medium-AHF zone and high-AHF zone increased, while the area of the zero-AHF zone and low-AHF zone decreased. During the study period, the zero-AHF zone decreased significantly, which accounted for a decrease from 63.70% to 56.69%, reduced by 103.06 km 2 , and the period of the fastest decrease was during 2004 to 2009, from 63.70% to 58.86%, with the area decreased by 71.09 km 2 . For the low-AHF zone, its area increased first and then decreased: during the period 2004 to 2009 and 2009 to 2014, the area of low-AHF zone was increasing steadily, with the proportion increased from 35.85% to 40.86% and the area increased by 73.79 km 2 . During the period 2014 to 2020, the area of low-AHF decreased significantly, from 40.86% to 26.56%, with the area decreased by 210.41 km 2 . During the period 2004 to 2020, the proportion of the medium-AHF zone increased significantly, from 0.46% to 16.48%, with the area increased by 235.65 km 2 , and the fastest growth period was during 2014 to 2020, with the area increased by 213.96 km 2 . The proportion of the high-AHF zone increased slightly, from 0.001% to 0.28%, with the area increased by 4.05 km 2 . Therefore, during period of 2004 to 2009 and 2009 to 2014, the transfer-out of the zero-AHF zone and the transfer-in of the low-AHF zone were the main changes, while during the period of 2014 to 2020, The transfer-out of the zero-AHF zone and low-AHF zone and the transfer-in of the medium-AHF zone were the main changes. In all, the level types of AHF in the central urban area from 2004 to 2020 changed extremely significantly, and the overall situation is unbalanced.
According to the transfer matrix (Table 7), the direction and quantity of changes among different level types of AHF could be analyzed. It can be seen that the various level types of AHF were transferred in and out during the period 2004-2020, and the transfer volume between types became more complicated.   The low-AHF zone was the largest of the transferred-out areas, with 217.06 km 2 , which mainly transferred to the medium-AHF zone. The next was the zero-AHF zone, with an area of 132.19 km 2, which mainly converted to the low-AHF zone and medium-AHF zone, with transfer rates of 8.41% and 5.57%, respectively. The smallest transfer-out area was the high-AHF zone, with area only of 0.01 km 2 .
The largest of the transfer-in areas was the medium-AHF zone, with an area of 237.64 km 2 , and the main contribution was from the zero-AHF zone and low-AHF zone. The second was the low-AHF zone with an area of 80.43 km 2 , and the main contribution came from the non-built-up area. The smallest was the high-AHF zone of 4.06 km 2 .
The absolute number of transfer-in and transfer-out areas of the low-AHF zone was relatively large; thus, the changes in low-AHF were the most intense, while the number of transfer-in and transfer-out areas in the high-AHF zones was small, which indicated that the high-AHF zones were relatively stable.

Spatial Migration of Gravity Center
The gravity center of different levels of AHF in the central urban area of Guangzhou is shown in Figure 6 and can be used to explore the migration trace of different levels of AHF and reveal the temporal and spatial evolution of anthropogenic heat. This study used the center of the study area as a reference. Spatially, the gravity center of the low-AHF zone moved to the northeast during 2004-2020, but the gravity centers of the low-AHF zone in the four years were all located in the southwest near the center of the study area, which indicated that the distribution of the low-AHF zone was uniform, and the southwest is relatively denser than the northeast. The gravity center of medium-AHF zone migrated to the northwest, from the southeast of the study area, and the gravity centers of medium-AHF zone in the four years were all located in the southern part of the center of the study area. This indicated that the southern part of the medium-AHF zone was denser than the north. The gravity center of the high-AHF zone migrated to the southwest. In 2004 and 2009, the gravity centers of the high-AHF zone were located in the northeast of the study area. And then it migrated to the north, located in the southeast of the center of the study area in 2014, while it was located in the southwest of the center of the study area in 2020.
The migration distances of different levels of AHF are shown in Table 8, which were significantly different. The shortest migration distance was the low-AHF zone with 1.31 km during 2004 to 2020, and its migration distances were 0.32 km in 2004-2009, 0.59 in 2009-2014, 0.55 km in 2014 to 2020, respectively. The next was the medium-AHF zone with 13.19 km, and its migration distances were 3.26 km, 5.3 km, and 5.7 km, respectively. The longest migration distance was the high-AHF zone with 19.06 km, and its migration distances were 2.31 km, 12.2 km, and 16.24 km, respectively. Based on Table 6 and Figure 5, the area and distribution of low-AHF played the main control role in the study area during the period 2004-2020. Therefore, the gravity center of migration distance of the low-AHF zone was relatively short. The medium-AHF was mainly distributed in the southeast of study area in 2004, while the medium-AHF gradually became the main type in 2020, concentrated in the southwest and northwest, and scattered in the north. Therefore, the gravity center of medium-AHF had a large shift. The high-AHF zone had a small change in area but a large change in space. The high-AHF zone was mainly located in the east in 2004. However, it expanded to the north and southeast in 2020. To sum up, the low-AHF zone had relatively stable spatial changes during the period from 2004 to 2020, while the medium-AHF zone and high-AHF zone changed more drastically.

Validation of the Estimated Heat Fluxes
The surface heat balance mainly depends on the surface material and local aerodynamic conditions. Since Meteorological parameters vary with location and time, more and accurate spatial meteorological data can reduce the error of anthropogenic heat flux estimation. However, due to the difficulty in obtaining multi-station meteorological data, we only used the data from one meteorological reference station to estimate anthropogenic heat flux. Moreover, there is a lack of anthropogenic heat flux measurement data of the study area. Here, we test the accuracy of heat flux estimated in this study by comparing the results of other studies in the same season (Table 9). Table 9. Heat fluxes ratios to net radiation between the present study and previous studies [19,38,39,44]. In the previous studies, the ratios of H/R n and L/R n varied from 0.08 to 0.46 and from 0.00 to 0.10, respectively [19,38,39,44]. The ratios of H/R n and L/R n in the present study are in the range of above studies, varying from 0.1 to 0.19 and from 0.03 to 0.07. The values of H/R n in this study are similar to the results of Nagoya, Japan in 2004 and Xiamen Island, China in 2009, but lower than the results of Nagoya, Japan in 2000 and Tokyo, Japan in 2002, which is mainly due to the incident solar radiation and the different characteristics of the surface [39]. The values of L/R n in this study are similar to the results of Tokyo, Japan in 2002, but higher than the results of Nagoya, Japan in 2000 and 2004 and Xiamen Island, China in 2009, which are related to the proportion of impervious surfaces [44]. We also compared the distribution of the mapped AHF in 2009 with the result by using the energy-consumption inventory based on top-down approach in 2010. Both distribution patterns are roughly the same, showing a trend of low AHF value in the central city and high AHF value in the edges in tiny spots [30]. Through the above comparative analysis, it is found that the results of this study are consistent with the results of other existing studies. Therefore, the anthropogenic heat flux estimated is credible to some extent.

Analysis of Influencing Factors of AHF
From 2004 to 2020, the total population of Guangzhou's central urban area increased from 3.45 million to 9.67 million, with a net increase of 6.22 million and an increase rate of 389,000 per year. Population is one of the most active factors in the urban structure [45] and has an important impact on AHF. Population data can be used to express the intensity of human activities; thus, the intensity of population activities is more intense in densely populated areas. A correlation analysis between the population data of each district and its AHF in 2004, 2009, 2014, 2020 was carried out to explore the relationship between AHF and population in various districts (Figure 7a). It can be seen that AHF in each district has a positive correlation with the total population (R 2 = 0.3046), and correlation is significant at the 0.01 level (two-tailed). As reported by Dong Y et al. [23], on a global scale, the areas accommodating a dense population have a relatively higher values of AHF than the areas with low population density. The increase in the urban population will inevitably be reflected in the strengthening of urban residents' demand for housing, transportation, and public facilities [46]. Therefore, the urban population can most directly affect anthropogenic heat, and the increase in the total population will lead to an increase in the intensity of anthropogenic heat.
In 2004, the total GDP of the central urban area of Guangzhou was 213.65 billion RMB. It grew to 1693.78 billion RMB by 2020, with a growth rate of 92.502 billion RMB/a. Here, we also used the GDP of each district and its AHF in 2004, 2009, 2014, 2020 to take a correlation analysis (Figure 7b). It can be seen that the intensity of AHF in each district has an obvious linear relationship with the total GDP of districts (R 2 = 0.3425), and correlation is significant at the 0.01 level (two-tailed). Previous studies have also found that the mid-latitudes of the northern hemisphere contribute the main global anthropogenic heat [23]. Countries with high GDP are mainly distributed in this region, such as the United States, China, Japan, Germany, and the United Kingdom. There is a big difference in anthropogenic heat between east and west in China [29], and space distribution of China's anthropogenic heat is similar to economic distribution characteristics, especially in some economically developed large cities, in which the intensity of anthropogenic heat is significantly higher than that in nearby areas [47]. The development of the economy and increased income has also stimulated the emission of anthropogenic heat from residential areas and industrial production activities.
In Figure 4b, the total area of industrial land in the central urban area of Guangzhou was 262.30 km 2 in 2020, and it accounts for 17.83% of the study area. An 8 × 8 grid covering the study area was established by ArcGIS, and the proportion of industrial land and the average AHF value in each grid were used to analyze the correlation relationship between the AHF and industrial land in this study (Figure 7c). It can be seen that AHF of each grid is linearly related to the proportion of industrial land area (R 2 = 0.4273), and correlation is significant at the 0.01 level (two-tailed). It indicates that the industrial land has a positive influence on AHF, and AHF increases as the area of industrial land increases. Previous studies showed that more than 70.0% of anthropogenic heat emissions in Mainland China come from industrial production [29], and the anthropogenic heat is estimated to have spatial heterogeneity [15]. The distribution of industrial land could contribute to the distribution of anthropogenic heat, and the type and density of industrial land have an influence on the intensity of anthropogenic heat.

Limitations
The anthropogenic heat was effectively estimated based on the energy balance model. However, there are some limitations in this study. First, February 2020 was in a

Limitations
The anthropogenic heat was effectively estimated based on the energy balance model. However, there are some limitations in this study. First, February 2020 was in a period of the COVID-19 pandemic when human activity and production had been limited. AHF is affected by human activities. Therefore, the estimation of the anthropogenic heat in 2020 may be lower than expected due to the influence of the COVID-19 pandemic. Second, the aerodynamic parameters used in this study are empirical values from previous studies [19,39]. However, the aerodynamic parameters vary with location and time, which makes the experimental results have certain limitations. A more accurate estimation of anthropogenic heat would be obtained if the spatial interpolation of aerodynamic parameters can be realized in real time. Moreover, the spatial pattern of the estimated anthropogenic heat lacks some detail contrasts with multiple publicly datasets, which might be a method of verifying the estimated anthropogenic heat. Validation of the estimated heat fluxes would be more comprehensive if the we can be further validated with existing publicly available anthropogenic heat datasets [48,49].

Conclusions
In this study, the central urban area of Guangzhou City was taken as an example to estimate the temporal and spatial variation of anthropogenic heat in different urbanization periods. The energy balance model was used based on remote sensing images. The transfer matrix and the migration of the gravity center was used to analyze the change characteristics of the anthropogenic heat and the drivers of anthropogenic heat change from population, GDP, and industry aspects were also explored. The main conclusions are as follows: The overall change trend of anthropogenic heat in the central urban area of Guangzhou was enhanced during the period from 2004 to 2020, and the degree of its enhancement was related to the types of urban functional land. The central area of Guangzhou is in a stage of rapid urbanization, which has stimulated a substantial increase in energy demand and supply, and the intensity of anthropogenic heat has also increased. The increase in the intensity of anthropogenic heat in industrial areas and airport areas is greater than that of commercial and residential areas.
There are obvious differences both in area and space for different types of anthropogenic heat with urban development. In terms of area expansion, the areas of the zero-AHF zone, low-AHF zone, and medium-AHF zone changed drastically, with the zero-AHF zone and low-AHF zone transferring out and medium-AHF zone transferring in as the main form, while the area change of the high-AHF zone was relatively stable. The urbanization of Guangzhou's central urban area has led to an increase in the demand for urban energy consumption, which is also the main reason for the unbalanced state of the single conversion of anthropogenic heat from low-AHF areas to high-AHF areas. Spatially, the distance change of the low-AHF zone is relatively stable; the migration distance is 1.31 km 2 . The distance changes of the medium-AHF zone and the high-AHF zone are more intense, with migration distances of 13.19 and 19.06 km, respectively. Some large-scale industrial areas were identified in the study area in 2020 with the improvement of transportation equipment, while industrial areas were mainly distributed in port areas in 2004. This is the main reason for the large changes of space in the medium-and high-AHF areas.
The expansion and enhancement of the influence of human activities, such as the increase in urban population, rapid economic development, and increased industrial production activities, have stimulated the emission of anthropogenic heat, which contributes to the positive impact on the intensity of anthropogenic heat and the occurrence of significant urban anthropogenic heat change. Analysis of the characteristics of temporal and spatial changes in anthropogenic heat and studies on the laws and influencing factors of anthropogenic heat changes will provide references for urban planning and construction in mitigating urban heat islands.

Data Availability Statement:
The data presented in this study are available from author upon reasonable request.