Sustainable Land Urbanization and Ecological Carrying Capacity: A Spatially Explicit Perspective

: Rapid urbanization has become a common occurrence all over the world, particularly in developing countries, and has thus resulted in various eco-environmental problems. In China, urban land has expanded at an unprecedented rate in the past several decades, and sustainable land urbanization has become an important issue in promoting sustainable development. Hence, scholars have proposed ecological carrying capacity (ECC) as a solution to balance socio-economic development and the ecosystems for achieving sustainable development. In the current work, we explored the spatial inﬂuence of ECC on land urbanization and its driving mechanism, using the Wuhan agglomeration as a case study. In the ﬁrst step, we calculated the ECC at the county level using the ecological footprint method. Then, we applied a combination of kernel density and the “densi-graph method” on the basis of points of interest, in order to identify urbanized areas and to measure land urbanization rates. Finally, we devised spatial models with ECC-based spatial weight matrices to examine the potential spatial interactions or constraints and the inﬂuencing factors. Results indicate the following. (1) Land urbanization rates in most counties increased, whereas the average ECC per capita in the Wuhan urban agglomeration decreased from 2010 to 2015; (2) China’s land urbanization is primarily driven by socio-economic development, in which ﬁxed asset investments and urban income present positive inﬂuences and agricultural outputs show a negative inﬂuence; (3) Spatial interaction was formulated through ECC during the land urbanization process. However, this effect was attenuated in 2010–2015. The ﬁndings are beneﬁcial for understanding the regional spatial inﬂuence of ECC on urban land urbanization. They should also facilitate the formulation of relevant policies for protecting, restoring, and promoting the sustainable use of terrestrial ecosystems to ultimately achieve coordinated and balanced regional development.


Introduction
Urbanization is the process of expansion in terms of the urban population and urban scale with relevant economic and social changes [1]. The urban population of the world increased rapidly from 751 million in 1950 to 4.2 billion in 2018 [2]. Tremendous demographic migration and socio-economic The remainder of this paper is organized as follows. Section 2 constructs an ECC-based spatial regression model to gauge the spatial interactions from ECC in land urbanization. The "densi-graph" method is employed to derive urbanized land using the Wuhan urban agglomeration as an example. Section 3 provides the results. Section 4 presents the discussion. Section 5 draws the conclusions.

Data and Study Area
The Wuhan urban agglomeration is situated in the eastern part of Hubei Province (upper left longitude 112 • 30 E and latitude 29 • 05 N, lower right longitude 116 • 10 E and latitude 31 • 50 N) [27]. The Wuhan urban agglomeration includes the capital of Hubei Province (i.e., Wuhan) and eight other prefecture-level cities, namely, Huangshi, Ezhou, Xiaogan, Xianning, Xiantao, Huanggang, Qianjiang, and Tianmen ( Figure 1). This agglomeration is composed of 48 counties and covers a total area of approximately 58,000 km 2 , which is less than one-third the size of Hubei Province [28]. Nevertheless, the Wuhan urban agglomeration accounts for more than 50% of the population and gross domestic product (GDP) of Hubei Province [29]. The Wuhan urban agglomeration includes Jianghan Plain, and the Yangtze River runs from west to east. The water area is 6344 km 2 , accounting for 11% of the total area in 2015 according to Landsat TM images. The Wuhan urban agglomeration was selected in this study because it is one of the largest and fastest growing urban agglomerations in China with a strategic position [28]. Similar to its central city, namely, Wuhan, the theme of which is "Wuhan, Different Every Day", the Wuhan urban agglomeration is highly representative of urbanization [27]. The remainder of this paper is organized as follows. Section 2 constructs an ECC-based spatial regression model to gauge the spatial interactions from ECC in land urbanization. The "densi-graph" method is employed to derive urbanized land using the Wuhan urban agglomeration as an example. Section 3 provides the results. Section 4 presents the discussion. Section 5 draws the conclusions.

Data and Study Area
The Wuhan urban agglomeration is situated in the eastern part of Hubei Province (upper left longitude 112°30′ E and latitude 29°05′ N, lower right longitude 116°10′ E and latitude 31°50′ N) [27]. The Wuhan urban agglomeration includes the capital of Hubei Province (i.e., Wuhan) and eight other prefecture-level cities, namely, Huangshi, Ezhou, Xiaogan, Xianning, Xiantao, Huanggang, Qianjiang, and Tianmen ( Figure 1). This agglomeration is composed of 48 counties and covers a total area of approximately 58,000 km 2 , which is less than one-third the size of Hubei Province [28]. Nevertheless, the Wuhan urban agglomeration accounts for more than 50% of the population and gross domestic product (GDP) of Hubei Province [29]. The Wuhan urban agglomeration includes Jianghan Plain, and the Yangtze River runs from west to east. The water area is 6344 km 2 , accounting for 11% of the total area in 2015 according to Landsat TM images. The Wuhan urban agglomeration was selected in this study because it is one of the largest and fastest growing urban agglomerations in China with a strategic position [28]. Similar to its central city, namely, Wuhan, the theme of which is "Wuhan, Different Every Day", the Wuhan urban agglomeration is highly representative of urbanization [27].  The data used included interpreted land use maps from Landsat TM images in 2010 and 2015, with a spatial resolution of 30 m. Thematic land use maps ( Figure 1) with eight categories, namely, arable land, forest, grassland, water area, urban area, rural area, other construction land, and unused land, were generated on the basis of the National Land Use Classification System (GB/T21010-2007) and the image processing results. The point of interest (POI) map of China was provided by Baidu open platform (http://lbsyun.baidu.com/index.php?title=lbscloud/poitags), and the POI for the Wuhan urban agglomeration was extracted ( Figure 1). The income, population, GDP, revenue, agricultural output, consumption, and fixed asset investment of all the counties were collected from the Wuhan City Yearbook (2010-2016), Hubei Yearbook (2010-2016), and Chinese City Yearbook (2010-2016). Not all the indices were available for the 48 counties in the agglomeration area. Thus, regressions and predictions were used as substitutes for the missing data.

Calculation of ECC
ECC can be calculated using various methods, with the ecological footprint being the most widely applied one [30,31]. The ecological footprint analysis method proposed by Rees and Wachernagel in 1992 is used to calculate ECC according to biological physical quantity. This method balances the imperfections of single carrying capacity models for resources or the environment and helps obtain necessary data [32]. Owing to its straightforwardness and effectiveness, this method has been employed by many scholars worldwide to evaluate ecological environmental quality and sustainable development [32][33][34][35][36]. Furthermore, most studies transformed the calculation of ECC into measurements of ecologically productive land areas through the quantification of a set of indicators [37,38], which have been widely used to assess ECC at different levels [39,40]. Specifically, the calculation equation for ECC is as follows: where ECC is the total ECC, N represents the population, ecc is the per capita ECC, 12% represents 12% of the land area deducted from the total area since the beginning of biodiversity conservation, j is the type of land use, a j is the ecological productive land area per capita, r j is the equilibrium factor, and y j is the yield factor. The data of the equilibrium factors were based on the "Living Planet Report" published by the World Wildlife Fund [41], and the data of the yield factors were adjusted on the basis of the output level in the study area [42] (Table 1). Urban development leads to urban land expansion, sprawl, and regeneration, thereby providing a basis for the formation of land urbanization. The boundary of an urbanized land area is defined as the urban spatial extent that has actually been constructed and developed. It has the essential urban infrastructure and facilities to meet the needs of daily life and production activities of residents.
Several studies used data from arbitrarily delineated areas obtained by the government, land surveying, or land use classification from remote sensing to represent urban land areas for calculating land urbanization [43]. However, this approach may induce bias or limitation because urbanized land is neither defined by the government nor confined to its superficial or spectral attributes. It is closely related to human activities and the physical attachments on the land. In the intersection of urban and rural areas, suburban areas, or urban fringes, the socio-economic development is high level. Moreover, anthropogenic activities are closely linked to the city. In this sense, these intersection area should also be regarded as part of the urbanized area. The urban-rural boundary threshold is a key element in determining the boundary of urban built-up areas, whereas the threshold in the existing literature often depends on human judgment [44,45].
In fact, a new method that identifies the boundaries of urbanized areas on the basis of the kernel density distribution of the POI data has become popular [24,25]. POI data refer to points, polylines, or polygons that describe the spatial distribution of primary infrastructure and facilities, such as shopping centers, hospitals, schools, and restaurants. POIs have various classifications according to regional characteristics and research objectives [46]. We classified POI data by considering the POIs that are used to derive urbanized land as follows ( Table 2). Green space and square land Parks, squares, and scenic spots Note: These data are based on the "Code for classification of urban land use and planning standards of development" (2011).
Kernel density and the densi-graph method are combined to obtain the boundary of the urbanized area. Kernel density estimation enables high-quality density estimation of point data without being affected by grid size and grid position [47], and it has been widely used in spatial analysis, especially in identifying urban-rural fringes [48,49]. The spatial distribution differences of POI density reflect the varying regional development levels, and the density in urbanized areas is generally higher than that in rural areas. Specifically, the estimated density of each point is the weighted average density of all points in the area. The kernel density P i of an arbitrary point i in space is defined as the highest density and exhibits an outward decreasing function at the center (2): where K j is the weight of the target j, D ij is the distance between neighboring point i and target j, R is the bandwidth of the selected regular region (D ij < R), and n is the range of the bandwidth R within the quantity of study object j. The selection of bandwidth R exerts a remarkable effect on the results of kernel density analysis [50], and a certain amount of bandwidth keeps the density center stable [51]. R is 2000 in 2015 and 1000 in 2010 according to the estimation of POI density distribution in different years ( Figure 2). Then, the densi-graph method is proposed to identify the actual urbanized areas using the contour map of the kernel density of the POI [25].The basic idea of the densi-graph is to obtain the urbanized area through the discovery of the point with abrupt change through the contour map generated by the kernel density estimation of POIs. In the first step, we produced the contour lines with consecutively increasing radii according to the kernel density estimation of POIs in each county. Then, we calculated the closed areas S d and the corresponding theoretical radius √ S d at each density point in each county. Subsequently, we produced the densi-graph with the density value d as the x-axis and the theoretical radius increment ∆ √ S d as the y-axis. Finally, we induced the critical point in the densi-graph, which showed the density curve undergoing an abrupt increase through differentiation and normalization; this curve was regarded as the threshold to generate the boundary of the urbanized and non-urbanized areas. We used the urbanized land percentage (ULP) to represent land urbanization (3).
where ULP i is the urbanized area in the ith county, UL i represents the urban land area of the ith county, and S i represents the land area of the ith county. induced the critical point in the densi-graph, which showed the density curve undergoing an abrupt increase through differentiation and normalization; this curve was regarded as the threshold to generate the boundary of the urbanized and non-urbanized areas. We used the urbanized land percentage (ULP) to represent land urbanization (3).
where ULPi is the urbanized area in the ith county, ULi represents the urban land area of the ith county, and Si represents the land area of the ith county.

Driving Factors and Spatial Interactions in Land Urbanization
The most widely accepted driving factors of land urbanization, urban sprawl, and urban land expansion can be categorized into socioeconomic, proximity, physical, accessibility, and neighborhood aspects on the basis of the review of existing studies [52]. By considering empirical studies, regional characteristics, and data accessibility, we selected nine potential factors, namely, disposable personal income for rural residents (DPIR), disposable personal income for urban residents (DPIU), proportion of contribution of second industries to GDP (PSI), proportion of contribution of tertiary industries to GDP (PTI), fixed asset investment per area (PFAI), government revenue per area (PGV), agricultural output value per capita (PAU), foreign trade export per area (PFT), and retail sales of consumer goods per capita (PRSC) [53][54][55]. We used "per capita" or "per area" to transform these explanatory variables into ratio form for comparison. Then, we performed correlation and regression analyses to find highly correlated factors and eliminate multicollinearity.
Spatial dependence and influence, which reflect the necessity to incorporate spatial effects into operational models, are widely considered in statistical calculation and regional science [56]. Spatial regression has been proven to be a suitable method for spatially identifying the urban expansion mechanism [57]. This approach provides a theoretical basis and spatial analysis technique for exploring underlying driving forces. In our study, the spatial interactions among ecological counties were represented by the mutual transfer process of urban ECC [58]. The spatial interactions among counties decline due to the distance among them, as described by Newton's law of universal gravitation. The gravity between two ecological counties can be described as follows Equation (4): where Eij is the ecological gravity between ecological counties i and j, ki and kj represent the ECC per capita that the two ecological counties can provide, dij is the distance between the two ecological counties, and r is the constant for the interaction of the ecological counties and is r = 1 in this study [58].

Driving Factors and Spatial Interactions in Land Urbanization
The most widely accepted driving factors of land urbanization, urban sprawl, and urban land expansion can be categorized into socioeconomic, proximity, physical, accessibility, and neighborhood aspects on the basis of the review of existing studies [52]. By considering empirical studies, regional characteristics, and data accessibility, we selected nine potential factors, namely, disposable personal income for rural residents (DPIR), disposable personal income for urban residents (DPIU), proportion of contribution of second industries to GDP (PSI), proportion of contribution of tertiary industries to GDP (PTI), fixed asset investment per area (PFAI), government revenue per area (PGV), agricultural output value per capita (PAU), foreign trade export per area (PFT), and retail sales of consumer goods per capita (PRSC) [53][54][55]. We used "per capita" or "per area" to transform these explanatory variables into ratio form for comparison. Then, we performed correlation and regression analyses to find highly correlated factors and eliminate multicollinearity.
Spatial dependence and influence, which reflect the necessity to incorporate spatial effects into operational models, are widely considered in statistical calculation and regional science [56]. Spatial regression has been proven to be a suitable method for spatially identifying the urban expansion mechanism [57]. This approach provides a theoretical basis and spatial analysis technique for exploring underlying driving forces. In our study, the spatial interactions among ecological counties were represented by the mutual transfer process of urban ECC [58]. The spatial interactions among counties decline due to the distance among them, as described by Newton's law of universal gravitation. The gravity between two ecological counties can be described as follows Equation (4): where E ij is the ecological gravity between ecological counties i and j, k i and k j represent the ECC per capita that the two ecological counties can provide, d ij is the distance between the two ecological counties, and r is the constant for the interaction of the ecological counties and is r = 1 in this study [58].
Spatial autocorrelation analysis is also performed to test whether the observations of a variable with a spatial position are closely associated with the observations at adjacent spatial points. Moran's I is a commonly used indicator that can be divided into two local spatial autocorrelation aspects [59]. We used global spatial autocorrelation to measure the degree of urban decentralization, which is defined as Equation (5): where n represents the number of observations, which is 48 in the Wuhan urban agglomeration; x i is the value of ULP in the ith county; x j is the value of ULP for county y j ; x is the mean value; and w ij is the value of the spatial weight matrix [60].    Figure 2 illustrates the kernel estimation of bandwidths with 500, 2000, and 4000 m in 2015. The kernel density distribution is too scattered in the 500 m bandwidth and is too aggregated in the 4000 m bandwidth. As a result, we selected 2000 m as the appropriate bandwidth to estimate the kernel density and produce a contour line map, as presented in Figure 3a. The densi-graph curve was generated (Figure 3b) to identify the point of abrupt change. Using the combined kernel density estimation and the densi-graph method in the 48 counties in the Wuhan urban agglomeration, we determined that the captured differentiated and normalized radius values with abrupt changes were 20% in 2010 and 5% in 2015. Twelve counties had a threshold of 2% in 2010, and nine counties had a threshold value of 5‰ in 2015.

Land Urbanization in the Wuhan Urban Agglomeration
We also compared the results of the urbanized land area extracted using the densi-graph method and the urban land area obtained from remote sensing images. Figure 4 further illustrates the results of urban land extraction through remote sensing images and the results of urbanized land area obtained through the combined kernel density and densi-graph method. The urbanized land area is slightly higher than the urban area because we included suburban areas or urban fringes with considerable urban infrastructure and great development potential. In this sense, the measurement of land urbanization is accurate given the focus on the "real" construction and terrestrial attributes of urbanized land. In terms of ULP, 45 counties exhibited an evident increase from 2010 to 2015, and the highest increment (24.64%) was recorded in Jiangan District. The land urbanization degrees in Qingshan District, Tieshan District, and Xisaishan District declined, but the decrease was less than 10%. The highest ULP values were observed in the central and western areas, particularly in the urban district of Wuhan. By contrast, the counties in the south, southeastern, and northeastern areas presented the least urbanization values. Qingshan District, Tieshan District, and Xisaishan District declined, but the decrease was less than 10%. The highest ULP values were observed in the central and western areas, particularly in the urban district of Wuhan. By contrast, the counties in the south, southeastern, and northeastern areas presented the least urbanization values.
(a) (b)   Qingshan District, Tieshan District, and Xisaishan District declined, but the decrease was less than 10%. The highest ULP values were observed in the central and western areas, particularly in the urban district of Wuhan. By contrast, the counties in the south, southeastern, and northeastern areas presented the least urbanization values.
(a) (b)    Table 4   After performing a regression analysis using different combinations of the significant variables in 2010 and 2015, we finally selected fixed asset investment (PFAI i (yuan/m 2 )), agricultural output (PAU i (yuan/cap)), and income (DPIU (yuan/m 2 )) as the potential driving factors of land urbanization. These three factors were calculated as the ratio of fixed asset investment to the total land area, the ratio of agricultural output to the total population, and the per capita income of urban residents.  Table 4   After performing a regression analysis using different combinations of the significant variables in 2010 and 2015, we finally selected fixed asset investment (PFAIi (yuan/m 2 )), agricultural output (PAUi (yuan/cap)), and income (DPIU (yuan/m 2 )) as the potential driving factors of land urbanization. These three factors were calculated as the ratio of fixed asset investment to the total land area, the ratio of agricultural output to the total population, and the per capita income of urban residents. The scatterplots of PFAI, PAU, and DPIU with land urbanization in 2010 and 2015 are illustrated in Figure 5a Table 5 presents the spatial interaction of ecological carrying capacities per capita amongst counties calculated through gravity model in the Wuhan urban agglomeration in 2010 and 2015. The strongest spatial interaction was observed between Caidian Disrtict and Hannan District, which are suburban districts of Wuhan city (2.9191 m 2 /cap 2 in 2010 and 3.6507 m 2 /cap 2 in 2015). From 2010 to 2015, the average gravity value of spatial interaction of ecological carrying capacities per capita amongst counties decreased.  Table 6 summarizes the spatial regression results in 2010 and 2015. A spatial relation was observed in the land urbanization process because the spatial coefficients in 2010 and 2015 were significant. The spatial interactions among counties at different ECCs per capita were high in 2010 and 2015. The spatial regression results provide a high degree of explanation given that R 2 is higher than 0.9. Note: ** and *** indicate 5% and 1% significance levels, respectively.

Discussion
Land urbanization is a spatio-temporal process that involves capital investment, infrastructure construction, social and government support, and ecological adjustment [61][62][63]. During this process, urban land expansion is inevitable, and ECC is proved to be an effective way to prevent eco-environmental problems and balance the socio-economic development and sustainable ecosystems [20,64]. At the regional level, the contribution of these factors to the expanded urbanized land and the effects of the spatial interactions generated from eco-constraints on this process are of vital importance for the comprehension and optimization of sustainable land urbanization [65,66].
The primary contributions of our study lie in two aspects, namely, the exploration of the spatial effects of ECC on land urbanization and the quantification of land urbanization using the densi-graph method. On the one hand, land urbanization has been closely related to urban expansion and is also generally accompanied by a decrease in ecological land [67][68][69]. In China, the proportion of ecological land space decreased by 0.85%, whereas the proportion of urban land space increased by 1.53% from 1984 to 2012 [17]. A high proportion of ecological land predictably yields high ecological footprint and ECC [70]. Therefore, the implicit and crucial relationship between ECC and land urbanization is plausible. However, previous studies primarily focused on the effects of land use change on ECC or the constraints and policy implications of ECC for the optimization of land use locally [20,71]. Spatial interaction, that is, the influences of ECC in one spatial unit on its neighboring spatial units, have seldom been explored in the context of ongoing land urbanization. In the contemporary socio-economic strategy, regional development has been underlined, and the evolution of the eco-environment has demonstrated increased regional characteristics [72,73]. Therefore, we used a gravity model to quantify the spatial interactions generated by ECC and embedded it into the spatial regression model to explore the importance and magnitude of this spatial influence in the Wuhan agglomeration. The average ECC per capita decreased, whereas the spatial interactions among the counties were substantial and exhibited a slight decline in 2010−2015. The results confirm that the spatial influence of ECC is a factor that cannot be neglected in promoting land urbanization. In fact, ECC reflects the ability of accommodating life and human production in one region, and the spatial influence among counties embodies resource allocation, transfer, compensation, and optimization [74]. Counties with high gravity induced by ECC are inclined to achieve a high degree of effective allocation and sustainable use of resources through various and vibrant trade-offs, which help coordinate regional development to achieve sustainable land urbanization [75]. Meanwhile, the quantification of urbanized land is essential in our study. We applied the densi-graph method based on POIs to obtain the urbanized land area. The reason behind this choice is that urban infrastructure on urbanized land should be sufficient and that the spatial distribution of these POIs should conform to certain density patterns [76]. We also compared the urbanized land area extracted through the densi-graph method and the urban area extracted from remote sensing images. The comparison showed that the urbanized land area is slightly higher than the urban area because we included suburban areas or urban fringes with considerable urban infrastructure and great development potential. In this sense, the measurement of land urbanization is accurate given the focus on the "real" construction and terrestrial attributes of urbanized land.
PFAI, PAU, and DPIU were found to be correlated with land urbanization, with PFAI being the most powerful. In general, the government tends to make investments in regions with favorable economic environment to enhance urban infrastructure and consequently accelerate the rapid expansion of urban built-up areas [77]. In addition, improved disposable personal income and consumption are also supposed to urge the government to expand urban land areas for commercial use, public facilities, and residence [78]. An increasing number of rural residents are also supposed to migrate to urban areas for job opportunities instead of engaging in agricultural activities [79]. Thus, a negative relationship between agricultural output and urbanization rate was observed. The use of strategic opportunities is advisable to achieve a new type of urbanization in which agricultural development, social opportunities, and economic investment are the key factors.

Conclusions
This study explored the spatial influence of ECC on land urbanization and the underlying driving forces using the Wuhan urban agglomeration as a case example. We devised spatial models with ECC-based spatial weight matrices for examining the potential spatial interactions or constraints and the influencing factors of land urbanization. Specifically, the densi-graph method based on the kernel density of POIs was used to derive the "real urban land," and the ecological footprint analysis method was applied to calculate ECC. The results revealed the following; (1) Land urbanization in most counties increased, whereas the average ECC per capita decreased in the Wuhan urban agglomeration; (2) Fixed asset investment and urban income exerted positive influences on land urbanization, whereas agricultural output was negatively correlated with land urbanization; (3) Spatial interaction was formulated through ECC during the land urbanization process; however, this effect was attenuated slightly in 2010−2015.
In fact, spatially explicit models are preferable in exploring spatial interactions at the regional level, especially in the context of spatio-temporal processes such as land urbanization. However, spatial influences are not confined to ECC; other social, economic, political, and technological factors, such as social network and transportation, can also been considered as potential candidates to gauge spatial effects on land urbanization. In addition, we applied the POI-based densi-graph method to identify urbanized land at the county level. This approach improved accuracy to a certain extent. However, we applied fixed thresholds, which are incapable of capturing the local differences in different counties, when kernel density was analyzed. Thus, the adjusted threshold is also a potential tool to improve our method. In the future, other integrated, dynamic, and spatially heterogeneous methods for exploring the driving mechanisms of land urbanization should be proposed to achieve sustainable urban development.