Proximate Causes of Land-Use and Land-Cover Change in Bannerghatta National Park: A Spatial Statistical Model

Abstract

Land change modeling has become increasingly important in evaluating the unique driving factors and proximate causes that underlie a particular geographical location. In this article, a binary logistic regression analysis was employed to identify the factors influencing deforestation and simultaneous plantation driven reforestation in Bannerghatta National Park, located at the periphery of one of the fastest growing cities in India, i.e., Bangalore. Methodologically, this study explores the inclusion of different sub-regions and statistical population to address spatial autocorrelation in land change modeling. The results show negative relationship between deforestation and protected area status and edge of previous forest clearing. In addition, the deforestation models found differences in the processes that are affecting forest clearing in our two sub-periods of 1973–1992 and 1992–2007. The plantation driven reforestation in the region were attributed to distance to major towns, Bangalore city, rural centers and major and minor roads suggesting the importance of accessibility to market for heavy cash crops such as coconut palm and eucalyptus. Finally, the inclusion of different sub-regions and statistical population facilitated a better understanding of varying driving factors in different zones within the overall landscape.

1. Introduction

One of the major research themes of human dimensions of global change, in the recent past, has been to explain various driving factors of land-use and land-cover change (LULCC). Case study evidence at the local level (landscape level) has been argued to better support the varying drivers behind LULCC [1], even though LULCC are studied at various “spatiotemporal scales” [2]. Our study region, Bannerghatta National Park (BNP) and its surroundings, situated in the southern part of India, provides an interesting case study for building a spatially explicit model to understand the interactions between LULCC and its driving factors. BNP is situated just 22 km south of the Bangalore urban area, which is one of the fastest growing cities in India. In spite of the close proximity to the city of Bangalore, rapid suburbanization and population growth, BNP and its surroundings have undergone recovery of forest cover inside the park after an initial deforestation trend between 1973 and 1992. Simultaneously a spatially distributed and patchy native forest regeneration and tree plantation driven reforestation trend is observed in the overall landscape [3]. Moreover, the rapid horizontal expansion of Bangalore urban area also indicates future pressure on the park and the ecosystem services its forest cover provides to the region. These localized LULCC patterns and trends certainly raises questions, e.g., what factors are important and how they are related to these changes.

Spatial-statistical models of LULCC are one of the important tools for quantifying and explaining the drivers of LULCC and to predict future potential LULCC [4,5]. Empirical models have been widely used for detailed analysis of case studies and to explore key driving factors affecting LULCC [6,7]. Further, building these empirical models for describing LULCC is made possible by data derived from satellite remote sensing and geographic information systems (GIS) which provides quantitative data at different spatial and temporal scales [7,8]. A commonly used empirical model to explain the spatial and temporal patterns of LULCC is a binary logistic regression model, which explains the probability of occurrence of a particular land cover category [7,9]. Logistic regression models have been used to identify proximate and underlying causes of urban development [10,11,12], agricultural expansion [5,13], development along rural–urban gradients [14], wildlife habitat studies [15,16], and deforestation [17,18,19,20].

Spatial-statistical models have been widely applied to land change studies, however, it has been criticized for its inability to incorporate spatial variability in the processes affecting LULCC [5]. LULCC in most regions of the world are very spatially heterogeneous. Along with the heterogeneity of LULCC comes the spatial variability in various biophysical, socio-economic, and political driving factors which affect LULCC at the macro level (regional, national and global level) [5]. This spatial heterogeneity in LULCC and its driving factors gives rise to spatial variability in its proximate causes affecting LULCC at the local level. Although the fundamental idea behind modeling LULCC processes is “to transcend the complexity of context, seeking to identify broad and universally applicable forces of change that crosscut the circumstances of place and period” [21], it is crucial that we understand the unique relationship between LULCC and its driving factors operating at a particular geographic location and time period [17,21]. Aggregating the spatial variability of LULCC and their driving factors or its proxy variables in one model would obscure the predictability of the model and produce weak links between the two [5]. To understand the uniqueness of the different LULCC processes, it is thus necessary to find the optimal spatial entity for each LULCC by choosing a relatively homogenous region where a particular LULCC is occurring; it would also require different model parameterizations [17].

Methodologically, this study would expand on previous modeling work on LULCC in two ways. First, to address the optimal spatial entity the research identifies different statistical population to represent each LULCC, addressed in literature theoretically but not applied statistically. Secondly, this research addresses spatial autocorrelation eminent in a pixel based regression model of LULCC, which is a common practice in geographical and ecological studies [22,23,24]. Spatial autocorrelation has been addressed in various spatial models by appropriate sampling schemes, scale, quantification of spatial patterns and various other statistical methods for use with spatial data [19,25,26,27]. Strength of this modeling approach is the inclusion of sub-regions as binary explanatory variables after careful study of residual cluster.

The overall goal of this research is fourfold: (1) to elucidate the relationship between LULCC in and around BNP and its various proximate causes and to infer about the underlying driving factors; (2) to determine whether there are differences in processes that drive land cover changes in different eras; (3) to address spatial autocorrelation and determine whether the explanatory variables that are driving land cover changes are different in different sub-regions within the BNP and its surroundings; and (4) to use different statistical populations for different LULCC, as addressed in our Section 3.1.1. The complexity of the driving factors of land change processes makes it challenging to capture time (e.g., temporal variations in road connectivity and changes in population density), space (effect of different spatial scales and different sub-regions) and human decision making (land tenure and policy changes) in a single model and thus the abundance of various modeling approaches to address different set of variables [7,28,29]. Our study focuses on distance based variables such as distance to roads, towns, villages, forest edge, and water bodies; attributes of the physical environment, such as elevation and slope; and zoning policies, i.e., protected area status in BNP and its surroundings. To account for spatial variability in proximate causes as well as address spatial autocorrelation, this research constructs various zones within the overall landscape after careful analysis of spatial pattern of LULCC in the region and preliminary statistical analysis. A combination of remote sensing and GIS techniques are used for data extraction and analysis.

2. Site Description

The total study area (638 km²) encompasses BNP (109 km²) and a zone of 5 km buffer around the park. It is situated in the southern part of Karnataka State, 22 km south of the city of Bangalore (Figure 1). BNP consists mainly of dry deciduous and scrub forests under the Terminalia–Anogeissus latifolia–Tectona series [30] with scattered patches of eucalyptus plantation and moist deciduous forest. The park is one of the oldest habitats of Asiatic elephants [31]. The study area is undulating with few rocky hillocks. The elevation within BNP ranges between 700 and 1000 m. Overall, the eastern part of the study area has a higher elevation than the western part of the study area and could be seen broadly as two zones of high and low elevation. The study area gets most of its rainfall during June–November and varies between 625 and 700 mm annually and the mean annual temperature is around 27 °C. There are six village communities inside BNP, which form three major enclosures inside the park. The main economic activity of the village communities inside the park is farming and animal grazing.

The 5 km buffer area includes Talli Reserve Forest and Bilikal Reserve Forest, which are situated in the southeastern and southern part of BNP, respectively, and are contiguous to the national park. The studies buffer distance is considered a suitable limit since most of the interaction between the people and park such as illegal logging and non-timber resource extraction took place within a 5 km buffer (personal communication and interviews). For our statistical analysis, we considered BNP, Talli and Bilikal Reserve Forest as one single entity of protected area and did not differentiate between the different levels of protection that is given to a national park and a reserve forest in India. To the western and eastern parts of BNP lie agricultural fields and several hundred village communities, village market areas and small towns. To the north and northeastern part of BNP lies the dense built area of suburban Bangalore, one of the fastest growing metropolitan areas in India. Farming, grazing, illegal sand mining, stone mining, and factory employment in the larger Bangalore’s Metropolitan area and surrounding towns are the most important economic activities of the village communities in the study area. These varied economic activities are reflected in the land-cover and land-use of agricultural land (includes the grassland and the fallow land), plantation, stone waste/mines, built area, water bodies and bare area in the buffer area. Furthermore, the park is very well connected with the city of Bangalore, surrounding village communities in the buffer area and various small towns through both paved and unpaved roads. Many of these paved and unpaved roads also run through BNP and the village communities inside the park, making the park accessible to people living in the buffer area.

BNP and its surroundings have experienced two major LULCC transitions between 1973 and 2007 details of which are provided in a previous study [3]. The transition between forest cover to non-forest cover (deforestation) and non-forest cover to tree plantation (plantation driven reforestation) are the two major land-use and land-cover changes that have occurred in the region. Most of the deforestation has occurred between 1973 and 1992 and most of the tree plantation increase has taken place between 1992 and 2007 (Figure 2). Deforestation is located mostly on the north central part of the study area outside BNP. Some deforestation has also occurred on the edges of the national park and near the village enclosures. The eastern and western parts of the study region outside BNP have recently experienced an increase in tree plantations which are mostly agroforestry crops in this region. Some of these tree plantations are coconut palm (Cocos nucifera), mango (Mangifera indica), teak (Tectona grandis) and eucalyptus (Eucalyptus cinerea). Many of these tree plantation areas were previously used for commercial and subsistence agricultural crops such as rice, raggi (Finger millets), vegetables, banana, mulberry and flowers.

3. Methods

3.1. Spatial Data

Land cover maps of the study site were created using Landsat Multi-Spectral Scanner (MSS) (United States Geological Survey [32]), Thematic Mapper (TM) (United States Geological Survey [32]), and Indian Remote Sensing Satellite Linear Imaging Self-Scanning (IRS LISS) (National Remote Sensing Center, Hyderabad, India) III images for 1973, 1992, and 2007, respectively. For mapping land cover, we used a hybrid classification approach of combining results from ISODATA clustering method and a supervised Gaussian maximum-likelihood classification method. All the satellite images were geometrically corrected and projected to UTM WGS 84. The resampling was done using the nearest neighborhood algorithm with a root mean square error of less than 0.5 pixels. Additionally a rule based classification was used to get the land-cover classes of native forest, tree plantation and non-forest. Our study focused on a simplified spatial extension of forest from a bird’s eye point of view of a multispectral satellite image, hence, limited our capability to map the complex forest structure.

A post-classification change analysis was applied to image pairs 1973–1992 and 1992–2007 by multiplying across the two dates for the three classes of native forest, tree plantation and non-forest. Out of the 27 trajectories, nine trajectories were chosen to show changes representing more than 1% of landscape. The land-cover maps as well as the changes in the landscape were validated using field data i.e., training samples used for running accuracy assessment, interviews on historical LULCC and other ancillary data (topographical maps). The details of image dates, change trajectory results and interview method have been discussed in a previous study [3,33]. The LULCC map layers were created using Erdas Imagine 9.3 (Hexagon Geospatial, Madison, AL, USA). Digital elevation model (DEM) was acquired from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) (NASA’s Earth Observing System, Washington, DC, USA) with a spatial resolution of 30 m. Slope was calculated from the elevation map using ArcGIS (ESRI, San Diego, CA, USA).

The road networks (major and minor roads), towns, villages, and water bodies (stream/rivers and lakes) were manually digitized from the 1:50,000 Survey of India topographic maps. Any unpaved roads, village roads and forest pathways were grouped as minor roads. All paved roads such as national and state highways were grouped as major roads. Park boundary was mapped from the topographic maps and the park management plan from the Karnataka Forest Department, Bangalore. Other important data layers were created from the land-cover maps of 1973, 1992 and 2007 i.e., distance from the edge of the forest. The distance layers (distance to road, water, town and villages) were created using the Euclidean distance tool of ArcGIS which calculates for each pixel the Euclidean distance to the closest source. A study area mask was created and used in the model to avoid the calculation of the background value. All data were resampled to a spatial resolution of 30 m as our land-cover change data (dependent variable) were at that scale.

3.1.1. Dependent Variables

Two dominant changes were observed in the 30-year time period from 1973 to 2007. The first dominant change in the region is forest loss in the Bannerghatta National Park (BNP) and its surroundings. This forest loss is mostly on the edges of the BNP and on the north-central part of the region outside the BNP. The forest loss is higher in the time period 1973–1992 as compared to the time period 1992–2007. The second dominant LULCC has been a massive increase in tree plantation in the region outside the BNP. The quantity of tree plantation gain is more in the second time period (1992–2007) than in the first time step (1973–1992). One spatial distinction could be made between the tree plantation on the eastern and the western part of the region. The eastern part of the study area has more eucalyptus plantation while the western part has more coconut palms and mango plantation with small patches of teak plantation. Further, the western part of the study area has a higher tree plantation gain than the eastern part.

These two dominant changes are modeled in this paper and make our dependent variables (

Table 1. GIS database.

	Variables	Type	Unit	Abbreviations
Dependent Variable
	Deforestation, 1973–1992	Binary	0–1
	Deforestation, 1992–2007	Binary	0–1
	Plantation Gain 1973–1992	Binary	0–1
	Plantation Gain 1992–2007	Binary	0–1
Independent Variables
Relief Related Variables	Elevation	Continuous	Meter	ELEV
	Slope	Continuous	Degrees	SLOPE
Proximity Variables	Distance to roads    All Roads    Major Roads    Village Roads
		Continuous	Meter	ALLRDDIST
				MJRDDIST
				MIRDDIST
	Distance to edge of the forest
	1973	Continuous	Meter	FORDIST73
	1992			FORDIST92
Zoning Policy	Distance to Bangalore	Continuous	Meter	BANGDIST
	Distance to towns	Continuous	Meter	MJTOWNDIST
	Distance to villages	Continuous	Meter	VILLDIST
	Distance to water	Continuous	Meter	WATERDIST
	Presence or absence of protected area	Binary	0–1	PA

), forest loss and tree plantation gain, both of which will be unidirectional models. The forest loss would be termed as deforestation in the absence of data availability to give a distinction between deforestation and degradation. Deforestation is defined in this paper as any area that was native forest on the earlier date and non-forest in the recent date. Similarly, tree plantation gains are areas which were non-forest in the first time step but became tree plantation in the second time step. The binary outcome variables of change and no change i.e., areas which were deforested as change (1) and rest of the landscape as no change (0) were created by masking all change trajectories but forest loss trajectory. Within the 5 km buffer zone, areas where high density of major roads and towns intersected in the north and northeastern parts of the study area were clipped out for better predictability of our deforestation model. To simplify the threshold determination of the clipped area, both the density maps were classified using natural breaks and a threshold of more than 1.2 km/km² (road density) and 0.01 km/km² (major town density) were considered appropriate threshold for elimination for our statistical population. This area is also part of the suburban extension of Bangalore city and never had any forest cover. Similarly, our binary outcome variable of tree plantation gain as change (1) and rest of the landscape as no change (0) was mapped out. For our tree plantation gain model, BNP, Talli and Bilikal Reserve Forest were clipped out and not included in the statistical analysis. Tree plantations in this region are agricultural crops. It has occurred mostly outside these protected areas. The factors that would drive tree plantation gain inside a protected area would be much different, e.g., reforestation strategies by various forest management policies. However, tree plantation outside the park should be driven by factors other than reforestation policies of political institutions.

3.1.2. Independent Variables

For BNP and its surroundings, the independent variables are elevation, slope, distance from the edge of the forest 1973, 1992 and 2007, distance from the towns, villages, Bangalore, water bodies, and presence and absence of protected area (Table 1). It is hypothesized that deforestation and tree plantation gain are influenced by the following mentioned variables.

Elevation and slope: Prior studies found topography to influence the spread and extent of deforestation [34,35]. Steeper slopes are least favored for land uses such as agricultural practices, infrastructure development and residential and commercial development [14,35,36]. Further, studies have found deforestation and slope gradient to be negatively correlated [25,37], as logging activities may be limited in the upper slopes because of inaccessibility. Tree plantation is an agricultural crop in this region and thus we would expect it to be on flatter slopes. Further, forest cover should remain intact at higher elevations and steeper slopes as steeper slopes are not preferred for land uses such as agricultural practices, road building, and residential and commercial development.

Distance to road: Various empirical studies have found deforestation to increase with greater access to forest and markets, with roads, rivers and railroads facilitating this access [38]. Prior studies of spatial regression models found strong relationships between roads and deforestation [17,36,39,40]. Thus, this variable is important as closer proximity to roads may encourage deforestation. Agriculturists favor closer proximity to roads as it provides access to market [36] and thus this variable is important for our tree plantation gain model. Distinctions between major and minor roads were made to assess the impact of each of them on deforestation and tree plantation gain separately. Further, we kept all roads as an added explanatory variable as they represented better the landscape level flow of the road network.

Distance to edge of the forest: Higher deforestation should be taking place closer to the edges of BNP. Distance to edge of the forest is not important for tree plantation gain as most of the forest cover is inside the national park and tree plantation gain is unlikely to increase inside the national park, as it is an agricultural crop in the region.

Distance to Bangalore and other towns (smaller urban centers): We used both statutory and census towns to represent towns in our study [41]. Statutory towns are places with a municipality, corporation, cantonment or notified town area. Census towns are defined as urban areas with: (a) minimum population of 5000; (b) population density of at least 400 per km²; and (c) 75% male population working on non-agricultural sector. This variable is an important explanatory variable as BNP is situated just 22 km south of the Bangalore city, which is growing at a fast rate as result of Information Technology related company development and increasing built area development [42]. Distance to Bangalore was calculated as a series of buffers of 1 pixel expanding from the center of Bangalore. A second variable, distance to all towns, was also created using the Euclidean distance tool of ArcGIS which calculates the Euclidean distance to the closest source (town). Distance to towns and Bangalore city are not an important explanatory variable for deforestation model as most of the deforestation took place far away from town. However, these two variables could be important for tree plantation as these towns provide market for these agricultural crops.

Distance to villages (rural centers): Villages for this study area represents the revenue villages of India [41]. This variable is important for both the deforestation and tree plantation gain model. There are seven villages inside the national park, which forms three enclosures inside BNP. The buffer area has more than 100 villages where population varies from 500 to 4000 people. No distinction was made between different sizes of the villages because of lack of village level population data. Closer proximity to villages is preferred by agriculturists for tree plantation crops. Further, people–park conflicts and various social-economic activities practiced by village communities are described as one of the major driving factors behind deforestation in heavily populated countries like India where many parks are situated in densely populated areas [43,44], which makes distance to villages an important variable for the deforestation model as well.

Distance to water: Closer proximity to water is valued by agriculturalists for tree plantation crops such as coconut palm, mango and eucalyptus. Additionally, disturbances to forest cover are expected to be higher near water bodies as wild animals and human beings need regular access to water. Streams, rivers, lakes and artificial water holes are mapped as water.

Space as a variable: We used the term “space” to denote zones/sub-regions in our larger study area to capture the optimal spatial entity at which the explanatory variables are driving the land cover change. The two dummy variables mentioned below are our space variables. Presence or absence of protected area: This is a dummy variable where presence of park is denoted by 1 and absence by 0. This variable is important for our deforestation model, as protection to forest should discourage deforestation inside the park. This variable is not used for our tree plantation gain model as our tree plantation gain is mostly located outside the park. Even though in the future there could be an increase in tree plantations inside the park, the driving factors that would be operating inside the park would be different from the driving factors outside the park.

Higher Plantation Gain Area: This is a dummy variable (DummyWest) (Figure 3) where the study area is divided into two different zones of East and West after examining the spatial distribution of residuals of the main effect explanatory variable models. Having higher tree plantation gain in the western part of the study area as compared to the eastern part justified this dummy variable creation. The digital form of only these variables, i.e., elevation, slope, distance to forest edge, road, water, Bangalore, towns and villages and presence or absence of protected areas, were available. Precipitation although an important variable for modeling LULCC is not considered for this study because of the small area of our study site, therefore, not much spatial variability of precipitation pattern is expected. Other variables such as population, land tenure, soil moisture is important to model deforestation and tree plantation gain in the present landscape but are not available in spatial form and should be considered as a data limitation of our study.

3.2. Sampling Procedure

For our present study, a data sample was selected from our spatially explicit LULCC (dependent) and proximate (independent) variables. Random sampling procedure was used to select our observations due to the presence of spatial auto-correlation in the data. The sample included 10% of the total study area for each model, which resulted in approximately 55,000 to 65,000 observations for tree plantation and deforestation models respectively, with an unequal number of 0 (no change) and 1 (change) observations of the dependent variables (

Table 2. Samples.

Sample Pixels	Deforestation (1973–1992)	Deforestation (1992–2007)	Plantation Gain (1973–1992)	Plantation Gain (1992–2007)
Absent	30,642	30,861	n/a	23,129
Present	34,826	34,531	n/a	34,516
Total	65,468	65,392	n/a	57,645

Absent and present represents the two strata of “no change area” and “change area”, respectively, from where samples were collected. n/a = not applicable.

). Unequal sampling of the dependent variables do not affect the coefficient estimates of the explanatory variables in logit models but only affect the constant term [5,45]. The entire data sample was imported into SPSS statistical software (IBM Analytics, Armonk, NY, USA) for our logistic regression analysis (Statistical Package for the Social Sciences, 2010).

Continuous independent variables that did not show linear behavior were transformed using either a logarithmic or a square root function. All independent variables were standardized using the formula below, prior to running the logistic regression analysis.

$$Z=\frac{X-\overline{X}}{S}$$

where Z is the standardized variable, X is the value of the original variable, $\overline{X}$ is the mean, and S is the standard deviation. The unstandardized logit coefficients measure the absolute contribution of each variable in determining the probability of occurrence of an event and thus could be misleading to interpret as a unit change in a variable as this is not equal from variable to variable. As such, the independent variables were standardized to reduce the disparities in scale of measurement and variance because of different units of the variables [13]. This process ensured that our logit coefficients and odd ratios are standardized, thus, it is easier to evaluate the contribution of each independent variable in our models.

3.3. Test of Multicollinearity

The independent variables in our models were tested for multicollinearity, as it is a prerequisite of any statistical method for the explanatory variable to be independent of each other so that the importance of each explanatory variable could be ascertained individually. We tested for Pearson’s correlations among all independent variables (

Table A1. Pearson Correlation Coefficients.

		FORDIST73	FORDIST92	ALLRDDIST	MJRDDIST	MIRDDIST	VILLDIST	DWATER	MJTOWNDIST	BANGDIST	SLOPE	ELEV	PRO_UNPRO	DummyWest
	FORDIST73	1.000	0.497	−0.032	−0.127	0.067	−0.031	0.020	−0.066	0.085	−0.201	−0.184	−0.173	0.043
	FORDIST92	0.497	1.000	−0.083	−0.189	−0.032	−0.151	−0.128	−0.137	−0.023	−0.234	−0.157	−0.417	0.098
	ALLRDDIST2	−0.032	−0.083	1.000	0.299	0.280	0.416	0.224	−0.112	0.316	0.205	−0.011	0.161	−0.027
	MJRDDIST2	−0.127	−0.189	0.299	1.000	−0.111	0.370	−0.022	0.039	0.213	0.263	−0.120	0.377	−0.061
	MIRDDIST2	0.067	−0.032	0.280	−0.111	1.000	0.296	0.677	0.036	0.227	0.042	0.039	−0.012	0.008
	VILLDIST	−0.031	−0.151	0.416	0.370	0.296	1.000	0.179	0.044	0.156	0.169	−0.152	0.245	−0.033
	STWDIST	0.020	−0.128	0.224	−0.022	0.677	0.179	1.000	0.158	0.319	0.084	0.276	0.175	−0.016
	MJTOWNDIST	−0.066	−0.137	−0.112	0.039	0.036	0.044	0.158	1.000	−0.080	0.088	−0.244	0.221	−0.043
	BANGDIST	0.085	−0.023	0.316	0.213	0.227	0.156	0.319	−0.080	1.000	0.067	−0.038	0.125	0.042
	SLOPE	−0.201	−0.234	0.205	0.263	0.042	0.169	0.084	0.088	0.067	1.000	0.050	0.255	−0.059
	ELEV	−0.184	−0.157	−0.011	−0.120	0.039	−0.152	0.276	−0.244	−0.038	0.050	1.000	0.083	−0.047
	PRO_UNPRO	−0.173	−0.417	0.161	0.377	−0.012	0.245	0.175	0.221	0.125	0.255	0.083	1.000	−0.114
	Dumy_P7392	0.043	0.098	−0.027	−0.061	0.008	−0.033	−0.016	−0.043	0.042	−0.059	−0.047	−0.114	1.000

Correlation is significant at the 0.01 level (2-tailed).

). A critical value of 0.80 was used to eliminate variables from our model [46]. All independent variables were used in our logistic regression analysis as the coefficients were all below 0.50 with the exception of distance to water and minor roads. We further checked the Variance Inflation Factor (VIF) and tolerance value to diagnose collinearity. A critical VIF value of 10 and above [47] and a tolerance value of less than 0.1 [46] is considered a cause of concern. None of the variables showed a VIF of higher than 2.5 and tolerance value lower than 0.1.

3.4. Logistic Regression Model

The relationship between the LULCC (forest loss and tree plantation gain) and environmental and proximate variables were tested using logistic regression analysis. For a logistic model, the dependent variable is categorical (binary in our case) presence or absence of any event. In this study, forest loss = 1 and other = 0 and tree plantation = 1 and other = 0 for the periods 1973–1992 and 1992–2007. The independent variable can be either continuous or categorical. For the present study, all the independent variables are continuous layers except for the presence or absence of protected area. The binary logistic technique produces coefficients for each independent variable based on a sample of data.

The logistic function gives the probability values that can be quantitatively expressed as a function of the explanatory variables stated in the following form:

$$p=E(Y)\frac{\exp (\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2+\dots {\beta }_n{X}_n)}{1+\exp (\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2+\dots {\beta }_n{X}_n)}$$

where p is the probability of forest loss or tree plantation gain in a pixel, E(Y) is the expected value of the binary dependent variable Y, α is the intercept, β_n a coefficient to be estimated for each independent variable X_n. The logistic function can be transformed into a linear response with the logit transformation.

$$\log (p)=\log \left[\frac{p}{1-p}\right]=\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2+\dots {\beta }_n{X}_n$$

Odd ratios are used to enable interpretation of logistic regression models [46]. Odds ratio is an indicator of the change in odds resulting from a change in one unit of the predictor and thus is a measure of association of how much more likely (or unlikely) it is for an outcome to be present for a set of values of independent variables [48]. Estimated odds values are computed as the exponential of the parameter estimate values [48,49]. The probability, logit and the odds are three different ways of expressing the same things [46].

$$\mathit{odds}\left(p\right)=\exp (\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2+\dots {\beta }_n{X}_n)$$

Model predictability was assessed using the model Wald chi-square statistic (Wald χ²), log likelihood statistics, parameter estimates (B), estimates of standard errors of the parameter estimates (SE), and the significance of the probabilities for parameter estimates. Wald χ² indicates if the parameter estimates of the independent variable is significantly different from zero and thus indicates the predictive power of each independent variable in the prediction of an event [47]. The log likelihood statistics provided an estimate of the variance unaccounted for after the model had been fitted and thus, larger values of log likelihood statistics indicates a poorly fitted model [47]. Parameter estimates (B) provided by the models were used to measure the association of independent variables with deforestation and tree plantation gain. A positive value of the parameter estimate shows that the likelihood of occurrence of an event increases as the predictor variable increases. Similarly, a negative value of the parameter estimate shows that an increase in predictor variables will decrease the likelihood of occurrence of an event. A variable was selected if it was statistically significant at p < 0.05. The goodness of fit of a logistic model is tested using Nagelkerke’s R² values in SPSS, which is an adjusted Cox and Snell coefficient and varies between 0 and 1 [50]. This R² is also called as pseudo R² as it is computed differently from a regular R² in a linear regression analysis [39]. A pseudo R² is interpreted similar to regular R² where a Nagelkerke’s R² value of 1 would indicate that the model predicts the outcome variable perfectly. Nonetheless, unlike the conventional regression analysis, a Nagelkerke R² of 0.2 and above is considered as a good fit of a logistic model because of the binary response variable [5,51]. Classification accuracy also indicates model fitness, with 100% indicating a perfect model [19]. Overall model fitness is measured by a significant Wald Chi-square value, Nagelkerke’s R² of 0.2 and above and a high correct classification percentage. Lastly, variables and their parameter estimates were also assessed to see if they made sense ecologically. The presence of outliers was measured using Cook’s distance, leverage and standardized residuals. Cases that had a standardized residual value more than 2 were excluded from the model [47].

In this study, several runs of backward stepwise logistic regression procedure were performed in SPSS 18 to select the best set of predictor variables (Figure 4), a widely used approach in land change modeling [52,53]. Independent variables were entered individually to assess the bivariate logistic relationship, entered all simultaneously, and entered in various combinations to assess predictability of the models. All these steps were run to maximize the predictability of the model. All the models were rerun after excluding cases which had a standardized residual value of more than 2. Spatial autocorrelation was run on the rest of the models residual by calculating Moran’s I in ArcGIS. Presence of spatial autocorrelation results in inefficient parameter estimates and inaccurate measures of statistical significance [54]. The results for all the models indicated the presence of clustering in our residuals. Since all the models indicated spatial autocorrelation, cluster and outlier analysis using Anselin Local Moran’s I was calculated using ArcGIS to check for the spatial distribution of the clustering [55]. Output of this analysis is a Local Moran’s I index, Z-score, p-value and cluster type code (COType). A high positive Z score for a feature indicates that the surrounding features have similar values (either high values or low value). The COType field indicates HH for a statistically significant (0.05 level) cluster of high values and LL for a statistically significant (0.05 level) cluster of low values. A low negative Z score for a feature indicates a statistically significant (0.05 level) spatial outlier. The COType field indicates if the feature has a high value and is surrounded by features with low values (HL) or if the feature has a low value and is surrounded by features with high values (LH). Dummy variables were created and were interacted with our main explanatory variables after careful study of the residual clusters (

Table 3. Logistic regression parameters estimated from deforestation.

Variables	B	S.E.	Wald	Sig.	Exp(B)
1 Deforestation Model 1973–1992
Constant	−0.555	0.020	775.018	0.000	0.574
DummyPA	−0.628	0.026	600.853	0.000	0.534
SLOPE	0.583	0.021	806.010	0.000	1.791
ELEV	0.864	0.018	2386.039	0.000	2.373
FORDIST73	−3.862	0.037	10,639.472	0.000	0.021
ZVILLDIST	−0.073	0.014	27.842	0.000	0.929
MIRDDIST	−0.067	0.018	14.768	0.000	0.935
MJRDDIST	0.157	0.013	145.453	0.000	1.170
ALLRDDIST	0.099	0.013	55.421	0.000	1.104
WATERDIST	−0.482	0.019	659.504	0.000	0.617
DummyPA × SLOPE	−0.775	0.025	968.927	0.000	0.461
DummyPA × ELEV	−0.421	0.029	214.057	0.000	0.656
DummyPA × MIRDDIST	−0.384	0.031	155.111	0.000	0.681
DummyPA × WATERDIST	0.613	0.028	474.414	0.000	1.846
2 Deforestation Model 1992–2007
Constant	−3.796	0.053	5136.465	0.000	0.022
DummyPA	−0.082	0.031	6.915	0.009	0.922
SLOPE	0.257	0.028	84.626	0.000	1.293
ELEV	0.138	0.025	31.656	0.000	1.148
FORDIST92	−11.209	0.110	10,439.912	0.000	1.6 × 10−5
VILLAGEDIST	−0.203	0.017	145.991	0.000	0.817
MIRDDIST	0.084	0.024	12.623	0.000	1.088

¹ Deforestation Model 1973–1992: Significant at 95% confidence level, R² = 0.602, Classification Accuracy = 80.2%. ² Deforestation Model 1992–2007: Significant at 95% confidence level, R² = 0.761, Classification Accuracy = 87%.

and

Table 4. Logistic regression parameters estimated from tree plantation gain.

Variables	B	S.E.	Wald	Sig.	Exp(B)
1 Tree Plantation Gain Model 1992–2007
Constant	0.061	0.020	8.932	0.003	1.063
SLOPE	−0.119	0.017	48.191	0.000	0.887
ELEV	−0.411	0.029	201.656	0.000	0.663
VILLDIST	−0.230	0.024	90.603	0.000	0.794
WATERDIST	−0.058	0.018	10.410	0.001	0.944
MIRDIST2	−0.329	0.020	279.259	0.000	0.720
MJRDIST2	−0.175	0.023	60.370	0.000	0.839
ALLRDDIST2	−0.073	0.012	40.599	0.000	0.929
MJTOWNDIST	−0.182	0.031	35.083	0.000	0.834
BANGDIST	−0.265	0.016	270.242	0.000	0.767
DUMMYWEST × SLOPE	−0.070	0.021	11.345	0.001	0.932
DUMMYWEST × ELEV	−1.039	0.056	342.470	0.000	0.354
DUMMYWEST × VILLDIST	−0.191	0.029	43.586	0.000	0.826
DUMMYWEST × WATERDIST	−0.070	0.024	8.745	0.003	0.932
DUMMYWEST × MIRDIST	0.805	0.027	876.347	0.000	2.238
DUMMYWEST × MJRDIST	0.169	0.025	46.284	0.000	1.184
DUMMYWEST × MJTOWNDIST	0.437	0.034	163.255	0.000	1.548
DUMMYWEST × BANGDIST	−0.288	0.030	90.645	0.000	0.749

¹ Tree Plantation Gain Model 1992–2007: Significant at 95% level, R² = 0.223, Classification Accuracy = 68.7%.

). Dummy variable DummyPA (Presence of protected area = 1 and absence of protected area = 0) were interacted with other independent variables for deforestation models and DummyWest (Dominant tree plantation in west of protected area = 1 and east of protected are = 0) were interacted with main effect variables in tree plantation model (Table 3). Specifics of dummy interaction variables for each model are explained in the Result Section 4.1 and Section 4.2. Our exploratory models were rerun with these interaction variables to assess each variable’s strength and significance and overall predictability of our models. Our models with interaction variables were compared with our main effect explanatory variable models. Final models of deforestation and tree plantation gain were run after careful study of the variables that were significant. Spatial distribution of residual clusters was calculated once again using Anselin Local Moran’s I on our remaining residual. The Nagelkerke R² and classification accuracy indicated our model’s overall performance.

4. Results

4.1. Deforestation Model

Two models were created that described the relationship between the explanatory variables and deforestation occurring between 1973–1992 and 1992–2007. All the deforestation models had a Nagelkerke R² more than our critical value of 0.2 and thus showed a good model fit. For the first deforestation model, interaction between DummyPA and elevation, slope, distance to minor roads, major roads, all roads and water were included for time period one after looking at the spatial distribution of the residuals. Distance to villages was excluded as there are only seven villages inside the park and most of the high values of residual clustering (under prediction) were away from the villages. For our second deforestation model, interaction between DummyPA and slope, elevation and minor roads were included in the model. The remainder of the interaction terms was not relevant for this time period. Our DummyPA variable also tests the effectiveness of protected area status on deforestation as compared to the overall landscape. Including the interaction variables in our model improved the model fit indicated by both an increase in Nagelkerke R² and classification accuracy of both time period deforestation models. The signs of estimated coefficients of most of the explanatory variables made ecological sense for this landscape as a result of inclusion of the dummy interaction variables in both the deforestation models and thus only the full models with interaction effects have been presented. However, for the second deforestation model, there were some signs of estimated coefficients for predictor variables which were not as expected. Both the models of deforestation were statistically significant at the p = 0.005 level.

Table 3 shows the results of the logistic regression model analysis of deforestation. The relationship between the explanatory variables and deforestation vary between the two time periods (Table 3 and

Table 5. Logistic regression signs of parameters estimated.

Independent Variables	Deforestation	Tree Plantation Gain
	1973–1992	1992–2007	1992–2007
	R2 = 0.602	R2 = 0.761	R2 = 0.223
SLOPE	+	+	−
ELEV	+	+	−
ALLRDDIST	+	+	−
MIRDDIST	−	+	−
VILLDIST	−	−	−
WATERDIST	−	−	−
MJRDDIST	+	n/s	−
FORDIST73	−	−	n/a
FORDIST92	n/a	−	n/a
TOWNDIST	n/a	n/a	−
BANGIDST	n/a	n/a	−
DummyPA	−	−	n/a
DummyPA × SLOPE	−	−	n/a
DummyPA × ELEV	−	−	n/a
DummyPA × MIRRDDIST	−	+	n/a
DummyPA × WATERDIST	+	n/a	n/a
DummyWest × SLOPE	n/a	n/a	−
DummyWest × ELEV	n/a	n/a	−
DummyWest × VILLDIST	n/a	n/a	−
DummyWest × WATERDIST	n/a	n/a	−
DummyWest × MIRDIST	n/a	n/a	+
DummyWest × ZMJRDIST	n/a	n/a	−
DummyWest × ALLRDDIST	n/a	n/a	n/s
DummyWest × MJTOWNDIST	n/a	n/a	−
DummyWest × BANGDIST	n/a	n/a	−

All predictors are significant unless indicated. n/s = not significant (was dropped from Backward LR model), n/a = not applicable.

). In the overall landscape in first time period, distance to forest edge (−) and protected area status (−) are the dominant drivers of deforestation, whereas, in second time period, distance to forest edge (−) and all roads (+) are the dominant drivers. Proximity to water bodies (−) also encourages deforestation in both time periods in overall landscape. Models also show a distinction between the driving factors inside and outside the park. Deforestation in both the eras in overall landscape has a positive relationship with elevation and slope. However, inside the protected area, deforestation is on higher elevation (+) and gentler slope (−) in first time period and lower elevation (−) and gentler slope (−) in the second time period. Further, proximity to minor roads is an important determinant of deforestation inside the park. Deforestation occurred closer to minor roads in first time period and farther away from minor roads in second time period, which is logical, as, in the previous time period, forests have already been cleared from areas closer to minor roads.

4.2. Tree Plantation Gain Model

The relationship between tree plantation gain and various explanatory variables were explained by the two models representing two periods (1973–1992 and 1992–2007). DummyWest was interacted with the main effect explanatory variables for both the models. Inclusion of DummyWest explained the explanatory variables dominant in the western part of our study site. Inclusion of DummyWest also increased our Negelkerke R² and classification accuracy.

The models for tree plantation gain have a lower explanatory power than the deforestation models (Table 4). The first model (1973–1992) for tree plantation gain is not significant (Nagelkerke R² = 0.045) and also has a low classification accuracy (60.2%) (

Table A2. Logistic regression parameters estimated from tree plantation gain.

Variables	B	S.E.	Wald	Sig.	Exp(B)
Tree Plantation Gain 1973–1992
Constant	0.566	0.022	662.487	0.000	1.762
ELEV	−0.378	0.032	140.957	0.000	0.685
VILLAGEDIST	−0.086	0.018	21.905	0.000	0.918
WATERDIST	0.115	0.018	42.848	0.000	1.122
MJTOWNDIST	−0.287	0.026	118.837	0.000	0.751
BANGDIST	0.477	0.015	974.628	0.000	1.612
MIRDIST	−0.052	0.014	13.843	0.000	0.950
MJRDIST	−0.265	0.023	130.639	0.000	0.767
ALLRDDIST	0.046	0.011	18.137	0.000	1.047
DummyWest × SLOPE	0.079	0.019	16.599	0.000	1.082
DummyWest × ELEV	0.796	0.055	213.047	0.000	2.218
DummyWest × VILLAGEDIST	0.155	0.023	46.491	0.000	1.167
DummyWest × MJTOWNDIST	0.300	0.029	105.724	0.000	1.350
DummyWest × MIRDIST	0.050	0.022	5.099	0.024	1.051
DummyWest × MJRDIST	0.216	0.025	73.707	0.000	1.240
DummyWest × WATERDIST	−0.359	0.022	262.580	0.000	0.698
DummyWest × BANGDIST	−0.145	0.027	29.215	0.000	0.865

Significant at 95% level, R² = 0.045, Classification Accuracy = 60.2% (1973–1992).

), however, the second model (1992–2007) has a reasonably good fit (Nagelkerke R² = 0.223, classification accuracy = 68.7%) (Table 4). The signs of the coefficient of estimates make ecological sense for the second model for tree plantation gain for explanatory variables. In the overall landscape elevation (−), distance to minor roads (−), distance to villages (−) and distance to Bangalore (−) play a significant role in explaining tree plantation gain. Tree plantation gain is likely to occur in lower elevation, closer to major and minor roads, villages and Bangalore. Tree plantation gain is also most likely to occur on flatter slopes (−). In the western part of our study area, as our DummyWest interaction variables indicate, tree plantation gain occurred on flatter slope, lower elevation, and closer to villages. Minor roads showed a positive relationship in this sub-region which in ecologically unexpected. However, a positive relationship between tree plantation gain and major towns in the sub-region is logical, as most of the towns are located on the eastern part of the study site.

4.3. Residual Analysis of Models

Cluster and outlier analysis for the present study was done using the Anselin Local Moran’s I module of ArcGIS 9.3 (ESRI, San Diego, CA, USA). The spatial distributions of residuals of deforestation and tree plantation gain models for both the time periods show a high residual inside the national park even after the inclusion of interaction variables between DummyPA and relevant explanatory variables. This suggests missing explanatory variables from our deforestation and tree plantation models, which could be spatially structured and thus causing the residual clusters. It is also possible that some of our explanatory variables are spatially interacting with each other.

5. Discussion

Spatial-statistical tools can be very powerful tools to examine the relationship between LULCC and its proximate causes and infer about the underlying driving factors. Ideally, spatial-statistical tools provide a means to quantify the impact of various biophysical, social and economic driving factors (explanatory variables) by mean of analyzing proxy variables affecting LULCC (outcome variable). Spatial statistics is useful in identifying the type (sign) of relationship and ranking the relative importance of these explanatory variables by estimating the magnitude (value) of effect brought by each of them on LULCC. These statistical estimates can further be used to predict future land-cover changes such as forest cover change [25,56,57,58].

The relationship between land-cover change and its various explanatory variables are affected by geographical scale of the study, different sub-regions within the overall system and different time periods considered [5,25,26,39,54,59,60]. Although our study is done at one scale (landscape level), different zones and eras are considered in explaining the relationship between explanatory variables and land-cover change. Our results show differences in the processes affecting deforestation between 1973–1992 and 1992–2007 and within protected and non-protected areas. For tree plantation gain, a bad model fit (1973–1992) did not allow us to compare between different time periods. However, as a result of the presence of residual clustering, we divided the study area between eastern and western zones to check for variation in relationship between explanatory variables and tree plantation gain in different sub-regions.

In both eras, we found distance to forest edge (indication of previous deforestation) to be the most important explanatory variable explaining deforestation. Prior studies have shown the spread effect of deforestation, where deforestation in one period occurs in proximity to previously deforested areas [25,61,62]. Protection to forest cover is also an important determinant of deforestation in both time periods, where areas inside BNP, Talli and Bilikal Reserve Forest were less likely to be deforested shown though the negative correlation between DummyPA and deforestation. Similar studies analyzing the relationship between various explanatory variables and deforestation found a negative correlation between protected areas and deforestation [26,54]. Globally, policy makers heavily depend on establishing protected areas to conserve forest resources [63,64], even though effectiveness of protected areas in conserving biodiversity has been questioned [44]. However, in India, it is generally agreed that protected areas are working in reducing deforestation [44,65]. This negative correlation between deforestation and BNP, Talli and Bilikal Reserve Forest throws light on the underlying driving factor of policy measures at the national level (establishing protected areas) and its impact on local level land-cover changes, here deforestation.

The relationship between topography and deforestation in our overall landscape is contradictory to many other studies that have found slope and elevation to be negatively correlated with deforestation. Lack of a definite relationship or a positive relationship between deforestation and topography is not uncommon as found by several studies [19,66]. In both time periods, the relationship of deforestation with elevation (+) and slope (+) (Table 3) is logical as most of the forest cover in this area is located in higher elevation and steeper slope areas (Figure 2). Further, the largest patch of deforestation in the first time period is centralized in the east central part of BNP which is comparatively in higher elevation, which explains the positive relationship of topography with deforestation. Conversely, the negative relationship of topography and deforestation inside our protected area is on par with the general idea that deforestation occurs in flatter slopes and lower elevation because of easy accessibility of these regions for illegal loggers and encroachers.

Deforestation occurring in both higher elevation and steeper slopes as well as in lower elevation and flatter slopes brings our attention towards two interlinked ideas. First, deforestation is governed by different driving factors inside and outside the protected area. Secondly, deforestation in the east central part of the study area, most of which is outside BNP may not be meditated by topography. Something else that would have overcome limitations posed by elevation and slope may have caused deforestation in this region. The most probable cause, as supported by our field interviews, is people–park conflict, occurring because of unresolved land ownership. Topography, even though it is not a constraint for deforestation, was found to be the most important determinant for tree plantation gain in the overall landscape as well as for the western zone. It is clear from this that land uses such as agriculture are strongly meditated by topography (see Section 4.1).

Looking at the results, proximity variables in BNP and its surroundings behave largely as expected for both deforestation and tree plantation gain models, but we find diversity in direction and magnitude of effect. Proximity to water is an important determinant of deforestation in this region as deforestation occurred closer to water bodies. Distance to water, however, plays a much more important role in explaining deforestation in the first than in the second time period. Water bodies, especially, rivers provide access to loggers, cultivators and encroachers [34]. In this landscape, water bodies are streams, lakes and ponds, which may not provide access but attract encroachers and illegal loggers as these are the main sources of water for the local village communities. Relatively small effect of proximity to water on tree plantations gain may be explained by the fact that these tree plantations (eucalyptus and coconut palm) do not entirely depend on direct irrigation as opposed to other field crops. However, finding tree plantations closer to the water bodies may be explained by the cultivation of field crops near and under coconut palm plantations (field observations).

Roads in general, and rural roads in particular, facilitate deforestation by opening the forested areas to loggers, encroachers and agriculturalists [34,36,67,68]. Thus, the relationship of deforestation with minor roads (−) and villages (−) in the first time period is likely a result of forest access that the minor roads provided to villagers. Major roads are not a significant indicator of forest clearing both the eras. As our previous study suggested, most of the deforestation in the region occurred as a result of domestic fuel gathering, animal grazing and illegal sand mining rather than large scale logging activity that would require highways for transportation of woody products [3]. This could explain why major roads such as highways are not an important determinant of deforestation in the region. We, however, find all roads to have significant relationship with deforestation (+) in the second time period. Overall connectivity and better flow of network provided by all roads might have been important in explaining deforestation in the second time period.

Negative relationship between tree plantation gain and town, Bangalore, villages, major and minor roads in the overall landscape is expected and shows the importance of accessibility to market in explaining agricultural expansion. Cash crop cultivation is promoted by access to road which connects them to urban markets [66]. Bangalore as one of the fastest growing cities in India, acts as a major market center for agricultural products from surrounding areas especially cash crops like coconut palm and eucalyptus. Looking at the beta coefficients, distance to Bangalore is more important than the distance to villages for tree plantations. As tree plantations such as eucalyptus and coconut palm trees are considered easy maintenance crops, it reduces the need for farmers to travel between their home and field. However, minor roads are more important than major roads as these are bulky agricultural products and minor roads provide access to the major roads for these crops transportation to the market.

Even though, in the overall landscape, the negative relationship of these distance variables with tree plantations makes ecological sense, within the western zone, a positive relationship of tree plantation gain with minor roads is unexpected. Further, all roads in the western zone and our accessibility to town variable do not show a significant relationship with tree planation. We would expect all roads to significantly affect our tree plantations in the western zone. These discrepancies could be because we have considered distance to different types of roads, towns, villages as individual entities and not part of an overall system of network. A majority of spatial regression model studies of LULCC define accessibility as a measure of only straight-line distance from roads or markets [25,26,39,40,68,69]. However, the importance of road should vary based on a combination of how connected or inter-connected any individual road is, how these roads connect different types of settlements to market, how the connectivity of roads change as a result of new road development, and whether these roads are access roads or major roads. It is thus important to remember that roads and various populations centers interact spatially. Few recent studies have used advanced measures of accessibility and found them to better represent their landscape than rather just the use of distance in kilometer terms [69]. An improved measure of accessibility that includes various network indices that describes the importance of roads would be a better representation of accessibility in a LULCC model.

Spatial regression models are widely used models in LULCC in explaining a variety of LULCC processes particularly deforestation and agricultural expansion. Our models provided substantial information about the relationship between various proximate causes and LULCC in this region. Additionally, we can also infer about some of the underlying root causes of LULCC and relate it to some well-established theoretical frameworks by examining the proximate causes as has been done by prior studies [5,13,69]. For example, our DummyPA explanatory variable highlights the importance of protected area establishment and its impact on local level LULCC. Additionally, deforestation not constrained by topography outside BNP may simply be an artifact of presence of forest cover in a densely populated area “outside” the PA boundary. However, inclusion of land tenure variables may be necessary to infer about driving forces such as land-use policies. In addition, the relationship of explanatory variables with the tree plantations may suggest a von Thünen like model where tree plantations, a bulky agroforestry product in this region, may be controlled by distance to market, a proxy to transportation cost [70]. Even though our study includes spatial variability of different proximate and biophysical causes, the present study, specifically our tree plantation models may be weak in predicting future LULCC as it does not include temporal variations in our explanatory variables. Including additional explanatory variables such as dynamic road network and changes in their connectivity, changes in importance of individual roads as a result of new emerging roads, changes in population density, changing land tenure regimes, and new emerging population centers in future research endeavors would be essential for predicting future LULCC. Finally, we tried to address spatial autocorrelation in our model by including dummies representing different parts of our landscape. Although this improved our model fit and provided coefficient estimates which made ecological sense, it did not remove our spatial autocorrelation altogether from our residuals.

6. Conclusions

This study integrates remote sensing and spatially explicit data to develop a statistical model for LULCC in Bannerghatta National Park (BNP) and its surroundings. Relationships between deforestation and tree plantation gain and their proximate and biophysical driving factors were quantified using a logistic regression model. Our results reaffirms the importance of protected area establishment in managing forest resources in densely populated regions and highlights the importance of including past land-cover changes (previously deforested areas) to explain present land-cover changes. We found proximity to roads, towns, villages and water sources as one of the important determinants of plantation driven reforestation in the region. The results also showed contradictory signs of the estimated parameters (not as expected in this landscape) emphasizing the need to incorporate weighted measures of accessibility for LULCC model building. Running models for different eras (1973–1992 and 1992–2007) rather than a single timespan (1973–2007) provided significant insights into the varying driving factors affecting the LULCC for different eras. The study also incorporated spatial variability of different proximate and biophysical causes to account for spatial autocorrelation by creating sub-regions as space variables such as the protected versus unprotected area and the DummyWest zone. The different directions and magnitude of effects of these proximity variables in different sub-regions provided us with valuable insights on the different driving factors for different sub-regions, for example, how the processes that affect deforestation are different inside and outside the park. Consequently, this research encourages future work on LULCC modeling using space as an explanatory variable to account for spatial autocorrelation.