Spatial modeling on land use change in regional scale of Java Island based-on biophysical characteristics

. This study discusses a biophysical-based spatial modeling for land use change in Java Island considering neighborhood interactions between land use types and the change area. These neighborhood characteristics used in logistic regression model to estimate the probability of the change events occurrence. Moreover, the future role of land use change is then projected using the Markov model based on the annual land use changes map. The results indicate that paddy rice with irrigation system (double cropping), especially in upland areas has a high positive spatial autocorrelation with the change areas. Residential area, paddy rice, and upland with intensive cropping have a high effect to the probability of change occurrence. Meanwhile, barren lands/dry land, bush-shrub and mixed garden give a negative correlation to the change occurrences in agricultural lands. In the case of forestland, the results show some land use types such as upland with intensive cropping and plantation have positive contribution to the change of land. The accuracy of the model was also assessed through comparison the projection with the actual area. The results indicate that the future role of each land use type is different based on the trend period predictor in the model.


INTRODUCTION
Understanding the role of land-use in regional and global environmental changes requires a historical reconstruction of past change conversions and projection of future trends (Lambin, 1997).Land use changes are the result of the complex interactions between human and biophysical-environment, which affect a wide range of temporal and spatial scales.Such changes can be measured and analyzed, including determination of their driving factors, however, it is difficult to aggregate these changes regionally and globally.On these scales, the aggregation level of the data is related to the simplification of the variability over geographic and socioeconomic conditions, and it might cause inaccuracy in measuring or estimating the effects of land use change (Verburg et al., 2004a).On the other hand, developing the spatial model which can identify critical areas of land use changes and provide insight on the changes pattern is a challenge in the modeling of regional and global land use changes.
Land use changes are driven by many and diversified underlying processes (Geist and Lambin, 2002;Firman, 2004;Turner et al., 2007).Anthropogenic drivers (demographic, technology and socio-economic conditions) and biophysical constraints such as soil, climate and topography, determine the spatial pattern of land use (Skole and Tucker, 1993;Turner et al., 1994), and their relationships are often used in the spatial land use models to allocate land use changes (Verburg et al., 1999).The spatial integration of socio-economic and biophysical-environmental data, which provides a mechanism to explore the relationship between human activities, landscape conditions and land use changes, is an important step to improve understanding of the land use change and its future role (Evans and Moran, 2002).However, in more complex land use patterns, such as in Java Island, land use change is the result of many, non-linear, interactions between socio-economic and cultural conditions, biophysical constraints and land use history.Such level of complexity has limited the scope of spatial modeling in this study.
Many of the biophysical processes related to land use change can be represented as interactions between different land use types, their biophysical conditions and the spatial elements of the landscape (Evans and Moran, 2002).These interactions need to be considered in order to understand the role of land use change on a regional scale.This study will explore the relationship between change in land use and the biophysical aspects of land, which will facilitate further analysis and understanding of the dominant process of land use change allocation.
This paper demonstrates a spatial analysis to investigate the possibility that land use changes in one area may be influenced by a land use type and/or other land use changes in surrounding areas.In other words, the spatial dependency will be considered to be the functional relationship between what happens with one land use type in space and what happens elsewhere (Anselin, 1988).Furthermore, considering the integrated analysis of the change area, neighborhood characteristics, and biophysical environmental factors, a spatial model will be developed in order to understand some important facts and trends related to land use change on Java Island.

Overview of Existing Land Use Models
A large number of studies developed many models to simulate the trends of land use change (Lambin, 1994;Himiyama, 1999).Differences among those modeling techniques often relate to differences in the purpose or the scale of study.Some explorative models were developed to design alternatives for present land use so as to optimize land use allocation based on the socio-economic and biophysical environment (Stoorvogel, 1995;Schipper, 1996).Other research explored possible changes of land use in the near future as a function of many driving factors (Costanza et al., 1990;Hall et al., 1995).
Meanwhile, in a global and/or regional scale, many of the model studies are regarded as development of global X models, such as the General Circulation Models (GCM), Global Vegetation Model (GVM) or IMAGE model.The IMAGE model is a system dynamic model considering both social and physical aspects of the world system (Alcamo et al., 1996).This model calculates the allocation of land use changes in the near future based on setting assumptions on population, economy and economic activities, and then computes future changes in consumption of energy, food and timber.Himiyama (1994) mentioned that these global models usually rely heavily on statistics, as well as on complicated mathematical functions, but have difficulty in quantitatively testing the results.Then, Robinson et al. (1994) proposed basic concepts of global land use/cover model (GLCM), in order to promote a global model for LULCC.
Moreover, a dynamic land use model called "Conversion of Land Use and its Effect (CLUE) model" was developed to simulate land use change in Costa Rica (Veldkamp and Fresco, 1996).It simulated land use change in space and time as a result of interacting biophysical and human drivers.The land use changes allocated by considering the complex interactions between historic and present land use, socio-economic conditions and biophysical constraints.The CLUE model was applied to explore land use changes for a scenario of further urbanization in Java (Verburg et al., 1999).The results of this study show that the land use change will especially occur in the lowland areas, either directly through construction or indirectly through the demand for higher value crops, meanwhile, the upland areas will stay primarily rural for the considered scenario.The most obvious pattern is the decrease in paddy rice fields due to increases in housing areas in the northern coastal plain of Java and an increase in estate crops in many upland areas.Moreover, the model also predicted some land use changes in presently undeveloped areas in the western and southern part of Java.

Land use patterns and neighborhood characteristics
In analogy with the theories used in landscape ecology (Forman and Godron, 1986), the pattern of land use change might be related with the interactions between the spatial elements of the landscape (Verburg et al., 2004b).For example, the influence of a residential area on the conversion of agricultural land uses.

Spatial Dependency and Spatial Weight Metric
In this study, the term of spatial dependency (spatial autocorrelation) can be considered to be a functional relationship between what happens at one point in space and what happens elsewhere (Anselin, 1988).It gives a means to the importance of distance and space, as mentioned by the Tobler's first law of geography that "everything is related to everything else, but near things are more related than distant things" (Tobler, 1970).Cliff and Ord (1973) present a general approach to express the interaction between two spatial areal units by using a combination of distance measures (inverse distance, or negative exponentials of distance) and a measure of the length of their common border.The function is as follows: where dij stands for the distance between spatial unit i and j, βij denotes the proportion of the interior boundary of unit i in contact with unit j,  and  are parameters.Spatial areal units are typically suited to have their interaction expressed by equation ( 1), both the distance between their centers and the relative importance of their common border are taken into account.The specific process through which spatial statistics integrates space and spatial data into the mathematics function above is explained by a spatial weight matrix (Getis and Aldstadt, 2004).It is a quantification of the spatial relationships that exist among the features in a data set.
Conceptually the spatial weights matrix is an N × N table (N is the number of features in the data set).There is one row for every feature and one column for every feature.The cell value for any given row/column combination is the weight that quantifies the spatial relationship between those row and column features.At the most basic level, there are two strategies for creating weights to quantify the relationships among data features: binary or variable weighting.By the binary strategies (such as fixed distance, K-nearest neighbors, and contiguity), a feature is either a neighbor (1) or it is not (0).Meanwhile, the weighted strategies (inverse distance or zone of indifference), neighboring features have a varying amount of impact and weights that are computed to reflect that variation.
In this study, the Generate Spatial Weights Matrix tool module in the ArcGIS version 10.3 (ESRI, 2016) was used to create a binary file defining the relationships among features in the dataset.These relationships are utilized in the mathematics the spatial statistics tools.

Logistics Regression
The analysis of neighborhood characteristics alone does not indicate what extent the spatial pattern land use changes.An analysis of the explanatory capacity of each factor is helpful for determining the relevance of neighborhood interactions in a particular case study (Verburg et al., 2004b).For this study, we used the neighborhood interaction in a logistic regression model to relate the change event with the explanatory capacity of many factors (independent variables).As statistical modeling, the logistic regression model estimates the probabilities of the occurrence of change events as a dichotomous dependent variable.
In logistic regression, the probability of conversion of a grid cell (P) is described as a function of a set of predictor variables following: where the independent variables (Xk) are a set of k predictor variables and βk are the coefficients to be estimated with a maximum likelihood estimation for k predictor variables.The number of independent variables included in the equation can be very large considering many interactions with different specification can be included in the model.The selection of interactions included in the model depends on the theoretical considerations of the researcher and the complexity of research.
In this study, some factors were selected from the biophysical environments to include in the analysis, such as road density, elevation and slope, and the value of spatial dependency that expresses the neighborhood characteristics.The resulting probability was compared with the locations that actually changed to determine the goodness of fit of the regression model.The goodness of fit is the measure of the amount of spatial variation of land use change that can be explained by neighborhood characteristics.
The logistic regression model was applied to Java Island to predict the future process of deforestation (Prasetyo et al., 2009).The study indicated that the highest deforestation occurred in the eastern part, followed by western of Java.Moreover, the result of statistical tests on the model shows that population density, road density and household income from the agricultural sector were significant parameters to predict the deforestation processes on Java Island.

Land Use Change Trends by Markov Model
The Markov model has been used to model changes in land use and land cover on a variety of spatial scales (Weng, 2002).Markov analysis in a small area of less than a few hectares tends to focus on changes in land cover/vegetation type (Baker, 1989).Meanwhile, land use studies using Markov chain models tend to focus on a much larger spatial scale, and involve both urban and non-urban covers (Jahan, 1986;Muller and Middleton, 1994).
In the case of land use change simulation, a basic assumption of the Markov model is to regard land use change as a stochastic process that move in a sequence of steps through a set of states (Stewart, 1994).This process gives a value at time t, Xt, and depends only on its value at time t-1, Xt−1, and not on the sequence of values Xt−2, Xt−3.... X0 that the process passed through in arriving at Xt−1.It can be expressed as: Moreover, it is convenient to regard the change process as one which is discrete in time (t = 0, 1, 2....).
The P {Xt = aj || Xt−1 = ai}, known as the one-step transitional probability, gives the probability that the process makes the transition from state ai to state aj in one time period.When steps are needed to implement this transition, the P {Xt = aj || Xt−1 = ai}is then called the step transition probability, Pij.If the P′ij is independent of times and dependent only upon states a i , a j , and then the Markov chain is said to be homogeneous.The treatment of Markov chains in this study will be limited to first order homogeneous Markov chains.In this event: where Pij can be estimated from observed data by tabulating the number of times the observed data went from state i to j, nij, and by summing the number of times that state ai occured, ni. Then; As the Markov chain advances in time, the probability of being in state j after a sufficiently large number of steps becomes independent of the initial state of the chain.When this situation occurs, the chain is said to have reached a steady state.Then the limit probability, Pj, is used to determine the value of P′ij: Where The transition area predicted by the trend of change on each period will then be assessed using the actual condition in 2007.The change model with more accurate prediction results relatively, will be used to predict the near future role of land use change pattern.

Neighborhood Characteristics of Land Use Changes
Neighborhood characteristics explored in this study focused on the spatial autocorrelation between the change areas and their neighborhood.A neighborhood weight value represents the effects of each land use type to the occurrence of change.Consequently, if the value for a land use type is higher, it means that the change area is positively correlated with the existence of the land use type.The characteristics of the neighborhood for each land use type from 2001 to 2007 are shown in Figure 3. Based on the results, the most significant neighborhood effects are found between a spatial weight of 0.1 and 0.25.The neighborhood effects of some land use types can be identified: timber forest, rubber and oil palm plantation areas show a slightly higher relationship with the change areas and indicate that the change area in many sites has been affected significantly by their existence.As the change mechanisms explained in Setiawan et al. (2014) indicate most of these changes rely on land management factors over a large area.Moreover, some specific agricultural lands such as paddy rice with irrigation system (double cropping), especially in upland areas, have a high positive spatial autocorrelation with the change areas.The positive effect of this agricultural land suggests a land is allocated to other uses, including paddy rice field.Moreover, Figure 3 indicates that the change area has also a significant spatial correlation with the change area itself, which means that the change in an area influences change in other areas over some distance.However, the result indicates the neighborhood characteristics of settlements (residential area) are not quite taking place at an equal pace in the change areas since the smaller change areas surrounding settlements cannot be defined by the method.

Logistic Regression Based on Neighborhood Characteristics
As mentioned in Setiawan et al. (2014), the main types of land use change in Java are recognized, as follows: (1) agricultural development, including some trajectories such as non-agricultural lands converted into intensive agricultural lands, (2) change in forestlands, either by natural factors, such as forest fire and volcanic activities, or human factors which is converted into agricultural land (upland/plantation), (3) change in plantation area, for example rubber plantation changed to oil palm plantation, and (4) urban development, e.g. a conversion from agricultural land into settlement/other facilities.Each change category has specific change mechanisms or processes, which then reveals the specific spatial model for each type.However, the datasets used allowed me to consider the first two change categories in developing a spatial model for the conversion probabilities using logistic regression.The coefficients of the regression equation for the conversion probability are given in Table 3.
Table 3 Coefficients of a logistic regression model explaining the spatial distribution of change area in agriculture and forest land by neighborhood characteristics, topography and road network (n= 6.557) The coefficients result above show that the change in agricultural land is related positively with paddy rice with double-irrigated cropping system in upland areas.Several lands use types, namely: settlement, upland with intensive cropping, mixed garden and paddy rice with triple-irrigated cropping system indicate a high possibility of change, but their contributions are not significant.Meanwhile, topographic factors, elevation and slope, have a negative impact on the occurrence of change area in agricultural land.In other words, the lands with agricultural uses in lowland or flat areas convert easily to other uses.Moreover, barren lands/dry land, bush-shrub and mixed garden give a negative impact on the change occurrences in agricultural lands.Based on the land survey results, most of these land use types are located in dry areas with little water.Moreover, distance to the road makes a significant contribution to the model and the negative coefficient result indicates that the changes are greater near the road than at some distance away from the road.
In the case of forestland, the results indicate that some land use types such as upland with intensive cropping, mixed garden with bush and oil palm have positive contribution to the change of land, however, most of them are not significant by statistics test (Table 3).The negative contribution to the occurrence of change in forestland is significantly provided by paddy rice field, mixed garden, forest plantation, and rubber plantation.
In this study, we included all significant variables of the biophysical-environment parameters in the logistic regression model.Based on this model, the conversion probability was calculated for the agricultural land as well as forestland that could potentially be changed.Locations with high probabilities for agricultural land use and forest are shown in Figure 4 and Figure 5, respectively.

Land Use Change Prediction on Regional Scale Using Markov Model
The land use transition predicted in Markov matrix from each period was assessed by predicting the actual land use in 2007.In order to simplify the complexity, the trends among 25 specific land use types, they were categorized into: 1) upland land use types (agricultural land in dryland areas), 2) paddy rice field, 3) forested lands, 4) plantations including timber forest plantation and 5) settlement/residential area.In Figure 6, the expected area for each category based on Markov model on each period is compared with the actual land use.The results indicated that the area of upland (agricultural land in dry land areas) can be projected regarding their trends in period 2001-2002, 2002-2003 and 2003-2004, with  Similar with the trends in agricultural upland, the paddy rice field in Java is also projected decreasing in the future (Figure 6).The Markov model's result shows that the pattern of change during 2001-2002 gives the best prediction for the transition areas in 2007, which is only 0.22% higher than the actual area in 2007.Instead of period 2005-2006, the predicted areas based on period 2002-2003, 2003-2004 and 2004-2005

CONCLUSION
This paper introduces a biophysical spatial modeling for land use change on Java Island considering neighborhood interactions between land use types and the change area.Then, these neighborhood characteristics used in logistic regression model to estimate the probability of the change events occurrence.Moreover, the future role of land use change is then projected using the Markov model based on the annual land use changes map from 2001 to 2007.
The results indicate that paddy rice with irrigation system (double cropping), especially in upland areas has a high positive spatial autocorrelation with the change areas.Residential area, paddy rice, and upland with intensive cropping have a high effect to the probability of change occurrence.Meanwhile, barren lands/dry land, bush-shrub and mixed garden give a negative impact to the change occurrences in agricultural lands.In the case of forestland, the results show some land use types such as upland with intensive cropping and plantation have positive contribution to the change of land.Even though the trend of land use pattern was projected, but the results of logistic regression model indicate that the interaction on the change area differ in different parts over Java Island, due to different environmental and socio-cultural conditions.
The future role of land use change was projected using the transition area and probability of Markov model.The accuracy of the model was also assessed through comparison the projection with the actual area in 2007.The results indicate that the future role of each land use type is different based on the trend period predictor in the model.For example, agricultural land in dry land areas/upland can be projected regarding their trends in period 2001-2002, 2002-2003 and 2003-2004, with different areas are 0.65%, 0.38% and 0.35%, respectively.

Figure 1 Figure 2
Figure 1 Distribution of land use/land cover change as detected by temporal vegetation dynamics change in Java (Source: Setiawan and Yoshino, 2014)

Figure 3
Figure 3 Spatial weight matrix of the change area corresponding to each land use type

Figure 4
Figure 4 Probability map of agricultural land use change based on logistics regression model

Figure 6
Figure 6 Transition area of land use expected based on trend of each period compared with the actual area in 2007 different areas are 0.65%, 0.38% and respectively.Meanwhile, the considerably increasing on this land use type from 2004 to 2005 caused the model overestimates the expected area in 2007.In the other side, the result of the model underestimates the area when the projection is based on period 2005-2006.Such condition indicates the dynamics change on the lands that occurred significantly from 2004 to 2006.Regarding on the Figure 6, the agricultural lands (uplands) in Java Island is in the direction of decreasing with time.Most of the change area of this land use type is expected distributed in the areas with high conversion probability on Figure 4.
are lower than the actual area in 2007.In the case of forested lands, the area in 2007 can be predicted regarding their trends either in period2002-2003 (0.29% lower) or 2003-2004 (0.49% higher).Meanwhile, these predicted areas are significantly lower then the actual area when predictions regarding the trend area in 2001-2002 (2.95% lower), 2004-2005 (2.08% lower) and 2005-2006 (4.59% lower).Figure6shows the forest timber plantation, rubber and oil plantation have similar trends with the forested area.Their future role can be predicted based on their trends area in period2002-2003 and 2003-2004, which are 0.07% lower and 0.13% lower, respectively.Different with others categories, settlement (residential area) can be projected based on most of trends in periods, since the prediction result from all periods are less than 1% different with the actual area in 2007.Even the change pattern during 2002-2003 gives the best prediction for the transition areas in 2007 (0.06% higher).

Table 1
The list of land use types as a summary of the analysis of EVI pattern (Setiawan et al., 2013) Number of LU Type Land Use Type LU-1 Mixed garden LU-2