Predicting Land Cover Change in the Mamminasata Area, Indonesia, to Evaluate the Spatial Plan

Abstract

The spatial plan program for Makassar City and the surrounding area called Mamminasata (Makassar, Maros, Sungguminasa, and Takalar) was created by the Indonesian Government. The program regulates the proportion of land cover, but predictions about land cover changes were not considered. Therefore, in this study, we predict what the land cover may be in 2031 using the multi-layer perceptron neural network and the Markov chain methods. For this purpose, image composite, support vector machine classifier, and change detection were applied to a time series of satellite data. Visual validation showed the hot-spots of land cover changes related to population density, and statistical validation scored 0.99 and 0.78 in no information kappa and grid-cell level location kappa, respectively. The model was performed to predict land cover in 2031, and the predicted result was then compared with the spatial plan using an overlapping method. The results showed that built-up area, dryland agriculture, and wetland agriculture occupied two, twenty, and eight percent of the protected zone, respectively. Meanwhile, fifteen percent of the development zone was covered by forest, mainly in the eastern part of Mamminasata. The result can be used to help the Government decide future plans for the Mamminasata area.

1. Introduction

The concept of sustainable development was established at Agenda 21 in Rio de Janeiro [1], leading to many countries, especially developing countries, striving to achieve the goals of the concept. The core concept comprised of three pillars, namely, social, economic, and environmental [2]. Many issues that affect sustainable development are geospatially interrelated and can be analyzed, modeled, and mapped within a geographic context [3]. Indonesia, as a developing country and a newly industrialized country that flourished in the mid-1990s [4], has applied the concept of sustainable development in spatial plan programs. These programs contain regulations related to the proportion of land cover area in a region. In the Mamminasata area, South Sulawesi province, the local government has established the spatial plan since 2006. Several issues related to the condition of land cover and spatial plan in this area had been addressed [5], such as: (a) the existing spatial plan was paying less attention to its implementation, especially for the conservation, protection and hazardous areas; (b) balance between the current land cover and future development were not fully taken into account under the existing spatial plan; and (c) lack of coordination among government agencies and institutions. Another important issue was that the land cover prediction has not been considered in the spatial plan.

Several researchers have studied the relationship between sustainable development and land cover change. García-Ruiz et al. [6] linked the effects of land use change on sustainable development in the Pyrenees mountain areas, Spain. Musa et al. [7] modeled the development of cities based on land cover changes for sustainable development in the Niger Delta Region, Nigeria. Salazar et al. [8] combined urban planning, natural conservation, and risk areas to build a prediction model for land use change as a guide for sustainable development in the Quito Metropolitan District, Ecuador. Land cover change is an important part of achieving sustainable development.

To simulate and predict land cover changes, a variety of models have been applied, mainly by combining remote sensing (RS) and a geographic information system (GIS) [9,10]. Losiri et al. [11] classified those models into four types: (a) empirical and statistical models such as the Markov chain and logistic regression; (b) dynamic models such as cellular automata, agent-based model, genetic algorithm, artificial neural network, and system dynamic; (c) integrated models such as conversion of land use and its effects at small regional extent (CLUE-S) and Dyna-CLUE; and (d) hybrid models such as Metronamica, land transformation model, land change modeler (LCM), and slope, land use, exclusion, urban extent, transportation, and hill-shade (SLEUTH).

The multi-layer perceptron neural network and the Markov chain (MLPNN + MC) in LCM was found to produce a considerably higher prediction accuracy than other models [11,12,13,14,15]. This higher prediction accuracy stems from neural networks being able to express changes in various land cover types more adequately than single probabilities such as the weights of evidence method [13]. The three layers (input, hidden, output) in MLPNN are used to express complex non-linear relationships between land cover changes and driving factors that produce transition potential, while the transition potential matrix is calculated by using MC to predict future land cover changes [10]. Based on the work of Pickard et al. [16], LCM has four specific characteristics and advantages: (a) flexible input data in pre-processing and the data collection process; (b) probability surface, exclusion layer, and regional stratification in the calibration process; (c) stochastic modeling to generate the model projections process; and (d) dynamic variables allowed and scenario analysis in extrapolating future scenarios. This model has been used in several studies on urban growth simulation [14,17,18], spatial dynamics of deforestation and fragmentation [19], land use changes following green infrastructure policies [20], and so on.

To obtain land cover classification, the support vector machine (SVM) with its non-linearity and multidimensional capabilities is a good choice because this method has a powerful image classification technique for handling multispectral satellite imagery [21]. For two or more periods, change detection can be applied to monitor land cover dynamics. Change detection methods are grouped into two categories: unsupervised direct comparison (UDC) and supervised post-classification comparison (SPCC) [22]. UDC is relatively simple and straightforward, and does not require a manually labeled training set [23], whereas SPCC uses an image object as the basic unit in change analysis and simultaneously detects the area and type of change [24]. Another simple change detection method is change vector analysis (CVA), which computes spectral change vectors using given multi-dates and compares their magnitudes with a specified threshold criterion [25]. Land cover classification is an important part of predicting land cover change, and is a common research tool that is often used for determining actual [26], previous [27,28,29,30], and future [9,31,32] land cover. Land cover classification and change detection are also supported by availability and access, as well as its continuity and affordability of medium- to low-resolution remote sensing data that can be provided free of charge [33].

The purpose of this research was to predict land cover in 2031 in the Mamminasata area, Indonesia, by using the multi-layer perceptron neural network and the Markov chain. The support vector machine and change vector analysis were used to obtain land cover classification in 2006, 2011, and 2016. Seven vector maps made by the Government were used as driving factors. The prediction model was validated using visual and statistical methods. The predicted land cover in 2031 was compared and evaluated with the Government’s spatial plan for the period from 2011 to 2031 using the overlapping method.

2. Research Site and Materials

2.1. Research Site

The Mamminasata Urban Area in Sulawesi Island, Indonesia, has an area of 247,000 ha and consists of Makassar, Maros, Sungguminasa, and Takalar municipalities (Figure 1). The area was established by a decree of the Governor of South Sulawesi province in 2003 because of the development of Makassar City, the capital of the South Sulawesi province. Makassar City has grown rapidly, resulting in agglomeration with the other three municipalities. In addition, the Mamminasata area is the economic center area in Eastern Indonesia and contributes a significant proportion of the gross regional domestic product of the South Sulawesi province: 39% in 2006, 44% in 2011, and 45% in 2016 [34,35,36].

The Government has created a spatial plan for 2011–2031 to make the Mamminasata Urban Area the center of international-scale services in Eastern Indonesia with the aim of achieving sustainable development. To bolster the plan, the Indonesian Government established the Mamminasata Urban Area as a national strategic area, with the aims of having this area act as a model metropolitan area and as an exemplar of urban development for Indonesia [37].

2.2. Materials

2.2.1. Satellite Data

Landsat 7 Enhanced Thematic Mapper plus (ETM+) and Landsat 8 Operational Land Imager (OLI) were used in this study because of the availability of long-term data free of charge. The observation dates are listed in

Table 1. Observation dates of Landsat in three periods.

Year	Satellite and Sensor	Observation Dates
2006	Landsat 7 ETM+	15/05/2005	21/07/2006
		11/02/2006	07/09/2006
2011	Landsat 7 ETM+	10/03/2010	21/09/2011
		11/04/2010	31/03/2012
		27/04/2010	09/10/2012
		29/05/2010	18/03/2013
		02/09/2010	22/06/2013
		13/03/2011	08/07/2013
		16/05/2011	25/08/2013
		03/07/2011	26/09/2013
		05/09/2011
2016	Landsat 8 OLI	15/02/2016	24/07/2016
		05/05/2016	10/09/2016

. Two consecutives scenes, 116/63 and 116/64 for path/row, were downloaded from the United States Geological Survey (USGS) [38]. Six bands from 1 to 7 except band 6 were used from Landsat 7 ETM+, and corresponding bands from 2 to 7 were used from Landsat 8 OLI. Map projection was in the Universal Transverse Mercator (UTM) Zone 50 South, and spatial resolution was 30 meters. Data were converted into top of the atmosphere reflectance in 8 bits. We selected less cloudy scenes from two seasons in Indonesia (wet season from November to May and dry season from June to September). In addition, multiple scenes were used for each season to avoid the remaining cloud and its shadows, and to fill the scan gaps in Landsat 7 ETM+.

2.2.2. Other Spatial Data

Vector data of capital, river and lake, road and settlement were used as driving factors; land use as the training data; and administration area as a boundary map of Mamminasata. In addition, elevation and slope were used based on a 30-meter resolution digital elevation model (DEM) generated by the Shuttle Radar Topography Mission (SRTM). The population density of the Mamminasata area in 2011 was used to calculate a population density per pixel as a driving factor. All data were converted into raster format at a resolution of 30 meters to apply LCM. These data were obtained from the database of the Centre for Regional Development and Spatial Information (WITaRIS), Hasanuddin University [37].

3. Methods

The integrated method used in this research required three mains steps. In Step 1, a land cover map for 2011 was derived using composite method and SVM classifier, and land cover maps for 2006 and 2016 were derived by change detection. In Step 2, land cover change was modeled by MLPNN + MC using three land cover maps, several driving factors, and a parameter set. MLPNN was carried out to obtain transition potential between 2006 and 2011, while MC was applied to predict land cover in 2016 and 2031. In Step 3, predicted land cover in 2031 was overlaid with the Government’s spatial plan. The result was used to evaluate the spatial plan and to help guide the future development of Mamminasata.

3.1. Land Cover Mapping

We used two schemes to obtain the land cover classification maps (Figure 2). First, we generated a land cover map for 2011 through image classification (left side). Second, land cover maps for 2006 and 2016 were generated through change detection with 2011 (right side). The right side of Figure 2 shows the case for 2011 and 2006, but the flow is the same for 2011 and 2016.

3.1.1. Land Cover in 2011

Areas covered by cloud and shadow are obstacles in optical remote sensing. The cloud masking method is an important solution to this problem [39]. We used cloud mask included in Landsat data to detect cloud and its shadow. In addition, image processing is often necessary prior to land cover classification, especially with the presence of striping scan gaps as occurred in the Landsat-7 ETM+ images [40] because of the failure of the Scan Line Corrector (SLC) since May 2003, which causes systematic data gaps on the captured imagery and eliminated the capacity to provide spatially continuous fields [41].

To fill-in the cloud and shadow areas and scan gaps in the Landsat ETM+ data, a two-step composite was applied prior to the land cover classification as follows. In the first step, two cloudless scenes were selected from images observed in the same season to minimize seasonal change (

Table 2. List of composite scenes.

	Scene 1	Scene 2	Scene 3
Step 1	11/04/2010 27/04/2010	03/07/2011 22/06/2013	05/09/2011 21/09/2011
Step 2	16/05/2011 13/03/2011 31/03/2012 10/03/2010 18/03/2013 29/05/2010	08/07/2013 02/09/2010 25/08/2013	09/10/2012 26/09/2013 02/09/2010 25/08/2013

). A pixel-based composite was created using these images and based on the following rules: (1) no data were selected if both pixels were cloud, shadow, or scan gap; (2) if one of the two images was clear data, it was selected; and (3) if both data were fine, the difference between the near infrared band (Band 4) and red band (Band 3) was calculated, and the pixel with the larger difference was selected. The last rule was adopted to avoid a thin cloud and its shadow, and the smoothness of the image could be kept by selecting neighboring pixels from the same image.

Although most of the pixels were filled-in during the first step, blank pixels remained if both image pixels were cloud, shadow, or scan gap. Blank pixels are unfavorable for land cover classification, so the second step was applied only to those blank pixels. Candidate images were selected from a similar season, although these images were from longer range of observation dates compared with the first step because fewer cloud-free scenes were available for our site. If the first composite image was a blank pixel, it was filled-in using the following procedure: (1) the 5 × 5 pixel block with the blank pixel at the center was selected from the composite image; (2) the sum of absolute difference between the composite image and each candidate image for all bands was calculated by using clear data; and (3) the pixel with the minimum difference was selected from the image.

The two-step procedure was adopted because a blank pixel has small patches, and a smooth image was generated by using a pixel block for the comparison. We started with a small number of candidate images and increased the number of images by checking the composite result. Using this composite method, three scenes were generated for different seasons to achieve high classification accuracy. We could not generate better images for other seasons because of heavy cloud cover and fewer available images.

A support vector machine (SVM) was adopted as a land cover classifier because it has been evaluated as a high-performance machine-learning algorithm and has been investigated in a number of studies [26,30,42,43]. In this study, the pixel-based classification was applied using six bands of three composite scenes. The pixel-based technique is generally used for land cover classifications with medium-resolution imagery at the regional level [44]. In addition, SVM is better at recognizing subtle patterns in complex datasets compared with other machine-learning methods [45]. This classifier is a kernel-based supervised-learning algorithm, which combines machine-learning theory, optimization algorithms from operations research, and kernel techniques from mathematical analysis [46].

The training area was extracted from the land use map that was produced by the Government in 2011 [37]. The original map contains 29 land use types, which we merged into six land cover types (built-up area, dryland agriculture, forest, shrub, waterbody, and wetland agriculture) for our purpose to apply the land cover map for spatial prediction modeling. The 30-meter-wide boundary of land use polygons in the training map was excluded to eliminate mixture of land cover in the training samples. Spectral reflectances of Landsat data (6 bands x 3 scenes) were extracted with respect to each land cover by overlapping the training map and Landsat images. These processes were carried out using ArcGIS software (ESRI, Redlands, CA, USA). The reflectances of all land cover types were subsequently input to the Statistic and Machine Learning Toolbox of MATLAB software (MathWorks, Natick, MA, USA). All the available samples were used as training data for the SVM classifier.

Because almost the whole area except the boundary was covered by the training data, the classification map was similar to the Government map. Therefore, accuracy was assessed only in areas that had different land cover types between the training data and the SVM classification. We picked 100 check points for each land cover class from a Google Earth image. A confusion matrix was generated and producer’s accuracy (PA), user’s accuracy (UA), overall accuracy (OA), agreement of chance (AoC), and kappa coefficient (KC) were calculated using the following formulas [47]:

$$PA=\frac{Ca}{Cg·sum}$$

(1)

$$UA=\frac{Ca}{Ci·sum}$$

(2)

$$OA=\frac{Cc·sum}{C·sum}$$

(3)

$$AoC=\frac{PM}{{\left(C·sum\right)}^2}$$

(4)

$$KC=\frac{OA-AoC}{1-AoC}$$

(5)

where

Ca	:	number of points where classification is correct,
Cg·sum	:	total ground truth points per class,
Ci·sum	:	total classification points per class,
Cc·sum	:	total number of points where classification is correct,
C·sum	:	total number of points,
PM	:	sum of the multiplication of ground truth points and classification points with respect to each land cover.

3.1.2. Land Cover in 2006 and 2016

Land cover classification maps in 2006 and 2016 were obtained through change detection with 2011, and land cover was identified only in the pixels where change was detected. Asokan and Anitha [48] classified the change detection approach into five types: algebra based, transform based, classification based, GIS, and an advanced model. Hussain et al. [44] summarized the differences between several change detection techniques into three types: pixel based, object based, and spatial data mining. Based on these studies, we used a pixel-based technique through an algebra-based approach (image differencing). This technique is based on the direct comparison between two images; the advantages are simple and easy to interpret [44]. In addition, this technique is suitable for use in medium-resolution imagery such as that from Landsat [44]. A common procedure in algebra-based change detection is the threshold selection to find the changed area [48].

In this change detection, the same method was applied to two periods (2011–2006 and 2011–2016) individually using four pairs of images listed in

Table 3. List of image pairs.

2006	2016
12/02/2012–11/02/2006 16/05/2011–15/05/2005 03/07/2011–21/07/2006 05/09/2011–07/09/2006	12/02/2012–15/02/2016 16/05/2011–05/05/2016 03/07/2011–24/07/2016 05/09/2011–10/09/2016

. The degree of change was calculated using the following equation:

$$D_n\left(i,j\right)=\sqrt{\frac{1}{P\times 6}{\displaystyle \sum }_{p=1}^P{\displaystyle \sum }_{b=1}^6{\left(\frac{D_{p,b}\left(i,j\right)-\overline{D_{p,b}}}{{\sigma }_{D_{p,b}}}\right)}^2}$$

(6)

where

$D_n\left(i,j\right)$	:	normalized difference for the image pixel at location $i$ (line) and $j$ (column),
$P$	:	number of valid image pair (maximum is 4 but it is decreased in case of cloud, shadow, and scan gap),
$D_{p,b}\left(i,j\right)$	:	difference in digital number between 2011 and 2006 (or 2016) at pixel $\left(i,j\right)$ ($b$ and $p$ are index for band and image pair, respectively),
$\overline{D_{p,b}}$	:	average of $D_{p,b}\left(i,j\right)$ for all pixels in a scene,
${\sigma }_{D_{p,b}}$	:	standard deviation of $D_{p,b}\left(i,j\right)$ for all pixels in a scene.

For Band (b), we used six bands from 1 to 7 except band 6 from Landsat 7, and Bands from 2 to 7 from Landsat 8. This equation means that the same weight was given for all bands because the variation in difference by band was normalized using average and standard deviation of the pixel difference by each band. It should be noted that the difference could be calculated if both image pixels have valid data. If one of the pair is cloud, shadow, or scan gap, the pixel is excluded from the calculation of D_n(i,j). In order to minimize the difference of TM and ETM+ spectral responses, the normalized difference was calculated using top of the atmosphere reflectance, not the digital number of the image. In addition, the residual difference could be eliminated by subtracting the average of difference ($\overline{D_{p,b}}$) in normalization process. After calculating this normalized difference, the land cover change was detected by applying the threshold, which was selected manually by comparing the detected area with a time series of Google Earth fine-resolution satellite images.

The new land cover was identified for the changed area based on the following method: (1) sample data (20,000 pixels) were selected from unchanged areas in 2006 (or 2016) images for all six land cover types; (2) for each changed pixel, the differences with the above sample data were calculated using Equation (7); (3) these differences were sorted in ascending order; and (4) by evaluating the first 20,000 samples, land cover was identified by the most common land cover among the six types.

$$D\left(c,s\right)=\sqrt{\frac{1}{P\times 6}{\displaystyle \sum }_{p=1}^P{\displaystyle \sum }_{b=1}^6{\left(\frac{R_{p,b}-R_{c,s,p,b}}{{\sigma }_{R_{c,p,b}}}\right)}^2}$$

(7)

where

$D\left(c,s\right)$	:	normalized difference for the image pixel at location $i$(line) and $j$(column),
$P$	:	number of valid image pair (maximum is 4 but it is decreased in case of cloud, shadow, and scan gap),
$R_{p,b}$	:	digital number of changed pixel ($b$ and $p$ are index for band and image pair, respectively),
$R_{c,s,p,b}$	:	digital number of sample data as same band and period with $R_{p,b}$ for sample data $s$ in land cover class $c$,
${\sigma }_{R_{c,p,b}}$	:	standard deviation of $R_{c,s,p,b}$ for all samples.

This method assigns the changed pixel to the nearest class based on the spectral reflectance of each land cover sample extracted from an unchanged area. Therefore, the method allows change between the same classes, e.g., from dryland agriculture to dryland agriculture, even if the pixel was regarded as a changed pixel. We adopted this “change detection and classification” scheme in order to cope with the variations of spectral reflectance due to sensor and season. These biases were compensated by comparing the reflectance of changed pixel with samples extracted from same scene. Change detection and subsequent determination of land cover type were carried out using self-coded software.

3.2. Land Cover Change Model

The land change modeler (LCM) module embedded in TerrSet software (Clark Labs, Worcester, MA, USA) was used in this study. The module developed by Clark Labs analyzes land cover changes, evaluates driving factors, calculates transition potential from one class to another, and builds a model prediction of land cover changes empirically, which can be tested for accuracy [10,18,49,50]. This module has been recognized as an effective tool in predicting land use change, with an easy-to-understand user interface [10].

To model land cover change, LCM requires three maps: maps from two previous periods are used to optimize the prediction model for the third period by comparing with the map of the third period [51]. We used the 2006 and 2011 maps to build the model by referring to the 2016 map. Figure 3 shows a land cover change modeling flowchart. Land cover maps in 2006 and 2011 were used to predict land cover in 2016 by applying the MLPNN method coupled with driving factors. Predicted land cover for 2016 was generated by applying the MC method, and the prediction result was validated visually and statistically. If the accuracy value is acceptable, the model can be applied to the prediction for 2031. Consequently, land cover prediction for 2031 was overlaid with the Government’s spatial plan map for 2011–2031. The result was used as an evaluation map for future development of the Mamminasata area.

3.2.1. Driving Factors

A driving factor is a variable that forces a change in land cover. Seven driving-factor maps (distance from capital, distance from river and lake, distance from road, distance from settlement, elevation, slope, and population density per pixel) were used to model land cover prediction. These driving factors were selected due to the availability of these data. The slope was calculated in the Slope module using DEM. The four distance maps were processed in the Distance module using vector maps of capital, river and lake, road and settlement, respectively. The population density per pixel was obtained using the following formula:

$$Kp=p\times A\times P\times C$$

(8)

where

Kp:

population density per pixel,

p:

population density of Mamminasata (= 981.29 people/km² [52]),

A:

population distribution area (= 3.14 $\times$ (2 km)² = 12.56 km²),

P:

population proportion expressed by Equation (8),

C:

conversion factor, from 1 km² to 1 pixel

 

(= 30 m $\times$ 30 m = 900 m² = 9 $\times$ 10⁻⁴ km² = 9 $\times$ 10⁻⁴ pixels).

The population proportion formula is as follows:

$$P=0.2402\times \exp \left(-0.9464\times \left(\frac{distance from settlement}{1000}\right)\right)$$

(9)

The map for population density per pixel was created assuming the population spreads radially with a 2 km radius from settlements, and the population increases as it approaches the center of its distribution [53].

3.2.2. Parameter Setting

Land cover maps for 2006 and 2011 were used as the basic maps to predict land cover. In this study, the cross-tabulation method in the change analysis stage was applied to 2006–2011 and 2011–2016. This method was used to identify and quantitatively assess the change in land cover from one class to another by producing a gains and losses graph.

In parameter setting, the basic maps (2006 and 2011) and the seven driving-factor maps were first fed into the model to produce a transition potential map, which indicates the probability of change from one class to another by using the MLPNN method [10]. We determined the land cover change by choosing a transition sub-model, which can consist of a single land cover transition or a group of transitions that are thought to have the same underlying driving factors [50]. In this study, we chose seven sub-models with a high transition potential and actual condition of our research site: from “dryland agriculture to built-up area”, “shrub to built-up area”, “wetland agriculture to built-up area”, “shrub to dryland agriculture”, “shrub to forest”, “waterbody to wetland agriculture”, and “wetland agriculture to dryland agriculture”. The selection of these sub-models could regulate the land cover changes unnecessary for the evaluation purpose, that is, we could remove the natural changes such as seasonal change in water surface.

As the second step in parameter setting, each driving-factor map was determined as a static or dynamic type. Both types are based on their behavior of change over time [10]. A static variable is a constraint or inhibiting factor for the transition, and does not change over time, while a dynamic variable, such as a development program and infrastructure [10,50], is the driver for transition and is updated during the model run. In this study, distance from the road was considered a dynamic variable and the other six driving factors as static variables.

Cramer’s V was used to test the relationship between the seven driving factors and land cover changes in 2006–2011. Cramer’s V is a common chi-square-based measure of nominal association, the most suitable measure of association, and a complete agreement formula between two nominal variables [11,19]. Values range from 0 to 1 and a higher score is better [10,49,54].

As the final step, we ran the model using the above-mentioned parameter setting. For other parameters, we used default values. In this process, the model generates transition potential maps based on the procedure that half of the pixel samples are used to run the model and the other half are used to validate the performance in predicting land cover changes [10,11]. This method is sophisticated, especially for complex problems and non-linear relationships [17].

The MC method was applied to predict land cover in 2016. The potential land cover transition is used by MC to describe actual changes and to predict the change in the following year [55]. The land cover prediction map for 2016 was validated with the actual land cover map. We used two validation methods: visual and statistical. Visual validation shows the prediction result using four categories: hit, null success, false alarm, and miss. Hit and null success represent the correctness of the model; false alarm and miss represent errors that result from the model as disagreement between predicted map and actual map [49].

To strengthen the accuracy test, the prediction result was validated statistically. Kappa index statistics were calculated for two types of target area: whole research site and different areas. For different areas, only areas that experienced changes between actual and predicted land cover maps for 2016 were tested. The accuracy result was evaluated using four kappa indices: kappa for no information that indicates the overall accuracy of the simulation run [56], kappa for grid-cell level location that indicates how well the grid cells are located on the landscape, kappa for stratum-level location that indicates how well the grid cells are located within the strata, and kappa for standard that indicates the proportion assigned correctly versus the proportion that is correct by chance [57]. These kappa values range from 0 to 1; the greater the value obtained, the more successful the prediction of land cover change. After the accuracy of land cover prediction in 2016 was accepted, the model can be used to predict land cover changes in 2031 with the same parameter settings.

3.3. Comparison with the Spatial Plan

The Indonesian Government created a spatial plan that consists of a spatial pattern map for 2011–2031. The spatial pattern map was constructed based on the actual land cover (when the map was made), the number and needs of the people (actual and future), Government policies, and analysis of the carrying capacity of the region. A spatial pattern is the distribution of land allocation in an area that includes a development zone and a protected zone [58]. A development zone is designated to be developed based on the conditions and potential of natural, human, and artificial resources, while a protected zone is designated to protect environmental sustainability, which includes natural and man-made resources [59].

The land cover prediction approach was not included in the spatial plan by the Government. Therefore, the land cover prediction for 2031 was overlaid with the Government’s spatial pattern map to provide a detailed land cover distribution in both the development and protected zones. This result could be helpful to the Government in considering future development in the Mamminasata area.

4. Results

4.1. Land Cover Maps of 2006, 2011, and 2016

Figure 4 shows the three composite images used for the classification in 2011, together with single-date images before the composite. Cloud, shadow, and scan gaps before the composite were reduced, and smooth and clear scenes were produced by the composite.

Most of the pixels were composed of two base scenes that have the smallest seasonal change in step 1, and a small fraction of the pixels were filled by other scenes in step 2 (based on Table 2). Therefore, the appearances were smoother than the simultaneous composite of all the scenes. However, there were still white areas (no-data, marked by red circles) because there was no valid data in any of the scenes due to the clouds and scan gaps.

The center panel in Figure 5 is the land cover in 2011 generated by SVM classifier, showing that the obstacles (cloud, shadow, and scan-gap areas) can be overcome because the differences in the reflectance were small and dissolved through statistical grouping by the SVM classifier.

Table 4. Confusion matrix.

		Classification Class						Total Row	PA
		BU	DA	Fo	Sh	Wb	WA
Ground Truth	BU	82	2	0	1	1	0	86	0.95
	DA	7	67	0	0	1	0	75	0.89
	Fo	11	22	98	25	3	4	163	0.60
	Sh	0	7	2	63	0	2	74	0.85
	Wb	0	2	0	5	93	2	102	0.91
	WA	0	0	0	6	2	92	100	0.92
Total Column		100	100	100	100	100	100	600
UA		0.82	0.67	0.98	0.63	0.93	0.92
OA		0.83
AoC		0.17
KC		0.79

Note: BU—built-up area; DA—dryland agriculture; Fo—forest; Sh—shrub; Wb—waterbody; WA—wetland agriculture; PA—producer’s accuracy; UA—user’s accuracy; OA—overall accuracy; AoC—agreement of chance; and KC—kappa coefficient.

shows the overall accuracy of land cover classification in 2011 was 0.83 with 600 points of total sample. The result indicated that the accuracy of land cover classification in 2011 was acceptable and can be used for land change modeling.

The land cover maps in 2006 and 2016 were obtained based on the change detection with 2011. The thresholds of change for normalized difference were 2.5 and 3.5 for respective years, determined by the visual comparison using a time series of Google Earth fine-resolution satellite images. As a result, total change areas were 4556 ha during 2006 to 2011, and 1301 ha during 2011 to 2016. The typical changed area of the final images can be seen in Figure 5. On the left side, some waterbody and wetland agriculture changed into built-up area; on the right side, some waterbody and shrub changed into dryland and wetland agriculture.

4.2. Land Cover Change Model

4.2.1. Land Cover Change Analysis

Figure 6 shows the gains and losses of land cover for two periods from 2006 to 2011 and from 2011 to 2016. In the analysis, there was no possibility of transition from built-up area to other classes, whereas other types of land cover had transitions between one another. Between 2006 and 2011, around 195 ha of wetland agriculture changed into built-up area, while dryland agriculture increased mainly from wetland agriculture (893 ha), shrub (746 ha), and waterbody (122 ha). In the same period, forest area increased, which was mostly obtained from shrub (347 ha). It can be seen that around 1472 ha shrub, 620 ha waterbody, and 1394 ha wetland agriculture changed into other types. However, there was an increase of around 843 ha of wetland agriculture from dryland agriculture, shrub, and waterbody. In the period between 2011 and 2016, the built-up area increased by 295 ha from wetland agriculture, 194 ha from dryland agriculture, and 130 ha from waterbody. This change is most likely because the Government program started during this period.

4.2.2. Land Cover Change Modeling

The seven driving factors are shown in Figure 7. Based on the relationship between driving factors and land cover in 2006–2011 using Cramer’s V test (

Table 5. Cramer’s V values of driving factors.

Variable	Overall Cramer’s V
Distance from Capital	0.3230
Distance from River and Lake	0.1103
Distance from Road	0.1965
Distance from Settlement	0.2542
Elevation	0.3263
Slope	0.2888
Population Density per Pixel	0.4920

), the highest and lowest influential factors were population density per pixel and distance from river and lake, respectively. It is reasonable that population density has a greater effect on land cover change, especially for agricultural and residential types of land cover. It is also natural that distance from water has less effect because many small rivers densely covered this area (as shown in Figure 7b), and enough water is available throughout the site.

The transition potential was obtained with 33.56% accuracy by running the MLPNN method from a combination of the basic maps (2006 and 2011) and the seven driving factors. This transition potential was used to simulate the land cover prediction for 2016 by using the MC method, and the resulting map is shown in Figure 8b.

After the land cover for 2016 was predicted, the resulting map was compared with the actual land cover map of 2016. Figure 8 shows the actual and predicted land cover for 2016 and the comparison of land cover classes between the two maps expressed as percentages. The two maps had almost the same values for all classes. The area of different land cover was 3125 ha.

Figure 9a shows the visual validation of the land cover prediction for 2016, displaying the four types of correctness and error: 0.00% hits, 98.89% null successes, 0.86% false alarms, and 0.25% misses. The correctness value (hits and null successes) indicated the same type of land cover between the actual map and the prediction map for 2016 (98.89%). Conversely, the error value (false alarms and misses) indicated the different types of land cover between the two maps (1.11%).

Statistical validation was carried out by comparing the whole area and different areas (areas of false alarms and misses).

Table 6. Accuracy test of change prediction.

Type of Kappa	Kappa Value
	Whole Area	Different Areas
No information	0.9961	0.9973
Grid-cell level location	0.9974	0.7821
Stratum-level location	0.9974	0.7821
Standard	0.9925	0.6544

shows the kappa values for the validation of the whole area and different areas. Pontius and Millones [60] revealed that the values of kappa for no information and grid-cell level location were somewhat more helpful than the values of kappa for stratum-level location and standard. Based on that, the results of accuracy can be accepted, and the model can be run for land cover prediction for 2031.

4.2.3. Land Cover Prediction for 2031

Figure 10a is the predicted land cover for 2031. The MC method was used with the same rules, seven driving factors, and basic maps (2006–2011) to predict land cover changes for 2031. Figure 10b shows the increases predicted for built-up area (252 ha), dryland agriculture (5578 ha), and forest (73 ha), and the decreases for shrub (2829 ha), waterbody (1282 ha), and wetland agriculture (1791 ha). Based on this prediction model, all types of land cover would have a transition potential between one another, including the built-up area, which did not experience a transition to other types of land cover in the period 2006–2016.

4.3. Comparing with the Spatial Plan

In the spatial pattern map from the Government, the total area of Mamminasata is divided into 80% development zone and 20% protected zone (Figure 11a). This map was overlaid with the land cover prediction map for 2031 to evaluate the Mamminasata spatial pattern (Figure 11b). The result shows that agricultural areas and waterbody occupy 70 and 7% of the development zone, respectively (Figure 12a), while 2% of the protected zone would become built-up area (Figure 12b).

The evaluation map also shows the small islands of Mamminasata called Spermonde, which are in the protected zone predicted to become built-up area (Figure 11c). Figure 11d shows waterbody in the development zone. The forest area occupies 15% of the development zone, mainly in the eastern and north-eastern parts of Mamminasata (Figure 11e,f). These areas are distant from the developed center, so special care is necessary to achieve sustainable development. Each land cover type in the development zone has a high possibility of changing into another land cover type in the future because the development zone is designed to be developed [59].

5. Discussion

In this section, we would like to address a few points based on our research aims. First, we have described a simple composite method to overcome the SLC problem of Landsat ETM+ coupled with the presence of cloud and shadow areas. This method proved to be useful in overcoming those obstacles in image processing. The land cover classification using the SVM classifier produced a precise land cover map with acceptable accuracy. The land cover classification map for 2011 using the SVM classifier in this study was helped by the availability of land use maps from the Government. We used this data as a training sample with the purpose of creating a land cover classification map without making new training data based on the visual interpretation of Landsat images.

As explained by Hussain et al. [44], the pixel-based technique using an algebra-based approach (image differencing) in this study has proven to be simple and easy to interpret. This technique is suitable for medium-resolution images such as those provided by Landsat. Based on the land cover maps of 2006 and 2016 shown in Figure 5 and the statistics in Figure 8, this technique helped us to detect an appropriate change in land cover. Our change detection method could be applied for another geographic scales and spatial resolutions if multispectral satellite data are available because the methodology is quite simple and general. Land cover is considerably complex with small patches in Asian countries, hence finer spatial resolution is preferable in order to reduce the mixture of land cover within a pixel. MultiSpectral Instrument (MSI) onboard the Sentinel-2 satellite would be a possible alternative because it has been operated continuously since 2015.

Although MLPNN showed a lower accuracy of transition potential than 80%, which is proposed by several researchers [10,11,50], our land cover prediction model for 2016 using the MC method showed the acceptable accuracy validated visually and statistically. Therefore, we would like to emphasize that MLPNN is a powerful method for modeling transition potential of land cover. It should be noted that our model only used current population and physical factors as driving factors. The economic factor was not included due to the unavailability of this data. Several researchers have proposed the economic factor in their land cover change model [11,61,62]. For further research, the economic factor should be selected into the model. All type of driving factor data (physical and socio-economic factors) could be included in the LCM. Furthermore, the influence order of driving factors will be determined using Cramer’s V.

Land cover prediction is not included in the current spatial plan in the study area. Therefore, the result of this study might be informative for current and future perspectives of the site. The importance of our result lies in the evaluation of the spatial pattern plan based on the land cover prediction. The evaluation map helped us to generate the detailed land cover distribution in both the development zone and the protected zone for 2031. Our model predicted that 30% of the protected zone would become built-up and agricultural areas, while forest accounted for 15% of the development zone. These conditions should be controlled by regulations in consideration of an unpredictable increase in population.

In reality, investment in land plots by both companies and individuals is also increasing rapidly with the establishment of Mamminasata as a national strategic area. Many agricultural, shrub, and waterbody (pond) areas have been converted into built-up areas, mainly residential and industrial areas, to meet the needs of the increasing population. Although the increase in built-up area cannot be avoided, adequate density and distribution should be considered together with environmental conservation.

Some regulations have already been enacted, for example, Indonesian Regulation No. 41 in 2009 (Protection for Sustainable Agriculture Land). Such regulations are needed because as an agrarian country, Indonesia needs to guarantee the provision of sustainable agricultural land as a source of work and decent livelihood for people by promoting the principles of togetherness, sustainability, environmental insight, etc., together with maintaining balance, progress, and national economic unity [63]. Indonesian Regulations (such as No. 41 in 1999 for Forestry, No. 17 in 2019 for Water Resources, and No. 26 in 2007 for Spatial Planning) are worthwhile and should be effectively put into practice. The evaluation map of the spatial pattern in this research could be used as guidance for the Government to apply these regulations to protect the conservation and hazardous areas.

In addition, local governments of the four municipalities in Mamminasata could cooperate to achieve common goals for development by balancing the development and the protected zones. Moreover, the National Government should strictly control investment in land plots to avoid the negative effects of the expansion of built-up areas. Managing the use of the Mamminasata area must be directed by zoning regulations, licensing, setting incentives and disincentives, and imposing sanctions for violations committed by any party. These actions should help to achieve sustainable development of the Mamminasata area.

6. Conclusions

Land cover of Mamminasata, Indonesia, in 2031 was predicted by the land cover change model using remote sensing and other geospatial data, and the resulting map was used to evaluate the spatial plan. The combination of predicted land cover and spatial plan helped to create a detailed map of future land cover in the development and protected zones. Significant land cover changes were expected in all four municipalities along with the increase in population density. The built-up area in the protected zone accounts for 2% of the total and agricultural areas would be 28% in 2031.

In terms of methodology, we found that the composite and change detection methods were simple and easy to follow. In addition, we showed that SVM is a powerful land cover classifier and useful even for medium-resolution images such as those provided by Landsat, and MLPNN + MC provides a good result in modeling land cover change prediction. It should be noted that the lower accuracy of transition potential in MLPNN could produce a good prediction of land cover, which was shown by the acceptable accuracy from the MC method (land cover prediction for 2016) validated by the actual land cover in 2016. This result is likely because the number of sub-models can affect the accuracy. Moreover, the transition potential accuracy in MLPNN should work well in most instances using default setting. As a suggestion, by experimenting the number of hidden layer nodes and the learning rate values in MLPNN as the two most critical parameters, it will affect the accuracy substantially.

Another topic that should be considered in further research is the prediction of population because land cover change is related to the needs of increased population, and the number of people can affect the carrying capacity of a region. These three topics are important in a spatial plan.

From our research, we suggest that the protected zone should be maintained or evenly expanded because Mamminasata is the economic center area in Eastern Indonesia coupled with its function as a national strategic area, which could result in an increase in population. Therefore, land cover prediction should be helpful in formulating a spatial plan, especially a spatial pattern plan to achieve goals for sustainable development.