Remote Sensing , Gis and Cellular Automata for Urban Growth Simulation

Cities are complex spatial systems and modeling their dynamics of growth using traditional modeling techniques is a challenging task. Cellular automata (CA) have been widely used for modeling urban growth because of their computational simplicity, their explicit representation of time and space and their ability to generate complex patterns from the interaction of simple components of the system using simple rules. Integrating GIS tools and remote sensing data with CA has the potential to provide realistic simulation of the future urban growth of cities. The proposed approach is applied to model the growth of the City of Montreal. Land use/land cover maps derived from Landsat data acquired in 1975 and 1990 were used to train a CA model which was then used to project the land use in 2005. A comparison of the projected and actual land uses for 2005 is presented and discussed.


Introduction
There are many characteristics of cities that place them among the most complex structures created by humans.From the local-scale behaviour of many individual objects, structured and ordered patterns emerge in the aggregate, such as peak-hour traffic congestions, and large spatial clustering of socioeconomic groups by residence (Torrens, 2000).Cities also exhibit several characteristics of complexity such as fractal dimensionality and self-similarity across scales (Batty & Longley, 1994).Modeling the dynamics and growth of cities may be a challenging task without the use of tools that are well suited for modeling complex systems.
Urban models and simulations are used for testing theories about spatial location, the interaction between land uses and related activities, and to improve the understanding of urban evolution so as to predict change.Although based on solid theories, traditional urban modeling has a number of weaknesses such as their poor handling of time-space dynamics, their implementation difficulties, and the failure of their top-down approaches to reproduce realistic simulations for the urban process (Torrens & O'Sullivan, 2000).Urban modeling has moved from theories and structures that describe land use and change in aggregate static terms, to more dynamic models of individual behaviour from which complex spatial structure emerges.Current urban simulation approaches are taking advantage of progress in information technology, data availability and complexity theory (such as Cellular Automata (CA), and artificial intelligence) to address the problems of earlier urban dynamics models (Batty, 2009).New approaches provide a rich environment for understanding the dynamics of the system, as well as offering some innovative simulation techniques such as agent-based and cellular automata models.These models are based on simulating aggregate urban dynamics, and are characterized by both simplicity in representing the structure of cities and the ability to embrace their complexity (Batty, 2009).
There have been many efforts (Batty & Xie, 1994a;Couclelis, 1985;Couclelis, 1988;Tobler, 1970;Wu & Webster, 1998;Alkheder and Shan, 2008) to combine different approaches for modeling and representing the change patterns of such dynamic phenomena.Geographic Information Systems (GIS) have been used for the storage and management of spatial data, and also for improving the mapping and visualization of spatial processes (Aronoff, 1989).However, most GISs have a limited capability for handling the temporal component of spatial data.Dynamic approaches, such CA and complexity theory, are becoming increasingly used to model the dynamics of spatial complexity of urban systems because of their simple structure, and their capabilities for generating complex patterns that from the local interaction of simple components of the system (Torrens, 2000).
Remote sensing imagery are an ideal data source for CA modeling because of their historical record, their ability to be easily converted to various land use/land cover categories, and their inherent cellular data structure.This paper implements a cellular automata prototype based on Landsat TM data.The principal objective is to provide a realistic simulation of the dynamics of urban growth for the City of Montreal, Quebec over the last three decades.

CA-Based Modeling of Urban Growth
The standard CA consists of a regular lattice in which the CA components exist and evolve over time.The lattices of CA can be any dimension; however, the grids used to model geographic space are generally two-dimensional.The lattice consists of identical cells with local interaction, and each cell can be in one of the states from a finite set of defined states.The cell states evolve synchronously in discrete time steps according to a set of local rules (Wolfram, 1984), in which the cell state at time (t +1) is a function of the cell state and the states of its neighboring cells at time (t).The neighborhood of a cell consists of the cell itself and all the cells at certain proximity of the cell, from which the cell draws input.Transition rules are generally in the form of <ifthen-else> statements, and are applied on input from a neighborhood template to evaluate the value of each cell (Torrens, 2000).Transition rules enable complex systems to be modeled using simple components that drive their dynamics (Batty, Xie, & Sun, 1999).Despite their simple construction, CA are capable of complex universal behaviour (Wolfram, 1984).
Cellular automata are dynamic and discrete models for complex natural systems that are studied in computability theory, mathematics, physics, and theoretical and microstructure modeling.CA are cell-based methods that can model two-dimensional space.Because of this feature, some geographers have adopted CA to simulate in land use or other geographical phenomena (White & Engelen, 2000).In particular, CA have become especially useful tools for modeling land use dynamics (Tobler, 1979).
Waldo Tobler was the first to consider using explicit cellular automata for modeling geographical phenomena (Tobler, 1979).Tobler proposed a cell-space model simulating the urban growth of Detroit region (Tobler, 1970).Applying his first law of geography that "everything is related to everything else but near things are more related than distant things" (Tobler, 1970), his model related the growth of a cell to the population of the same and neighboring cells in the immediately preceding time period.Following this research, Tobler started investigating the use of cellular automata for other geographic systems.Accordingly, he classified land use change models into five types (I, II, III, IV, V).In his study, Tobler used a type V model, which he classified as a dynamic model, to show that the land use of a given location g ij at time step t+Δt is a function F of land use at that location at time t and of a measure of all the land uses in neighbourhood n of that location.The neighbourhood used was Von Neumann neighbourhood and was defined as the geographical domain of influence.Transition rules were defined similar to the standard CA definition.After the pioneering work of Tobler, many research efforts explored CA and its application to geography.For example, Coulelis (Coulelis, 1985;Coulelis, 1988;Coulelis 1989) provided theoretical and practical frameworks for using CA in modeling the dynamics of geographical spatial phenomena.Further, Batty (Batty, Xie & Sun, 1999) developed the Dynamic Urban Evolutionary Model (DUEM) to model land use change through transition rules.They defined various decision rules that embedded distance and direction, density thresholds, and transition or mutation probabilities into the model's dynamics.White and Engelen (White & Engelen 1993;White and Engelen, 1994) integrated models for socio-economic and natural systems to predict the future demand for the various land uses.Clarke, Hoppen and Gaydos (Hoppen & Gaydos, 1997) coupled two CA models -an urban growth model and a land use change model -to predict urban growth as part of a project for estimating the regional and broader impact of urbanization.Their model started as a basic urban growth model using weighted maps as inputs and applied a set of CA rules that determined whether or not cells will change from nonurban to urban for each year in the sequence.This was coupled to a second cellular automaton that received the number of new urban cells in each time period and created and enabled change among land uses other than the urban class, for example from wetlands to agriculture.Wu and Webster (Wu & Webster, 1998) integrated multicriteria evaluation into CA simulation to define non-deterministic, multidimensional, and multilevel transition rules.Chen (Chen et al., 2002) developed CA model to investigate future growth scenarios for the city of Beijing from 1976 to 1997.In their model, they used an adaptive Monte Carlo method for the calibration of the factor weights used in CA transition rules.Shan (Shan et al., 2008) proposed the use of genetic algorithms for the calibration of CA transition rules.
The popularity of urban CA models is in part due to the weaknesses of traditional urban models, such as their poor handling of time-space dynamics, their static aggregate representation of change, and their unrealistic simulation of space (Torrens, 2000).Also, CA models offer several advantages that make them well recognised for the modeling and simulation of complex processes, such as urban growth dynamics.Some of these advantages are listed as follows (Torrens, 2000): The most appealing characteristic of cellular models is their simplicity, In CA models, realistic and complex patterns emerge from simple rules and the interaction of simple components at the local level.In CA models complexity emerges from within the model; it is not a part of the model design.CA models are also used as explanatory tools for real-world applications.CA are explicitly spatial and make implicit use of spatial complexity.The parallel processing approach of CA makes them decentralized.
The decentralized approach treats the system dynamics as nonlinear interactions among its individual components.This mirrors the ways in which real cities decentralize, and move away from monocentric to polycentric city systems.CA have intrinsic affinity with raster data, and seem to be well suited to remotely sensed data and GIS.Further, both CA and spatial information systems represent attributes in a layered fashion (e.g., class and spectral layers in remote sensing, themes in GIS, and cell state-spaces in CA) and manipulate that information with spatial operators (e.g., texture measures in remote sensing, neighbourhood operators in GIS, and transition rules in CA.The capability of CA to handle fine-scale dynamics with computational efficiency make them suitable for detailed simulation.This stems from the progress made in operational land-use and transport models, but with an explicit attention to detail.CA models have the ability to represent system dynamics.Traditional urban models are relatively static and treat urban dynamics very poorly.CA models represent a significant advance in the treatment of time, since they are interactively dynamic and move in short time steps.This makes them capable to approximate real-time urban dynamics.Also, CA allow multiple time scales to be represented facilitating their suitability for modeling the activities and events of cities that vary temporally, such as long-term economic cycles, daily activities, and hour -by-hour social interactions.CA models simulate the emergence of large-scale patterns from the interaction of local components.CA has much to offer for urban simulation, city is complex system composed of subsystems that interact with each other at different scales.CA models can connect these subsystems by allowing large-scale patterns to emerge from the interaction of local components.The highly visual simulation environments of CA models enhance the engagement and interaction of users with the model.Also, the visual dynamics of CA models allows the evolution of the system to be displayed as it changes through time.
In order to use CA models for urban simulation they must be modified from their basic standard form.The urban CA model consists of a rectangular grid of square cells.CA cells are considered as pixels of a land use map derived from a classified remote sensing image, except that each cell is able to process information and visualize the evolution of its state.As such, each cell is associated with vectors of suitabilities, one for each land use taking part in the dynamics.
The neighbourhood space of the urban automata is defined as a circular region around the central cell with a radius of eight cells.A neighbourhood effect is calculated for each cell to represent the effects of the various land uses and land covers within its neighbourhood.Distant cells will have a smaller effect in the evaluated central cell.
A cell can be in one of the states represented by the different land uses classes in the area.The state of a cell (i) at consecutive time t+1 is a function of its state (S), its neighbourhood effect (N), and set of transition rules (T) at an initial time t.Using this function, transition rules are applied to each cell to determine what state it should change to during a time transition.In land use change modeling, these rules represent how change occurs in the real world.The transition rules change each cell to state for which it has the highest potential.The cell state can be mathematically represented as: (Barredo, Kasanko, McCormick, & Lavalle, 2003)

Study Area
The CA model is applied to model the urban growth of the city of Montreal, Québec.The study area (Figure 1) covers the Greater Montreal Metropolitan Area and is approximately 3600 km 2 .The city got its name from the most prominent geographical feature in the island, The Mount Royal Mountain.According to statistics Canada, Montreal's census Metropolitan Area (CMA) (land area 4,258.97km 2 ) is recognized as the second largest populated region, with a population of 3,451,027 in the year 2001, and a population of 3,635,571 as of 2006 census, with growth rate of 5.3% (Statistics Canada, 2006 Census of Population).The population growth had significant affect in the rate of urban sprawl in Montreal in the last few decades (Statistics Canada).NAD 1983) and resampled to 100 m grid cell spacing.All image analysis was completed using PCI Geomatica V10.3.A 100 m grid spacing was used because: (i) we were not interested in land use changes at resolutions smaller than this; and (ii) this setting balanced the differences in spatial resolutions between the sensors.
Canadian Digital Elevation Data (CDED) was obtained from the Canada Centre for Topographic Information (Natural Resources Canada, 2000).The elevation data were mosaicked and registered to the Universal Transverse Mercator projection (UTM Zone 18, NAD 1983) and resampled to a 100 m grid cell spacing to match the satellite imagery.The percent slope for each grid cell was calculated for use in the CA model.Road and street data were obtained from Statistics Canada.These data had been pre-registered to a UTM projection.
Only main and major road classes were used because these were thought to have the greatest impact on urban growth.

Image Classification and Accuracy Assessment
The three Landsat images were independently classified using the maximum likelihood algorithm into six land use/land cover classes: water, residential, commercial-industrial, vegetation, cropland, and woodland.The residential and commercial-industrial classes represent the primary classes of interest for this study of urban growth.Land use maps were produced for the years 1975, 1990, and 2005.Figure 2 shows the land use maps for the years 1975, 1990, and 2005.
In order to use the produced maps for further analysis, an accuracy assessment was carried out.Accuracy assessment compares the agreement of the classified image/map with the reference data or ground truth.
One of the most commonly used forms of classification accuracy is the error, or confusion matrix.Many accuracy statistics are derived from the error matrix, such as the overall accuracy and the Kappa Index of Agreement (KIA) (Cohen, 1960).
KIA is a multivatiate analysis technique appropriate for the analysis of categorical discrete data such as the classification data (Conglaton et al.,1983, Conglaton, 1991).KIA measures the difference between the actual agreement between classified data and reference data and the chance agreement between the reference data and a random classifier.The Kappa index ranges between 0 and 1, where a value of 1 indicates a perfect agreement and a value of 0 indicates the level of agreement is due to chance agreement.
The overall accuracies of the maps of 1975, 1990, and 2005 were 94%, 80%, and 94% respectively; the overall Kappa statistics were 0.92, 0.75 and 0.92.The achieved accuracies were deemed to be adequate for effective and reliable land use change analysis and modeling.

CA Modeling of Urban Land Use: Model Description and Implementation
In this study we applied an integrated approach that combined cellular automata with Markov analyses (Bell, 1974;Brown et al, 2000;Muller & Middleton, 1994) for modeling land use change and project the future urban growth in the study area.We used the CA-Markov algorithm as implemented in the IDRISI Tagia GIS.Land use maps from different time periods were used to develop the urban probability statistics for land use change for a future time period.Although a stochastic land cover map can be created from the conditional probability images using stochastic choice decision model that evaluates the conditional probabilities for each land cover at each pixel.The result lacks spatial knowledge of the distribution of each class.This missing spatial information specific to the prediction of land use change is introduced through cellular automata.The CA allows the transition probabilities of one pixel to also be a function of its neighbouring pixels as defined by a contiguity filter.This filter is spatially explicit contiguity-weighting factor applied to the set of suitability maps created for each land use that express the suitability of each pixel for each of the considered land use types.A pixel close to specific land use is more likely to change to that category than a distant one.The filter reduces from the suitability of areas away from existing area of that type.A suitability map was created for each kind of land use (e.g., Figure 4 shows the suitability map for residential land use).

Figure 3. Conditional probability image for class 3
A multi-criteria approach was applied to create the land use suitability maps (Eastman, 2009).The approach combines various criteria into a single index that indicates the suitability of each location in the study area for that specific land use.Criteria are of two types: factors and constraints.Constraints exclude or constrain areas from being considered for a type of land use and factors are criteria that assign a relative degree of suitability for each location for the land use under consideration.For example, water bodies are not considered as suitable for residential development under any condition and proximity to roads is a factor that enhances the suitability of land for residential development.Constraints are represented by Boolean images, where 0 is assigned to areas that are not suitable and 1 is assigned to suitable areas.Factors define a degree of suitability and are represented by continuous scales / continuous images; different factors may have different units -e.g., distance from roads in meters and slope of land in degrees.In order to be able to combine factors, they are standardized to a scale of 0 to 255 using fuzzy functions (Eastman, 2009).Non-suitable areas are represented by 0 and the most suitable areas are represented by the value 255.The criteria applied to create the suitability maps were selected from the existing literature.For example, the criteria used here for residential land-use were land use type, distance from roads, slopes, distance from water bodies and distance from developed areas. In

Results and Discussion
A visual comparison between the actual map for the year 2005 and the predicated map for the same year was carried out to assess the accuracy of the predicated map.
The water areas were correctly predicted due to the constraint applied in the model (i.e., the water areas were considered as out of bounds for the growth of any other type of land use).The residential class was over-predicted in some areas and under-predicted in others.This is probably due to the effect of applying the contiguity filter.The filter increased the suitability of a pixel to become residential if the local neighbourhood pixels were residential, and the suitability value of the pixel to become residential (or remain residential) was decreased if the pixel was away from existing residential land.The same observation applies for the commercial-industrial class.The cropland appeared to take over the commercial-industrial land in areas where there were only isolated areas of commercial-industrial development and large contiguous areas of cropland.The vegetation land use class was noticeably over-predicted and suffered the same smoothing effect (growth of areas that were homogeneously vegetation) indicated above; in the actual map, vegetated areas tend to appear as small pockets scattered over the map area.Small areas of forest were lost (went to other classes) and larger areas had small inclusions of other land uses removed.An error matrix (Table 2) was created for the real and predicted maps for the year 2005.The real classes represented by the columns of the matrix are compared to the classes in the predicted map represented by rows.
In addition to giving an overall accuracy measure of the classification, error matrix provides more detailed information about the individual classes' accuracies.Having the accuracy of each category represented in the map can be very valuable when assessing the effectiveness of the classified map for a specific application (Story & Conglaton, 1986).
Classes' accuracies are calculated by dividing the number of correctly classified cells by the number of total of pixels in that class (as defined by the corresponding row or column totals).The producer's accuracy for a category is obtained by dividing the number of correctly classified pixels of that category by the total number of the category pixels in the reference data (column total).This measure is an indication of percentage of the ground data for certain category that was classified as such.On the other hand, the user's accuracy (related to errors of commission -) is obtained by dividing the number of correctly classified pixels in each category by the total number of pixels in the reference data classified into that category.The resulting percentage is an indication of the probability that what was classified as a certain category in the map is in fact that category on the ground (Story & Conglaton, 1986;Conglaton & Green, 2008) Derived from the error matrix (shown in Table 2), the user's and producer's accuracies for each class were computed (shown in Table 3).We can see that only 29% (producer's accuracy) of the commercial-industrial class was predicated correctly as that class and that the largest components of the error of omission (E O ) for the class are cropland (50%) and residential (16%); the largest portion of the error of commission (E C ) was residential (12%).For the residential class, 53% (producer accuracy) was predicted correctly as residential and the largest components of E O are cropland (24.5%) and commercial-industrial (12%); the largest component of E C was commercial-industrial (16% To evaluate the model and assess the overall degree of similarity of the simulated and actual land use maps for 2005, the two maps were compared using the Kappa Index of Agreement (KIA).Different Kappa indices were obtained: Kappa for no information (Kno) and Kappa for location (Klocation).Kno indicates the proportion classified correctly relative to the expected proportion classified correctly by a simulation with no ability to specify accurately quantity or location.Klocation is defined as the success due to a simulation's ability to specify location divided by the maximum possible success due to a simulation's ability to specify location perfectly.The overall accuracy of the predicted map was 59.4%, Kno was 0.51, and Klocation was 0.66.However, combining the two urban classes (residential and commercial-industrial) into a single urban class increases the overall accuracy of the prediction to 65%, Kno to 0.56, and Klocation to 0.73.
In addition to the effect of the contiguity filter on the predictive accuracy of the model, as discussed above, changes in economic and social factors between the two time periods (1975-1990 and 1990-2005) (the year for which land use is being predicted) may also contribute towards prediction error; for example, differences in the level of economic prosperity / activity and changes in development policies between the two time periods may lead to the calculation of probabilities that don't accurately reflect the growth in the later period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005).
Another major aspect, which may have considerable impact on the simulation results, is the availability of the data layers that represent the drivers of change (social and economic factors) in land uses and are necessary for the creation of more accurate suitability maps.

Conclusion
The CA prototype is anticipated to be able to provide a good representation for the dynamics of urban growth of cities.However, as can be seen from the results obtained from the model used, some additional refinement may be required to achieve better predictive accuracy.Such refinement may include the addition of more data layers, which may be a problem in some cases due to the lack of availability of the data.Additionally, other modifying factors such as economic performance may better help simulate the rate of growth of urban areas.Another factor that may be important is some refinement to the application of the contiguity filter, so that it more realistically handles areas of heterogeneous land use and avoids the tendency to smooth those areas away, creating large patches of homogenous land use.Also, the availability of reference data can better help assess the accuracy of the classification and simulation results.
The most challenging task is specifying the criteria that are subsequently applied to create the land use suitability maps.The work on this paper motivated the development of the ontology-based information extraction system, for the domain of land use suitability analysis, described in (Al-Ageili, 2014), to automate the extraction of the criteria and their values from bylaw and regulation documents related to the geographic area of interest.

Figure 1 .
Figure 1.Location of Study Area: Greater Montreal Area, Quebec, Canada

Figure 4 .
Figure 4. Residential land use suitability map this study, land use maps from the years 1975 and 1990 were input to the CA model.The model produced a transition probability matrix for the land use classes between 1975 and 1990, as well as a set of conditional probability images for each class of the land use classes defined (e.g., Figure 4 shows the conditional probability image for being commercial/industrial land use class).The transition probability matrix, the collection of suitability images (one image for each kind of land use) and the CA contiguity filter were used to project the land use change for the year 2005 (Figure 5 illustrates Projected Land Use Map for the Year 2005).

Table 1 .
). Landsat Images All the images were geometrically corrected to the Universal Transverse Mercator projection (UTM Zone 18, ).For vegetation, 62% was predicted correctly, with the largest components of E O cropland (22%) and forest (13%); the largest components of E C were cropland (29.5%) and forest (13%).Pixels that are cropland in the real map were the most accurately predicted (at 75%), with the largest component of E O being vegetation (13%); the largest components of E C were commercial-industrial (50%), residential (24.5%), and vegetation (22%).Forestland was 59% correctly predicted, with the largest component of E O also being vegetation (29.5%); the largest component of E C was vegetation (13%).It is apparent from the above that the majority of the problems with the predicted map are with respect to residential, commercial-industrial, and cropland, with cropland responsible (errors of commission) for most of the predictive error with respect to the two urban classes.

Table 2 .
Error Matrix