Spatial-Temporal Evolution and Relationship between Urbanization Level and Ecosystem Service from a Dual-Scale Perspective: A Case Study of the Pearl River Delta Urban Agglomeration

: Since the beginning of the 21st century, urbanization has brought about dramatic changes in population, life, and economy, while having a signiﬁcant impact on the distribution of ecosystem service. As research on the relationship between urbanization and ecosystem service has gradually speciﬁed, we decided to explore it at different scales. In this paper, we quantiﬁed and mapped the spatial–temporal evolution and relationship between urbanization and ecosystem service value in the Pearl River Delta Urban Agglomeration from 2000 to 2019 based on a dual-scale perspective of county and 3 km × 3 km raster. Our results show that the overall trend of urbanization level and ecosystem service value was increasing. Urbanization and ecosystem service value at the county scale showed a negative spatial correlation, while it was not signiﬁcant at the raster scale. The “high–high” agglomeration was more concentrated, while the other three “low–low”, “low–high” and “high–low” agglomerations were more dispersed. Our ﬁndings suggest it is crucial to identify the key factors of small urban areas to grasp the development mechanism in the urbanization process and maintain the balance of the ecosystem.


Introduction
Since the beginning of the 21st century, urbanization has been accelerating, which has led to ecological damage and other problems while promoting rapid development of society, economy, life, and culture. Urbanization in a narrow sense refers to the transfer of rural population to urban areas; urbanization in a broad sense includes the urbanization process of population, economy, society, natural ecology and land [1], and the urbanization discussed in this paper focuses on the broad sense. Ecosystem Service (ES) is the direct or indirect benefits that humans obtain from ecosystem [2,3], including tangible and intangible services that contribute to human survival and quality of life [3], which is important for maintaining urbanization development. After studying the impact of urbanization on the structure and function of wetland forests, Faulkne [4] found that urbanization development changes regional ecosystem service functions. Kreuter [5] obtained from a coupled perspective that urbanization changes the ecosystem service by changing land cover and natural land use. Therefore, it is important to understand the spatial-temporal relationship between urbanization and ecosystem service to coordinate regional sustainable development.
At present, the research on quantitative assessment of ecosystem service value (ESV) is improving, mainly including the comprehensive indicator assessment method [6], the value quantity method [7] and the energy value analysis method [8], among which the an area of about 56,500 km 2 . Since the reform and opening up, the urbanization of the PRD has developed rapidly, with the urbanization rate of core area reaching 86.28%, per capita GNP of 160,000 yuan and a population density of 1212 people/km 2 at the end of 2019. As the most densely populated urban agglomeration in the country, the ecosystem of PRD has been fragmented by urbanization and the differences in population density and socioeconomic development have exacerbated the uneven development of the region, presenting significant spatial differences. The coordination of urbanization and ecosystem service has become increasingly prominent.
Three types of spatial and statistical data are used in this study. of the Pearl River, it belongs to the southern subtropical maritime monsoon climate and covers an area of about 56,500 km 2 . Since the reform and opening up, the urbanization of the PRD has developed rapidly, with the urbanization rate of core area reaching 86.28%, per capita GNP of 160,000 yuan and a population density of 1212 people/km 2 at the end of 2019. As the most densely populated urban agglomeration in the country, the ecosystem of PRD has been fragmented by urbanization and the differences in population density and socioeconomic development have exacerbated the uneven development of the region, presenting significant spatial differences. The coordination of urbanization and ecosystem service has become increasingly prominent.
Three types of spatial and statistical data are used in this study. First, land use/cover data for five years of 2000, 2005, 2010, 2015 and 2019 include cropland, forest land, grassland, water body, construction land and unused land with resolution of 30 m × 30 m, which are obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). Secondly, demographic, social and economic data are obtained from Guangdong Provincial Statistical Yearbook, China County Statistical Yearbook and related city and county statistic yearbook. Finally, the third type of data is the nighttime lighting data from the dataset of An extended time-series (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) of global NPP-VIIRS-like nighttime light data which were obtained from Harvard Dataverse Platform (https://dataverse.harvard.edu/dataset.xhtml?persisten-tId=doi: 10.7910/DVN/YGIVCD) and has a ground resolution of 30 arc second. Considering the feature of China, all data are projected to WGS_1984_UTM_Zone_49N, clipped by the shape file of Pearl River Delta Urban Agglomeration and created to 3 km × 3 km raster with ArcGIS 10.2 to accumulate the level of urbanization and ecosystem service value. See Figure 1.

Caculation of Value of Ecosystem Service
The quantitative assessment of ecosystem service value is important for maintaining regional ecological security and promoting the coordinated development of regional economy and environment [27]. In this paper, the value of ecosystem service in the PRD Note: The map is based on the standard map (review number GS (2016 2556)) downloaded from the standard map service website of the National Bureau of Surveying, Mapping and Geographic Information, with no modifications to the base map.

Caculation of Value of Ecosystem Service
The quantitative assessment of ecosystem service value is important for maintaining regional ecological security and promoting the coordinated development of regional economy and environment [27]. In this paper, the value of ecosystem service in the PRD is measured by the improved unit area value equivalent factor method of Xie [28]. Among them, the value of one standard unit ESV equivalent factor is equivalent to 1/7 of the average economic value of grain yields in the PRD (Equation (1)), which is calculated as where: C = the value of one standard unit of ESV equivalent factor (yuan/hm 2 ); P = the average price of grain in the PRD; Q = the grain yield per unit area in the PRD (kg/hm 2 ).
To eliminate the interference of price fluctuations on value changes, this study introduces the Consumer Price Index (CPI) to adjust the average food price data of each year to the price level of 2019 [23,29], and uses the correction method of Xu [30] and Xie [31] for the equivalent factor table of the ecosystem per unit area. The economic value of the service function is corrected (Equations (2) and (3)) to derive the economic value of food production service provided per unit area of the farmland ecosystem, and finally the corrected ecosystem service value coefficient of the PRD. The calculation formula is: where: λ = the equivalent factor correction factor; Q = the average grain yield per unit area in the PRD from 2000 to 2019 (kg/hm 2 ); Q 0 = the national annual average grain yield per unit area (kg/hm 2 ); θ = the biomass factor of farmland ecosystems in different provinces in China proposed by Xie Gaodi; E i = the modified equivalent factor of land use type i; E 0 = the equivalent factor of land use type determined by Xie. The specific process is as follows: The average grain production in the PRD from 2000-2019 is 5291.3 kg/hm 2 according to the Guangdong Provincial Statistical Yearbook. Based on the grain price statistics from the South China Grain Exchange, the 2019 PRD grain purchase price was estimated to be 2.55 yuan/kg. Using the CPI to adjust the average grain price data for each year to the 2019 price level, the value of 1 standard unit of ESV equivalent factor in the PRD in 2000, 2005, 2010, 2015 and 2019 was calculated 1278. 28, 1369.8, 1873.29, 1928.07 and 1873.48 yuan/hm 2 , respectively. Based on the national agricultural land unit area grain yield of 5004.65 kg/hm 2 and Guangdong Province land use type equivalent factor of 1.4, the correction coefficient of the PRD ecological service value equivalent factor table was obtained as 1.48, and then the land use/cover of the PRD was linked to the ecosystem type, in which the ecological service value generated by construction land was not considered for the time being, arable land corresponded to agricultural land, forest land corresponds to forest, beaches and mudflats in water bodies correspond to wetlands, and the corresponding value of desert is taken for unused land and multiplied by the ecological service value per unit equivalent to obtain a table of ecosystem service value factors suitable for the study area (Table 1). It shows that ecosystem service value coefficient per unit area is on the rise during the period of 2000-2015 and then decrease. Water body contributes more than any other land use types, followed by wetland and woodland, with over 300 thousand, 100 thousand and 50 thousand yuan, respectively. It was found through the literature that demographic, economic, life, social, ecological and land urbanization are the main focus in urbanization evaluation studies. However, as ecological urbanization is more relevant to the ecosystem service involved in this study, the overlap between society and life is high but society does not cover as much as life, and land urbanization covers more aspects of life, economy, and population. To avoid contradictions in the later analysis and to consider the availability of data, this paper selects three aspects of life, economy and population to build a comprehensive urbanization level evaluation system for the Pearl River Delta Urban Agglomeration, which contains 8 specific indicators ( Table 2).  [26] + 0.24 1, 2 "+" represents a positive indicator and "−" a negative indicator.

Weighting Measurement Method
The entropy method is an objective assignment method, which calculates the index weights according to the dispersion degree of each index, the greater the dispersion degree of the index, the greater the influence of the index on the comprehensive evaluation, we adopt the entropy method to determine the index weights, and the calculation steps are as follows. 1 Standardization of original data.
Positive indicators: Negative indicators: x kij = the value of the jth indicator in the ith district and county in year k; 2 Given r (k = 1, 2, ..., r) counties, m (I = 1, 2, ..., m) years, and n (j = 1, 2, ..., n) indicators, calculate the weight y kij of the value of the jth indicator in the ith year of kth county: 3 Calculate the information entropy value e j for the jth indicator: 4 Calculation of indicator weights W j : 5 Calculate the index score F j : The index scores F j are calculated by Formula (9) to reflect the comprehensive urbanization level, life urbanization level, economy urbanization level and population urbanization level of each district or county, and the larger the value is, the higher the urbanization level is.

Spatialization of Comprehensive Urbanization Level
The level of urbanization is a comprehensive reflection of urbanization in terms of population, human activities, industry and economy, etc. The nighttime light data are a combined reflection of the results of the interaction of these factors, therefore, the integrated processed NPP-VIIRS-like nighttime light data are selected to correct the comprehensive urbanization level in this study, and the spatialization path of the comprehensive urbanization level is shown in Figure 2.

1
Modified construction of the light index reflecting the level of regional urbanization The lighting index L i is constructed from two attributes of regional night lighting distribution to correct the comprehensive urbanization level: (i) regional average lighting intensity; (ii) regional lighting area. The two attributes are defined and calculated by the following indexes, and then linearly weighted to form the lighting index (L i ), which is defined as a relative value to remove the influence of scale differences and to facilitate year-to-year comparisons.
where L i = the lighting index of region i; I i = the average light intensity indicator of region i, which characterizes the proportional relationship relative to the maximum possible light intensity; Sustainability 2021, 13, 8537 7 of 20 S i = lighting area index of region i, which is defined as the total area of all lighting areas in the region (light intensity is not zero area) accounted for the proportion of the entire area, reflecting the extensibility of the light in space; u, v (u + v = 1) are the weights of I i and S i , which will be determined by the correlation analysis later on; D j = the light intensity of the jth light area in region i; M j = the area of the jth light area in region i; D max = the maximum light intensity in the study area; ∑M j = the total area of all lighting areas in region i; A i = the total area of region i. 1. Modified construction of the light index reflecting the level of regional urbanization The lighting index is constructed from two attributes of regional night lighting distribution to correct the comprehensive urbanization level: (i) regional average lighting intensity; (ii) regional lighting area. The two attributes are defined and calculated by the following indexes, and then linearly weighted to form the lighting index ( ), which is defined as a relative value to remove the influence of scale differences and to facilitate year-to-year comparisons.
where = the lighting index of region i; = the average light intensity indicator of region i，which characterizes the proportional relationship relative to the maximum possible light intensity; = lighting area index of region i, which is defined as the total area of all lighting areas in the region (light intensity is not zero area) accounted for the proportion of the entire area, reflecting the extensibility of the light in space; , ( + = 1) are the weights of and , which will be determined by the correlation analysis later on;

Correlation analysis of the light index and the comprehensive urbanization level of cities
In this study, the correlations between different nighttime lighting indices and comprehensive urbanization levels were calculated to determine the optimal lighting indices for the sub-districts, and modeling regression analysis was performed using the econometric software STATA16.0. Experimentally, it was proved that the correlation and fitting effect was best after logarithmic treatment of the variables, so the regression model as in Equation (13) was established after logarithmic transformation of the variables first.
where: UR i = the comprehensive urbanization level of region i; (10); ε i = the random error term in the model; and a and b are the regression model coefficients.
Regression analysis of the weight combinations of different lighting indices yielded the following results for the fit (R 2 ). As seen in Table 3, separate regressions were performed from the (0.1, 0.9) combination of weights to the (1, 0) combination, and the results showed the best fit at the (0.9, 0.1) combination, with an R 2 of 0.81. The formula for the lighting index L i was thus determined as

Weighted Combinations
Note: All passed the significance test at the 0.01 level.

Spatialization of comprehensive urbanization level
The comprehensive urbanization level of the raster is assigned to the corresponding raster after the coefficients are calculated based on the regression. Finally, the actual comprehensive urbanization level of each district and county is linearly adjusted with the predicted comprehensive urbanization level to correct the value of each raster, and the corrected comprehensive urbanization level of the PRD at the raster scale is generated with the following equations.
Equation (15) is transformed from Equation (13), where: UR c i = the comprehensive urbanization level of the raster after correcting for the actual comprehensive urbanization level; UR i = the predicted comprehensive urbanization level of each raster; UR i = the comprehensive urbanization level of the district and county to which the raster belongs; UR * i = the predicted comprehensive urbanization level of the district and county to which the raster belongs.

Pearson Correlation Analysis
Analyzing the correlation between single and comprehensive urbanization level and ecosystem service value can better reflect the mutual influence of urbanization and ecosystem service. In this study, the Spearman coefficients were calculated by STATA to identify the correlation between the comprehensive urbanization level (Compr_ur), population urbanization level (Pop_ur), economy urbanization level (Eco_ur), life urbanization level (Life_ur) and ecosystem service value (ESV) from 2000 to 2019, respectively. It can also be used to analyze the sensitivity of the correlation and explore the extent of its impact.

Bivariate Spatial Autocorrelation Analysis
Bivariate spatial autocorrelation analysis is a measure of whether multiple variables are clustered in a spatial distribution [35] and explains the correlation between the value of spatial unit attributes and the value of other attributes in its adjacent space, including two parts: global spatial autocorrelation and local spatial autocorrelation [25]. The global spatial autocorrelation is used to describe the degree and significance of spatial correlation between variables in the study area, and the local spatial autocorrelation can discern the aggregation effect of an element in space on the surrounding elements to achieve the study of the dynamics of local changes in space [36], and the calculation equations [37] are: where I = the global Moran's I index; I i = the local Moran's I index; n = the number of spatial cells; x i , x j = the observed values of cell i and cell j, respectively; (x i − x) = the deviation of the observed value on the ith spatial cell from the mean, W ij = the spatial weight value between cell i and cell j; and S 2 is the variance.
In this paper, we used GeoDa1.14, a spatial analysis software, to establish a spatial distance weight matrix to analyze the spatial autocorrelation between comprehensive urbanization level and ecosystem service value in the PRD, and used LISA (Local indicators of Spatial association) distribution maps to explore the spatial correlation between the two and their spatial heterogeneity characteristics. Since 2005, the highest growth rate value of each district and county has gradually decreased, and the difference has become smaller and smaller, especially between the highest and lowest values of its growth rate between year 2015 and 2019 is less than 20%, indicating that the overall trend tends to urbanization maturity. −11.24% in Longhua District during 2015-2019. Since 2005, the highest growth rate value of each district and county has gradually decreased, and the difference has become smaller and smaller, especially between the highest and lowest values of its growth rate between year 2015 and 2019 is less than 20%, indicating that the overall trend tends to urbanization maturity.

Spatial-Temporal Evolution of the Raster Scale
As seen in Figure 4:  Table 4.
x FOR PEER REVIEW 12 of 21

Spatial-Temporal Evolution of Urbanization Level at County Scale
As seen in Figures 5 and 6 (1) the comprehensive urbanization level of the PRD rose from 0.1532 to 0.5178 between 2000 and 2019, a 3.38-fold rise in rapid development over 20 years. In terms of individual urbanization levels, economy urbanization grew explosively from 0.0319 in 2000 to 0.3562 in 2019, an 11.2-fold increase. Within the same development stage, the population urbanization level increased from 0.0327 to 0.0563, and the life urbanization level increased from 0.0885 to 0.1053, showing a slow growth trend. Among them, the economy urbanization level of the PRD urban agglomeration is the fastest growing, and the comprehensive urbanization level ranks second, while the other two types of urbanization levels are relatively stable within the study period.

Spatial-Temporal Evolution of the Raster Scale
As seen in Figure 7: (1)

Spatial-Temporal Evolution of the Raster Scale
As seen in Figure 7: (1)

Correlation Analysis of Urbanization Level and ESV
As shown in Table 5: Ecosystem service value shows negative correlation with comprehensive urbanization, population urbanization and economy urbanization, and positive correlation with life urbanization, respectively, and all correlations are significant at the level of 0.01. The degree of influence of urbanization development on ESV is strong. For each 1% increase in comprehensive, population and economy urbanization levels, ESV will decrease by 43.9, 53.9 and 28.1 ten thousand yuan, respectively, with population urbanization level having the greatest impact; for each 1% increase in life urbanization, ESV increases by 32.6 ten thousand yuan. The high life urbanization level motivated by the rising energy efficiency, people's awareness of environment protection and advanced technology has promoted the environment protection which will greatly improve the ecosystem value. However, how specific factors affect ecosystem need more complex models to deduce, such as system dynamic, neural networks, etc. Table 5. Pearson relating analysis between urbanization level and ESV.

Spatial Autocorrelation Analysis of Comprehensive Urbanization Level and ESV
1. Spatial-temporal changes at the county scale As seen in Table 6 2) In 2000, the areas with "high-high" agglomeration of integrated urbanization level and ecosystem service value were Chancheng District, Pengjiang District and Jianghai District; the areas with "low-high" agglomeration were Deqing, Guangning County, Enping City, Nansha District and Huicheng District; the areas with "high-low" agglomeration were Tianhe District, Yuexiu District and Haizhu District of Guangzhou City and Nanshan District, Luohu District and Futian District of Shenzhen City. "High-Low" clustering areas are Tianhe District, Yuexiu District and Haizhu District of Guangzhou City and Nanshan District, Luohu District and Futian District of Shenzhen City, and "Low-Low" clustering does not reach 95% confidence level. In 2005, the "high-high" agglomeration area was Nansha District, the "low-high" agglomeration area increased by three-Pengjiang District, Jianghai District and Xiangzhou District-and decreased by one-Nansha District-the "high-low" agglomeration The overall change of "high-low" agglomeration is not significant, with the addition of Baiyun District and the reduction of Yantian District, and the "low-low" agglomeration still does not reach the 95% confidence level; in 2010, the "high-high" agglomeration did not reach the 95% confidence level. The overall number of "low-high" agglomeration areas remained unchanged, among which Pengjiang District and Jianghai District in 2005 were changed to Fengkai County and Kaiping County, and the number of "high-low" agglomeration areas increased by Baoan District and Longhua District and decreased by Baiyun District. The "low-low" agglomeration areas are Baiyun District and Yantian District. In 2015, the "high-high" agglomeration did not reach the 95% confidence level, and the "low-high" agglomeration expanded extensively. (3) From 2000 to 2005, the "low-low" agglomeration areas in the PRD did not reach the 95% confidence level, and the only significant agglomerations were "high-high," "low-high," and "high-low". However, since 2010, the "high-high" agglomerations have disappeared and have been replaced by "low-low" agglomerations, indicating that the level of urbanization in the PRD has been increasing. Intense human activities have seriously damaged the ecological environment and resources, resulting in a significant reduction in the value of regional ecosystem service. (4) The "high-low" agglomeration areas are mainly located in the city centers of Guangzhou and Shenzhen, where Tianhe, Yuexiu and Haizhu districts of Guangzhou and Nanshan and Futian districts of Shenzhen have been in the "high-low" agglomeration in the past 20 years. These areas are in the plains, with a high concentration of humans and a high proportion of economic production. The strong human activities have led to the continuous improvement of the infrastructure construction in population, economy and life, and the increase of the comprehensive urbanization level, while destroying the local green space and water environment to a large extent, resulting in low ecosystem service value.(5) The "low-high" gathering areas are mainly located in Guangning County and Deqing County in the northwest of Zhaoqing City, Enping County and Kaiping County in Jiangmen City in the southwest, and Huicheng District in Huizhou City in the northeast, etc. These areas are high, mountainous and rich in forest resources, which have obvious ecological advantages. Moreover, the topographic environment with many mountains and little land leads to a small distribution of local population and inconvenient transportation, which to a certain extent hinders the gathering of people and development of industries, and makes it difficult for urban construction land to expand outward, which inhibits the improvement of the comprehensive urbanization level and presents a spatial gathering pattern of low urbanization and high ecological service value.  2. Spatial-temporal variation at the raster scale As seen in Table 7 and Figure 9: (1) the global Moran's I index of comprehensive urbanization level and ecosystem service value in the PRD at raster scale was positive from 2000 to 2015, and the two were positively correlated. Places with high urbanization levels were more likely to gather high ecosystem service value, which changed to negative in 2019, and the two were negatively correlated, but their absolute values all converged to zero, indicating that the two spatially differ greatly and there is no significant spatial correlation. (2) From 2000 to 2005, the distribution of "high-high" agglomeration in the PRD was mainly concentrated in Foshan City, Duanzhou District in southeastern Zhaoqing City, Sanshui District and southeastern Sihui City, northwestern Zhongshan City, the

2.
Spatial-temporal variation at the raster scale As seen in Table 7 and Figure 9: (1) the global Moran's I index of comprehensive urbanization level and ecosystem service value in the PRD at raster scale was positive from 2000 to 2015, and the two were positively correlated. Places with high urbanization levels were more likely to gather high ecosystem service value, which changed to negative in 2019, and the two were negatively correlated, but their absolute values all converged to zero, indicating that the two spatially differ greatly and there is no significant spatial correlation.
(2) From 2000 to 2005, the distribution of "high-high" agglomeration in the PRD was mainly concentrated in Foshan City, Duanzhou District in southeastern Zhaoqing City, Sanshui District and southeastern Sihui City, northwestern Zhongshan City, the junction of western Dongguan City with Panyu District and Nansha District, and the junction of northeastern Jiangmen City with Foshan City. Since 2010, Foshan City changed to "high-low" agglomeration and the "high-high" agglomeration area spread from Foshan City to Sanshui District, Dinghu District, Duanzhou District and the southeast of Sihui City in the southeast of Zhaoqing City, and to the middle of Nansha District, Dongguan City and the junction of Panyu District and Nansha District in the southeast of Dongguan City. and Nansha District, southward to the junction of Jiangmen City and Foshan City, Jianghai District, Jinwan District, central and northeastern Xinhui District, and western and northern Zhongshan City, the proportion of "high-high" clusters in the past 20 years showed a trend of first increasing and then decreasing, with an overall relative increase, and the spatial pattern gradually growing from a small ram's horn layout to a large ram's horn.
(3) In the past 20 years, the spatial pattern of "low-low" agglomeration areas has changed less, mainly in the western part of Conghua District, Huizhou City, the northeastern part of Boluo County, the eastern and southern part of Huidong County, the western part of Huaiji County, Zhaoqing City, the northern and southern part of Enping City, Jiangmen. (4) The spatial pattern and numerical changes of the "low-high" agglomeration areas in the past 20 years are relatively small, mostly concentrated in Sanshui District, Dinghu District, Xinhui District, Doumen District, Jinwan District and other "high-high" agglomeration areas and the southern part of Taishan City and the eastern part of Huidong County. (5) From 2000 to 2005, the "high-low" agglomeration areas in the PRD were concentrated in Guangzhou, Shenzhen, central Dongguan and southwestern Conghua District, and from 2010, the "high-low" agglomeration areas added Foshan City on top of the previous ones, except for a slight expansion north of Shenzhen and a significant shift from the original "low-low" to "high-low" in northeastern Boluo County, and the proportion of agglomeration areas continued to increase in the past 20 years. Table 7. Percentage of each cluster in the Pearl River Delta in 2000-2019.

Percentage of Each Cluster/% Global Moran's I Not Significant
High-High Low-Low Low-High High-Low Neighborless

Calculation of the Ecosystem Service Value in the PRD
In this paper, the ESV of the PRD urban agglomeration is assessed using a modified unit area value equivalent factor method based on correction coefficients, which is easy to operate compared to market value, shadow engineering and alternative cost methods, and is applicable to regional, national and global scales [38]. In recent years, more studies on the PRD have accounted for ESV with the help of this method; Xiao [39] assessed the ESV of the PRD from 2000 to 2015 with Guangzhou-Foshan-Zhaoqing as an example, which were 1402.68, 1440.98, 1413.88, and 139.909 billion yuan, showing a trend of first increasing and then decreasing year by year. The assessment results after area conversion are similar to ours and the change trend is consistent; Xu [40] assessed the ESV of the PRD urban agglomeration from 2000 to 2015 to be about 300 billion yuan, and the assessment results are also similar, so the conclusions of this study are reliable.

Calculation of the Ecosystem Service Value in the PRD
In this paper, the ESV of the PRD urban agglomeration is assessed using a modified unit area value equivalent factor method based on correction coefficients, which is easy to operate compared to market value, shadow engineering and alternative cost methods, and is applicable to regional, national and global scales [38]. In recent years, more studies on the PRD have accounted for ESV with the help of this method; Xiao [39] assessed the ESV of the PRD from 2000 to 2015 with Guangzhou-Foshan-Zhaoqing as an example, which were 1402.68, 1440.98, 1413.88, and 139.909 billion yuan, showing a trend of first increasing and then decreasing year by year. The assessment results after area conversion are similar to ours and the change trend is consistent; Xu [40] assessed the ESV of the PRD urban agglomeration from 2000 to 2015 to be about 300 billion yuan, and the assessment results are also similar, so the conclusions of this study are reliable.

Calculation of Urbanization Level in the Pearl River Delta
There are many methods for calculating the urbanization level [18,32], and this paper draws on existing studies to characterize the comprehensive urbanization level of the region by selecting indicators from three aspects: population, life, and economy. On the one hand, the dimensions affected by the urbanization process are far richer than the urbanization level system constructed in this study, and it is more difficult to provide an integrated description of the urbanization level due to data limitations. On the other hand, no scholars have yet spatialized the comprehensive urbanization level using nighttime lighting data, and more scientific and accurate and realistic adjustments will be made in the future through relevant literature studies and parameter corrections.

Calculation of Urbanization Level in the Pearl River Delta
There are many methods for calculating the urbanization level [18,32], and this paper draws on existing studies to characterize the comprehensive urbanization level of the region by selecting indicators from three aspects: population, life, and economy. On the one hand, the dimensions affected by the urbanization process are far richer than the urbanization level system constructed in this study, and it is more difficult to provide an integrated description of the urbanization level due to data limitations. On the other hand, no scholars have yet spatialized the comprehensive urbanization level using nighttime lighting data, and more scientific and accurate and realistic adjustments will be made in the future through relevant literature studies and parameter corrections.

Spatial-Temporal Evolution Patterns of Urbanization Level and ESV
By analyzing the spatial-temporal evolution of urbanization level and ESV in the PRD from 2000 to 2019, it is found that: high-value areas with comprehensive urbanization level are concentrated in the middle of the PDR, mainly in Guangzhou and Shenzhen, with a large expansion of high-value area starting in 2005 in particular. Guangzhou and Shenzhen, as the first cities to develop, have driven the development of the surrounding areas, and the urbanization spillover effect is obvious, showing a mature urbanization pattern in the region. The total value of ecosystem service shows a trend of first increasing and then decreasing, and the value added mainly comes from the ecological value contribution of water bodies, forest land and other land resources, while the decrease of ESV year by year after 2010 reflects the negative effect of urbanization changing land use cover has emerged, which means the expansion of building land has caused more serious damage to ecosystem service. ESV from the spatial distribution of high in the belly and low around the PRD in 2000 to low in Guangzhou, Shenzhen and Foshan and high in Foshan with northwest, southeast and south in 2019, with enhanced spatial heterogeneity and uneven distribution. Therefore, it is necessary not only to increase the overall efforts to protect natural resources in the future, but also to pay more attention to the spatial equalization of ecological service value, especially to adjust and optimize the proportion of land types in areas with high urbanization levels.

Spatial Relationship between Urbanization Level and ESV
The spatial interaction between urbanization level and ESV generally has a spillover effect, and a large number of studies have shown that they are negatively correlated [26,41], and by studying the spatial relationship between the comprehensive urbanization level and ESV in the PRD from 2000-2019, it was found that there is a negative correlation between the two at the county scale, and there is no significant spatial correlation at the raster scale, indicating that the spatial The spatial correlation between the two was found to be negative at the county scale, but not at the raster scale, indicating that the spatial correlation between the two was different at different scales, with obvious scaledependent characteristics [37].At the spatial level, the raster scale shows obvious spatial heterogeneity compared with the county scale, and the "high-high" agglomeration was mostly concentrated in and around Foshan City from 2000 to 2005, but from 2010 onward, the region was influenced by the urbanization level of Guangzhou City and changed to "high-low" agglomeration, which was mainly distributed in "As one of the cities with the most drastic transformation in the past 20 years, Foshan City should focus on ecological environmental protection during its urbanization development, avoiding the negative impact of the spillover effects of arable land operation and urban construction and development on the natural ecology, and at the same time, it is necessary to set up as many eco-friendly land types in the natural areas at the edge of the city as possible [42], for instance add water body and wetland in the edge of Foshan City

Conclusions
This study measured the urbanization level and ecosystem service value of the PRD Urban Agglomeration from 2000-2019 from a dual-scale perspective, and analyzed their spatial-temporal evolution and spatial relationships, which are important for fine-grained analysis of the development of the region. The results of the study indicate that:

•
In the past 20 years, the overall trend of urbanization level and ecosystem service value in the PRD has been on the rise. Among them, the regional differences in urbanization level are significant, and the urbanization spillover effect in the central region drives the development of neighboring regions with centripetal nature; the spatial heterogeneity of ecosystem service value is strong, with the largest change between 2005 and 2010 and a gradual decline after 2010. • At the numerical level, both single and comprehensive urbanization are significantly correlated with ESV. ESV is negatively correlated with comprehensive, population, and economy urbanization levels, and positively correlated with life urbanization. • At the spatial level, comprehensive urbanization and ESV in the PRD from 2000 to 2019 are spatially negatively correlated at the county scale, while the spatial correlation is not significant at the raster scale. From 2000 to 2005, the distribution of "high-high" agglomeration is mainly concentrated in and around Foshan City; from 2010, the "high-low" agglomeration area expands from Tianhe District to Foshan City and occupies Foshan; from 2015 to 2019, the "high-high" and "high-low" agglomeration distribution is dominant. • Foshan City is the area that has changed most dramatically in the past 20 years, from a "high-high" agglomeration to a "high-low" agglomeration, and the rapid urbanization has damaged to Foshan's waters and green areas. Tianhe, Yuexiu District in Guangzhou City and Nanshan District in Shenzhen City, have always maintained the "high-low" clustering, indicating that the construction land has encroached on the local natural land cover and the ecological resource environment is poor.
At the same time, there are shortcomings in this study. First, with the deepening of the concept of "new urbanization", the indicators used in this study cannot accurately reflect the new urbanization level of the region, and a more combined and complete measurement system is needed. Secondly, as the factors influencing urban agglomerations are too complex, covering population, economic and living conditions as well as ecology, smaller urban areas will be focused on for specific exploration. How to identify the key factors of smaller urban areas, reveal their development mechanisms and maintain a balanced ecosystem therefore proves to be necessary in the urbanization process. Data Availability Statement: New data generated is shared through this article. All other sources of data are cited throughout the paper.