Understanding Spatio-Temporal Patterns of Land Use/Land Cover Change under Urbanization in Wuhan, China, 2000–2019

Abstract

Exploring land use structure and dynamics is critical for urban planning and management. This study attempts to understand the Wuhan development mode since the beginning of the 21st century by profoundly investigating the spatio-temporal patterns of land use/land cover (LULC) change under urbanization in Wuhan, China, from 2000 to 2019, based on continuous time series mapping using Landsat observations with a support vector machine. The results indicated rapid urbanization, with large LULC changes triggered. The built-up area increased by 982.66 km² (228%) at the expense of a reduction of 717.14 km² (12%) for cropland, which threatens food security to some degree. In addition, the natural habitat shrank to some extent, with reductions of 182.52 km², 23.92 km² and 64.95 km² for water, forest and grassland, respectively. Generally, Wuhan experienced a typical urbanization course that first sped up, then slowed down and then accelerated again, with an obvious internal imbalance between the 13 administrative districts. Hanyang, Hongshan and Dongxihu specifically presented more significant land dynamicity, with Hanyang being the active center. Over the past 19 years, Wuhan mainly developed toward the east and south, with the urban gravity center transferred from the northwest to the southeast of Jiang’an district. Lastly, based on the predicted land allocation of Wuhan in 2029 by the patch-generating land use simulation (PLUS) model, the future landscape dynamic pattern was further explored, and the result shows a rise in the northern suburbs, which provides meaningful guidance for urban planners and managers to promote urban sustainability.

1. Introduction

Land use/land cover (LULC) serves as a basic indicator for multiple interdependencies between human society and natural environment, which has a profound sense for human life and earth ecology [1,2,3]. Over the past decades, intensified anthropogenic processes have greatly accelerated LULC changes, which has caused a series of problems for the ecosystem, climate, biodiversity, food security and so forth [4,5,6]. A central focus is on urbanization, which increasingly alters the global landscape, especially for developing countries [7,8]. The changes from natural covers to artificial surfaces bring negative social and ecological consequences and cause great concerns for ecology safety, human settlements and public health in addition to numerous beneficial aspects. Hence, understanding LULC dynamics is critical for promoting urban sustainability and ecological balance.

With the rapid development of Earth observation technology, remote sensing has shown increasing superiority in LULC investigation and provides an effective and spatially explicit alternative to urbanization monitoring due to its advantages of wide coverage, short revisit cycle and a strong view of the current situation, which is emphasized as the “unique view” of urban dynamics [9]. A large number of studies concerning urbanization have been carried out at various scales based on remote sensing techniques. For example, Schneider et al. mapped the global urban extent with Moderate Resolution Imaging Spectroradiometer (MODIS) images [10]. Mertes et al. investigated the urban land dynamics based on MODIS data at the continental scale [11]. Yu et al. monitored urban sprawl at a regional scale by using Landsat images [12]. Additionally, many research studies were committed to investigating the consequent effects of LULC change in the context of urbanization, including environmental deterioration [13], climate change [14], urban heat islands [15], food security [16] and so on.

Time series mappings have been reported to be helpful for exploring spatio-temporal patterns of LULC change and urban evolution law due to its abundant information [17,18]. For instance, Seto et al. explored the landscape dynamics of four Chinese cities based on time series LULC maps [19]. Li et al. investigated urban sprawl and LULC changes in Des Moines and Ames from 1985 to 2015 by time series mapping [20]. Zheng et al. analyzed the land dynamics of 30 megacities in the world in the period 2012–2019 by using time series nightlight remote sensing data [21]. Therefore, it is conducive to take advantage of time series data to characterize the LULC change pattern and urban development mode.

As the largest developing country in the world, China has experienced unprecedented urbanization and significant landscape changes during the past decades of reform and opening-up [17,22]. The urbanization rate of China increased sharply from 17.92% to 59.58% during the period 1978–2018 [23]. Furthermore, many Chinese cities have shown accelerated urbanization pace because of the financial revenue reliance on real estate and land transactions [24]. In recent years, based on time series LULC maps, extensive studies have focused on the urbanization in China and its related LULC changes at various scales, including the whole country [25,26], typical regions such as Zhujiang Delta [27] and Guangdong-Hong Kong-Macau Greater Bay Area [28] and typical mega cities, e.g., Beijing [29], Shanghai [30], Shijiazhuang [31], etc., which greatly improves the understanding of the universality and regionality of urbanization in China and contributed to urban sustainability and ecology protection.

Driven by the rise strategies of Central China, Wuhan has experienced rapid development and urbanization, with frequent landscape changes as an accessory. A set of studies has been carried out to explore the LULC changes in Wuhan and their consequent influences. Kabba et al. simply analyzed the landscape changes in Wuhan and their ecological impacts from 1987 to 2005 based on the classification results of two phases [32]. Li et al. investigated the LULC changes in Wuhan and its underlying physical and socioeconomic driving factors in the period 1990–2010 on the basis of the LULC maps of five target years [33]. Liu et al. explored the land use dynamics in Wuhan and the driving forces in the context of urban-rural construction based on the land investigation results in 1996 and 2009 and found a large difference between the inner city and outer city [34]. Ji et al. analyzed the LULC changes in Wuhan from 2006 to 2012 and evaluated the subsequent ecological security effects [35]. Liang et al. tried to understand the drivers of urban sprawl in Wuhan with the land patch-generating land use simulation (PLUS) model [36]. However, most of these past studies focused on the macro change trends of different land uses, with the evolvement pattern exploration being temporally or spatially coarse-grained. Undoubtedly, this hinders a more conceptual and complete understanding of the development mode of Wuhan, and certain uncertainties and challenges remain for urban land management and planning.

In view of this, in order to better understand the Wuhan development mode, this research study investigates the spatio-temporal patterns of LULC change in Wuhan in detail during the past 19 years in the context of urbanization. The main objectives and contributions of this paper are four-fold. Firstly, the continuous time series LULC maps at a spatial resolution of 30 m in the period 2000–2019 are produced, based on Landsat observations and machine learning with the support vector machine (SVM). Secondly, a detailed temporal change analysis is conducted by a land use theme to find out the change characteristics of major land uses. Thirdly, based on five typical years (2000, 2004, 2009, 2014 and 2019), we comprehensively explore the landscape dynamic pattern of Wuhan from multiple perspectives, including single land use dynamicity, landscape activity and the urbanization mode. Lastly, the future land use structure of Wuhan in 2029 is predicted by the PLUS model [36] in order to identify the landscape dynamic pattern over the next 10 years, which provides a reference for urban planning and policy making.

2. Study Area and Data Resources

2.1. Study Area

Wuhan (29°58′–31°22′N, 113°41′–115°05′E) is located in the east of Hubei, China, and is the capital of the Hubei Province, as shown in Figure 1. It lies in the middle reaches of the Yangtze River, at the intersection of the Yangtze River and its largest tributary, i.e., the Han River. It belongs to the Jianghan Plain, with the plain area accounting for more than 80% of the total area of the city. There are a few mountains distributed in the northeast, northwest and central parts of Wuhan. In addition, Wuhan is known as “the city of hundreds of lakes”, with the Yangtze River and Han River meeting in the center of the city and hundreds of lakes crisscrossing on both sides, forming a complex water network. It belongs to the subtropical monsoon climate, characterized by abundant rainfall, a hot summer and a cold winter, with an annual average temperature of 15.8–17.5 °C and annual precipitation of 1150–1450 mm.

Wuhan has 13 administrative districts in total, including seven metropolitans, i.e., Jiang’an, Jianghan, Qiaokou, Hanyang, Wuchang, Qingshan and Hongshan, and six suburbs, i.e., Caidian, Jiangxia, Huangpi, Xinzhou, Dongxihu and Hannan, with a total area of 8569.15 km². As the economic and geographical center of central China and the core of the “1 + 8” city circle and Yangtze River economic belt, Wuhan has the honor of “the thoroughfare of nine provinces”. Meanwhile, it is also one of the three major intelligence intensive areas in China. Wuhan has developed industry, commerce, technology and tertiary industry and possesses an innovation and entrepreneurship center with global influence, i.e., China Optics Valley. It had a permanent resident population of 12.3265 million and a GDP of 1.56 trillion yuan in 2020.

2.2. Data Resources

In this study, the Landsat observations, which has served for over 40 years and provides wall-to-wall coverage of most areas of Earth with considerable spatial resolution of 30m and a revisit cycle of 16 days [37], were utilized to produce the continuous time series LULC maps for Wuhan, according to the path/row geographies of Landsat Worldwide Reference. Specifically, the level-1 products at path/row of 122/39, 123/38 and 123/39 were downloaded for each year in the period of 2000–2019 from the United States Geological Survey (USGS) Resources Observation and Science Data Center (https://earthexplorer.usgs.gov/) (accessed on 3 August 2020) under the principle of cloud-free or low-cloud coverage, which has been systematically corrected by the providers to exclude the systematic errors of sensors. There were 60 images collected in total, with detailed information listed in

Table 1. Description of Landsat images utilized in this study.

Year	Sensor	Acquisition Date (/Month-Day/) (Path/Row: 122/39, 123/38, 123/39)	Utilized Bands
2000	Landsat 5 TM	9-22/9-13/9-13	B/G/R/NIR/SWIR-I/SWIR-II
2001	Landsat5 TM	7-10/9-3/9-3	B/G/R/NIR/SWIR-I/SWIR-II
2002	Landsat5 TM	5-26/5-1/5-1	B/G/R/NIR/SWIR-I/SWIR-II
2003	Landsat5 TM	9-17/9-24/9-24	B/G/R/NIR/SWIR-I/SWIR-II
2004	Landsat 5 TM	7-18/9-11/9-11	B/G/R/NIR/SWIR-I/SWIR-II
2005	Landsat5 TM	6-19/4-7/4-7	B/G/R/NIR/SWIR-I/SWIR-II
2006	Landsat5 TM	7-8/7-31/7-31	B/G/R/NIR/SWIR-I/SWIR-II
2007	Landsat5 TM	3-4/4-28/4-28	B/G/R/NIR/SWIR-I/SWIR-II
2008	Landsat5 TM	7-13/9-6/9-6	B/G/R/NIR/SWIR-I/SWIR-II
2009	Landsat5 TM	11-5/11-12/11-12	B/G/R/NIR/SWIR-I/SWIR-II
2010	Landsat5 TM	8-20/6-8/6-8	B/G/R/NIR/SWIR-I/SWIR-II
2011	Landsat5 TM	5-10/5-17/5-17	B/G/R/NIR/SWIR-I/SWIR-II
2012	Landsat7 ETM+	8-9/9-17/9-17	B/G/R/NIR/SWIR-I/SWIR-II
2013	Landsat8 OLI	10-15/10-6/10-6	C/B/G/R/NIR/SWIR-I/SWIR-II
2014	Landsat8 OLI	9-22/9-13/9-13	C/B/G/R/NIR/SWIR-I/SWIR-II
2015	Landsat8 OLI	10-18/10-25/10-25	C/B/G/R/NIR/SWIR-I/SWIR-II
2016	Landsat8 OLI	8-1/7-23/7-23	C/B/G/R/NIR/SWIR-I/SWIR-II
2017	Landsat8 OLI	4-30/7-26/8-27	C/B/G/R/NIR/SWIR-I/SWIR-II
2018	Landsat8 OLI	4-17/4-8/4-8	C/B/G/R/NIR/SWIR-I/SWIR-II
2019	Landsat5 TM	8-18/8-25/8-25	C/B/G/R/NIR/SWIR-I/SWIR-II

Note: C (Coastal); B (Blue); G (Green); R (Red); NIR (Near Infrared); SWIR-I (Shortwave Infrared-I); SWIR-II (Shortwave Infrared-II).

. In addition, the administrative area shape file of Wuhan City is also collected for LULC mapping.

Furthermore, 15 natural and socioeconomic driving factors were collected as the auxiliary data for land use simulation of Wuhan. These driving factors involve point data, line data and area data, which contribute to the simulation from different aspects. Specifically, the driving factors include population, GDP, soil, DEM, annual mean temperature, annual precipitation, the distance from highway, railway, trunk/primary/second road, high-speed railway stations, towns, provincial centers and prefectural centers. In practice, the distance can be extracted by Euclidean distance analysis. Due to the similarity of the study area, most factors were same with Reference [36]. Meanwhile, the distances from the administrative centers were refined in this research study, with some factors with less influences in Reference [36] removed. The detailed information of driving factors is provided in

Table 2. The introduction of the driving factors.

Category	Data	Original Resolution	Data Resource
Socioeconomic driving factors	Population	1000 m	http://www.geodoi.ac.cn/WebCn/Default.aspx (accessed on 3 March 2021)
	GDP
	The distance from highway	30 m	OpenStreetMap (https://www.openstreetmap.org/) (accessed on 3 March 2021)
	The distance from railway
	The distance from trunk road
	The distance from primary road
	The distance from second road
	The distance from high-speed railway stations
	The distance from provincial centers
	The distance from prefectural centers
	The distance from towns
Natural driving factors	Soil	1000 m	HWSD v1.2 (http://westdc.westgis.ac.cn/data/844010ba-d359-4020-bf76-2b58806f9205) (accessed on 4 March 2021)
	DEM	30 m	NASA SRTM1 v3.0
	Annual mean temperature	1000 m	WorldClim v2.0 (http://www.worldclim.org/) (accessed on 4 March 2021)
	Annual precipitation

.

3. Method

In this study, the continuous time series LULC maps for Wuhan during the period 2000–2019 were firstly produced with Landsat images, followed by a detailed analysis to learn the spatio-temporal patterns of LULC change and the Wuhan development mode. Meanwhile, the land use structure in 2029 was predicted by the PLUS model, with the future landscape dynamic patterns further explored. In summary, the workflow of our study can be depicted in Figure 2.

3.1. Data Pre-Processing

Systematic pre-processing was conducted before LULC mapping, and change analysis was conducted to eliminate errors and establish direct connections between different data. Firstly, the projection systems of all Landsat images and the administrative area shape file were geographically unified to the Universal Transverse Mercator (UTM) project system with the 1984 World Geodetic System (WGS) datum. Due to atmospheric influences, the digital numbers generally deviated from the real spectral responses of ground objects, which hinders the fine information extraction. Hence, an atmosphere correction was conducted for each image by using the FLAASH module of ENVI 5.3. Then, the overall study area was generated by mosaicking three adjacent images (Path/Row: 122/39, 123/38, 123/39) and masking the result with the administrative division.

Although the images were selected under the principle of cloud-free or low-cloud coverage, cloud contamination was still inevitable for some special years, which blocks the underlying surface information and degrades the mapping quality to a large degree. In order to relieve this problem, the potential clouds contained in the image were firstly masked by the general cloud index [38]. It takes full advantages of cloud reflection properties to construct the following two indices:

$$C{I}_1=\frac{B_{NIR}+2\times B_{SWIR-1}}{B_B+B_G+B_R},C{I}_2=\left(B_B+B_G+B_R+B_{NIR}+B_{SWIR-1}+B_{SWIR-2}\right)/6$$

(1)

where B_B, B_G, B_R, B_NIR, B_SWIR-1 and B_SWIR-2 denote the blue, green, red, near infrared, short wave infrared-I and short wave infrared-II band, respectively. Then, with the recommended thresholds for CI₁ and CI₂ [38], the potential cloud pixels can be effectively identified. Then, based on the cloud detection results, the cloud pixels were filled with the corresponding clear-sky pixels obtained from time series observations, according to the time proximity.

Furthermore, due to different sources, the 15 driving factors have various geographical systems and spatial resolutions and the cloud is not directly utilized for land use simulation. Hence, they were geographically unified with the Landsat images in order to possess the same geographical system and projection system, i.e., WGS84 and UTM, with the spatial resolution uniformly resampled to 30 m.

3.2. LULC Mapping

3.2.1. Feature Extraction and Classification

A six-type land system was adopted for this study, including built-up land, cropland, water, forest, grassland and unused land, according to the first classification standard of the Chinese Academy of Sciences. The spectral signatures of pixels were utilized as basic features. Specifically, six multispectral bands were utilized for Landsat TM/ETM+ images, with seven bands utilized for OLI images, as listed in Table 1. Moreover, three typical spectral indices, i.e., normalized difference vegetation index (NDVI), normalized difference water index (NDWI) and normalized difference built-up index (NDBI), were also extracted to enhance the separability between different classes, which have been verified to have strong indications for special land types, i.e., vegetation, water and impervious surfaces [39,40].

Supervised classification was employed to perform a pixelwise classification for each year over the study period. A famous machine learning algorithm, i.e., support vector machine (SVM) [41,42], was utilized for this study. It performs classification by searching an optimal hyperplane between two classes with a maximum margin, as shown in Figure 3a. The hyperplane is generally determined by only a few support vectors, which results in outstanding performance in the case of limited training samples. As a result, SVM has been a powerful tool in remote sensing classification and applications [43,44]. Due to the complexity of imaging environment, remote sensing pixels are commonly nonlinearly separable. Hence, the nonlinear SVM model with a soft margin and kernel is utilized.

Specifically, take the binary classification problem as an example. The SVM learns an optimal hyperplane as a decision boundary as follows:

$$\underset{w,b}{\min }\frac{{\Vert w\Vert }^2}{2}+C{\displaystyle \sum _{i=1}^n{\xi }_i}s.t.y_i\left[w^T\phi \left(x_i\right)+b\right]\ge 1-{\xi }_i,{\xi }_i\ge 0$$

(2)

where $x_i\in {ℝ}^D$ and $y_i\in \left\{+1,-1\right\}$ denote to the i-th sample of D dimension and its corresponding class label; w and b are hyperplane parameters; ${\xi }_i$ is a slack variable, measuring the distance of misclassified samples to the hyperplane; C is a penalty parameter; n is the number of training samples; and $\varphi$ is a kernel function, defined as $\varphi :\kappa \left(x_i,x_j\right)=\varphi {\left(x_i\right)}^T\varphi \left(x_j\right)$, which maps features from the original space into a much higher space in order to obtain an approximately linearly separable problem, as shown in Figure 3b. This study utilized the radial basis function (RBF) as the kernel function, i.e., $e^{-\delta {\Vert x_1-x_2\Vert }^2}$, where $\delta$ is the kernel width. Based on the binary classification model above, SVM can be extended to the multi-class version through “one-against-all” or “one-against-one” strategies [45].

In this study, the sample database for classification was firstly constructed through visual interpretation by experienced professionals using the region of interest (ROI) tool to label adequate representative samples for each class, with the high-resolution images and public LULC products, e.g., land use/cover data [46], as references. Then, 20% of the labeled samples were randomly selected for training, with the rest utilized for testing. With the help of SVM, a primary LULC result cloud was obtained for each year over the study period 2000–2019.

Then, two common quantitative metrics, i.e., overall accuracy (OA) and Kappa statistic, were utilized in order to assess mapping accuracy. The OA describes the proportion of correctly classified pixels in the result, and Kappa statistic considers the bias of different categories in classification and provides the consistency evaluation which can better measure the classification quality. Both these two metrics are based on the classification confusion matrix and can be calculated according to the following rules:

$$OA=\frac{p_c}{n},kappa=\frac{p_0-p_e}{1-p_e}$$

(3)

where p_c denotes the number of correctly classified samples, n is the total number of samples, p₀ is the sum of the number of correctly classified samples for each category divided by the total number of samples and $p_e=\left(a_1\times b_1+a_2\times b_2+\cdots a_k\times b_k\right)/\left(n\times n\right)$, where a_i and b_i are the number of real samples and the number of predicted samples for each category and k is number of classes.

3.2.2. Post-Processing

The SVM classification results were then post-processed to correct misclassifications in order to obtain the refined LULC maps before change analysis. The post-processing mainly involved temporal consistency analysis and spatial coherence analysis. Firstly, based on the time-series maps, the LULC conversion logic was verified for each pixel to correct the irrational labels. As we know, a continuous closed-loop change is generally unreasonable. For example, a cropland pixel is not permitted to change into built-up land but then it changed back to cropland in the next year; this is most likely to be a misclassification. Based on this fact, the irrational labels were corrected according to the following rule:

$$\mathrm{if}{y}_{t-1}=i\mathrm{and}{y}_t=j\mathrm{and}{y}_{t+1}=i,\mathrm{then}{y}_t=is.t.i,j\in \left\{1,2,3,4,5,6\right\}$$

(4)

where $y_t$ denotes the label of the target pixel in the t-th year, with 1, 2, 3, 4, 5 and 6 standing for built-up land, water, forest, cropland, unused land and grassland, respectively. Due to the inconsistency of image size, two time series were constructed for temporal consistency analysis in this study, i.e., {2000–2009, 2013–2015, 2017} and {2010, 2011, 2018, 2019}. In addition, considering the peculiarity in 2016 when a serious flood broke out, the LULC result of 2016 was reserved separately.

Then, the spatial coherence check was conducted for each LULC map. According to the first law of geography, land use should gradually change in the space. That is to say that the pixels in a small neighborhood belong to the same land use type in a large probability. It means that the isolated heterogeneous pixels are commonly caused by misclassifications. Specifically, majority voting was employed to eliminate such pixels and to enhance the homogeneity of the LULC results, which is a common strategy in remote sensing applications [47,48]. Specifically, it opened a small spatial window for each pixel and took the land use type with the highest frequency as the final class label of the central pixel. Considering the relatively low spatial resolution of Landsat images, i.e., 30m, a local window at a size of 3 × 3 was selected as the unit for majority voting.

3.3. LULC Change Analysis and Modeling

Based on the produced continuous time series LULC maps, we comprehensively explored the spatio-temporal patterns of LULC change in Wuhan both qualitatively and quantitatively, with several typical metrics utilized to model the landscape dynamics, including single land use dynamicity, landscape activity, urban gravity center and transition distance.

Specifically, the single land use dynamicity describes the change extent and speed of certain land use over a given period [34], which is widely used in LULC change analysis. It can be calculated according to the following formulation:

$$K=\frac{A_j-A_i}{A_i}\times \frac{1}{T}\times 100\%$$

(5)

where K denotes the change velocity of the target land use with a positive value representing an increasing rate and a negative value representing a decreasing rate; T is the time interval; and A_i and A_j stand for the area of the target land use in the i-th and j-th year, respectively. In this study, this index was calculated on the scale of the whole city and the scale of the inner administrative district to comprehensively uncover the landscape dynamic characteristics of Wuhan.

A landscape activity index was defined to explore the internal imbalance and regional difference between the 13 administrative districts. It describes the total dynamicity of a certain region over all land uses, i.e., the ratio of the changed area to the total area of a given district, which can be formulated as follows:

$$D_i=\frac{A_c}{A_{p,i}}\times 100\%$$

(6)

where D_i denotes the landscape activity of the i-th district, with a larger value indicating a more active state; and A_c and A_p,i stand for the changed area and the total area of the i-th district in the p-th period, respectively. Compared with single land use dynamicity, this index removes the influence of base area of different land uses and comprehensively describes the landscape activity and vitality of each district.

The urban gravity center describes the equilibrium point of the spatial distribution of built-up land, which is an important index for urban evolution analysis [49]. Generally, it is defined as the weighted average of coordinates of urban block geometric centers and can be calculated as follows:

$$X_t=\frac{{\displaystyle {\sum }_{i-1}^N\left(A_{ti}\times X_i\right)}}{{\displaystyle {\sum }_{i-1}^N{A}_{ti}}},Y_t=\frac{{\displaystyle {\sum }_{i-1}^N\left(A_{ti}\times Y_i\right)}}{{\displaystyle {\sum }_{i-1}^N{A}_{ti}}}$$

(7)

where X_t and Y_t denote to the coordinates of urban gravity center in the t-th year; X_i and Y_i stand for the coordinates of geometric center of the i-th urban block; and A_ti is the area of the i-th urban block, with N being the number of urban blocks. With the help of urban gravity centers, the urban sprawl direction and historical evolution track can be effectively inferred. Moreover, the transition distance of urban gravity center was also calculated to quantitatively analyze the urbanization intensity in a given period, which can be formulated as follows:

$$d=\sqrt{{\left(X_{t2}-X_{t1}\right)}^2+{\left(Y_{t2}-Y_{t1}\right)}^2}$$

(8)

where d denotes to the transition distance of urban gravity center in period t₁–t₂; X_t1, Y_t1 and X_t2, Y_t2 represent the coordinates of urban gravity centers in the t₁-th and t₂-th year, respectively.

3.4. LULC Dynamic Simulation

In order to better understand the land dynamic system of Wuhan and explore its future dynamic pattern, the recently developed advanced PLUS model [36] was employed to predict the land allocation in 2029. It is derived from the well-known cellular automata (CA) model. As is well known, CA is a dynamic model focusing on the spatio-temporal interactions and potential rules. With bottom-up modeling, powerful parallel computing and prominent spatio-temporal dynamic expression abilities, CA is very suitable for simulating the complex geographical evolution process and has been an important and mature tool in land use simulation studies [50,51]. The PLUS model enhances the prediction ability of CA by better integrating the impacts of various spatial factors with the dynamics of geographical cells. It introduces a new rule-mining strategy and a patch-generating mechanism to better learn the nonlinearity in LULC changes. Consequently, it generally results in higher accuracy than other state-of-the-art models [36]. Specifically, the mechanism of the PLUS model fits the task in this research well to simulate the land use dynamics by integrating multiple spatial driving factors based on multi-temporal LULC results. Therefore, the PLUS model was employed to predict the land use pattern of Wuhan in 2029.

The PLUS model involves two basic modules, i.e., a land expansion analysis strategy (LEAS)-based rule-mining module and a multi-type random patch seeds-based CA (CARS) module. The LEAS module excavates the impact of each driving factor on the expansion of geographical units based on random forest and derives the final development potential of each land use in the space. The CARS module is a special CA model, which combines the influences of “top-down” and “bottom-up” mechanisms on the land system. It simulates the local micro land use competition to drive the total land use to meet the macro future demand under the comprehensive action of adaptive coefficient, neighborhood effect and development probability. The neighborhood effect of a land use unit is determined by the proportion of the number of units in its own category and the priori diffusion coefficient. The adaptive coefficient is influenced by the difference between the number of units in its own category and the future demand, while the development probability is derived from the mining results of LEAS. For detailed information, please refer to [36]. Generally speaking, the PLUS model takes advantage of transition analysis and pattern analysis to provide consideration to effectiveness and efficiency, and the software can be obtained from http://www.urbancomp.net/2020/07/25/plus/ (accessed on 11 March 2021).

Specifically, the PLUS model needs two LULC maps as basic input for land use simulation, together with multiple concerning driving factors to learn the land change rules. In this study, 15 different factors were collected to assist LULC simulation, as described in Sub-Section 2.2. The simulation accuracy was firstly tested over a past time interval, i.e., 2009–2019, with 2009 and 2014 as basic dates for an equal interval prediction. Then, the LULC in 2029 was predicted for future land use dynamic analysis.

4. Results

4.1. LULC Maps and Accuracy Assessment

The continuous time series LULC maps of Wuhan in the period of 2000–2019 are given in Figure 4. The LULC results depicted rapid urbanization in Wuhan over the past 19 years, with large landscape changes taking place. The urban area expanded rapidly similar to a large flatbread, with a large account of cropland occupied, which greatly changed the land use structure of Wuhan. In addition, the forest and grassland also showed some changes, especially for the northwestern and central parts of Wuhan. The change in the area covered by water was relatively small over the study period, except for 2016.

Table 3. Accuracy assessment of the LULC maps of Wuhan during the period 2000–2019.

Year	Original Classification		Temporal Consistency		Spatial Coherence
	OA (%)	Kappa	OA (%)	Kappa	OA (%)	Kappa
2000	86.17%	0.8218	86.17%	0.8218	88.31%	0.8490
2001	79.94%	0.7415	82.97%	0.7789	84.86%	0.8032
2002	86.96%	0.8318	88.19%	0.8469	89.54%	0.8641
2003	85.25%	0.8142	86.19%	0.8250	87.86%	0.8459
2004	84.31%	0.7938	86.85%	0.8250	89.04%	0.8539
2005	90.36%	0.8687	91.27%	0.8806	92.87%	0.9022
2006	87.61%	0.8368	89.89%	0.8663	91.58%	0.8885
2007	89.11%	0.8561	90.46%	0.8731	91.91%	0.8923
2008	87.44%	0.8354	88.76%	0.8520	90.49%	0.8744
2009	87.66%	0.8369	88.23%	0.8438	89.66%	0.8626
2010	82.88%	0.7763	82.88%	0.7763	85.52%	0.8098
2011	84.09%	0.7943	83.61%	0.7877	85.77%	0.8153
2012	88.56%	0.8530	88.56%	0.8530	89.44%	0.8638
2013	85.58%	0.8122	86.50%	0.8229	88.04%	0.8427
2014	89.27%	0.8611	90.21%	0.8727	91.66%	0.8913
2015	88.40%	0.8510	87.54%	0.8384	89.27%	0.8606
2016	87.02%	0.8261	87.02%	0.8261	88.46%	0.8449
2017	84.67%	0.8001	84.67%	0.8001	86.56%	0.8244
2018	87.09%	0.8247	86.73%	0.8194	88.45%	0.8425
2019	89.60%	0.8632	89.60%	0.8632	91.09%	0.8826

gave the quantitative evaluations of the LULC mapping results. For all years, the OA and Kappa exceeded 84% and 0.80, with most of them being larger than 88% and 0.84, which meets the requirements of the subsequent analysis according to the 80% rule of OA in remote sensing applications [52]. In addition, it can be found that the mapping accuracy was improved by temporal consistency analysis for most years, i.e., 2001–2009, 2013 and 2014, and was improved by spatial coherence analysis for all years, which demonstrates the effectiveness of the post-processing operators adopted in this study.

Moreover, Figure 5 intuitively shows the function of post-processing with two typical scenes. From the figure, it can be observed that misclassifications caused by imaging environment and seasonal changes, e.g., illumination change and agricultural fallow, were effectively corrected to a large extent by using temporal consistency analysis, and the scattered heterogeneous pixels were smoothed out to obtain a more homogeneous result through spatial coherence analysis. Accordingly, more reasonable LULC results were obtained after post-processing.

4.2. Temporal Change Analysis by Land Use Theme

Based on the continuous time series LULC maps, yearly change of each land use was firstly analyzed to discover the change trend and characteristics. For convenient analysis, the highly correlated land uses were explored together.

4.2.1. Built-Up Land Expansion and Cropland Consumption

Over the past 19 years, the built-up land and cropland have both greatly changed under the background of urbanization. Figure 6 provides the detailed changes in areas of built-up land and cropland from which an obvious growth trend can be observed for built-up land, with a declining trend for cropland.

Specifically, the built-up area grew year by year, increasing from 430.62 km² in 2000 to 1413.28 km² in 2019, with an increment of 982.66 km² (228%). That is to say, the urban scope tripled during the past 19 years, which indicates significant urbanization. In addition to the advantages, the large-scale growth of impervious surface was bound to exacerbate urban heat island and imposes pressure on public health. As the most direct consumption of urban sprawl, the cultivated land area decreased from 6078.85 km² in 2000 to 5361.71 km² in 2019, with a reduction of 717.14 km² (12%). Undoubtedly, this results in a certain decline of the local land productivity and threatens food security to some degree.

In addition, an obvious abnormal phenomenon could be noted for cropland in 2016. The plunge was caused by the serious flood caused by the concentrated heavy rainfall in summer. The breakwater of Tangxun Lake and many other lakes resulted in the inundation of large areas of farmland, which can be obviously observed from Figure 4. With the cast off of flood, the farmland area rose rapidly in 2017.

4.2.2. Natural Habitat Dynamics

The natural habitat is an important component of ecological system and has a critical impact on local climate and environment. In this study, the natural habitat involves water, forest and grassland, with yearly changes depicted in Figure 7.

From the figure, it can be observed that the natural habitats were also affected by urbanization to some extent. Different from built-up land and cropland, there is no monotonous change trend for any kind natural habitat, instead, it alternately increases and decreases. Generally speaking, the natural habitats shrank to some extent over the study period. Specifically, the water area reduced from 1052.98 km² in 2000 to 870.46 km² in 2019, with 182.52 km² (17.33%) occupied by other land uses. There was an obvious anomaly in 2016 where the water area reached an unprecedented high value because of the flood disaster. The forest area decreased from 678.43 km² to 654.51 km² over the past 19 years, with the grassland area decreasing from 226.01 km² to 161.06 km². All these reductions created a threat to the local ecological environment and hindered urban sustainable development to some degree.

4.2.3. Land Use Conversions

In order to further explore the interactions among different land types and to quantitatively analyze their conversions induced by urbanization over the past 19 years, the land use conversion matrix in period 2000–2019 is provided in

Table 4. Land use conversion matrix of Wuhan in period 2000–2019.

	2019	Built-Up Land	Water	Forest	Cropland	Unused Land	Grassland
2000
Built-up land		305.08	7.73	5.46	99.05	4.83	8.40
Water		55.13	707.46	20.92	248.92	1.92	17.94
Forest		42.13	42.75	229.46	348.52	4.61	10.74
Cropland		947.37	103.45	349.69	4463.39	100.04	113.86
Unused land		34.77	7.64	4.24	61.10	5.42	2.54
Grassland		28.72	0.97	44.57	139.74	4.82	7.14

.

From Table 4, it can be observed that certain conversions have taken place under the background of urbanization during the past 19 years. There were 947.37 km² cropland changing to built-up land as the consumption of urban sprawl. By comparison, the conversion area of other land use to built-up land was relatively small. Specifically, the areas of water, forest, unused land and grassland converted to built-up land were 55.13 km², 42.13 km², 34.77 km² and 28.72 km², respectively. In addition, it can be observed that the conversions between forest and cropland were also significant. Generally speaking, all these conversions resulted in a large change of land use pattern in Wuhan.

4.3. Spatio-Temporal Pattern Analysis of LULC Change

Five typical years (2000, 2004, 2009, 2014 and 2019) with representative LULCs were selected as time nodes to generate four typical periods, i.e., 2000–2004, 2004–2009, 2009–2014 and 2014–2019, at approximately five-year intervals. Based on this, we explored the spatio-temporal patterns of LULC change in Wuhan from multiple views.

4.3.1. Single Land Use Dynamicity

Firstly, the single land use dynamicity of Wuhan City was explored to characterize the overall dynamic characteristics of the landscape during the past 19 years, with the results provided in

Table 5. Single land use dynamicity of Wuhan City in different periods.

Land Use Dynamicity (%)	2000–2004	2004–2009	2009–2014	2014–2019	2000–2019
Built-up land	7.01	5.15	9.09	8.03	12.01
Cropland	−0.28	−0.71	−0.66	−0.87	−0.62
Water	0.34	0.44	−2.12	−2.16	−0.91
Forest	4.03	0.66	−1.70	−2.42	−0.19
Grassland	−13.20	5.84	0.81	2.45	−1.51
Unused land	−12.12	−2.83	21.29	3.02	0.27

.

From Table 5, it can be observed that different land uses presented distinct dynamic patterns. Specifically, the built-up land area was constantly increasing, with an average annual growth rate of 12.01%, which resulted in a typical urbanization course that sped up first, then slowed down and then accelerated again. The cropland area continuously declined with an aggravated velocity. However, due to the large base area, its dynamicity stuck at a low level. Generally, the cropland was consumed at an average annual rate of 0.62%. For natural habitats, i.e., water, forest and grassland, the area presented a fluctuating state but reflected a decreasing trend in general, with annual declining rates of 0.91%, 0.19% and 1.51%, respectively.

Then, the single land use dynamicity was explored at the administrative district level to uncover the regional differences of LULC changes in Wuhan, with the results provided in Figure 8. For convenient analysis, the forest and grassland, which have similar ecological functions, were merged into one class, i.e., vegetation. Affected by the base area, the dynamicity was small for cropland while it was large for vegetation.

Generally, the built-up land dynamicity occupied the dominant position for most districts over the whole study period, i.e., 2000–2019, in the context of urbanization, as shown in Figure 8e. However, the land use dynamic pattern varied at different stages. Specifically, Hongshan, Hanyang, Jiangxia and Dongxihu showed a larger built-up land dynamicity at most stages (the early, middle-to-late and late stages), i.e., 2000–2004, 2009–2014 and 2014–2019. Especially for Hongshan, owing to the preferential policies, abundant resources and extensive space, this district developed rapidly, with the built-up land dynamicity maintaining at a high level at all stages over the study period. Additionally, Huangpi, Xinzhou, Caidian and Hannan showed a relatively high dynamicity for built-up land as well at certain stages, but it was more due to the smaller base area. In particular, Huangpi, Caidian and Hannan showed a high built-up land dynamicity at the middle-to-late and late stages, i.e., 2009–2014 and 2014–2019, while Xinzhou showed a high built-up land dynamicity at the early and middle-to-late stages, i.e., 2000–2004 and 2009–2014. However, for most old metropolitan districts, such as Wuchang, Jiang’an, Jianghan, Qiaokou and Qingshan, the built-up land dynamicity was not significant because of the greater base area. By comparison, the dynamicity of cropland and natural habitats was more obvious.

4.3.2. Comprehensive Landscape Activity

The landscape activity index was utilized to further evaluate the dynamic pattern of the 13 administrative districts in order to better understand the development mode of Wuhan, with the results provided in Figure 9.

From Figure 9, it can be clearly observed that there exists an obvious internal imbalance and regional difference between the 13 administrative districts in the development process of Wuhan. From the whole research stage, Hanyang, Hongshan, Jiang’an and Dongxihu presented a more active landscape, with Hanyang being the active center. Among them, most are metropolitans. By comparison, the suburbs, such as Huangpi, Xinzhou and Hannan, lacked vitality and had a less active landscape.

On the other hand, it can be observed that the active pattern varied a lot at different stages, with the active districts frequently changing. Specifically, at the early stage, i.e., 2000–2004, Jiang’an, Qingshan and Dongxihu owned the largest activity, followed by Hanyang, Hongshan, Caidian, Huangpi and Xinzhou. At the middle stage, i.e., 2004–2009, Hanyang started to be the active center, with the landscape of Jiang’an, Qingshan and Dongxihu still being very active, while the landscape activity of Huangpi and Xinzhou obviously declined. At the middle-to-late stage, i.e., 2009–2014, the landscape of Hanyang, Dongxihu and Caidian was still very active, followed by that of Wuchang and Hongshan. At the late stage, i.e., 2014–2019, the large activity remained in these three districts, while the landscape of Hannan and Jiangxia started to be active.

4.3.3. Land Dynamics of Typical Districts

In this section, four representative districts with different characteristics were selected, i.e., Jianghan, Hongshan, Jiangxia and Huangpi, in order to further uncover the landscape dynamic difference between metropolitans and suburbs, with an old metropolitan, a new metropolitan, an urban fringe suburb and a suburb included. The LULC results of these four districts at five typical time nodes (2000, 2004, 2009, 2014 and 2019) are provided in Figure 10.

From the figure, it can be observed that the land use structure and major land dynamics were quite different for metropolitans and suburbs. For metropolitans, i.e., Jianghan and Hongshan, the built-up land area accounted for a considerable proportion, and its dynamicity was significant and dominant over the study period. Especially for the rising metropolitan, i.e., Hongshan, the built-up land dynamicity was conspicuous, with its proportion increased from 5.75% in 2000 to 32.91% in 2019. On the other hand, for suburbs, i.e., Jiangxia and Huangpi, most of the area was covered by the arable land and natural habitats, with built-up land change submerged into other changes. By contrast, as an urban fringe area and the natural extension of urbanization, Jiangxia showed a larger development potential than Huangpi, with the built-up area proportion increased from 1.25% to 9.71% (nearly eightfold) over the past 19 years.

In addition, the urbanization pattern was quite different for metropolitans and suburbs. For metropolitans, urbanization presented a rapid multi-directional expansion in the form of a big pie at anywhere that was available. While for the suburbs, in addition to the slow sprawl of original towns, urbanization mainly took place close to the metropolitans, such as the south of Huangpi and the northwest of Jiangxia.

4.3.4. Urbanization Mode and Land Consumption

In order to discover the land consumption of urbanization in Wuhan, we explored the situation of urban land crowding out other land uses in different periods, with the results provided in Figure 11.

From the figure, it can be observed that the main consumption of urban expansion is the cultivated land, accounting for the vast majority in each period. For the other land uses, the consumed land types were different in various periods. For example, there was a relatively obvious consumption for water during the period 2000–2014, while it was obvious for unused land during the period 2014–2019. For the natural habitats, i.e., forest, grassland and water, the consumed part mainly concentrated in the regions nearer to the city, which suggested that their growth was greatly affected by human activities.

In order to better learn the urbanization mode of Wuhan, we comprehensively explored the urban sprawl pattern, gravity center transition and development direction, with the results given in

Table 6. Urban gravity center transition distance during different periods.

Period	2000–2004	2004–2009	2009–2014	2014–2019
Transition distance (km)	4.10	0.38	1.73	4.46

and Figure 12 and Figure 13.

From Figure 12 and Table 6, it can be clearly observed that the urban sprawl pattern and gravity center transfer route varied in different periods. Specifically, during the period 2000–2004, the urban sprawl in the east and south was obvious. As a result, the gravity center transferred towards the southeast, with a large transition distance of 4.10 km. During period 2004–2009, the urbanization course was not significant. Consequently, the gravity center only slightly moved. During period 2009–2014, urban expansion began to accelerate, with the gravity center shifted to the east at a relatively large step, i.e., 1.73 km, while during period 2014–2019, the urbanization in the southwest was prominent. Accordingly, the gravity center transferred in a quite different direction, i.e., the southwest, with a large transition distance of 4.46 km.

From Figure 13, it can be observed that the principal development direction varied at different stages, but with similar trends as above. Specifically, during period 2000–2004, the principal direction was the southeast. During period 2004–2009, the development direction was not obvious, with the southeast slightly dominant. During period 2009–2014, the principal direction was the northeast, while in 2014–2019, the principal direction changed to the southwest. On the whole, Wuhan mainly developed towards the east and south.

4.4. Simulated LULC Results

Figure 14 provides the land use simulation results of the PLUS model. From the figure, it can be observed that the simulation result in 2019 (Figure 14b) was basically consistent with the reference map (Figure 14a), which demonstrates the effectiveness of the PLUS model. The composition and distribution of various land uses approximately coincide with the reality. As a result, an OA of 74.69% and a figure of merit (FOM) of 0.18 were obtained, which reaches a good simulation standard and can guarantee the subsequent analysis [36].

From the predicted LULC result in 2029 (Figure 14c), it can be observed that a large landscape change will occur in the next 10 years. The land use dynamics will still be significantly induced by urbanization. Specifically, the built-up land obviously expands outside, while the cropland has a large proportion of reduction as the consumption of urbanization, especially for the regions close to the metropolitan, further exacerbates the threat to food security. The unused land has no obvious change, and the forest, grassland and water seem to have some incremental change. For instance, there is an obvious growth of forest in the northeastern part of Wuhan.

Figure 15 provided the importance analysis of the 15 driving factors utilized in this research study for urban expansion. From the figure, it can be observed that DEM contributed most to built-up land expansion. This is because topographic factors determined the difficulty and cost of urban construction. Generally, urban sprawl will avoid mountainous areas with large topographic fluctuations. From Figure 15a, it can be observed that the expanded built-up areas were mostly distributed in the relatively flat part of Wuhan, which is in line with the general law of urban development. In addition, it can be observed that the distance from prefectural centers and the distance from primary/second roads were also dominant driving factors in the growth of built-up land. The reasons could be that the distribution of administrative centers was generally an important factor in urban planning. Therefore, the distance from prefectural centers will naturally affect the direction and extent of urban expansion, while road is an important part of traffic and is an important prerequisite for urban development. Thus, the distance from primary/second roads had a large impact on the growth of built-up land.

Figure 16 and Figure 17 provided future land use dynamics and urbanization mode based on the predicted LULC result in 2029. From the figures, it can be observed that the urban sprawl will continue in the next 10 years, with the “converted into” area much larger than the “converted from” area. The built-up land continues to rapidly expand outward along the two rivers, especially for Hongshan, Hanyang, Jiangxia and Caidian. Unlike before, the urban area will mainly develop toward the northeast in the future, with more significant urbanization in Huangpi and Xinzhou. Accordingly, the urban gravity center will transfer to the northeast, with a relatively large transition distance of 2.62 km. As a result of urbanization, a large quantity of farmland will continue to be occupied. Unlike the past 19 years, the natural habitats, i.e., forest, grassland and water tend to increase in the future, which is a good signal for ecological environment improvement and urban sustainability.

5. Discussion

5.1. Land Use Structure and Dynamicity Analysis

5.1.1. Land Use Structure Characteristics

The land use structure of Wuhan is very distinct. Specifically, the urban area is mainly distributed along the two rivers, i.e., the Yangtze River and Han River, forming a typical pattern of three towns (Wuchang, Hankou and Hanyang) across the rivers. The built-up lands are mainly concentrated in the seven metropolitan districts in central Wuhan, with a scattered distribution in the outer suburbs. Over the past 19 years, the urban area constantly spread to the periphery, similar to a pie. As a result, great changes took place in the urban scope and form, from concentration in the hinterland of the Yangtze River during the initial stage to all-around development along the two rivers to four banks during the current stage. However, as a central city, the dominant land use of Wuhan is still arable land, accounting for nearly 60%, which is mainly distributed in the six suburban districts, such as the south of Huangpi, most of Xinzhou, Dongxihu and Hannan, the northwest of Caidian and the west of Jiangxia, and occupied the vast majority of their area.

On the other hand, Wuhan enjoys a good reputation of “the city of hundreds of lakes”. Water occupies a considerable area, with large amounts of rivers, lakes, ponds and ditches in different sizes in addition to the Yangtze River and Han River, which is embedded in the separated area of the two rivers and forms a developed water network. For the other two natural habitats, forest and grassland had a relatively concentrated distribution, mainly in the north of Hangpi, the east of Xinzhou and the middle of Hongshan, Jiangxia and Caidian, and were less affected by human activities.

5.1.2. Land Dynamic Characteristics

Driven by a series of national strategic measures, such as the rise of central China, Wuhan experienced unprecedented urbanization over past 19 years, which can be characterized by a distinct process that sped up first, then slowed down and then accelerated again (Table 5). Consequently, great landscape changes were triggered. The urban area constantly increased at an annual rate of 12%, with an increment of 228%. The immediate consequence of the rapid urban sprawl was a great reduction in cultivated land, i.e., 12%, which poses a threat to local land productivity and food security. Meanwhile, certain fluctuations and reductions have took place in natural habitats, which brings some negative effects to urban sustainability.

On the other hand, there was an obvious internal imbalance and regional difference between the 13 administrative districts in the development process of Wuhan (Figure 8 and Figure 9). On the whole, Hanyang, Hongshan and Dongxihu owned a more active landscape than others, with Hanyang being the active center. Meanwhile, the built-up land dynamicity was significant and dominant for most districts over the study period.

Specifically, for most old metropolitan districts, such as Wuchang, Jianghan, Qiaokou and Qingshan, the built-up land dynamicity was not obvious and the landscape lacked activity. This was because the urban construction had already been basically completed at the end of last century, with limited preparation space for urban development. The land dynamics mainly came from land structure upgrading, urban greening, infrastructure renewal and old city reconstruction. By comparison, the dynamicity of natural habitats and cropland was more obvious. However, as an old metropolitan, Hanyang was an exception. It was because of the backward foundation and large development space. Driven by a series of preferential policies and the real estate, Hanyang had joined in the main battlefield of urban construction over the past 19 years and became the active center of Wuhan. In addition, as an emerging metropolitan, Hongshan demonstrated large potential and stamina owing to the advantages of policy, resources and space. Since 2000, with a series of policies enacted, such as the planning of China Optical Valley and the construction of high-tech industrial bases in the Donghu new technology development zone, Hongshan ushered in an era of rapid development. The landscape has been very active, with the built-up area risen by nearly six-fold.

For the suburb districts, i.e., Huangpi, Xinzhou, Dongxihu, Jiangxia, Caidian and Hannan, due to the small base area, their built-up land dynamicity was obvious. This was especially true for the urban fringe districts, such as Dongxihu, Caidian and Jiangxia, which were natural extensions of urban development. Frequent conversions happened between built-up land and farmland or other land uses. As a result, the landscape of Dongxihu and Caidian was very active in most stages over the study period.

5.2. Urbanization Mode Analysis

5.2.1. Analysis of Urban Land Crowding out Other Land Uses

Under the background of urbanization, great landscape changes have taken place in Wuhan over the past 19 years, which brought about frequent land use conversions. Here, we focused on the land expenditure of urbanization, i.e., the situation of urban land crowding out other land uses (Figure 11).

Overall, the cropland to built-up land conversion dominated during various periods, with a small amount of natural habitat to built-up land conversions. In particular, the arable land around the city accounted for the land expenditure of urbanization to the greatest extent, which is in line with the general law of China. Specifically, the resulting consumption of arable land reached as high as 200.10 km², 234.34 km², 360.45 km² and 569.17 km² during the periods 2000–2004, 2004–2009, 2009–2014 and 2014–2019, respectively, with a total expenditure of 947.37 km² over the study period, which indicates the cultivated land occupation development mode for Wuhan. In addition, there was also a small number of natural habitats crowded by urban land, especially for the places that are close to human activities.

5.2.2. Urban Sprawl Pattern and Gravity Center Transition

The urban scope and outline changed a lot over the past 19 years in the context of urbanization, resulting in a distinct urban sprawl pattern and gravity center transfer route for Wuhan (Figure 12 and Table 6). Generally, the urban are constantly expanded along the rivers and roads, with various velocities and patterns at different stages. Specifically, during the early stage, i.e., 2000–2004, the urban area rapidly expanded to the periphery, especially towards the east and south, such as Hongshan and Jiangxia. In the next stage, i.e., 2004–2009, the urban area only showed a slight sprawl in small pieces in the periphery. Relative to the urban scope in 2009, many new urban areas appeared in the fringe areas in 2014, which suggests an acceleration of urbanization, while in the late stage, i.e., 2014–2019, the urban scope had large-scale growth along the two rivers to four banks. All of these results demonstrated the typical urbanization course of Wuhan once more.

On the other hand, the urban gravity center ceaselessly shifted throughout the whole study period. Generally speaking, it mainly appeared in the central zone of Wuhan near the Yangtze River. Specifically, during the period 2000–2004, the gravity center transferred towards the southeast as a result of the rapid urbanization. It transferred from the northwest of Jiang’an to the junction of Wuchang and Jiang’an, with a very large transition distance of 4.10 km. During the period 2004–2009, the gravity center further moved to the junction of these two districts, but with a small transition distance of only 0.38 km. During the period 2009–2014, the gravity center mainly shifted to the east and transferred to the junction of Wuchang and Qingshan, with a relatively large transition distance of 1.73 km. During the period 2014–2019, the gravity center moved in a quite different direction, i.e., the southwest, and back to the interior part of Jiang’an, with a large transition distance of 4.46 km. The above transfer route illustrated that Wuhan mainly developed towards the east and the south over the past 19 years.

5.2.3. Urban Development Direction

The urban development direction (Figure 13) indicated the development priorities of Wuhan at different stages, which is an important component of the urbanization mode.

Generally speaking, the principal development direction of Wuhan varied at different stages, but with similar trends as above. Specifically, during the early stage, i.e., 2000–2004, the urban development in the southeast was particularly prominent, with much weaker urbanization in the north and west, which resulted in an obvious principal direction. In the middle stage, i.e., 2004–2009, the urbanization slowed down, with the principal direction being not obvious. By contrast, the development in the southeast seemed more prominent. During the middle-to-late stage, i.e., 2009–2014, the development in the east showed a certain advantage, especially for the northeast. During the late stage, i.e., 2014–2019, a vastly different trend could be observed: the principal direction changed into the southwest. On the whole, over the entire study period, the urban development principal direction of Wuhan was the east and the south, especially for the southeast.

5.3. Future Land Use Dynamic Pattern Analysis

5.3.1. Future Land Use Dynamics

In order to explore the relationship between future land use conversions and past trends, the future dynamic characteristics and conversion trajectories of major land uses was further analyzed (Figure 16).

The results suggested that the “converted into” area was much larger than the “converted from” area for built-up land, which demonstrates a further urban sprawl over the next 10 years. The built-up area will reach 1794.27 km² in 2029, with an increment of 380.99 km² (26.96%) compared with that of 2019. On the other hand, it can be observed that the urbanization will still mainly take place in agricultural regions, at an expense of a reduction of 378.49 km² (7.06%) for cropland. The continuous farmland occupation sounds an alarm for the red line of cultivated land protection and food security. On the plus side, for all natural habitats, i.e., water, forest and grassland, the “converted into” area is larger than the “converted from” area, which heralds an improvement of the ecological environment in the future. It may attributed to the related national strategic measures, such as “green water and mountain”, “beautiful China”, etc., and suggests that the future development mode of Wuhan tends to be eco-friendly.

5.3.2. Future Urbanization Mode

In order to explore the continuity of Wuhan development mode in the future, the urbanization sprawl pattern, gravity center transition, development direction and the situation of urban land crowding out other land uses in the next 10 years were further analyzed (Figure 17).

The results indicate a very different development direction of Wuhan in the next 10 years, i.e., the northeast. In addition to the urbanization in the periphery of the metropolitans, a large number of new built-up areas appeared in the northern suburbs, i.e., Huangpi and Xinzhou, and mainly distributed along the main traffic lines, including expressways, national/provincial roads and other trunk roads, where convenient transportation is provided for urban construction. Accordingly, the urban gravity center transferred towards the northeast but is still inside the Jiang’an district, with a large transition distance of 2.62 km. On the other hand, it can be observed that farmland will still be the major land expenditure of urbanization, which further threatens food security and hinders urban sustainable development. This should arouse the specific attention of urban planners and managers.

6. Conclusions

This paper carried out an in-depth study on spatio-temporal patterns of LULC change in Wuhan under the background of urbanization over the past 19 years, i.e., 2000–2019, on the basis of the continuous time series LULC maps produced by SVM by using Landsat observations. The results indicated a typical urbanization course that sped up first, then slowed down and then accelerated again, with the land use structure changing significantly. The built-up area tripled, with an increment of 982.66 km² (228%). Consequently, a large amount of cropland was occupied, with a reduction of 717.14 km² (12%), which poses a threat to food security. Meanwhile, the water, forest and grassland area had certain reductions of 182.52 km², 23.92 km² and 64.95 km², respectively. On the other hand, an obvious internal imbalance was found between the 13 administrative districts, with Hanyang, Hongshan and Dongxihu owing a more active landscape relative to others and Hanyang being the active center. Accordingly, the urban gravity center constantly changed, transferring from the northwest to southeast of Jiang’an district, with different principal development directions at various stages. Generally, Wuhan mainly developed toward the east and the south. Lastly, based on the simulated LULC result in 2029 by the PLUS model, the future landscape dynamic pattern was further explored, and the results illustrated a quite different urbanization mode in the next 10 years, mainly toward the northeast, which is meaningful for urban planning and management.

The large change of the underlying surface, especially for the rapid increase in the impermeable layer, will result in certain impacts on the local climate of Wuhan. For example, the rapid urbanization broke the balance of the original heat flow mode and further intensified the urban heat island effect of Wuhan. In addition, the large number of natural covers transforming into artificial surface in Caidian may account for the tornado that broke out in May 2021 in this area to some degree. Furthermore, the rapid urbanization in the northern suburbs, i.e., Huangpi and Xinzhou, in the next 10 years, may increase the risk of natural disasters. In the next work, we will explore the climate changes, ecological impacts and urban heat island effect under different development scenarios in Wuhan to promote urban sustainable development.