Validation of CA-Markov for Simulation of Land Use and Cover Change in the Langat Basin, Malaysia

Validity of CA-Markov in land use and cover change simulation was investigated at the Langat Basin, Selangor, Malaysia. CA-Markov validation was performed using validation metrics, allocation disagreement, quantity disagreement, and figure of merit in a three-dimensional space. The figure of merit, quantity error, and allocation error for total landscape simulation using the 1990-1997 calibration data were 5.62%, 3.53%, and 6.13%, respectively. CA-Markov showed a poor performance for land use and cover change simulation due to uncertainties in the source data, the model, and future land use and cover change processes in the study area.


Introduction
Globally, land use change impacts biodiversity, water and radiation budgets, emission of green house gases, carbon cycling, and livelihoods.The study of Land Use and Cover Change (LUCC) and its dynamics is crucial for environmental management, especially with regard to sustainable agriculture and forestry.Different models, in terms of structure and application, have been used to understand LUCC dynamics [1,2].
In this study, CA-Markov was chosen to simulate land use change based on the following requirements: dynamic simulation capability, high efficiency with data scarcity, simple calibration, ability to simulate multiple land covers and complex patterns.
CA-Markov has the ability to simulate land use changes among multiple categories and combines the CA 1 and Markov chain procedures [6].Markov analysis does not account for the causes of land use change and it is insensitive to space.However, CA-Markov using the CA ap-proach relaxes strict assumptions associated with the Markov approach and explicitly considers both spatial and temporal changes [7].CA-Markov also enables a more comprehensive simulation as compared to other LUCC models such as GEOMOD and CLUE [8].However, CA-Markov calibration is operationally based on a single period of time, which renders difficulty in simulating land cover dynamics on a temporal scale [9].In addition, the assumptions underlying CA-Markov method tend to be somewhat simplistic when looking at the micro-level changes in land use [10].Pontius and Malanson (2005) [11] have demonstrated that the predictive power of CA-Markov is higher for cases where it concentrates on the major signal of land changes and ignores noises.Pontius and Malanson (2005) [11] compared CA-Markov and GEOMOD in terms of simulation power and suitability for different applications in Central Massachusetts, USA.They applied a three stage method to measure simulation power.At first, calibration process was separated from validation process.Secondly, the accuracy was assessed at multiple resolutions.Finally, the calibrated model was compared to a null model that simulates pure persistence.Their investigation showed that the added complexity of spatial contiguity rule in CA-Markov was of no benefit.Paegelow and Camacho (2005) [9] studied the potential and limitation of prospective GIS 2 -based LUCC modeling.Their approach comprised the Markov chain for temporal simulation, MCE 3 , MOLA 4 , and CA to perform spatial contiguity on the simulated land cover scores.Results showed three distinct limitations; the first was caused by complex variability within the land cover categories, the second was caused from using only two land cover maps for calibration, and the third was caused by the fact that MCE, MOLA, and CA only affected spatial distribution, and not temporal distribution, of simulated scores.In another study, Samat (2009) [12] applied CA-Markov to evaluate urban spatial growth in Seberang Perai, Malaysia.Validation analysis was performed using Kappa Agreement Index (KIA).In this work, the poor performance of CA-Markov in the development modeling of commercial/public facilities and industrial activities mainly resulted from: Model development based on physical factors, model assumption based on the uniform spatio-temporal growth of urbanized area, and inability of the model in recognizing new development.Samat et al. (2011) [13] showed that the prediction accuracy of CA-Markov decreased when the model tried to predict for a longer period of time, possibly due to the fact that a uniform transition rule was used by the model throughout the simulation period.According to Jokar Arsanjani et al. (2011) [14], the critical issue in employing CA-Markov is to combine the social, human and economic dynamics in the simulation, which can be realized in agent-based modeling systems.However, they tried to solve this problem using the integration of logistic regression and CA-Markov.
Araya and Cabral (2010) [15] modeled and analyzed urban land use change using CA-Markov and landscape metrics in Portugal.They reported that CA-Markov was validated successfully, with Klocation and Kquantity of 87% and 83%, respectively.However, they did not analyze the model capability for change simulation using a three-dimensional approach.Wang et al. (2010) [16] studied land use change based on vector data source using CA-Markov in China.Although the vector data improved simulation accuracy, validation analysis was performed only based on Kappa values.
Recent validation techniques proposed by Pontius and Millones (2011) [17] and Pontius et al. (2011) [18] allow the understanding of complex relationships caused by multiple simulation parameters.Additionally, these techniques facilitate three-dimensional spatial and temporal comparisons of LUCC, and eliminate the weakness of Kappa values for validation analysis.
This work was aimed at evaluating the capability of CA-Markov in simulating LUCC in a tropical catchment using the validation techniques proposed by Pontius and Millones (2011) [17] and Pontius et al. (2011) [18].

Study Area
In recent decades, the Langat Basin has undergone rapid urbanization, industrialization and agricultural development [19,20].The Langat Basin is also a main source of drinking water for surrounding areas, a source of hydropower and plays an important role in flood mitigation.Over the past four decades, the Langat Basin has served approximately 50% of the Selangor State population.However, the Selangor State is currently facing water shortage problems, especially in urban areas [21][22][23].Hydrometeorologically, the Langat basin is affected by two types of monsoon, i.e. the northeast (November-March) and the southwest (May-September) monsoons.Average annual rainfall is about 2400 mm.The wettest months are April and November with average monthly rainfall exceeding 250 mm, while the driest month is June with average monthly rainfall not exceeding 100 mm.Topographically, the Langat basin can be divided into three distinct areas with reference to the Langat River: the mountainous area in the upstream, undulating land in the centre and flat flood plain in the downstream (Figure 1).The basin has a rich diversity of landforms, surface features and land cover [23,24].Dominant soil types in the basin are steepland and Rengam-Jerangau soil series with sandy clay loam and clay textures, respectively [20].
Due to the national importance of the Langat Basin, three sub basins (upstream of the Langat River), i.e.Lui, Hulu Langat, and Semenyih, with a total area of 694.13 km 2 were selected for investigation of land use change (Figure 1).

Data Set
Land use maps dated 1990, 1997 and 2002 and soil maps, with the scale of 1:50,000 were collected from the Soil Resource Management and Conservation Division, Department of Agriculture, Malaysia.Ancillary data such as road and stream networks, contour lines, and population centers were extracted from the topographic maps, with the scale of 1:50,000.Landslide points, geology map (scale: 1/100,000), and location of the water treatment plants were obtained through the library of Universiti Putra Malaysia.Sawmill and oil refinery location maps were extracted from the agriculture yellow pages of Malaysia (www.yellowpages.com.my) and Google Earth.Daily TRMM 5 radar data during 1998-2009 were obtained through the TOVAS and utilized to produce rainfall map.TRMM Online Visualization and Analysis System (TOVAS), developed by the NASA, provides a user friendly web-based interface for visualization and analysis of TRMM gridded rainfall products and other precipitation data [25].In this work, land use classification and map resolution (i.e. 30 m) were outlined based on the objective of a broad study, which was aimed at hydrological analysis of the Langat Basin during 1984-2020.The optimal pixel size was calculated based on the complexity of terrain method [26].

CA-Markov
A random process, X(t), is a Markov process if for any moment in time, [27].The random process satisfies (see Equation ( 1) below) where: t n is the present time so that t n + 1 represents some points in the future and t 1 , t 2 , ••• , t n − 1 represent various points in the past.Based on the present data, future is independent of the past.In other words, future of random process depends neither only on where it is now nor on how it got there [28].
If X[k] is a Markov chain with the states {x 1 , x 2 , x 3 , •••}, then probability of transition from the state i to the state j in one time instant is, If Markov chain has a finite number of states, i.e. n, transition probability matrix can be defined as follows: 1,1 If transition probabilities vary with the time, the above matrix need to be written explicitly as a function of k (e.g.p i, j, k ) [28].
CA-Markov in IDRISI program involves two techniques, i.e.Markov chain analysis and cellular automata [15].Cellular automata underlie dynamics of the change events based on proximity concept so that the regions closer to existing areas of the same class are more probable to change to a different class.A cellular automaton is a cellular entity that independently varies its condition based on its previous state (according to a Markov transition rule) and adjacent neighbors [6].The following 5 × 5 contiguity filter was utilized in this work:

Modeling
Firstly, three calibration periods, i.e. 1990-2002, 1997-2002 and 1990-1997 were considered and correspondent simulation results were pre-analyzed in terms of goodness of fit.Results showed that the period 1990-1997 was more capable to simulate the future changes of land.Therefore, in this study, calibration data over the period 1990-1997 was used to extract transition probability matrix.Then two types of criteria (constraints and factors) were developed to determine which lands are suitable for future development.The constraint factors (dam reservoir and stream) were standardized into a Boolean value (0 and 1), while the parameter maps were standardized into a suitability continuous scale from 0 (minimum suitability) to 255 (maximum suitability) (Figures 2 and  3).Three types of fuzzy membership functions, i.e. sigmoidal, symmetrical, and user defined were used to rescale parameter maps into the range 0 -255.
Analytic Hierarchy Process (AHP) was used to determine the weights of driving factors.AHP is a measurement theory founded on expert judgment to drive priority scales using pair wise comparisons [29].In this method, comparisons are made based on a scale of absolute judgment that shows how much more one element prevails over the other for a given attribute (Table 1).However, the judgment may be inconsistent [30].Considering n elements to be compared (C  for all i, j and k.In this condition, Aω = λω, ω is an eigenvector (of order n) and λ is an eigenvalue.For a consistent matrix, λ = n.Human judgments are inconsistent to a greater or lesser degree, therefore in matrices concerning human judgment, the condition ij jk ik   will not be reliable.In such a case, Aω = λ max × ω and λ max ≥ n.The difference between λ max and n is an indicator for inconsistency of the judgment.Consistency Index (CI) can be measured by  where: RI is the Random Consistency Index which is a function of number of elements.If the CR exceeds 0.1, set of judgments may be too inconsistent to be reliable and the CR of zero means that the judgments are perfectly consistent [31].SINMAP6 approach [32] under ArcView GIS was  applied to derive stability and saturation indices.Stability index was used in the extraction of suitability maps for agriculture, bare land, mining land, and urban areas, while saturation and roughness indices were used in suitability mapping for marshland and water bodies.Runoff accumulation map was produced based on impervious areas and Digital Elevation Model (DEM).This map was utilized to determine suitable areas for marshland and water bodies.Sediment transport capacity index (topography factor in USLE 7 model) was used to extract suitability maps for grassland, marshland, and agriculture (Figure 2).Distance to sawmill and oil refinery centers were important parameters in the extraction of suitability maps for forest and oil palm stands, respectively.Distance to Water Treatment Plant (WTP) was used as a driving factor in AHP and MCE processes due to its impact on development of agriculture and urbanization (Figure 2).
tion, which was calibrated using the data of 1990-1997.Disagreement parameters determine the disagreement between simulated map and observed map [17,18,33].Quantification error (quantity disagreement) happens when the quantity of cells of a category in the simulated map is different from the quantity of cells of the same category in the reference map.Location error (allocation disagreement) occurs where location of a class in the simulated map is different from location of that class in the reference map [17].These metrics were calculated based on the methods presented by Pontius and Millones (2011) [17] as follows.

Disagreement Components
These calculations and indices were derived directly from Table 2. J refers to the number of categories and number of strata in a typical stratified sampling design.Each category in the comparison map is indexed by i, which ranges from 1 to J. The number of pixels in each stratum is denoted by N i .Each observation is recorded based on its category in the comparison map (i) and the reference map (j).The number of these observations is summed as the entry n ij in row i and column j of the contingency matrix.Proportion of the study area (P ij ) that is category i in the simulated map and category j in the observed map is estimated by the following equation [17]: As illustrated in Figure 2, Weighted Linear Combination (WLC) method was used in multi criteria evaluation of the drivers and constraints.WLC multiplies each standardized factor map by its factor weight and then sums the results.This result is then multiplied by each constraint to mask out unsuitable areas [6].

Model Validation
Land use map dated 2002 was used in the model valida-


Reference Total The Equation ( 6) calculates quantity disagreement (q g ) for an arbitrary category g. 1 1 The overall quantity disagreement (Q) which incorporates all J categories is computed by Equation ( 7): Allocation disagreement (a g ) for an arbitrary category g is calculated using Equation (8).The first argument within minimum function is the omission of category g while the second argument is the commission of category g. 2min , The overall allocation disagreement (A) is computed by the following equation: The proportion of agreement C is estimated by Equation (10): Based on Equation ( 11), total disagreement D is a sum of overall quantity of disagreement and overall allocation of disagreement:

Figure of Merit
The figure of merit equals to the intersection of observed change and simulated change divided by the union of observed change and simulated change.The figure of merit ranges from 0, i.e. no overlap between real and predicted changes, to 100%, i.e. perfect overlap between real and predicted changes.

Figure of Merit
where: A is the area of error due to observed change simulated as persistence, B is the area of correct due to observed change simulated as change, C is the area of error due to observed change simulated as change to wrong category, and D is the area of error due to observed persistence simulated as change [33].

Land Use Change
Calibration of CA-Markov was done based on the changes in land cover over the period 1990-1997 (Figure 4).
According to transition probability matrix (Table 3), the biggest change occurred on grassland.Field and desk investigations demonstrated that grassland with bare land and mining land are transit categories and hence exposed to more changes during the time.Agriculture land 8 , oil palm and rubber stands also changed over the time.Agriculture land shifted mostly to oil palm and rubber plantings and urban areas while older rubber stands were transformed to agricultural and urban areas.Bare land changed mostly to grassland, forest and oil palm plantings.Forest stand primarily changed to grassland, mining land, urban areas, and rubber planting.Grassland and marshland mainly shifted to rubber planting, agricultural and urban areas.Mining land and oil palm stand changed  primarily to forest and rubber plantings, bare land, and urban areas.In this region, water bodies did not change much.However, there were some transitions of water bodies to marshland and urban areas.However, some of these changes in land use as indicated in the land use maps appear as a mismatch in reality.For example, the transitions of bare land and grassland to forest over a period of 7 years seem unreasonable.Also, the transition of urban areas to agricultural and rubber plantings seem unlikely under the Malaysian policy of land development at the time [34,35].This mismatch is not caused by the model, but rather it could have been caused by the use of different land use classification criteria or inconsistent resolution of satellite imagery over the 12-year period.
As shown in Table 4, twenty factors were entered into AHP to determine the weights proportional to their importance in land allocation process.Consistency ratio in all land cover categories is lower than 0.1, which infers a consistent and reliable judgment.Figure 3 depicts the transition suitability of different land uses within the Langat Basin.

Validation
Based on error analysis (Tables 5 and 6), change simulation of bare land and grassland categories resulted in the lowest figure of merit, i.e. 47.39% and 49.45%, respectively.This could be due to the transit nature of these categories.Conversely, simulation of forest and water bodies resulted in the highest figure of merit, i.e. 94.39% and 80.32%, respectively.This could be caused by high percentage of persistence in these categories, i.e. 97% and 99% for forest and water bodies, respectively.Due to protection policy of the Hulu Langat Forest Reservoir, transition rate in the forest area was very low, which caused a high persistence.Figure of merit in change simulation of marshland and agriculture land was 51.03% and 57.42%, respectively.In simulation of agriculture land, quantity error was nearly similar to allocation error.This indicates that simulation of agriculture land is impacted by both types of errors.However, in simulation of marshland, quantity error was significantly higher than the allocation error.In this work, simulations of marshland and agriculture land were mostly impacted by transition inconsistency between calibration and validation periods.Figure of merit in change simulation of mining land and oil palm was 60.85% and 66.01%, respectively.In simulations of mining land and oil palm stands, transition consistency between calibration and validation periods was relatively high, which resulted in lower quantity error than allocation error.Figure of merit in change simulation of urban area and rubber stands was 76.37% and 67.71%, respectively.Change simulation of urban area was more affected by allocation error rather than quantity error.However, in simulation of rubber stands, due to transition inconsistency between calibration and validation periods, quantity error reached 2.39%.
Table 5 shows total quantity error, total allocation error, and total figure of merit calculated using a validation analysis over three time periods which are: reference T1, reference T2, and simulated T2.Using the same approach, components of agreement and disagreement as a fraction of landscape area are given in Table 6.
Total figure of merit, total allocation error, and total quantity error were 5.62%, 6.13%, and 3.53%, respectively.According to Table 6 and Figures 5 and 6, the components "persistence simulated as change" and "change simulated as persistence" were 5.63% and 3.46%, respectively.The model was able to correctly simulate only 0.57% of the changes.Simulated change is the sum of "persistence simulated as change", "change simulated as change to wrong category" and "change simulated correctly."Observed change is the sum of "change simulated as change to wrong category", "change  When the main signal of landscape is persistence, it is important that the model focus on the most important signal of change in the landscape.In this region, there were some noises resulting from the transitions of small patches specifically in agriculture and urban categories.These noises were projected by CA-Markov.However, studies have shown that the predictive power of CA-Markov is higher for cases where it focuses on major signal and ignores noises [11].
Land use maps utilized in this work were generated from satellite imagery, topographic maps, field surveys, and manual digitization.The digitized land use maps clearly showed some mismatch among the land use polygons, especially that of grassland, bare land and ur-ban area, due to different interpretation of category definition.This led to some errors in the digitization process and consequently affected the simulation accuracy.
In this study, as confirmed by Pontius and Neeti (2010) [36], there were three sources of uncertainty, viz. the data, the model, and future land change processes.Calibration and validation processes in this study required land use maps in three points of time.The duration of time intervals was dependent on data availability.As mentioned previously, there were some inconsistencies between the calibration and validation periods in terms of change intensity that was mostly due to data scarcity and lack of freedom in data selection.This study also showed that CA-Markov was unable to simulate accurately based on marginal changes in land use over the given time period.CA-Markov simulations using the 1990-1997 calibration data set had an overall Kappa of 89% and a standard Kappa of 85%.These Kappa agreement statistics infer high capability of the model in land use change simulation based on the given calibration period.However, the results of three-dimensional validation approach as proposed by Pontius et al. (2011) [18] demonstrated that these Kappa accuracies are mostly caused by high land persistence over time, as shown in this study (Figure 6 and Table 6).Therefore, the Kappa agreement statistics are unable to measure the strength of CA-Markov in simulating land use change.Typically, Kappa statistics compare accuracy to a randomness baseline.According to Pontius and Millones (2011) [17], however, randomness is not a logical option for mapping.In addition, several Kappa indices suffer from basic theoretical errors.Thus, it is recommended that quantity disagreement and allocation disagreement be used for accuracy assessment, instead of Kappa statistics.Each LUCC model has its own capabilities and limitations.A single model is not capable to capture all of key processes in land use change projection [5].In this study, the apparent differences in land use change between urban and non-urban areas were not well projected by CA-Markov.Therefore, according to the results of this study and other recent works [37,38], and with regard to data scarcity and the limitation of the applied resolution in this work, the following refinements are suggested as improvement to the simulation accuracy by CA-Markov: 1) Modeling urban and non-urban areas separately; 2) Integration of logistic regression and CA-Markov to prepare more realistic suitability maps; 3) Incorporation of more deterministic socio-economic parameters in land change simulation, such as population growth.

Conclusions
In this work, the validity of CA-Markov was studied using the 1990-1997 calibration period.
CA-Markov did not accurately simulate land use change dynamics in this area which resulted in the high values of quantity and allocation disagreements and low value of figure of merit computed using a three-dimensional approach.
Results showed that the new patches of grassland and bare land did not expand from existing patches.CA-Markov imposed spatial dependency by the contiguity rule, which caused inaccurate simulation of land use change for bare land and grassland.Moreover, CA-Markov simulation included noises which were dominant in agriculture and urban area categories.However, the projecting power of CA-Markov was higher in situations where major signals dominated over noises.Mismatching among land cover polygons due to different interpretations of land use definitions contributed to inaccuracies in the simulation of land use change.Data scarcity and lack of freedom in data selection caused some transition inconsistencies between calibration and validation intervals.In addition, CA-Markov exhibited difficulty in change simulation due to marginal changes in this area over time.
This study demonstrates the utility of three-dimensional analysis and disagreement components in validating land use change models.

Figure 1 .
Figure 1.Geographic locations of the Lui, Hulu Langat and Semenyih sub basins within the State of Selangor, Peninsular Malaysia.

Figure 2 .
Figure 2. Computational framework of the study.

Figure 3 .
Figure 3. Transition suitability maps for different land covers generated using MCE.

Figure 4 .
Figure 4. Observed land use map versus simulated land use map.8 Refers to cultivation of annual field crops and horticultural crops, and does not include perennial industrial crops, i.e. oil palm and rubber.

Figure 5 .
Figure 5. Components of agreement and disagreement.

Figure 6 .
Figure 6.Map representation of agreement and disagreement components. 1 ••• C n ), a ijdenotes the relative weight or priority of C i over C j .A = (a ij ) is a square matrix of order n with the constraints, 1