Simulating multiple urban land use changes by integrating transportation accessibility and a vector-based cellular automata: a case study on city of Toronto

ABSTRACT The accessibility provided by the transportation system plays an essential role in driving urban growth and urban functional land use changes. Conventional studies on land use simulation usually simplified the accessibility as proximities and adopted the grid-based simulation strategy, leading to the insufficiencies of characterizing spatial geometry of land parcels and simulating subtle land use changes among urban functional types. To overcome these limitations, an Accessibility-interacted Vector-based Cellular Automata (A-VCA) model was proposed for the better simulation of realistic land use change among different urban functional types. The accessibility at both local and zonal scales derived from actual travel time data was considered as a key driver of fine-scale urban land use changes and was integrated into the vector-based CA simulation process. The proposed A-VCA model was tested through the simulation of urban land use changes in the City of Toronto, Canada, during 2012–2016. A vector-based CA without considering the driving factor of accessibility (VCA) and a popular grid-based CA model (Future Land Use Simulation, FLUS) were also implemented for comparisons. The simulation results reveal that the proposed A-VCA model is capable of simulating fine-scale urban land use changes with satisfactory accuracy and good morphological feature (kappa = 0.907, figure of merit = 0.283, and cumulative producer’s accuracy = 72.83% ± 1.535%). The comparison also shows significant outperformance of the A-VCA model against the VCA and FLUS models, suggesting the effectiveness of the accessibility-interactive mechanism and vector-based simulation strategy. The proposed model provides new tools for a better simulation of fine-scale land use changes and can be used in assisting the formulation of urban and transportation planning.


Introduction
Urbanization is a global trend unique to the past few centuries, featured with continuous urban land expansion in urban fringes Yeh, Li, and Xia 2021) and multitype urban land use changes within cities . As suggested by many previous studies, urban expansions and urban land use changes are closely linked to the development of transportation system (Wegener 2021;Soteropoulos, Berger, and Ciari 2019). In fact, the urban transportation network not only provides opportunities for socioeconomic activities within/between cities but also contributes to shape the city forms and urban land use patterns over the long term (Iacono, Levinson, and El-Geneidy 2008;Wee 2004). Specifically, the transport supply from the transportation network determines the attractiveness of location to people's activities, for example, the choices for places to live, work, sport, and leisure, which therefore drives the urban growth and the changes the functional land use in urban areas (Bertolini 2017;Rodrigue, Comtois, and Slack 2016). Thus, urban transportation system has long been regarded as the fundamental driving force of the socioeconomic developments and urban land use changes, especially in metropolitan areas (Glaeser, Kahn, and Rappaport 2008;Zhang, Zhang, and He 2021).
Given that urban land use and transportation system are interacted together as a synthesis, a series of land use and transport interaction models have been designed to investigate the urban land use changes driven by accessibility provided by transportation infrastructure (Chang 2006;Ortuzar and Willusen 2001;Wegener 2021). Previous models can be roughly generalized into three categories: method based on regression analysis (Iacono and Levinson 2009), agent-based modeling (Salvini and Miller 2005;Wise, Crooks, and Batty 2016), and Cellular Automata (CA) land use simulations (Abolhasani et al. 2016). Among these models, the CA-based approaches can provide advantage in comprehending the complex dynamics of the urban growth and its interactions with the local environment. Conventional CA models are composed of five components (Santé et al. 2010;Verburg et al. 2004;Batty, Xie, and Sun 1999): a space containing a collection of homogeneous cells, a sequence of discrete time steps, a definition of neighborhood of each cell, a set of states describing the dynamics of cells, and a set of transition rules determining the changes of state of each cell. Because of its simplicity of computations and assumptions, CA models are immensely popular for simulating the spatiotemporal dynamics of land use changes (Wu and Martin 2002;Mohan Rajan and Loganathan 2021).
In the context of urban expansion and urban functional land use change, accessibility, which is known as a measure of easiness with which people can reach or interact with destinations (Hansen 1959), plays an essential role in connecting the transportation system and urban land use dynamics (Kelobonye et al. 2019;Liu and Zhu 2004). A series of studies suggested that the dynamic couple of accessibility into the land use model is a feasible way to investigate the evolution of urban land use driven by the transportation (Aljoufie et al. 2013;Liu and Zhu 2004;Kasraian, Maat, and van Wee 2017). In fact, many existing studies on urban change simulation have already taken into account the impact from transportation system in CA modeling by adopting transportation-related variables, such as distance to road networks and bus stations Verburg et al. 2004). Although great progress has been made previously (Chang 2006;Ortuzar and Willusen 2001), there are still limitations on the integration of accessibility in CA modeling, which needs further improvements for a better simulation of urban land use change driven by the transportation system.
The first limitation is that most traditional CA models usually interpretate the accessibility as simplified variables, for example, proximity to transportation facilities (Euclidean distance to road network or bus stations) (Feng and Tong 2018;Liu et al. 2010;Zhang et al. 2020). There is no doubt that the urban development is closely related to the proximity of land patch to roads and other travel facilities. However, these simplified proximity variables are incapable of reflecting actual people's movement and traffic flows, whereas the urban land use dynamics, especially the changes among different urban functional types (e.g. transition from residential to mixed land use or commercial use), are more directly linked to this kind of travel information. Some studies also pointed out that the routinely collected travel information, such as travel time and traffic volumes on the transportation network, is a more intuitive and rational representation in the modeling of interaction between urban transportation and land use (Luna, Miller, and Shoshanna 2018;Rad and Alimohammadi 2021). Unfortunately, there is limited research on the integration of such detail traffic information into the land use simulation, probably due to the difficulty of traffic data collection.
The second limitation is that conventional CA models adopt the grid-based representation (raster data format) in the urban land use change simulation because of its simplicity of data organization and convenience of computation (Aljoufie et al. 2013). These CA models typically take grid cells as their basic modeling unit without considering the fine-scale shape of ground objects. Therefore, they are usually considered inadequate for applications with fine resolution spatial data (Pinto, Antunes, and Roca 2017;Moghadam, Karimi, and Habibi 2018) and thus unable to simulate the spatial configuration of urban patches (Chen et al. 2019;Yao et al. 2017). In addition, the simulation of land use patterns from the grid-based CA model is usually sensitive to the cell size parameter and neighborhood configuration (Chen et al. 2014;Dahal and Chow 2015). Moreover, the subtle urban land use changes, which are closely related to human daily travel behaviors, are barely able to be captured and simulated with the regular geometry of grid cell in the CA simulation Stevens and Dragićević 2007). In fact, people living in urban areas usually schedule their destination and routine of activities according to the location's accessibility provided by the transportation network (Badoe and Miller 2000;Li et al. 2021). When introducing the accessibility into the land use simulation, the complexity of accessibility itself (e.g. travel time on the road network) requires fine-scale shape of land patches to better represent the actual travel cost in the transportation network rather than simplified Euclidean distances between the origins and destinations. Compared with the grid-based representation, the vector format data is more capable of capturing fine-scale ground objects with irregular polygons and thus can avoid the bias caused by the grid-based representation Chen et al. 2019;Barreira-González, Gómez-Delgado, and Aguilera-Benavente 2015). Recently, a series of patch-based and vector-based CA models were proposed to simulate fine urban land use changes and proved to be more realistic for revealing the temporal dynamics of ground object (Abolhasani et al. 2016;Dahal and Chow 2015;Barreira-González, Gómez-Delgado, and Aguilera-Benavente 2015;Chen et al. 2014;Lu, Cao, and Zhang 2015). Thus, the incorporation of vector-based simulation strategy into the modeling is expected to better characterize the interaction between the transportation system and land use dynamics, especially the subtle urban changes driven by the accessibility.
The motivation of this study is to overcome the two abovementioned limitations and develop a better model for the simulation of urban land use change driven by the transportation system. To achieve this, we proposed an Accessibility-interacted vector-based CA (A-VCA) model to simulate the subtle urban land use changes among different functional types. Different from previous models, we took the accessibility derived from actual urban traffic data as a key driver to investigate the urban land use changes driven by transportation infrastructures. In order to simulate the urban land use changes among different urban functional types, for example, transition from residential to mixed land use, we adopted the Random Forest Algorithm (RFA) to calibrate the A-VCA model since it is popular and effective for highdimensional data modeling (Belgiu and Drăguţ 2016;Yao et al. 2017). Specifically, the accessibility variables at both local and zonal scales were estimated based on using actual auto-in-vehicle travel time collected by the Transportation Tomorrow Survey (TTS) and was processed with the EMME transportation network analysis software (https:// www.inrosoftware.com/en/products/emme/). Then, the estimated accessibility variables were incorporated into the vector-based CA model to explore the transition rules and derive the land use transition probabilities. The proposed A-VCA model was applied and tested in the City of Toronto, Canada, during 2012-2016. Note that we mainly focused on the fragmented processes of urban land use change among different functional types within the urban area. The number of such land use changes is usually not large, and most of them mainly occur in discrete and fragmented land parcels within the urban area. Thus, simulations for a relatively shorter time span, that is, from 2012 to 2016 in this study, are more suitable to test the model's capability of simulating these fragmented and subtle urban land use changes among different functional types. Hopefully the model is expected to accurately simulate the subtle urban land use changes in the City of Toronto with good morphological features.
The rest of the article is organized as follows: In Section 2, we give a brief description of the study area (City of Toronto, Canada) and the data used to conduct the simulations and analyses. Section 3 presents the development of the proposed A-VCA model, including how to define and derive the accessibility (at both local and zonal scales), and how to embed the accessibility variable into a vector-based CA model. In Section 4, we apply the proposed A-VCA model in the study area and evaluate the performance of the simulation. A comparison of the proposed and existing models is also presented in this section. Finally, some conclusive remarks and limitations are presented in the Section 5. The results and findings in this study can deepen our undersetting of the interaction between land use changes and transportation, and provide useful knowledge for urban growth management and transportation planning.

Study area and data
The proposed model was tested through the simulation of urban land use changes in the City of Toronto, Canada ( Figure 1). The City of Toronto is the most populous city in Canada (2.96 million in 2018) and the fourth largest city in North America (630.2 km 2 ). In recent decades, Toronto has been among the fastest-growing large metropolitan areas in high-income countries and became the principal commercial center in Canada. Like most megacities in North America, Toronto formed a monocentric urban structure with large urban sprawl in the suburban area after urbanization. The downtown area is densely built with commercial and high-rise residential land use. The monocentric urban structure and longterm urban growth indicate that Toronto is in the phase of steady urbanization stage, which led to the rapid urban land use changes within the study area. Thus, the City of Toronto is a suitable study area for testing the performance of the proposed A-VCA model.
The travel data used to estimate the accessibility were acquired from the 2011 TTS. The TTS is a 1-day telephone-based travel survey that is sampled at a rate of 5% of households within the Greater Toronto Area (Data Management Group 2013). In this study, we randomly selected the auto in-vehicle travel time records from the TTS on a weekday to estimate the accessibility variables. The auto in-vehicle travel time records are archived as an origin-destination matrix of average travel hours from 625 TTS Traffic Analysis Zones (TAZs) of the City of Toronto. The accessibility of a land parcel depends not only on the traffic demands from the origin to the destination but also on the width of the road network. Since the TTS travel time records are synthesis observations from the real traffic flow in the transportation network, these traffic time data should have already implicitly taken both the traffic demand and road width into consideration. The time period we used in this study was the morning rush hour on a weekday.
The land use data used in the simulation were obtained from the Geospatial Competency Centre, City of Toronto, which originally includes the detailed land use information in 50 subcategories. The vector-format land use data set in the study area of the City of Toronto consists of totally 487,509 parcels. To simplify the initial data set, the land use data were generalized and then reclassified into 10 categories according to Table 1, including high-rise residential, low-rise residential, commercial, industrial, government, institution, or community (G/IC), open space, mixed, utilities, vacant, and water. The initial year of vector land use data we used is 2012, and the simulation was performed from 2012 to 2016. Because of the availability of land use data in 2016, we conducted the model calibration only based on land use data in the downtown Toronto, where detailed spatial information of multifarious land use changes is available.
The auxiliary geospatial variables data, which are contributed to drive the urban land uses from transforming into others, were used to estimate land use suitability with the RFA regression model. The details of the geospatial variables selected in our study are summarized in Table 2. The auxiliary spatial variables can be categorized into three major categories, including natural condition variables, traffic location variables, and urban environment variables. Figure 2 shows the visualizations of the spatial variables after using spatial analysis tools (e.g. Euclidean distance, Kernel density) via ArcGIS toolboxes (ESRI Company, CA, USA; https://www.esri.com/software/arcgis). All the data used in this study were registered and reprojected into the Universal Transverse Mercator zone 17 north coordinate system.

Methodology
In this study, we developed a A-VCA model to simulate the land use changes among different urban subtypes driven by the transportation system. Figure 3 illustrates the flow of the land use change simulation using the A-VCA model. Firstly, the selected geospatial driving factors and urban land use maps were used to calibrate a random forest regression to estimate the land use suitability (probability-ofoccurrence) for each urban land use subtype. Secondly, the transportation data collected from TTS were processed to derive the local and zonal accessibilities for vector and TAZ scales, respectively.
Thirdly, a vector-based CA model was built, taking into account the land use suitability (probability-ofoccurrence), accessibility (combination of local and zonal accessibilities with the Cobb-Douglas function), neighborhood environments, development constraints, and unexpected random factors. Finally, the land use changes of difference urban subtypes were iteratively simulated with the vectorbased simulation strategy.

Accessibility
The accessibility to urban land use types was generated from the travel time derived from the actual trip generation and conducted from the production and attraction per vector for each specific land use category. The transportation data from TTS were based on the TAZ level, and the trip generation could be simply calculated from the average number of trip origins and trip destinations per TAZ. However, the trip generation results calculated in this manner could be biased and bring the inhomogeneity since they only focus on the TAZ center rather than the global extent. In fact, the actual trip generations for the rest of the TAZ are heterogeneous and may be greatly different from the average situation because of some traffic condition such as congestions. Thus, we designed a twoperspective accessibility (i.e. local and zonal) to better describe trip generation at a finer scale: where A k;i denotes the accessibility for urban type k in urban land parcel i, A local denotes the local accessibility calculated at the vector scale, and A zonal denotes the zonal accessibility calculated at TAZ scale. We applied the Cobb-Douglas production function to numerically combine these two accessibilities together, and Equation (1) can be rewritten as where μ, a, and b are parameters set according to empirical knowledge. The Cobb-Douglas production function has long been popular for measuring and regressing factors accounting for the differences in productivity (Lobo et al. 2013;Montolio and Solé Ollé 2009) because of its simple formulation and empirical estimation mechanisms, and can be generalized in urban studies and geography (Santías, Cadarso-Suárez, and Rodríguez-Álvarez 2011). The local accessibility in Equation (2) is defined as a function to describe the accessibility for urban land parcels relative to the road network: where M k;r is calibrated parameter for differenct land use type k and road type r (we considered two road types in this study: highways and main roads), d i is the nearest distance from current land use k in land parcel i to the road network, and D r is the distance decay parameter to control the proximity effect for road type r. For current land parcels located near the highway interchanges within a specific distance, the accessibility for people to get to the highways is superior to other roads (main roads). Thus, we built a set of buffer rings around the highway interchanges and calculated the d i for these land parcels according to the buffer distance without considering other road types. For current land parcels located outside the buffer rings, the d i were calculated through finding the nearest Euclidean distance to the road network. The buffer rings distance and D r can be calibrated in calibration procedure.  The zonal accessibility in Equation (2) represents the accessibility that can be reached by all land use types from all TAZs. We assumed that the zonal accessibility for the current land use types in a TAZ is influenced by all other land uses within the study area through a distribution function. The gravity model (Erlander and Stewart 1990) was adopted to calculate the zonal accessibility of all land parcels from all other TAZs: where z and z' denote the origin and destination of TAZs, respectively, t z;z 0 denotes the in-vehicle travel time from origin to the destination zone, L k counts the current iterated land use amount for land use type k in zone z', and α and β are the parameters of the gravity model that vary from cases to cases and calibrated and set in calibration procedure.

A-VCA model
Different from the conventional grid-based CA models, the proposed A-VCA model uses land parcel as the basic modeling unit for the simulation process. Thus, the land use data should be divided into small land parcels beforehand in order to simulate the subtle urban land use changes. Based on previous studies (Vanegas et al. 2009;Yao et al. 2017), we adopted an interactive dichotomy strategy to spit large land parcels into subparts and maintain the along-the-road direction after the land parcel subdivision. For the detailed process of land parcel subdivision, please refer to Yao et al. (2017).
In general, the transition rules in a classic CA model typically consist of four components (Li and Yeh 2002;Liu et al. 2017): the land use suitability to describe the suitability for specific land use type to convert according to the essential nature of land use, the neighborhood interaction to represent the spatial autocorrelation, the constrain to denote the difficulty of change from one land use type to another one for a specific land parcel, and the randomness to account for unexpected random factors. In this study, we introduced the accessibility estimated from the previous section as an additional independent component to derive the urban land use transition rule within the proposed A-VCA model: where the subscripts k and i represent land use type k and land parcel i, respectively; P k;i denotes the final transition probability; S k;i denotes the land use suitability; N k;i is neighborhood interactions; A k;i represents the processed accessibility; the stochastic factors R s can be further expressed as 1 þ ðÀ ln γÞ α , where γ and α are uniform random variables where γ ranges from 0 to 1 and α ranges from 0 to 10, respectively; and C represents the constraints of the conversion. The land use suitability S k;i in Equation (5) was calibrated with a series of geospatial variables listed in Table 2. We adopted the RFA regression model to calibrate and estimate the land use suitability since it is capable of overcoming the multicollinearity problem among spatial variables and efficient for highdimensional fitting tasks (Palczewska et al. 2014). For each land parcel, the land use suitability can be expressed as where nTtee denotes total number of decision tree in the RFA, I � ð Þ is the indicator function, X denotes the high-dimensional vector of auxiliary spatial variables listed in Table 2, and Pre t X ð Þ denotes the prediction land use type of the tth decision tree for X.
The neighborhood interaction in the CA model is an important component to describe the spatial autocorrelation of complex geographical phenomena (Dahal and Chow 2015;Liu et al. 2017). Traditional neighborhood definitions in the grid-based CA model, such as the Moore neighborhood and Von Neumann neighborhood, are not applicable for the vector-based CA models because of the irregular geometry of land parcels. In this study, we adopted a land parcel area-weighted method (Abolhasani et al. 2016;Yao et al. 2017) to define the neighborhood effect in the A-VCA model: where Ω t i;k denotes the neighborhood effect of land use type k for land parcel i at time t; d ij is the distance between the centroids of land parcel i and j; d is the buffer range from land parcel i; S j and S i represent the geometric areas of land parcel j and i, respectively; S max and S min are the maximum and minimum geometric areas of land parcel within the entire study area. The Ω t i;k is updated in each simulation iteration and only valid when d ij < d.
The constraints of the conversion C in Equation (5) represent the easiness (or cost) of a particular land use type to change to another. Here, we set land use of water-type-restricted areas that cannot be changed during the simulation. The conversion costs among other land use types were estimated based on an analysis of the historical land use data in the study area and regional expert opinions. During the simulation, we iteratively estimated the final transition probability for each land use type in every land parcels according to Equation (5). A conversion was taken when the highest transition probability at a land parcel exceeds the conversion thresholds for different land use types. The threshold for a specific conversion (e.g. from residential to commercial) was set to be the mean probability of all transition probabilities of this specific conversion estimated from the entire study area. The simulation will keep iterating until all land use changes reach the land use areas in the terminate year.

Model assessment indicators
In this study, three indicators were selected to assess the simulation performance of the A-VCA model, including the kappa indicator, Figure of Merit (FoM) indicator, and the Cumulative Producer's Accuracy (CPA). Firstly, we used the kappa indicator, a traditional and reliable approach for evaluating the accuracy of CA land use simulations, especially for multiple land uses (Sim and Wright 2005;Van Vliet, Bregt, and Hagen-Zanker 2011). After generating a confusion matrix, the kappa indicator can evaluate the simulation results against the actual land use pattern: where p k+ represents the proportion of land parcels in the confusion matrix simulated in land use type k, p +k represents the proportion of land parcels in the confusion matrix that are land use type k in the actual land use pattern, and p kk represents the proportion of land parcels that were simulated as land use type k and are actually land use type k in the actual land use pattern.
Previous studies suggested that using the kappa indicator might be unstable in simulations in which the unchanged land parcels dominate (Pontius et al. 2008;Chen et al. 2014;Liu et al. 2017). Thus, we also adopt the FoM indicator (Pontius et al. 2008) to assess the performance of the simulation results. The FoM is often used in the CA simulation estimations to reveal the relevance of the ground truth land use changes and the simulated urban land use changes. It is a ratio of four subvariables, namely, the denominator (A), the numerator (B), and the errors C and D: Specifically, A represents the error due to observed land use change and is simulated as the persistence, B represents the number of samples that are observed and simulated correctly, C is the error caused by the observed changes and is meanwhile simulated as an incorrect gaining category, and D represents the error caused by observed persistence and is simulated as changed.
The last indicator is the producer's spatial accuracy (PA, and its cumulative form, namely, CPA). The PA indicator is commonly applied to evaluate the performance of simulating correctly and stands for the proportion of correctly simulated cells to all cells (Lu et al. 2019): where C area represents the count of correctly modeled transformed cells, and TC area is the total number of all transformed cells. Note that the PA indicator can also be calculated through four subvariables in the FoM indicator, that is, B/(A + B + C). However, the PA calculated in this manner is more reasonable for assessing the accuracy of grid-based CA simulations (He et al. 2015). In this study, the PA indicator is calculated specifically for vector-based CA by taking the area of transition land parcels into consideration. Furthermore, we designed the cumulative form of PA, namely, CPA, to assess the stability of the simulation model through multiple repeating simulations: where M is the repeating times of simulation.

Calibration of the A-VCA model
In this study, the auto in-vehicle travel time was calculated using EMME V4.2 software, and the urban land use changes was simulated using the proposed A-VCA model. Through model calibrations, we calibrated all parameters for the calculation of the gravity modelbased zonal accessibility, the distance decay effectbased local accessibility, the Cobb-Douglas function, the RFA formulation, and the neighborhood distance. In local accessibility calibration, M k;r is defined as the weight parameter for road type r for land use k, and D k;r is the impedance parameter related to distance decay effect for road type r and land use k. Table 3 lists the results of the calibrated parameters of local accessibility. All other certain calibrated parameters are listed in Table 4, where the value of β in zonal accessibility calibration refers to a previous study on Toronto (Kasraian, Maat, and van Wee 2017). Besides, the data sets for the A-VCA calibration were randomly divided into two groups, that is, training samples (60%) and validation samples (40%). Thus, the training samples used to calibrate the A-VCA model are completely independent from the validation samples used to evaluate the model performance.

Simulation of urban land use changes
The estimation of accessibility is the first step of the implementation of the proposed A-VCA model for the land use simulation in this study. As mentioned in Section 3, the accessibility used in this study consists of two separate components, that is, the zonal accessibility, generated by gravity model with actual in-vehicle travel time from the TTS records, and the local accessibility estimated based on the hierarchical distance decay effect. Figures 4 and 5 illustrate the zonal accessibility results for nine land use types and the local accessibility results toward different land use and road hierarchies, respectively. Generally, the higher zonal accessibility for urban land types of residential, commercial, G/IC, mixed, and utilities is mainly located in the central area of Toronto. The higher zonal accessibility for industrial land use is found in the west, whereas the higher zonal accessibility for open space and vacant is in the east. For local accessibility, the higher values are found mainly around the highways (Ontario 401 Express, Don Valley Express, and Gardiner Express) since they have higher priority in the estimation of local accessibility.
We used the training data set that contains the urban land use changes and auxiliary geospatial variables to calibrate the RFA regression model and estimate the land use suitability for each urban land use type. Figure 6 visualizes the spatial patterns of the land use suitability for nine urban land use types after a minimum-maximum normalization. They describe the suitability for specific land use type to covert according to the essential nature of land use. This pattern of land use suitability will be combined with the estimated accessibility, the neighborhood effect, the constraints of conversion, and stochastic factor to update the final conversion probability in every iteration of the simulation process according to Equation (5).
The proposed A-VCA model was applied to simulate the land use changes among different urban land use types from 2012 to 2016. The simulation will keep where T is the total iteration time, K is the number of urban functional types, N is the total land use parcels in the study area, and J is the number of land use parcels within the neighborhood buffer range. The proposed A-VCA model was implemented in C++ programming language. We carried out the simulations on a local workstation equipped with Intel Core i9-9900 K CPU, NVIDIA GeForce RTX 2080 (8 GB graphics memory), and 16 GB dual-channel DDR4 RAM. The computational time cost for the entire simulation process was 37.32 ± 3.42 s in 10 times repeating simulations.
We presented the actual (Figure 7(a)) and simulated (Figure 7(b)) patterns in 2016 to show a visual comparison at global scale. Overall, the simulated pattern agrees well with the actual land use across the entire study area. The zoom-in view of the simulation result in downtown Toronto is also illustrated in Figure 8 to present more spatial details of land use changes and the spatial agreement of the simulation result. From the zoom-in view of the downtown Toronto (Figure 8(a,b)), we can see that obvious land use changes were taking place during this period, especially in areas marked with dashed circles. Favorably, most of the changed land parcels were well simulated with the proposed A-VCA model (Figure 8(c)).
The regions marked with dashed circles (Figure 8) are three typical areas where urban land use changes were taking place during this period. Region 1 is an aggregated commercial zone near the Yonge Street. The main change of land use in this area was the transition from vacant to commercial land use. In region 2, many vacant land parcels were transformed into G/IC buildings during this period because this area possesses very high accessibility to multilevel roads and various land uses. Region 3 is the south cabbage town of Toronto with many low-rise residential buildings aggregated (e.g. semi-detached house, town house). Some vacant land parcels were developed into functional buildings like G/IC and mixed uses (residential and commercial). Our model generated a good result to show the details of land use changes in these three regions.

Assessments of the urban land use simulation
To quantitatively evaluate the performance of the simulation, we first summarized the proportions of actual and simulated areas of land use changes for each land use type and examined their agreement (Table 5). We found that the average bias of all land use types is +0.56% for the entire City of Toronto, and the largest bias of the simulation is the vacant type (−2.52%). As to the downtown Toronto, relatively higher average bias was observed for all land use types (+1.37%), and the largest bias was found in the commercial land type (−3.84%). This might be due to the fact that there are much more fragmented land parcels in the downtown Toronto than in the entire city. Nevertheless, the small biases of the simulated land use pattern in Table 5 suggest that the proposed A-VCA model is capable of simulating the changes of multiple urban land use types at least from the quantity perspective.
In order to quantify spatial agreement of the simulation against actual land use pattern, we constructed the confusion matrix of the simulation result and estimated three indicators described in Section 3, that is, the kappa, FoM, and CPA indicators. The results of these three indicators are summarized in Table 6. From the assessment results shown in the table, the proposed A-VCA model yields satisfactory pattern agreement when applied in the City of Toronto with both values of kappa and overall accuracy higher than 0.90. Moreover, the FoM indicator, which focuses on the change simulation and can avoid bias introduced by persistent unchanged land parcels, yields a value of 0.283 for the study area. According to the recent studies on traditional urban development simulations (from nonurban to urban), a range of FoM values from 0.10 to 0.30 have been reported at various spatial scales (Chen et al. 2014;Thapa and Murayama 2011;Xu, Zhang, and Chen 2020). Thus, our simulation result achieves similar and even slightly better than the results reported previously, considering the simulation was conducted among different subtypes of urban land use changes. Moreover, we conducted multiple repeating simulations to estimate the CPA and its standard deviation that take the transition area into consideration. A CPA value of 72.83% ± 1.535% was derived after 10 times repeating simulations, indicating that the proposed A-VCA model is capable of describing the urban land use change with satisfactory accuracy and stability.
To demonstrate the outperformance of our proposed A-VCA model, we compared the A-VCA model with three existing models that are capable of simulating multitype land use changes. These comparative models include the original VCA model Figure 5. The normalized local accessibility in the study area. ) without the accessibility component and two grid-based CA models, that is, the ANN-CA model (Li and Yeh 2002) and the future land use simulation (FLUS) model . The purpose of the comparison with the VCA model is to test the effectiveness of the accessibility-interactive mechanism in the vector-based simulation strategy; and the comparisons with the two grid-based CA models are to examine the superiority of the vectorbased strategy in simulating fragmented land use  changes at a fine scale. All these comparative models were calibrated with the same geospatial variables listed in Table 2, and the transition probabilities were estimated without the accessibility variables. As to the grid-based ANN-CA and FLUS models, the urban land use change was calibrated and simulated at spatial resolution of 30 m, which is a commonly used resolution in many previous gridbased CA studies. The comparison results are summarized in Table 6. We can observe from the table that the simulated result from the VCA model is apparently not as good as the result simulated from the A-VCA model. All three assessment indicators show decline in varying degree (the second row in Table 6) compared with those of the A-VCA model. This confirms the importanct role of accessibility in driving the urban land use dynamics, and the cooperation of accessibility in the CA model can better simulate the subtle land use changes in urban area. The comparisons against the two grid-based CA models also show lower kappa, FoM, and cumulative PA values of the ANN-CA and FLUS models (the last two rows in Table 6). Note that the cumulative PA values for these two grid-based CA models were estimated based on grid unit and thus may not be comparable to the value of the A-VCA model. It is worth mentioning that the assessment indicators for the two grid-based CA models are even lower than that of the VCA model, which suggests the effectiveness of the vector-based strategy in the simulation of fractal land use changes among multiple urban subtypes at a finer spatial scale.  The size of land parcel is the key variable for determining the details of urban land use from both perspectives of spatial detail and urban functional types. Thus, the sensitivity of the simulation performance to the size of land parcel needs to be examined. In order to conduct the simulation on different sizes of land parcel, we manually broke down each of the original land parcel into two (1/2), three (1/3), and four (1/4) subparcels, and maintained the land use types of the subparcel consistent to the original ones. Within each original land parcel, the areas of the broken subparcel were kept as equal as possible. Then we conducted the simulations on the broken subparcels to test the sensitivity of the model performance to the land parcel size. The quantitative comparison among different land parcel sizes is summarized in Table 7. From the table, we can observe that all the assessment indicators show a slight decrease as the size of parcel becomes smaller. Such a minor performance degradation might probably be due to the reason that the increased number of land parcels introduces much more complex conflicts among different urban functional types. As expected, the computational cost increased dramatically (from 37.32 to 202.57 s) as the size of land parcel became smaller since the total number of land parcels was greatly multiplied in each iteration and the neighborhood effect involved much more land parcels within the buffer range.

Morphological pattern of Toronto urban land use changes
The visualization of simulation results shows that most of the urban land use changes were mainly taking place in downtown area during the study period. Most of the high-rise residential land uses are located in the downtown area, which possesses the highest land price in Toronto and is intensively built. The extent of high-rise residential land parcels was spreading because of higher housing demands in the urban center. The newly converted high-rise residential land parcels are found around the existing ones and scattered alongside main road networks, which is mainly driven by the convenient transportation facilities and high accessibility around this area. Both Figures 4 and 5 give evidence that the downtown Toronto has superior local and zonal accessibilities toward most of the land uses and better access to various facilities, motivating tenants to live in downtown and inspiring investors and government to develop more high-rise residential land uses in this area. On the contrary, most of the transitions into low-rise residential land use were simulated outside downtown Toronto, most of which were formed with old communities and rowed by townhouses or semidetached houses.
The G/IC and open space land uses were simulated closed with residential (both high-rise and row-rise) land parcels since the deployments of public facility have to take the range of geographic service areas into consideration (Donnelly 2013). This can be partially supported by the zonal accessibility maps in Figure 4 that high zonal accessibility of G/IC, open space, and residential land uses share a large overlapping area. The transition to mixed land use (commercial and high-rise residential) was mainly simulated to be near to the main road in the downtown area (e.g. Yonge Street of Toronto, Figure 8) because of the demand for both commercial and residential land use types in the densely built downtown area. The industrial land use was simulated around the original industrial estates along expressways (404 and 401) and in the southeast lake shore area of Toronto near to the harbor where many traditional factories are located. Both of those industrial estates have superior accessibility for promoting productive elements and reducing transportation costs.

The integration of accessibility and vector-based land use simulation
As previous studies suggest, the transport system is one of the crucial drivers for urban growth and urban land changes by providing accessibility and economic opportunity to the surrounding lands and social activities (Iacono, Levinson, and El-Geneidy 2008;Wee 2004). Accessibility is the typical indicator for describing travel behavior in urban area (Luna, J, and Shoshanna 2018). Thus, in this study, we introduced the accessibility derived from actual travel time data rather than simplified Euclidean distances to simulate fine-scale urban land use changes driven by transportation facilities. The accessibility at both scales, that is, the zonal accessibility calculated at TAZ scale and the local describing proximity of land parcel to the multilevel road network, was considered in this study and was numerically combined into one through the Cobb-Douglas function. In most previous CA simulation studies (Chen et al. 2014;Yao et al. 2017), the accessibility described as simplified proximity usually has a relatively lower priority and is only applied for estimating the basic land use suitability. In this study, the coupling mechanism of accessibility and CA simulation considers the accessibility as a coequal component with the land use suitability for the estimation of the final transition probability. The comparison results listed in Table 6 clearly show the outperformance of this proposed coupling mechanism in the simulation of subtle urban land use changes among different functional types. Note that, apart from the travel time, the information of road width is supposed to be another key variable of the transportation system that contributes to driving the urban land use changes. Incorporation of the road width variable into the A-VCA model is potentially helpful for better simulation of urban land use changes. This hypothesis will be tested in future studies as long as accuracy measurements of road width are available.
The vector-based land use simulation strategy was adopted to establish the A-VCA model and applied for the simulation of urban land use changes in the City of Toronto. Compared with traditional grid-based representation, previous studies suggested that the vectorbased strategy was more suitable for fine-scale land use simulations (Abolhasani et al. 2016;Dahal and Chow 2015;Barreira-González, Gómez-Delgado, and Aguilera-Benavente 2015;Chen et al. 2014;Lu, Cao, and Zhang 2015). When introducing the accessibility to CA simulation, the fine-scale shape of land parcel is more suitable for describing subtle land use changes driven by accessibility (e.g. actual travel cost in the transportation network). The vector format data is capable of capturing the changes of fine-scale land parcels with irregular polygons and avoid the bias caused by regular grid cells Chen et al. 2019;Barreira-González, Gómez-Delgado, and Aguilera-Benavente 2015). The application of the proposed A-VCA model in urban land use simulation of the City of Toronto shows satisfied accuracies and good morphological patterns. The superiority of vectorbased strategy in such fine-scale simulation was also confirmed in this study (Table 6) by the comparison of the proposed A-VCA and the popular grid-based FLUS model.

Conclusions
This study proposed an accessibility-interacted vectorbased land use simulation model, namely, A-VCA model, to investigate the urban land use changes driven by transportation facilities. The accessibility at both local and zonal scales derived from auto-in-vehicle travel time data was considered as a key driver of fine-scale urban land use changes among different urban subtypes. The proposed A-VCA model was applied and tested in the City of Toronto, Canada, to simulate the fine-scale urban land use changes during 2012-2016. The results show that our model is capable of accurately simulating the subtle urban land use changes driven by the accessibility with good morphological features. The satisfactory performance of the A-VCA model was testified by the quantitative assessment with both values of kappa and overall accuracy higher than 0.90 and FoM larger than 0.28. The effectiveness of the accessibility-interactive mechanism was confirmed by the comparison with the original VCA model. Besides, the superiority of the vector-based strategy in the simulation of subtle urban land use change among different functional types has been demonstrated through the comparisons with the two popular grid-based CA models.
However, due to the unavailability of detailed urban functional land use data, we only calibrated and evaluated the proposed model in the same period (2012)(2013)(2014)(2015)(2016). In fact, calibrating the model in one period, and then applying the calibrated model for the simulation and evaluating it in another period can better justify the performance of the proposed model. Besides, whether or not the proposed A-VCA model can be generalized and directly applied for other cities has not been tested and still needs further investigations. Nevertheless, the proposed model provides new tools for fine-scale simulation of subtle land use change among different functional types in urban area, which are rather useful for many urban-related studies such as job-housing matching and urban facility layout optimization. Additionally, the results and findings can deepen our understanding of the interaction between land use changes and transportation, and provide knowledge for urban growth management and transportation planning.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
with the School of Geography and Planning, Sun Yat-Sen University. His research interests include image processing, artificial intelligence, and geographical simulation.
Jinpei Ou received his PhD degree in geographic information science from Sun Yat-Sen University, Guangzhou, China, in 2016. He is currently an associate professor with the School of Geography and Planning, Sun Yat-Sen University. His research interests are spatial simulation and optimization, remote sensing processing, and carbon emission.
Xinxin Wu received her MS degree in geographic information science from Sun Yat-Sen University, Guangzhou, China, in June 2020. She is currently a PhD candidate in cartography and geography information system of Sun Yat-Sen University, Guangzhou, China. Her research interests include multi-source spatial data mining and urban morphology understanding.