www.mdpi.com/journal/ijgi/ Modeling Urban Land Cover Growth Dynamics Using Multi-Temporal Satellite Images: A Case Study of Dhaka,

The primary objective of this research is to predict and analyze the future urban growth of Dhaka City using the Landsat satellite images of 1989, 1999 and 2009. Dhaka City Corporation (DCC) and its surrounding impact areas have been selected as the study area. At the beginning, a fisher supervised classification method has been applied to prepare the base maps with five land cover classes. In the next stage, three different models have been implemented to simulate the land cover map of Dhaka city of 2009. These have been named as “Stochastic Markov (St_Markov)” Model, “Cellular Automata Markov (CA_Markov)” Model and “Multi Layer Perceptron Markov (MLP_Markov)” Model. Then the best-fitted model has been selected by implementing a method to compare land cover categories in three maps: a reference map of time 1, a reference map of time 2 and a simulation map of time 2. This is how the “Multi Layer Perceptron Markov (MLP_Markov)” Model has been qualified as the most appropriate model for this research. Later, using the MLP_Markov model, the land cover map of 2019 has been predicted. The MLP_Markov model extrapolates that built-up area increases from 46% to 58% of the total study area during 2009–2019.


Background of the Research
Like many other cities in the world, Dhaka, the Capital of Bangladesh is also the outcome of spontaneous rapid growth.As the growth of population in Dhaka is taking place at an exceptionally rapid rate, it has become one of the most populous Mega Cities in the world.
Dhaka City has undergone radical changes in its physical form, not only in its vast territorial expansion, but also through internal physical transformations over the last decades.These have created entirely new kinds of urban fabric.In the process of urbanization, the physical characteristics of Dhaka City are gradually changing as plots and open spaces have been transformed into building areas, open squares into car parks, low land and water bodies into reclaimed built-up lands, etc.
Dhaka is now attracting a huge amount of rural-urban migrants from all over the country due to well-paid job opportunities, better educational, health and other daily life facilities.This kind of increasing and over population pressure is putting adverse impacts on Dhaka City like unplanned urbanization, extensive urban poverty, water logging, growth of urban slums and squatters, traffic jams, environmental pollution and other socio-economic problems [1].
If this situation continues then Dhaka would soon become an urban slum with the least livable situation for the city dwellers.Therefore, the primary objective of this paper is to forecast the future urban land cover changes of the selected study area within the greater Dhaka City.

Existing Works
With the advancement of technology, reduction in data cost, availability of historic spatio-temporal data and high resolution satellite images, Remote Sensing (RS) and Geographic Information System (GIS) techniques are now very useful for conducting researches like land cover change detection analysis and predicting the future scenario [2].
Many researchers have conducted number of researches to detect the land use/land cover change pattern over time and predict the future growth of urban areas.They have introduced and applied different techniques to achieve the research objectives.Among them, Griffiths et al. have approached to map the urban growth of Dhaka megacity region (1990 to 2006) using multi-sensoral data.They have used a Support Vector Machine (SVM) classifier and post-classification comparison to reveal Spatio-temporal patterns of urban land-use and land-cover changes [3].
Dewan and Yamaguchi have tried to evaluate land cover changes and urban expansion in greater Dhaka, between 1975 and 2003 using satellite images and socio-economic data.A supervised classification algorithm and the post-classification change detection technique in GIS have been implemented by them.They have found the accuracy of the Landsat-derived land cover maps ranged from 85% to 90% [4].
Emch and Peterson have quantified mangrove forest cover change in the Sundarbans of south-west Bangladesh from 1989 to 2000 using Landsat Thematic Mapper (TM) satellite imagery.They have used three image processing techniques: Normalized Differential Vegetation Index (NDVI), maximum likelihood classification and sub-pixel classification [5].
Kashem has implemented SLEUTH urban growth model to simulate the historical growth pattern of Dhaka Metropolitan Area.SLEUTH model incorporates Slope, Landuse, Exclusion layer (where growth cannot occur), Urban, Transportation and Hill-shade data layers.SLEUTH uses a modified Cellular Automata (CA) to model the spread of urbanization [6].Lahti has predicted the urban growth of Sydney, Australia till 2106.The modeling has been performed using the CA model Metronamica, developed by the Research Institute for Knowledge Systems (RIKS) in the Netherlands [7].Cabral and Zamyatin have implemented three land change models to forecast the urban dynamics in Sintra-Cascais municipalities of Portugal, for the year 2025.The models are CA Markov chain model (CA_Markov), CA_Advanced and Geomod.They have used image segmentation and texturing procedures to classify the Landsat images of 1989, 1994 and 2001 [8].
Wang and Mountrakis have developed a GIS-based modeling framework titled Multi-Network Urbanization (MuNU) model, which integrates multiple neural networks, to predict the urban growth of Denver Metropolitan Area, CO, USA [9].Cheng and Masser have used a spatial data analysis method that comprises exploratory data analysis and spatial logistic regression technique, to seek and model major determinants of urban growth in the period1993-2000 of Wuhan City in PR China [10].

Study Area
Dhaka is located in central Bangladesh at 23°43′0″N, 90°24′0″E, on the eastern banks of the Buriganga River (Figure 1(a)).Dhaka City is located in Dhaka District that is surrounded by rivers.The proposed study area for this research is Dhaka City Corporation (DCC) and its surrounding impact areas (Figure 1(b)).This area has increased potential to face massive urban growth in near future based on the current trend of rapid urbanization.The study area covers the oldest organic core part of Dhaka city (old Dhaka), the planned areas and even the unplanned new generation organic areas that are called "Informal Settlements".This selected study area almost covers the biggest urban agglomeration and is the central part of Bangladesh in terms of social and economic aspects [11].

Remote Sensing Data
To prepare the base maps for analysis purpose and applying the different methods to achieve the research objectives, the Landsat satellite images (1989, 1999 and 2009) have been collected from the official website of US Geological Survey (USGS).Landsat Path 137 Row 44 covers the whole study area.Map Projection of the collected satellite images is Universal Transverse Mercator (UTM) within Zone 46 N-Datum World Geodetic System (WGS) 84 and the pixel size is 30 m. Five land cover types have been identified for this research (Table 1).

Land Cover Type Description
Built-up Area All residential, commercial and industrial areas, villages, settlements and transportation infrastructure.

Water Body
River, permanent open water, lakes, ponds, canals and reservoirs.

Vegetation
Trees, shrub lands and semi natural vegetation: deciduous, coniferous, and mixed forest, palms, orchard, herbs, climbers, gardens, inner-city recreational areas, parks and playgrounds, grassland and vegetable lands.

Low Land
Permanent and seasonal wetlands, low-lying areas, marshy land, rills and gully, swamps, mudflats, all cultivated areas including urban agriculture; crop fields and rice-paddies.

Fallow Land
Fallow land, earth and sand land in-fillings, construction sites, developed land, excavation sites, solid waste landfills, open space, bare and exposed soils.
The Landsat satellite images used for analysis are of different qualities and dates (Table 2).For the purpose of ground information, several base maps of Dhaka City (for the year of 1987, 1995 and 2001) have been collected from the Survey of Bangladesh (SoB).Google Earth was another option to get some ideas about the recent land cover pattern of Dhaka city.The detail land use map (2009) of Dhaka

Base Map Preparation
Supervised classification relies on the a priori knowledge of the study area [12].Therefore, for this research, a supervised classification method has been used.

Training Site Development
Training sites are the areas defined for each land cover type within the image.The chosen color composite is used for digitizing polygons around each training site for similar land cover.Then a unique identifier is assigned to each known land cover type [12].At the beginning, the training sites have been developed based on the collected reference data and ancillary information.

Signature Development
This is the stage of creating the spectral signature for each type of land cover.This is done by analyzing the pixels of the training sites.When the digitization of training sites is finished, the statistical characterizations of each land cover class are needed.These are called signatures [12].

Classification
After developing signature files for all the land cover types the images have been classified using a hard classifier called "Fisher Classifier".Fisher classifier uses the concept of the linear discrimination analysis [12].Fisher Classifier performs well when there are very few areas of unknown classes and when the training sites are representative of their informational classes [12].This is why fisher classifier is appropriate for this particular research, because most areas for the classes are known.

Generalization
After image classification, sometimes many isolated pixels may be found [12].These isolated pixels may belong to one or more classes that differ from surrounding pixels.Therefore it is necessary to generalize the image and remove the isolated pixels.Filtering is the solution for this type of problem.Therefore a 3 × 3 mode filter has been applied to generalize the fisher classified land cover images.This post-processing operation replaces the isolated pixels to the most common neighboring class.Finally the generalized images are reclassified to produce the final version of land cover maps for different years (Figure 2).It is important to remember that the supervised classification in every year is independent of the classification of the other years, because this method can lead to substantial overestimates of land change in the post-classification change analysis.

Accuracy Assessment
The next stage of image classification process is accuracy assessment.It is not typical to ground truth each and every pixel of the classified image.Therefore some reference pixels are generated [13].The randomly selected points within the classified image list two sets of class.The first set of class values represents the land cover type in the map.The second set of class values are known as reference values.These reference values are input by the researcher that is based on ground truth data.Therefore 250 reference pixels have been generated for each classification image for this research to perform accuracy assessment.Then the collected base maps have been used to find the land cover types of the reference points.

User's, Producer's and Overall Accuracy
User's accuracy for category K is the percent of category K in the reference information, given that the map shows category K. Producer's accuracy for category K is the percent of category K in the map, given that the reference information shows category K [14].The producer's and user's accuracy for all the years are found ranging approximately from 71% to 100% (Table 3-5).The overall accuracy represents the percentage of correctly classified pixels [14].It is achieved by dividing the number of correct observations by the number of actual observations.The overall accuracies for 1989, 1999 and 2009 are found 85.20%, 86.80% and 91.60% respectively.

Map Error vs. the Amount of Differences among the Maps
There is 14.8%, 13.2% and 8.4% reported map error in the maps of 1989, 1999 and 2009 respectively.These errors in the maps are clearly not negligible.It is a very common problem in land change science that the amount of error in the maps is nearly as large as the amount of change on the ground.Thus it is necessary to compare the error in the maps with how much difference there is during 1989-1999 and 1999-2009.
The amount of land change ranges from approximately 2% to 20% for 1989-1999, while the figure is approximately 1% to 17% for 1999-2009 (Tables 6 and 7).Therefore, on an average, it can be stated that the errors in the maps are not much larger than the amount of land change between the two points in time (1989-1999 and 1999-2009).Table 6.Amount of differences (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999).For older images it is not possible to visit the field to find out the actual land cover types.These things can be improved by recent field visit to collect Global Positioning System (GPS) data for land cover verification/ground truthing purpose.To solve the problems with older satellite images, the historical base maps of similar years should be collected.
The overall accuracies of the base maps can be achieved better if the above mentioned limitations can be reduced.

Land Cover Change Detection Analysis
In remote sensing, "Change Detection" is defined as the process of determining and monitoring the changes in the land cover types in different time periods.It provides the quantitative analysis of the spatial distribution in the area of interest [15].To perform change detection analysis, class-level metrics have been used.Class-level metrics are integrated over all the patches of a given type (class) [16].For this research purpose, 8-cell patch neighborhood rule has been applied.The 8-cell rule considers all 8 adjacent cells, including the 4 orthogonal and 4 diagonal neighbors.

Number of Patches (NP)
NP equals the number of patches of the corresponding patch type (class).Number of patches of a particular class or land cover type is a simple measure of the extent of subdivision or fragmentation of the class [16].The numbers of patches (urban blocks in this case) of built-up area and fallow land types have increased over time.Again the decrease of water body and vegetation is prominent (Figure 3).This reveals the gradual development of urban infrastructures converting water bodies, vegetation and low lands.

Gains and Losses in Land Cover Types
In case of built-up area, the core southern part of Dhaka city has remained the same.While the north-east and south-west parts have converted to built-up areas.The northern part of Dhaka city has gained water body followed by a massive decrease in the south-east and south-west parts (Figure 4).No particular pattern on gains or losses is found for vegetation.In cases of low land the changes are evident in eastern and western parts.Fallow land has decreased markedly and the losses are clear in north-western and mid parts (Figure 4).

Stochastic Markov Model
The first model that has been implemented is given the name as "Stochastic Markov Model (St_Markov)", because this model combines both the stochastic processes as well Markov chain analysis techniques [17].This kind of predictive land cover change modeling is appropriate when the past trend of land cover changing pattern is known [12].
A Markov chain is a stochastic process (based on probabilities) with discrete state space and discrete or continuous parameter space [18].In this random process, the state of a system s at time (t+1) depends only on the state of the system at time t, not on the previous states.

Markov Property
In a Markov chain the probability of the next state is only dependent upon the current state.This is called Markov property and stated as [19]: The probability of a Markov chain ξ 1 , ξ 2 ,........can be calculated as [19]: The conditional probabilities: | These are called the "Transition Probabilities" of the Markov chain [19].

Transition Matrix for a Markov Chain
Let's consider a Markov chain with n states s 1 , s 2 ,.......,s n .Let p ij denote the transition probability from state s i to state s j , i.e., | The transition matrix (n × n) of this Markov process is then defined as [19]: Predictions of the future state probabilities can be calculated by solving the matrix equation [19]: With increasing time steps, a Markov chain may approach to a constant state probability vector, which is called limiting distribution [19]:

Analysis with Stochastic Markov Model
At the beginning, Markov chain produces a transition matrix (Table 8), a transition areas matrix (Table 9) and a set of conditional probability images by analyzing two qualitative land cover images (Figure 5) from two different dates 1989 and 1999 [12].

Built-up Area Water Body Vegetation Low Land Fallow Land
The matrix of transition probabilities (Table 8) shows the probability that each land cover category will change to other categories in 2009.Table 9 represents the number of cells/pixels (30 m × 30 m) that will be transformed over time from one land cover type to other types.Markov Chain Analysis also produces related conditional probability images (Figure 5) with the help of transition probability matrices.Each conditional probability image shows the possibility of transitioning to another land cover class.
After analyzing Figure 5, it is clear that most areas will convert into built-up areas.The Markovian conditional probability of being built-up area ranges up to 0.66 which is the highest among all other land cover types.This probabilistic prediction is dependent upon the past trend of the last ten years (1989-1999).

Future Prediction Using Stochastic Markov Model
The next step is to make one single land cover map for future prediction aggregating all the Markovian conditional probability images.This prediction is performed by a stochastic choice decision model.Stochastic choice creates a stochastic land cover map by evaluating and aggregating the conditional probabilities in which each land cover can exist at each pixel location against a rectilinear random distribution of probabilities [12].
Stochastic Choice creates a stochastic land cover map by evaluating and aggregating the conditional probabilities in which each land cover can exist at each pixel location against a rectilinear random distribution of probabilities [12].The Stochastic Markov (St_Markov) predicted land cover map of 2009 is shown in Figure 6.

Cellular Automata Markov Model
The second model that has been implemented is named as "Cellular Automata Markov Model (CA_Markov).CA_Markov combines the concepts of Markov Chain, Cellular Automata (CA) [20], Multi-Criteria Evaluation [21] and Multi-Objective Land Allocation [12].
Stephen Wolfram defined CA as follows: "Cellular automata are simple mathematical idealizations of physical systems in which space and time is discrete, and physical quantities take on finite set of discrete values" [22].

The Elements of Cellular Automata
The components that comprise an elementary cellular automaton are as follows [23]:

Mathematical Notation of Cellular Automata
A CA model represents discrete dynamic system consisting of four elements [20]: where, L = the discrete lattice or the physical environment; ∑ = the set of possible states, where each i th cell of the lattice at time step t has a state σ i (t) ∑; = the neighborhood of a cell automaton, which is defined as all cells that fall within a radius r around the actual cell.It denotes the idea of neighborhood template [24]: , , , where , ; = the relative index of all neighbors of a particular cell; δ = the local transition rule, which is denoted as follows [20]: It shows that the state of the i th cell at the next time step t+1 is computed by δ based on the states of all the cells in its neighborhood at the current time step t.
Here, i (t) = the associated neighborhood with i th cell at time t; | | = the number of cells in the neighborhood.
The local transition rule is given by a rule table where given the sizes of ∑ and , the total number of possible rules equals [20]: | | (10) where each of the |∑ N | possible configurations of a cell's neighborhood is mapped to the number of possible states a cell can be in.Now considering the ordered set of all the states of all cells collectively at time step t, a CA's global configuration can be denoted as follows [20]: (11) Now applying the local transition rule to all the cells in the CA's lattice, the next configuration of the CA can be computed by its induced global map [20]: (12) In brief, the standard CA can be generalized as follows [24]: where, S = the set of all possible states of the cellular automata; N = a neighborhood of all cells providing input values for the function f and f = a transition function that defines the change of the state from t to t+1.

Simulation with Cellular Automata Markov Model
CA_Markov is useful for modeling the state of several categories of a cell based on a matrix of Markov transition areas; transitional suitability images and a user defined contiguity filter.A Markov model applies contiguity rule like a pixel near to an urban area is most likely to be changed into urban area [12].In this research, a 3 × 3 mean contiguity filter has been used (Figure 7) for modeling purpose.This kind of filtering is applied on suitability images for each land cover class.This is a defined neighborhood.The suitability of a pixel is determined by the pixel values within this defined filtering kernel.The more pixels of the same category of land cover exist in the neighborhood, the more the suitability value for that particular land cover type increases.Otherwise the pixel value remains the same [12].

Suitability Maps for Land Cover Classes
The suitability maps determine which pixels will change as per the highest suitability of each land cover type.The higher the suitability of a pixel, the possibility of the neighboring pixels to change into that particular class is higher.
Preparing a suitability map for each land cover type is difficult in terms of data and information availability.It is not possible to incorporate all types of factors or constraints that exist within the study area.Therefore a simple assumption has been assumed for fuzzy factor standardization.
The basic assumption for preparing suitability images is the pixel closer to an existing land cover type has the higher suitability.It means a pixel that is completely within vegetation has the highest suitability value (255) and pixels far from existing vegetation pixels will have less suitability values.The farthest pixels from vegetation will show the lowest suitability values.Here the suitability decreases with distance.Though this idea is not perfect always, like in case of water body and vegetation the scenario can be different.In reality, a pixel near to forest is more likely to convert into built-up area rather in forest.But this can be considered to be the research limitation.
Therefore a simple linear distance decay function is appropriate for this basic assumption.It serves the basic idea of contiguity.The land cover maps have been standardized (Figure 8) to the same continuous suitability scale (0-255) using fuzzy set membership analysis process [12].

Future Prediction Using Cellular Automata Markov Model
At the end, the Markov transition area matrix (Table 9), all the suitability images (Figure 8), the 3 × 3 CA contiguity filter (Figure 7) and the base map of Dhaka city (1999) have been used to predict the land cover map of 2009.The CA_Markov predicted final land cover image (2009) of Dhaka city is illustrated in Figure 9.

Multi Layer Perceptron Markov Model
The term "Artificial Neural Network (ANN)" has been inspired by human biological nervous system [25].In a typical ANN model, simple nodes are connected together to form a network of nodes.Some of these nodes are called input nodes; some are output nodes and in between there are hidden nodes [26].Multi Layer Perceptron (MLP) is a feed-forward Neural Network with one or more layers between input and output layers.The great advantage of using MLP perceptron neural network is that it gives the opportunity to model several or even all the transitions at once [12].

The Feed-Forward Concept of Multi Layer Perceptron Neural Network
MLP neural network uses the back propagation (BP) algorithm.The calculation is based on information from training sites [12].Back propagation involves two major steps, forward and backward propagation.The input that a single node receives is weighted as: ∑ where, w ij = the weights between node i and node j; O i = the output from the node i The output from a given node j is computed as [26]: f = a non-linear sigmoid function that is applied to the weighted sum of inputs before the signal passes to the next layer This is known as "Forward Propagation".Once it is finished, the activities of the output nodes are compared with their expected activities.In normal circumstances, the network output differs from the desired output (a set of training data, e.g., known classes).The difference is termed as the error in the network [26].The error is then back-propagated through the network.Now the weights of the connections are corrected as follows [12]: η = the learning rate parameter; δ j = an index of the rate of change of the error; α = the momentum parameter.
The process of the forward and backward propagation is repeated iteratively, until the errors of the network minimized or reaches an acceptable magnitude [26].The purpose of training the network is to get proper weights both for the connection between the input and hidden layer, and between the hidden and the output layer for the classification of unknown pixels [12].Several factors affect the capabilities of the neural network to generalize [26].These include:

Number of Nodes
In general, the larger the number of nodes in the hidden layer, the better the neural network represents the training data [26].The number of hidden layer nodes is estimated by the following equation [12]: (17) where, N h = the number of hidden nodes; N i = the number of input nodes; N o = the number of output nodes

Number of Training Samples and Iterations
The number of training sample also affects the training accuracy.Too few samples may not represent the pattern of each category while too many samples may cause overlap.Again too many iterations can cause over training that may cause poor generalization of the network [12].Over training can be prevented by early stopping of training [25].The acceptable error rate is evaluated based on the Root Mean Square (RMS) Error [25]: where, N = the number of elements; i = the index for elements; e i = the error of the i th element; t i = the target value (measured) for i th element; a i = the calculated value for the i th element.

Multi Layer Perceptron Markov Modeling
The basic concept of modeling with MLP neural network adopted in this research is to consider the change in built-up area over the years.In general, it means other land cover types are primarily contributing to increase the built-up area.At this stage, the issue of which variables affect the change to built-up area (1989-1999) has been considered.Therefore, only the transitions from "water body to built-up area", "vegetation to built-up area", "low land to built-up area" and "fallow land to built-up area" have been considered for model simulation.These four transitions have been termed as "All" here.Figure 10 exhibits the transition from all to built-up area.It is logical that new areas will be converted to built-up area where there are existing built-up areas.Therefore six driver variables have been selected for MLP_Markov modeling.These are (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999): Distance from all to built-up area, distance from water body, distance from vegetation, distance from low land, distance from fallow land and empirical likelihood image).
The empirical likelihood transformation is an effective means of incorporating categorical variables into the analysis (Figure 11).It has been produced by determining the relative frequency of different land cover types occurred within areas of transition (1989 to 1999).The numbers (legend) indicate the likelihood of changing into built-up area.The higher the value the likeliness of the pixel to change into the built-up cover type is more.Now it is important to test the potential explanatory power of each variable.The quantitative measures of the variables have been tested through Cramer's V [27].It is suggested that the variables that have a Cramer's V of about 0.15 or higher are useful while those with values of 0.4 or higher are good [12].
After getting satisfactory Cramer's V values for all the driving variables, now the turn is to run MLP neural network model.For this purpose, 10,000 iterations have been chosen.The minimum number of cells that transitioned from 1989 to 1999 is 4,794.Therefore, the maximum sample size has been chosen as 4,794.For each principal transition particular weights have to be obtained.The RMS error curve has been found smooth and descent after running MLP neural network.After all these combinations, the MLP running statistics gives a very high accuracy rate of 91.36% (this accuracy is a measure of calibration, not validation).Based on these running statistics the transition potential maps have been produced (Figure 12).These maps depict, for each location, the potential it has for each of the modeled transitions [12].These are not the probability maps where the sum of values for a particular pixel location will not be 1.The reason behind this is because the MLP neural network outputs are obtained by applying fuzzy set to the signals into values from 0 to 1 with activation function (sigmoid).Here the higher values represent a higher degree of membership for that corresponding land cover type [12].

Future Prediction Using Multi Layer Perceptron Markov Model
Using this kind of MLP neural network analysis it is possible to determine the weights of the transitions that will be included in the matrix of probabilities of Markov Chain for future prediction.The transition probabilities are shown in Table 10.Based on all these information from MLP neural network, the final land cover map of 2009 (Figure 13) has been simulated through Markov chain analysis.The whole procedure, for predicting the land cover map by this way, has been termed as "MLP_Markov" model.

Model Validation and Selection
Now the task is to select the most appropriate model.In general sense, model validation refers to comparing the simulated and reference maps.But the traditional way of validating a model or comparing two maps, using Kappa statistics, is now out-of-date [28,29].
Therefore, a method of comparing three maps (a reference map of time 1, a reference map of time 2 and a simulation map of time 2) has been implemented for model validation [30].In this case; the base map of 1999, the base map of 2009 and the simulated maps of time 2009 (St_Markov, CA_Markov and MLP_Markov) have been used.

Components of Agreement and Disagreement
The three maps comparison method consists of two components of agreement and three components of disagreement.The components of agreement are persistence simulated correctly and change simulated correctly [30].The components of disagreement are change simulated as persistence (the entries where reference t1 matches simulation t2 but does not match reference t2), persistence simulated as change (the entries where reference t1 matches reference t2 but does not match simulation t2) and change simulated as change to wrong category (the entries where all three maps disagree) [30].

Validation Results
Among the three models, the percentages of disagreement components are lowest (28.066%) while the percentages of agreement components (71.934%) are found highest for MLP_Markov (Table 11).This is how, implementing the three maps comparison method, it is found that the simulated map of "MLP_Markov 2009" is showing the best results in terms of the percentages of disagreement and

Simulating the Land Cover Map of 2019
The base maps of 1999 and 2009 have been used to predict the land cover map of 2019 (Figure 15).The procedure followed here is the same as stated in the MLP_Markov modeling section.

Analysis of the Final Predicted Map (2019)
The predicted map of 2019 reveals that 58% of the total area will be occupied by the "built-up area" cover type (Figure 16).On the other hand, water body (6%) and fallow land (11%) types are going to decrease in a notable way.

Conclusions
This research has revealed that the increase in built-up area is prominent in Dhaka City over the years .Later this paper has presented three different methods to simulate the land cover map of 2009 being persistent with the inherent changing characteristics.The methods have been named as "Stochastic Markov (St_Markov)", "Cellular Automata Markov (CA_Markov)" and "Multi Layer Perceptron Markov (MLP_Markov)" model.Then the MLP_Markov model has been found most appropriate, based on the three map comparison method, for future prediction (the land cover map of 2019).
Our hope might be realized if the error in the base maps is reduced to the point where the error becomes smaller than apparent change in land.We hope the interpretation of depicting the future scenario in quantitative accounts, as demonstrated in this research, will be of great value to the urban planners and decision makers, for the future planning of a modern and habitable Dhaka City.Moreover, it is our belief that this kind of research has a high potential to contribute towards the sustainable urban growth both at any local and regional level in the world.

Figure 1 .
Figure 1.(a) Location of Dhaka City in Bangladesh, (b) location of the study area within greater Dhaka City; Map 1 (a) Source: Banglapedia, National Encyclopedia of Bangladesh, 2011.
collected from the Dhaka City Corporation (DCC).These reference data have been used for training site selection and accuracy assessment of the final land cover maps.

Figure 2 .
Figure 2. Land cover maps of the study area.

Figure 3 .
Figure 3. Number of patches over the years.
(a) The physical environment or the space represented by an array of cells, on which the automaton exists (its lattice).(b) The cell in which the automaton resides that contains its state(s).(c) The neighborhood around the automaton.(d) Transition rules that describe the behavior of the automaton.(e) The temporal space in which the automaton exists.

Figure 12 .
Figure 12.Transition potential maps from all to built-up area (1989 to 1999).

Figure 16 .
Figure 16.Percentages of presence of land cover types over the years (1989-2019).

Table 1 .
Details of the land cover types.

Table 2 .
Details of the Landsat satellite images.(Source: US Geological Survey, 2011).
Class NameReference Totals Classified Totals Number Correct Producer's Accuracy User's Accuracy

Table 7 .
Amount of differences(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009).Limitations in Base Map Preparation and Accuracy Assessment Finally it can be stated that few misclassifications have been observed in the classified land cover maps of Dhaka city.The reasons may be as follows:i.The same spectral characteristics of some land cover types.For example, in case of 1989 base map, certain built-up areas were misclassified as fallow land.Again, in most cases it was really difficult to separate water bodies and low/cultivable lands categories.The reasons may be the seasonal variations of the satellite images for different years and the similar spectral properties of land covers in some cases.The images collected for 1999 (November) and 2009 (October) are from the same winter season.But the image of 1989 (February) is Therefore other seasonal images can be important evaluating the land cover change pattern of this kind of highly dynamic urban environment like Dhaka.iii.Again the spatial resolution of the images is important.For this research purpose, Landsat satellite images have been chosen that are only commercially available but can be found in free public-domain.The main problem of working with Landsat images is low resolution.The spatial resolution of Landsat Image is 30 m. IKONOS, QuickBird or other satellite images with higher resolution can be better option.iv.The next limitation regarding this research is the collection of reference data or maps.The reference data are necessary for ground truthing purpose of the base maps (1989, 1999 and 2009) that have been prepared from the Landsat satellite images.But reference maps of the respective years(1989, 1999 and 2009)are not available.Therefore the base maps of Dhaka city of the years 1987, 1995 and 2001, collected from Survey of Bangladesh (SoB), have been used for referencing purpose.Google Earth images (2010) are used as reference data for ground truthing the base map of 2009.v.This kind of research needs extensive field visit for image classification and assessing the accuracies of the satellite images.Another point is the verification of the older images.
from another season, summer.This kind of variation creates problems while preparing base maps for analysis.ii.Moreover, less image spectral resolution has directed to spectral mixing of different land cover types.This has caused spectral confusions among the cover types.It is also important to mention that the images of 1999 and 2009 represent winter season while the image of 1989 represents spring season.

Table 10 .
Transition probabilities grid for Markov chain (1989 to 1999) in MLP modeling.

Table 11 .
Components of agreement and disagreement for model validation.