A spatial zoning approach to calibrate and validate urban growth models

ABSTRACT Calibration and validation of models predicting urban growth have been largely developed using internal variables. Further investigation is required to improve model’s calibration and validation mixing internal and external variables. To reach this objective, a spatial zoning approach simulating long-term expansion of Mashhad, the second largest city of Iran, was presented in this study. Spatial zoning approaches distinguish local-scale urban dynamics in districts with different socioeconomic characteristics. Thiessen polygons were used to identify districts with different morphology and functional attributes. Urban growth was subsequently simulated for each district using a Multi-Layer Perceptron (MLP) neural network and Markov chains (MC) analysis. MLP and MC algorithms were respectively used to derive transition maps from non-urban to urban use of land and to determine spatial evolution of built-up areas at the metropolitan scale. Results of simulations based on spatial zoning were compared with outcomes of traditional urban growth models. Spatial zoning improved significantly model’s accuracy in respect to more traditional simulation modes. The approach proposed here is appropriate when simulating land-use changes under discontinuous urban expansion.


Introduction
To ensure urban sustainability, decision makers and planners need precise information on recent urbanization trends, especially around large cities and metropolitan regions (Jiang and Yao 2010). Recent studies (Hu and Lo 2007, Huang et al. 2009, Dubovyk et al. 2011) emphasize the importance of spatiotemporal analysis of urban expansion and the relevance of integrated approaches monitoring changes in peri-urban areas, with special regard to land take and soil sealing (Poelmans and Van-Rompaey 2009). Understanding patterns and processes of land-use change is a key issue when implementing dynamic models of urban growth (Jokar et al. 2013). Future scenarios of urban growth are important for planners, policy makers and stakeholders interested in land-use change (Lambin and Geist 2006).
Landscape transformations involve complex processes that are influenced by dynamic, nonlinear human-nature interactions, which are difficult to be monitored through standard approaches (Vega et al. 2012, Kolb et al. 2013. Various approaches have been implemented to simulate land-use change and predict urban expansion in different socioeconomic contexts (Cao et al. 2014, Mozumder and Tripathi 2014, Olmedo et al. 2015. Urban growth models are usually classified in three categories: (i) models implementing Markov chain (MC) or general regression approaches (Bell 1974, Clark 1965; (ii) dynamic models implementing Cellular Automata (CA), agent-based models, weights of evidence, artificial neural networks (ANN), System Dynamic (SD) models and colony optimization approaches (ACO); (iii) integrated and mixed approaches (Matthews et al. 2007, Shen et al. 2007, Rafiee et al. 2009, Thapa and Murayama 2011. To improve model's calibration and validation, methodologies integrating different algorithms have been suggested (e.g. Liu and Li 2007, Yang and Li 2007, Qiu and Chen 2008. Hybrid models have been recently proposed considering together empirical, statistical and dynamic approaches, such as ANN and MC (Thapa and Murayama 2012), CA and linear regression (Jokar et al. 2013), MC and linear regression (Geneletti 2012), CA and fuzzy logic (Wu 1998), conversion of land-use and its effect model (CLUE) and SD (Luo et al. 2010).
On the one hand, CA has demonstrated to have potential in simulating the spatiotemporal characteristics of complex systems (Li and Yeh 2002); Clancy et al. 2010). On the other hand, ANNs investigate nonlinear relationships and multiple (changing) thresholds and can be conceived as complex mathematical functions that convert input data to a desired output (Yang et al. 2008, Pijanowski et al. 2014, Eastman 2015. Although ANNs may answer questions such as 'which pixels are changing' based on land transformation matrices, they are less effective in simulating intensity and direction of change. Therefore, land-use changes have been better estimated using empirical and dynamic models such as MC or SD (Shen et al. 2007, Luo et al. 2010. In this sense, MC analysis is a convenient tool for modeling land-use when landscape modifications are latent or too difficult to be assessed with a traditional change detection analysis (Benito et al. 2010). MC analysis serves as a method assessing direction and magnitude of change thanks to the intrinsic capability of prediction of both linear and nonlinear trends in a multiple set of variables, such as landscape metrics.
Although integrated models may accurately simulate land-use changes, they often consider relevant variables as space-invariant. In fact, urban growth patterns vary largely in a given landscape according to the spatial variation of certain drivers of change. In this sense, our study introduces a methodology based on spatial zoning with the aim of improving calibration and validation of urban growth models. Spatial zoning operates partitioning a given space into homogeneous units according to specific criteria (Moreno-Regidor et al. 2012). The idea underlying the notion of spatial zoning is that each district is homogeneous for urban patterns, landscape composition and rate of change in the use of land along a defined time interval. These place-specific attributes influence model's capability in providing accurate scenarios of future urbanization. Appropriate discrimination of urban zones based on homogeneous criteria contributes to stable and reliable urban models. Various techniques were applied to zoning natural or urban areas including automated zoning procedure (Williams 2002), mathematical programming models (Shirabe 2005), graph partitioning models (Assunção et al. 2006, Tavares-Pereira et al. 2007, cluster analysis (Ochoa et al. 2009), artificial intelligence such as genetic algorithm (Demetriou et al. 2013), artificial immune systems (Liu et al. 2011), ant colony optimization (Li et al. 2011), simulated annealing (Santé et al. 2016), multiobjective hybrid metaheuristic approaches (Wei and Chai 2004) and spatial multi-criteria analysis (Habtemariam and Fang 2016). Many of the mentioned models are efficient only when applied to problems with a low number of zones (M) and/or basic spatial units (N), or with a set of less restrictive spatial conditions, as in cases where the centers of a given regions are predetermined (Moreno-Regidor et al. 2012). Thiessen (Voronoi) polygon is a geometric method that has been largely used to partitioning geographical space into different zones thanks to the joint use of spatial information and computational geometry: Casaer et al. (1999) applied this method with the aim to estimate wildlife spaceuse patterns; Morillo et al. (2016) used Voronoi polygons to zoning Madrid's districts in order to estimate housing prices; Moreno-Regidor et al. (2012) designed agricultural districts using adaptive additively weighted Voronoi diagrams. This method is particularly effective in partitioning a heterogeneous space into homogeneous urban zones or districts.
Zoning techniques have been marginally considered in land change modeling. To investigate urban growth at the county level in the United States, Jantz et al. (2010) have implemented a regional model of land-use change using the SLEUTH CA approach. Tayyebi et al. (2011) have modeled expansion of Tehran, Iran, partitioning the study area in three regions and developing separate models of growth.
Our study proposes an original framework integrating a spatial zoning approach (Thiessen polygons) into a growth model based on MC and ANN: land-use changes were analyzed implementing a transition matrix based on a MC approach; transition rules were based on a Multi-Layer Perceptron Neural Network (MLP-NN). This approach was used to simulate long-term expansion of Mashhad, the second largest city in Iran. Mashhad is a representative example of compact cities expanding under intense demographic pressure until early 1990s, and then shifting to a progressively dispersed morphology. The paper is organized as follows: Section 2 gives a brief description of the study area and the applied methodology, introducing the proposed model. Section 3 illustrates model's outcomes for Mashhad expansion providing a scenario for future land-use changes covering 12 years (2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025). Section 4 offers a comparative discussion of the results gathered in this study and provides conclusions and indications for future research.

Study area
To test the model proposed here, medium-term (2013-2025) urban growth has been simulated for Mashhad, the second largest city in Iran, situated in the northeast part of the country (Figure 1). Resident population in Mashhad amounted to 2.4 million inhabitants in 2006 and increased to 2.7 million inhabitants in 2011. Internal migration alimented urban growth producing uneven population gaps with other big cities like Isfahan, Tabriz and Shiraz. Mashhad is a tourism-specialized city with nearly 27 million visitors every year, 90 percent domestic. More recently, about 0.7 million refugees from Afghanistan settled in Mashhad creating new neighborhood units.

Data and input variables
Land change modeling used in this study was illustrated in Figure 2. To investigate landuse changes in Mashhad, three remote-sensed images (28 April 1986, 19 April 1999 and22 April 2013) derived from Thematic Mapper (Landsat 5) and OLI/TIRS (Landsat 8) sensors were considered. The maximum likelihood supervised classification algorithm has been used to extract and classify seven land-use classes in the area: (i) urban landuse (buildings and streets), (ii) vegetation (garden and agriculture lands), (iii) urban green space, (iv) barren land, (v) mountainous and rocky land, (vi) water surfaces and (vii) sedimentary surfaces. 'Overall accuracy' and 'Kappa coefficient' were used to evaluate classification accuracy. Kappa coefficient is a measure of agreement between model simulation and reality; it determines if the values contained in an error matrix may represent a spatial pattern significantly better than a random matrix (Wang et al. 2012). Results were compared with sample areas extracted from historical land-use maps dated 1986 and 1999 and obtained from Mashhad's municipality. Google Earth high resolution imagery was used to assess accuracy of 2013 classification. Kappa coefficients for classification results were respectively 0.89, 0.88 and 0.91 for 1986for , 1999for . Overall accuracy was 0.903, 0.896 and 0.930, respectively, for 1986for , 1999 Additional variables influencing urban growth were selected based on research  Table 1).

Spatial zoning
A problematic issue in urban growth models is model's calibration and validation is a problematic issue in urban growth simulations (Silva and Clarke 2002, Jantz et al. 2010, Wang et al. 2011, 2012. To improve accuracy in model's calibration and validation, the present study integrates a spatial zoning approach with a traditional land-use change model. Using a spatial zoning approach is justified from the fact that Mashhad urban districts display distinct growth patterns based on their socioeconomic characteristics, territorial potentials and biophysical constraints. Spatial zoning is intended to improve model calibration and validation and to obtain higher accuracy in the simulation of future urban growth. Spatial zoning implies that the study area can be segmented into homogeneous districts. For this purpose, administrative divisions or geographic units defined with specific criteria can be used. Earlier studies have simulated small-scale urban expansion at enumeration district, census track or cadastral parcel spatial levels (Pinto and Antunes 2010). In some cases, simulating urban expansion at detailed spatial scales has provided more realistic urban scenarios compared with those assured by standard aspatial models (Barreira- González et al. 2015b). However, administrative districts are often heterogeneous as far as urban growth patterns are concerned.
Thiessen polygons were used in this study to classify the study area into homogeneous districts as far as settlement density and dispersion are concerned. This methodology was first suggested by Thiessen (1911) as a way of spatial interpolation of biophysical variables. Thiessen polygons have been used in spatial partitioning problems across several disciplines, including geospatial science (Moreno-Regidor et al. 2012). This approach was also applied in economic analysis and urban planning, e.g. in market zone analysis (Okabe et al. 2000), retail trade zone design (Boots and South 1997) and service location problems (Okabe and Suzuki 1997). Thiessen spatial partitions make use of geometric distance to determine neighborhoods, producing a polygon around each target location that identifies all places matching with that target.
To classify the study area into homogeneous districts, pixels of each land-use class at time 1 (T1) and 2 (T2) were determined (30 m grid) with changes over time being calculated between T1 and T2. The required points for creation of Thiessen polygons were determined using the ArcGIS spatial analyst tool and visual analysis, considering density and dispersion of the changed pixels. To facilitate point selection, density of changed areas was calculated using a kernel density function. Locations with high and low density of settlements were selected as candidate district's centers. District boundaries were thus determined using a Thiessen polygons function based on selected points. To ensure the appropriateness of Thiessen polygons, visual analysis was applied to correct few nodes not representative of Mashhad's urban structure. Thiessen polygons separated dense urban zones from districts with dispersed settlements and discontinuous urban expansion. The spatial distribution of changed pixels was assessed using Moran's I index of spatial autocorrelation. In each zone, clustered or dispersed patterns indicate appropriate discrimination from neighboring polygons. Adjacent zones with similar land-use distribution and urban growth rates were finally merged.

Model calibration and simulation
Urban growth was simulated separately for each district at T3 (Time 3: 2013) based on land-use data at T1 (1986) and T2 (1999). To validate model accuracy, urban growth over 1999-2013 was simulated based on 1986-1999 data. Urban growth was subsequently simulated for 2025 (T4) based on 1999-2013 data. Urban expansion was supposed to advance proportionally to the growth rate observed in the previous time period at each district. Local-scale growth drivers were identified using Cramer's V coefficient analysis. Cramer's coefficients were used to quantify the association between variables and land-use classes. A high value for Cramer's V indicates a strong relationship between two variables. Variables with a low Cramer's V can be removed from subsequent analysis (Mozumder and Tripathi 2014); in this way, variables poorly associated with urban growth were excluded from the model. This would facilitate model's calibration and improves the accuracy of the final model.
Simulation was done using the land change modeler (LCM) provided by IDRISI (TerrSet). Urban growth in T3 (2013) was assessed using two variables: the transition potential map (TPM) and the amount of land-use change. Taken as one of the most effective stochastic processes to assess changing-state probability (Jokar et al. 2013), a MC process was used to quantify the amount of change for each land-use class and time interval. MC models are considered effective tools for modeling land-use transformations, providing indicators for direction and magnitude of change (Benito et al. 2010, Eastman 2015. MC analysis elaborates a transition probability matrix (A): ; (1) where A ij (with both i and j = 1, 2, . . ., n) is the area changing from the i-th state to the j-th state during the investigated time interval, and n is the number of land-use classes.

Multi-Layer perceptron (MLP) neural network
A Multi-Layer Perceptron (MLP) artificial neural network was used to determine the potential of land to transition and spatial location of transitioning cells. An artificial neural network consists of a connected network of processing units that are modeled on basic properties of neurons in the human brain. Neural networks are nonlinear and can be conceived as a complex mathematical function that converts input data to a desired output. The MLP approach using Back-Propagation (BP) learning algorithms is one of the most widely used neural network models. A typical MLP network contains one input layer, one or more hidden layers and one output layer (Eastman 2015). These layers are used for data input, data processing and data output, respectively (Hu and Weng 2009). In our study input layers represent areas experiencing land-use changes and the related criteria. The hidden layer was used to identify the relations between changed pixels and criteria. The output layer is based on a TPM. Functioning of a MLP neural network involves two major steps, forward and backward propagation, to accomplish its modification of the neuron connection weights. During training, each sample is fed into the input layer and the receiving node sums the weighted signals from all nodes to which it is connected in the preceding layer. The following algorithm was used to calculate the input that a single MLP node j received: where net j refers to the input that a single node j receives; w ij represents the weight between node i and node j and O i is the output from node i. The output from a node i was calculated as follows: The function ƒ is usually a nonlinear sigmoidal function that is applied to the weighted sum of inputs before the signal passes to the next layer. Once the forward signal is finished, the activities of the output nodes are compared with the related expected activities. Except unusual circumstances, the actual output will differ from the desired outcome, the difference being associated with errors in the network. Errors are then propagated backward with weights for relevant connections corrected via a relation known as the delta rule: where η is learning rate, α is momentum factor and δ is computed error. The forward and backward passes continue until the network has learned the characteristics of all land-use classes. The purpose of network training is to get the proper weights for the connections between input and hidden layers, and between hidden and output layers. The number of the hidden layer nodes is estimated by the following equation: where N h , N i and N o are the number of the hidden, input and output layer nodes respectively (Eastman 2015).

Land-use change analysis
Land-use maps for the study area are shown in Figure 3. These contain seven classes:  1986-1999, 1999-2013and 1986-2013. Other classes except open land (1999-2013 have experienced negative rates of change.

Spatial zoning
To partition the study area into homogeneous districts, land-use maps for 1986 and 1999 were overlaid computing change-detection maps for conversion from non-urban classes to built-up areas; spatial density of changed areas were estimated using a kernel density function. District centers were selected using visual analysis and were used as input for creation of Thiessen polygons. To improve simulation accuracy, urban growth was predicted separately for each Thiessen polygon. Point selection and creation of Thiessen polygons was iterated several times until the appropriate zoning was obtained.
Figure 4(a) shows the density of urbanized areas for 1986 and 1999 along with the selected points and Thiessen polygons. The spatial distribution of transformed pixels is expected to differ between urban districts. For this purpose, we computed the Moran's I coefficient of spatial autocorrelation in the use of land based on the polygons created above (Figure 4(b)). Some adjacent districts have similar coefficients and were merged to reduce simulation runtime, obtaining a total of 18 districts with Moran's I coefficients illustrated in Figure 4(c) and Supplementary Materials, Table 2. Peripheral districts showed higher Moran's I coefficients compared with inner districts, outlining relevant differences in growth pattern and amount of change among urban zones.

Simulating urban expansion for 2013
Urban expansion was simulated separately for all districts using MC analysis and a MLP neural network ( Figure 5(b)). Table 2 provides a summary of the Markov probability matrix used for simulating the transition between initial and final cell state. The highest probability of conversion to built-up areas was observed for open land (4270 ha), vegetation (1685 ha), mountainous and rocky land (503 ha) and urban green space (412 ha). Since our aim is simulation of urban growth, changes to urban land-use were specifically considered in the analysis. Supplementary Materials Table 3 shows   0.14 absolute rates of change from non-urban classes to urban land-use calculated for each district. Intensity of change in the use of land for each district was different at the end of the simulation interval (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) in respect with the rate of change in previous time intervals. For example, the highest amount of open land converted to built-up area was observed in district 5 (1198 ha). The highest amount of land converted to built-up area respectively for vegetation, mountainous land, sedimentary land and urban green was observed in district 6 (211.2 ha), 14 (338.5 ha), 16 (34 ha) and 10 (82.7 ha).

Model calibration and change allocation
Model calibration for each district was carried out selecting important target variables and creating TPMs. Cramer's coefficients were used to quantify the level of association between selected variables and land-use class and to identify the most relevant factors of growth in each district. According to Mozumder and Tripathi (2014), values of Cramer's coefficient ranging between 0.15 and 0.4 are satisfactory; significant coefficients are marked in bold (Table 3). Distance to urban areas, settlement density and evidence likelihood (probability of change to urban land-use) were identified as key variables influencing urban growth in all districts. Distance to sedimentary land and slope were moderately associated to urban land-use and, based on low Cramer's coefficients, these variables were excluded from the model. Settlement density was an important variable in central districts (1,2,3,5,7,10,17). Distance to roads was relevant in both external (6,8,11) and internal city districts (4,9,12). Urban green space was a relevant variable in districts 1, 10, 12, 14, 15, 16 and 17. Distance to river was relevant only in districts surrounded by water bodies (18, 5, 16, 13, 10 and 11). Distance to vegetation influenced urban growth in the inner city and external districts where built-up areas are generally located adjacent to patches with natural vegetation. Elevation was recognized as a relevant factor shaping urbanization in districts situated in mountainous areas south and southeast of Mashhad (1,2,6,14,15,17 and 18).
After identifying important variables of urban growth for each district, MLP parameters were specified for extraction of TPMs. Number of input layers were chosen in accordance with effective variables for each district; number of hidden layers was considered equal to 20 and TPMs were generated as output layers; momentum factor was considered equal to 0.5. Criteria to stop training were based on a 0.01 root mean square error (RMSE), with 50,000 iterations and 100% accuracy. A sample of pixels (changed or unchanged pixels in the previous time interval, i.e. calibration period) was also selected for MLP training stage. At the end of each training, algorithm performance including RMSE and accuracy was evaluated and best training was selected. Accuracy >80% indicated a good MLP performance; values below 75% were considered not reliable (Eastman 2015). An accuracy index above 75% was obtained for all districts except 12 and 15.
Determination of a TPM for each class completes the training process. TPMs for both non-zoning and zoning simulation modes were achieved by merging all TPMs (transition maps from non-urban classes to urban class). High map values indicate higher probability of change compared to neighboring pixels (Eastman 2015). In the non-zoning mode ( Figure 5(a1)), pixels with the highest probability of change were situated close to inner city. In the zoning mode ( Figure 5(b1)) the highest density of pixels with the highest probability of change was found in district 5 and in sparse external city patches. The simulation process was based on the TPMs developed from the MLP neural network and the intensity of change derived from the MC analysis. Figure 5(a2,b2) show results of urban expansion simulation in 2013 in non-zoning and zoning modes, respectively.

Validation of MC analysis
To validate MC estimated changes in the use of land, actual changes in land area (conversion to built-up area) were derived based on cross tabulation of classification data (1999 and 2013) in each district (Table 4). Differences between simulation modes' output and reality were subsequently computed. This allows comparing the estimated quantity of changes using MC analysis in each district under zoning and non-zoning mode. The highest accuracy was reached with zero-differences; negative and positive values respectively indicate that simulation output was less or more than reality. Although the observed difference for non-zoning mode (−4.07) is less than the difference found using the zoning mode (−6.13), zoning mode performed better than nonzoning mode in all districts except 1, 2, 7, 8 and 11; differences observed in districts 1, 2 and 7 were quite moderate.

Kappa Index of Agreement (KIA)
The Kappa Index of Agreement (KIA) indicates the degree of agreement between simulated and actual changes (Pontius and Suedmeyer 2004). KIA was used to evaluate simulation accuracy and indicates the extent to which two maps agree in terms of location (Pontius and Cheuk 2006). A 0.69 score was observed in the study area with the non-zoning mode increasing to 0.75 when the zoning mode was applied.

Validation through map overlay
To evaluate spatial accuracy, simulation results were compared with reality through map overlay. Figure 5(a3,b3) show validation maps for non-zoning and zoning modes, respectively, in 2013 where Hits (green) indicate correct simulated change, Misses (red) indicate non-simulated change and False alarms (yellow) indicate false simulated change. Amount of Hits was 3295 ha in zoning mode and 2750 ha in non-zoning mode (Table 5), improving accuracy by 19.8%. The highest improvements were observed in districts 16, 14, 10 and 5. Ratios of Hits to actual changes and Hits to simulated changes were respectively 0.406 and 0.352 in non-zoning mode, increasing to 0.484 and 0.423 under zoning mode. In central districts (1,5,7,10,13,14,16), accuracy (Hits to actual changes ratio) under zoning mode was higher than under non-zoning mode. In the inner city districts (4, 9 and 15) and in the external city districts (2, 8 and 11) both ratios were found higher under non-zoning mode compared with zoning mode. Taken together, simulation for districts in the external city (8, 18, 6, 11 and 2) was less accurate than in other districts. By contrast, districts 16, 17, 12, 5, 14, 9 and 1 had more accurate simulations.

Simulating urban expansion in 2025
Spatial zoning was used for simulating expansion of Mashhad in 2025. The study area was partitioned into 88 zones based on density of changed areas over 1999-2013 ( Figure 6(a)). Adjacent zones with similar distribution of changed pixels were subsequently merged, obtaining a total of 25 districts. The shift from 18 to 25 districts is  motivated with changes in urban boundaries and socioeconomic characteristics during 1986-1999 and 1999-2013, reflecting disting growth rates and land-use change patterns. Figure 6(b) shows TPM derived from MLP neural network at the district scale. According to this map, adjacent land to mountainous areas southwest of the city and free land spatially separated from consolidated urban areas have more potential for urbanization. Figure 6(c) illustrates simulation of the city expansion in 2025. According to simulation, 5818 ha will add to urban land-use in the 2013-2025 time intervals. The major growth has occurred in patches physically separated from consolidated city. Urban land-use will expand mostly on open land (vegetation and agriculture land) along north, northeast, east and southeast directions. Mountainous land in southern and southwestern parts of the metropolitan region acts as a natural barrier for future urban expansion.

Discussion
While appropriate identification of urban districts allows improving model accuracy, effective methodologies classifying districts into homogeneous groups are still required. Many of the models and algorithms used to partitioning space into homogeneous zones (e.g. cluster analysis, graph partitioning, mathematical programming) are efficient only when applied to problems in which the number of zones (M) and/or basic units (N) is low, or with a set of less restrictive spatial conditions, as in cases where centers of a given region are predetermined (Moreno-Regidor et al. 2012). In this study, we used a Thiessen polygons method based on density and distribution of changed pixels in the calibration period. Thanks to the use of spatial information and computational geometry in the zoning process, determining number and boundaries of urban districts according to settlement morphology proved to be more flexible than using administrative divisions. However, regularity in polygon boundaries may prevent a precise matching of polygons with urban districts' boundaries. Using improved versions of Voronoi diagrams, e.g. power Voronoi diagrams, additively weighted Voronoi diagrams (AWVDs) or multiplicatively weighted Voronoi diagrams (MWVDs) by smoothing boundary of zones, may contribute to an effective identification of homogeneous districts. Another limitation in the use of Thiessen polygons in urban growth simulations is that changing polygons could impact final model's outcome; adding or removing a point in a specific location would influence the state of surrounding polygons, impacting simulation results in a given district. In this sense, points have to be chosen based on the spatial structure of the changed area. Unsuitable merging of neighboring districts can negatively influence model accuracy. To overcome this problem, chosen districts should have similar growth patterns that can be verified using a Moran's spatial autocorrelation analysis.
Because cities are expanding and their boundaries are changing over time, the position of district's centers may change. Differential point selection and polygon design according to urban dynamics would help predicting realistic urban spatial structures. Accuracy is another important issue in land change simulation when evaluating overall model performance. To evaluate performance of Thiessen polygons, overall accuracy should be compared with other approaches after modeling land change at districts generated using different zoning systems. Kappa indexes comparing observed accuracy relative to a baseline accuracy expected due to randomness, were used in this study. However, in practical applications, randomness is an uninteresting, irrelevant and/or misleading baseline. Also the standard Kappa and its computational variants result rather complicated, difficult to understand and unhelpful to interpret (Pontius et al. 2011). In this sense, allocation agreement and disagreement was regarded as useful criterion to validate simulation results. Allocation indexes are computed from map overlay (or cross-tabulation matrix) of reference maps and simulated maps. Overlay produces maps that show three types of information, namely Hits (agreement), False alarms (disagreement) and Misses (Olmedo et al. 2015). Simple indexes can be derived by linear composition of these variables with the aim at validating different approaches in land change simulation.
Based on empirical evidences illustrated in our study, Mashhad's districts were classified into three categories based on land-use patterns and amount of change: (i) fringe districts adjacent to open land and vegetation, displaying rapid urbanization rates; (ii) external city districts with moderate urbanization rate and a lower simulation accuracy compared with fringe districts; (iii) inner city districts experiencing moderate urban changes due to restricted availability of buildable land.

Conclusions
Land-use simulation is important for a variety of planning and management issues as well as for academic research (Deng et al. 2008). In this study, a spatial zoning approach was used to improve calibration and to increase accuracy of a model simulating urban growth in Mashhad. Under the hypothesis that quantity of change is proportional to urbanization rate in earlier time intervals, use of spatial zoning increased the accuracy of MC analysis, containing the impact of heterogeneity in land-use and urban growth patterns on simulation outcomes and contributing in the identification of place-specific, relevant drivers of urban growth.

Disclosure statement
No potential conflict of interest was reported by the authors.