Modeling Urban Expansion in Bangkok Metropolitan Region Using Demographic–Economic Data through Cellular Automata-Markov Chain and Multi-Layer Perceptron-Markov Chain Models

Abstract

Urban expansion is considered as one of the most important problems in several developing countries. Bangkok Metropolitan Region (BMR) is the urbanized and agglomerated area of Bangkok Metropolis (BM) and its vicinity, which confronts the expansion problem from the center of the city. Landsat images of 1988, 1993, 1998, 2003, 2008, and 2011 were used to detect the land use and land cover (LULC) changes. The demographic and economic data together with corresponding maps were used to determine the driving factors for land conversions. This study applied Cellular Automata-Markov Chain (CA-MC) and Multi-Layer Perceptron-Markov Chain (MLP-MC) to model LULC and urban expansions. The performance of the CA-MC and MLP-MC yielded more than 90% overall accuracy to predict the LULC, especially the MLP-MC method. Further, the annual population and economic growth rates were considered to produce the land demand for the LULC in 2014 and 2035 using the statistical extrapolation and system dynamics (SD). It was evident that the simulated map in 2014 resulting from the SD yielded the highest accuracy. Therefore, this study applied the SD method to generate the land demand for simulating LULC in 2035. The outcome showed that urban occupied the land around a half of the BMR.

1. Introduction

In the past three decades, the cities in developing countries have experienced a rapid increase in the rate of the population growth. According to the world urbanization prospects 2014 revision, the trend in urbanization recently has been focusing on the population transition between urban and rural [1]. Developing land for an urbanization and migration of people from rural areas is a global phenomenon. The statistical report endorses that about 53 percent of the population live in urban areas in 2015. About 77 percent of those people living in urban areas come from developed countries and 48 percent come from developing countries [2]. This change in urban population is continuously increasing year by year. By 2035, the world’s population living in urban areas will reach 5000 million, principally in African and Asian countries [3]. This situation will result in land use, economic, and environmental deterioration problems if each country does not have a suitable approach for better urban planning.

Bangkok Metropolis (BM) is one of the largest primary cities in Southeast Asia [4]. BM has expanded rapidly in both size and population since the 1960s. Local organizations have tried to improve the management to keep up with the changes; however, it fails to stop the growth. The issue of the primacy of the city of Bangkok is even more pronounced. The main cause is that BM has defined roles in the National Social and Economic Development Plan (NESDP) virtually in every issue for the prosperity of the country in all forms [5]. It is the administrative center of the country, as well as the center of investment in real estate development, various business services, and industrial center. This results in the expansion of the built-up area in its vicinity and the characteristics of a metropolitan region have become more noticeable.

There is a policy to develop major and minor cities in each region of the NESDP to decentralize BM. In contrast, the NESDP continues to promote the development of the metropolitan region that combines BM and five surrounding provinces called Bangkok Metropolitan Region (BMR). The metropolitan role of the BM focuses on supporting the diversification of production and serves as a residential area to prevent the evacuation of the population. The urbanization pattern in the BMR was initially found along the Chao Phraya River and then expanded itself along the transportation network and centers of minor cities. Nowadays, the major transportation network plays as an important factor in extending the urban areas [6].

Presently, land use and land cover (LULC) analysis has been extensively used to monitor urban expansion forms in specific space and time. Many satellite images have been used by scientists to comprehend the spatial-temporal evolution of an urban area [7]. However, they could not account the need for prediction of future LULC without an incorporation of land use models [8,9]. These models are powerful techniques that can support the policy making and future land demand evaluation, and simulate the spatial pattern of the land use based on the driving factors [10]. The LULC simulation models in an urban study are developed from the theories of urban morphology and dynamic process of LULC to forecast the urban expansion in different patterns and scales [11,12]. Those models can be categorized into four types: (i) empirical and statistical models such as Markov Chain (MC), logistic regression, etc. [13,14]; (ii) dynamic models such as cellular automata (CA), agent-based model (AGB), genetic algorithm (GA), artificial neural network (ANN), system dynamic (SD), etc. [15,16,17,18,19,20,21]; (iii) integrated models such as conversion of land use and its effects at small regional extent (CLUE-S) and Dyna-CLUE [22]; and (iv) hybrid models, like Metronamica, land transformation model (LTM), land change modeling (LCM), SLUETH, etc. [23,24,25,26,27,28,29]. They are composed of the affirmation models, which consider the past changes in the calibration step, and a sigmoidal pattern of the goodness of fit in the change potential function [30,31,32].

To understand the urban expansion phenomena and LULC in BMR, this current study tries to address the complex phenomenon of the urban context by integrating physical, demographic, and economic data into the LULC model to calibrate and simulate the future land use change from the base year 2014. This study also employs a multi-criteria evaluation and multilayer perception to determine the driver factors from the CA-MC and LCM, respectively. Moreover, the projected population and economic growth rates are applied to compute the transition probability by extrapolation, MC, and SD to simulate the urban LULC in 2035.

2. Study Area

BMR is known as the biggest urban agglomeration in Thailand. It consists of BM and five bordering provinces, Nakhon Pathom, Pathum Thani, Nonthaburi, Samut Prakan and Samut Sakhon, which are located in the lower part of Chao Phaya Delta. The total area covers about 7762.00 km² (Figure 1) [33].

The registered population in the BMR is about 11 million in 2014. More than half of the population lives in BM with the population density of 5258 people per km². Samut Prakan is the second most populous province and it is also a Southeast Asian air transportation hub. The third province is Nonthaburi, which plays an important role as the government center and the major residential zone of the region. Pathum Thani is the center for education and research centers for science and technology as well as industrial and residential zones supporting the population in the Northern part. The other two provinces, Nakhon Pathom and Samut Sakhon, are the least populated. Both cities have a role in supporting the settlement of the population in the East as well as a space for agricultural activities and fishing [34]. Approximately 75 percent of the land in BMR is utilized for agriculture. The second largest LU is urban and built-up area, accounting for approximately 20 percent. The other LUs are water, miscellaneous, and forest, sequentially. In the case of economics, BMR has a major role in driving the Gross Domestic Product (GDP) of Thailand because it is the center of economic activities. Around 97 percent of Gross Regional Product (GRP) in BMR comes from non-agriculture sections, which BM serves as the primary source [35].

3. Data and Methodology

In order to simulate the future urban LULC of BMR, this study implemented four major processes. The first process was the extraction of LULC data from the satellite images for several years. The second process was that of data input preparation for urban LULC modeling. This was the most important step to prepare the data set for each simulation period. In the third process, an analysis of the urban LULC change was taken into account to identify the major classes that had changed in each specific period. The last process was to model the urban LULC change using IDRISI Selva software packet version 17 developed by Clark Labs, Clark University [36]. The details of these four processes are explained in the following sections, and the overall operation of the methodology is shown in Figure 2.

3.1. Preparation of Satellite Image and LULC Classification

Extracting LULC information from the satellite images is one of the most important processes in data preparation. Multiple dates of the Landsat data, acquired from the U.S. Geological Survey (USGS) via the USGS Global Visualization Viewer site, were utilized in this study. The five periods of satellite images were generated using the Landsat 5-TM images in November of 1998, 2003, 2008, and 2011 and by Landsat 8-OLI image for 2014 of path/row 129/50 and 129/51. After the mosaic process, a coupler of image in the same year applied a geometric correction and re-projection into a common UTM zone 47 North. The first 2014 Landsat image was geo-reference using topographic map from the Royal Thai Survey Department (RTSD). Then the other images were corrected through an image-to-image rectification method based on the corrected 2014 image. The total root mean square error (RMSE) of each image was assessed to confirm the geometric error, which was less than one pixel. Furthermore, the Landsat images were classified into five classes by the maximum likelihood classification. This method applies the Bayesian theory [37] to categorize the image data into the interested LULC map. The sets of training data were used to predict the probability of each pixel being classified in each class. Moreover, an accuracy assessment was applied to confirm the quality of the resulted LULC maps by the overall accuracy and Kappa coefficient. The reference data to be used in this process were very high resolution images of IKONOS and QuickBird, Landsat RGB, and ground truth data provided by Land Development Department (LDD) or gathered directly from the field survey. All classified results have an overall accuracy greater than 80% and also a kappa value is higher than 0.8, which shows perfective agreement.

3.2. Data Input for the Urban LULC Modeling

The affirmation satellite data pre-processing was an important step to extract the LULC in each year of study. To model the urban LULC, however, there were many input data to collect, such as demographic, economic data, and other geo-database of physical and environment related to urban features including topography, transportation network, stream, and administrative boundary. These data were created and stored in different formats and scaled according to the data owner. The details of the input data are shown in

Table 1. Input data for the urban land use and land cover (LULC) modeling.

Data	Year	Scale/Resolution	Sources
Demographic data: density of population and housing, growth rate of population and housing	1993, 1998, 2003, 2008, 2011, 2014	1:50,000	Department of Provincial Administration
Economic data: Gross regional and provincial product	1988, 1993, 1998, 2003, 2008, 2011, 2014	1:50,000	Office of the National Economic and Social Development Board
Geo-data base of physical and environment related to urban features	n.d.	1:25,000 to 1:50,000	Bangkok Metropolis Administration GIS center
Digital Elevation Model (DEM)	n.d.	30 m	Bangkok Metropolis Administration GIS center

.

Demographic data at the sub-district level and economic data at the provincial level were converted through rasterization, performed using ArcGIS 10.3 (ESRI, Redlands, CA, USA), to produce the raster input maps. Some distance maps were generated such as distance to river, transportation network, agriculture, build-up, forest, and miscellaneous to use in the data analysis step. Finally, all raster inputs were resized into 100-meter resolution.

3.3. Change Analysis

To investigate the changes in an urban area, change analysis is an important process to identify the areas that have changed between the two time periods of interest. Change analysis is based on the historical changes that can be assessed between time 1 and time 2 of the LULC maps. The changes identified are transitions from one LULC state to another within the specific time period. It is likely that with many LULC classes the potential combination of transitions can be complex. In this study, a cross-tabulation analysis from the differential times was applied to determine the LULC change over five periods: 1988–1993, 1993–1998, 1998–2003, 2003–2008, and 2008–2011.

3.4. Urban LULC Modeling

In this section, the study utilized the LULC model to quantify the urban and future LULC changes as shown in Figure 3. On the whole, the analysis workflow is the following: (i) Cellular Automata-Markov Chain; (ii) Multi-Layer Perceptron-Markov Chain; (iii) Model validation; and (iv) Model application. The details of each part are described in the following subsections.

3.4.1. Cellular Automata-Markov Chain (CA-MC)

It is a combination of a sophisticated method between cellular automata (CA) and Markov chain (MC). Firstly, CA is a discrete model from disciplines of physics and biology [38]. It adopts a natural adaptation from the John Conway’s Game of Life [36] to simulate and understand a complex behavior [39]. The CA consists of four important elements: (i) a grid cell, which is a spatial unit in a finite cellular space; (ii) the state, each cell also has an individual state in the system; (iii) a neighborhood, which is a set of surrounding cells in the grid system; and (iv) a transition rule, which applies a change to the cell state based on specific properties of the cell and its neighborhood. Secondly, MC is a stochastic process and widely useful technique in examining LULC changes [40]. MC uses the knowledge based probability to predict the entire state of a given time (t2) from an earlier time (t1) without historical consideration. A transition probability and area matrix are calculated to illustrate the likelihood of the changes in the pixel status and total area for the next time step [36]. Moreover, a set of condition probability images is also created to identify the probability that areas will be changed in the next generation. The CA-MC, built in IDRISI Selva version 17 [36], requires two different periods of LULC maps. The time interval for the future LULC prediction is based on a number of iterations determined by a user. The 5 × 5 mean filter of the contiguity constraint of the CA was used for the future land allocation.

3.4.2. Multi-Layer Perceptron-Markov Chain (MLP-MC)

This method uses the combination between multi-layer perceptron (MLP) and MC to predict the future LULC. The MLP-MC performed in Land Change Modeler (LCM) utilizes three stepwise processes: change analysis, transition potential modeling and change prediction [28,36]. The details of each process are examined in detail.

• Change Analysis

This step comes with two times of LULC maps like change analysis process. It uses these to identify the major transitions that can be modeled separately like sub-model in the next step.

• Transition Potential Determination

The transition potential analysis is also an important part to determine driver factors that cause a conversion between a pair of LULC classes [31,41,42]. Some affirmation data were passed through the Euclidean distance function to generate distance maps. Statistical data such as slope, population, and economic data are also used to calculate the density and growth rate maps. All of these factors are shown in

Table 2. Definition of driving factors of LULC conversion in each specific time.

Data	Detailed Description
Distance to agriculture	Distance to any agricultural activity
Distance to BTS and MRT	Distance to BTS and MRT line
Distance to built-up area	Distance to any man-made construction
Distance to forest	Distance to any forest
Distance to main road	Distance to any permanent road
Distance to miscellaneous	Distance to any open land and grassland
Distance to river	Distance to permanent flowing water and surface sources
Distance to urban center	Distance to city center
Growth rate of economics	Growth rate of GPP
Growth rate of housing	Annual growth rate of housing
Growth rate of population	Annual growth rate of population
Housing density	Housing unit density per square kilometer
Population density	Persons per square kilometer
Slope	Mean slope in degrees

.

In the driving factors identification, Cramer’s V, a popular nominal scale measurement, based on qui-square, was applied to identify the most associated factors between two categorical maps. Equation (1) shows the formula of Cramer’s V statistic where φ (phi) is the coefficient of contingency and t is the smaller of the number of rows minus one or the number of columns minus one [41].

$$V=\sqrt{\frac{{\phi }^2}{\mathrm{t}}}$$

(1)

According to the literature, if the Cramer’s V value is about 0.15 or higher, it means that the association is useful. Moreover, if its value is greater than 0.4, the association is good [38].

Once all the explanatory factors were tested in each land transition, the LCM adopted MLP to generate the transition potential maps. The MLP is a feed-forward artificial neural network in which all data are passed from input layers to output layers. Each neuron calculates a combined value that is the result of the multiplication of values in the nodes of the previous layer by the connection weights between the nodes. This combined value is then the input to an activation function of the node, typically a non-linear activation function such as a sigmoid function [43,44]. Further, a back propagation training algorithm is typically used for the supervised training of the neural network. The algorithm iteratively changes the weights of the neural connections between the nodes of the different layers in such a manner that the error between the nodes’ outputs and their desired outputs is minimized [45]. In each land transition sub-model of this study, the MLP generates a random set of sample cells that have and have not experienced the transition. Therefore, each sub-model, consisting of one transition and two examples of classes, including one transition class and one persistence class, was fed. The MLP constructs a network of neurons between two example classes and driving factors, together with a web of connections that consist of sets of weights. Then, the sample cells are divided into two groups. The first 50% of the sample cells are used for training and the second 50% for validation [31,45,46]. This process is to develop the relationship between explanatory factors and activation levels of the transition potential. Thus, the weight has been adjusted to minimize the error to improve accuracy [46,47]. When the accuracy rate is more than 80% or reaches the maximum iteration (1000 epochs), the transition potential map of each sub-model, which is an output of the MLP, yields the suitability of LULC classes for the future period.

• LULC Change Prediction

The LCM uses a change prediction process to determine the quality of the future LULC in specific time modeled by MC. It produces two types of results. The first is a hard prediction model which results from a projected LULC map through a land allocation function [36]. The second is a soft prediction which is a comprehensive vulnerability map based on the set of transition. Generally, the values of the soft prediction outputs are between 0 and 1 [28].

3.4.3. Model Validation

The validation is a process to assess the quality of the simulated LULC results by comparing with the reference maps from the classified images in the years of 1998, 2003, 2008, and 2011. This study used the cross tabulation based accuracy assessment to calculate an overall accuracy of each predicted map [48]. Moreover, the Kappa analysis was also used to examine the quality of results [49].

3.4.4. Model Application

Once all predicted LULC maps were assessed for acceptable accuracy, the model was used to predict the future LULC in 2014 and 2035. Although the demographic and economic data were introduced as inputs of the model already, much literature [50,51,52,53] also suggests that it will improve the good performance of the model to incorporate the population and economic growth rate for future predictions. This study considers that reason to integrate both important data to transition rate calculation through the statistical extrapolation and SD.

Statistical extrapolation is a method to estimate the land demand (transition probability) based on the future population. The probabilities of land changes were predicted from an existing LULC and population change. The predicted values were generated by simple linear regression [31]. In this study, the change in population throughout the specific time period was compared with the change in each LULC category to develop a ratio of the future change using ordinary forecast and trend functions.

SD is a non-spatial simulation technique developed in the late 1950s by Jay W. Forrester of the Massachusetts Institute of Technology. SD was originally applied in the engineering field called “Industrial Dynamics” and also applied in the urban planning field, namely “Urban Dynamics”. It has attracted adequate attention as a tool for system characteristics investigation rather than predictions based on numerical simulations [54]. To calculate the transition probability of land demand using SD, this study applied the existing and projected population proposed by DOPA and NESDB in 2014 and 2035, respectively, to compute annual growth rates. The economic growth rates, calculated from the GPP, were also used in the model. The process in SD for land demand calculation can be divided into two parts: the first is driving factors module and the second is a land use module. The driving factors module handled the impact of population and economic growth while the land use module paid attention to the interaction and conversions of land use classes that were driven by population and economic data. For example, the demand for the urban class would be increased if the population and GPP growth rates increased. On the other hand, the agriculture and forest classes would decrease because of those factors. This study calculated the land demands using STELLA software from isee system. It was used to design the stocks (land use classes) and flows (driving forces) and to automatically generate relevant equations according to the feedback loop of the land demand in the SD environment. Once all were completed, the model estimated the land demands for each LULC class of the interested year, which were then converted into the transition probability matric.

4. Results and Discussion

In order to understand the LULC changes and dynamics of urban expansion in the BMR, this study divides the results into four sections: (i) LULC classification and accuracy assessment from the satellites images; (ii) Analysis of LULC change in the study periods; (iii) LULC prediction and validation; and (iv) LULC model application.

4.1. LULC Classification and Accuracy Assessment

The LULC maps for the six years of the Landsat images, 1988, 1993, 1998, 2003, 2008, and 2011, were analyzed through the maximum likelihood algorithm and improved the map qualities using auxiliary datasets (Figure 4). The classification results from the six LULC maps indicated that the agricultural class was the predominate land in the BMR throughout the study period. However, its share decreased dramatically. In contrast, urban land seems to be increased continuously due to the land demand for residential and industrial purposes. There was a slight fluctuation in the amount of miscellaneous class. Water and forest areas were relatively stable.

This study employed the advantages of total accuracy value and the kappa value for evaluating the classified LULC maps, as listed in

Table 3. Accuracy Assessment of the LULC maps for the years 1988, 1993, 1998, 2003, 2008, and 2011.

Years	Overall Accuracy	Kappa
1988	85.19%	0.80
1993	86.25%	0.82
1998	90.12%	0.87
2003	86.95%	0.83
2008	87.65%	0.83
2011	89.11%	0.86

. The overall accuracy values of the LULC maps were observed to be 85.19, 86.25, 90.12, 86.95, 87.65, and 89.1 percent for the six time periods, respectively. Furthermore, the kappa values of the maps were 0.80, 0.82, 0.87, 0.83, 0.83, and 0.86, respectively. All values of the overall accuracy and the kappa coefficient were greater than 85 percent and 0.80, which indicated a strong perfect agreement between the classified maps and the ground truth data. Therefore, the accuracy of the LULC maps was satisfactorily accurate for investigating the future urban LULC simulations.

4.2. Change Analysis

To quantify the changing in an urban area, the LULC is the most important. A cross tabulation process was applied to identify the major changes between two LULC maps of the specified time periods. The LULC change maps together with a number of changes in each class from the five periods are displayed in Figure 5.

It can illustrate evidently that, firstly, the conversion between agriculture and urban was the most distinctive change of the study period. The amount of change increased gradually from approximately 300 km² in 1988–1993 to about 850 km² in 2008–2011. The second was the transition from agriculture to miscellaneous, which could be found from the analysis. The percentage of its change was less than 10 percent from the first four periods. On the other hand, in the last period (2008–2011), a significant increase in this transition was observed. The third was the conversion from miscellaneous to agriculture, which could be detected from the change analysis with a lower rate of change when compared to the second. Moreover, the conversion between miscellaneous land and urban was another influence that could be seen in the LULC change. Nevertheless, the change analysis indicated other transitions in each LULC class, which also had effects on the urban growth and LULC changes in the BMR.

4.3. Urban LULC Prediction

4.3.1. Identification of Driving Factors

After detecting the changes in particular time periods, the identification of the driving factors was investigated regarding the changes in the specific time duration. To determine the association between the changes and the factors for the LULC modeling and simulation, Cramer’s V values were calculated. The factors that had a Cramer’s V value higher than 0.15 were selected to generate the transition potential map of land transition in each sub-model based on a root mean square error (RMSE) value less than 0.1 within 10,000 iterations of the MLP neural network parameters.

Table 4. Sub-modules and Cramer’s V value for simulation in the years of study.

Year	Factor	Overall V	A	W	U	F	M
1998	Distance to agriculture	0.4533	0.6302	0.9938	0.2662	0.9176	0.1941
	Distance to built-up area	0.3843	0.8002	0.7752	0.2400	0.8207	0.0654
	Distance to center of city	0.1725	0.2229	0.3752	0.1268	0.3429	0.0918
	Distance to forest	0.1863	0.3178	0.2590	0.1170	0.2409	0.3127
	Distance to miscellaneous	0.1514	0.1373	0.1381	0.0618	0.1085	0.0806
	Distance to road	0.1349	0.1758	0.2737	0.0717	0.2902	0.0338
1998	Distance to water	0.2699	0.3649	0.3191	0.5405	0.2517	0.0658
	Economic growth rate	0.1582	0.5003	0.3029	0.0994	0.3112	0.1283
	Elevation	0.2913	0.1914	0.1896	0.6409	0.0989	0.0349
	Slope	0.0652	0.0151	0.0532	0.1273	0.0303	0.0216
2003	Distance to agriculture	0.4498	0.7677	0.9949	0.2285	0.9328	0.1356
	Distance to built-up area	0.4526	0.4701	0.9383	0.1378	1.0000	0.0389
	Distance to center of city	0.1698	0.2121	0.3648	0.1209	0.3362	0.0958
	Distance to forest	0.1790	0.9365	0.2334	0.1086	0.2193	0.3051
	Distance to miscellaneous	0.4520	0.3931	0.1998	0.0646	0.1342	0.0283
	Distance to road	0.1462	0.2646	0.3036	0.0633	0.3181	0.0284
2003	Distance to water	0.2481	0.7880	0.2448	0.5162	0.1779	0.0759
	Economic growth rate	0.1547	0.8363	0.3090	0.0788	0.3102	0.1182
	Elevation	0.2722	0.3918	0.1625	0.6009	0.0945	0.0247
	Housing density	0.2401	0.6925	0.5083	0.1114	0.5019	0.0555
	Housing growth rate	0.1643	0.9972	0.3150	0.0954	0.3137	0.1450
	Population density	0.2337	0.9773	0.4975	0.0964	0.4939	0.0362
	Population growth rate	0.1997	0.6992	0.4047	0.0927	0.3997	0.1077
	Slope	0.0776	0.1496	0.0582	0.1694	0.0297	0.0180
2008	Distance to agriculture	0.3668	0.5380	0.7965	0.1635	0.7814	0.0837
	Distance to BTS and MRT	0.2466	0.3331	0.5356	0.0902	0.5235	0.0696
	Distance to built-up area	0.3918	0.7218	0.8035	0.2060	0.8427	0.0534
	Distance to center of city	0.1734	0.1824	0.3712	0.1161	0.3462	0.0753
	Distance to forest	0.2033	0.5839	0.3141	0.2883	0.2908	0.2760
	Distance to miscellaneous	0.1347	0.6038	0.2719	0.0159	0.2650	0.0557
	Distance to road	0.1651	0.0806	0.3463	0.0751	0.3564	0.0322
	Distance to water	0.2943	0.6357	0.2366	0.6190	0.1979	0.0859
	Economic growth rate	0.1802	0.3877	0.3675	0.1105	0.3495	0.1177
	Elevation	0.2585	0.3879	0.1397	0.5676	0.1075	0.0219
	Housing density	0.2430	0.6503	0.5161	0.1292	0.5009	0.0691
	Housing growth rate	0.2061	0.3140	0.4045	0.1720	0.3869	0.1101
	Population density	0.2384	0.4323	0.5101	0.1147	0.4940	0.0687
	Population growth rate	0.2079	0.3740	0.4047	0.1504	0.3915	0.1454
	Slope	0.0661	0.1865	0.0284	0.1447	0.0260	0.0137
2011	Distance to agriculture	0.4496	0.5449	0.9966	0.1789	0.9238	0.1283
	Distance to BTS and MRT	0.2466	0.5991	0.5356	0.0902	0.5253	0.0696
	Distance to built-up area	0.4723	0.7801	0.9306	0.3010	1.0000	0.0656
	Distance to center of city	0.1734	0.7794	0.3712	0.1161	0.3462	0.0753
	Distance to forest	0.2027	0.4043	0.2934	0.1734	0.2912	0.2918
	Distance to miscellaneous	0.2338	0.4055	0.3398	0.0567	0.2914	0.0384
	Distance to road	0.1651	0.0109	0.3463	0.0751	0.3564	0.0322
2011	Distance to water	0.3646	0.5748	0.2449	0.7866	0.1878	0.0899
	Economic growth rate	0.1801	0.0298	0.3674	0.1104	0.3495	0.1172
	Elevation	0.2585	0.6823	0.1397	0.5676	0.1075	0.0219
	Housing density	0.2404	0.5694	0.5053	0.1330	0.4951	0.0776
	Housing growth rate	0.2322	0.2842	0.4554	0.1634	0.4384	0.1609
	Population density	0.2370	0.4840	0.5058	0.1209	0.4909	0.0648
	Population growth rate	0.2136	0.7338	0.4244	0.1551	0.4086	0.1311
	Slope	0.0661	0.1697	0.0284	0.1447	0.0260	0.0137
2014	Distance to agriculture	0.4491	0.6516	0.9981	0.1619	0.9456	0.1127
	Distance to BTS and MRT	0.2272	0.3075	0.4932	0.0771	0.4896	0.0627
	Distance to built-up area	0.4559	0.7176	0.9484	0.2107	0.4552	0.0821
	Distance to center of city	0.1681	0.4045	0.3639	0.0916	0.3325	0.0939
	Distance to forest	0.2069	0.6778	0.2977	0.1547	0.3009	0.3128
	Distance to miscellaneous	0.2781	0.4165	0.0854	0.2331	0.1068	0.0794
	Distance to road	0.1545	0.2389	0.3302	0.5810	0.3345	0.0318
	Distance to water	0.2663	0.4229	0.1508	0.0660	0.0980	0.0786
	Economic growth rate	0.1700	0.1255	0.3366	0.0718	0.3415	0.1443
	Elevation	0.2764	0.0618	0.1358	0.6100	0.1003	0.0275
	Housing density	0.2251	0.2134	0.4769	0.0991	0.4732	0.1076
	Housing growth rate	0.1847	0.1065	0.3495	0.1267	0.3514	0.1670
	Population density	0.2218	0.2052	0.4712	0.0880	0.4671	0.1105
	Population growth rate	0.2075	0.3831	0.4078	0.1104	0.3978	0.1905
	Slope	0.0814	0.0361	0.0595	0.1719	0.0431	0.0182

Note: A = agriculture; W = water; U = Urban; F = forest; M = miscellaneous.

shows the influence factors that were used to generate a transition potential map for each specific period for the simulation.

4.3.2. Transition Probability Analysis

To simulate the future LULC and urban expansion, this study considered a sophisticated method to estimate the rate of land transition in each specific time using the conventional scheme of the MC and application of the projected population to calculate the transition probabilities from statistical extrapolation and SD approach. The details of this analysis are discussed in two separate sections.

• Transition probability calculation from the MC

The MC was used to quantify the change between a pair of LULC maps. It created a transition probability matrix for each LULC category to be converted to any other categories (

Table 5. Transition probabilities of LULC change.

Class		Agriculture	Water	Urban	Forest	Miscellaneous
1998	Agriculture	0.9486	0.0025	0.0418	0.0032	0.0039
	Water	0.1121	0.8490	0.0358	0.0005	0.0026
	Urban	0.0000	0.0041	0.9880	0.0000	0.0079
	Forest	0.3487	0.0073	0.0160	0.6236	0.0044
	Miscellaneous	0.6223	0.0180	0.3080	0.0000	0.0516
2003	Agriculture	0.9465	0.0022	0.0451	0.0013	0.005
	Water	0.1745	0.7905	0.0269	0.0034	0.0047
	Urban	0.0000	0.0032	0.9905	0.0002	0.0061
	Forest	0.6245	0.0016	0.0186	0.3524	0.0029
	Miscellaneous	0.7160	0.0272	0.1657	0.0067	0.0843
2008	Agriculture	0.9112	0.0022	0.0682	0.0042	0.0141
	Water	0.1277	0.8167	0.0446	0.0051	0.0059
	Urban	0.0000	0.0030	0.9795	0.0011	0.0164
	Forest	0.294	0.0010	0.0805	0.6245	0.0000
	Miscellaneous	0.5427	0.0069	0.2602	0.0024	0.1878
2011	Agriculture	0.9357	0.0025	0.0458	0.0015	0.0144
	Water	0.1349	0.8521	0.0126	0.0005	0.0000
	Urban	0.0000	0.0021	0.9922	0.0002	0.0055
	Forest	0.2310	0.0155	0.0463	0.6426	0.0646
	Miscellaneous	0.3111	0.0051	0.4806	0.0000	0.2033
2014	Agriculture	0.8221	0.0097	0.1586	0.0051	0.0045
	Water	0.1542	0.7630	0.0650	0.0158	0.0021
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.3185	0.0159	0.0494	0.6148	0.0015
	Miscellaneous	0.5736	0.0077	0.3688	0.0142	0.0357
2035	Agriculture	0.1780	0.0140	0.8023	0.0046	0.0010
	Water	0.2400	0.0849	0.6624	0.0111	0.0015
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.2504	0.0213	0.7122	0.0147	0.0015
	Miscellaneous	0.1327	0.0105	0.8525	0.0036	0.0008

). The contingency matrix, which displays the relative frequencies of change in a certain time step, was utilized to calculate the transition probability. The resulted matrices were then used with the transition potential map in the next step of the change prediction.

• Integration of the future population and economic demand for the transition probability of urban expansion in 2014 and 2035

To quantify the future LULC in 2014, this study considered the existing population data collected by DOPA, Ministry of Interior. Based on the population growth calculation between 2011 and 2014, the annual growth rate of population was observed as 0.01 per year. As reported by NESDB, the future population in the BMR would reach up to 16,632,220 in 2035. The increase in the rate of the population when compared with the base year 2014 was 0.15 per year. Moreover, the economic growth rate obtained from the GPP change in the same period would approach a value of 0.12 per year. Here, both data were used to derive the transition probability using simple extrapolation and the non-spatial technique of SD.

Table 6. Transition probability rate in each LULC class calculated by extrapolation and system dynamic model for 2014.

Classes		Agriculture	Water	Urban	Forest	Miscellaneous
Extrapolation	Agriculture	0.8426	0.0099	0.1393	0.0030	0.0052
	Water	0.1541	0.7854	0.0481	0.0125	0.0000
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.2861	0.0234	0.0326	0.6255	0.0324
	Miscellaneous	0.0809	0.0030	0.7959	0.0056	0.1146
SD	Agriculture	0.8995	0.0055	0.0896	0.0029	0.0026
	Water	0.2016	0.7976	0.0007	0.0002	0.0000
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.4011	0.0001	0.5986	0.0002	0.0000
	Miscellaneous	0.0079	0.0001	0.0051	0.0002	0.9868

and

Table 7. Transition probability rate in each LULC class calculated by extrapolation and system dynamic model for 2035.

Classes		Agriculture	Water	Urban	Forest	Miscellaneous
Extrapolation	Agriculture	0.5941	0.0000	0.3860	0.0022	0.0177
	Water	0.1717	0.7269	0.0728	0.0286	0.0000
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.1458	0.0862	0.6602	0.1000	0.0078
	Miscellaneous	0.7859	0.0136	0.0000	0.0009	0.1996
SD	Agriculture	0.5365	0.0253	0.4130	0.0133	0.0118
	Water	0.1535	0.7644	0.0647	0.0158	0.0017
	Urban	0.0000	0.0000	1.0000	0.0000	0.0000
	Forest	0.2276	0.0012	0.7668	0.0042	0.0002
	Miscellaneous	0.1366	0.0018	0.7703	0.0034	0.0879

show the estimated transition values in 2014 and 2035, respectively. Then, the LULC transitions were applied in the CA-MC and MLP-MC environments to generate the transition potential maps in each sub-model for the simulation.

4.3.3. LULC Prediction and Validation

• LULC Prediction from CA-MC and MLP-MC

This study considered two different schemes for simulation based on the conventional model to explore the urban LULC and urban expansion in the BMR. The simulated LULC maps in the four time periods 1998, 2003, 2008, and 2011 were generated, as shown in Figure 6 and

Table 8. Accuracy assessment of simulated LULC maps from the MC in CA and MLP environments.

Types	Years
	1998		2003		2008		2011
	CA-MC	MLP-MC	CA-MC	MLP-MC	CA-MC	MLP-MC	CA-MC	MLP-MC
Overall accuracy	92.3728	92.1738	88.4218	90.5811	89.4320	89.4969	85.5354	85.5690
K standard	0.9097	0.9108	0.8748	0.8979	0.8913	0.8924	0.8570	0.8571
Kno	0.9374	0.9383	0.9088	0.9257	0.9168	0.9171	0.8859	0.8863
KlocationStrata	0.9148	0.9146	0.9054	0.9298	0.9145	0.9072	0.9359	0.9356

. The validation procedure was investigated to confirm the reliability and efficiency of the CA-MC and MLP-MC for the LULC prediction. The overall accuracy and several kinds of Kappa coefficients were applied to validate the simulated maps with the LULC maps generated by the maximum likelihood classification. After the superimposing, all simulated maps yielded overall accuracies greater than 80 percent. Furthermore, Kappa coefficients such as Kstandard, Kno, and KlocationStrata were also calculated to assure the accurate result of the simulation, all of which had values higher than 0.8 (Table 8). However, the accuracy assessment from the simulated LULC maps proved that the MLP-MC method contributed higher accuracy values than CA-MC. Since the MLP-MC deliberated the influence in physical, demographic and economic factors on the LULC conversion, the simulated maps from this method showed a higher accuracy when compared with the CA-MC results (

Table 9. The simulated results of LULC maps from the MC transition probability in CA and MLP.

LULC Classes	Years
	1998		2003		2008		2011
	CA-MC	MLP-MC	CA-MC	MLP-MC	CA-MC	MLP-MC	CA-MC	MLP-MC
Agriculture	6359.45	6339.11	6075.25	6067.47	5399.02	5478.06	5027.78	5101.09
Water	74.96	81.25	77.86	76.58	79.76	82.47	93.99	97.51
Urban	1184.09	1192.15	1468.49	1489.37	2031.3	2055.99	2429.72	2424.47
Forest	41.03	52.42	36.83	15.77	52.53	53.54	29.2	33.65
Miscellaneous	35.36	28.32	36.46	44.06	131.87	23.19	114.66	36.53

and Figure 6). CA-MC does not, in contrast, consider the change spatially. This method is not always appropriate to allocate future LULC in suitable locations [14,55].

4.3.4. LULC Model Application

• Comparative study on the LULC prediction from the conventional MC and demographic–economics data to estimate the transition probability

According to the previous section, the accuracy assessment confirms that the MLP-MC method yielded a strong level of perfect agreement between the simulated maps and reference data. The model application in this study adopted the transition rates calculated by MC in CA-MC and MLP-MC methods to simulate the LULC in 2014. The transition probabilities calculated from the population and economics via the extrapolation and SD were also utilized in the MLP as a platform to simulate in 2014. The LULC map in 2011 was used as a base map for the future simulation. The results from all schemes are illustrated in Figure 7. All values of the overall accuracy and Kappa coefficients illustrated and certify the performance of all schemes that were stratified to simulate the future LULC (

Table 10. Accuracy assessment of simulated LULC maps in 2014 from four schemes.

Types	CA-MC	MLP-MC	Extrapolation	SD
Overall accuracy	90.4473	90.8636	93.5340	94.6711
Kstandard	0.8908	0.8949	0.9212	0.9324
Kno	0.9113	0.9148	0.9364	0.9451
KlocationStrata	0.9706	0.9634	0.9588	0.9754

). The transition rate calculated using the SD displayed the highest accuracy when compared with others. The SD scheme projected the urban area slightly lower than the existing LULC from the Landsat image; however, it showed a strong performance in the simulation (

Table 11. The actual and simulated area of LULC maps in 2014 from Landsat 8 and the transition probability calculated by MC, the extrapolation, and SD.

LULC Classes	Landsat 8	CA-MC	MLP-MC	Extrapolation	SD
Agriculture	4368.03	3715.57	3715.6	3814.48	4265.74
Water	123.56	138.47	147.58	148.48	123.01
Urban	3179.87	3761.32	3750.99	3677.12	3254.88
Forest	8.76	56.44	57.38	49.66	43.38
Miscellaneous	13.03	21.45	21.70	3.51	6.24

). From the extrapolation scheme, its result yielded an overestimation in the urban class and underestimation in agriculture class similar to the results from the conventional MC scheme when performed in CA and MLP environments.

• Simulation of the future urban LULC and urban expansion in 2035.

Because of the strong performance of the demographic–economics data, the transition probabilities from the statistic extrapolation and SD schemes were applied to simulate the LULC and urban area. Both transition probabilities were calculated from the estimated population and annual growth rate of 2035. The simulation, in this part, was based on the transition potential map created by the MLP of the LCM model. Figure 8 illustrates the simulated results from both methods. It can visualize evidently that, by using the SD scheme, the predicted urban area is 3750.77 km² or 48.75 percent of the whole area of the BMR. Agricultural land is almost equal to the urban class. The lands for water, forest, and miscellaneous classes occupy the smallest amounts of land, sequentially (

Table 12. The simulated area of LULC maps in 2035 from the transition probability calculated by the extrapolation and SD.

LULC Classes	Extrapolation		SD
	km	%	km	%
Agriculture	2478.6	32.22	3735.72	48.56
Water	217.19	2.82	147.67	1.92
Urban	4896.27	63.64	3750.77	48.75
Forest	98.55	1.28	57.45	0.75
Miscellaneous	2.64	0.03	1.64	0.02

). In terms of the simulation from the extrapolation, all of the LULC have a similar trend with the SD scheme. However, the proportion of the LULC is different, especially in the urban class, which occupies around 64 percent of land.

To analyze the urban expansion of the BMR in 2035, this study used the simulated map from the SD scheme to examine the urban extent (Figure 9). The conversion between agriculture and urban areas plays as a primary role in the LULC change. The large urban agglomeration is located along the Chao Phraya River where the urban and built-up areas have increased dramatically since 1988. The urban area expanded itself over administrative boundaries, particularly in the connected zone between BM and Pathum Thani, Nonthaburi, and Samut Prakan. Those provinces play an important role as residential areas supporting the BM population and also as lands for industrial purposes based on the fundamental function of the city. Moreover, the urban area expands in a linear form along the main roads that connect neighboring provinces to the BM. On the other hand, the center of the BMR remains as most population dene zone where all built up areas will be constructed in response to land demands for business and commercial enterprises. Since the construction of the Bangkok Mass Transit System (BTS) and Mass Rapid Transit (MRT), many condominiums with medium and high classes were built along those lines to respond to the demand for urban lifestyles with reduced commutes to work. Furthermore, there is a different pattern of urban expansion that occurred in Nakhon Pathom. Its urban area has expanded from the center of the province and will continuously extend until 2035. There is also a scattered random pattern in the built-up area, in which farmers occupy the land near agricultural areas. This situation can be found predominantly in Nakhon Pathom. If the spatial planners do not consider this, it will create urban sprawl, which is already a crucial problem in that area.

5. Conclusions

Recently, Thailand faces with the problem of urban expansion due to the national plan policy that was established to confirm BM as the center for all aspects of the country. Moreover, industrialization and innovation in the globalized world have turned BM into a modern city where everyone would like to live. However, BM does not grow by itself, but rather brings about spatial planning problems in the vicinity, which collectively can be called the BMR. This study used scientific methodology to explore the urban expansion and LULC changes in a specific time period based on the innovation of satellite imagery. Landsat images were utilized to generate LULC maps using the sophisticated maximum likelihood classification in different times. Then, the change detection technique was investigated to develop a pair of LULC maps in order to identify the land transformations. From this study, it can noticeably evident that the major land conversion in this area is agriculture to urban. To explore the future LULC and urban expansion in the future, this study applied the integrated methods of CA-MC and MLP-MC to generate transition probabilities and transition rates for the simulations. The simulated results in the years of 1988, 2003, 2008, and 2011 were validated to identify the best method for the predictions. The accuracy of the simulations showed that MLP-MC yielded higher accuracy than CA-MC. Therefore, the simulation in the years of 2014 and 2035 used the MLP as a platform for the prediction.

To project the LULC to the year 2035, even though the physical and economic factors were considered as input layers for the engendering the transition potential in each conversion, the future population and economic growth were applied to control the transition rate of the land conversions. The results from the study evidently confirm that in 2014 the predicted LULC maps from the SD and statistics extrapolation shows reliable results with robust accuracy, particularly the projected result from the SD method. Therefore, the transition probability rates, proposed by the SD method, within the MLP environment were used to explore the urban expansion in 2035. The simulation illustrates that urban and built-up area are the main LULC types and occupy about half of the land in the BMR. Most of the agricultural land would be converted to urban or built-up area if proper measures are not taken. This study shows that integration of the LULC modeling approach, remote sensing data, and demographic and economic data can be professionally used to simulate the future LULC change and, hence, the urban expansion phenomena in spatial distribution, direction, and time can be monitored. Future LULC is crucial information for planners and organizations to allocate important infrastructures and manage sustainable land use in response to the diverse requirements of the people in this region.