GEO-INFORMATION BASED SPATIO-TEMPORAL MODELING OF URBAN LAND USE AND LAND COVER CHANGE IN BUTWAL MUNICIPALITY , NEPAL

Unscientific utilization of land use and land cover due to rapid growth of urban population deteriorates urban condition. Urban growth, land use change and future urban land demand are key concerns of urban planners. This paper is aimed to model urban land use change essential for sustainable urban development. GI science technology was employed to study the urban change dynamics using Markov Chain and CA-Markov and predicted the magnitude and spatial pattern. It was performed using the probability transition matrix from the Markov chain process, the suitability map of each land use/cover types and the contiguity filter. Suitability maps were generated from the MCE process where weight was derived from the pair wise comparison in the AHP process considering slope, land capability, distance to road, and settlement and water bodies as criterion of factor maps. Thematic land use land cover types of 1999, 2006, and 2013 of Landsat sensors were classified using MLC algorithm. The spatial extent increase from 1999 to 2013 in built up , bush and forest was observed to be 48.30 percent,79.48 percent and 7.79 percent, respectively, while decrease in agriculture and water bodies were 30.26 percent and 28.22 percent. The predicted urban LULC for 2020 and 2027 would provide useful inputs to the decision makers. Built up and bush expansion are explored as the main driving force for loss of agriculture and river areas and has the potential to continue in future also. The abandoned area of river bed has been converted to builtup areas.


INTRODUCTION
The urban sector has been given high priority by Nepalese policy makers since efforts for planned development began.This strategy is essential because urban sector contributes 60 percent of GDP and 4.5 times higher of average income than rural income (Karki 2004) and beside this, cities have been focal points of economic activity and engines of economic growth; centers of excellence for education, health care, innovation, entrepreneurship, business, commerce, industry, culture and social services; large markets for a wide range of products, goods and services; nodal points for transportation, telecommunications and information technology systems; and the primary location for employment and livelihood opportunities ( UN 2002 ).Urban sector development needs in Nepal to span the full range of activities, infrastructure, services, and facilities required to improve the efficiency of cities as places in which to invest and work, and to improve the condition of cities as places in which to live and prosper.
The government of Nepal has recognized that targets of poverty reduction and achievement of the MDGs have not been kept pace with increasing population growth.To meet this short fall and to fulfill rising expectation about the standard of living, urban growth and its sustainable management must be addressed.Sustainable urban growth can be achieved by providing urban facilities to growing future urban population particularly future urban land demanded.
With this in mind, the government has declared more urban centers or municipalities and their planning authorities as Town Development Committee, 1988 (TDC) and given mandates for planning and physical development of cities to Municipalities by the Local Self Government Act, 1999and Rural Urban Partnership Program, 2001(ADB 2006).
The urban economy in Nepal is growing at faster rate than the rural economy or the overall economy.As a result, predominantly rural activities are shrinking while predominantly urban activities like industry and services are growing more rapidly.The urban economy is estimated to have grown twice the rural economy per annum (ADB, 2000).Though urban population in Nepal has increased from 9.2 percent in 1991 to 17 percent in 2011 (CBS 2011), the current urban population is very low even with reference to the south Asian level (Haub 2007).
In Nepal, the major problems of sustainable urban planning and management are the lack of accurate and up-to-date information about urban expansion and changes in urban areas and its magnitude and pathways and future urban land demand for ever-growing urban population.This paper aims to model the urban land use change and its future spatial prediction using Markov and CA-Markov in Butwal municipality, Nepal.

Study Area
Geometrically, Butwal is located from 83˚22' 52.70"-83˚30' 22.61" E and 27˚39' 56.65"-27˚44' 55"N and on both banks of Tinau River that follows from north to south (Fig. 1).Historically, Butwal was known as Batauli and it was designated as municipality in 1959 AD.Average annual growth of urban population in Butwal is 4.62 percent between 2001-2011(MPRC 2013/14) that is relatively lower as compared to previous inter-censual rate of 5. 47 percent, 6.96 percent and 5.83 percent between 1991-2001, 1981-1991 and1971-1981 respectively and also slightly less than of national level (6.4%) (Sharma 2005).It is ranked 10 th on specialization of non-agriculture activities.The spatial extent of Butwal Municipality is 6822.81ha with a maximum eastwest distance of about 12 km and north-south length of about 11 km.It is located on strategic position of the crossroads between the east west and the north-south Siddhartha highway Functional base of municipality is expanding stronger, as manufacturing industries, construction and repair works, banking, insurance, education and catering services, and tourism.Over 80 percent of the economically active population is engaged in non-primary production .works,but the hideous part of Butwal is its growing squatter settlement having about one-third of the total households.

Data Set
A time series of remote sensing data including vector layers and satellite images was used to generate land use and cover maps.Topographic vector layers of 1993 at a scale of 1:25000 were employed for producing suitability maps required for land use cover maps.Landsat-5 TM, Landsat-7 ETM+ and landsat-8 OLI images (path 142, row 41) with 30 m spatial resolution dated 13-12,1999 , 5-10,2006 and 5-10,20013 downloaded at ftp://glcf.umiacs.umd.edufrom Global Land Cover Facility of the University of Maryland, USA were used in this study.Aster DEM downloaded at http://www.gdem.aster.ersdac.or.jp/ having 30 m spatial resolutions was used to derive slope map.Scanned toposheet with the scale of 1:25000 were used for image rectification at 0.5 accuracy level.The affine transformation and nearest neighborhood methods were used.

Image Classification and Accuracy Assessment
All satellite images including Aster DEM were rectified to Modified Universal Transverse Mercator (MUTM) projection with at least 20 well distributed ground control points matching with the vector layers.At first geometric correction was carried out using scanned toposheets no 098-12, 098-16 and 099-13 with less than a pixel of root mean square error.Then rectified toposheets were used to register all satellite images including Aster DEM through map to image registration using the first degree polynomial equation of image transformation.All images were resembled to 20 m pixel size using nearest neighbor resampling method in order to avoid altering the original pixel values of the image data.The digital vector layers of the municipality as land capability, settlement road, river and landuse were used for determing factors of suitability map.The accuracy of the ground referenced data required for image classification was assessed by using combining GPS, Google-earth, ETM+ image, aerial photographs, and toposheets during the field Training areas of 20-30 on the image for each class were selected as the sample signature using different band combinations for obtaining appropriate indication of discrimination among the classes.Maximum likelihood classifier (MLC), the most commonly and widely used algorithm, was run to classify the pixels.Smoothening of the resultant classified images to remove the salt and pepper appearance was achieved using a 3x3 majority filter.Producer's accuracy, user's accuracy and overall accuracy and kappa coefficient for 1999 and 2006 image classification were achieved (Mandal 2013).Overall accuracy and kappa statistics were found at 91.3 and 0.89 for 2013 image.These values are considered above the minimum value (85%) required in the identification of land use and cover types from imagery data (Murteira 1990).The detailed concepts and equations for quantitative measures of accuracy assessment applied in this research were mentioned in author's work (Mandal 2011).

Markov Chain Process
Markov chain framework was the first introduced by Brown in geography for the study of the innovation of diffusion and has since been followed by many scholars (Brown 1963cited in Collins 1975).Markovian chain is a stochastic process in which the probability, pij of a random variable X being in a state j at any point of time t+1 depends only on the state i; it has been at t, but not on states at previous points of time (Thomas and Laurence 2006).In this analysis, Markov Chain was produced the transition probability matrix from the two period images and based on this , a set of condition probability images for each land use/cover types were generated by analyzing with two period land use land cover images.Land Change Modular (LCM) in IDRISI is a change prediction modeler, which allows the transition likelihood of one pixel to be a function of nearest pixels, (Pontius & Chen, 2006).Markov Chain Process analyzes a pair of land cover images and outputs a transition probability matrix, a transition areas matrix, and a set of conditional probability images.The transition probability matrix is a text file that records the probability that each land cover category will change to every other category.The transition areas matrix is a text file that records the number of pixels that are expected to change from each land cover type to each other land cover type over the specified number of time units.In both of these files, the rows represent the older land cover categories and the columns represent the newer categories.The transition area matrix obtained from two time period was used as the basis for predicting the future LULC scenario.The detailed mathematical expression of Markov Chain equation is given in author's work (Mandal 2013).Markov chain process adapted in this study is shown in Figure 2  (Soe and Le 2006).The method is based on probability that a given piece of land will change from one mutually exclusive state to another (Thomas and Laurence 2006).These probabilities are generated from past changes and then applied to predict future change.However, a stochastic Markov model is not appropriate because it does not consider spatial knowledge distribution within each category and transition probabilities are not constant among landscape states; so it may give the right magnitude of change but not the right direction (Boerner et al. 1996).

CA-Markov Model
Cellular Automata (CA)-Markov incorporates the spatial component (Soe and Le 2006) and thereby adds direction to modeling.It has the ability to change its state, based on a rule that relates the new state to the previous state and those of its neighbors (Clarke and Gaydos 1998).It is implemented in LULC models that are able to simulate multiple land use types (Thomas and Laurence 2006).Hence, Cellular Automata Markov (CA-Markov) allows any number of categories and can simulate the transition from one category to another.CA-Markov is standard hybrid approach to model both spatial and temporal changes) the Markov process controls temporal dynamics among the LULC types through the use of transition probabilities ii) spatial dynamics are controlled by local rules through a CA mechanism considering either neighborhood configuration and transition probabilities (Sylvertown et al. 1992)

Multi-criteria Evaluation (MCE)
Multi-criteria evaluation(MCE) approach was adopted, in the situations in which a single decision-maker is faced with a multiplicity of usually incompatible criteria or in which a number of decision-makers must consider criteria, each of which depends on the decisions of all the decision-makers (Ademiluyi and Otun 2009).MCE is a decision support tool aiding a choice to be made between alternatives.The basis for a decision is known as a In a Multi-Criteria Evaluation, an attempt made to combine a set of criteria to achieve a single composite basis for a decision according to a specific objective.A decision need to be made about what areas are the most suitable for specific land use type development.
In this analysis, criteria or factors affecting suitability of specific land use types as agriculture, forest, built up and bush include bio-physical factors as slope, land capability, distance to road, distance to river and distance settlement.
The physical proximity to an existing land-use class is thought to be a driver of change into a particular land-use class in the future.Suitability maps for built-up, agriculture, forest and bush were generated from the MCE process in which weight was derived from the pair wise comparison in the AHP process (Figure 3).Also some restrictions such as developing an urban area into agricultural area were also taken into consideration for preparation of land-use suitability maps.
Figure 3. Suitability maps for agriculture, built up, forest and bush

Analytical Hierarchy Process (AHP)
AHP is a systematic procedure for developing factor weights required for preparing suitability map.The weights assigned to different factors were obtained by analytical hierarchy process (AHP).The larger the weight, the more important is the criterion in the overall suitability (Malczewski 1999).Factors are not hard rule like constraints; they allow the analyst to determining the degree of suitability from very low to high.The low to high suitability in LCM can be in the Table1.The weights computed from AHP process.
range value either from (0.0 to 1.0 range) in the real number or (0 to 255 ranges) in byte .Generally areas near an existing land-use class are likely to change into that class compared to areas that are far from that class.So, inverse J-shaped monotonic decrease function was used to determine the relative suitability of an existing land-use to change into that same land-use.The weight for the factors for preparing suitability map of agriculture, forest and built up were shown in Table 1.
In AHP, a pair wise comparison matrix is created by setting out one row and one column for each factor (Satty and Vargas 2001).Since the matrix is symmetrical, only the lower triangular half actually filled in (figure 4).In developing the weights, an individual factor compared with every other possible pairing, entered the ratings into a pairwise comparison matrix.In matrix, each cell is to be considered in terms of the relative importance of the row factors to its corresponding column factors.Factors or criteria were rated according to the following 9-point scale.Table 2. Factor weight computation from AHP process.
First few ratings were illustrated for considering the suitability of agriculture.It was observed that land use types of 2013 was very important than slope, and thus received a rating of 3 (Table2).This process were continued for all criteria or factors required for transition suitability maps for all land use types considered and all of the cells in the lower triangular half of the matrix were filled.The final factor weights were obtained along with consistency ratio less than 0.1 indicating good and then final weight were assigned to specify the relative importance of each factor in determining the aggregate suitability map.

Contiguity Rule in CA
CA-Markov is useful for modeling the state of several categories of a cell based on a matrix of Markov transition areas; transitional suitability images and a user defined contiguity filter.CA model applies contiguity rule like a pixel near to an urban area is most likely to be changed into urban area.Here, 5 × 5 mean contiguity filter was used on suitability images for each land cover class to define neighborhood relationship.The suitability of each pixel is determined by the pixel values within this defined filtering kernel and the more pixels of the same category of land cover exist in the neighborhood, the more the suitability value for that particular land cover type increases.Otherwise, the pixel value remains the same (Eastman, 2009).

Urban Land Use and land Cover Change Dynamics
With the increase of urban population significantly at about 44 percent within seven years from 2001 to 2008 and 145 percent from 1991 to 2008 ( Mandal 2013) Land use land cover change statistics using Markov and CA-Markov along with change percentage was given in Table 4.The results as shown in Table 4 showed that urban built-up area and forest were increased by 48.3% and 7.79 % between the period 1999 and 2013, whereas agriculture and water bodies were decreased by 30.26% and 28.22% .respectively.Unexpectedly, bush was also significantly increased between the period 199-2013 but it was decreased between 2006-2013.It is because of the fact that most of the shrub lands in the north found to be added into bush.The urban expansion has confined along the Siddhartha highway (north-south) and the Mahendra highway (east-west), as well as in the south of municipality, consuming agriculture land around the river bank
The transition area matrix obtained (table 6), was used as the basis for future LULC change prediction for year 2020.6 showed the probable area that might convert from one class to another in 2020.The results of predicted LULC scenario showed significant increase in built up by 6.18 percent and decrease agriculture by almost 20 percent respectively for the year 2020 (figure 4 and Table 4.The decrease in forest land by 11.53% municipality during 2013 to 2020 were high(

CONCLUSION
Unscientific utilization of land use and land cover due to rapid growth of urban population deteriorates urban condition.Urban growth, land use change and future urban land demand are key concerns of urban planners and decision makers essential for urban development and environmental management.It is aimed to predict urban land use change over time for better sustainable management of urban growth and development.Geo-information science has been used to study the urban change dynamics using Markov Chain and CA-Markov and predicted the magnitude and spatial pattern of the future urban scenario, based on past trend in an urban unit, Butwal Municipality, of Nepal.Using CA-Markov, spatio-temporal modeling of urban landuse change was performed using the probability transition matrix from the Markov chain process the suitability map of each land use/cover images and the contiguity filter.Suitability maps for built-up, agriculture, forest and bush were generated from the MCE process in which weight was derived from the pair wise comparison in the AHP process.In AHP process, land use types, slope, land capability, distance to road, settlement and water bodies as criterion of factor maps were evaluated with pair wise comparison in decision support system.The AHP weight was obtained consistency ratio for built-up, agriculture, forest and bush are 0.09, 0.07, 0.08 and 0.09 respectively considering acceptable threshold limit of 0.1.Thematic land use land cover types of base images as 1999, 2006, and 2013 of, TM, ETM & OLI sensors of Landsat, were classified using MLC algorithm of supervised classification accepting overall accuracies more than 90 and kappa statistics more than 0.87.The spatial extent increase from 1999 to 2013 in built up , bush and forest was observed to be 48.30percent,79.48percent and 7.79 percent, respectively, while decrease in agriculture and water bodies were 30.26 percent and 28.22 percent.The predicted urban land use land cover scenario for the year 2020 and 2027, with reasonably good accuracy would provide useful inputs to the urban planners for optimum management of the urban landscape.The study revealed the fact that the built up and bush expansion are the main driving force for loss of agriculture and river areas in Butwal Municipality and has the potential to continue in future.The abandoned area of river bed has been converted to built-up areas.The study utilizes three time period changes to better account for the trend and the modeling of urban dynamics.
The application of geo-information based integrated model that combines Markov and CA models for modeling, analyzing and predicting the changes in municipality LULC dynamics is demonstrated satisfactorily and effectively.The main challenges of the modelling and prediction work out is the suitability rating using multi-criteria decision-making.And analytical hierarcy process(AHP).The LULC reasonable states could be predicted by integrating biophysical and socioeconomic factors with the current LULC change in a municipality.

Fig 1 .
Fig 1.Location map of study area survey.Collection of ground reference data were carried out for land use and land cover classification.Four categories of land use and cover were extracted from TM, ETM + and OLI sensors.Supervised classification was used to train the algorithm to find similar pixels based on the established statistical parameters.Training areas of 20-30 on the image for each class were selected as the sample signature using different band combinations for obtaining appropriate indication of discrimination among the classes.Maximum likelihood classifier (MLC), the most commonly and widely used algorithm, was run to classify the pixels.Smoothening of the resultant classified images to remove the salt and pepper appearance was achieved using a 3x3 majority filter.Producer's accuracy, user's accuracy and overall accuracy and kappa coefficient for 1999 and 2006 image classification were achieved(Mandal 2013).Overall accuracy and kappa statistics were found at 91.3 and 0.89 for 2013 image.These values are considered above the minimum value (85%) required in the identification of land use and cover types from imagery data(Murteira 1990).The detailed concepts and equations for quantitative measures of accuracy assessment applied in this research were mentioned in author's work(Mandal 2011).

Figure 2 :
Figure 2: Flow chart of Markov chain process , iii) GIS and remotely sensed data is used to define initial conditions, to parameterize CA-Markov model, to calculate transition probabilities and to determine the neighborhood rules (Wang and Zhang 2001).CA -Markov is composed of quadruple elements as defined in the following equation (1) (White & Engelen, 2003): S) is the state of each cell representing any variable considered, in this case, it is land use land cover types.The state transition of Cell, S in CA is defined by following equation (of time periods St = initial state of a Cell at time, t S t +1 = new states of a Cell at next time t+1 P = transition probability matrix Spatio-temporal modeling of urban land use for year 2013 was predicted using CA-Markov process.For this, transition area matrix file required was derived from the LULC 1999 and 2006 using Markov Chain module in Idrisi Taiga and it was validated with existing LULC map 2013.As the validation results has sufficient reliable limit of kappa index of agreement (KIA), LULC maps for year 2020 was predicted using transition area matrix file obtained from 2006 to 2013.The 2013 LULC map was used as the base map for estimating future LULC scenario for the year 2020.

Figure 4 .
Figure 4. Predicted land use land cover

Table 3 .
Quantification of land use types (Area in ha).

Table 4 .
Land use land cover change statistics (Area in ha).

Table 5 .
Table 5 depicts the transition probability used for prediction of the land use type for 2020.Transition probability matrix

Table 6 .
Table 6 showed that 210 ha of agriculture land area has the probability of converting into the forest type.Similarly 316,131 and 65 ha area of agriculture have the probability to get converted into bush, built up and water bodies class from 2013 to 2020 respectively.Transition area matrix (2006-2013) for prediction of 2020.The statistical values in Table Table 4) .It is because of the fact that community forest program has not been so successfully that is started in the early 1980 to conserve the forest.