Report submitted to EPA

July 2011

Spatial Modeling of Future Oil Palm 

Expansion in Indonesia, 2000 to 2022

Nancy Harris, Sean Grimland and Timothy Pearson
2 SPATIAL MODELING OF OIL PALM EXPANSION IN INDONESIA, 2000-2022 DRAFT REPORT

TABLE OF CONTENTS

SUMMARY.............................................................................................................................................................................5

1.0 BACKGROUND...............................................................................................................................................................6

2.0 METHODS........................................................................................................................................................................6

2.1 MODEL DESCRIPTION ......................................................................................................................................................7

2.2 DATA INPUTS ...................................................................................................................................................................8

2.3 CREATION OF FACTOR MAPS.........................................................................................................................................11

2.3.1 Elevation................................................................................................................................................................12

2.3.2 Precipitation ..........................................................................................................................................................14

2.3.3 Temperature...........................................................................................................................................................15

2.3.4 Slope ......................................................................................................................................................................16

2.3.5 Soil Type ................................................................................................................................................................17

2.3.6 Land Cover in 2001 ...............................................................................................................................................20

2.3.7 Infrastructure (Roads, Rivers, Railroads, Settlements) .........................................................................................21

2.3.8 Distance from Existing Plantations.......................................................................................................................24

2.3.9 Distance from Palm Oil Mills................................................................................................................................25

2.3.10 Peat Soils.............................................................................................................................................................26

2.3.11 Land Management Zones.....................................................................................................................................27

2.3 CREATING SUITABILITY-FOR-CHANGE MAPS ................................................................................................................28

2.4 EVALUATING SUITABILITY MAPS USING THE ROC STATISTIC.......................................................................................30

2.5 EVALUATING MODEL PERFORMANCE USING THE KAPPA STATISTIC..............................................................................31

3.0 RESULTS AND DISCUSSION .....................................................................................................................................32

3.1 EVALUATION OF SINGLE-FACTOR SUITABILITY MAPS FOR PREDICTING CHANGE.........................................................33

3.2 EVALUATION OF MULTI-FACTOR SUITABILITY MAPS FOR PREDICTING CHANGE..........................................................36

3.3 FINAL SUITABILITY MAP ...............................................................................................................................................37

3.4 PALM OIL EXPANSION, 2000-2022 ................................................................................................................................39

3.4.1 Expansion by State.................................................................................................................................................39

3.4.2 Expansion by Soil Type..........................................................................................................................................40

3.4.3 Expansion by Land Cover Type .............................................................................................................................41

4.0 CONCLUSIONS.............................................................................................................................................................42

5.0 REFERENCES ...............................................................................................................................................................43
FIGURES

Figure 1. Components of the land cover change model GEOMOD. From Pontius and Schneider (2001). 8

Figure 2. Representative dynamics of palm oil plantations in Indonesia between 2000 and 2009. 10

Figure 3. Candidate locations for palm oil expansion. Areas of expansion were assumed to occur only in the 

non-oil palm plantation areas. 11

Figure 4. Stratified factor map for elevation. 13

Figure 5. Stratified precipitation factor map for Indonesia. 15

Figure 6. Stratified temperature factor map for Indonesia. 16

Figure 7. Stratified slope factor map for Indonesia, with slope categorized into 5% bins. 17

Figure 8. Stratified soil factor map for Indonesia. 19

Figure 9. Stratified 2001 land cover factor map for the Indonesian islands of Sumatra and Kalimantan. 21

Figure 10. Stratified factor map for (A) roads; (B) rivers; (C) railroads; and (D) settlements for Indonesia. 23

Figure 11. Stratified distance from plantations factor map for Indonesia. 25

Figure 12. Stratified distance from palm mills. 26

Figure 13. Stratified peat soil factor map. 27

Figure 14. Management Zones reclassified 28

Figure 15. Example suitability map created from one or several factor maps, where blue is least suitable and red 

is most suitable. 29

Figure 16. Example ROC curve used to validate suitability maps for simulating land use change in Indonesia. 31

Figure 17. Ranking of ROC scores. 34

Figure 18. Ranking of single-factor suitability map ROC scores. 35

Figure 19. Locations of 2009 palm oil plantations located in areas of low suitability with respect to the elevation 

factor map. 36

Figure 20. Suitability map that incorporates all 14 factors equally weighted. 38

Figure 21. Predicted expansion of palm oil from 2000-2022 in the equal weight scenario. 38

Figure 22. Distribution of palm oil expansion among Indonesia provinces. 39
TABLES

Table 1. Total area (ha) of Sumatra and Kalimantan, Indonesia palm oil plantations in 2000 and 2009 according 

to IAI 2011. 9

Table 2. Factor data considered in the GEOMOD analysis. 12

Table 3. Elevation ranges of palm oil plantations (first ten bins) 13

Table 4. Rainfall suitability classes for growing oil palm. From Paramananthan (2003). 14

Table 5. Temperature suitability classes for growing oil palm. From Paramananthan (2003). 15

Table 6. Slope suitability classes for growing oil palm (first ten bins). 16

Table 7. Soil types in Indonesia and the reclassification scheme used to create the soil type factor map. 18

Table 8. Soil types for growing palm oil in Indonesia. 19

Table 9. Land cover area (ha) based on MODIS 2001 land cover for palm oilproduct. 20

Table 10. Distance from roads (first ten bins)for palm oil plantations. 21

Table 11. Distance from rivers (first ten bins) for palm oil plantations. 22

Table 12. Distance from railroads (first ten bins) for palm oil plantations. 22

Table 13. Distance from settlements (first ten bins) for palm oil plantations. 22

Table 14. Distance from 2000 oil palm plantations (first ten bins). 24

Table 15. Distance from palm oil mills (first ten bins). 25

Table 16. Proportion of peat or non-peat soil. 26

Table 17. Indonesian Land Management Zones 27

Table 18. Comparison of Kappa indices for model runs using single-factor suitability maps. As model 

performance increases, Klocation increases and disagreement for location decreases. 33

Table 19. Weighting schemes for combined-factor GEOMOD runs where all factors were included in the model.

36

Table 20. Ranking of ROC and Kappa scores for multi-factor suitability maps. 37

Table 21. Predicted area of palm oil expansion per Indonesian province between 2000 and 2022 under a control 

case and the additional expansion resulting from a palm-only policy scenario. 39

Table 22. Predicted palm oil expansion per soil type between 2000 and 2022. 41

Table 23. Predicted palm oil expansion per land cover type between 2000 and 2022. 41
SUMMARY

The US Environmental Protection Agency (US EPA) has commissioned Winrock International to conduct an 

analysis of land use change in the Indonesian islands of Sumatra and Kalimantan to predict locations of future 

palm oil expansion until the year 2022. To achieve this goal, we used the land use change simulation model 

GEOMOD. The model evaluates the spatial distribution of existing palm oil plantations relative to various 

biogeophysical factors to predict the most likely locations for future expansion. For this analysis, fourteen 

factors underlying palm oil plantation expansion were evaluated: elevation, precipitation, temperature, slope, 

soil type, land cover type in 2001, distance to roads, distance to rivers, distance to railroads, distance to 

settlements, distance to palm oil mills, peatland soils, Indonesian land management zones, and distance to 

existing plantations. The factors were used both individually and in combination to calibrate the model, and 

the results of each candidate model run were validated against a 2009 map of oil palm plantations. 

The most successful model in predicting expansion, as measured by two statistics, used all fourteen

biogeophysical factors. The final model was run with all factors weighted equally. The quantity of palm oil 

expansion predicted to occur between 2000 and 2022 was simulated under two policy scenarios, which EPA 

calls the “control” case and the “palm-only” case. In one scenario, the “control” case estimates the area of 

expansion assuming that all RFS2 fuels are produced at their mandated levels in 2022, and it does not include 

the use of any palm oil biofuels in the United States. This is then compared to a second scenario, the “Palm 

Only” case which is similar to the control case, but with a change in the amount of land needed to supply 

roughly 400 million additional gallons of palm oil biodiesel and renewable diesel to the United States in 2022. 

For both scenarios, the locations were compared with respect to the area of change in each Indonesian island

as well as the area of change by land cover type and by presence of peat soil, both which affect the magnitude 

of expected GHG emissions.

Approximately 6.2 million hectares of palm oil were expected to exist in Indonesia in 2022, with 63,000 ha of 

additional expansion occurring as a result of the palm-only scenario compared to the control case. Our results 

indicate that nearly 88% of this additional expansion will occur on the island of Sumatra. There are eight 

provinces in Sumatra; the majority of the expansion on Sumatra is predicted to occur in the provinces of Riau

(29%), Sumatera Utara (18%) and Sumatera Selatan (16%). Under the palm-only scenario, palm oil will expand

into histosol soils by 4,894 ha more than under the control case. The vast majority (>90%) of additional 

expansion predicted to occur as a result of the palm-only policy was on two land cover types: forest and

cropland/natural vegetation mosaic.
1.0 BACKGROUND 

Winrock International has been working with US EPA to estimate the lifecycle greenhouse gas (GHG) emissions 

from indirect land use change (iLUC) associated with the increased production of biofuels under the 2007 

Energy Security and Independence Act’s Renewable Fuel Standard (RFS) program. In an initial analysis, 

Winrock estimated the extent of recent land cover change at national and sub-national scales for all countries 

using land cover products derived from moderate resolution satellite imagery, and also developed GHG 

emission factors for various land cover conversions using IPCC methodologies. The U.S. Renewable Fuel 

Standard 2 (RFS2) final rule includes factors for indirect land use change (iLUC) as part of the lifecycle GHG 

calculation of biofuels, which are calculated based on the results of EPA’s modelling which includes the partial 

equilibrium FAPRI-CARD model of international agricultural markets. 

Because palm oil is potentially an important biofuel feedstock, US EPA commissioned Winrock to conduct a 

more detailed analysis of land cover change in Indonesia and Malaysia, where close to 90% of global palm oil is 

produced. Therefore, this analysis supports EPA’s efforts to update the assumptions and data used in 

economic and spatial models to estimate the lifecycle greenhouse gas (GHG) emissions associated with palm 

oil biofuels.

EPA contracted with Integrated Applications Incorporated (IAI), to develop spatial data on the locations of 

existing palm oil plantations for two years, 2000 and 2009 (IAI 2011). The predicted quantity (area) of future 

palm oil expansion comes as output from the FAPRI-CARD model. In this model, a “control” case estimates the 

area of expansion assuming that all fuels under the RFS program are produced at their mandated levels in 

2022. This is then compared to a case similar to the control case, but with a change in the amount of land 

needed to supply roughly 400 million additional gallons of palm oil biodiesel and renewable diesel to the 

United States in 2022.

Given these data inputs, it became possible to use these data as well as various decision rules in a spatial 

model to predict where future expansion is likely to occur across the landscape.

2.0 METHODS

We used information about the locations of existing palm oil plantations on the Indonesian islands of Sumatra 

and Kalimantan in 2000 and 2009 with the GEOMOD land use change simulation model to predict potential 

locations of future palm oil expansion until the year 2022. First, we created maps of factors that could

potentially influence where future expansion occurs, such as elevation, slope, and proximity to roads.

Second, we compared the factor maps against a map of existing palm oil plantations in 2000 to construct a 

series of suitability maps. This is the calibration stage. For each suitability map, the model assigns higher 

suitability values to locations that have a combination of characteristics similar to the land already cultivated in 

palm oil and low suitability values to locations that are less similar to existing palm oil areas. Each resulting 
suitability map shows one possible way in which humans change the landscape, with pixels with the largest 

suitability values being changed first. 

Third, each candidate suitability map was overlain with a map of existing plantations in the year 2009. This is 

the validation stage. The Relative Operating Characteristic (ROC) statistic was used to evaluate whether the 

locations with the highest values in the suitability map were also the locations of actual change to palm oil 

between 2000 and 2009. The best performing suitability map was then used to simulate palm oil expansion 

between 2000 and 2009. The level of agreement between actual and predicted locations of palm oil in 2009

was evaluated using Kappa indices. 

Finally, the suitability map that performed best in predicting the location of change between 2000 and 2009

was chosen to simulate locations of future palm oil expansion across Sumatra and Kalimantan until the year 

2022.

2.1 Model Description

GEOMOD is a spatially explicit simulation model of land cover change that uses maps of biogeophysical 

attributes and maps of existing land cover to extrapolate the known pattern of land cover from one point in 

time to other points in time (Pontius et al. 2001). The goal of GEOMOD is to generate potential scenarios of 

land cover change and evaluate which scenario has the best correlation to real land cover change1; this 

scenario is then used to predict land cover change either forwards or backwards in time. The information 

needed to run GEOMOD includes reference land cover maps containing only two categories (e.g., palm oil vs. 

non-palm oil or forest vs. non-forest) for at least two points in time and information on the quantity of change 

that occurred between the two points in time. Figure 1 shows the components of a GEOMOD simulation. 

GEOMOD has been used to model land cover changes across the world in many different ecosystems such as 

in Costa Rica (Pontius et al. 2001), Indonesia (Harris et al. 2008), India (Rashmi and Lele 2010), the 

Mediterranean (Geri et al. 2011), Massachusetts, USA (Pontius and Schneider 2001), Puerto Rico (Murphy et 

al. 2010), Chile (Echeverria et al. 2008), the Hudson River Basin, USA (Hong et al. 2007), and Bolivia (Hall et al. 

2007). Pontius et al. (2001) give the most complete peer-reviewed description and application of GEOMOD.

 

1 Real land cover change in this case is derived from spatial data on oil palm plantations for 2000 and 2009
Figure 1. Components of the land cover change model GEOMOD. From Pontius and 

Schneider (2001). 

2.2 Data inputs 

The locations of oil palm plantations in the years 2000 and 2009 were supplied in a geographic information

system (GIS) format to EPA by IAI. For information about the maps supplied to Winrock refer to (IAI 2011). 

The 2000 map was used for model calibration and the 2009 map was used for model validation.

Calibrating a model involves the adjustment of model parameters such as the number of strata or how the 

strata bins are created, in order to meet known criteria. The criteria for this analysis was to find the best 
model run using a combination of one or more factors which produced the best possible ROC score and Kappa 

for Location statistic. Model validation is the process of evaluating the model results against a set of 

independent data, the 2009 plantations, to either accept or reject the model. If the model is rejected because 

either or both the ROC score and Kappa is unacceptably low, then the model is re-calibrated and process run 

again. 

Table 1 shows the total area of palm oil in 2000 and the additional areas recorded in 2009 according to the (IAI 

2011) spatial data. Of the 3.9 million hectares of palm oil present in 2009 according to the spatial data, close to 

1.7 million hectares were represented also in the older 2000 spatial data, with the remaining 2.2 million ha 

representing areas of new palm oil expansion between 2000 and 2009. However, approximately 1.5 million ha 

of palm oil areas in 2000 were no longer classified as palm oil by 2009 (Figure 2), indicating a highly dynamic 

landscape. This could be explained in part by the planting cycle for oil palm because the plantations replanted 

between 2000 and 2009 may have been immature in 2000. Immature plantations were not classified as oil 

palm on the 2009 spatial data provided by (IAI 2011). Other explanations may include land abandonment and 

classification errors in (IAI 2011). Assuming that no areas were retired in this time period gives a total area of 

palm oil in 2009 of 5,526,896 ha an increase of 40%. However, at this time to be conservative we operate the 

modeling based on the lower numbers represented by the actual net increase recorded between 2000 and 

2009 (Table 1).

Table 1. Total area (ha) of Sumatra and Kalimantan, Indonesia palm oil plantations 

in 2000 and 2009 according to IAI 2011.

Source 2000 ha 2009 ha 2000-2009 

Increase

Spatial data (IAI 2011) 3,297,543 3,984,520 20%
Figure 2. Representative dynamics of palm oil plantations in Indonesia between 

2000 and 2009. 

Green areas represent areas that were classified as palm oil in both 2000 and 2009 (persistence). Blue areas 

represent areas of palm oil expansion between 2000 and 2009. Yellow areas represent areas that were 

classified as palm oil in 2000 and as non-palm oil in 2009.

Before proceeding with the GEOMOD model runs, we constrained the locations where future palm oil 

expansion could occur in three ways. 

 Plantations could not expand into locations where palm oil already exists, so the existing palm oil 

plantations in 2000 were excluded from the model’s spatial domain. This was done to prohibit the 

model from projecting expansion onto areas that were already replanted and occupied with immature 

oil palm plantations.

 Inland water bodies were excluded due to the obvious reason that palm oil will not grow in these 

locations. 

 Areas designated as protected by the government, according to World Database of Protected Areas 

(WDPA 2010) were also excluded. Figure 3 shows the spatial domain of the non-palm oil areas where 

future expansion could potentially occur.
Figure 3. Candidate locations for palm oil expansion. Areas of expansion were 

assumed to occur only in the non-oil palm plantation areas.

2.3 Creation of Factor Maps

Several different factors can influence the location of future land use change, many of which vary spatially 

over the landscape. The factors used in the GEOMOD analysis were chosen based on both data availability as 

well as past experience; certain factors (e.g., elevation) have been shown in prior analyses (e.g., Harris et al. 

2008 for Indonesia) to correlate well with patterns of land use change.

For this analysis, fourteen factor maps were created: elevation, precipitation, temperature, slope, soil type, 

land cover type in 2001, distance to roads, distance to rivers, distance to railroads, distance to settlements, 

distance to palm mills, peat soil location, land allocation, and distance to existing plantations. Details on each 

factor, including data source and spatial scale, are given in Table 2 and explained in further detail in the 

sections below. For use in GEOMOD, the range of each factor’s values must be stratified into discrete ‘bins’ so 

that categorical values of suitability can be created. For example, elevation data was stratified into 100 m 

‘bins’, with each bin containing pixels with values of 0-100 m, 100-200 m, etc. 
Table 2. Factor data considered in the GEOMOD analysis.

Biogeophysical Attribute Note / Source

Elevation

SRTM 90m digital elevation data v 4

Jarvis, A., H.I. Reuter, A. Nelson, E. Guevara (2008)

Precipitation

WorldClim database v 1.4

Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis (2005)

Temperature

WorldClim database v 1.4

Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis (2005)

Slope Percent slope derived from SRTM

Soil type Harmonized World Soil Database v 1.1

FAO/IIASA/ISRIC/ISSCAS/JRC ( 2009)

Land cover in 2001

MODIS Land Cover Type (MCD12Q1)

These data are distributed by the Land Processes Distributed Active Archive 

Center (LP DAAC), located at the U.S. Geological Survey (USGS) Earth 

Resources Observation and Science (EROS) Center (lpdaac.usgs.gov)

Roads WRI (2008)

Rivers

Global Map v1.0

Source: http://www.iscgm.org/gmd/

Railroads Global Map v1.0

Source: http://www.iscgm.org/gmd/

Settlements Combined 2 datasets, builtupp and mispopp both from Global Map v1.0

Source: http://www.iscgm.org/gmd/

Palm Mills Digitized by Winrock from existing maps created by Wilmar International 

Ltd. 2010.

Peatland Wetlands International (2008)

Land Allocation WRI (2008)

Existing Palm Oil Plantations 2000 IAI (2011)

Protected Areas World Database on Protected Areas

IUCN and UNEP ( 2010) 

Water bodies Global Map v1.0

Source: http://www.iscgm.org/gmd/

Country Boundary GADM database of Global Administrative Areas, version 1.0

http://www.gadm.org/

2.3.1 Elevation 

The elevation data originated from NASA’s Shuttle Radar Topography Mission (SRTM) 90 meter digital 

elevation model (http://srtm.csi.cgiar.org/). Table 3 shows the distribution of land area by elevation across 
the Indonesian islands of Sumatra and Kalimantan and across the areas covered by palm oil plantations. The 

elevation factor map was created by stratifying the range of elevation values into forty 100 meter bins (Figure 

4). 

Table 3. Elevation ranges of palm oil plantations (first ten bins)

Area in each elevation bin as a percentage of area in:

Elevation Bin Bin Range (m) Palm Oil Sumatra and 

Kalimantan 2000 2009

1 0 - 100 94.44% 94.60% 88.52%

2 101 - 200 4.45% 3.87% 3.14%

3 201 - 300 0.48% 0.78% 1.50%

4 301 - 400 0.25% 0.30% 1.05%

5 401 - 500 0.21% 0.17% 0.84%

6 501 - 600 0.13% 0.10% 0.76%

7 601 - 700 0.03% 0.08% 0.66%

8 701 - 800 0.01% 0.07% 0.58%

9 901 - 1000 0.03% 0.53%

10 1001 - 1100 0.01% 0.52%

Figure 4. Stratified factor map for elevation.
2.3.2 Precipitation 

Annual precipitation data originated from the WorldClim global climate dataset2. The precipitation map was

stratified into five suitability classes (0-4) based on Paramananthan (2003) resulting in the stratified 

precipitation map shown in Figure 5. Note that Figure 5 shows only three suitability classes; according to the 

WorldClim data, no land area in either of the two Indonesian islands receives less than 1,540 mm or more than 

4,120 mm of precipitation per year (Table 4).

Table 4. Rainfall suitability classes for growing oil palm. From Paramananthan 

(2003).

Area in each precipitation bin as a 

percentage of area in:

Climatic 

Suitability

Growth 

Limitations

Mean Annual 

Precipitation (mm)

Precipitation 

Bin

Palm Oil Sumatra and 

Kalimantan

2000 2009

Unsuitable Very Severe < 1,250

> 6,000 0 0% 0% 0%

Suitable 1 None 2,501 - 3,500 1 64.53% 64.45% 62.82%

Suitable 1 Slight 1,701 - 2,500

3,501 - 4,000 2 35.18% 35.27% 35.38%

Moderate Moderate 1,451 - 1,700

4,001 - 5,000 3 0.17% 0.17% 1.80%

Marginal Severe 1,250 - 1,450

5,001 - 6,000 4 0% 0% 0%

 

2

http://www.worldclim.org/, Hijmans, R.J., 2005
Figure 5. Stratified precipitation factor map for Indonesia.

2.3.3 Temperature 

Average annual temperature data originated from the WorldClim global climate dataset3. Figure 6 shows the 

distribution of land area by temperature across the Indonesian islands of Sumatra and Kalimantan. Table 5 

shows the distribution of temperature by strata across the two Indonesian islands versus the areas covered by 

palm oil plantations. The temperature map was stratified into five suitability classes (1-5) based on 

Paramananthan (2003) . 

Table 5. Temperature suitability classes for growing oil palm. From Paramananthan 

(2003).

Area in each temperature 

suitability class as a proportion of 

area in:

Climatic

Suitability

Growth

Limitations

Mean Annual

Temperature (C)

Temperature 

Bin

Palm Oil Sumatra and 

Kalimantan

2000 2009

Suitable 1 None 25 – 29 1 99.37% 99.31% 76.63%

Suitable 2 Slight 22 - 25

29 – 32 2 0.51% 0.56% 14.53%

Moderate Moderate 20 - 22

32 – 35 3 0% 0.01% 5.33%

Marginal Severe 10 - 20

35 – 37 4 0% 0% 3.50%

 

3

http://www.worldclim.org/, Hijmans, R.J., 2005
Unsuitable Very Severe < 10

> 37 5 0% 0% 0.002%

Figure 6. Stratified temperature factor map for Indonesia.

2.3.4 Slope 

Slope is the steepness, incline, or grade of a surface. The slope is (in the simplest of terms) defined as the ratio 

of the "rise" divided by the "run". A higher slope value indicates a steeper incline. There are two common ways 

to describe how steep a slope is. One is by the angle in degrees, and the other is in a percentage. For this 

analysis a slope map, calculated as a percentage, was derived from the elevation dataset and was stratified 

into 5% bins. Table 6 shows the distribution of land area by slope across the Indonesian islands of Sumatra 

and Kalimantan and the areas covered by palm oil plantations. The stratified slope map for Indonesia is shown 

in Figure 7.

Table 6. Slope suitability classes for growing oil palm (first ten bins).

Area in each slope suitability bin as a 

percentage of area in:

Slope Bin Bin Range (%) Palm Oil Sumatra and 

Kalimantan

2000 2009

1 0 - 5 98.70% 98.67% 66.92%

2 5 - 10 1.13% 1.10% 10.56%

3 10 - 15 0.13% 0.16% 6.98%

4 15 - 20 0.03% 0.04% 5.23%

5 20 - 25 0.01% 0.01% 3.84%
Area in each slope suitability bin as a 

percentage of area in:

Slope Bin Bin Range (%) Palm Oil Sumatra and 

Kalimantan

2000 2009

6 25 - 30 0% 0.01% 2.64%

7 30 – 35 0% 0% 1.67%

8 35 – 40 0% 0% 0.98%

9 40 – 45 0% 0% 0.54%

10 45 - 50 0% 0% 0.29%

Figure 7. Stratified slope factor map for Indonesia, with slope categorized into 

5% bins. 

2.3.5 Soil Type

Soil type data originated from the Harmonized World Soil Database4and was stratified by grouping soils into 

major soil families (Table 7). In 2009, almost all plantations were located on one of three soil types: Acrisols, 

70%, Histosols, 13%, and Gleysols, 11% (Table 8). 

 The acrisol proportion in plantations is close to the proportion across the country, while the other two

soil types are more common in plantations (Table 8). This indicates that soil type may be a targeted 

consideration in the location of oil palm plantations. Acrisols are clay-rich, relatively infertile soils that 

 

4

http://www.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML, FAO/IIASA/ISRIC/ISSCAS/JRC, 2009
contain toxic levels of aluminum that pose limitations to agricultural use. They are associated with 

humid, tropical climates and often support forested areas. 

 In the USA soil taxonomy, a histosol is an organic soil with an organic carbon content by weight of 12 

percent or more5. Histosols are sometimes referred to as organosols in other countries. When palm 

oil is cultivated on a histosolic soil, the water must be drained as part of site preparation activities 

which results in significant GHG emissions. 

 A Gleysol is a wetland soil that, unless drained, is saturated with groundwater for long enough periods 

to develop a characteristic gleyic colour pattern.

Table 8 shows the distribution of land area by soil strata across the Indonesian islands of Sumatra and 

Kalimantan and across the areas covered by palm oil plantations. The stratified soil type factor map for 

Indonesia is shown in Figure 8.

Table 7. Soil types in Indonesia and the reclassification scheme used to create 

the soil type factor map.

Soil Name Reclass Name Reclass ID

Humic Acrisols Acrisols 1

Ferric Acrisols Acrisols 1

Ferric Acrisols Acrisols 1

Orthic Acrisols Acrisols 1

Eutric Gleysols Gleysols 2

Lithosols Lithosols 3

Eutric Regosols Regosols 4

Orthic Acrisols Acrisols 1

Orthic Acrisols Acrisols 1

Orthic Acrisols Acrisols 1

Ferralic Cambisols Cambisols 5

Rhodic Ferralsols Ferralsols 6

Dystric Gleysols Gleysols 2

Dystric Gleysols Gleysols 2

Eutric Fluvisols Fluvisols 7

Thionic Fluvisols Fluvisols 7

Thionic Fluvisols Fluvisols 7

Calcic Luvisols Luvisols 8

Orthic Luvisols Luvisols 8

Dystric Nitosols Nitosols 9

Dystric Histosols Histosols 10

Dystric Histosols Histosols 10

Dystric Histosols Histosols 10

Gleyic Podzols Podzols 11

 

5

http://soils.usda.gov/technical/classification/taxonomy/
Soil Name Reclass Name Reclass ID

Albic Arenosols Arenosols 12

Table 8. Soil types for growing palm oil in Indonesia.

Area in each soil type as a percentage 

of area in:

Reclass Name Palm Oil Sumatra and 

Kalimantan

2000 2009

Acrisols 48.76% 33.65% 46.94%

Gleysols 17.36% 21.88% 11.31%

Lithosols 11.83% 13.62% 11.20%

Cambisols 2.84% 2.71% 6.42%

Ferralsols 1.34% 1.97% 3.97%

Fluvisols 6.63% 6.44% 3.70%

Luvisols 2.13% 6.27% 3.43%

Nitosols 0.29% 0.10% 2.98%

Histosols 8.76% 13.28% 10.06%

Figure 8. Stratified soil factor map for Indonesia.
2.3.6 Land Cover in 2001

Land cover data originated from the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover 

Type version 5 for the year 2001.

6 MODIS is a sensor that is on two satellites, the Terra (EOS AM) and Aqua 

(EOS PM) satellites operated by NASA. The Land Cover dataset is produced on an annual basis starting in 2000. 

The dataset uses the International Geosphere Biosphere Programme (IGBP) global vegetation classification 

scheme which contains 17 categories of land cover. The primary land cover scheme includes 11 natural 

vegetation classes, 3 developed and mosaiced land classes, and three non-vegetated land classes. Table 9

shows the distribution of land area by land cover across the Indonesian islands of Sumatra and Kalimantan and 

across the areas covered by palm oil plantations. The year 2001 was chosen to represent the land cover 

condition of areas before conversion to oil palm between 2000 and 2009. Data were stratified by aggregating 

the most common land cover types into a Boolean map of two categories (Figure 9) by assigning broadleaf 

evergreen forest and natural mosaic classes a value of 1 because these two categories account for over 86% of 

the land cover where existing palm oil plantations occur. All other land cover categories were assigned a 

value of 0. 

Table 9. Land cover area (ha) based on MODIS 2001 land cover for palm oilproduct.

Area in each land cover type as a percentage 

of area in:

MODIS 2001

Palm Oil Sumatra and 

Kalimantan

2000 2009

Evergreen Needleleaf Forest 0.01% 0.00% 0.05%

Evergreen Broadleaf Forest 51.52% 57.53% 73.94%

Deciduous Needleleaf Forest 0.00% 0.00% 0.03%

Deciduous Broadleaf Forest 0.16% 0.24% 0.33%

Mixed Forests 0.97% 0.83% 0.65%

Closed Shrublands 0.02% 0.03% 0.08%

Open Shrublands 0.06% 0.07% 0.04%

Woody Savannas 1.41% 1.47% 0.88%

Savannas 4.99% 4.30% 1.24%

Grasslands 0.95% 0.94% 0.38%

Permanent Wetlands 0.70% 2.10% 4.27%

Croplands 4.03% 2.88% 2.84%

Urban And Built-Up 0.13% 0.14% 0.22%

Cropland/Natural Vegetation Mosaic 34.96% 29.33% 14.98%

Barren Or Sparsely Vegetated 0.01% 0.03% 0.07%

 

6

https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/land_cover/yearly_l3_global_500_m/mcd12

q1
Figure 9. Stratified 2001 land cover factor map for the Indonesian islands of 

Sumatra and Kalimantan.

2.3.7 Infrastructure (Roads, Rivers, Railroads, Settlements)

Roads, rivers, railroads and settlements data originated from Global Map version 17. The datasets were

created and are maintained by the International Steering Committee for Global Mapping (ISCGM). Each data 

layer was stratified into Euclidean distances where each class width is 1000 meters (e.g. category 1= 0-1000m; 

category 2 = 1000-2000m, etc). This resulted in “distance from” factors for the model ( i.e. distance from 

roads, distance from rivers, etc.) (Figure 10), with low values close to the feature and high values far from the 

feature.

Table 10. Distance from roads (first ten bins)for palm oil plantations.

Area in each distance from roads bin as a 

percentage of area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil Sumatra and Kalimantan

2000 2009

1 0 - 1 79.4% 64.2% 36.6%

2 1 - 2 10.9% 14.1% 17.1%

3 2 - 3 4.1% 7.3% 10.0%

4 3 - 4 2.0% 4.5% 6.9%

5 4 - 5 1.2% 3.0% 4.9%

 

7

http://www.iscgm.org/gmd/
6 5 - 6 0.8% 2.1% 3.8%

7 6 - 7 0.5% 1.5% 3.4%

8 7 - 8 0.3% 0.8% 2.6%

9 8 - 9 0.2% 0.6% 2.1%

10 9 - 10 0.1% 0.3% 1.7%

Table 11. Distance from rivers (first ten bins) for palm oil plantations.

Area in each distance from rivers bin as a 

percentage of area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil

Sumatra and Kalimantan

2000 2009

1 0 - 1 11.8% 10.6% 19.1%

2 1 - 2 14.7% 13.5% 18.1%

3 2 - 3 14.0% 13.3% 14.2%

4 3 - 4 12.4% 12.1% 11.2%

5 4 - 5 10.0% 10.1% 8.2%

6 5 - 6 8.0% 8.3% 6.2%

7 6 - 7 7.0% 7.7% 5.2%

8 7 - 8 5.0% 5.7% 3.8%

9 8 - 9 3.8% 4.2% 2.8%

10 9 - 10 2.8% 3.1% 2.0%

Table 12. Distance from railroads (first ten bins) for palm oil plantations.

Area in each distance from railroads bin as a 

percentage of area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil

Sumatra and Kalimantan

2000 2009

1 0 - 1 1.2% 1.1% 0.5%

2 1 - 2 1.4% 1.2% 0.5%

3 2 - 3 1.3% 1.0% 0.4%

4 3 - 4 1.1% 1.0% 0.4%

5 4 - 5 1.1% 1.0% 0.4%

6 5 - 6 1.1% 0.9% 0.4%

7 6 - 7 1.3% 1.0% 0.4%

8 7 - 8 1.2% 0.9% 0.4%

9 8 - 9 1.2% 0.9% 0.4%

10 9 - 10 1.0% 0.8% 0.4%

Table 13. Distance from settlements (first ten bins) for palm oil plantations.
Area in each distance from settlements bin as a 

percentage of area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil

Sumatra and Kalimantan

2000 2009

1 0 - 1 0.2% 0.2% 0.3%

2 1 - 2 0.9% 0.8% 1.2%

3 2 - 3 1.6% 1.5% 1.9%

4 3 - 4 2.4% 2.1% 2.5%

5 4 - 5 3.2% 2.8% 2.9%

6 5 - 6 4.1% 3.7% 3.4%

7 6 - 7 5.4% 4.8% 4.3%

8 7 - 8 6.0% 5.4% 4.7%

9 8 - 9 6.0% 5.4% 4.7%

10 9 - 10 5.7% 5.2% 4.5%

A) B)

C) D)

Figure 10. Stratified factor map for (A) roads; (B) rivers; (C) railroads; and (D) 

settlements for Indonesia.
2.3.8 Distance from Existing Plantations

The distance from 2000 oil palm plantations data originated from the 2000 data on oil palm plantations. 

First, we examined the distances between 2000 plantations and new 2009 plantations and found that 60% of 

new plantations were located within a kilometer of a 2000 plantation. This would indicate that new

plantations are being created near existing plantations, likely in order to leverage existing infrastructure. 

The plantations data were stratified into Euclidean distances where each class width is 1000 meters (e.g. 

category 1= 0-1000m; category 2 = 1000-2000m, etc). This results in a distance from plantations factor for the 

model (Figure 11). Table 14 shows the distribution of land area by distance from existing oil palm plantations 

across the Indonesian islands of Sumatra and Kalimantan and across the area covered by palm oil plantations.

Table 14. Distance from 2000 oil palm plantations (first ten bins).

Area in ‘distance from 2000 palm oil plantation’ bins as a 

percentage of area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil

Sumatra and Kalimantan

2000 2009

1 0 - 1 100% 52.2% 5.2%

2 1 - 2 5.0% 2.1%

3 2 - 3 3.3% 1.8%

4 3 - 4 2.9% 1.7%

5 4 - 5 2.5% 1.6%

6 5 - 6 2.4% 1.6%

7 6 - 7 2.6% 1.7%

8 7 - 8 2.4% 1.6%

9 8 - 9 2.1% 1.6%

10 9 - 10 1.8% 1.5%
Figure 11. Stratified distance from plantations factor map for Indonesia.

2.3.9 Distance from Palm Oil Mills

Palm oil mill locations were created by digitizing a series of map documents provided to the EPA by Wilmar 

International Ltd. which showed the locations of palm oil mills across Sumatra and Kalimantan (Wilmar 2010). 

In total there were 557 mills as of August 2010, according to the data provided. The majority (445) are located 

in Sumatra. The density of palm oil mills is one mill for every 96,839 hectares in Sumatra. In Kalimantan the 

density is one mill for every 470,611 hectares. 

The map of palm oil mills was stratified into Euclidean distances where each class width is 1000 meters (e.g. 

category 1= 0-1000m; category 2 = 1000-2000m, etc). This results in a distance from palm oil mills factor for 

the model (Table 15; Figure 12). Sixty-two percent of plantations in 2000 and 56% of plantations in 2009 were 

within 10km of a palm oil mill. This compares to just 12% of the area of Indonesia that is within the same 

distance.

Table 15. Distance from palm oil mills (first ten bins).

Area in ‘distance from palm oil mills’ bins as a percentage of 

area in:

‘Distance 

from ‘

Bin

‘Distance 

from’

Range (km)

Palm Oil

Sumatra and Kalimantan

2000 2009

1 0 - 1 1.5% 1.1% 0.2%
2 1 - 2 5.1% 4.1% 0.6%

3 2 - 3 7.1% 5.9% 0.9%

4 3 - 4 8.2% 7.0% 1.1%

5 4 - 5 7.7% 6.8% 1.2%

6 5 - 6 7.6% 6.8% 1.4%

7 6 - 7 7.7% 7.1% 1.6%

8 7 - 8 6.8% 6.6% 1.7%

9 8 - 9 5.6% 5.9% 1.7%

10 9 - 10 4.6% 5.1% 1.6%

Figure 12. Stratified distance from palm mills.

2.3.10 Peat Soils

A dataset of detailed peat soil characteristics was used as an additional factor. The data originated from 

Wetlands International and came to Winrock via World Resources Institute (WRI). The data originally 

described peat soil depth but for this study was reclassified into peat or non-peat soils (Figure 13). The 

majority of existing oil palm plantations occurred on non-peat soil, though that percentage decreased 

between 2000 to 2009 (Table 16). The area of plantations on peat soils in 2009 was over 25% compared with 

only 13% of soils in this classification across Indonesia.

Table 16. Proportion of peat or non-peat soil.

Area of each stratum as a percentage of 

area in:
Strata 

Description

Strata

Class

Palm Oil Sumatra and 

Kalimantan

2000 2009

Non-Peat Soil 1 82.17% 74.74% 86.83%

Peat Soil 2 17.76% 25.16% 13.17%

Figure 13. Stratified peat soil factor map.

2.3.11 Land Management Zones

Land management zones are designated and regulated by the government of Indonesia. The data were 

provided to Winrock by World Resources Institute (WRI). The data represents land management zones in the 

year 2008. Zones that are conducive to conversion from forest to plantation (conversion & production) were 

aggregated into one stratum; all others (conservation, protection, non-forest area and no-data) were not 

included in the model. Table 17 shows the distribution of land area by management zone strata across the 

Indonesian islands of Kalimantan and Sumatra and across the area covered by palm oil plantations.

Table 17. Indonesian Land Management Zones

Area of each management zone as a percentage 

of area in:

Reclass ID

Management 

Zone

Palm Oil Sumatra and 

Kalimantan

2000 2009

0 Conservation 0.89% 0.72% 15.65%

1 Conversion 19.14% 23.23% 10.55%

0 No Data 0.34% 0.25% 0.39%

0 Non-forest area 51.16% 44.07% 30.92%

1 Production 28.33% 31.54% 37.04%

0 Protection 0.13% 0.19% 5.45%
Figure 14. Management Zones reclassified

2.3 Creating Suitability-for-Change Maps

After factor maps were created, they were used individually and in combination to create a series of suitability 

maps. Each model run uses one candidate suitability map to simulate future change, placing change in cells 

that have the highest suitability values. If the suitability map were perfect, the order of the suitability values 

would match the order in which humans change the landscape, with the largest suitability values being 

changed first (Pontius and Schneider 2001). For validation, the model run result, a map of simulated future 

change, is compared against the real world map of 2009 plantations. Several model runs can generate a series 

of simulated maps, each with a different location of change, hence each with a different level of agreement 

with the map of reality.

In the simplest case, a single factor such as elevation can be used to create a suitability map. In this case, for 

example, a pixel that falls in the 0-100m bin of the stratified elevation map is compared to the percentage of 

all pixels that fall into the same 0-100 m bin in locations of existing palm oil. The suitability map is generated 

by assigning the highest suitability score to the pixels that have an elevation value that is found most 

commonly in the existing plantations, and the lowest suitability score is assigned to the pixels that have an 

elevation value that is found least commonly in the existing plantations. The result is a suitability map where 

reds (high values) are the most suitable areas and blues (low values) are the least suitable (Figure 15).
Figure 15. Example suitability map created from one or several factor maps, where 

blue is least suitable and red is most suitable.

A suitability map can also be created by combining multiple individual factor maps together and weighting 

each factor equally or unequally. The weights for the factor maps are applied when GEOMOD creates the 

suitability map by summing the weights for each grid cell from all the factor maps. Hence, the suitability in 

each cell is calculated according as:

 ( ) [∑* 

( )+

 

 

] [∑ 

 

 

]

Where:

R(i) = suitability value in cell(i),

a = a particular factor map,

A = the number of factor maps,

Wa = the weight of factor map a, and

Pa(i) = percent-developed in category ak of factor map a, where cell(i) is a member of category ak.

For this analysis, 34 different suitability maps were created. In the first stage, we created fourteen suitability 

maps, one from each single factor, to test the performance of a single factor in predicting the location of palm 

oil expansion between 2000 and 2009. After single factors were tested, we used this information to create 

suitability maps from several combined factors and using different weighting schemes. Results from the 

validation procedures of each scenario were used to refine subsequent simulations until a simulation model 
achieved the best validation results possible given available time and resources. When the best model was 

identified based on the validation exercises, we used this model to simulate expansion of oil palm plantations 

from 2000 to 2022.

2.4 Evaluating Suitability Maps using the ROC statistic

The validity of each suitability map in predicting change was tested using a statistic called the Relative 

Operating Characteristic (ROC), which has been used widely to assess model accuracy (Kim 2010, Pontius et al.

2001, Manson 2002). The ROC is a quantitative score that helps to evaluate the ability of a suitability map to 

distribute change accurately across the landscape. It results from a comparison of the suitability map with a 

reference map and is a way to assess, with a simple figure, the quality of a suitability map by showing clearly 

how agreement differs across a variety of change scenarios. The ROC score is a better indicator of map 

accuracy than a simple ‘percent correct’ statistic because in many cases, the ‘percent correct’ value may be 

high simply because the model predicts no change where there is no change. The more important question in 

most cases is whether the model correctly predicts change where there is change in reality. 

The first step in computing an ROC score for a suitability map is to reclassify it into a map with several 

suitability percentile groups. For example, a possible reclassification rule could assign 10% of the cells that 

have the largest suitability values to group 1, another 10% of the cells that have the next largest suitability 

values to group 2, etc. Next, the modeler overlays the map of suitability groups with the map of reality8. A 

tabular result shows the proportion of grid cells of each suitability group that are categorized correctly as 

change vs. non-change in reality. This exercise is an efficient way to show results for 11 scenarios of change 

using the same map of reality and the same map of suitability; however, the quantity of change in each 

scenario is different. In the first scenario, the quantity of change is zero and therefore none of the 10 

suitability groups are classified as change in that scenario. In the second scenario, change is simulated only in 

suitability group 1 (top 10% of suitability values), the next scenario simulates change in both suitability groups 

1 and 2 (top 20% of suitability values), and so on until the last scenario simulates change in all 10 groups. For 

each scenario, the percent of grid cells that each scenario classifies correctly as non-change and change is 

calculated. If a grid cell is simulated correctly as change in a model scenario, it is a ‘positive’. Therefore, a ‘true 

positive’ is a cell that is categorized as change in both reality and the modeled scenario. A ‘false-positive’ is a 

cell that is incorrectly categorized in the model as change but that is non-change in reality. 

The ROC statistic is calculated as the area under the blue curve that connects the plotted points for percent 

‘true-positive’ and ‘false-positive’ (Figure 16). If the sequence of the suitability values matches perfectly the 

sequence in which real land-cover change has occurred, then ROC would be equal to 1, corresponding to the 

area under the purple curve, because as the amount of modeled change increases from 0% to 100% 

(horizontal axis), there are no (0%) false positives and all (100%) true positives. If the order of suitability values 

 

8

For this analysis – map of reality is the map of oil palm plantations in 2009
were assigned at random locations across the landscape, then the expected value of the ROC would be 0.5, as 

shown by the red line in Figure 16. 

Figure 16. Example ROC curve used to validate suitability maps for simulating land 

use change in Indonesia.

What constitutes a “good” ROC? Any ROC above 50% (red line) is better than random. In other fields, ROCs 

have ranged from 71 to 89% for weather forecasting, 75 to 97% for library information retrieval, 81 to 93% for 

medical imaging diagnosis, 68 to 93% for material strength testing, and 55 to 98% for polygraph lie detection 

(Pontius and Schneider 2001). 

2.5 Evaluating Model Performance using the Kappa statistic

The ROC statistic assesses how well a candidate suitability map assigns locations of land cover change under 

different scenarios. Kappa indices of agreement are a collection of seven related statistical measures used to 

assess how well model predictions agree with reality shown on the spatial dataset when a single quantity of 

change is simulated. In other words, the ROC score is an evaluation of the suitability map, and the Kappa 

indices are an evaluation of the prediction map that results when the suitability map is used to predict real

land cover change. 

It has become customary in the remote sensing field to report the standard kappa index of agreement,

especially for purposes of accuracy assessment, because it is generally thought to be a more robust measure 

than simple percent correct since kappa takes into account random agreement. The Kappa is usually 

attributable to Cohen (1960), and is sometimes referred to as Cohen’s kappa, but kappa can be traced to 

Galton (1892). The use of kappa indices in map comparison was first proposed by Congalton et al. (1983), later 
advocated by Monserud and Leemans (1992), Congalton and Green (1999), Smits et al. (1999), and Wilkinson 

(2005).

For all Kappa statistics, 0 indicates that the level of agreement is equal to the agreement due to chance

(randomness) and 1 indicates perfect agreement. The kappa indices of agreement related to location are:

 Kappa for no information (Kno): indicates overall agreement between the prediction map and the 

reality map relative to random agreement. 

 Kappa for location (Klocation): indicates the extent to which the two maps agree in terms of location 

of each category (in this case, palm oil vs. non-palm oil).

Kappa indices are also expressed as percentages indicating how much each index contributes to the result of 

the model simulation:

 Agreement due to chance (always 50% in two-category scenarios)

 Agreement due to quantity

 Disagreement due to quantity (always 0% for this analysis, because quantity is predetermined)

 Agreement due to location: indicates the percent correct in terms of location between prediction map 

and map of reality

 Disagreement due to location: indicates the additional percent correct obtainable if the prediction

map were to have perfect location given no change in its quantities of each category. 

For this analysis, the important kappa indices are the ones concerned with location. Since the change quantity 

was predetermined, the statistics for quantity are irrelevant and are not discussed in detail here. The higher 

the Kappa for location and the lower the disagreement due to the location, the better the model performed at 

predicting correctly the locations of land cover change.

3.0 RESULTS AND DISCUSSION

A single model showed the best predictions of palm oil expansion in Indonesia between 2000 and 2009. It

incorporated all fourteen biogeophysical factors and weighted all factors equally with no neighborhood

constraint. The suitability map for the model had an ROC score of 0.88, a better score than any other singlefactor suitability map and tied with one other combination-factor suitability map. The Knoscore of 0.90, and 

Kappa for location scores of 0.46, were the best Kappa scores of all models runs including single-factor and 

combination-factor runs. In the final model, disagreement due to location was 4.9%. All of these metrics 

indicate good model performance. The resulting suitability map was used to simulate locations of future palm 

oil expansion to the year 2022 under two policy scenarios: a control case and palm only case. In each scenario, 

the quantity of palm oil expansion was predetermined using output from the FAPRI-CARD model. In both 

model variations, nearly 88% of the additional expansion that occurred as a result of the palm-only scenario 

occurred on the island of Sumatra. Under the palm-only scenario, palm oil expanded into histosol soil areas by 

4,894 ha more than under the control case, depending on which model variation was selected. The vast 
majority (79%) of additional expansion predicted to occur as a result of the palm-only policy was on two land 

cover types: forest and cropland/natural vegetation mosaic.

3.1 Evaluation of Single-Factor Suitability Maps for Predicting Change

The results of fourteen model runs, each which used a single biogeophysical factor to create a suitability map,

offered insight into how well each factor individually could predict the expansion of oil palm. The factor with 

the best ROC score was distance from existing palm oil mills, followed by distance from existing plantations 

and elevation (Figure 17). This result was to be expected, as palm oil mills are dependent on being in close 

proximity to palm plantations. The factors of elevation and slope were nearly identical in both their ROC and 

Kappa scores, which was not surprising as slope is associated with elevation and there was little variability in 

slope across existing oil palm areas. 

The distance from roads and temperature factor maps were the next best factors in predicting palm oil 

expansion. Distance from railroads and distance from rivers were also good explanatory variables for oil palm 

expansion. Precipitation, land cover type, distance from settlements and land management were poor 

predictors of expansion when used individually, likely because these factors do not vary much across the 

landscape. Overall, the factors lined up in similar rankings when ranked either by ROC or by Kappa statistics 

(Figure 17, Error! Reference source not found.). 

For some factors with the highest ROC scores (e.g., elevation, slope), Figure 18 shows that the factor 

performed worse than randomly for the first 10% of modelled change in the landscape. One reason for this is 

the fact that some plantations exist in locations that are outliers compared to where most other existing 

plantations occur. Figure 19 shows 2009 plantations overlain onto an elevation suitability map and shows 

clearly that several plantations exist in the lowest suitability areas for the factor. GEOMOD would not allocate 

change in these areas until all of the other higher suitability areas have transitioned. If the number of high 

suitability pixels exceeds the quantity of pixels predicted to change, GEOMOD would miss these locations of 

change. This effect is seen in the graphs for the initial 10% of change (Figure 18). As a result, it is impossible to 

achieve perfect location of change in the validation stage as long as these plantations are located in areas of 

low suitability. Thus, ROC is a valuable tool for to answer the question ‘How well is the new oil palm

concentrated at the locations of relatively high suitability?’ while Kappa indices are useful for answering the 

question ‘How well did my model do in predicting locations of new oil palm compared to actual locations of 

new oil palm?’ 

Table 18. Comparison of Kappa indices for model runs using single-factor 

suitability maps. As model performance increases, Klocation increases and 

disagreement for location decreases. 

Factor Kno Klocation Agreement Due to 

Location

Disagreement Due 

to Location

Elevation 0.896 0.424 3.848 5.224
Factor Kno Klocation Agreement Due to 

Location

Disagreement Due 

to Location

Landcover 2001 0.895 0.421 3.822 5.249

Land Management 0.895 0.422 3.828 5.244

Distance from Palm Mills 0.899 0.442 4.013 5.059

Peatland 0.896 0.425 3.857 5.214

Precipitation 0.895 0.421 3.817 5.255

Distance from Rail 0.896 0.424 3.849 5.222

Distance from Rivers 0.895 0.422 3.828 5.244

Distance from Roads 0.895 0.423 3.836 5.236

Distance from Settlements 0.895 0.422 3.825 5.247

Slope 0.895 0.423 3.839 5.233

Soil Type 0.895 0.423 3.834 5.238

Temperature 0.895 0.424 3.843 5.229

Distance from 2000 Plantations 0.900 0.448 4.063 5.008

Figure 17. Ranking of ROC scores.

Dist. From Palm Mills ROC = 0.842 Dist. From Plantations ROC = 0.724 Elevation ROC = 0.697
Slope ROC = 0.670 Dist. From Roads ROC = 0.646 Temperature ROC = 0.618

Dist. From Railroads ROC = 0.612 Dist. From Rivers ROC = 0.603 Soil Type ROC = 0.563

Peatlands ROC = 0.554 Dist. From Settlements ROC = 0.549 Land Allocation ROC = 0.525

Precipitation ROC = 0.520 Land Cover Type 2001 ROC = 0.501

Figure 18. Ranking of single-factor suitability map ROC scores.
Figure 19. Locations of 2009 palm oil plantations located in areas of low 

suitability with respect to the elevation factor map.

3.2 Evaluation of Multi-Factor Suitability Maps for Predicting Change

After running single-factor scenarios to evaluate which factors performed best in predicting plantation 

expansion, we attempted to improve model performance further by generating suitability maps that included

multiple factors. First, we created a suitability map from the combination of all fourteen factors where each 

factor was weighted equally. Second, we created a suitability map from the two factors with the best ROC 

scores, distance from palm mills and distance from plantations, equally weighted. Third, we added the next 

two strongest ROC scoring factors, elevation and distance from roads. Next different stratification methods 

were tried on all the ‘distance from’ factors. This involved using alternative methods to group the distance

values into bins. The methods were 500 meter equal intervals instead of 1000, quantile which limits the 

variance between bins, and standard deviation groupings. The newly stratified single factor was run 

individually in GEOMOD. If it performed better than the original stratification it was used in as a replacement 

in a combination run. 

After all the trials, the combined 14-factor suitability map showed the most significant improvement over any 

single-factor map and over any combination map, achieving ROC scores of 0.879 for an equal-weighted 

scenario. The other combination scenarios also performed better than single-factor scenarios, but to a lesser 

extent (Table 20).

Table 19. Weighting schemes for combined-factor GEOMOD runs where all factors were 

included in the model.
Factor Equal Weight 

Scenario

Elevation 0.0714

Land cover data in 2001 0.0714

Land management 0.0714

Distance from Palm Mills 0.0714

Peatland 0.0714

Precipitation 0.0714

Distance from Rail 0.0714

Distance from Rivers 0.0714

Distance from Roads 0.0714

Distance from Settlements 0.0714

Slope 0.0714

Soil order 0.0714

Temperature 0.0714

Distance from 2000 Plantations 0.0714

1.0000

Table 20. Ranking of ROC and Kappa scores for multi-factor suitability maps. 

Weighting Scenario ROC 

Score Kno Klocation

Agreement 

due to 

location

Disagreement 

due to 

location

All Factors (equal weights) 0.879 0.901 0.457 4.145 4.926

Two Best ROC scores 0.879 0.901 0.455 4.124 4.947

Four Best ROC scores 0.876 0.901 0.457 4.144 4.928

Four Best ROC scores (alternative 

stratification) 0.856 0.901 0.456 4.140 4.932

5 Best ROC scores (alternative stratification) 0.872 0.901 0.457 4.142 4.929

All Factors (alternative stratifications) 0.873 0.901 0.457 4.142 4.929

All Factors (alternative stratifications) 0.877 0.901 0.457 4.142 4.929

3.3 Final Suitability Map

The suitability map shown in Figure 20 corresponds to the scenario of all factors included and equally 

weighted; this scenario achieved the highest combination of, ROC score (0.879), Kappa for no information 

(Kno) (0.901), Kappa for location (Klocation) (0.46), agreement due to location (4.15), and disagreement due to 

location (4.93). Zooming in to specific areas of the final prediction map (Figure 21) shows slight variations in 

where additional palm oil expansion is predicted to occur depending on which case outcome is used in the 

simulation. 
Figure 20. Suitability map that incorporates all 14 factors equally weighted. 

Figure 21. Predicted expansion of palm oil from 2000-2022 in the equal weight 

scenario. 

Green pixels show predicted expansion in both control and palm-only scenarios, red pixels show additional 

expansion in the palm-only scenario.
3.4 Palm Oil Expansion, 2000-2022

3.4.1 Expansion by State

The area of predicted palm oil expansion between 2000 and 2022 per Indonesian Province is shown in Table 

21. Of the 63,453 ha of additional expansion predicted as a result of a palm-only policy in 2022, nearly 88% is 

expected to occur on the island of Sumatra, with a majority of that occurring in the three provinces of Riau

(29%), Selatan (18%) and Utara (16%) (Figure 22). 

Table 21. Predicted area of palm oil expansion per Indonesian province between 

2000 and 2022 under a control case and the additional expansion resulting from a 

palm-only policy scenario.

Island Province Control 

Case (ha)

Additional Expansion,

Palm-Only Case (ha)

Kalimantan Kalimantan Barat 350,559 1,610

Kalimantan Kalimantan Selatan 342,488 1,374

Kalimantan Kalimantan Tengah 496,441 3,113

Kalimantan Kalimantan Timur 206,695 1,438

Sumatra Aceh 246,020 2,833

Sumatra Bengkulu 112,245 1,739

Sumatra Jambi 295,585 6,526

Sumatra Lampung 262,012 2,469

Sumatra Riau 1,511,197 18,675

Sumatra Sumatera Barat 205,407 2,318

Sumatra Sumatera Selatan 814,737 11,313

Sumatra Sumatera Utara 1,334,833 10,046

6,178,220 63,453

Figure 22. Distribution of palm oil expansion among Indonesia provinces.
3.4.2 Expansion by Soil Type

The area of predicted palm oil expansion between 2000 and 2022 per soil type is important primarily for 

determining how much expansion occurs on histosols. In both the FAO soil classification and the USA soil 

taxonomy, a histosol is an organic soil with an organic carbon content by weight of 12 percent or more.

Histosols are sometimes referred to as organosols. When palm oil is cultivated on a histosolic soil, the water

must be drained as part of site preparation activities and results in significant GHG emissions. Approximately 

9% (531,366 ha) of palm oil expansion are predicted to occur in histosol soils with an additional 4,894 ha in the 

palm only model. The majority of the expansion occurs on acrisols that are found commonly throughout 

Indonesia (Table 22).
Table 22. Predicted palm oil expansion per soil type between 2000 and 2022. 

Soil Type Control 

Case (ha)

Additional Expansion, 

Palm-Only 

Case (ha)

Acrisols 3,095,142 33,487

Gleysols 926,102 8,479

Lithosols 855,243 10,411

Histosols 531,366 4,894

Fluvisols 354,165 2,855

Cambisols 193,794 2,340

Luvisols 106,149 236

Ferralsols 91,123 751

Nitosols 20,393 0

6,173,476 63,453

3.4.3 Expansion by Land Cover Type

The area of predicted palm oil expansion between 2000 and 2022 per land cover type is important primarily 

for determining emission factors for various land use transitions. When forest land is cleared for oil palm, 

significantly more GHG emissions result than when grassland is replaced with palm oil, for example. When 

locations of predicted palm oil expansion between 2000 and 2022 were overlain with MODIS land cover type 

for the year 2001, most predicted expansion of palm oil between 2000 and 2022 occurred in the forest and 

mixed (cropland/natural vegetation mosaic) land cover types (Table 23).

Table 23. Predicted palm oil expansion per land cover type between 2000 and 2022. 

Land cover type Control Case (ha)

Additional Expansion, 

Palm-Only 

Case (ha)

Evergreen Broadleaf Forest 2,937,218 26,489

Cropland/Natural Vegetation Mosaic 2,249,086 23,956

Savannas 376,018 5,753

Croplands 334,631 4,379

Woody Savannas 91,488 880

Mixed Forests 53,214 472

Permanent Wetlands 53,107 601

Grasslands 52,398 623

Urban And Built-Up 13,137 64

Deciduous Broadleaf Forest 8,715 86

Open Shrublands 3,220 21

Water 2,576 86

Closed Shrublands 1,782 21

Barren Or Sparsely Vegetated 1,052 0

Evergreen Needleleaf Forest 429 0

Deciduous Needleleaf Forest 150 21

6,178,220 63,453
4.0 CONCLUSIONS

The GEOMOD model was used to predict locations of future expansion of palm oil across Sumatra and 

Kalimantan, Indonesia until the year 2022. Fourteen spatial factors were tested for their ability to predict the 

distribution of plantation expansion between 2000 and 2009 and were combined to create a suitability map. 

The best single factor suitability map had an ROC score of 0.842. The best combined factor suitability map had 

an ROC score of 0.879 indicating that most new palm oil plantations in 2009 were concentrated in locations 

that were assigned high suitability scores. This suitability map was then used to simulate locations of 

expansion between 2000 and 2009, and the simulated results were compared (validated) with actual 2009

palm oil locations. The simulation resulted in an overall Kappa for no information (Kno) statistic of 0.9 with a 

Kappa for location of 0.46, indicating good model performance. After this model validation step, the model 

was run to the year 2022 to predict locations of future palm oil expansion under a control case and a palmonly case. Additional expansion due to the palm-only scenario occurred primarily on Sumatra, in the provinces

of Riau, Utara and Selatan, and 8% of the additional expansion occurred on histosolic soils. 

The largest source of uncertainty in the analysis was the quality of the input data for palm oil plantation

locations. A comparison of the 2000 vs. 2009 spatial data suggests highly dynamic change between the two 

years. Although there was a net expansion of 686,963 ha, there was only 32% persistence of plantations during 

that time. The result was that the areas modeled as suitable for transition were not always the same areas 

that were converted by 2009, which affected the validation of the model by increasing the disagreement due 

to location. 

 
5.0 REFERENCES

Cohen, J. 1960. A coefficient of agreement for nominal scales. Educational and Psychological Measurement, 

20, pp. 37-46.

Congalton, R.G., Oderwald, R.G. and Mead, R.A. 1983. Assessing Landsat classification accuracy using discrete 

multivariate analysis statistical techniques. Photogrammetric Engineering and Remote Sensing, 49, pp. 1671-

1678.

Congalton, R.G. and Green, K. 1999. Assessing the accuracy of remotely sensed data: principles and practices, 

pp. 183 (Boca Raton, FL: Leweis)

Echeverria, C., Coomes, D., Hall, M., Newton, A. 2008. The application of GIS-based modeling techniques to 

analyze the patterns of forest loss and fragmentation between 1976 and 2020 in southern Chile. Ecological 

Modeling 212: 439-449.

FAO/IIASA/ISRIC/ISSCAS/JRC. 2009. Harmonized World Soil Database (version 1.1). FAO, Rome, Italy and IIASA, 

Laxenburg, Austria. http://www.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/

Galton, F. 1892. Finger Prints, pp. 216 (London: Macmillan)

Geri, F., Amici, Valerio., and Rocchini, D.2011. Spatially-based accuracy assessment of forestation prediction in 

a complex Mediterranean landscape, Applied Geography, Volume 31, Issue 3, July 2011, p.881-890.

Hall, M.H., Dushku, A., and Brown, S. 2007. Scale issues in developing a deforestation baseline for the region of 

the Noel Kempff Mercado Climate Action Project, Bolivia, pp 280-302 in G. LeClerc and C.A.S. Hall (Eds.) 

Making Development Work: Scientific Alternatives to Neoclassical Economic Theory, University of New Mexico 

Press, Albuquerque, NM.

Harris, N. L, S. Petrova, F. Stolle, and S. Brown. 2008. Identifying optimal areas for REDD intervention: East 

Kalimantan, Indonesia as a case study. Environmental Research Letters 3: 035006.

Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis. 2005. Very high resolution interpolated climate 

surfaces for global land areas. International Journal of Climatology 25: 1965-1978. http://www.worldclim.org/

Hong, B., Limburg, K.E., Hall, M.H. and Erickson, J.D. 2007. Scenario analysis of economy-ecology interactions 

in the Hudson River Basin, pp. 97-112 in Erickson, J.D., F. Messner, and I. Ring (Eds.), Advances in the 

Economics of Environmental Resources Vol 7: Ecological Economics of Sustainable Watershed Management. 

JAI Press, Elsevier, Amsterdam.
Integrity Applications Incorporated (IAI). 2010. High Resolution Land Use Change Analysis of Oil Palm in 

Sumatra and Kalimantan Circa 2010. BPA-09-03. September 20, 2010.

Integrity Applications Incorporated (IAI). 2011. High Resolution Land Use Change Analysis for Sumatra and 

Kalimantan Circa 2000. BPA-09-03. April 8, 2011.

IUCN and UNEP. 2010. The World Database on Protected Areas (WDPA). UNEP-WCMC. Cambridge, UK.

http://www.wdpa.org/

Jarvis, A., H.I. Reuter, A. Nelson, E. Guevara, 2008, Hole-filled SRTM for the globe Version 4, available from the

CGIAR-CSI SRTM 90m Database: http://srtm.csi.cgiar.org.

Kim, Oh S. 2010. An Assessment of Deforestation Models for Reducing Emissions from Deforestation and 

Forest Degradation (REDD). Transactions in GIS 14, no. 5 (10): 631-654. doi:10.1111/j.1467-9671.2010.01227.x.

Manson, S. M. 2002. Integrated assessment and projection of land-use and land-cover change in the southern 

Yucatan peninsular region of Mexico. Diss. Clark University.

Monserud, R.A. and Leemans, R. 1992. Comparing global vegetation maps with the Kappa statistic. Ecological 

Modelling, 62, pp. 275-293.

Murphy, D.J.R., Hall, M.H., Hall, C.A.S., Heisler, G., Stehman, S., and Anselmi-Molina, C. 2010. The relation 

between land cover and the urban heat island in northeastern Puerto Rico. International Journal of Climate, 19 

Apr 2010 (Articles on line in advance of print).

Paramananthan, S. 2003.Chapter XIV. Land Selection for oil palm in: T.H. Fairhurst and R. Härdter (eds). Oil 

Palm: Management for Large and Sustainable Yields. IPNI and IPI, Singapore.

Petrova, S. 2011. Methodology for predicting oil palm expansion for Indonesia and Malaysia. Preliminary 

report prepared for EPA. 

Pontius Jr, R G, J Cornell and C Hall. 2001. Modeling the spatial pattern of land-use change with Geomod2: 

application and validation for Costa Rica. Agriculture, Ecosystems & Environment 85(1-3) p.191-203.

Pontius Jr, R G and L Schneider. 2001. Land-use change model validation by a ROC method for the Ipswich 

watershed, Massachusetts, USA. Agriculture, Ecosystems & Environment 85(1-3) p.239-248.

Pontius Jr, R. G, and S. H Petrova. 2010. Assessing a predictive model of land change using uncertain data. 

Environmental Modelling & Software 25, no. 3: 299–309.
Pontius Jr, R. G, and H. Chen. 2006. GEOMOD modeling. Clark University.

Pontius Jr, R. G, D. Huffaker, and K. Denman 2004. Useful techniques of validation for spatially explicit landchange models. Ecological Modelling 179, no. 4: 445–461.

Rashmi, M. and Lele, N. 2010. Spatial modeling and validation of forest cover change in Kanakapura region 

using GEOMOD. Journal of the Indian Society of Remote Sensing p.45-54. 

Silva, E. and K. Clarke. 2002. Calibration of the SLEUTH urban growth model for Lisbon and Porto, Portugal. 

Computers, Environment and Urban Systems 26. pp. 525-552.

Smits, P.C., Dellepiane, S.G. and Schowengerdt, R.A. 1999. Quality assessment of image classification 

algorithms for land-cover mapping: a review and proposal for a cost-based approach. International Journal of 

Remote Sensing,20, pp. 1461-1486.

Wetlands International. 2008. Peatland Distribution in Sumatra and Kalimantan

http://www.wetlands.or.id/PDF/Atlas%20Review.pdf

Wilkinson, G.G. 2005. Results and implications of a study of fifteen years of satellite image classification 

experiments, IEEE Transactions on Geoscience and Remote Sensing, 43, pp. 433-440.

World Resources Institute, Indonesia Ministry of Forestry, Sekala and South Dakota State University. 2008. 

Interactive atlas of Indonesia's forests. World Resources Institute: Washington, DC.
For questions or comments contact:

Nancy Harris, PhD Carbon and Land Use Specialist, office +1.703.302.6519, nharris@winrock.org

Sean Grimland GIS Analyst, office +1.703.302.6504, sgrimland@winrock.org

Winrock International | 2121 Crystal Drive, Suite 500 | Arlington, VA 22202, USA | www.winrock.org
47 SPATIAL MODELING OF OIL PALM EXPANSION IN INDONESIA, 2000-2022 DRAFT REPORT