### R eport submitted to EPA

July 2011



# Spatial Modeling of Future Oil Palm Expansion in Malaysia , 200 3 to 202 2

Nancy Harris and Sean Grimland2 SPATIAL MODELING OF OIL PALM EXPANSION IN MALAYSIA, 2003

- 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

................................ ................................ ................................ ................................ ................................ ... 9

2.3 C REATION OF F ACTOR M APS

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

2.3.1 Elevation

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

................................ ................................ ................................ ................................ .......................... 13 2.3.3 Temperature

................................ ................................ ................................ ................................ ........................... 14 2.3.4 Slope

................................ ................................ ................................ ................................ ................................ ...... 15 2.3.5 Soil Type

................................ ................................ ................................ ................................ ................................ 16

2.3.6 Landcover in 2001

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

................................ ................................ ......................... 20 2.3.8 Distance from Existing Plantations

................................ ................................ ................................ ....................... 24

2.4 C REATING S UITABILITY

- FOR

- C HANGE M APS

................................ ................................ ................................ ................ 26

2.5 E VALUATING S UITABILITY M APS USING THE ROC STATISTIC

................................ ................................ ....................... 27

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

................................ ................................ .............. 29

3.0 RESULTS AND DISC USSION

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

3.1 E VALUATION OF S INGLE

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

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

3.2 E VALUAT ION OF M ULTI

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

................................ .......................... 34

3.3 F INAL S UITABILITY M AP

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

3.4 P ALM O IL E XPANSION

, 2003

- 2022

................................ ................................ ................................ ................................ 38 3.4.1 Expansion By State

................................ ................................ ................................ ................................ ................ 38 3.4.2 Expansion by Soil Type

................................ ................................ ................................ ................................ .......... 39 3.4.3 Exp ansion by Land Cover Type

................................ ................................ ................................ ............................. 40

4.0 CONCLUSIONS

................................ ................................ ................................ ................................ ............................. 42

5.0 REFERENCES

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

Figure 1. Components of the land cover change model GEOMOD. Factor maps are combined with a land cover map from one point in time to generate a series of suitability maps, which are then used to generate land cover change scenarios. Each scenario is evaluated against a land cover map from a second point in time and validation statistics are computed. From Pontius and Schneider (2001). 8
Figure 2. Representative dynamics of palm oil plantations in Malaysia between 2003 and 2009. Green areas represent areas that were classified as palm oil in both 2003 and 2009 (persistence). Blue areas represent areas of palm oil expansion between 2003 and 2009. Yellow areas repr esent areas that were classified as palm

oil in 2003 and as non

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

n on

- oil palm plantation areas (brown). Constraints to the model included existing oil palm plantations in 2003 (green) and inland water bodies and protected areas (black). 11
Figure 4. Stratifie d factor map for elevation 13
Figure 5. Stratified precipitation factor map for Malaysia. 14
Figure 6. Stratified temperature fact or map for Malaysia. 15
Figure 7. Stratified slope factor map for Malaysia, with slope categorized into 5% bins. 16
Figure 8. Stra tified soil factor map for Malaysia. 18
Figure 9. Stratified landcover factor map for Malaysia. 20
Figure 10. Stratified factor ma p for roads. 22
Figure 11. Stratified factor map for rivers. 23
Figure 12. Stratified factor map for railroads. 23
Figure 13. Stratified factor map for settlements. 24
Figure 14. Stratified distance from plantations factor map for Malaysia. 25
Figure 15. Example suitability map created from one or several factor maps, where red is most suitable and blue

is least suitable. 26
Figure 16. E xample ROC curve used to validate suitability maps for simulating land use change in Malaysia. 28
Figure 17. Ranking of single

- factor suitability map ROC scores. 32
Figure 18. Ranking of Kappa for Location scores. 33
Figure 19. Locations of 2009 palm oil plantations located in areas of low suitability with respect to the elevation factor map. Oil palm areas located in low suitability regions (red) may cause the ROC score to be less than 0.50 for the first 10% of modelled change in the landscape. 34
Figure 20. S uitability map that incorporates all 11 factors weighted equally. 36
Figure 21. Predicted expansion of palm oil from 2003

- 2022 in the (A) Equal weight scenario and (B) ROC weight scenario. Gree n pixels show predicted expansion in both control and palm

- only scenarios, pink pixels show additional expansion in the palm

- only scenario. The locations of additional expansion are slightly different depending on which suitability map is used for the simu lation. 37
Figure 22. Distribution of palm oil expansion in the control and palm only scenarios according to the prediction

model. 39
Figure 23. Malaysian states. 39TABLES

Table 1. Total area (ha) of Malaysian palm oil plantations in 2003 and 2009 according to two data sources. 9
Table 2. Factor data considered in the GEOMOD analysis. 11
Table 3. Elevation ranges of palm oil plantations (first ten bins) 12
Table 4. Rainfall suitability classes for growing oil palm. From Paramananthan (2003). 13
Table 5. Temperature suitability classes for growing o il palm. From Paramananthan (2003). 15
Table 6. Slope suitability classes for growing oil palm (first ten bins). 16

Table 7. Soil types in Malaysia and the reclassification scheme used to create the soil type factor map. 17
Table 8. Comparison of soil types in Malaysia vs. in palm oil plantations. 18
Table 9. Areas (ha) of each land cover type according to MODIS 2001 land cover data within the oil palm plantations and across Malaysia as a whole. 19
Table 10. Distance from roads (first ten bins). 20
Table 11. Distance from rivers (first ten bins. 21

Table 12. Distance from rail roads (first ten bins). 21
Table 13. Distance from settlements (first ten bins). 21
Table 14. Distance from 2003 oil palm plantati ons (first ten bins). 24
Table 15. Comparison of Kappa indices for model runs using single

- factor suitability maps. As model performance increases, Klocation increases and disagreement for loca tion decreases. 33
Table 16. Weighting schemes for combined

- factor GEOMOD runs where all factors were included in the model.

35

Ta ble 17. Ranking of ROC and Kappa scores for multi

- factor suitability maps. 35
Table 18. Predicted area of palm oil expansion per Malaysian State between 2003 and 2022 under a baseline (control case) and the additional expansion resulting from a palm

- only policy scenario. Results from two models are shown: the ‘equal weights’ model uses all factors weighted equally and the ‘varied weights’ model uses all factors weighted by their ROC scores. 38
Table 19. Predicted palm oil expansion per soil type between 2003 and 2022. Results displayed as for Table 18

above. 40

Table 20. P redicted palm oil expansion per land cover type between 2003 and 2022. Results displayed as for

Table 18 above. 41

# SUMMARY

The U nited S tates Environmental Protection Agency (EPA) has commissioned Winrock International to conduct

an analysis of land use change in Malaysia 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 pa lm oil plantations relative to various biogeophysical factors to predict the most likely

locations for future expansion. For this analysis, eleven factors were evaluated: elevation, precipitation,

temperature, slope, soil type, land cover type in 2001, dis tance to roads, distance to rivers, distance to

railroads, distance to settlements, and distance to existing plantations . The factors were used both individually

and in combinat ion to calibrate the model, and the results of each candidate model run were va lidated against

a 2009 map of oil palm plantations.

The most successful model s in predicting expansion used all eleven biogeophysical factors weighted equally

. The final model was run in two variations, first with all factors weighted equally and the s econd with different

weights assigned to different factors. The equal weights scenario scored slightly better on statistics that

measure how well the model predictions agree with the spatial data set for 2009. The model with all factors

weighted equally w as then used to project the quantity of palm oil expansion predicted to occur between

2003 and 2022 under two policy scenarios provided by EPA , a control case and a palm

- only case, and the

locations were compared with respect to the area of change in each Malaysian state as well as the area of

change by land cover type and by presence of peat soil , both of which affect the magnitude of expected GHG

emissions.

Based on the projections provided by EPA, a pproximately 5.2 million hectares of palm oil were exp ected to

occur in Malaysia in 2022, with 40,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 60% of th is additional expansion w ill

occur in the three states of Joh or (23%), Pahang ( 20%) and Sabah (16%) . Under th e palm

- only scenario, palm

oil will expand into histosol soils by 3,500 more hectares than under the control case , representing
approximately 9% of the additional expansion . The vast majority (> 90%) of additi onal 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 BACKG ROUND

Winrock International has been working with EPA to estimate the lifecycle greenhouse gas (GH G) 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 rece nt 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. T he RFS2 final rule, published

in March 2010 (EPA 2010), includes factors for indirect land use change (iLUC) as part of the lifecycle GHG

calculation of biofuels, which are calc ulated internationally based on the results of the FAPRI

- CARD partial

equilibri um model (EPA 2010 a )

.

Because palm oil is an important biofuel feedstock, EPA commissioned Winrock to conduct a more detailed

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

produced

. Therefore, t his analysis supports EPA’s efforts to use the best available assumptions and data in

economic and spatial models for biofuel lifecycle greenhouse gas (GHG) analysis related to iLUC.

EPA obtained spatial data from the Malaysian Government on the locatio ns of existing palm oil plantations for

two years, 2003 and 2009. 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 RFS2

f uels are produced at their mandated levels in 2022 . This is then compared to a “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

.

Given these data inputs, it became possible to use 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 pa lm oil plantations in Malaysia in 2003 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 futur e 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 2003 t o construct a series 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 values to locations that are

less similar to existing palm oil areas. Third, e a ch candidate suitability 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 sui tability map were also the

### locations of actual change to palm oil between 2003 and 2009. The best performing suitability map was then

used to simulate palm oil expa nsion between 2003 and 2009. T he level of agreement between actual and

predicted locations o f palm oil in 2009 was evaluated using Kappa indices

. Finally, t he suitability map that

performed best in predicting the location of change between 2003 and 2009 was chosen to simulate locations

of future palm oil expansion across Malaysia 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 poin ts 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; this

scenario is then used to predict land cover change either forw ards 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 ( Po ntius 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 Ba sin, USA (Hong et al. 2007), and Bolivia (Hall et al.

2007)

. Pontius et al. (2001) contains the most complete peer

- reviewed description and application of

GEOMOD.Figure 1 . Components of the land cover change model GEOMOD. Fac tor maps are

combined with a land cover map from one point in time to generate a series of suitability maps, which are then used to generate land cover change scenarios. Each scenario is evaluated against a land cover map from a second point in time

and va lidation statistics are computed. From Pontius and Schneider (2001).

### 2 .2 Data inputs

The location s of oil palm plantations in the years 2003 and 2009 w ere supplied in shapefile format to EPA by

the G overnment of Malaysia. The shapefi l e s contained no metadata and the accuracy of the maps is

unknown. The shapefiles included an age attribute (immature, mature or old ) tha t was not accounted for in

our analysis 1 

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

Table 1 shows the total area of palm oil in 2003 and 2009 according to the spatial data supplied by the

Malaysian government. Of the 5.2 million hectares of palm oil present in 2009 according to the spatial data,

close to 3 millio n hectares were represented also in the older 2003 spatial data, with the remaining 2.2 million

ha representing areas of new palm oil expansion between 2003 and 2009 . However, approximately 8 66 ,000 ha

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

landscape. Therefore, the spatial data indicate a net increase in oil palm area between 2003 and 2009 of

1,386,025 ha

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

replanted between 2003 and 2009 may have been immature in 2003. Other explanations may include land

abandonment and classification errors in the data. Assuming that no areas were retired in this time period

gives a total area of palm oil in 2009 of 6,084,280 ha an increase of 59 % from 2003 . 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 (3 7 %) Table 1

.

Table 1

. Total area (ha) of Malaysian palm oil plantations in 2003 and 2009 according to two data sources.

Source 200 3 ha 2009 ha 200 3

- 2009

Increase

Spatial data (IAI 2011) 3, 831,270 5,217,295 36 %

1 Both 2003 and 2009 sh apefiles had approximately the same proportion of mature to immature areas (3:1).Figure 2

. Representative d ynamics of palm oil plantations in Malaysia between 2003

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

betwe en 2003 and 2009. Yellow areas represent areas that were classified as palm

oil in 2003 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 la ntations could not expand into locations where palm oil already exists, so the existing palm oil

plantations in 2003 were excluded from the model’s spatial domain.

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

locations.

 A reas protected by the Malaysian 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 (brown). Constraints to

the model include d existing oil palm plantation s in 2003 (green) and inland water bodies and protected areas (black).



### 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 GEO MOD analysis were chosen based on both data availability as well as past experience; certain factors (e.g., elevation) have been shown in other analyses (e.g., Harris et al.

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

.

For this analysis, e leven 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, and

distance to existing plantations

. Detai ls 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’ s o that categorical values of suitability can be created . For
example, elevation data was stratified into 100m ‘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

Elevation SRTM 90m digital elevation data v 4
Jarvis, A., H.I. Reuter, A. Nelson, E. Guevara (2008)

Precipitation WorldClim database v 1.4Biogeophysical Attribute Note / Source

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 cov er 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)

R oads Global Map v1.0
Source: http://www.iscgm.org/gmd/

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

S ource: http://www.iscgm.org/gmd/

Existing Palm Oil Plantations 2003 Data delivered to Winrock by EPA

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

Malaysia and across palm oil plantations. The elevation factor map was created by stratifying the range of

elevation values into 40 100

- m 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 Malaysia

2003 2009

1 0

– 100 88.08% 86.11% 46.08%2 101

- 200 8.74% 9.95% 13.08%

3 201

- 300 1.70% 1.93% 8.35%

4 301

- 400 0.81% 0.73% 6.58%

5 401

- 500 0.39% 0.43% 5.47%

6 501

- 600 0.16% 0.39% 4.33%

7 601

- 700 0.09% 0.17% 3.44%

8 701

- 800 0.02% 0.16% 2.81%

9 901

- 1000 0.01% 0.12% 2.24%

10 1001

- 1 100 0.01% 0.01% 1.84%

Figure 4

. Stratified factor map for elevation



### 2 .3.2 Precipitation

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



### . Error! Reference source not

ound. shows the distribution of land area by annual precipitation across Malaysia and across palm oil

plantations. The precipitation map was stratified into five suitability classes (0

- 4) based on Paramananthan

(2003) ( Table 4 ) , resulting in the stratified precipitation map shown in Figure 5

. Note that Figure 5 shows only

three suitability classes . A ccording to the WorldClim data, no land area in Malaysia receives less than 1, 450

mm or more than 5,000 mm of precipitation per year.

Table 4

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

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

2 http://www.worldclim.org/

, Hijmans, R.J., 2005Climatic

Suitability Gro wth

Limitations Mean Annual

Precipitation (mm) Precipitatio

n Bin Palm Oil Malaysia

200 3 2009

Unsuitable Very Severe < 1,250
> 6,000 0 0% 0% 0%

Suitable 1 None 2, 501

- 3, 500 1 44.92% 45.29% 43.73%

Suitable 1 Slight 1,701

- 2, 500
3, 501

- 4,000 2 53 .92% 52.30% 47.39%

Moderate Moderate 1,451

- 1,700
4,001

- 5,000 3 0.97% 2.09% 8.87%

Marginal Severe 1,250

- 1,450 5,001

- 6,000 4 0% 0% 0%

Figure 5 . Stratified precipitation factor map for Malaysia.



### 2 .3.3 Temperature

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



### . Table 5

. Temperature suitability classes for growing oil palm . From Paramananthan (2003). Error! Reference source not

found. shows the distribution of land area by temperature across Malaysia and across palm oil plantations by

strata . The temperature map was stratified into five suitability classes (1

- 5) based on Paramananthan (2003),

resulting in the stratified temperature ma p shown in Figure 6

.

3 http://www.worldclim.org/

, Hijmans, R.J., 2005Table 5 . Temperature suitability classes for growing oil palm . From Paramananthan

(2003).

Area in each temperature suitability class as a proportion of

area in:

Climatic

Suitabil ity Growth

Limitations Mean Annual

Temperature (C) Temperature

Bin Palm Oil Malaysia

200 3 2009

Suitable 1 None 25

– 29 1 99.16% 98.35% 72.77%

Suitable 2 Slight 22

- 25

29

– 32 2 0.65% 1.32% 20.34%

Moderate Moderate 20

- 22

32

– 35 3 0.01% 0.01% 4. 82%

Marginal Severe 10

- 20

35

– 37 4 0% 0% 2.07%

Unsuitable Very Severe < 10

> 37 5 0% 0% 0.01%

Figure 6 . Stratified temperature factor map for Malaysia.



### 2 .3.4 Slope

Slope is the steepness, incline, or grade of a surface. T he 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 i n 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

. The stratified slope map for Malaysia 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 Malaysia

2003 2009

1 0

- 5 89.35% 87.52% 51.30%

2 5

- 10 7.01% 8.46% 1 4.76%

3 10

- 15 2.08% 2.45% 10.76%

4 15

- 20 0.87% 0.95% 8.25%

5 20

- 25 0.38% 0.39% 5.98%

6 25

- 30 0.18% 0.14% 3.83%

7 30

– 35 0.08% 0.05% 2.12%

8 35

– 40 0.03% 0.02% 1.02%

9 40

– 45 0.01% 0.% 0.44%

10 45

- 50 0.% 0.% 0.19%

Figure 7 . S tratified slope factor map for Malaysia , 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 F amilies ( 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 acro ss 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 in fertile 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 a n organic soil with an organic car bon content by weight of 12

percent or more 5



### . 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 signifi cant GHG emissions. For this study, histosol soils will be used as a proxy for peat

soils because spatial data for peat soils in Malaysia was not available.

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

to develop a characteristic gleyic colour pattern.

The stratified soil type factor map for Malaysia is shown in Figure 8

.

Table 7 . Soil types in Malaysia and the reclassification scheme used to cr eate the soil type factor map.

S oil Name R eclass Name R eclass 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 Fl uvisols 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

Albic Arenosols Arenosols 12

5 http://soils.usda.gov/technical/classification/taxonomy/Table 8 . Comparison of soil types in Malaysia vs. in palm oil plantations.

Area in each soil type as a
percentage of area in:

R eclass Name Palm Oil M alaysia

2003 2009

Acrisols 70.11% 67.27% 79.77%

Ar enosols 0.46% 0.51% 0.35%

Cambisols 0% 0.03% 0.01%

Ferralsols 1.25% 1.02% 0.36

Fluvisols 6.01% 5.84% 4.01%
Gleysols 11.08% 9.09% 4.16%

Histosols 8.38% 13.32% 8.04%

Lithosols 0.14% 0.36% 1.27%

Luvisols 2.43% 2.35% 1.71%

Nitosols 0.12% 0.09% 0.05%

P odzols 0.02% 0.01% 0.01%
Regosols 0.01% 0.10% 0.26%

Figure 8 . Stratified soil factor map for Malaysia.



### 2 .3. 6 Landcover in 2001

Landcover 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

6 

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

q1

### (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 landcover. The primary land cover scheme includes 11 natural

vegetation classes, 3 developed and mosaicked land classes, and three non

- vegetated l and classes ( Table 9 ).

The year 2001 was chosen to represen t the land cover condition of areas before conversion to oil palm

between 2003 and 2009. Approximately 92 % of palm plantations in 2003 and 2009 were locate d in areas that

were classified as broadleaf evergreen forest or natural mosaic classes in 2001. Data were stratified ( Figure 9 

### )

by assigning broadleaf evergreen forest and natural mosaic classes a value of 1 and all other catego ries a value

of 0.

Table 9

. Areas (ha) of each land cover type according to MODIS 2001 land cover

data within the oil palm plantations and across Malaysia as a whole.

Area in each land cover type as a percentage of area in:

MOD IS 2001 Palm Oil Malaysia

2003 2009
barren or sparsely vegetated 0.01% 0.01% 0.02%

closed shrublands 0.02% 0.02% 0.03% cropland/natural vegetation mosaic 31.61% 30.53% 13.93% croplands 3.30% 2.96% 1.88%

deciduous broadleaf forest 0.24% 0.27% 0.25%

deciduous needleleaf forest 0.00% 0.00% 0.00% evergreen broadleaf forest 60.51% 61.45% 78.38% evergreen needleleaf forest 0.00% 0.00% 0.05% grasslands 0.71% 0.72% 0.30%

mixed forests 0.59% 0.64% 0.55% open shrublands 0.11% 0.09% 0.05%

permanent snow and ice 0.00% 0.00% 0.00% permanent wetlands 1.01% 1.52% 2.82%

savannas 1.10% 1.03% 0.30%

urban and built

- up 0.23% 0.23% 0.85%

water 0.03% 0.05% 0.23% woody savannas 0.52% 0.49% 0.34%Figure 9 . Stratified landcover factor map for Malaysia.



### 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 International Steering Committee for Gl obal 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.), with low values close to the feature and h igh values far from the feature.

Table 10

. Distance from roads (first ten bins).

Area in each ‘ distance from roads ’ bin as a percentage of area in:

‘Dista nce from ‘

Bin ‘Distance

from’ Range (km) Palm Oil

Malaysia

2003 2009

1 0

- 1 24.80% 21.94% 14.44%

2 1

- 2 19.84% 17.66% 10.72%

3 2

- 3 13.58% 12.23% 7.57%

4 3

- 4 9.79% 9.03% 5.97%

5 4

- 5 6.76% 6.55% 4.70%

6 5

- 6 4.72% 4.95% 3.85%

7 6

- 7 4.09% 4.48% 3.77%

8 7

- 8 2.86% 3.40% 3.05%

9 8

- 9 2.12% 2.69% 2.59%

7 http://www.iscgm.org/gmd/10 9

- 10 1.55% 2.05% 2.19%

Table 11

. Distance from rivers (first ten bins )

.

Area in each ‘ distance from rivers ’ bin as a percentage of area in:

‘Dist ance from ‘

Bin ‘Distance

from’ Range (km) Palm Oil Malaysia

2003 2009

1 0

- 1 28.10% 27.24% 29.32%

2 1

- 2 23.98% 24.49% 25.45%

3 2

- 3 15.05% 16.10% 16.75%

4 3

- 4 9.52% 10.39% 10.47%

5 4

- 5 5.64% 6.09% 6.04%

6 5

- 6 3.71% 3.77% 3.64%

7 6

- 7 2.83% 2.75% 2.51%

8 7

- 8 2.12% 1.91% 1.57%

9 8

- 9 1.62% 1.42% 1.02%

10 9

- 10 1.35% 1.12% 0.71%

Table 12

. Distance from rail roads (first ten bins).

Area in each ‘ distance from railroads ’ bin as a percentage of

area i n:

‘Distance

from ‘

Bin ‘Distance

from’ Range (km) Palm Oil Malaysia

2003 2009

1 0

- 1 1.49% 1.34% 1.05%

2 1

- 2 1.73% 1.54% 1.07%

3 2

- 3 1.56% 1.41% 0.99%

4 3

- 4 1.66% 1.36% 0.98%

5 4

- 5 1.63% 1.34% 0.93%

6 5

- 6 1.53% 1.28% 0.89%

7 6

- 7 1.73% 1.48% 1.01%

8 7

- 8 1.72% 1.47% 0.94%

9 8

- 9 1.66% 1.47% 0.91%

10 9

- 10 1.57% 1.38% 0.85%

Table 13

. Distance from settlements (first ten bins).Area in each ‘ distance from settlements ’ bin as a percentage of

are a in:

‘Distance

from ‘

Bin ‘Distance

from’ Range (km) Palm Oil Malaysia

2003 2009

1 0

- 1 1.52% 1.45% 1.90%

2 1

- 2 5.98% 5.82% 6.46%

3 2

- 3 9.22% 9.20% 8.70%

4 3

- 4 10.82% 10.84% 9.47%

5 4

- 5 10.22% 10.10% 8.60%

6 5

- 6 9.34% 9.39% 7. 93%

7 6

- 7 9.04% 8.99% 7.57%

8 7

- 8 7.61% 7.55% 6.32%

9 8

- 9 6.02% 6.11% 5.08%

10 9

- 10 4.91% 4.97% 4.11%

Figure 10 . Stratified factor map for roads

.Figure 11 . Stratified factor map for ri vers.

Figure 12

. Stratified factor map for railroads.Figure 13

. Stratified factor map for settlements.



### 2 .3. 8 Distance from Existing Plantations

The distance from 2003 oil palm plantations data origi nated from the 2003 shapefile of oil palm plantations.

First, we examined the distances between 2003 plantations and new 2009 plantations and found that 68% of

new plantations were located less than 1.5 km from a 2003 plantation, and more than 90% were l ocated less

than 2 km from a 2003 plantation.

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 facto r for t he

model ( Figure 14 ).

Table 14

. Distance from 200 3 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 Malaysia

2003 2009

1 0

- 1 100% 72.03% 19.15%

2 1

- 2 7.72% 6.92%

3 2

- 3 3.98% 4.74%

4 3

- 4 2.72% 3.84%

5 4

- 5 1.92% 3.18%6 5

- 6 1.47% 2.82%

7 6

- 7 1.33% 2.89%

8 7

- 8 0.97% 2.57%

9 8

- 9 0.75% 2.27%

10 9

- 10 0.57% 2.01%

Figure 14 . Stratified distance from plantations factor map for Malaysia.

### 2.4 Creating Suitability

- for

- Change Maps

After factor maps were created, they were used individually and in combin ation 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 largest suitability values. If the suitability map were perfect, the order of the suitability v alues

would match the order in which humans change the landscape, with the largest suitability values being
changed first (Pontius and Schneider 2001). For validation, a map of simulated future change is compared to a

map of recent real land

- cover change. Several model runs , each producing a unique suitability map

, 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 in to the same 0

- 100 m bin in locations of exist ing 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 ha ve an

elevation value that is found least commonl y in the existing plantations. The res ult is a suitability map where

red (high values) are t he most suitable areas and blue (low values) are the least suitable ( Figur e 15 )

.

Figur e 15 . Example s uitability map created from one or several factor map s, where

red is most suitable and blue is least suitable.



### A suitability map can also be created by combining multiple individual factor maps together and weightin g 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

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, 18 different suitability maps were created. In the first stage, we created eleven suitability

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

oil expansion between 2003 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 we re 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 expans ion of oil palm plantations

from 2003 to 2022.

2.5 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 o f 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 tha n 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 ch ange 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 t hat

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 reality. A

tabular result shows the propor tion 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; how ever, 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 a ll 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.

To define the ROC, the rate of true

- positives (change in both reality and model ) is plotted on t he vertical axis

and the rate of false

- positives (change in model, no change in reality) on the horizontal axis for each of the 11

scenarios ( Figure 16 ). The ROC statistic is calculated as the area under the blue curve that connec ts the plotted

points. 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 c hange increases from 0% to 100% (horizontal axis), there are no (0%) false

positives and all (100%) true positives (purple lines in Figure 16 ) . If the order of suitability values were assigned

at random locations across the landsc ape, 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 Malaysia

. 

### 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% f or material strength testing, and 55 to 98% for polygraph lie detection

( Pontius and Schneider 2001 ).

### 2.6 Evaluating Model Performance using the Kappa statistic

The ROC statistic assesses how well a candidate suitability map assigns locations of land co ver 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 when a single quantity of change is simulated. In other

words, the ROC scor e 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.

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

especially for purposes of accuracy asses sment

, because i t 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 C ohen (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 Gr een (1999), Smits et a l. (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 loc ation 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 ea ch category (in this case, palm oil vs. non

- palm oil)

.

Kappa indices can also be expressed as percentages, which allow for an easy review of how much each index

contributes to the results 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 rea lity

 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.

F or this analysis, the important kappa indices are the ones co ncerned 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 mode l performed at

predicting correctly the locations of land cover change.

# 3 .0 RESULTS AND DISCUSSION

Two models resulted in the best predictions of palm oil expansion in Malaysia between 2003 and 2009. Both

incorporated all eleven biogeophysical factor s , but one variation weighted all factors equally and the other

weighted factors by their individual ROC score s

. The ROC weighted scenario resulted in the highest ROC score
but the equal weights scenario generated better scores for all four Kappa indices ( Table 15 

### ). Based on the

Kappa scores we determined that the equal weights scenario produced a superior prediction map. For this

reason the equal weights scenario 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 model

. The model projected n early 60% of the

additional expansion that occurred as a resu lt of the palm

- only scenario occurr ed in the three states of Johor

(23%), Pahang ( 20%) and Sabah (16%)

. Under the palm

- only scenario, palm oil expanded into histosol soils by

3,500 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.

3 .1 Evaluation of Single

- Factor Suitability Maps for Predicting Change

The results of eleven model runs , each whic h 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 plantations (DEP), followed by elevation and slope ( Figure 17).

This result was to be expected

, as DEP is not independent to the initial map and therefore may violate the

model’s requirement for using independent factors to predict future land use change

. It is also not surprising

that areas adjacent to existing plantations would tend to be favourable for expansion. The factors of elevation

and slope we re nearly identical in both their ROC and Kappa scores, which was not surprising as slope is

derived from elevation and there was little variability in slope across existing oil palm areas.

The roads and railroads factor map s were the next best factors in predicting palm oil expansion. The

performance of the railroads factor map was somewhat surprising, as the data were limited to main land

Malaysia with a single rail line in Sabah. Temperature and soil type were also good explanatory v ariables for

oil palm expansion. Precipitation, land cover type, and distance from settlements and rivers wer e poor predictors of expansion when used in dividually, perhaps 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 18



### , Table 15 

### ).

For some factors with the highest ROC scores (e.g., elevation, slope), Figure 17 

### 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 in Sarawak exist far inland and at high elevations, and are outliers compared to

where most other existing plantations occur. Figure 19 

### shows 2009 plantations overlai n onto the elevation

suitability map and show c learly that several plantations exist in the lowest suitability areas for the elevation

### 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 exce eds 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 17 

### ). 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 category of interest (new oil palm) concentrated at the locations of relatively h igh 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?’Dist. From Plantations ROC = 0.799 Elevation ROC = 0.75 0 Slope ROC = 0.723

Dist. From Roads ROC = 0.677 Dist. From Railroads ROC = 0.657 Temperature ROC = 0.649

Soil Type ROC = 0.595 Dist. From Settl ements ROC = 0.57 0 Precipitation ROC = 0.54 0

Dist. From Rivers ROC = 0.517 Land C over Type 2001 ROC = 0.492

Figure 17

. Ranking of single

- factor suitability map ROC score s

.Figure 18 . Ranking of Kappa for Location scores.

Table 15 . Comparison of Kappa indices for model runs using single

- factor suitability maps

. As model performance increases, Klocation increases and disagreement for location d ecreases.

Kappa index E levation P recipitation Temperature Slope Soil Land Cover

Kno 0.7315 0.7194 0.7239 0.7292 0.7324 0.7178

Klocation 0.5311 0.5100 0.5178 0.5270 0.5326 0.5071 Agreement location 15.2040 14.6011 14.8244 15.0881 15.2470 14.5192 Disagr eement location 13.4249 14.0279 13.8046 13.5409 13.3820 14.1098

Kappa index Roads Rivers Railroads Settlements Distance to Existing Plantations

Kno 0.7265 0.7235 0.7362 0.7204 0.7625

Klocation 0.5223 0.5172 0.5392 0.5116 0.5852 Agreement location 14.9 541 14.8060 15.4380 14.6474 16.7538 Disagreement location 13.6749 13.8230 13.1909 13.9816 11.8752Figure 19 . Locations of 2009 palm oil plantations located in areas of low

suitability with respect to the elevation factor map. Oil palm areas located in low suitability regions (red) may cause the ROC score to be less than 0.50 for the first 10% of modelled change in the landscape.



### 3

. 2 Evaluation of Multi

- Factor Suitability Maps for Predicting Change

After running single

- facto r scenarios to evaluate which factors performed best in predicting plantation

expansion, we improve d model performance further by generating suitability maps that inc luded multiple

factors. First, we created a suitability map from the combination of all el even factors where each factor was

weighted equally

. Second , we created a suitability map from all eleven factors, but different weights were

assigned to different factors in proportion to individual ROC score s ( Table 16 

### ). Third, we removed from

consideration redundant factors (slope ; redundant with elevation ) and poorly performing factors ( distance

from settlements, precipitation, distance from rivers, and land cover type , see Figure 17 

### ), included the

rem aining six fact ors (distance from existing plantations, elevation, distance from roads, distance from

railroads, temperature, and soil type) and tested different weighting scenarios

. The 11

- factor suitability maps

showed a significant improvement over any single

- factor map, achieving ROC scores of 0.881 and 0.879 for

### ROC

- weighted and equal

- weighted scenarios, respectively. The other six

- factor scenarios also performed

better than single

- factor scenarios, but to a lesser extent ( Table 17 

### ).

Table 16 . Weighting schemes for combined

- factor GEOMOD runs where all factors were

included in the model

.

Factor Weight (Equal Weight Scenario) Weight ( ROC Weight ed Scenario)

Elevation 0.0909 0.1100

Precipitation 0.0909 0. 0773

Temperature 0.0909 0.0929

Slope 0.0909 0.1035

Soil Type 0.0909 0.0851

Land Cover 2001 0.0909 0.0704

Distance from Roads 0.0909 0.0969

Distance from Rivers 0.0909 0.0740

Distance from Railroads 0.0909 0.0940

Distance from Settlements 0.0909 0.0 816

Distance from Existing Plantations 0.0909 0.1143

Total 1 .0000 1 .0000

Table 17 . Ranking of ROC and Kappa scores for multi

- factor suitability maps.

Weighting Scenario ROC

Score Kno Klocation Agreement

due to

location Disagree ment

due to location

All Factors, Weighted by ROC score 0.881 0.7687 0.5960 17.0638 11.5651

All Factors, Equal Weights 0.879 0.7699 0.5981 17.1232 11.5058

6 Factors, Weighted by Standard

Deviation 0.874 0.7644 0.5920 16.9485 11.6805

6 Factors, Weighted by Pairwise Comparison v1 0.874 0.7633 0.5867 16.7956 11.8334

6 Factors, Weighted by ROC score 0.874 0.7663 0.5918 16.9421 11.6869

6 Factors, Equal Weights 0.873 0.7662 0.5917 16.9401 11.6889

6 Factors, Weighted by Pairwise Comparison v2 0.856 0.7664 0 .5920 16.9485 11.6805

### 3 .3 Final Suitability Map

The suitability map shown in Figure 20 

### corresponds to the scenario of all factors included and weighted

equally ; this scenario achieve d the second highest ROC score ( 0.879 ) a nd best scores for all four Kappa

metrics

. The equal weights scenario was determined to be superior based on the combination of second best

ROC score and best Kappa scores. T he suitability map for the next best scenario (all factors weighted

according to ROC scroe ) looks virtually identical to Figure 20 

### at the country scale. However, zooming in to

specific areas of the prediction map ( Figure 21 

### ) shows slight variations in where additional palm oil expansi on is predicted to occur in the palm

- only case depending on which suitability map is used in the simulation

.

Figure 20

. Suitability map that incorporates all 11 factors weighted equally

.Figure 21

. Predicted expansion of palm oil from 2003

- 2022 in the (A) Equal weight

scenario and (B) ROC weight scenario . Green pixels show predicted expansion in both control and palm

- only scenarios, pink pixels show additional expans ion in the palm

- o nly scenari o. The locations of additional expansion are slightly different depending on which suitability map is used for the simulation. (A)

( B )

### 3 .4 Palm Oil Expansion , 2003

- 2022

3 .4.1 Expansion By State

The area of predicted palm oil expansion between 2003 and 2022 per Ma lays ian State is shown in Table 18



### . Of

the 40,764 ha of additional expansion predicted in the palm

- only scenario , nearly 60% is expected to occur in
the three states of Johor (23%), Pahang (20%) and Sabah (16%) (Figure 22)

. Sabah is project ed to have the

greatest area of expansion in the control case but Johor has more additional expansion in the palm only case.

Figure 23 is a map of the Malaysian states.

Table 18

. Predicted area of palm oil expansion per Malaysian State between 2 003

and 2022 under a baseline ( control case ) and the additional expansion resulting from a palm

- only policy scenario. Results from t wo models are shown: the ‘ equal weights’ model uses all factors weighted equally and the ‘varied we ights’ model uses all factors weighted by their ROC scores.

EQUAL WEIGHTS ROC WEIGHTS

State Control Case

(ha) Additional

Expansion,

Palm

- Only Case (ha) Control Case

(ha) Additional

Expansion,

Palm

- Only Case (ha)

Johor 909,766 9,317 917,988 9,960

Kedah 123,472 1,159 124,545 880

Kelantan 98,958 386 99,044 493

Melaka 82,214 923 84,425 1,374

Negeri

Sembilan 221,270 3,070 222,665 2,490

Pahang 585,911 8,200 590,032 7,427

Perak 571,658 4,851 569,489 4,036

Perlis 773 43 751 22

Pulau Pinang 30,932 451 31, 769 516

Sabah 1,351,126 6,333 1,336,765 6,011

Sarawak 834,400 2,726 832,146 3,091

Selangor 259,394 1,352 257,483 1,696

Trengganu 131,650 1,954 134,419 2,769

Total 5,201,523 40,764 5,201,523 40,764Figure 22 . Distribution of palm oil expansion in the control and palm only

scenarios according to the prediction model.

Figure 23 . Malaysian states.



### 3 .4. 2 Expansion by Soil Type

The area of predicted palm oil expansion between 2003 and 2022 per soil typ e is important primarily for determining how much expansion occurs on histosols

. In this study, histosol soils are used as a proxy for peatsoils because a more detailed map of peat soil was not available. It is a logical assumption as according to the

USA soil taxonomy, peat soil is a type of histosol 8



### . As mentioned above (Section 2.3.5), 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

.



### When palm oil is cultivated on peat, the peat must be drained as part of site preparation activities and results

in significant GHG emis sions. Approximately 10% of the total expansion between 2003 and 2022 is projected

to occur on peat soils. Approximately 9 % (3,477 ha) of the additional expansion in the palm

- only scenario is

predicted to occur on peat soils, with the majority of the remaining expansion occurring on acrisols that are

found commonly throughout Malaysia ( Tabl e 19 

### ).

Tabl e 19

. Predicted palm oil expansion per soil type between 2003 and 2022

. Results displayed as for Table 1 8 above.

EQUAL WEIGHTS VARIED WEIGHTS

Soil Type Control Case (ha) Additional Expansion,

Palm

- Only Case (ha) Control Case

(ha) Additional

Expansion,

Palm

- Only Case (ha)

Acrisols 3,378,663 30,009 3,400,794 30,418

Arenosols 26,596 193 26,704 171

Ferralsols 45,443 666 44,606 258

Fluvisols 391,838 3,799 389,026 3,048
Gleysols 663,682 1,395 653,357 1,438

Histosols 525,678 3,477 521,835 4,272

Lithosols 11,463 0 11,484 43

Luvisols 140,859 1,052 136,566 987

Nitosols 7,213 0 7,105 43

Podzols 923 21 880 0
Regosols 2,361 0 2,361 0
No soil type information 6,805 150 6,805 86

Total 5,201,523 40,764 5,201,523 40,764



### 3 .4.3 Expans ion by Land Cover Type

The area of predicted palm oil expansion between 2003 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 existing annual cropland is simply replaced with palm oil.

In fact replacing grasslands or savannah with oil palm can actually result in net carbon sequestration over time

because of the carbon stored in oil palm trees and related biomass. When l ocations of predicted palm oil

8 http://soils.usda.gov/technical/cla ssification/taxonomy/

### expansion between 2003 and 2022 were overlain with MODIS land cover type s for the year 2001

, most

predicted expansion of palm oil ( between 2003 and 2022 ) occurred in the forest and mixed (cropland/natu ral

vegetation mosaic) land cover types ( Table 20 

### ). The land cover results are based on the assumption that land

cover in 2001 predicts the type of land cover that will be present at the time of conversion , between 2003 and

2022.

Table 20 . Predicted palm oil expansion per land cover type between 2003 and 2022. Results displayed as for Table 1 8 above.

EQUAL WEIGHTS ROC WEIGHTS

Land Cover Type Control Case (ha) Additional

Expansion,

Palm

- Only Case (ha) Co ntrol Case

(ha) Additional

Expansion,

Palm

- Only Case (ha) barren or sparsely vegetated 665 0 623 64

closed shrublands 880 21 1,073 0 cropland/natural vegetation

mosaic 1,847,997 14,404 1,837,650 12,278 croplands 140,558 2,232 159,878 3,520

deciduous b roadleaf forest 11,398 21 11,849 21

deciduous needleleaf forest 86 0 86 0 evergreen broadleaf forest 2,968,429 22,045 2,952,265 22,174 evergreen needleleaf forest 236 0 258 0 grasslands 32,371 150 33,272 279

mixed forests 29,838 129 30,696 193

open s hrublands 4,121 43 4,594 86 permanent snow and ice 86 0 86 0 permanent wetlands 73,070 966 70,623 558

savannas 44,585 193 47,096 601

urban and built

- up 18,053 258 20,929 580

water 4,443 86 4,293 64 woody savannas 24,707 215 26,253 343

TOTAL 5,201,52 3 40,764 5,201,523 40,764

# 4.0 CONCLUSION S

The GEOMOD model was used to predict locations of future expansion of palm oil across Malaysia until the

year 2022. Eleven spatial factors were tested for their ability to predict the distribution of plantation expansion

between 2003 and 2009 and were combined to create a suitability map. The two best suitability map s had ROC

score s of 0.881 and 0.879 , indicating that most new palm oil plantations in 2009 were concentrated in

locations that were assigned high su itability scores. Th ese suitability map s were then used to simulate

locations of expansion between 2003 and 2009, and the simulated results were compared with actual 2009

palm oil locations. The simulation with all factors equally weighted was chosen as th e best model because it

resulted in the best overall Kappa scores

. 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 in the states of Johor, Pahang and Sabah, and 9 %

of the ad ditional expansion occurred on peat soils.

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

land cover categories . A comparison of the 2003 vs. 2009 spatial data suggests highly dynamic change in oil

palm area between the two years. Although there was a net expansion of 1,386,025 ha, there was only 49%

persistence of plantations during that tim e. The result was that the areas modeled as suitable for transition

we re not always the same areas that were converted by 2009, which affected the validation of the model by

increasing the disagreement due to location. Using a static land cover map from 2001 to project land cover

impacts in 2022 also introduces uncertainty because land cover is a dynamic spatial feature.

There are several possibilities for further improvements to the model that could be incorporated at a later

stage . The first would b e to verify the locations of existing oil palm plantations

. Additional improvements

could include:

 Obtaining data for additional factors that could influence the location of palm oil expansion such as

locations of palm oil processing facilities

.

 Obtaining higher resolution maps of peat soils and peat soil depths.

 Obtaining more recent land cover maps with higher resolution and greater verification

.

# 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.

Co ngalton, 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.

G alt on , F. 1892

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

Geri, Francesco, Valerio Amici, Duccio Rocchini.2011. Spatially

- based accuracy assessment of forestation

prediction in a complex Mediterranean landscape, Applied Geography, Volume 31, Issue 3, July 2 011, 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 W ork: 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 s tudy. 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.

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: Eco logical Economics of Sustainable Watershed Management.

JAI Press, Elsevier, Amsterdam.

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

Forest Degradation (REDD). Transactions in GIS 14, no. 5 (10): 6 31

- 654. doi:10.1111/j.1467

- 9671.2010.01227.x.

Manson, Steven 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 i n 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, Silvia. 2011. Methodology for predicting oil palm expansion for Indonesia and Malaysia . Prepared for

EPA.

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

- use change w ith Geomod2:
application and validation for Costa Rica. Agriculture, Ecosystems & Environment 85(1

- 3) p.191

- 203.

Pontius Jr, R G & L Schneider. 2001. Land

- use change model validation by a ROC method for the Ipswich

watershed, Massachusetts, USA. Agricultu re, 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, & H. Chen. 2006. GEOMOD m odeling. Clark University.

Pontius Jr, R. G, et al . 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 cover c hange 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.

United States Environmental Protection Agency (EPA). 2010. Regulation of Fuels and Fuel Additives: Changes

to Renewable Fuel Standard Program; Final Rule. March 26, 2010. 75 FR 14760.

United States Environmental Protection Agency (EPA), 2010 a . Renewable Fuel Standard Program (RFS2)

Regulatory Impact Analysis, February 2010, EPA

- 420

- R

- 10

- 006

.

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.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, offi ce +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 MALAYSIA, 2003

- 2022 DRAFT REPORT