Report submitted to EPA

July 2011

Spatial Modeling of Future Oil Palm 

Expansion in Malaysia, 2003 to 2022

Nancy Harris and Sean Grimland
2 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 MODEL DESCRIPTION ......................................................................................................................................................7

2.2 DATA INPUTS ...................................................................................................................................................................9

2.3 CREATION OF FACTOR MAPS.........................................................................................................................................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 CREATING SUITABILITY-FOR-CHANGE MAPS ................................................................................................................26

2.5 EVALUATING SUITABILITY MAPS USING THE ROC STATISTIC.......................................................................................27

2.6 EVALUATING MODEL PERFORMANCE USING THE KAPPA STATISTIC..............................................................................29

3.0 RESULTS AND DISCUSSION .....................................................................................................................................30

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

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

3.3 FINAL SUITABILITY MAP ...............................................................................................................................................36

3.4 PALM OIL EXPANSION, 2003-2022 ................................................................................................................................38

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

3.4.2 Expansion by Soil Type..........................................................................................................................................39

3.4.3 Expansion by Land Cover Type .............................................................................................................................40

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

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

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 represent 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 

non-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. Stratified factor map for elevation 13

Figure 5. Stratified precipitation factor map for Malaysia. 14

Figure 6. Stratified temperature factor map for Malaysia. 15

Figure 7. Stratified slope factor map for Malaysia, with slope categorized into 5% bins. 16

Figure 8. Stratified soil factor map for Malaysia. 18

Figure 9. Stratified landcover factor map for Malaysia. 20

Figure 10. Stratified factor map 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. Example 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. Suitability 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. Green 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 simulation. 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. 39
TABLES

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 oil 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 railroads (first ten bins). 21

Table 13. Distance from settlements (first ten bins). 21

Table 14. Distance from 2003 oil palm plantations (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 location decreases. 33

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

35

Table 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. Predicted palm oil expansion per land cover type between 2003 and 2022. Results displayed as for 

Table 18 above. 41
SUMMARY

The United States 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 palm 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, distance to roads, distance to rivers, distance to 

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

and in combination to calibrate the model, and the results of each candidate model run were validated against 

a 2009 map of oil palm plantations. 

The most successful models 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 second 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 was 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, approximately 5.2 million hectares of palm oil were expected 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 nearly 60% of this additional expansion will

occur in the three states of Johor (23%), Pahang (20%) and Sabah (16%). Under the 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 additional expansion predicted to 

occur as a result of the palm-only policy was on two land cover types: forest and cropland/natural vegetation 

mosaic.
1.0 BACKGROUND 

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

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

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

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

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

emission factors for various land cover conversions using IPCC methodologies. The 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 calculated internationally based on the results of the FAPRI-CARD partial 

equilibrium model (EPA 2010a). 

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, this 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 locations 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 

fuels 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 palm 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 future expansion 

occurs, such as elevation, slope, and proximity to roads. Second, we compared the factor maps against a map 

of existing palm oil plantations in 2003 to construct a series of suitability maps. This is the calibration stage. 

For each suitability map, the model assigns higher suitability values to locations that have a combination of 

characteristics similar to the land already cultivated in palm oil and low suitability values to locations that are 

less similar to existing palm oil areas. Third, each candidate suitability map was overlain with a map of existing 

plantations in the year 2009. This is the validation stage. The Relative Operating Characteristic (ROC) statistic 

was used to evaluate whether the locations with the highest values in the suitability map were also the 
locations of actual change to palm oil between 2003 and 2009. The best performing suitability map was then

used to simulate palm oil expansion between 2003 and 2009. The level of agreement between actual and 

predicted locations of palm oil in 2009 was evaluated using Kappa indices. Finally, the suitability map that 

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

of future palm oil expansion across Malaysia until the year 2022.

2.1 Model Description

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

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

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

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

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

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

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

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

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

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

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

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

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

GEOMOD.
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). 
2.2 Data inputs 

The locations of oil palm plantations in the years 2003 and 2009 were supplied in shapefile format to EPA by 

the Government of Malaysia. The shapefiles contained no metadata and the accuracy of the maps is 

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

our analysis1. 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 million 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 866,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 (37%) Table 1. 

Table 1. Total area (ha) of Malaysian palm oil plantations in 2003 and 2009 

according to two data sources.

Source 2003 ha 2009 ha 2003-2009 

Increase

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

 

1 Both 2003 and 2009 shapefiles had approximately the same proportion of mature to immature areas (3:1).
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 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. 

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

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

locations. 

 Areas 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 included existing oil palm plantations in 2003 (green) and inland water 

bodies and protected areas (black).

2.3 Creation of Factor Maps

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

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

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

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

For this analysis, eleven 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. 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 100m ‘bins’, with each bin containing pixels with values of 0-100 m, 

100-200 m, etc. 

Table 2. Factor data considered in the GEOMOD analysis.

Biogeophysical Attribute Note / Source

Elevation SRTM 90m digital elevation data v 4

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

Precipitation WorldClim database v 1.4
Biogeophysical 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 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 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

Source: 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 - 1100 0.01% 0.01% 1.84%

Figure 4. Stratified factor map for elevation

2.3.2 Precipitation 

Annual precipitation data originated from the WorldClim global climate dataset2. 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. According 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., 2005
Climatic 

Suitability

Growth 

Limitations

Mean Annual 

Precipitation

(mm)

Precipitatio

n Bin

Palm Oil Malaysia

2003 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 

Average annual temperature data originated from the WorldClim global climate dataset3. 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 map shown in Figure 6. 

 

3

http://www.worldclim.org/, Hijmans, R.J., 2005
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 Malaysia

2003 2009

Suitable 1 None 25 – 29 1 99.16% 98.35% 72.77%

Suitable 2 Slight 22 - 25

29 – 32 20.65% 1.32% 20.34%

Moderate Moderate 20 - 22

32 – 35 30.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. The slope is (in the simplest of terms) defined as the ratio 

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

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

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

into 5% bins. 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% 14.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. Stratified 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 Database4and was stratified by grouping soils into 

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

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

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

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

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

 

4

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

humid, tropical climates and often support forested areas. 

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

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

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

which results in significant GHG emissions. 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 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

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:

Reclass Name Palm Oil Malaysia

2003 2009

Acrisols 70.11% 67.27% 79.77%

Arenosols 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%

Podzols 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 dataset is produced on an annual basis starting in 2000. 

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

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

vegetation classes, 3 developed and mosaicked land classes, and three non-vegetated land classes (Table 9).

The year 2001 was chosen to represent 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 located 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 categories 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:

MODIS 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)

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

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

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

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

roads, distance from rivers, etc.), with low values close to the feature and high 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:

‘Distance 

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:

‘Distance 

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 railroads (first ten bins).

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

area in:

‘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 

area 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 rivers.

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 originated 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 located less 

than 2 km from a 2003 plantation.

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

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

model (Figure 14).

Table 14. Distance from 2003 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 combination to create a series of suitability 

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

that have the largest 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, 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 into the same 0-100 m bin in locations of existing palm oil. The suitability map is generated 

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

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

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

red (high values) are the most suitable areas and blue (low values) are the least suitable (Figure 15).

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

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

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

each cell is calculated according as:
 ( ) [∑* 

( )+

 

 

] [∑ 

 

 

]

Where:

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

a = a particular factor map

A = the number of factor maps

Wa = the weight of factor map a

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

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 were used to refine subsequent simulations until a simulation model 

achieved the best validation results possible given available time and resources. When the best model was 

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

from 2003 to 2022.

2.5 Evaluating Suitability Maps using the ROC statistic

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

suitability groups are classified as change in that scenario. In the second scenario, change is simulated only in
suitability group 1 (top 10% of suitability values), the next scenario simulates change in both suitability groups 

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

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

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

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

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

To define the ROC, the rate of true-positives (change in both reality and model) is plotted on the 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 connects 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 change 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 landscape, then the expected value of the ROC would be 0.5, as shown by the 

red line in Figure 16. 

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

use change in Malaysia.

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

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

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

(Pontius and Schneider 2001). 
2.6 Evaluating Model Performance using the Kappa statistic

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

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

assess how well model predictions agree with reality when a single quantity of change is simulated. In other 

words, the ROC score is an evaluation of the suitability map, and the Kappa indices are an evaluation of the 

prediction map that results when the suitability map is used to predict real land cover change. 

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

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

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

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

Galton (1892). The use of kappa indices in map comparison was first proposed by Congalton et al. (1983), later 

advocated by Monserud and Leemans (1992), Congalton and Green (1999), Smits et al. (1999), and Wilkinson 

(2005).

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

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

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

reality map relative to random agreement. 

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

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

Kappa indices 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 reality

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

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

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

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

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

predicting correctly the locations of land cover change.
3.0 RESULTS AND DISCUSSION

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

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

weighted factors by their individual ROC scores. 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 nearly 60% of the 

additional expansion that occurred as a result of the palm-only scenario occurred 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 which used a single biogeophysical factor to create a suitability map,

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

the best ROC score was distance from existing 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 were 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 maps 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 mainland 

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

oil palm expansion. Precipitation, land cover type, and distance from settlements and rivers were poor 

predictors of expansion when used individually, 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 overlain onto the elevation 

suitability map and show clearly 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 exceeds the quantity of pixels predicted to change, 

GEOMOD would miss these locations of change. This effect is seen in the graphs for the initial 10% of change 

(Figure 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 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?’ 
Dist. From Plantations ROC = 0.799 Elevation ROC = 0.750 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 Settlements ROC = 0.570 Precipitation ROC = 0.540

Dist. From Rivers ROC = 0.517 Land Cover Type 2001 ROC = 0.492

Figure 17. Ranking of single-factor suitability map ROC scores.
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 decreases. 

Kappa index Elevation Precipitation 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

Disagreement 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.9541 14.8060 15.4380 14.6474 16.7538

Disagreement location 13.6749 13.8230 13.1909 13.9816 11.8752
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. 

3.2 Evaluation of Multi-Factor Suitability Maps for Predicting Change

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

expansion, we improved model performance further by generating suitability maps that included multiple 

factors. First, we created a suitability map from the combination of all eleven 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 scores (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 

remaining six factors (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

Weighted 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.0816

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

Disagreement 

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 achieved the second highest ROC score (0.879) and 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. The 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 expansion 

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 expansion in the 

palm-only scenario. 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 Malaysian 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 projected 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 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.

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,764
Figure 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 type is important primarily for 

determining how much expansion occurs on histosols. In this study, histosol soils are used as a proxy for peat 
soils 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 histosol8. 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 emissions. 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 (Table 19).

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

Results displayed as for Table 18 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 Expansion 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 locations of predicted palm oil 

 

8

http://soils.usda.gov/technical/classification/taxonomy/
expansion between 2003 and 2022 were overlain with MODIS land cover types for the year 2001, most 

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

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 18 above. 

EQUAL WEIGHTS ROC WEIGHTS

Land Cover Type Control Case 

(ha)

Additional

Expansion,

Palm-Only Case 

(ha)

Control 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 broadleaf 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 shrublands 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,523 40,764 5,201,523 40,764
4.0 CONCLUSIONS

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 maps had ROC 

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

locations that were assigned high suitability scores. These suitability maps 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 the best model because it 

resulted in the best overall Kappa scores. After this model validation step, the model was run to the year 2022 

to predict locations of future palm oil expansion under a control case and a 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 additional expansion occurred on peat soils. 

The largest source of uncertainty in the analysis was the quality of the input data 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 time. The result was that the areas modeled as suitable for transition 

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

increasing the disagreement due to location. 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 be 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.

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

pp. 183 (Boca Raton, FL: Leweis)

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

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

Modeling 212: 439-449.

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

Laxenburg, Austria.

Galton, 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 2011, p.881-

890.

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

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

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

Press, Albuquerque, NM.

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

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

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

surfaces for global land areas. International Journal of Climatology 25: 1965-1978.

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.

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): 631-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 in northeastern Puerto Rico. International Journal of Climate, 19 

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

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

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

Petrova, 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 with 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. 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, & H. Chen. 2006. GEOMOD modeling. 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 change in Kanakapura region 

using GEOMOD. Journal of the Indian Society of Remote Sensing p.45-54. 
Silva, E. and K. Clarke. 2002. Calibration of the SLEUTH urban growth model for Lisbon and Porto, Portugal. 

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

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

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

Remote Sensing,20, pp. 1461-1486.

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), 2010a. 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, office +1.703.302.6504, sgrimland@winrock.org

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