Computers, Environment and Urban Systems 35 (2011) 25–34

Contents lists available at ScienceDirect

Computers, Environment and Urban Systems journal homepage: www.elsevier.com/locate/compenvurbsys

Urban growth modeling of Kathmandu metropolitan region, Nepal Rajesh Bahadur Thapa *, Yuji Murayama 1 Division of Spatial Information Science, Graduate School of Life and Environmental Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, Japan

a r t i c l e

i n f o

Article history: Received 31 March 2010 Received in revised form 21 July 2010 Accepted 22 July 2010

Keywords: Land cover change Bayesian approach Weight of evidence Cellular automata Urbanization LUCC

a b s t r a c t The complexity of urban system requires integrated tools and techniques to understand the spatial process of urban development and project the future scenarios. This research aims to simulate urban growth patterns in Kathmandu metropolitan region in Nepal. The region, surrounded by complex mountainous terrain, has very limited land resources for new developments. As similar to many cities of the developing world, it has been facing rapid population growth and daunting environmental problems. Three time series land use maps in a fine-scale (30 m resolution), derived from satellite remote sensing, for the last three decades of the 20th century were used to clarify the spatial process of urbanization. Based on the historical experiences of the land use transitions, we adopted weight of evidence method integrated in cellular automata framework for predicting the future spatial patterns of urban growth. We extrapolated urban development patterns to 2010 and 2020 under the current scenario across the metropolitan region. Depending on local characteristics and land cover transition rates, this model produced noticeable spatial pattern of changes in the region. Based on the extrapolated spatial patterns, the urban development in the Kathmandu valley will continue through both in-filling in existing urban areas and outward rapid expansion toward the east and south directions. Overall development will be greatly affected by the existing urban space, transportation network, and topographic complexity. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Urban growth is recognized as physical and functional changes due to the transition of rural landscape to urban forms. The time– space relationship plays an important role in order to understand the dynamic process of urban growth. The dynamic process consists of a complex nonlinear interaction between several components, i.e., topography, river, land use, transportation, culture, population, economy, and growth policies. Many efforts have been made to improve such dynamic process representation with the utility of cellular automata (CA) coupling with fuzzy logic (Liu, 2009), artificial neural network (Almeida, Gleriani, Castejon, & Soares-Filho, 2008; Li & Yeh, 2002), Markov chain with modified genetic algorithm (Tang, Wang, & Yao, 2007), weight of evidence (Soares-Filho et al., 2004), non-ordinal and multi-nominal logit estimators (Landis, 2001), SLEUTH (Clarke, Hoppen, & Gaydos, 1997; Jantz, Goetz, Donato, & Claggett, 2010) and others (Batty, Couclelis, & Eichen, 1997; White & Engelen, 1997). Models based on the principles of CA are developing rapidly. CA approach provides a dynamic modeling environment which is well suited to modeling complex environment composed of large num* Corresponding author. Tel.: +81 29 853 4321; fax: +81 29 853 4211. E-mail addresses: [email protected], [email protected] (R.B. Thapa), [email protected] (Y. Murayama). 1 Tel./fax: +81 29 853 4211. 0198-9715/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compenvurbsys.2010.07.005

ber of individual elements. The land use change and urban growth process can be compared with the behavior of a cellular automaton in many aspects, for instance, the space of an urban area can be regarded as a combination of a number of cells, each cell taking a finite set of possible states representing the extent of its urban development with the state of each cell evolving in discrete time steps according to local transition rules. Therefore, CA based urban models usually pay more attention to simulating the dynamic process of urban development and defining the factors or rules driving the development (Batty et al., 1997). Different CA models have been developed to simulate urban growth and urban land use/cover change over time. The differences among various models exist in modifying the five basic elements of CA, i.e., the spatial tessellation of cells, states of cells, neighborhood, transition rules, and time (Liu, 2009). CA models have demonstrated to be effective platforms for simulating dynamic spatial interactions among biophysical and socio-economic factors associated with land use and land cover change (Jantz et al., 2010). While new urban models have provided insights into urban dynamics, a deeper understanding of the physical and socioeconomic patterns and processes associated with urbanization is still limited in developing countries in South Asia. Although, emerging geospatial techniques are bridging the spatial data gap recently, empirical case studies are still very few (Thapa & Murayama, 2009). This research aims to simulate urban growth in Kathmandu metropolitan region in Nepal using weight of evidence technique

26

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

incorporating with CA. As the result of population growth and migration from rural to urban areas, urbanization has been recognized as a critical process in metropolitan areas of Nepal (Bhattarai & Conway, 2010; Haack, 2009; Portnov, Adhikari, & Schwartz, 2007). The Kathmandu metropolitan region, capital and major tourist gateway, has been facing rapid urbanization over the last three decades. Recently, it has an estimated population of 2.18 million with an annual growth rate of 5.2% (Thapa & Murayama, 2010). Such urbanization pressure results rapid changes in the urban landscape pattern of the region adding more constructions and the loss of natural lands. Kathmandu, the capital of Nepal, has long history of development and exhibits a typical city surrounded by complex mountain terrains in the Himalayan region. History has witnessed its development as a strategic center of power, politics, culture and commerce (Thapa, Murayama, & Ale, 2008). However, along with the establishment of modern transportation infrastructures bringing easy access to the city, the agglomeration of rural settlements of Kathmandu valley into the city began in the early 1960s. The predominantly agricultural landscape gradually changed to an urban landscape with increasing human settlement in the 1960s and 1970s. The land changing process has escalated since the 1980s. Spatial diffusion of urban/built-up areas has spread outward from the city core and along the major roadways. Agricultural encroachment in rural hills and mountain peripheries and urbanization in the valley floor area are identified as the most common phenomenon in the valley (Bhattarai & Conway, 2010; Haack, 2009; Thapa & Murayama, 2009). Several urban land use development planning and policy initiatives for the valley have been made by the government in the past decades (Thapa et al., 2008). A latest planning document ‘Long Term Development Concept for Kathmandu Valley’ (Kathmandu Uptyakako Dirghakalin Bikas Avadharana) was released in 2002 (KVUDC, 2002). This document as planning reference conceptual-

izes scenarios to develop the Kathmandu metropolitan region by 2020. This long-term plan recommends the promotion of the tourism led service sector; guided urban development seeking compact urban form and the conservation of agricultural land; infrastructure development coordinated with land use; a new outer ring road to connect the traditional settlements in the metropolitan region; and rigorous regulation of areas defined as environmentally sensitive. All these policy recommendations eventually affect the future spatial pattern of urbanization. 2. Method 2.1. Study site The study area (Fig. 1) selected to apply the urban growth model follows the watershed boundary, which is derived from 20 m digital elevation data. Topography rises to elevation of 1100– 2700 m above the sea level that forms a bowl shaped valley. As most of the areas outside the watershed boundary are having high mountains, forest, shrubs land, and very low human settlements, therefore, urban expansion outside this boundary is largely restricted by these natural barriers. The valley is drained by the Bagmati river system which is the main source of water for drinking and irrigation in the valley (Thapa & Murayama, 2009). The study area covers 685 km2 whereas 14% of the land is defined as urban area that includes five urban centers, Kathmandu, Lalitpur, Bhaktapur, Kirtipur, and Madhyapur Thimi. In addition, the region consists of 97 sub-urban and rural villages. 2.2. Database preparation In this study, we used data from various sources for modeling, calibrating, and validating of urban growth in the Kathmandu

Fig. 1. Study area – Kathmandu valley, Nepal. Data source: ICIMOD/UNEP (2001)

27

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

metropolitan region. Three land use maps at 30 m spatial resolution for the years 1978, 1991, and 2000 were processed. These maps were acquired from Thapa (2009), which were created using remote sensing techniques. Due to the heterogeneity and complexity of the landscape in urban regions, for example, sub-urban residential areas forming a complex mosaic of trees, lawns, roofs, concrete, and asphalt roadways, require land use and land cover classification techniques that combine more than one classification procedure to improve remote sensing-based mapping accuracies. Therefore, a hybrid approach with series of processing steps was framed to create multi-temporal thematic maps of the Kathmandu valley. The detailed discussion on method adopted to create these maps and mapping accuracies can be found in Thapa and Murayama (2009). In this paper, it is assumed that these reference maps are the most accurate maps available for the particular points in time, so they serve as the basis to measure the accuracy of the prediction. The maps were further generalized to reduce the complexity in urban growth modeling. The land use types, i.e., built-up areas, industrial areas, roads, airport, institutional areas, government secretariat area, and royal palace were aggregated into built category. The other land use types, agricultural area, forest, shrubs, water, and open space were kept in the same state. Protected areas were not separated from the land cover maps, they were merged according to their belongingness, i.e., forest cover, open space, water, or built category. Two historical transition matrices were calculated from the land use maps for the period of 1978–1991 and 1991–2000. Based on the transition matrices analysis, the water and open space areas were excluded in dynamic modeling as they represented less land areas and mostly remained static. Therefore, we restricted the four broad land cover categories, i.e., built, agriculture, forest, and shrubs for the dynamic modeling of urban growth. Furthermore, the land cover transition rates for each category were computed normalizing row sum equal to 1. As the model was set to run in yearly time steps, the transition rates were further converted to annual rates by simply dividing the year differences, i.e., 13 and 9 for the land cover transition period of 1978–1991 and 1991–2000, respectively. The transition rates per year were passed to the modeling as fixed parameters. The other data, i.e., elevation, urban and village boundaries, population census, were also prepared. After creating all the required input maps in ArcGIS software, a simulation model of urban growth was designed in DINAMICA, a spatially explicit CA based modeling software (http://www.csr.ufmg.br/dinamica/dinamica.html). The transition matrices are passed onto the DINAMICA model that allocates the changes across the landscape based on spatial data layers representing physical and socioeconomic conditions which are stored in a GIS environment. As interactions between landscape elements occur in different ways, depending on local characteristics and transition rates, DINAMICA model produces distinct spatial patterns of land cover change. 2.3. Determining the factors of urban growth The selection of the best set of input variables and internal software parameters in order to produce the best fit between the empirical data and the observable reality are very important aspects of urban growth modeling. After calculating the existent land transitions, we identified a different set of factors governing each change of four land cover categories. In general, the urban growth that causes land use changes is a result of the complex interaction between behavioral and structural factors associated with the demand, technological capacity, and the social relations affecting demand and capacity, ultimately straining the environment (Thapa & Murayama, 2010). The factors available for the modeling analysis do not always represent the set of necessary variables able to produce ideal simulation results. However, there are no universal

driving factors of the change. Although similar driving factors have been found in several studies, the degree to which they contribute to landscape change differs (Almeida et al., 2008; Jantz et al., 2010; Verburg, Ritsema van Eck, de Nijs, Schot, & Dijst, 2004). People, government plans and programs, landforms, landscape change processes, and available resources often cause differences in the importance of various factors. In Kathmandu Valley, similarities among these factors are apparent. In fact, people’s behaviors and daily interactions with the environment over time have caused observable changes in the valley’s landscape (Thapa & Murayama, 2010). Indeed there is a set of factors for land use transitions that substantially respond to landscape changes where these factors effectively guide the modeling experiment. Several maps with biophysical, infrastructure, and social factors have been generated on the basis of the information extracted from land use maps and other data sources (Table 1). Digital elevation model (DEM) map was created based on the elevation point data (ICIMOD/UNEP, 2001), slope map was derived from the DEM map, and annual population growth rate map was prepared based on the census data (CBS, 2001). In all cases, distances were calculated using the Euclidean distance method. The distance to existing built area was defined as dynamic, update each time of model iteration, because of its neighborhood characteristics of land use change attraction. And all the remaining factors were static. Spatial independency of the input factor maps was checked using a set of measures, i.e. the Cramer test and the Joint-Uncertainty Information (Bonham-Carter, 1994). Both tests present the value between 0 and 1 showing the degree of association as independent to full, respectively, between the maps compared. As a principle, correlated variables must be removed or combined into a third that will be used in the model. Among the factors selected, no significant spatial dependency was revealed where Cramer’s test and JointUncertainty Information were found <0.5 and <0.6, respectively. 2.4. Transition probability map calculation and model calibration The development of abstraction methods capable of adequately representing complex processes with respect to quantity and location is a great challenge (Godoy & Soares-Filho, 2008). Weight of evidence, entirely based on the Bayesian approach of conditional probability, is traditionally used by geologists to point out areas favorable for geologic phenomena, such as seismicity and mineralization (Bonham-Carter, 1994; Goodacre, Bonham-Carter, Asterberg, & Wright, 1993). This method can combine spatial data from variety of sources to describe and analyze interactions, provide evidences for decision making and make predictive models (Almeida et al., 2008; Soares-Filho et al., 2004). In our spatial Table 1 Definition of the land cover change driving factors. Land change factors (biophysical, infrastructure and social)

Year

Digital elevation model at 30 m spatial resolution Slope in degrees Distances to rivers

1995 1995 1978, 2000 1978, 2000 1978, 2000 1978, 2000 1978, 2000 1978, 2000 1978, 2000 1991,

Distances to industrial estates Distances to five urban centers (Kathmandu, Lalitpur, Kritipur, Bhaktapur, and Madhayapur Thimi) Distances to major roads and highways Distances to ring road Distances to feeder roads Distances to existing built-up surface Annual population growth rate

1991, 1991, 1991, 1991, 1991, 1991, 1991, 2000

28

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

context, this approach detects the favorability of a certain event, for example, an event of land cover change from agriculture to built surface in relation to potential evidences (proximity to urban centers, roads, water, etc.) often called driving factors of change. Weights are estimated from the measured association between the land cover change occurrences and the values on the driving factors maps which are to be used as predictors. In this research,

we employed weight of evidence method for selecting the most important variables needed for the land cover change analysis and quantifying their influences to each type of land cover transition event. The weights of evidence represent each variable influence on the spatial probability of a transition i ) j can be calculated as follows (Eq. (1)):

OfDjBg ¼ OfDg Table 2 Simulation internal parameters. Land cover

Shrubs Forest Built Agriculture

1978

1991

2000

MPS

Var

Iso

MPS

Var

Iso

MPS

Var

Iso

81 27 31 49

21,457 3433 3755 5682

1.0 1.0 1.0 1.0

62 41 7 20

13,849 8245 114 2573

1.1 1.1 1.1 1.1

69 40 6 9

15,842 7510 222 438

1.1 1.1 1.1 1.1

Note: MPS = mean patch size; Var = variance; Iso = isometry.

PfBjDg PfBjDg

ð1Þ

where O{D|B} is the odds ratio of occurring event D, given a spatial pattern B, O{D} is prior odds ratio of event D, and PfBjDg=PfBjDg is known as the sufficiency ratio. In weight of evidence, the natural logarithm of both sides of Eq. (1) is taken as following (Eq. (2)):

logfDjBg ¼ logfDg þ W þ

ð2Þ

where W+ is the weight of evidence of occurring event D, which is calculated from the data. The spatial probability of a transition

Fig. 2. Land covers transition (1978–1991 and 1991–2000).

29

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

i ) j, given a set of spatial data can be expressed as following equation (Eq. (3)):

 x; y P i)j ¼

v

P

Wþ kn

i)jðVÞx;y ek P þ P W 1 þ ij e K kni)jðV Þx;y

ð3Þ

where v is a vector of k spatial driving factors, measured at location þ þ x, y and represented by its weights W þ 1xy ; W 2xy ; . . . ; W nxy being n the number of categories of each factor k (for detail mathematical discussion, see Bonham-Carter, 1994; Soares-Filho et al., 2004). In this way, weights of evidence are assigned for categories of each factor represented by its spatial data layers. After creating local transition probabilities, CA simulation model is calibrated by internal parameters, which concern the average size and variance of patches and patch isometry (Table 2). These functions enable the formation of a variety of sizes and shapes of patches of change. The patch isometry varies from 0 to 2. The patches assume a more isometric form as this number increases. The sizes of change patches are set according to a lognormal probability distribution, therefore, it is necessary to specify the parameters of this distribution represented by the mean and variances of the patch sizes to be formed (Soares-Filho et al., 2004). The mean patch size and variances were determined from the source maps while the isometry was determined empirically. The calibration parameters were computed for three temporal years, i.e. 1978 for calibrating the predictive model of 1991, 1991 for the model of 2000, and 2000 for projecting the future land cover patterns in 2010 and 2020. Furthermore, the reference map of time 2 were not used during the model calibration, for example, the land cover map of 2000 was not used as input map while calibrating the model of 2000. 2.5. Model validation Validation of a landscape dynamics model is usually carried out by comparing the predicted result to the empirical map to determine the prediction ability of the model. In this paper, we used a three-map comparison approach for model validation which is recommended by Pontius et al. (2008). This validation technique considers the overlay of all three maps: the reference map of time 1, the reference map of time 2, and the prediction map of time 2. This three-map comparison approach allows one to distinguish the pixels that are correct due to persistence vs. the pixels that are correct

due to change. In this paper, two validation maps were created for the model 1991 and 2000. For further quantitative clarification of each model, sources of percent correct and percent error were analyzed by computing observed change, predicted change, figure of merit (Eq. (4)), producer’s accuracies (Eq. (5)), user’s accuracies (Eq. (6)), and overall accuracies (Eq. (7)). The figure of merit (FoM) is the ratio of the intersection of the observed change and predicted change to the union of the observed change and predicted change. The FoM can range from 0%, meaning no overlap between observed and predicted change, to 100%, meaning perfect overlap between observed and predicted change. Eq. (4) shows mathematically notation of the FoM.

FoM ¼

B AþBþCþD

ð4Þ

where A is the area of error due to observed change predicted as persistence, B area of correct due to observed change predicted as change, C area of error due to observed change predicted as wrong gaining category, and D area of error due to observed persistence predicted as change. The producer’s accuracy (PA) (Eq. (5)) shows the proportion of pixels that the model predicts accurately as change, given that the reference maps indicate observed change. The user’s accuracy (UA) (Eq. (6)) shows the proportion of pixels that the model predicts accurately as change, given that the model predicts change. The overall accuracy (OA) (Eq. (7)) provides overall agreement between the reference and predicted maps.

B AþBþC B UA ¼ BþCþD BþE OA ¼ AþBþCþDþE

PA ¼

ð5Þ ð6Þ ð7Þ

where E in Eq. (7) is the area of correct due to observed persistence predicted as persistence. 3. Results 3.1. Land covers transition analysis Land covers transition matrix provided an important basis to analyze the temporal and spatial changes of land cover, and to

Table 3 Land covers transition area matrix (1978–1991). 1978

1991 Shrubs

Shrubs Forest Built Agriculture

Forest

Built

Agriculture

Km2

Rate

Km2

Rate

Km2

Rate

Km2

Rate

78.040 3.140 0.000 0.000

0.645 0.019 0.000 0.000

6.320 132.220 0.000 0.420

0.052 0.814 0.000 0.001

0.570 1.200 32.100 22.920

0.005 0.007 1.000 0.065

36.090 25.970 0.000 328.520

0.298 0.160 0.000 0.934

Table 4 Land covers transition area matrix (1991–2000). 1991

2000 Shrubs

Shrubs Forest Built Agriculture

Forest

Built

Agriculture

Km2

Rate

Km2

Rate

Km2

Rate

Km2

Rate

71.070 0.540 0.000 0.080

0.876 0.004 0.000 0.000

0.920 132.630 0.000 0.450

0.011 0.957 0.000 0.001

0.950 1.490 57.060 27.150

0.012 0.011 1.000 0.069

8.210 3.860 0.000 363.420

0.101 0.028 0.000 0.929

30

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

examine the driving forces behind those changes in the Kathmandu metropolitan region. Fig. 2 shows the landscape transition maps for the two time periods, i.e. 1978–1991 and 1991–2000. The maps demonstrated substantial landscape transitions during the study period. Agricultural area gained a large amount of land at the expenses of shrubs and forest lands during the period of 1978–1991 (Table 3). A large proportions of shrubs (36 km2) and forest (26 km2) lands were transformed into agricul-

tural land in the surrounding rural mountain areas in the region. This can be observed mostly in the northeastern and southern parts of the region. The built area received 23 km2 from agricultural land and 1.2 km2 from forest. Development of ring road around the existing urban core during the 1970s and extension of major and feeder roads into rural areas in the 1980s had accelerated built area expansion at the expense of agricultural area. Shrubs and forest

Fig. 3. Actual vs. simulated land cover patterns ((A) Actual 1978. (B) Actual 1991. (C) Simulated 1991. (D) Actual 2000. (E) Simulated 2000.).

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

land cover also contributed to built area at lower rate. Depending on the locations, the land cover transitions in between the forest and shrubs lands were also observed. During the period of 1991–2000, a large amount of agricultural land (27 km2) was transformed into built area which was increased by 4 km2 as compared to previous 13-years period (Table 4). But the transformation of shrubs and forest land cover into agricultural lands noticeably decreased. The shrubs (1 km2) and forest (1.5 km2) lands were also changed to built areas due to the expansion of rural roads in the 1990s. Overall, agricultural encroachment in rural hills and mountain peripheries was occurred where shrubs and forest landscape in rural areas of the valley mostly changed to agricultural areas. While in the valley floor, conversion of agricultural lands into built surface was identified. A small land cover transition between the forest and shrubs land covers was also noticed in this period. 3.2. Urban growth model validation results By varying the parameters in each model iteration, various simulation results were produced, the best results that the model has generated are illustrated in Fig. 3.C and E, which are compared with the actual maps of land cover from 1978, 1991 and 2000 (Fig. 3.A, B and D). Similarity in the spatial patterns between the simulated and reference maps are very important. A visual comparison of the model’s simulated results with the actual maps from 1978–2000 shows that the results produced by the model matched well with the actual urban extent. However, a systematic validation method requires quantifying the degree of error of the simulation results. The Fig. 4 shows the validation results graphically that allow the reader to access the nature of the prediction errors which are in various colors. These are obtained by overlaying the reference map of time 1, reference map of time 2, and prediction

31

map of time 2. The medium pink pixels show where the model predicts change correctly. Purple pixels show where change is observed and the model predicts change, however, the model predicts a transition to the wrong category, which is a type of error that can occur in multiple land cover category models. Medium blue pixels show error where change is observed at locations where the model predicts persistence. Light blue pixels show error where persistence is observed at locations where the model predicts change. White pixels show locations where the model predicts persistence correctly. Fig. 5 presents a summary of the error analysis according to the logic of the legend for Fig. 4. Each bar is a rectangular Venn diagram where the two central segments with different gray scale represent the intersection of the observed change and the predicted change. The second segment from the left shows the change that the model predicts correctly. The union of the segments on the left and center portions of each bar represents the area of change according to the reference maps, and the union of the segments on the center and right portions of each bar represents the area of change according to the prediction map. As represented by the figure of merit, 26% of overlap in the observed change and the predicted change is found in 1991 while this overlap is decreased by 7% in 2000. Decreasing pattern is also noticed in the producer’s accuracy and the user’s accuracy. However, the overall agreement between the reference and the predicted maps as shown by the overall accuracy is 91% in 2000, which is increased by 8% compared to 1991. The Table 5 evidences the results by land cover type at quantitative level, i.e., actual vs. simulated with minor differences. 3.3. Forecasting urban dynamic patterns for the years 2010 and 2020 Using the same configuration of the 1991–2000 simulation model and input map of 2000 (time 1) and road network of

Fig. 4. Validation results ((A) 1991 and (B) 2000) obtained by overlaying the reference map of time 1, reference map of time 2, and prediction map of time 2.

32

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

Fig. 5. Sources of percent correct and percent error in the model validation.

Table 5 Actual vs. simulated land cover in percentages (1991–2000). Land cover

Shrubs Forest Built Agriculture

1991

2000

Actual

Simulated

Actual

Simulated

12.14 20.78 8.58 58.50

12.01 20.82 8.70 58.47

10.71 20.04 13.02 56.23

10.72 20.10 13.02 56.16

2000, we performed a simulation aiming to project the spatial patterns of urban growth in the metropolitan region for the year 2010 and 2020. The Fig. 6 shows the urban growth consistently expanding eastwards agglomerating the sub-urban villages and the two urban centers, i.e., Madhayapur Thimi and Bhaktapur. The built surface of the villages in the southeastern part also starts agglomeration by 2010. Current agricultural area between the Madhayapur Thimi and Kathmandu–Lalitpur urban centers will be converted into built in the 2010s. By 2020, all the urban centers will be aggregated into a greater metropolitan region in the valley. The Fig. 7 shows the quantitative results produced from the simulated maps by the selected land cover categories. The built area will be increased from 87 km2 to 148.7 km2, 59% of increase, until 2020. While other land covers are diminishing at different level. Most vulnerable seems to be shrubs land changing to either agricultural area or forest. The agricultural area will be converted mostly to the built area. However, the urban growth rate is decreasing in the upcoming decades (Table 6). Although varies at rate, the trend of negative growth of forest and agricultural areas

Fig. 7. Status of land covers in km2 (2000–2020).

Table 6 Annual urban growth ratea (2000–2020). Land cover

Shrubs Forest Built Agriculture a

Years 2000

2010

2020

1.298 0.391 5.755 0.427

1.318 0.374 3.020 0.473

1.309 0.381 2.213 0.490

Rate = ((present area  past area)/(past area)  100)/number of past years.

Fig. 6. Simulated land covers patterns (2010 and 2020).

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

is increasing. But the negative growth of shrubs land will be improved by 2010 at lower rate.

4. Discussion and conclusions Modeling urban growth has been the objective of urban research for many years. This article analyzed historical land cover transition and simulated future urban dynamics in the Kathmandu metropolitan region using the Bayesian approach incorporating with CA and GIS techniques. The historical evidences of land cover transition explained the rate of encroachment of urban areas on other land cover has been quite rapid, with scattered patches of urban development characterizing the urban sprawl in metropolitan region. It may be due to conversion of agricultural lands to built areas in the urban fringes, which forced the farmers to migrate in vicinities in one hand. On the other hand, due to road expansion and market accessibility to rural areas, farmers were encouraged to develop agricultural activities in the rural hills that enhanced to deforestation and encroaching shrubs lands. After the establishment of democracy in 1990, Kathmandu became the center of political power and hub of business activities (Thapa et al., 2008). The business and economic opportunities led to a population influx in the valley putting housing demands that eventually increased the built surfaces. Landscape fragmentation and heterogeneous land use development are recent phenomenon in the area (Bhattarai & Conway, 2010; Thapa & Murayama, 2009). This type of development around the city often described as economically inefficient and aesthetically unattractive (Cadwallader, 1996). The form of urban sprawl, with its uneven spatial directions, could have been promoted by weakness of the local government regulations as reported in Thapa et al. (2008). The land tenure system could also have contributed to this form of sprawl since private land has been subdivided and used for development of unplanned residential areas. The modeling framework adopted in DYNAMICA, which is used in this research, can be compared to state change model concept which was first used to model urban development patterns in CUF II model (Landis, 2001). The state change models project future land use and land cover using the information on land uses at two points of time to calibrate a statistical model which relates a set of independent variables to the observed land use changes at each location (Klosterman & Pettit, 2005). Current and projected values for the independent variables are then used to project future land use changes. Similar concept is applied in DYNAMICA differentiating with probability functions, i.e., the DYNAMICA uses weight of evidence while CUF II uses non-ordinal logit (Landis, 2001). Recently, the other models, i.e., new version of Clarke’s urban growth model, SLEUTH-3r (Jantz et al., 2010) and Land Change Modeler of IDRISI (Eastman, 2009) adapted similar state change concept by differentiating in probability function. In order to derive future spatial pattern of urban development in the valley, the calibration of the model was conducted systematically over time with the similarity of the model’s outcomes being compared to the actual urban development of the metropolitan region between 1978 and 2000. The model was applied to generate perspective views of the city for the next 20 years. Results presented dynamic patterns and composition of Kathmandu urban development up to the year 2020. The perspective views of Kathmandu generated by the model highlighted the requirement of urban planning controls, leading to the conservation of more forest and agricultural lands and the promotion of intensive development. This analysis suggests that urban development in the Kathmandu valley will continue through both in-filling in existing urban areas and outward expansion toward the east, south, and west directions in the future. Development will be greatly affected by the existing urban space and the transportation network. Unsur-

33

prisingly, under extremely unfavorable topographic conditions, close proximity to road network play a crucial role in urban development, much more than in cities with milder topographies and better infrastructure networks (Liu, 2009). However, this study assumed that transport systems would not change during the future simulation period. There are substantial uncertainties in simulating future changes of road systems since they are subject to relatively frequent changes and are often affected by urban transportation policies and land use planning. The map of the future transportation network is unavailable in the valley. The topographical constraints on this development will also be important, especially when the urban areas of Kathmandu extend further outward from existing urban areas. The topographical complexity and unavailability of future transportation network map could be the main cause to decrease urban growth rate in future. Furthermore, a concrete future land use plan for the Kathmandu valley is missing in this modeling. The spatial data related to future urban development planning, such as land use plans, are still lacked in the authority. Therefore, the model used in this research for simulating future growth scenario missed the planning guidance. The simulation estimates are based on extrapolation from historic processes which are not guaranteed to continue in the future but it mirrors spatial patterns of land cover in the metropolitan region if the historic process is not altered. In this case, the model has generated maps to show where and how the urban development in Kathmandu is heading in the next two decades from 2000, which may be a critical reference to make decisions for guiding future urban development and land management in the valley. This study has also demonstrated the usefulness of the data acquired from satellite remote sensing and CA based urban growth modeling in providing land use and land cover maps and change information, which are very valuable for planning and research. The approach adopted in this study can be used for the analysis of urban growth and land cover changes in developing countries where the amount and quality of geographic information and other ancillary data are very limited. Acknowledgements The authors wish to thank to the three anonymous reviewers for their creative comments and suggestions that helped us to improve this manuscript substantially. The financial support for this research from Japan Society for Promotion of Science (Grant #2109009) to study spatial process of urbanization and its impact on environment in Kathmandu is greatly acknowledged. References Almeida, C. M., Gleriani, J. M., Castejon, E. F., & Soares-Filho, B. S. (2008). Using neural networks and cellular automata for modeling intra-urban land-use dynamics. International Journal of Geographical Information Science, 22, 943–963. Batty, M., Couclelis, H., & Eichen, M. (1997). Urban systems as cellular automata. Environmental and Planning B, 24, 175–192. Bhattarai, K., & Conway, D. (2010). Urban vulnerabilities in the Kathmandu valley, Nepal: visualizations of human/hazard interactions. Journal of Geographic Information System, 2, 63–84. Bonham-Carter, G. (1994). Geographic information systems for geoscientists: modeling with GIS. New York: Pergamon. BS, C. (2001). Population of Nepal (selected data – central development region). Kathmandu: His Majesty’s Government of Nepal. Cadwallader, M. T. (1996). Urban geography: an analytical approach. New Jersey: Prentice-Hall. Clarke, K. C., Hoppen, S., & Gaydos, L. J. (1997). A self-modifying cellular automaton model of historical urbanization in the San Francisco Bay area. Environment and Planning B, 24, 247–261. Eastman, J. R. (2009). IDRISI Taiga, guide to GIS and remote processing. Guide to GIS and Remote Processing. Worcester: Clark University. Godoy, M. M. G., & Soares-Filho, B. S. (2008). Modelling intra-urban dynamics in the Savassi neighbourhood, Belo Horizonte city, Brazil. In M. Paegelow & M. T. C. Olmedo (Eds.), Modelling environmental dynamics (pp. 319–338). Berlin: Springer.

34

R.B. Thapa, Y. Murayama / Computers, Environment and Urban Systems 35 (2011) 25–34

Goodacre, C. M., Bonham-Carter, G. F., Asterberg, F. P., & Wright, D. F. (1993). A statistical analysis of spatial association of seismicity with drainage patterns and magnetic anomalies in western Quebec. Tectonophysics, 217, 285–305. Haack, B. (2009). A history and analysis of mapping urban expansion in the Kathmandu valley, Nepal. The Cartographic Journal, 46, 233–241. ICIMOD/UNEP (2001). Kathmandu valley GIS database. Nepal: Kathmandu. Jantz, C. A., Goetz, S. J., Donato, D., & Claggett, P. (2010). Designing and implementing a regional urban modeling system using the SLEUTH cellular urban model. Computers, Environment and Urban Systems, 34, 1–16. Klosterman, R. E., & Pettit, C. J. (2005). Guest editorial: an update on planning support systems. Environment and Planning B, 32, 477–484. KVUDC (2002). Long term development concept of Kathmandu valley. Kathmandu: Kathmandu Valley Urban Development Committee. Landis, J. (2001). CUF, CUF II, and CURBA: a family of spatially explicit urban growth and land-use policy simulation models. In R. K. Brail & R. E. Klosterman (Eds.), Planning support systems: integrating geographic information systems, models and visualization tools (pp. 157–200). Redlands: ESRI Press. Li, X., & Yeh, A. G. (2002). Neural-network-based cellular automata for simulating multiple land use changes using GIS. International Journal of Geographical Information Science, 16, 323–343. Liu, Y. (2009). Modelling urban development with geographical information system and cellular automata. Boca Raton, FL: Taylor and Francis Group. Pontius, R. G., Boersma, W., Castella, J., Clarke, K. C., de Nijs, T., Dietzel, C., et al. (2008). Comparing the input, output, and validation maps for several models of land change. Annals of Regional Science, 42, 11–47.

Portnov, B. A., Adhikari, M., & Schwartz, M. (2007). Urban growth in Nepal: does location matter? Urban Studies, 44, 915–937. Soares-Filho, B. S., Alencar, A., Nespad, D., Cerqueira, G. C., Dial, M., Del, C., et al. (2004). Simulating the response of land-cover changes to road paving and governance along a major Amazon Highway: the Santarem-Cuiaba corridor. Global Change Biology, 10, 745–764. Tang, J., Wang, L., & Yao, Z. (2007). Spatio-temporal urban landscape change analysis using the Markov chain model and a modified genetic algorithm. International Journal of Remote Sensing, 28, 3255–3271. Thapa, R. B. (2009). Spatial process of urbanization in Kathmandu valley, Nepal. PhD Dissertation. Graduate School of Life and Environmetnal Sciences, University of Tsukuba, Ibaraki, Japan. Thapa, R. B., & Murayama, Y. (2009). Examining spatiotemporal urbanization patterns in Kathmandu valley, Nepal: remote sensing and spatial metrics approaches. Remote Sensing, 1, 534–556. Thapa, R. B., Murayama, Y., & Ale, S. (2008). Kathmandu. Cities, 25, 45–57. Thapa, R. B., & Murayama, Y. (2010). Drivers of urban growth in the Kathmandu valley, Nepal: examining the efficacy of the analytic hierarchy process. Applied Geography, 30, 70–83. Verburg, P. H., Ritsema van Eck, J., de Nijs, T., Schot, P., & Dijst, M. (2004). Determinants of land-use change patterns in the Netherlands. Environment and Planning B, 31, 125–150. White, R., & Engelen, G. (1997). Cellular automata as the basis of integrated dynamic regional modeling. Environment and Planning B, 24, 235–246.