### Report submitted to EPA

July 2011



# Spatial Modeling of Future Oil Palm Expansion in Indonesia , 200 0 to 202 2

Nancy Harris

, Sean Grimland and Tim othy Pearson2 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 M ODEL D ESCRIPTION

................................ ................................ ................................ ................................ ...................... 7

2.2 D ATA INPUTS

................................ ................................ ................................ ................................ ................................ ... 8

2.3 C REATION OF F ACTOR M APS

................................ ................................ ................................ ................................ ......... 11

2.3.1 Elevation

................................ ................................ ................................ ................................ ................................ 12
2.3.2 Precipitation

................................ ................................ ................................ ................................ .......................... 14
2.3.3 Temperature

................................ ................................ ................................ ................................ ........................... 15
2.3.4 Slop e

................................ ................................ ................................ ................................ ................................ ...... 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 C REATING S UITABILITY

- FOR

- C HANGE M APS

................................ ................................ ................................ ................ 28

2.4 E VALUATING S UITABILITY M APS USING THE ROC STATISTI C

................................ ................................ ....................... 30

2.5 E VALUATING M ODEL P ERFORMANCE USING THE K APPA STATISTIC

................................ ................................ .............. 31

3.0 RESULTS AND DISC USSION

................................ ................................ ................................ ................................ ..... 32

3.1 E VALUATION OF S INGLE

- F ACTOR S UITABILITY M APS FOR P REDICTING C HANGE

................................ ......................... 33

3.2 E VALUATION OF M ULTI

- F ACTOR S UITABILITY M APS FOR P REDICTING C HANGE

................................ .......................... 36

3.3 F INAL S UITABILITY M AP

................................ ................................ ................................ ................................ ............... 37

3.4 P ALM O IL E XPANSION

, 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

................................ ................................ ................................ ................................ ............................... 43FIGURES

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. 39TABLES

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. S oil 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 se ttlements (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 perf ormance 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 c ontrol 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 Indonesia n 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 relati ve to various

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

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

soil type, land cove r 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 individual ly and in combinat ion 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 200 0 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 n early 88 % of th is additional expans ion w ill occur on the island of Sumatra

. There are eight

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

( 29 %), Sumatera Utara ( 18 %) and Sumatera Selatan (16%). Under th e 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 vegetati on mosaic.

# 1 .0 BACKG ROUND

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 Securit y 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 methodologie s. The U.S. Renewable Fuel

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

calculation of biofuels, which are calc ulated 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, t his analysis supports EPA’s efforts to update the assumptions and dat a 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, 200 0 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 R FS 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 diese l 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 Indonesia n 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. Firs t, we created maps of factors that could

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

Second

, we compared the factor maps against a map of existing palm oil plantations in 200 0 t o construct a

serie s of suitabi lity 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 value s 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, e ach candidate suit ability map was overlain with a map of existing plantations in the year 2009

. This is

the validation stage. The Relative Operating Characteris tic (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 200 0 and 2009

. The best performing suitability map was then used to simulate palm oil expa nsion

between 2000 and 2009 . T he level of agreement between actual and predicted locations of palm oil in 200 9

was evaluated using Kappa indices

.

Finally, t he 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

2 022

.

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 extrapola te 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 change 1 

### ; this

scenario is then used t o 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 t ime 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 us ed 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 GE OMOD.

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

Schneider (2001).



### 2 .2 Data inputs

The location s of oil palm plantations in the years 200 0 and 2009 w ere 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 200 0 map was used for model calibration and the 2009 map was used for model validation.

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

strata bins are created, in order to meet known criteria. The c riteria for this analysis was to find the best

### model run using a combination of one or more factors which produce d the best possible ROC s core and Kappa

for Location statistic. Model v alidation is the process of evaluat ing 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 200 0 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 dat a, close to

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

representing areas of new palm oil expansion between 200 0 and 2009

. However, approximately 1.5 million ha

of palm oil areas in 200 0 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

classificati on 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 represen ted by the actual net increase recorded between 2000 and

2009 ( Table 1 )

.

Table 1

. Total area (ha) of Sumatra and Kalimantan, Indonesia palm o il plantations

in 200 0 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 d ynamics of palm oil plantations in Indonesia between

200 0 and 2009

. Green areas represent areas that were classified as palm oil in both 200 0 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.

 P lantations could n ot expand into locations where palm oil already exists, so the existing palm oil

plantations in 200 0 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 oc cupied with immature

oil palm plantations.

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

locations.

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

(WDPA 20 10) 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

S everal different factors c an influence the location of future land use change, many of which vary spatially

over the landscape. The f actors 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, dist ance 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 val ues of 0

- 100 m, 100

- 200 m, etc.Table 2

. Factor data considered in the GEOMOD analysis

.

Biogeophysical Attribute Note / Source

Ele vation 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 Digitize d by Winrock from existing maps created by Wilmar International

Ltd. 2010.

Peatland Wetlands International ( 2008 )

Land Allocation W RI (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 Indonesia n islands of Sumatra and Kalimantan and across the areas covered by palm oil plantations. The

elevation factor map was created by st ratifying the range of elevation values into forty 100 m eter 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

A nnual precipitation data originated from the WorldClim global climate dataset 2



### . 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 a rea in either of the two Indonesia n islands receives less than 1, 54 0 mm or more than

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

.

Table 4

. R ainfall 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., 2005Figure 5 . Stratified precipitation factor map for Indonesia

.



### 2 .3.3 Temperature

A verage annual temperature data originated from the WorldClim global climate dataset 3



### . 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., 2005Unsuitable 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 slo pe 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 data set 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 . S tratified 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 Database 4 

### and 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%, Histosol s, 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

- Wor ld

- soil

- database/HTML

, FAO/IIASA/ISRIC/ISSCAS/JRC, 2009

### contain toxic levels of alumi num 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 a n organic soil with an organic carbon content by weight of 12

percent or more 5



### . Hi stosols 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 soi l 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

. S oil 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 C over 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 data set 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 y ear 2001 was chosen to represen t 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 oil product

.

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

q1Figure 9 . Stratified 2001 land cover factor map for the Indonesia n islands of

Sumatra and Kalimantan

.



### 2 .3. 7 Infrastructure ( Roads , Rivers, Railroads, Settlements)

R oads , rivers, railroads and settlements data originated from Global Map version 1 7 

### . The dataset s w ere

created and are maintained by t he Internation al Steering Committee for Global Mapping (ISCGM)

. Each data

layer was stratified into Euclidean distance s where each class width is 1000 meters (e.g. category 1= 0

- 1000m;

category 2 = 1000

- 2000m, etc) . This result ed in “ distance from ” factor s for the model ( i.e. distance from

roads, distance from rivers, etc.) ( Figure 10 ), with low values close to the feature and h igh 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 rail roads (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 200 0 oil palm plantations data originated from the 200 0 data on oil palm plantations.

First, we examined the distances between 200 0 plantations and new 2009 plantations and found that 6 0 % of

new plantations were located within a kilometer of a 200 0 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 distance s 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 t he

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 plantati ons.

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. wh ich 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. T he 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

d istance.

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 d istance 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 zo nes 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 h ave 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 diff erent 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 in to 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 commonl y in the existing plantations. The res ult is a suitability map where

reds (high values) are t he most suitable areas an d blues (low values) are the least suitable ( Figure 15 )

.Figure 15 . Example s uitability map created from one or several factor map s, 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,

W

a 

### = the weight of factor map a, and

P

a 

### (i) = percent

- developed in category a k 

### of factor map a, where cell(i) is a member of category a k



### .

For this analysis, 34 different suitability maps were crea ted. 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 200 0 and 2009 . After single factors were tested, we used this in formation 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 validati on 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 200 0 to 2022.

2 .4 Evaluating Suitabili ty Maps using the ROC statistic

The validity of e ach 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 qu antitative 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 ma y 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 suit ability 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 n ext largest suitability

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

### . A

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

change vs. n on

- 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 quan tity 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 a ll 10 groups. For

each scenario , the percent of grid cells that each scenario classifies correctly as non

- change and change is

calculat ed. 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 purp le 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 thi s 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 suitab ility maps for simulating land
use change in Indonesia

.


### What constitutes a “good” ROC? Any ROC above 50% (red line) is better than random. In othe r 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 P erformance 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 predic tion map that results when the suitability map is used to predict real

land cover change.

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

especially for purposes of accuracy asses sment

, because i t is g enerally 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 (1 892). 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 a l. (1999), and Wilkinson

(2005)

.

For all Kappa statistics, 0 indicates th at 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 ma p 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 quanti ty 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 hav e perfect location given no change in its quantities of each category.

F or 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 200 0 and 2009

. It

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

constraint. The suitability map for the model had an ROC score of 0.8 8

, a better score than any other single

factor suitability map and tied with one other combination

- factor suitability map. The Kno score of 0. 90 , and

Kappa for location scores of 0. 46

, were the best Kappa score s of all models runs including single

- factor a nd

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 tw o 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, n early 88 % of the additional expansion that occurred as a result of the palm

- only scenario

occurr ed 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

, offer ed 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 we re nearly identical in both their ROC and

Kappa scores, which was not surprising as slope is associated with elevation an d 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 wer e 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 a n elevation suitability map and shows

clearly t hat 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 o f 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 loc ations 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.

F actor Kno Klocation Agreement Due to

Location Disagreement Due

to Location

Elevation 0.896 0.424 3.848 5.224F actor 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.697Slope 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 score s

.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 Suitabilit y 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 inc luded

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 dist ance 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 me thods to group the distance

values into bins. The methods were 500 meter equal interval s instead of 1000

, q uantile 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 1 4

- 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 achieve d the highest combination of, ROC score (0.8 79 )

, Kappa for no information

(Kno ) ( 0. 901 ), Kappa f or 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 1 4 factors equally weighted.

Figure 21

. Predicted expansion of palm oil from 200 0

- 2022 in the e qual weight

scenario

. Green pixels show predicted expansion in both control and palm

- only scenarios, red pixels show additional expansion in the pal

m

- only scenario.

### 3 .4 Palm Oil Expansion , 200 0

- 2022

3 .4.1 Expansion by State

The area of predicted palm oil expansion between 2000 and 2022 per Indonesia n 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 t he three provinces of Riau

(2 9 %), Selatan ( 18 %) and Utara (16%) ( Figure 22 

### ).

Table 21

. Predicted area of palm oil expansion per Indonesian province between

2 00 0 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 200 0 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 so il

taxonomy, a histosol is a n 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 emis sions. 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. T he majority of the expansion occurs on acri sols that are found commonly throughout

Indonesia ( Table 22 

### ).Table 22

. Predicted palm oil expansion per soil type between 200 0 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 200 0 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 exa mple . When

l ocations of predicted palm oil expansion between 200 0 and 2022 were overlain with MODIS land cover type

for the year 2001

, most predicted expansion of palm oil between 200 0 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 CONCLUSION S

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 200 0 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 fac tor 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 200 0 and 2009

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

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

Kappa for location of 0. 4 6 , indicating good model performance. After this model validation step, t he model

was run to the year 2022 to predict locations of future palm oil expansion under a control case and a palm

only case. Additional expansion due to the palm

- only scenario occurred primarily on Sumatra, in the prov inces

of Riau

, Utara and Selatan , and 8 % of the ad ditional expansion occurred on histosolic soils.

The largest source of un certainty in the analysis was the quality of the input d ata for palm oil plantation

locations. A comparison of the 200 0 vs. 2009 s patial 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 we re not a lways 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 S ensing, 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 II ASA,

Laxenburg, Austria. http://www.iiasa.ac.at/Research/LUC/External

- World

- soil

- database/HTML/

G alton , F. 1892

. Finger Prints, pp. 216 (London: Macmillan)

Geri, F

. , Amici

, Valerio

. , and Rocchini , D .2011. Spatially

- based accuracy assessment of forestati on 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 Clim ate 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. Bro wn. 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 cli mate

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). 20 10. 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 Lan d Use Change Analysis for Sumat ra 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.2 010.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. Intern ational 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 p repared for EPA.

Pontius Jr, R G, J Cornell and C Hall. 2001. Modeling the spatial pattern of land

- use change with Geomod2:

appli cation 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 land

change models. Ecological Modelling 179, no. 4: 445

– 461.

Rashmi, M. and Lele, N. 2010. Spatial modeling and validation of forest c over 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 2 6. 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

experim ents, 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.org47 SPATIAL MODELING OF OIL PALM EXPANSION IN INDONESIA, 2000

- 2022 DRAFT REPORT