A Fractal Approach to Urban Boundary Delineation Based on Raster Land Use Maps: A Case of Shanghai, China

Abstract

Given the diverse socioecological consequences of rapid urban sprawl worldwide, the delineation and monitoring of urban boundaries have been widely used by local governments as a planning instrument for promoting sustainable development. This study demonstrates a fractal method to delineate urban boundaries based on raster land use maps. The basic logic is that the number of built-up land clusters and their size at each dilation step follows a power-law function. It is assumed that two spatial subsets with distinct fractal characteristics would be obtained when the deviation between the dilation curve and a straight line reaches the top point. The top point is regarded to be the optimum threshold for classifying the built-up land patches, because the fractality of built-up land would no longer exist beyond the threshold. After that, all the built-up land patches are buffered with the optimum threshold and the rank-size distribution of new clusters can be re-plotted. Instead of artificial judgement, hierarchical agglomerative clustering is utilized to automatically classify the urban and rural clusters. The approach was applied to the case of Shanghai, the most rapidly urbanizing megacity in China, and the dynamic changes of the urban boundaries from 1994 to 2016 were analyzed. On this basis, urban–rural differences were further explored through several fractal or nonfractal indices. The results show that the proposed fractal approach can accurately distinguish the urban boundary without subjective choice of thresholds. Extraordinarily different fractal dimensions, aggregation and density and similar average compactness were further identified between built-up land in urban and rural areas. The dynamic changes in the urban boundary indicated rapid urban sprawl within Shanghai during the study period. In view of the popularization and global availability of raster land use maps, this paper adds fuels to the cutting-edge topic of distinguishing the morphological criteria to universally describe urban boundaries.

1. Introduction

Today, more than 55% of the global population resides in urban areas, and this figure is projected to reach 68% by 2050 [1]. Due to the need of land for urban development, large areas of rural land have been converted into urban land for residences, transportation, industry or education [2]. In China, urban built-up land increased by 124.8% from 2000 to 2015, more than doubling the urban population growth rate [3]. Such rapid urban sprawl has raised diverse social and environmental problems of sustainability, including underutilized infrastructure, poverty and inequity, food insecurity, farmland loss, biodiversity decline and climate degeneration [4,5,6]. In this regard, how to restrict uncontrolled urban sprawl has become a hot debate in both scientific and political communities concerned with sustainability.

Aimed at curbing uncontrolled urban sprawl, Smart Growth was proposed in America in the 1990s. It emphasizes sustainable development of the environment, society and economy, which represents a compact, concentrated and efficient schema. In this context, the urban growth boundary (UGB) was applied by local governments as a planning instrument for promoting higher urban densities and protecting rural lands such as basic farmlands [7]. UGB has been implemented in many countries all over the world, such as the United States of America [8], Canada [9] and Great Britain [10]. In recent years, the former Ministry of Land and Resources and the Ministry of Housing and Urban-Rural Development in China have tried to establish the urban development boundary (UDB, similar to UGB) in 14 pilot cities, including Beijing and Shanghai. The Chinese Academy of Urban Planning and Design defines the UDB as a boundary in space that divides the urban administrative district into a construction expansion permitted zone and a construction expansion prohibited zone. Both UGB and UDB are considered to be effective in controlling urban sprawl; however, heated issues of their establishment and prediction remain [11].

To evaluate and improve existing UGB and UDB, many scholars concentrate on urban boundary delineation [12,13,14]. There is no consensus at present about the best approach to delineating the urban boundary in terms of criteria, thresholds or methods [15]. Several methodologies in different aspects have been proposed for the delimitation of an urban boundary, such as morphological delimitation, demographic delimitation, delimitation based on the economic and social structure, and functional delimitation [16]. Most use buildings extracted from digital maps and buffer with a series of number and distance thresholds to realize the discontinuity spatially. These methods may work well in their respective contexts, but they do not remove user experience such as a predefined threshold or continuity constraint [17], which seem to be hardly comparable and subjective. Although some reasonable methods have been proposed, such as a city clustering algorithm using census data [18] or street nodes extracted from OpenStreetMap [19], they exactly prescribe the clustering resolution, either the length of a square or the radius of a circle, depending on experience.

Fractal theory has received increasing popularity in the urban literature over recent decades. Fractal geometry, as an effective mathematical tool, has the potential to analyze chaotic and fragmentized geographic phenomena. Many studies have explored and discovered the characteristic of self-similarity in land use systems, which promotes the fractal theory applied to the explanation of laws in urban spatial distribution, the morphology of urban boundaries [20,21,22] and the deforestation process [23,24]. Tannier et al. proposed a common fractal approach to identifying urban boundaries based on vector-built maps [25] and then compared six archetypal theoretical cities with Belgium’s 18 largest towns [26]. The results indicate that such a fractal method is steady, comparable and fruitful. However, whether or not such a method works well for raster land use maps remains unknown. The varying data structure of raster data and vector data necessitates the computing procedure of fractal analysis for raster land use maps, so as to boost the comparability of such methods. More specifically, although the fractal method can be adapted to the raster data structure in theory, it may have difficulties when applied to raster data. One problem is that raster resolution cannot satisfy the needs of exponential growth of the buffer radius similarly to vector data. Another limit of such a method is the lack of approach distinguishing urban and rural clusters automatically from rank-size distribution, which reduces comparability of urban boundaries between different spatial or temporal dimensions. Additionally, we are aware that globally accessible land use maps recently released by international academic institutions are all recorded in raster structure. Hence, a fractal approach to urban boundary delineation based on raster land use maps should open the door for distinguishing the morphological criteria to universally describe urban boundaries.

Our study aimed to close such a research gap through proposing a fractal approach to deriving urban boundaries from raster land use maps. This approach was further applied to delineate urban boundaries for different years and explore their dynamic characteristics in Shanghai, the largest megacity in China. The contributions should be twofold: (1) methodologically, the reliability of fractal analysis for extracting urban boundary based on raster land use maps is evidenced; and (2) practically, spatiotemporal dynamics of urban boundaries in Shanghai are captured, with essential implications for land use planning.

2. Materials and Methods

2.1. Study Area

Located at 31°14′ N and 121°29′ E, Shanghai is one of the municipalities directly under the central government in China, and it is the main city in the Yangtze River delta (Figure 1).

Shanghai is close to the estuary of the Yangtze River and is a plain densely covered with rivers. In addition, Shanghai has a highly developed transportation system, which jointly leads to fragmented landscapes. Shanghai has approximately 6833 square kilometers of land area and more than 24 million permanent residents within 16 municipal districts. With regard to the overall spatial distribution of built-up land in Shanghai, the Huangpu River runs through the land and separates it into four parts in four directions together with another two narrower artificial rivers. This fragmented landscape determines the spatial pattern of urban morphology to some degree, which implies that the urban area is more likely composed of several agglomerations instead of a huge one. Since the 1990s, Shanghai has experienced rapid urbanization [27]. The problems caused by urbanization such as population explosion, land shortage, and worsening environmental pollution have kept Shanghai from modernization construction. In China, urbanization took a formal modality of industrial zones and satellite towns to maximize land income or to open up new spaces for alleviating the main city zone from the pressure of land and population growth by a city plan [3]. Especially in Shanghai, the city master plan has a significant impact on city growth and contributes to the operation of the local growth engine [28,29], which makes its spatial distribution of built-up land more polycentric. In this context, identifying the urban boundary and exploring the dynamic changes for years in Shanghai has a good demonstration effect and referential significance to other modern cities undergoing rapid urbanization. It should be explained that we only focused on the metropolitan area of Shanghai city and the surrounding islands were not taken into consideration, because these islands are generally rural spaces and follow a quite different urbanization process from the metropolitan area.

2.2. Data and Preprocessing

Raster built-up land maps (1994, 2000, 2003, 2006 and 2016) with 10 m × 10 m resolution were provided by the local government. These maps were converted from the official digital land use maps at 1:10,000 scale in vector format following the data confidentiality standard. What resolution to choose depends on the actual situation of built-up land distributions in different cities. Normally, it should not be too large to discover the distance threshold, nor be too small owing to algorithm time-consumption and real data resolution [27]. The 10 m resolution was chosen because this was accurate enough to reflect hierarchical differences in roads within this study area. According to the national standard, the built-up land includes residential land, public service land, commercial land, industrial land, warehouse land, transportation land, and squares. The premise of the fractal approach is hierarchical empty lanes between built-up land patches; therefore, we removed transportation land and squares to split the built-up land. The raster maps are shown in Figure 2. Detailed information is presented in Table 1.

2.3. A Fractal Approach to Urban Boundary Delineation

2.3.1. Theoretical Foundation

The term ‘fractal’ was first proposed by Mandelbrot, with the initial meaning of irregular and broken. It was defined as a rough geometric shape that can be subdivided into smaller parts, which are the copy of the whole [30]. Fractals are born in nature and often used to describe intricate natural phenomena, such as the length of a coastline, the shape of a cloud and the outline of a leaf. Fractal geometry has been applied in various fields, and has been widely adopted in geography for approximately 30 years. A number of publications have discussed the multiscale spatial organization resulting from urban growth [31,32,33]. Others explore the statistical scaling relationship between the size and number of urban clusters [25,34,35]. Fractals such as Fournier dusts have strict hierarchy between distances and numbers of empty lanes: their numbers will increase as their widths become narrower. Such characteristics could be recognized by the dilation method proposed by Minkowski [36], as shown in Figure 3. In dilation, every element is surrounded by a buffer of increasing width. It is obvious that the number of clusters will ultimately decrease to one, which is to say that there exists only one enormous cluster (Figure 3d). If a spatial organization resembles Fournier dusts, the number of clusters N and the size of buffer width ε will follow a power-law function. If we take the logarithm of both sides, a linear relationship with slope equaling to fractal dimension D would appear (Equation (1)).

(1)$\log \left(\mathrm{N}\right)=\log \left(\mathrm{m}\right)+\mathrm{D}·\log \left(\mathsf{\epsilon }\right)$

where m is a constant of proportionality. If we draw a log–log plot about N and ε, and use a curve to fit it, we can obtain a dilation curve reflecting the spatial scaling behavior of patterns.

Built-up land in an urban area (except transportation land and squares) (Figure 4b) is morphologically similar to Fournier dusts (Figure 4a). Both are separated by hierarchical streets (lanes) following the rule that the narrower, the more. In contrast, built-up land in a rural area (Figure 4d) is more random and has an irregular spatial organization (Figure 4c). Tannier built three theoretical urban patterns: (1) regular fractal city, (2) random fractal city, and (3) regular fractal city in a nonfractal environment (an appropriate substitute of a real city involving urban and rural areas) [25]. The dilation curves of these three patterns express different scaling characteristics. Only the third pattern shows a sharp change in its dilation curve. In the third theoretical urban pattern, the urban area is a fractal pattern which is well-bedded, hierarchical, high-density and agglomerative, whereas the rural area is amorphous, random, low-density and scattered. Different spatial organization leads to different scaling behavior, so that we may find a significant change in one city’s dilation curve, separating urban and rural areas.

According to the aforementioned findings, Tannier et al. developed a fractal-based approach to identify the urban boundary [25]. Although this method is adapted to both vector and raster data in theory, it may have difficulties when applied to raster data. One problem is that raster resolution cannot satisfy the need of exponential growth of buffer radius similarly to vector data. Another limit of this method is the lack of approach distinguishing urban and rural clusters automatically from rank-size distribution, which reduces the comparability of urban boundaries between different spatial or temporal dimensions. Therefore, this paper proposed an improved method to delineate urban boundaries.

2.3.2. Urban Boundary Delineation

The basic logic of the fractal approach to delineating urban boundaries is that the number of built-up land clusters and their size at each dilation step should follow a linearly shaped power-law function [25]. It is assumed that two spatial subsets with distinct fractal characteristics would be obtained when the deviation between the dilation curve and a straight line reaches the top point. The top point is regarded to be the optimum threshold for classifying the built-up land patches, because the fractality of built-up land would no longer exist beyond the threshold. After that, all the built-up land patches are buffered with the optimum threshold and the rank-size distribution of new clusters can be re-plotted. A challenge emerges in the last step to distinguish urban and rural clusters. Instead of artificial judgement, the focus of our approach is to automatically classify the urban and rural clusters. Figure 5 shows the flowchart of the proposed approach and details for each step are described as follows:

First, every built-up patch will be covered by a steadily widening buffer in every iteration, and the buffer radius ε and the number of built-up clusters N will be stored.

Second, to ensure that ε can exponentially grow, we adopt a method of interpolation: cubic hermite interpolation (CHI) [37,38]. Cubic hermite interpolation can ensure that both function value and 1st–nth derivatives remain the same between the interpolation polynomial and original function. Especially, we can assign derivatives of the interpolation polynomial in the node. In view of the distribution of log(ε) and log(N) monotonically decreasing, the closer two nodes are, and the more similar their 1st derivatives are. Therefore, we calculated every 1st derivative in the node with the following equation:

(2)$\mathrm{y}{\left({\mathrm{x}}_{\mathrm{i}}\right)}^{\prime }={\mathrm{k}}_1·\frac{{\mathrm{d}}_2}{{\mathrm{d}}_1+{\mathrm{d}}_2}+{\mathrm{k}}_2·\frac{{\mathrm{d}}_1}{{\mathrm{d}}_1+{\mathrm{d}}_2} (\mathrm{i}=1, 2,\dots , n-2)$

where n is the number of nodes, k₁ and k₂ are the slope from node i − 1 to node i and that from node i to node i + 1, respectively, whereas d₁ and d₂ are the distance from node i − 1 to node i and that from node i to node i + 1, respectively. The 1st derivatives of the first and last node are equal to slopes of their nearest node to themselves, respectively. In this paper, we set 0.1 as the sample interval, which meant that the position (log(ε₀) + 0.1 ∗ n) could be interpolated. Meanwhile, the first and last points will be reserved.

Third, log(N) and log(ε) were fitted with a polynomial from the 1st to 13th degree to obtain 13 estimated curves. Then, the best dilation curve was selected according to the Bayesian information criterion (BIC). BIC can evaluate the trade-off between complexity and goodness of fit for statistical models [39], as shown in Equation (3), where n is the number of data, p is the number of model parameters, and σ² is the mean-square error (MSE).

(3)$\mathrm{BIC}=n·\ln \left({\mathsf{\sigma }}^2\right) +p·\ln \left(n\right)$

This equation shows that only when the polynomial fits well and its complexity remains low does the BIC value reach the minimum. After we selected the best estimated dilation curve, we needed another indicator to catch the discontinuity in scaling behavior: curvature k (Equation (4)).

(4)$\mathrm{k}=\left(\frac{{\mathrm{y}}^{″}}{{(1+{\mathrm{y}}^{\prime })}^{3/2}}\right)$

where y′ is used to measure the velocity of the decreasing number of clusters, y″ is its accelerated velocity, and k is the ratio of y″ and y′. BIC: Bayesian information criterion; Existing method: the fractal approach [25].

In this study, the point of main curvature of a dilation curve with minimal BIC corresponds with the distance threshold. To avoid the point of main curvature arising from estimation artifacts, two extremities of the estimated curve should be removed before calculating the curvature [25]. If there is not a point of main curvature within a valid range, other dilation curves with secondary BIC values would not be selected until a meaningful curvature appeared. The data range can be determined by Equation (5), the same as the Morpholim software [40]:

(5)$\left\{\begin{array}{c}{\mathrm{R}}_{\min }={\mathrm{X}}_{\min }+({\mathrm{X}}_{\max }-{\mathrm{X}}_{\min })\times 0.2\\ {\mathrm{R}}_{\max }={\mathrm{X}}_{\max }-({\mathrm{X}}_{\max }-{\mathrm{X}}_{\min })\times 0.1\end{array}\right.$

where X means log(ε) and R means the valid range of X.

Next, built-up land was covered with a buffer radius equaling the distance threshold. Then, the rank-size distribution of new clusters was plotted. Instead of artificial judgement, we adopted a bottom-up clustering approach named hierarchical agglomerative clustering (HAC) [41] to distinguish urban and rural clusters. In consideration of the characteristics of the rank-size distribution of urban and rural clusters, we chose the single-linkage as the linkage criteria of HAC. This means the minimum distance between any two subsets which determines whether these two subsets are able to merge together or not. Urban clusters seem to be top-ranking, few and large, whereas rural clusters are always lower-ranking, abundant and much smaller; therefore, Euclidean distance is sufficient as a metric. All clusters will be divided into two clusters: the one including the first rank cluster must be the urban set, and the other one is the rural set. The clustering result will be exhibited in hierarchical clustering dendrograms. Finally, we extracted the urban boundary from urban clusters. All of the above steps were performed using Visual C++ programming.

2.3.3. Urban Boundary Characteristics

After we extracted the boundary of the clusters, the city was divided into two spatial forms: urban areas within the boundary and rural areas outside of the boundary. To explore the shape of the boundary and the difference between urban and rural area as the assessment of delineated boundary, we then calculated multidimensional indices [3]. The first one characterized fractal patterns or sets by quantifying their complexity as a ratio of the change in detail to the change in scale [42]. The second one measured the average morphological compact degree of spatial organization. The third one was an index reflecting aggregation effect, and the final one showed the proportion of one landscape in a scope.

• Fractal dimension

The fractal dimension is a ratio reflecting the space-filling capacity of a pattern in changing scales, which does not have to be an integer. If we used rulers of different lengths to measure the length of a coastline, we would not obtain a unique value, because coastlines show more details at higher resolution. This is why we need an index named fractal dimension instead of a perimeter to describe a fractal-like coastline. It can be calculated by different methods, such as the structured walk method [43], equipaced polygon method [44], hybrid walk method [45] and cell count method [46]. However, there is a small diversity in the calculated dimensions [47]. In this paper, we used the box-counting method [48,49] realized in a fractal analysis software package named Fractalyse [50] to calculate the fractal dimension of the boundary (D_B) and surface (D_S). This method uses the exponential growth box size to cover the whole picture and then counts the boxes containing objects. Then, the result is displayed on a log-log plot. If a pattern is fractal, points can be fitted by a straight line whose slope is its fractal dimension.

To be specific, the raster images in the five years were all extracted by the minimum enclosing rectangle of the urban boundary in 2016. These extracted images were converted into BMP format and added into the Fractalyse software, which was downloaded from the official address (https://sourcesup.renater.fr/www/fractalyse/, accessed on 3 March 2020). The upper left corner was used as the start point for box-counting. We chose six theoretical patterns (Figure 6) with known fractal dimensions to test the approximate accuracy of box-counting in Fractalyse: boundary and surface of circle, square, Sierpinski carpet and Fournier dusts. The box size was set following an exponential function whose exponent equaled 2 and the grid algorithm. The result presented a high correlation coefficient of size and number (near 1.0) and tiny deviation (smaller than 0.06), which demonstrates the high precision of the box-counting method in Fractalyse in calculating the fractal dimension (Table 2). Theoretically, the fractal dimension of a point is 0; a line is 1; and a surface is 2, which equals their own topological dimensions. If a square is narrow enough, its fractal dimension may be close to 1, because the space-filling capacity of this pattern finally degenerates to that of a line in large scales. If a line is complex and winding enough or widens, its space-filling capacity will be stronger in small scales than smoother lines, such that its fractal dimension may be close to 2. When we calculate the fractal dimension of sets of lines or squares, it reflects the space-filling capacity of all the elements in changing scales, not a part of them. In theory, the space-filling capacity of built-up land in urban areas is stronger than that in rural areas. With respect to the urban boundary, the minimum size was 2 pixels and maximum size was 8192 pixels. The minimum size was 1 pixel and the maximum size was 8192 pixels for both urban clusters and rural clusters.

• 2.. Area-weighted compactness index

The compactness index (CI) is used to evaluate whether the shape of one patch is compact or not, which in some places is named after the shape complexity. It is calculated using the following formula:

(6)$\mathrm{CI}=\frac{2·\sqrt{\mathsf{\pi }·\mathrm{A}}}{\mathrm{P}}$

where A is the surface area and P is the perimeter.

A circle has the most compact geometry, whose CI reaches a maximum equaling 1. If a patch is very compact, CI is close to 1. In contrast, if a patch has a long perimeter and small surface area, CI is close to 0. We expected to compare overall patches in urban areas to rural areas in terms of the compactness index; therefore, we designed the area-weighted compactness index (AWCI):

(7)$\mathrm{AWCI}={\displaystyle \sum }_{\mathrm{i}=1}^n\left(\frac{\mathrm{Ai}}{\mathrm{A}}\times {\mathrm{CI}}_{\mathrm{i}}\right)$

where A means the total area.

We calculated the AWCI of urban and rural built-up land to explore the difference in shapes between two such patterns.

• 3.. Aggregation index

The aggregation index (AI) is a global measurement regarding the aggregation of landscapes, which is mathematically described as the ratio of actual shared edges versus maximal possible shared edges [51]. It was calculated as follows:

(8)$\mathrm{AI}=\frac{{\mathrm{g}}_{\mathrm{ii}}}{\max \to {\mathrm{g}}_{\mathrm{ii}}}\times 100\%$

where g_ii is the number of similar adjacencies tallied using the single-count method, whereas the max g_ii is the maximum possible number of similar adjacencies corresponding to the whole landscape.

The single-count method means that if two neighboring pixels belong to the same landscape, it tallies only once. Hence, the maximal aggregation is achieved when the patch type consists of a single and compact patch, which is not necessarily a square patch. We calculated this index with FRAGSTATS v4 [52], which is a computer software program designed to compute a wide variety of landscape metrics for categorical map patterns.

• 4.. Density

Density in this paper is the ratio of built-up land surface area and the urban or rural area, divided by the urban morphological boundary we have defined. Its formula is:

(9)${\mathrm{D}}_{\mathrm{u}/\mathrm{r}}=\frac{{\mathrm{S}}_{\mathrm{u}/\mathrm{r}}}{{\mathrm{A}}_{\mathrm{u}/\mathrm{r}}}$

where D_u/r means the built-up land density in the urban or rural area, S_u/r means the built-up land area in urban or rural area, and A_u/r means the whole area in the urban or rural region.

By comparing D_u to D_r, it is easy to observe the difference in built-up land distribution in spatial scales between urban and rural areas, contributing to further explorations of the functions of such morphological boundaries.

2.4. Spatiotemporal Analysis

On the basis of urban area extracted through years, we further explore the spatiotemporal characteristics of urban growth in Shanghai. We quantified the sprawl extent over intervals of five years according to the scale of district/township. There were 197 district/townships within the study area. First, we counted the number of pixels transformed into urban areas from rural areas during these periods. Secondly, we calculated the proportion of transformed areas within every district/township. Finally, all proportions were reclassified into four classes by the Natural Break Jenks method: non-sprawl, light sprawl, medium sprawl, and high sprawl.

3. Results

3.1. Urban Boundary Delineation and Urban–Rural Differences

The values of box sizes for corresponding box numbers are shown in Table 3 for urban boundaries, Table 4 for the urban clusters, and Table 5 for the rural clusters. Table 6 presents detailed results about the best estimated dilation curve and the distance thresholds from 1994 to 2016. In 1994, the threshold reached 138.995 m, has clearly shortened since 2000, and then maintained around 20–40 m. The results of fitting data by CHI, the best estimated dilation curves, and thresholds for five years are shown in Figure 7a–e.

Figure 8a–e shows the rank-size distributions of built-up land covered by distance thresholds in each year, whereas the hierarchical clustering dendrograms (Figure 9a–e) reflect different quantities of urban clusters in different years. In 1994, only the largest cluster belongs to urban area, but the 2, 2, 2, and 3 largest clusters belong to urban area in 2000, 2003, 2006 and 2016, respectively.

The identified urban boundaries for five years are shown in Figure 10. The results of the qualitative analysis can be summarized as follows: (1) In 1994, built-up land clusters near the east and the west of Huangpu River merged into one large morphological agglomeration because of the long threshold radius of 140 m. Such a scope just conforms to the initial central urban area in Shanghai, which includes most of the range of Huangpu, Hongkou, Yangpu, Jingan, Putuo, Changning, and Xuhui, the east region of Baoshan and a small region in the northwest of Pudong. (2) In 2000, built-up land was separated into east and west regions by the Huangpu River, and the west area was larger than the east area, which corresponded to their development levels. In the west, the urban area sprawled over most of Minhang west of the Huangpu River and the middle area of Baoshan. In the east, the urban area rapidly covered the west part of Pudong, which illustrates the rapid development in Pudong in the late 1990s. (3) In 2003, the urban area in the east and west regions showed a striking contrast: the urban area in the west region grew quickly and began to invade Jiading and Songjiang, much more quickly than the east regions. Meanwhile, the west urban boundary seems highly tortuous, with a morphology of tentacular structures similar to Sierpinski carpets, which occurs where interstitial spaces are filled in along transportation axes. (4) In 2006, the urban area in Pudong sprawled sharply toward the south and the west part continued to sprawl toward the west. (5) In 2016, the north region of Shanghai had nearly complete “urbanization” because of the geographic restriction of the administrative division. The region of Minhang in the east of Huangpu River finally became an urban area and some space in Fengxian was transformed into an urban area because of the developing transportation industry. Notably, the third largest cluster R3 combined with regions of Jinshan and Fengxian was divided into an urban area. Although it was smaller than the largest two clusters, it clearly exhibited a big leap compared with the fourth cluster, which could not be ignored. In general, Shanghai seems to have transformed its spatial structure from a monocentric city to a polycentric one.

Table 7 shows the results of multidimensional indices of urban boundaries and urban–rural differences. From 1994 to 2016, D_B first increased (1994–2000), and then maintained around 1.4 (Figure 11a). The explanation of such a phenomenon may be as follows. In 1994, the urban area concentrated on the central region whose form was relatively compact so that D_B was not high. From 2000 to 2016, the urban boundaries showed similar space-filling capacity. Furthermore, as the urban area expanded, more and more rural areas are gradually surrounded by urban area, which mainly occurs around the urban fringe. In general, D_B was maintained between 1.3 and 1.41, which reflects the complexity of its morphology.

The average of D_S of the built-up land in urban areas was 1.743, whereas that in rural areas was 1.589, which means that the space-filling capacity in urban areas is much better than that in rural areas. In urban areas, D_S increased gently, close to a line (Figure 11b, y1) with a slope of 0.0051 (R² = 0.8368), whereas in contrast, D_S in rural areas had a nonsignificant downward trend with a slope of −0.0035 (R² = 0.1492) (Figure 11b, y2). This indicates that built-up land in urban and rural areas has different characteristics in terms of space-filling capacity and opposite development tendencies.

AWCI values for urban and rural built-up land seem indiscriminate but irrelevant at the level of the patch. They basically remain at 0.62, reflecting that the overall shape of urban and rural built-up land patches is neither too compact nor too fragmented. AI values for urban and rural built-up land are obviously high enough that all of them are larger than 90%. The AI value changes with map resolution and measurement resolution, but fixing them makes the results comparable. As shown in Table 7, all the urban AI values are larger than rural ones, which means that the degree of aggregation in urban areas is higher than that in rural areas. Hence, on the level of patterns, there exists a significant difference between urban and rural built-up land distinguished by the urban boundary. There is also a clear difference between urban and rural built-up land density, indicating that the area proportion of urban built-up land in urban scope is much higher than that of rural land in rural scope, as expected.

3.2. Spatiotemporal Changes of Urban Boundaries

Figure 12 reveals the sprawl extent in Shanghai at the level of District/township. Four periods basically follow the rule that urban sprawls are more severe at the urban boundary. Non-sprawling and light sprawling areas are mainly distributed in the city center and outskirts, whereas the majority of medium and high sprawling areas are more likely to appear near the boundary. Table 8 shows statistical results for the number and land area percentage of District/township in four kinds of sprawl and new urban area per year during four periods. In view of the number of Districts/townships in different sprawls, non-sprawl becomes the main force. New urban area per year reflects fast urban sprawl from 1994 to 2000, whereas since 2000, the speed of urban sprawl slows down gradually. However, considering the percentage of land area in different sprawls, it is obvious that the scope of urban sprawl is expanding, and the percentage of more severe urban sprawl is increasing, which implies that the sustainable relationship between development and nature has worsened even further.

4. Discussion and Conclusions

Most recent studies delineating urban area are based on data from remote sensing images, such as Landsat data or night-time stable light data (NTL), have developed many common methods such as a Bayesian approach, cluster-based method, quantile-based model and land-use information entropy [53,54]. However, these methods have several shortcomings: (1) the rough resolution of remote sensing images essentially results in an ambiguous urban boundary; (2) it is difficult for local government or small-scale countries to benefit from these models; and (3) in a similar context, it is difficult to determine the strengths of different methods.

In this paper, we improve the fractal-based method of [25,26] based on raster built-up land maps, which resolves two problems: (1) how to deal with the resolution problem when applying such a method to built-up raster maps; and (2) how to automatically distinguish urban and rural clusters from rank-size distribution. It is noteworthy that there are three major challenges to this method. First, the resolution of a built-up land map has a lack of standards. In theory, the higher the resolution, the more accurate the boundary would be, but computational efficiency must decline paradoxically. On the other hand, too low a resolution may ignore the necessary scaling behavior difference and cut down the amount of statistical data. Second, built-up land must remove roads, which raises difficulties in data acquisition. Finally, it is difficult to evaluate the accuracy of urban boundaries directly, owing to different official criteria for defining urban areas and setting urban boundaries [55]. Most official land use classification maps used to assess urban boundary extraction methods have been handled by researchers according to different standards; thus, past urban boundaries are not comparable.

In this study, multidimensional indices were applied to reveal the characteristics of urban morphological boundaries and the possible quantitative difference between urban and rural built-up land. Our results present extraordinarily different fractal dimensions (e.g., the number of boxes), aggregation and density and similar average compactness between built-up land in urban and rural areas. This indicates that the fractal method used in this research can effectively distinguish the boundary between two spatial distribution organizations in morphology, which may contribute to identifying the urban boundary in other modern cities.

Urban boundary delineation as the basis of urban sprawl measurement and analysis has been emphasized and used for policy assessment by many scholars [13,28]. Different from most urban sprawl analysis methods, this paper proposes an integrated urban sprawl analysis approach based on the urban boundary mentioned above. A new urban area visual map for different periods presents the extent of local sprawling all over the city. In addition to the urban boundary, urban–rural differences in multidimensional indices are developed. The urban/rural dichotomy seems increasingly insufficient for distinguishing urban and rural scope due to blurred urban–rural differences according to a world economic and social survey [55]. However, given the demand for sustainable development in the rapid urbanization context, the urban boundary still has significant functions. First, the urban boundary as defined in this study presents the spatial discontinuity of built-up land in scaling behavior, which is valuable in exploring the urban morphology and development history. In addition, if the proposed method could not identify the distance threshold, it means the scaling behavior of the whole built-up land seems homogeneous, implying that there hardly exists a morphological difference between urban and rural areas. Second, many policy assistance tools and assessment indices of sustainable development may be integrated, improved or reconstructed based on urban boundaries, such as urban–rural difference monitors, urban sprawl measurement and driving force analysis. In addition, urban boundaries can assist in decision making and facilitate urban–rural synergistic planning. For example, a local boundary where excessive changes appear needs enhanced land use monitoring by a local department. It also supplies a reference for making or evaluating UDB, which has been proposed in many pilot cities, such as Shanghai (Shanghai City Master Plan, 2017–2035). We thus argue that monitoring urban boundaries should be taken into account and put into practice, because it can reflect urban morphology and the extent of sprawl. Establishing a normative urban sprawl evaluation system to track rapid urbanization can be valuable for land use planning. Local surveying and mapping departments should unify the type and precision of land use, keeping the annual land use map comparable.

Different from most urban sprawl analysis methods, this paper proposes an integrated urban sprawl analysis approach based on the urban boundary mentioned above. A new urban area visual map for different periods presents the extent of local sprawling all over the city. Finally, establishing a normative urban sprawl evaluation system to track rapid urbanization is another necessary measure. Local departments of surveying and mapping should unify the type and precision of land use, keeping the annual land use map comparable. Monitoring urban boundaries should be taken into account and put into practice, because it can reflect urban morphology and the extent of sprawl. In addition to the urban boundary, urban–rural differences in multidimensional indices, a new urban area map and the driving force analysis proposed in this paper can also facilitate evaluating urban sprawl and decision-making.

However, there are several limitations to this study. First, this paper did not explore the meaning of changing distance thresholds by time. It is hard to ensure that a digital land use map is strictly accordant in the aspects of properties and precision. In addition, an interval of three years is a relatively high temporal resolution for the comparison between distance thresholds. In theory, urbanization brings a more compact built-up land pattern that results in a narrower distance threshold to a monocentric city. Secondly, such urban boundaries may cover the natural landscape such as rivers, lakes and mountains. Some geo-processing operations such as overlap or cutting will optimize the urban boundary. Thirdly, dilation in this study was not real or could not be developed given the sea and Yangtze River delta urban limit. Future research should focus on the following issues: (1) How and how much does the resolution influence the urban boundary delineation?; (2) How can one assess the accuracy and validity of the urban boundary with a more scientific approach?; and, (3) How can one predict complex urban boundaries in the future and then curb urban sprawl?

Author Contributions

Conceptualization, C.Z., Y.L. and M.W.; methodology, C.Z. and Y.L.; software, Y.L.; validation, C.Z. and M.W.; formal analysis, C.Z.; investigation, C.Z.; resources, C.Z.; data curation, M.W.; writing—original draft preparation, C.Z. and Y.L.; writing—review and editing, M.W.; visualization, Y.L.; supervision, M.W.; project administration, M.W.; funding acquisition, M.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China, grant number 2016YFB0501403.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

This research received data from the local government and the data are confidential.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location of the study area (a), the administrative division (b) and the areal imagery in 2020 (c).

Figure 2. Spatial distribution of built-up land in Shanghai.

Figure 3. Minkowski dilation process.

Figure 4. Theoretical and real urban and rural pattern.

Figure 5. Flowchart of urban boundary delineation.

Figure 6. Six experimental patterns.

Figure 7. Results of dilation curve (a–e) from 1994 to 2016. The red, blue, green and gray dotted lines represent 40%, 50%, 60%, and 70% of urban coverage, respectively.

Figure 8. Results of rank-size distribution (a–e) from 1994 to 2016.

Figure 9. Results of hierarchical clustering dendrograms (a–e) from 1994 to 2016.

Figure 10. Identified urban boundaries where the R1–3 mean the 1st, 2nd and 3rd rank cluster, respectively.

Figure 10. Identified urban boundaries where the R1–3 mean the 1st, 2nd and 3rd rank cluster, respectively.

Figure 11. (a,b) represent the fractal dimensions of boundary and built-up land respectively, where CC means the correlation coefficient of log(r) and log(N).

Figure 12. Urban sprawl extent within districts of Shanghai, China.

Figure 12. Urban sprawl extent within districts of Shanghai, China.

Table 1

Information of raster maps for five years.

Year	SR(m)	Type	Region (Pixel)	Proportion (%)
1994	10	Built-up land	10,651 × 9072	12.9
2000	10	Built-up land	10,651 × 9072	26.0
2003	10	Built-up land	10,651 × 9072	31.6
2006	10	Built-up land	10,651 × 9072	37.1
2016	10	Built-up land	10,651 × 9072	45.6

SR: spatial resolution; Proportion: percentage of built-up land area in Shanghai inland area (about 542,647.49 ha).

Table 2

Comparison between theoretical fractal dimensions and experimental values in Fractalyse.

Pattern	Type	Dt	Dp	Deviation	CC
Circle	B	1	0.95	0.0523	0.999988
Square	B	1	0.98	0.0151	0.999765
Circle	S	2	1.97	0.0260	0.999989
Square	S	2	1.97	0.0270	0.999956
Sierpinski carpet	S	1.89	1.88	0.0138	0.998266
Fournier dusts	S	1.51	1.52	0.0051	0.979444

Type B: boundary of pattern; Type S: surface of pattern; D_t: theoretical fractal dimension; D_p: practical fractal dimension; Deviation: the absolute value between theoretical and calculated fractal dimension; CC: correlation coefficient of log(r) and log(N).

Table 3

Number of boxes in different years and box sizes for urban boundaries.

	Size	2	4	8	16	32	64	128	256	512	1024	2048	4096	8192
Year
1994		30,328	15,010	7386	3551	1647	686	259	98	37	14	7	4	1
2000		101,274	49,300	23,204	10,030	3996	1488	541	201	75	26	8	4	1
2003		129,854	63,375	30,055	13,166	5277	1967	717	263	100	34	12	4	1
2006		155,394	75,690	35,480	15,318	6128	2334	869	319	123	43	15	5	2
2016		209,098	101,670	47,830	20,503	8173	3087	1141	417	152	54	21	8	3

Table 4

Number of boxes in different years and box sizes for the urban clusters.

	Size	1	2	4	8	16	32	64	128	256	512	1024	2048	4096	8192
Year
1994		3,078,062	799,932	215,369	60,876	17,044	4850	1396	400	120	38	14	5	2	1
2000		7,200,854	1,968,728	543,778	149,380	41,601	11,662	3272	928	273	84	29	11	6	2
2003		10,577,007	2,894,354	803,569	222,804	62,617	17,575	4878	1358	385	113	35	11	6	2
2006		14,231,067	3,879,508	1,072,089	295,691	82,752	23,133	6425	1762	491	143	45	12	6	2
2016		22,501,120	6,281,349	1,791,108	505,638	142,673	39,549	10,698	2867	776	217	64	21	8	3

Table 5

Number of boxes in different years and box sizes for the rural clusters.

	Size	1	2	4	8	16	32	64	128	256	512	1024	2048	4096	8192
Year
1994		3,925,754	1,051,851	298,938	93,472	33,036	13,449	5966	2373	759	220	67	20	8	3
2000		6,915,358	2,089,801	685,612	245,737	93,984	33,705	10,089	2830	786	223	66	20	8	3
2003		6,565,693	1,957,783	631,354	221,221	82,463	29,043	8824	2602	757	218	65	20	8	3
2006		5,912,252	1,739,515	551,198	189,490	69,890	24,942	7851	2414	728	216	65	20	8	3
2016		2,217,723	677,360	222,344	79,529	30,683	11,753	4245	1563	573	193	61	21	8	3

Table 6

Main results of distance threshold calculation.

Year	Range	Degree	BICmin	R2	Curvature	Log(r)	Threshold/m
1994	1.400–2.800	12	−159.439	0.999992	−2.033	2.143	138.995
2000	1.369–2.661	4	−95.537	0.999580	−1.638	1.406	25.468
2003	1.323–2.451	9	−109.335	0.999974	−1.029	1.595	39.355
2006	1.350–2.573	4	−81.763	0.999275	−1.304	1.393	24.717
2016	1.329–2.479	6	−82.586	0.999704	−0.786	1.580	38.019

Abbreviations: Range means the valid range of the logarithm of dilation radius, Degree means the polynomial degree corresponding with the minimal BIC, Log(r) means the logarithm of dilation radius corresponding with point of main curvature, and Threshold means the unique dilation radius distinguishing the urban and rural areas.

Table 7

Several indices to characterize urban and rural area.

		Urban				Rural
Year	DB	DS	AWCI	AI	Du	DS	AWCI	AI	Dr
1994	1.308	1.692	0.589	0.981	0.662	1.536	0.625	0.965	0.079
2000	1.409	1.701	0.630	0.953	0.797	1.655	0.637	0.898	0.153
2003	1.407	1.749	0.641	0.953	0.759	1.643	0.618	0.906	0.163
2006	1.391	1.779	0.601	0.955	0.799	1.624	0.640	0.914	0.162
2016	1.407	1.794	0.629	0.942	0.682	1.489	0.630	0.892	0.104

D_B: the fractal dimension of urban boundary; D_S: the fractal dimension of surface area; D_u: urban density; D_r: rural density.

Table 8

Statistical information about the extent of different sprawls.

	1994–2000	2000–2003	2003–2006	2006–2016
Unit number (LAP) of A	131 (71.33%)	142 (64.17%)	137 (52.12%)	114 (23.85%)
Unit number (LAP) of B	32 (15.41%)	28 (20.38%)	35 (25.98%)	46 (35.91%)
Unit number (LAP) of C	27 (11.00%)	19 (11.57%)	17 (15.02%)	22 (25.60%)
Unit number (LAP) of D	7 (2.26%)	8 (3.88%)	8 (6.88%)	15 (14.64%)
New urban area (km2/year)	76.0	172.8	153.4	132.6

A: non-sprawl; B: light sprawl; C: medium sprawl; D: high sprawl; LAP: land area percentage.