Spatiotemporal Change Analysis and Future Scenario of LULC Using the CA-ANN Approach: A Case Study of the Greater Bay Area, China

Land use land cover (LULC) transition analysis is a systematic approach that helps in understanding physical and human involvement in the natural environment and sustainable development. The study of the spatiotemporal shifting pattern of LULC, the simulation of future scenarios and the intensity analysis at the interval, category and transition levels provide a comprehensive prospect to determine current and future development scenarios. In this study, we used multitemporal remote sensing data from 1980–2020 with a 10-year interval, explanatory variables (Digital Elevation Model (DEM), slope, population, GDP, distance from roads, distance from the city center and distance from streams) and an integrated CA-ANN approach within the MOLUSCE plugin of QGIS to model the spatiotemporal change transition potential and future LULC simulation in the Greater Bay Area. The results indicate that physical and socioeconomic driving factors have significant impacts on the landscape patterns. Over the last four decades, the study area experienced rapid urban expansion (4.75% to 14.75%), resulting in the loss of forest (53.49% to 50.57%), cropland (21.85% to 16.04%) and grassland (13.89% to 12.05%). The projected results (2030–2050) also endorse the increasing trend in built-up area, forest, and water at the cost of substantial amounts of cropland and grassland.


Introduction
Landscape change is based on the natural ecosystem and the structure of human society. The process of landscape change results from changes in the physical environment and socioeconomic factors [1]. Fast economic development and population growth results in rapid urbanization and land-use changes. These changes significantly affect LULC dynamics and the cycle and structure of the ecosystem [2]. Anthropogenic influence on landscapes is one of the fundamental driving forces of regional LULC change mechanisms [3,4]. Human exploitation of the natural environment, such as urban expansion, fragmentation of agricultural land and green spaces, expressively disturbs the local environment [5]. Among these activities, rapid urban growth is considered the leading driver of cropland and green area losses [6], which may greatly influence climatic changes and human life [7,8]. The dynamics of socioeconomic development and industrialization have led to a drastic expansion of built-up areas in periurban regions. Such an increase in impervious surfaces affects the urban environment [9].
Global LULC trends are divided into two categories: intensification, such as deforestation, agricultural and urban expansion; and extensification, such as afforestation and land protection. Intensification refers to the rapid urban and agricultural expansion in the natural landscape, whereas extensification is the proposed strategies to control the 1.
Modeling spatiotemporal LULC patterns to analyze the magnitude and direction of change over the last 4 decades.

2.
Predicting future LULC under the influence of physical and socioeconomic factors. 3.
Identifying current LULC change intensity and potential impacts of LULC change on the spatial pattern. 4.
Analyzing the predicted LULC intensity scenario.

Study Area
The Guangdong Hong Kong Macao Greater Bay Area, also known as the Pearl River Delta (PRD), is located between 111°-115°-240° E and 22°-25° N ( Figure 1) and has a total area of approximately 56,000 km 2 , a population of over 71 million, and a GDP of over 10 trillion RMB [79]. The Guangdong Hong Kong Macao Greater Bay Area is a concept introduced by the Chinese government as a part of the 13th Five-Year Plan for Socioeconomic Development in 2017. The GBA consists of nine cities in Guangdong Province, including Zhaoqing, Foshan, Guangzhou, Dongguan, Jiangmen, Zhongshan, Zhuhai, Shenzhen, Huizhou, and the two Special Administrative Regions (SARs) of Hong Kong and Macao. As one of the largest bay areas in the world, and an economic and industrial symbol of China, the GBA experienced a dramatic LULC change phenomenon and socioeconomic development, especially in terms of the expansion of built-up land and rapid population growth. Rapid urban expansion is mainly attributed to population growth and socioeconomic development. Subsequently, the rapid increase in population and impervious surfaces resulted in many environmental changes in the region [80].

Data Collection
The LULC data for the Guangdong Hong Kong Macao Greater Bay Area from 1980,1990,2000,2010 and 2020, and the socioeconomic variables, i.e., GDP and population data , were downloaded from the Resources and Environment Science Data Center (RESDC), Chinese Academy of Sciences. China's land use remote sensing monitoring database is a multitemporal land use database covering the land area of the whole country. Data production is based on temporal Landsat (MSS/TM/ETM+/OLI) images generated through artificial visual interpretation. The data include six primary and 23 secondary

Data Collection
The LULC data for the Guangdong Hong Kong Macao Greater Bay Area from 1980,1990,2000,2010 and 2020, and the socioeconomic variables, i.e., GDP and population data , were downloaded from the Resources and Environment Science Data Center (RESDC), Chinese Academy of Sciences. China's land use remote sensing monitoring database is a multitemporal land use database covering the land area of the whole country. Data production is based on temporal Landsat (MSS/TM/ETM+/OLI) images generated through artificial visual interpretation. The data include six primary and 23 secondary land-use types, including arable land, woodland, grassland, water, residential and barren land, with a 30 m spatial resolution. To maintain the accuracy of LULC classification, random sampling methods, GPS field surveys and kappa coefficient tests are engaged. The accuracy rate of cultivated land, urban and rural areas is not less than 95%, and that of grassland, woodland and water area is not less than 90%; that of barren land is not less than 85%. Socioeconomic data are based on national subcounty statistics and the multifactor weight distribution method (www.resdc.cn, accessed on 28 July 2020). The digital elevation model (DEM) was obtained from WorldClim with a spatial resolution of 30 arc-seconds and was derived from the SRTM elevation data (https://www.worldclim.org/, accessed on 16 August 2020). The slope was calculated from the DEM, the road network was obtained from the Socioeconomic Data and Application Center (SEDAC NASA) from 1980 to 2010 (gROADS V1) and the stream network was downloaded from Diva-GIS (https: //www.diva-gis.org/, accessed on 16 August 2020). Proximity factors such as distance from roads, distance from streams and distance from the city center, were calculated using the Euclidean distance method in ArcGIS 10.4 (Table 1). After data collection, the study area was clipped and projected to WGS_1984_UTM_Zone_49N. Figure 2 represents the methodological workflow of the study.

Collection of Spatial Variables
Researchers mainly consider the physical and socioeconomic factors responsible for LULC changes, as their contribution to the LULC change mechanism is more significant. Physical factors such as topography and climate are considered the most influential factors that encourage anthropogenic activities [81]. Socioeconomic factors, such as population and GDP, have a potential influence on LULC change [82], and proximity factors, e.g., proximity to roads, distance from the city center and distance from stream network, help in analyzing the driving causes of the landscape pattern [83]. We used a collection of physical, socioeconomic, and proximity factors ( Table 1).
The MOLUSCE plugin offers some well-known methods, such as Pearson's correlation and Cramer's coefficient, to evaluate the correlation between LULC data and spatial variables. Cramer's V is a numerical measure of association that ranges from 0 to 1, where 1 indicates a 'perfect association' between LULC and the spatial driver and 0 indicates 'no association' [15]. The values are not considered ultimate, and they can only help determine whether to include transition potential modeling or not, but generally, a value higher than 0.1 is considered useful [5,84].

LULC Classification
After extraction and projection of the study area, we grouped the subtypes of land-use data to obtain five categories: urban land, rural settlements, and other construction land as built-up area; woodland, and shrubs as forest; paddy fields, dry land, and cultivated land as cropland; pasture, parks, and green spaces as grassland; and rivers, lakes, reservoirs, and canals as water (Table 2) using the reclassifying tool in ArcGIS 10.4. Furthermore, Google Earth high-resolution images from 2000 were used to interpret and verify the LULC categories. To use the data for transition potential modeling and prediction, we used resampling techniques in ArcGIS to fix the spatial resolution differences between the LULC data and spatial variables. Table 2. LULC classification scheme.

LULC Type Description
Built-up Area Impervious surface, residential, commercial and other infrastructure Forest All types of forest cover land Cropland Agricultural and farmland Grassland Parks, green spaces and pasture Water Rivers, lakes, ponds and dams

Change Analysis and Transition Potential Modeling
We used the Modules for Land-Use Change Simulation (MOULSCE) plugin within QGIS to determine the spatiotemporal changes and calculate the LULC transition between the study intervals (1980 to 1990, 1990 to 2000, 2000 to 2010, and 2010 to 2020) and produced four LULC change maps. Based on the LULC data and explanatory variables, we acquired an area change and transition probability matrix, which comprised rows and columns of landscape categories in the initial and final years. We employed the ANN multilayer perception method for transition potential modeling. The DEM, slope, GDP, population, distance from roads, distance from the city, and distance from streams were employed as explanatory variables (Table 1). These variables are commonly used for LULC change analysis and provide replicable information about the physical and anthropogenic impacts on LULC dynamics.

LULC Classification
After extraction and projection of the study area, we grouped the subtypes of landuse data to obtain five categories: urban land, rural settlements, and other construction land as built-up area; woodland, and shrubs as forest; paddy fields, dry land, and cultivated land as cropland; pasture, parks, and green spaces as grassland; and rivers, lakes, reservoirs, and canals as water (Table 2) using the reclassifying tool in ArcGIS 10.4. Furthermore, Google Earth high-resolution images from 2000 were used to interpret and verify the LULC categories. To use the data for transition potential modeling and prediction, we used resampling techniques in ArcGIS to fix the spatial resolution differences between the LULC data and spatial variables.

Change Analysis and Transition Potential Modeling
We used the Modules for Land-Use Change Simulation (MOULSCE) plugin within QGIS to determine the spatiotemporal changes and calculate the LULC transition between the study intervals (1980 to 1990, 1990 to 2000, 2000 to 2010, and 2010 to 2020) and produced four LULC change maps. Based on the LULC data and explanatory variables, we acquired an area change and transition probability matrix, which comprised rows and columns of landscape categories in the initial and final years. We employed the ANN multilayer perception method for transition potential modeling. The DEM, slope, GDP, population, distance from roads, distance from the city, and distance from streams were employed as explanatory variables (Table 1). These variables are commonly used for LULC change analysis and provide replicable information about the physical and anthropogenic impacts on LULC dynamics.

Prediction and Model Validation
The purpose of simulation models is to simplify the dynamics of composite urban structures and interpret them in an easily understandable way. We applied the CA-ANN approach within the MOLUSCE plugin for transition potential modeling and future simulation, as many researchers believe that the CA-ANN approach is more potent than linear regression [85,86]. The MOLUSCE plugin effectively computes land-use change analysis [87] and is well suited to analyzing spatiotemporal forest and land-use changes, transition potential modeling and simulating future scenarios [86].
Based on the LULC data of 2000 and 2010, explanatory variables, and transition matrix, we projected the LULC for 2020. To validate the model and prediction accuracy, the MOLUSCE plugin offers a kappa validation technique and comparison of actual and projected LULC images. In the ANN learning process, 100 iterations and a neighborhood value of 3×3 pixels, a learning rate of 0.001, 12 hidden layers, and 0.05 of momentum, were chosen to project the LULC for 2020. After obtaining satisfactory results from the model validation, we employed LULC data from 2010 and 2020 to forecast the LULC in 2030, the LULC of 2000 and 2020 for 2040, and the predicted data for 2030 and 2040 were used to forecast LULC for 2050 ( Figure 2).

Intensity Analysis
Intensity analysis is a statistical framework for expressing variations within various categories over multiple time intervals [88]. This approach analyzes a cross tabulation matrix and examines the transition variation among categories and the persistence of each category across the time intervals. Intensity analysis consists of three-level analysis: interval, category and transition. The interval-level analysis defines the time in which the overall annual change rate is comparatively fast or slow. The category-level analysis examines the changing intensity among LULC categories, and at this level, the intensity of gains and losses for each category are calculated and compared with the uniform intensity. The transition-level analysis compares the changes in terms of size and intensity during each time interval. It compares each category's gain and loss under the influence of other classes during each time interval. The interval level addresses the question of in which time intervals the annual rate of overall change is relatively fast versus slow, and the category level elaborates which categories are relatively dormant versus active in a given time interval. Based on the above two level analysis, the transition level answer the question of which transitions are intensively avoided versus targeted by a given category in a given time interval [70].
This study used an open-source Microsoft Excel program downloaded from the Intensity Analysis website (https://sites.google.com/site/intensityanalysis/, accessed on 10 August 2020). This program was developed by Safaa Zakaria Aldwaik and Robert Gilmore Pontius Jr., and is based on statistical equations for intensity analysis [70].
We further used the following equation to assess the spatiotemporal magnitude and rate of change in LULC categories: where ARC is the annual rate of change in LULC categories, A 1 and A 2 are the initial and final year areas, respectively, and t is the time interval.
During 2000-2010, the urban growth rate was 1.50%, and the decreasing rates for forest, cropland, grassland, and water were −0.05%, −0.27%, −0.24%, and −0.17%, respectively. During 2010-2020, the urban growth was 0.30%, the decreasing rate for forest and cropland was −0.04%, the grassland growth rate was 0.01%, and the water category's decreasing rate was −0.21%.   The LULC change analysis explored the spatial dynamic variations in the LULC pattern during the study period. The results from 1980 to 2020 showed a notable expansion in urban areas and a shrinking phenomenon in the forest and cropland. Figure 4 and Table 4 show the spatiotemporal area and percentage change in all LULC categories. Cropland, grassland, forest and water contributed 25.47%, 14.35%, 9.93%, and 5.09% to the built-up area, respectively. Figure 5 and Table 5 represent the inter-transition and contribution of LULC categories to changing phenomena from 1980 to 2020. Cropland was the largest contributor to the change from 1980-2020, with 25.47% to built-up area, 9.63% to water, and 1.79% to forest, while forest contributed 9.93% to built-up area, 1.58% to cropland, 12.07% to grassland and 0.73% to water. Grassland contributed 14.35% to built-up land and 7.02% to forest, and water contributed 5.09% to built-up land and 4.20% to cropland.  The LULC change analysis explored the spatial dynamic variations in the LULC pattern during the study period. The results from 1980 to 2020 showed a notable expansion in urban areas and a shrinking phenomenon in the forest and cropland. Figure 4 and Table  4 show the spatiotemporal area and percentage change in all LULC categories. Cropland, grassland, forest and water contributed 25.47%, 14.35%, 9.93%, and 5.09% to the built-up area, respectively. Figure 5 and Table 5 represent the inter-transition and contribution of LULC categories to changing phenomena from 1980 to 2020. Cropland was the largest contributor to the change from 1980-2020, with 25.47% to built-up area, 9.63% to water, and 1.79% to forest, while forest contributed 9.93% to built-up area, 1.58% to cropland, 12.07% to grassland and 0.73% to water. Grassland contributed 14.35% to built-up land and 7.02% to forest, and water contributed 5.09% to built-up land and 4.20% to cropland.

LULC Transition Analysis
The transition matrix plays an essential role in analyzing temporal changes within a set of LULC categories. The rows of the matrix table represent the categories in the initial year, while the columns show the same order of LULC categories in the final year. The diagonal entries show the size of class stability, and each off-diagonal entry represents the size of the transition from one class to different classes. Values close to 1 in diagonal entries represent the stability of a category. Researchers mostly use transition matrices to compare the temporal changes in different regions. To show how each LULC type was projected to change in our study area, we calculated the transition potential matrix with the help of the MOLUSCE plugin during the periods of 1980-1990, 1990-2000, 2000-2010, and 2010-2020 based on the existing LULC conditions and explanatory variables. Table 6 shows the transition potential matrix during 1980-1990, in which built-up area, forest, and water were the most stable categories with probabilities of 0.954, 0.982, and 0.949, respectively, while cropland and grassland decreased with transition probabilities of 0.914 and 0.937, respectively, and contributed 0.035 and 0.009, respectively, to built-up land. During 1990-2000 (Table 7), the transition values of built-up land and forest were 0.930 and 0.966, respectively, those of cropland and grassland were 0.843 and 0.859, respectively, and that of water was 0.890. Cropland and grassland again contributed 0.060 and 0.066 to the built-up area, respectively. Accordingly, during 2000-2010 (Table 8), built-up area and forest were stable with values of 0.918 and 0.949, respectively, and cropland, grassland, and water had values of 0.809, 0.761, and 0.736, respectively. Here, cropland, grassland, and water contributed 0.116, 0.123, and 0.119, respectively, of the transition to built-up land. During 2010-2020 (Table 9), built-up land and forest class areas stabilized again with transition probabilities of 0.988 and 0.984, respectively, and those of cropland and grassland were 0.958 and 0.961, respectively.
Moreover, as we used LULC data from 2000-2020 along with spatial factors for the prediction of 2040 and the transition probability matrix, the transition matrix between 2000 and 2020 (Table 10) showed that built-up area and forest were still stable with values of 0.931 and 0.936, respectively, while cropland, grassland and water had fragmentation values of 0.799, 0.742 and 0.706, respectively. During 2030-2040, cropland will still experience a decreasing phenomenon (Table 11), while other classes will show comparatively stable behavior. Table 12 represents the transition probability matrix during 2040-2050. Here, we observed that built-up land, forest, grassland, and water were the most stable categories, with transition values of 1.0, 1.0 1.0, and 0.99, respectively, while cropland will experience rapid fragmentation, with a value of 0.732. During the study period, along with the built-up area, the forest was stable because of the conservation, afforestation, and reforestation policy of the Chinese government, but the significant pressure was on cropland, which donated the largest share to other LULC categories. Table 13 and Figure 6 elaborate the gained and lost areas of each category over each time interval.

Selection of Spatial Variables
According to transition matrix analysis, we observed that the significant growth in built-up land was mainly due to cropland fragmentation. All these transitions were based on physical and socioeconomic driving factors. Researchers mostly use these spatial factors to investigate LULC dynamics. Table 14 shows the prospective Cramer's V value of each spatial variable. The Cramer's V value suggests that the variables are ideal for transition potential modeling, as their values are more significant. According to the values, the selection of physical and socioeconomic explanatory variables is more effectual, i.e., DEM: 0.35; slope: 0.33; population: 0.30; and GDP: 0.28. Figure 7 shows the spatial variables used in this study.

Selection of Spatial Variables
According to transition matrix analysis, we observed that the significant growth in built-up land was mainly due to cropland fragmentation. All these transitions were based on physical and socioeconomic driving factors. Researchers mostly use these spatial factors to investigate LULC dynamics. Table 14 shows the prospective Cramer's V value of each spatial variable. The Cramer's V value suggests that the variables are ideal for transition potential modeling, as their values are more significant. According to the values, the selection of physical and socioeconomic explanatory variables is more effectual, i.e., DEM: 0.35; slope: 0.33; population: 0.30; and GDP: 0.28. Figure 7 shows the spatial variables used in this study.

Transition Potential Modeling and Model Validation
The MOLUSCE plugin integrates some well-known algorithms for transition potential modeling, such as the ANN (multilayer perception), weights of evidence, multicriteria evaluation, logistic regression, and CA algorithm, for future simulation. The spatial variables for model calibration were chosen based on their relatively strong association with LULC, as measured by Cramer's coefficient.
We used the CA-ANN approach for transition potential modeling and prediction. We employed LULC data from 2000-2010 along with spatial variables to projected LULC for 2020 and obtained a validation kappa value of 0.76. After obtaining the projected LULC, we compared the actual LULC of 2020 with the projected data and obtained an overall accuracy of 96.25% and an overall kappa value of 0.94. Figure 8 and Table 15 show the actual and forecasted maps and statistics for 2020.

Transition Potential Modeling and Model Validation
The MOLUSCE plugin integrates some well-known algorithms for transition potential modeling, such as the ANN (multilayer perception), weights of evidence, multicriteria evaluation, logistic regression, and CA algorithm, for future simulation. The spatial variables for model calibration were chosen based on their relatively strong association with LULC, as measured by Cramer's coefficient.
We used the CA-ANN approach for transition potential modeling and prediction. We employed LULC data from 2000-2010 along with spatial variables to projected LULC for 2020 and obtained a validation kappa value of 0.76. After obtaining the projected LULC, we compared the actual LULC of 2020 with the projected data and obtained an overall accuracy of 96.25% and an overall kappa value of 0.94. Figure 8 and Table 15 show the actual and forecasted maps and statistics for 2020.

Prediction of LULC
After obtaining satisfactory results from model validation, we predicted the LULC for 2030, 2040 and 2050. We employed the temporal LULC data from 2010 and 2020, the spatial variables, and the transition probability matrix (Table 9) to predict the LULC for 2030 and obtained a kappa value of 0.95. Furthermore, the LULC of 2000 and 2020, including the explanatory variables and transition matrix (Table 10), were used for the prediction of 2040, and we obtained a kappa value of 0.78. Finally, we predicted the LULC for 2050 based on the projected data for 2030-2040 and the transition matrix (Table 11) and obtained a kappa value of 0.91. Figure 9 and Table 16 represent the map, area statistics, and overall validation kappa of the predicted LULC of 2030, 2040, and 2050.
Land 2021, 10, x FOR PEER REVIEW

Prediction of LULC
After obtaining satisfactory results from model validation, we predicted the for 2030, 2040 and 2050. We employed the temporal LULC data from 2010 and 20 spatial variables, and the transition probability matrix (Table 9) to predict the LU 2030 and obtained a kappa value of 0.95. Furthermore, the LULC of 2000 and 2020, ing the explanatory variables and transition matrix (Table 10), were used for the tion of 2040, and we obtained a kappa value of 0.78. Finally, we predicted the LU 2050 based on the projected data for 2030-2040 and the transition matrix (Table obtained a kappa value of 0.91. Figure 9 and Table 16

Intensity Analysis
After change detection, transition analysis and prediction, we analyzed the change in terms of intensity at three levels, i.e., interval level, category level, and transition level [70]. Table 17 and Figure 10 show the change intensity at the interval level in which we observed that the uniform intensity was 0.90. The annual changes during 1980-1990 and 2010-2020 were smaller than the uniform intensity, i.e., 0.53 and 0.40, indicating that the change intensity was slow compared with the uniform intensity. During the second and third intervals (1990-2000 and 2000-2010), the annual change intensity was greater than the uniform change intensity, which supported the fast change during these intervals.

Intensity Analysis
After change detection, transition analysis and prediction, we analyzed the change in terms of intensity at three levels, i.e., interval level, category level, and transition level [70]. Table 17 and Figure 10 show the change intensity at the interval level in which we observed that the uniform intensity was 0.90. The annual changes during 1980-1990 and 2010-2020 were smaller than the uniform intensity, i.e., 0.53 and 0.40, indicating that the change intensity was slow compared with the uniform intensity. During the second and third intervals (1990-2000 and 2000-2010), the annual change intensity was greater than the uniform change intensity, which supported the fast change during these intervals.  We analyzed each category's intensity in category-level analysis and compared the gain and loss intensity with the uniform intensity during each interval. Table 18 represents the gains and uniform intensity of each category during each interval. Except for the first interval, the built-up area had the largest size compared with uniform intensity and other categories, which indicated the fast gain in the built-up area as a result of the fragmentation of other categories. Table 19 shows the category-level loss intensity. During the first and second intervals, the loss intensity value of cropland was the highest (0.86 and 1.57, respectively), while during the third interval, cropland, grassland and water lost the most, as their lost intensities were larger than uniform intensity, with values of 1.91, 2.39, and 2.64, respectively. During the final interval, the cropland, grassland, and water intensities were greater than uniform intensity, which indicated a fast fragmentation in these categories and eventually the contribution to other categories, especially to built-up area, which was further described in the transition-level intensity section. We analyzed each category's intensity in category-level analysis and compared the gain and loss intensity with the uniform intensity during each interval. Table 18 represents the gains and uniform intensity of each category during each interval. Except for the first interval, the built-up area had the largest size compared with uniform intensity and other categories, which indicated the fast gain in the built-up area as a result of the fragmentation of other categories. Table 19 shows the category-level loss intensity. During the first and second intervals, the loss intensity value of cropland was the highest (0.86 and 1.57, respectively), while during the third interval, cropland, grassland and water lost the most, as their lost intensities were larger than uniform intensity, with values of 1.91, 2.39, and 2.64, respectively. During the final interval, the cropland, grassland, and water intensities were greater than uniform intensity, which indicated a fast fragmentation in these categories and eventually the contribution to other categories, especially to built-up area, which was further described in the transition-level intensity section.  Table 20 and Figure 11 show the transition of other categories to built-up area during the study period. Overall, cropland and grassland contributed the most considerable amounts to built-up area. During the first interval transition, the intensity of cropland was 0.35, and the uniform intensity was 0.13. Hence, the transition intensity was greater than the uniform intensity, indicating fast fragmentation of cropland to built-up area. During 1990-2000, cropland and grassland transition intensities were greater than those of other categories, i.e., 0.60 and 0.66, respectively, while the uniform intensity during this period was 0.45. During 2000-2010, the uniform intensity was 0.94, and the cropland, grassland, and water transition intensities were 1.16, 1.23, and 1.19, respectively. During the last interval, the uniform intensity was 0.30, and the transition intensities of cropland, grassland, and water were 0.40, 0.34, and 0.35, respectively. Therefore, the gain in built-up land from cropland means that cropland's transition intensity was greater than the uniform intensity during the study period.  Figure 11. Transition to Built-up Area. Table 21 and Figure 12 represent the size of transition intensity of other categories to forest against the uniform intensity, indicating that grassland contributed more than the other classes to forest during the first three intervals. In contrast, during the last interval, the built-up area donated more than the other categories. During the first three intervals, the grassland transition intensity was 0.41, 0.38, and 0.48, respectively, and the uniform intensity was 0.17, 0.22, and 0.31, respectively. During 2010-2020, the built-up land and   Figure 12 represent the size of transition intensity of other categories to forest against the uniform intensity, indicating that grassland contributed more than the other classes to forest during the first three intervals. In contrast, during the last interval, the built-up area donated more than the other categories. During the first three intervals, the grassland transition intensity was 0.41, 0.38, and 0.48, respectively, and the uniform intensity was 0.17, 0.22, and 0.31, respectively. During 2010-2020, the built-up land and grassland transition intensities were 0.03 and 0.02, respectively, while the uniform intensity was 0.02. Hence, forest gains would cause a transition in grassland as grassland donated more than the other categories.  Figure 11. Transition to Built-up Area. Table 21 and Figure 12 represent the size of transition intensity of other categories to forest against the uniform intensity, indicating that grassland contributed more than the other classes to forest during the first three intervals. In contrast, during the last interval, the built-up area donated more than the other categories. During the first three intervals, the grassland transition intensity was 0.41, 0.38, and 0.48, respectively, and the uniform intensity was 0.17, 0.22, and 0.31, respectively. During 2010-2020, the built-up land and grassland transition intensities were 0.03 and 0.02, respectively, while the uniform intensity was 0.02. Hence, forest gains would cause a transition in grassland as grassland donated more than the other categories.     Figure 13 show the transition intensity of other categories to cropland. We observed that the water transition intensity was more significant than the uniform intensity throughout the study period, i.e., 0.22, 0.37, 1.26 and 0.52, while the uniform intensity was 0.13, 0.22, 0.42 and 0.52, respectively.   Figure 13 show the transition intensity of other categories to cropland. We observed that the water transition intensity was more significant than the uniform intensity throughout the study period, i.e., 0.22, 0.37, 1.26 and 0.52, while the uniform intensity was 0.13, 0.22, 0.42 and 0.52, respectively.   Table 23 and Figure 14 show the transition to water from other categories. During the first three intervals, cropland's transition intensity was more significant than the uniform intensity, which indicated that the transition of cropland to water was faster than other categories, i.e., 0.32, 0.69 and 0.56, respectively, while the uniform intensity was 0.11, 0.25 and 0.23, respectively. During the last interval, as we discussed earlier, the water size decreased, so the other categories did not contribute to water.    Figure 14 show the transition to water from other categories. During the first three intervals, cropland's transition intensity was more significant than the uniform intensity, which indicated that the transition of cropland to water was faster than other categories, i.e., 0.32, 0.69 and 0.56, respectively, while the uniform intensity was 0.11, 0.25 and 0.23, respectively. During the last interval, as we discussed earlier, the water size decreased, so the other categories did not contribute to water.    Figure 15 represent the expected change intensity during 2030-2050. According to the results, during 2030-2040, the annual change intensity was greater than the uniform change intensity, i.e., 0.57 and 0.56, respectively, while during 2040-2050 the annual change intensity was smaller than the uniform intensity (0.54 and 0.56, respectively). The observed change during 2030-2040 was 5.72 and that during 2040-2050 was   Figure 15 represent the expected change intensity during 2030-2050. According to the results, during 2030-2040, the annual change intensity was greater than the uniform change intensity, i.e., 0.57 and 0.56, respectively, while during 2040-2050 the annual change intensity was smaller than the uniform intensity (0.54 and 0.56, respectively). The observed change during 2030-2040 was 5.72 and that during 2040-2050 was 5.39.    Figure 15 represent the expected change intensity during 2030-According to the results, during 2030-2040, the annual change intensity was greater the uniform change intensity, i.e., 0.57 and 0.56, respectively, while during 2040-205 annual change intensity was smaller than the uniform intensity (0.54 and 0.56, re tively). The observed change during 2030-2040 was 5.72 and that during 2040-2050 5.39.   Tables 25 and 26 show the gain and loss intensity at the category level. Accordi the projected results, our study area will experience a rapid expansion in built-up ar the future, i.e., the gain intensity of the built-up area will be greater than the uniform tensity during both intervals (1930-1940 and 1940-1950). The uniform intensity du  Tables 25 and 26 show the gain and loss intensity at the category level. According to the projected results, our study area will experience a rapid expansion in built-up area in the future, i.e., the gain intensity of the built-up area will be greater than the uniform intensity during both intervals (1930-1940 and 1940-1950). The uniform intensity during 2030-2040 and 2040-2050 will be 0.57 and 0.54, respectively, while the built-up area's gain intensity will be 1.85 and 2.12, respectively. Accordingly, cropland will still experience rapid fragmentation in the future. During 2030-2040, the loss intensity of cropland will be 1.92 against a uniform intensity of 0.57, while during 2040-2050, cropland will experience a loss intensity of 2.68 against a uniform intensity of 0.54.

Discussion
Unprecedented urban expansion has transformed the natural environment and landscape patterns worldwide, especially in the twenty-first century [89]. Physical and socioeconomic factors such as geography, population and economic growth are considered the most significant driving forces of urbanization [3]. However, socioeconomic development has a greater impact on urban expansion than population growth [90]. The extent and pace of urban expansion and fragmentation of landscape patterns have led to concerns about climatic changes, food security and scarcity of natural resources.
LULC changes are closely linked to geographical location and development policies. After the 'opening up' policy of China in the late 1970s, the economic reforms resulted in massive migration, immigration and urban expansion, especially in the southern part of the country. Based on the spatiotemporal LULC data and physical and socioeconomic driving factors, we analyzed the change from 1980 to 2020 and created a transition probability matrix for each interval using the MOLUSCE plugin within QGIS software. Furthermore, with the help of the CA-ANN multilayer perception approach within the MOLUSCE plugin, we projected the LULC for 2030, 2040, and 2050.
Our results indicate that physical and socioeconomic factors significantly affected landscape patterns during the study period. Generally, areas with lower elevations are associated with rapid LULC changes, as the geography of such areas is more suitable for anthropogenic activities. The highest changes occurred in the plain areas of the GBA, especially along the coast of the Pearl River, as the slope of this part is comparatively lower than that of other parts. The northern, eastern, and western portions, which consist of mountains, hills, and forest, do not experience rapid fragmentation.
The concept of GBA embodies broader goals to achieve market-oriented reforms and cooperation among the eleven cities in the socioeconomic sector, education, innovation, international business, and technological advancement. According to the "Outline Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area", by 2035, the GBA will become a globally competitive modern economic and industrial system; however, increasing the population and scarcity of land resources is one of the major challenges in the GBA [79]. Many studies have indicated that population growth and economic development are the driving forces of built-up expansion [91]. The increasing built-up area negatively impacts the natural environment, water quality and biodiversity [14]. According to our results, the GBA has experienced a dramatic conversion of LULC during the last four decades, especially in the rapid transition of cropland into built-up area. From 1980 to 2020, the built-up area increased from 4.75% to 14.75%, and cropland contributed 25.47%, grassland contributed 14.35%, forest contributed 9.93% and water bodies contributed 5.09%. Additionally, the future simulation results endorse that built-up land will continue to increase from 2030 to 2050 with percentages of 14.76% to 21.17%, and cropland will continue to decrease from 16.11% to 9.93%.
Furthermore, intensity analysis helps in analyzing and understanding the size and intensities of the transition at the interval, category, transition levels and the change patterns [70]. In this study, the intensity analysis reflected the rapid transition phenomenon in the LULC categories during past and future scenarios. During the first and fourth time intervals, the rate of annual change intensity was slower than that during the second and third time intervals, which indicated that the impacts of socioeconomic driving factors during the four decades were different. Furthermore, except for the first time interval, the gain intensity in the built-up area was more significant than the uniform intensity. The loss intensity of cropland was greater than the uniform intensity, which supports the result that the overall gain in the built-up area and the overall loss in cropland were greater than other categories. Similarly, the transition intensity of cropland to built-up area was greater than the uniform intensity. Hence, the gain in built-up area targets cropland more than other categories throughout the study period.
Ultimately, the drastic changes in LULC, especially urban sprawl and cropland fragmentation, could endanger natural resources, the environment, and food security. Therefore, the spatiotemporal and projected LULC simulation results will help policymakers analyze the change intensity of LULC and the impacts of socioeconomic factors and help with the promotion of environmental conservation and sustainable development strategies. Moreover, we considered only the physical and socioeconomic factors in LULC modeling and prediction. However, development policies and climatic factors can also be considered in future studies.

Conclusions
Based on four decades of data and physical and socioeconomic factors, we used the CA-ANN model within the MOLUSCE plugin in QGIS software to quantify the transition, patterns, and future simulation of LULC in the Greater Bay Area. Along with the CA-ANN approach, we used an intensity analysis technique to analyze the size and intensity of land change over time. The results show that the cropland and grassland in the GBA are facing enormous pressure and will face further pressure in the future. The spatial variables, including population, GDP, DEM, slope, and proximity factors, are the most critical driving factors, as they significantly affected the LULC change mechanism.
Due to the rapid urban growth and fragmentation of cropland, forest, and grassland, the GBA faces a shortage of agricultural land, environmental degradation and water quality depletion. These issues are creating challenges in maintaining regional development and environmental protection. In this study, we considered only physical and socioeconomic factors for modeling and prediction, but development policies, migration, immigration and climatic factors may also influence landscape patterns. Furthermore, agricultural policies and development strategies can be interconnected to promote sustainable urbanization. It is recommended that future studies consider the data and more variables to explore their impacts on landscape patterns.
Author Contributions: Z.A. designed the approaches, processed the data, performed the analysis, and wrote the manuscript. Y.Z. (Yuanjun Zhong), and Y.Z. (Yaolong Zhao) proposed the basic idea, designed the approaches, and supervised the study. G.Y. helped to modify the manuscript and provide useful suggestions. All authors have read and agreed to the published version of the manuscript.