Combining Radar and Optical Imagery to Map Oil Palm Plantations in Sumatra, Indonesia, Using the Google Earth Engine

Monitoring the expansion of commodity crops in the tropics is crucial to safeguard forests for biodiversity and ecosystem services. Oil palm (Elaeis guineensis) is one such crop that is a major driver of deforestation in Southeast Asia. We evaluated the use of a semi-automated approach with random forest as a classifier and combined optical and radar datasets to classify oil palm land-cover in 2015 in Sumatra, Indonesia, using Google Earth Engine. We compared our map with two existing remotely-sensed oil palm land-cover products that utilized visual and semi-automated approaches for the same year. We evaluated the accuracy of oil palm land-cover classification from optical (Landsat), radar (synthetic aperture radar (SAR)), and combined optical and radar satellite imagery (Combined). Combining Landsat and SAR data resulted in the highest overall classification accuracy (84%) and highest producer’s and user’s accuracy for oil palm classification (84% and 90%, respectively). The amount of oil palm land-cover in our Combined map was closer to official government statistics than the two existing land-cover products that used visual interpretation techniques. Our analysis of the extents of disagreement in oil palm land-cover indicated that our map had comparable accuracy to one of them and higher accuracy than the other. Our results demonstrate that a combination of optical and radar data outperforms the use of optical-only or radar-only datasets for oil palm classification and that our technique of preprocessing and classifying combined optical and radar data in the Google Earth Engine can be applied to accurately monitor oil-palm land-cover in Southeast Asia.


Introduction
Tropical forests are important ecosystems for biodiversity and carbon storage and they provide numerous material and non-material benefits to people [1]. However, these ecosystems are threatened as a result of fires and conversion for other land-uses, chiefly for the purpose of agriculture [1]. Recent global assessments on the drivers of forest degradation and loss showed that commodity-driven deforestation, which is the loss of forest cover for commodity crops, such as palm oil and soy, was the major driver for deforestation in the tropics [2,3]. This is especially pronounced in tropical Southeast Asia, where one of the major drivers for this land-use change is the expansion of the oil palm (Elaeis guineensis), an economically important vegetable oil crop for the purpose of food, oleochemical industries, and biofuel [4].
Indonesia is the largest producer of palm oil in the world, with 12 million ha of oil palm cultivated and 34 Mt exported in 80 countries, making up 48% of global palm oil export in 2017 [5]. The oil palm industry employs more than 2 million people in Indonesia and contributes to 1.5% of national gross domestic production [6]. However, rapid expansion of oil palm plantations has also led to large-scale deforestation, greenhouse gas emissions, and biodiversity loss, as well as land grabs and social conflicts in Indonesia [7][8][9], leading to many concerns about the negative effects of palm oil on people and the environment [4,[10][11][12]. Hence, it is crucial to map and monitor oil palm plantations for evaluating and managing their social and environmental effects [13][14][15].
Multiple techniques have been developed to map oil palm land-cover, using satellite remote sensing data across the tropics [16][17][18][19][20]. Optical methods for mapping oil palm land-cover rely on information derived from the phenology or image characteristics of the oil palm plantation [17]. Phenology-based methods, such as the use of temporal changes in vegetation greenness by Gutierrez-Velez et al. (2011) to detect the area deforested by large-scale oil palm expansion in the Peruvian Amazon. Image-based methods utilize spectral signatures, as well as textural information, to differentiate oil palm trees from their surroundings [21,22]. Oil palm plantations can be manually digitized from satellite images, based on the unique textural information of oil palm plantations (e.g., long rectangular blocks for industrial plantations, geometric shape of oil palm canopy, presence of roads), along with expert knowledge on the land system, or it could be automated based on spectral image analysis, which classifies pixels based on their spectral class thresholds [15]. However, it is challenging to separate oil palm plantations from other spectrally similar land-cover types (e.g., forests, rubber trees) [23], and the frequent presence of cloud cover in the tropics hinders image analysis [17].
To resolve challenges of mapping oil palm plantations in the humid tropics, radar data, which is all-weather and all-time capable, have been tested for mapping oil palm. Microwave energy scattered by vegetation depends on the size, and density, as well as the orientation and dielectric properties of elements that are comparable to the size of the radar wavelength, while optical energy reflected by vegetation depends on the leaf structure, pigmentation, and moisture [24]. Hence, radar data provide more information on the structural properties of the land, such as harvesting trails and roads, which are a prominent feature of oil palm plantations. Radar data are also able to distinguish closed canopy oil palm stands, with their large fronds supported by branchless trunks from other woody vegetation and tree cover [18]. Past studies have successfully used radar data to detect new fronts of deforestation [25], as well as tree-crop plantations in particular oil palm [17,18,26,27]. As a result, a combination of optical and radar data has been used for mapping oil palm across different geographical contexts and using different data combinations [28][29][30][31][32]. In a review of the performance of optical and radar data fusion for land-use mapping and monitoring, 28 out of 32 studies obtained higher accuracies when fusing the two datasets, compared to using only one of them [24]. More specifically for palm oil land cover classification, fusing optical and radar datasets has also resulted in higher accuracies than classification using only one dataset [33].
While the use of optical and radar imagery shows high accuracy for mapping plantation agriculture, automated mapping of oil palm plantations can be problematic and has been shown to underestimate extents of oil palm plantations, as compared to visual classification [15]. In Borneo, automated classification of oil palm was limited to closed canopy oil palm stands, while visual classification was able to detect all areas used for oil palm agriculture, regardless of the phase or condition of the oil palm plantation [15]. This difference is related to automated approaches that are only able to interpret physical characteristics of land-cover, while visual approaches integrate other land-use characteristics, using contextual information and interpreter knowledge of the area into the mapping of oil palm plantations [15].
We mapped oil palm plantations for the year 2015 to test the use of combined radar and optical data to map oil palm plantations in three provinces of Sumatra -Riau, Jambi, and South Sumatra. We compared the accuracy of our combined map with maps generated using only radar or optical data. We also compared the oil palm area of our highest accuracy map with two other oil palm land-cover products produced in 2015: Miettinen et al. (2016), which used Moderate Resolution Imaging Spectroradiometer (MODIS) imagery (250m) to conduct a semi-automated unsupervised land-cover classification with Sentinel-1 data to detect oil palm plantations, and Austin et al. (2017), which manually digitized the extents of oil palm planted using Landsat imagery (30 m) and official plantation boundaries provided by the Indonesian government. Through our study, we aim to address the following research questions: (i) does the use of radar and optical imagery for oil palm classification outperform radar-only and optical-only classification?; (ii) how does our 30 m semi-automated map of oil palm plantations compare with two maps a 250 m semi-automated map and a 30 m digitized map of oil palm plantations for the same year? The study proposes improved methods to map oil palm, using a combination of radar and optical imagery on the Google Earth Engine and a systematic comparison of previous oil palm maps to identify similarities to drive advancement of new methods.

Study Site
Our study site were the provinces of Riau (0.5333 • N, 101.4500 • E), Jambi (1.5833 • S, 103.6166 • E), and South Sumatra (2.9789 • S, 104.7584 • E), with a total area of 237,000 km 2 or 23,700,000 ha. These three provinces lie east of the Barisan mountain range, spanning an elevation gradient from 0 to 3736 m. Approximately 30% of the land is underlain by peatland close to the eastern coast, with peat depth up to 8 m [34]. These three provinces are also important oil palm producing provinces, with a total growing extent of~4 million ha, representing 36% of Indonesia's total oil palm cultivated extent in 2015 [5]. Riau has the largest extent of oil palm cultivated among the three provinces (2,400,876 ha), of which more than half are owned by smallholders (57% or 1,354,503 ha) and the rest (43% or 1,046,373 ha) are owned by industrial government or private estates. South Sumatra has the second largest extent of cultivated oil palm (952,083 ha), of which 45% (431,104 ha) are owned by smallholders and 55% (520,979 ha) are owned by industrial government or private estates. Jambi has the smallest extent of oil palm cultivated (714,399 ha), but the highest proportion of smallholder owned oil palm (63% or 450,075 ha), with the remaining 37% (264,324 ha) owned by industrial government or private estates.

Satellite Imagery
We chose optical and radar satellite data for the year 2015 to map oil palm plantations in Riau, Jambi, and South Sumatra. Mapping oil palm land-cover for the year 2015 would allow us to compare our extent of oil palm planted area with maps produced for the same year by Miettinen et al. (2016) and Austin et al. (2017), hereafter labelled as 'Miettinen2016' and 'Austin2017', respectively. Landsat 8 surface reflectance products were imported from the Google Earth Engine data catalogue for our three provinces in Sumatra. Our Landsat imagery consisted of 12 bands, including ultra-blue, blue, green, red, near infrared, and shortwave infrared 1 and 2. Advanced land observing satellite-2 phased array L-band synthetic aperture radar-2 (ALOS-2/PALSAR-2) imagery was imported from the Japan Aerospace Exploration Agency (JAXA) in Google Earth Engine's data catalogue. Data layers under ALOS-2/PALSAR-2 include horizontal transmit-vertical receive (HV) polarization and horizontal transmit-horizontal receive (HH) polarizations.

Workflow for Image Processing, Classification, and Analysis
Our workflow includes the following steps: image pre-processing, image classification, accuracy assessment, and comparison with other oil palm maps ( Figure 1). Image pre-processing for Landsat imagery removed cloud cover from the surface reflectance Landsat product to create the best available pixel image composite, while that of PALSAR imagery included speckle filtering and conversion of normalized radar cross-sections, as well as the calculation of texture measures and indices. Image classification involved the creation of image stacks using Landsat and PALSAR layers, drawing regions of interest and classifying the image using a machine learning classifier. Three stacks of data layers were created for image classification: PALSAR-only layers ('synthetic aperture radar (SAR) layer'), Landsat-only layers ('Landsat layer'), and a combination of Landsat and SAR layers ('Combined layer'). Accuracy assessments were conducted for SAR, Landsat, and Combined layers to evaluate which map produced the best oil palm extent map. We compared our best oil palm land-cover map based on overall classification accuracy with Miettinen2016 and Austin2017 (1) by examining areas of disagreement against high-resolution Google Earth Pro imagery and (2) by calculating pairwise Jaccard similarity for oil palm cover between the best map produced from our image classification and each of the two existing maps.
Remote Sens. 2020, 11, x FOR PEER REVIEW 4 of 18 drawing regions of interest and classifying the image using a machine learning classifier. Three stacks of data layers were created for image classification: PALSAR-only layers ('synthetic aperture radar (SAR) layer'), Landsat-only layers ('Landsat layer'), and a combination of Landsat and SAR layers ('Combined layer'). Accuracy assessments were conducted for SAR, Landsat, and Combined layers to evaluate which map produced the best oil palm extent map. We compared our best oil palm landcover map based on overall classification accuracy with Miettinen2016 and Austin2017 (1) by examining areas of disagreement against high-resolution Google Earth Pro imagery and (2) by calculating pairwise Jaccard similarity for oil palm cover between the best map produced from our image classification and each of the two existing maps. Overall workflow for mapping oil palm plantations in Riau, Jambi, and South Sumatra, which includes pre-processing of our phased array L-band synthetic aperture radar-2 (PALSAR-2) and Landsat 8 images, image classification using a random forest algorithm, developing accuracy assessments for our oil palm land cover map, and comparing our oil palm land cover map with two other oil palm map sources.

Image Pre-Processing for Landsat Imagery
To remove cloud cover in Landsat images, Landsat image composites were based on a time series, where the best available observation from 40   Overall workflow for mapping oil palm plantations in Riau, Jambi, and South Sumatra, which includes pre-processing of our phased array L-band synthetic aperture radar-2 (PALSAR-2) and Landsat 8 images, image classification using a random forest algorithm, developing accuracy assessments for our oil palm land cover map, and comparing our oil palm land cover map with two other oil palm map sources.

Image Pre-Processing for Landsat Imagery
To remove cloud cover in Landsat images, Landsat image composites were based on a time series, where the best available observation from 40 Landsat 8 images of the path 127 and rows 59 and 60 from the year 2015 were selected through pixel-based image compositing ( Figure 2). The image compositing the Google Earth Engine script was adapted from De Alban et al. (2018) and applied to the Landsat 8 surface reflectance product. The image compositing script selected the best available pixel in Landsat 8 images. We applied the median compositing method, which takes the median value of each band and therefore is very resistant to outliers. We used a hybrid cloud shadow masking method, including a function of mask for cloud detection, a cloud score function for cloud masking, and a temporal dark outlier mask (TDOM) for cloud shadow. The threshold for masking clouds was <5% and the shadows z-score = −1. We used the tier 1 surface reflectance Landsat product and Landsat 8 OLI/TIRS sensors.
Remote Sens. 2020, 11, x FOR PEER REVIEW 5 of 18 surface reflectance product. The image compositing script selected the best available pixel in Landsat 8 images. We applied the median compositing method, which takes the median value of each band and therefore is very resistant to outliers. We used a hybrid cloud shadow masking method, including a function of mask for cloud detection, a cloud score function for cloud masking, and a temporal dark outlier mask (TDOM) for cloud shadow. The threshold for masking clouds was <5% and the shadows z-score = −1. We used the tier 1 surface reflectance Landsat product and Landsat 8 OLI/TIRS sensors.

Image Pre-Processing for SAR Imagery
Our 2015 SAR mosaic was obtained from the seamless yearly mosaic available on the Google Earth Engine database and clipped around Sumatra, our area of interest. We applied a 3 × 3 Refined Lee filter on each channel, HH and HV, to smoothen the image, reduce the speckle effects in the raw SAR imagery, and conserve the information of the edges and linear features between two land-cover types [47]. The digital numbers of the HH and HV images were converted into sigma-naught values, With the processed Landsat image composite, we calculated five vegetation indices, including the normalized difference vegetation index (NDVI) [35,36] the enhanced vegetation index (EVI) [37,38], the soil adjusted total vegetation index (SATVI) [39], the normalized difference tillage index (NDTI) [40], and the normalized difference water index (NDWI) [41,42], also called the modified NDVI or land surface water index (LSWI; Table 1). The NDVI and LSWI can separate land-cover types, such as forests and croplands [43]. The NDTI, SATVI, and NDWI have been shown to discriminate plantations from forests [32]. Other studies integrating radar and optical satellite imagery have shown that the NDVI, NDWI, and EVI improve the detection of deciduous rubber plantations and forests [32,[44][45][46]. Table 1. Definition of vegetation indices calculated from Landsat imagery. Near-infrared band (NIR); red band (RED), blue band (BLUE); shortwave infrared band (SWIR).

Vegetation Indices Formulas Characteristics
Normalized Difference Vegetation Index (NDVI) NDVI = NIR -RED NIR+RED NDVI is widely used as a vegetation metric because it is sensitive to change in greenness and therefore green vegetation biomass.
EVI improves the vegetation signal especially in regions with dense green vegetation. Further, it is little influenced by the soil and atmosphere.
SATVI can detect both senescent and green vegetation.
NDTI is sensitive to tillage intensity and crop residue cover.
NDWI is sensitive to changes in vegetation liquid water content. It is higher for green than dry vegetation.

Image Pre-Processing for SAR Imagery
Our 2015 SAR mosaic was obtained from the seamless yearly mosaic available on the Google Earth Engine database and clipped around Sumatra, our area of interest. We applied a 3 × 3 Refined Lee filter on each channel, HH and HV, to smoothen the image, reduce the speckle effects in the raw SAR imagery, and conserve the information of the edges and linear features between two land-cover types [47]. The digital numbers of the HH and HV images were converted into sigma-naught values, following Equation (1).
σ 0 is the radar backscatter per unit area; DN is the digital number, and CF is the calibration factor with a given value of -83.0 dB for the SAR mosaic [47]. Speckle filtering and conversion to normalized radar cross-section was implemented in the Google Earth Engine based on a script adapted from the SNAP 3.0 toolbox of the European Space Agency (Figure 2).
To complement the information obtained through the SAR channels, HH and HV, we calculated texture measures, which are the relationship between neighbouring pixels and their grey tones in a defined window. The second-order texture measures are the image statistical distribution of grey tone and describe how often a grey tone appears in a specific spatial relationship to another grey tone. To extract textural information, it is necessary to calculate grey-level co-occurrence matrices (GLCMs). GLCM statistic values allow us to understand the organization of pixels and show variation among pixels in different parts of a picture [48]. They outline the angular relationships and distance between adjacent pixels [48]. Texture measures add information to discriminate land-covers and enhance classification accuracy of broad land-cover types [49]. They can predict vegetation properties, such as height, structure, and canopy cover [48].
The GLCM texture measures we calculated included angular second moment (ASM), which represents the orderliness of a picture, the average (AVG), which is the frequency of occurrence of one grey tone in combination with another certain grey tone, the contrast (CON), which is the difference between pixels, the correlation (COR), which is the predictability and linear relationship between two neighbouring pixels, entropy (ENT), which is the degree of disorder among pixels, the inverse difference moment (IDM), which inversely correlated with the contrast, and the variance (VAR), which is the dispersion around the average of pixel values within a defined window. Texture measures were calculated from the average of the directional bands within a 3x3 window using the glcmTexture function available in the Google Earth Engine [50].
To add layers of information and enhance classification, we calculated six indices (Table 2) from the raw synthetic aperture radar backscattering coefficient, including average (AVE), difference (DIF), and HH and HV channels' simple ratios (RAT1 and RAT2). These four indices yielded good accuracies when differentiating forest to non-forest and coconut palm [51], oil palm, and rubber [33]. In addition, they improve mapping of broad land-cover types, such as cropland, forest, built-up areas, and water [26,45]. The DIF showed a good result to discriminate plantations, such as oil palm and rubber [18]. The two other indices calculated were the normalized difference index (NDI), which is useful to detect new deforestation events [25], and the NL index, which can enhance overall accuracies for land-cover classification using L-band SAR data [52].

Creating Image Stacks
The final Landsat image stack contained 12 layers, comprising 5 indices (NDVI, EVI, SATVI, NDTI, and LSWI) and 7 bands (blue band (BLUE), green band (GREEN), near-infrared band (NIR), red band (RED), shortwave infrared bands 1 and 2 (SWIR1 and SWIR2), TIR). The final SAR image stack contained 24 layers (i.e., HH and HV with their corresponding textures, AVE, DIF, NDI, NLI, RAT1, and RAT2). The final Combined stack contained 36 layers, including indices, textures, and bands from the two previous stacks. Despite the importance of image features to increase class separability and classification accuracy, increasing the number of features does not necessarily translate in higher accuracies [53]. Each feature was chosen according to its potential to improve oil palm classification, as cited in the literature.

Regions of Interest and Sampling Design
We defined six land-cover classes, which included forest (FOR), forest shrub mosaic (FSM), oil palm (OPM), bare ground-crop (BGC), built-up area (URB), and water (WTR) (Supplementary Table  S1). We delineated regions of interest (ROIs) for these land-cover classes by visual interpretation of high-resolution images in Google Earth Pro from 01 January 2015 to 31 December 2015. We used features such as location, size, shape, tone, colour, shadow, texture, and pattern to identify and differentiate our land-cover classes in Google Earth Pro [54] (Supplementary Table S1). To ensure representability of our sample ROIs, we obtained between one and four ROIs for each land-cover class in each regency for each of our provinces of interest; Riau, Jambi, and South Sumatra. This sampling scheme ensured a spread of ROIs across our different satellite image tiles. The final sample of 788 polygons was imported as a feature collection in the Google Earth Engine.
We used the randomColumn and sampleRegions functions in the Google Earth Engine to separate the ROIs into testing and training polygons. The randomColumn function gives a random number between 0 and 1 to each polygon, which would be appointed as training data when they have a value of ≤0.70 and as testing data when they have a value of >0.70. The same selection process was repeated at the pixel level. Using the same Google Earth Engine functions, every pixel within each training and testing polygon was assigned a random number between 0 and 1. Pixels from testing polygons with an assigned value of ≤0.70 were segregated as testing pixels, while pixels from training polygons with an assigned value of ≤0.70 were appointed as training pixels. In total, the sample for the SAR classification had 191,663 training pixels and 108,816 testing pixels, the Landsat sample had 133,064 training pixels and 75,502 testing pixels, and the Combined sample had 133,039 training pixels and 75,501 testing pixels. This two-level random sampling design prevented training and testing data to be spatially dependent [55].

Image Classification
We used a non-parametric random forest (RF) machine learning classifier for supervised land-cover classification. RF is an ensemble of decision trees that, based on the data they receive, vote for the most popular class for an input pixel [56]. In general, non-parametric machine-learning classifiers, such as RF or support vector machine (SVM), yield higher accuracies than parametric classifiers. RF can improve accuracy of land-cover classification [57] and detect plantations like oil palm [54]. We used the RandomForest function included in the Google Earth Engine and trained the machine learning algorithm with the training data obtained from each of our three image stacks: Landsat, SAR, and Combined. The following parameters were used as inputs for the algorithm: 100 for the number of decision trees created per class, 10 for the minimal length of a terminal node, and 0.50 as the input fraction to bag per tree. Following classification, we applied a majority filter to homogenise the classified image and remove small clusters of pixels. Isolated pixels were replaced by another pixel value based on the majority of their contiguous neighbouring pixels.

Accuracy Assessments
We evaluated our classification results using standard accuracy assessment metrics as the error matrix, overall accuracy, user's, and producer's accuracies [58]. Since we were more interested in oil palm mapping, we focused on the producer's and user's accuracies for oil palm land-cover class. We identified the best map as the one with the highest overall, producer's, and user's accuracies. To calculate the metrics above, we used Google Earth Engine functions errorMatrix(), accuracy(), consumersAccuracy(), producersAccuracy(), and kappa().

Comparing Classified Oil Palm Extent with Other Oil Palm Maps for 2015
We extracted the oil palm mapped under our best classified map and compared this with two other oil palm land-cover products produced in 2015 from Miettinen2016 and Austin2017 using difference images. Our difference images allowed for the visualization of areas of agreement, where both maps detected oil palm or both maps detected no oil palm, as well as areas of disagreement, where one map detected oil palm but not the other [59]. We first used the function r.cross from GRASS GIS to create the two difference images, comparing our classified map with Miettinen2016 and Austin2017. The raster layers were polygonised to calculate the areas of agreement and disagreement in each map comparison. These values of agreement and disagreement were used to calculate the Jaccard similarity coefficient (J) using the following equations: Under Equation (2), the Jaccard coefficient (J) is calculated by dividing the areas shared by maps X and Y by the areas belonging to both or either maps. The Jaccard coefficient values vary between 0 and 1. A value of 0 means that the two compared maps have no areas in common, while a value of 1 means that the two maps are identical and share the same area.
To investigate how land-cover classes were classified differently under our extents of disagreement between our classified oil palm extent and Miettinen2016, as well as our classified oil palm extent and Austin2017, we extracted the areas of disagreement in each of our difference images. These areas of disagreement refer to pixels identified as oil palm in one map extent but not the other (i.e., when our map identified oil palm but Miettinen2016 or Austin2017 did not and where Miettinen2016 or Austin2017 identified oil palm where our map did not). For each extent of disagreement, we randomly selected 100 points.
We visually inspected each of these points in Google Earth Pro, using high-resolution images from 2015, and recorded the land-cover for each point. We counted the number of points for each land-cover class and grouped them according to whether the point was identified as oil palm under our map or either Miettinen2016 or Austin2017. Since the area of oil palm identified under Miettinen2016, Austin2017, and our maps differed, the random sample of 100 points from our disagreement extents reflected this difference in area. Hence, we weighted our count data of each land-cover class with the total number of random points drawn from each map. For example, 56 pixels were correctly identified as oil palm in the Combined map and 15 in Miettinen2016; however, after weighting these count data, both maps had 71% of their pixels correctly identified as oil palm. We present the proportion of points that were accurately identified as oil palm in our Google Earth Pro image and the proportion of points which fell under other land-cover classes, such as forest shrub mosaic or bare ground cover land-classes.

Comparing Landsat, SAR, and Combined Maps
Our Combined (Landsat + SAR) map had the highest overall accuracy score of 84%, followed by our Landsat map (78%) and SAR map (74%) (Table 3, Figure 3). The Combined map also had the highest producer's and user's accuracies for oil palm land-cover class (84% and 90%, respectively), as compared to Landsat (70% and 62%) and SAR (81% and 66%) maps ( Table 3). The error matrix for our Landsat classified image showed a high number of oil palm pixels misclassified as forest and forest shrub mosaic land-cover classes (Figure 4). Under our SAR classified image, oil palm pixels were confused with mostly bare ground-crop, followed by forest shrub mosaic and forest land-cover classes (Figure 4). The error matrix for our Combined classified image showed that oil palm pixels were confounded with bare ground-crop and forest shrub mosaic (Figure 4). Overall, combining SAR and Landsat imageries outperformed single imagery (SAR-only and Landsat-only) classification in mapping oil-palm cover due to the superior producer 's and user's accuracy of the former technique (Table 3).

Comparing Oil Palm Extents with Miettinen2016 and Austin2017
The extent of oil palm mapped in our Combined classified image for Riau, Jambi, and South Sumatra was 4,418,176 ha, which was 8.6% higher than the oil palm extent from the Indonesian Ministry of Tree Crop Statistics (4,067,358 ha) for the year 2015 [5]. Our mapped extent of oil palm was 47% higher than the extent of oil palm mapped by Miettinen2016 (2,997,943 ha) and 22% higher than the extent of oil palm mapped by Austin2017 (3,608,815 ha). Table 3. Producer's (%) and user's (%) accuracy of Combined, Landsat, and SAR classified maps for all land-cover classes. Producer's accuracy (PA); user's accuracy (UA). Values in bold represent the highest producer's and user's accuracy for oil palm land-cover and highest overall accuracy.    The Jaccard coefficient between our Combined map and the Miettinen2016 map was 0.55, which was higher than the Jaccard coefficient between our Combined map and the Austin2017 map (0.40).

Land-Cover Class Combined Landsat SAR
The difference images in Figure 5 show extents of agreement and disagreement for the Combined-Miettinen2016 and Combined-Austin2017 maps. These difference images show the locations of agreements, where oil palm was detected in both maps, and locations of disagreements, where oil palm was detected in one map but not the other. In the Combined-Miettinen2016 comparison, the total extent of disagreement was 2,155,972 ha, of which 83% (1,788,479 ha) stemmed from our Combined map, identifying oil palm in areas where Miettienen2016 did not, and 17% (367,496 ha) were from areas of disagreement where Miettinen2016 identified oil palm and where our Combined map did not. In the Combined-Austin2017 comparison, the total extent of disagreement was 3,444,846 ha, of which 62% (2,126,000 ha) came from our Combined map, identifying oil palm in areas where Austin2017 did not, and 38% (1,318,846 ha) were from areas of disagreement where Austin2017 identified as oil palm and where our Combined map did not. These extents were largely concentrated in South Sumatra. Approximately 62% (2,125,999 ha) of the extent of disagreement was a result of our Combined map detecting oil palm where Austin2017 did not, and these areas were found mostly on the eastern coasts of Riau and Jambi.

Comparing Oil Palm Extents with Miettinen2016 and Austin2017
The extent of oil palm mapped in our Combined classified image for Riau, Jambi, and South Sumatra was 4,418,176 ha, which was 8.6% higher than the oil palm extent from the Indonesian Ministry of Tree Crop Statistics (4,067,358 ha) for the year 2015 [5]. Our mapped extent of oil palm was 47% higher than the extent of oil palm mapped by Miettinen2016 (2,997,943 ha) and 22% higher than the extent of oil palm mapped by Austin2017 (3,608,815 ha). The Jaccard coefficient between our Combined map and the Miettinen2016 map was 0.55, which was higher than the Jaccard coefficient between our Combined map and the Austin2017 map (0.40).
The difference images in Figure 5 show extents of agreement and disagreement for the Combined-Miettinen2016 and Combined-Austin2017 maps. These difference images show the locations of agreements, where oil palm was detected in both maps, and locations of disagreements, where oil palm was detected in one map but not the other. In the Combined-Miettinen2016 comparison, the total extent of disagreement was 2,155,972 ha, of which 83% (1,788,479 ha) stemmed from our Combined map, identifying oil palm in areas where Miettienen2016 did not, and 17% (367,496 ha) were from areas of disagreement where Miettinen2016 identified oil palm and where our Combined map did not. In the Combined-Austin2017 comparison, the total extent of disagreement was 3,444,846 ha, of which 62% (2,126,000 ha) came from our Combined map, identifying oil palm in areas where Austin2017 did not, and 38% (1,318,846 ha) were from areas of disagreement where Austin2017 identified as oil palm and where our Combined map did not. These extents were largely concentrated in South Sumatra. Approximately 62% (2,125,999 ha) of the extent of disagreement was a result of our Combined map detecting oil palm where Austin2017 did not, and these areas were found mostly on the eastern coasts of Riau and Jambi. Of the 100 random sample points for our Combined-Miettinen2016 comparison, 79 points were drawn from the Combined map and 21 points were drawn from the Miettinen2016 map. Whereas, for our Combined-Austin2017 map, 58 points were drawn from the Combined map, while 42 points were drawn from the Austin2017 map. We calculated the weighted percentage of points for each Of the 100 random sample points for our Combined-Miettinen2016 comparison, 79 points were drawn from the Combined map and 21 points were drawn from the Miettinen2016 map. Whereas, for our Combined-Austin2017 map, 58 points were drawn from the Combined map, while 42 points were drawn from the Austin2017 map. We calculated the weighted percentage of points for each land-cover class using the total number of points identified under each map above ( Figure 6). The majority of the points identified as oil palm in Combined (70.8%) and Miettinen2016 (71.4%) were verified as oil palm under Google Earth Pro ( Figure 6). About 15.2% of the oil palm points under the Combined map were forest shrub mosaic, and 13.9% of these points were identified as bare ground crop with Google Earth Pro. For the Miettinen2016 map, 19% of the oil palm points were forest shrub mosaic, 4.8% were bare ground cover, and another 4.8% identified as urban extents under Google Earth Pro ( Figure 6). Under the Combined-Austin2017 comparison, the majority of the points identified as oil palm in Combined (65.5%) were verified as oil palm under Google Earth Pro, while 28.5% of oil palm points under Austin2017 were verified as oil palm under Google Earth Pro ( Figure 6). About 38.1% of the oil palm points under Austin2017 were forest shrub mosaic, 28.6% were bare ground cover, and 4.8% were forest under Google Earth Pro. Under the Combined map, 22.4% of the oil palm points were bare ground cover, while 12.1% were forest shrub mosaic under Google Earth Pro ( Figure 6).

Discussion
Our results showed that combining optical and radar datasets improved the accuracy of oil palm land-cover classification. This corroborates with results from other recent studies that use both optical and radar data to map oil palm plantations in Myanmar [32,33] and Indonesia [30]. Moreover, we developed a code to pre-process radar and optical datasets, as well as classify land-cover entirely in the Google Earth Engine. The improved accuracy for detecting oil palm plantations with the use of optical and radar data highlights the importance of multisensory and data fusion techniques, which integrate multiple data features to improve map detail and accuracy [60]. This is especially important in the tropics, where oil palm plantations expand most and where optical-based mapping is limited due to high extents of cloud cover. In our study, the use of only Landsat data underestimated the The majority of the points identified as oil palm in Combined (70.8%) and Miettinen2016 (71.4%) were verified as oil palm under Google Earth Pro ( Figure 6). About 15.2% of the oil palm points under the Combined map were forest shrub mosaic, and 13.9% of these points were identified as bare ground crop with Google Earth Pro. For the Miettinen2016 map, 19% of the oil palm points were forest shrub mosaic, 4.8% were bare ground cover, and another 4.8% identified as urban extents under Google Earth Pro ( Figure 6). Under the Combined-Austin2017 comparison, the majority of the points identified as oil palm in Combined (65.5%) were verified as oil palm under Google Earth Pro, while 28.5% of oil palm points under Austin2017 were verified as oil palm under Google Earth Pro ( Figure 6). About 38.1% of the oil palm points under Austin2017 were forest shrub mosaic, 28.6% were bare ground cover, and 4.8% were forest under Google Earth Pro. Under the Combined map, 22.4% of the oil palm points were bare ground cover, while 12.1% were forest shrub mosaic under Google Earth Pro ( Figure 6).

Discussion
Our results showed that combining optical and radar datasets improved the accuracy of oil palm land-cover classification. This corroborates with results from other recent studies that use both optical and radar data to map oil palm plantations in Myanmar [32,33] and Indonesia [30]. Moreover, we developed a code to pre-process radar and optical datasets, as well as classify land-cover entirely in the Google Earth Engine. The improved accuracy for detecting oil palm plantations with the use of optical and radar data highlights the importance of multisensory and data fusion techniques, which integrate multiple data features to improve map detail and accuracy [60]. This is especially important in the tropics, where oil palm plantations expand most and where optical-based mapping is limited due to high extents of cloud cover. In our study, the use of only Landsat data underestimated the extent of oil palm plantations due to confusion with spectrally similar land-cover classes, such as forest and forest shrub mosaic, while the use of only SAR data overestimated the extent of oil palm plantations, largely as a result of confusion between bare ground cover and oil palm (Figure 4). In general, classification using SAR data overlooks new oil palm plantations, because they have similar backscatter as herbaceous vegetation and bare land [17].
Our combined Landsat and SAR approach produced an oil palm extent of 4,418,176 ha, which was closest to the oil palm area reported by the Indonesian Ministry of Tree Crop Statistics (4,067,358 ha), as compared to Miettinen2016 (2,997,943 ha) and Austin2017 (3,608,815 ha) for the provinces of Riau, Jambi, and South Sumatra. However, disagreement maps between our combined approach and Miettinen2016 and Austin2017 showed large extents of non-overlap in the identification of oil palm land-cover, with our Combined map showing a higher level of agreement with the Miettinen2016 map as compared to the Austin2017 map. These differences are mainly the result of different methods used to define and classify oil palm land-cover in these studies.
Our semi-automated supervised 30 m Combined map was at a higher resolution compared to Miettinen2016, which used Sentinel-1 C-band radar data (20 m) that was resampled to 250 m to match MODIS images [30]. This difference in spatial resolution and Miettinen2016's definition of palm plantations as contiguous closed canopy plantations larger than 1 km 2 could be one reason for the difference in oil palm area (Combined map had 47% larger extent), as well as the differences in misclassification, where a higher proportion of oil palm points were forest land-cover under Miettinen2016, as compared to a higher proportion of oil palm points as bare ground cover under the Combined map ( Figure 6). While some pixels from Miettinen2016 could be coconut and sago palm plantations, the authors visually inspected and removed any potential misclassified areas to avoid classification errors. This was not carried out in our study. Given that coconut and sago plantations are prominent along the coasts of Sumatra (Lee J.S.H. pers. obs.), the large extents of disagreement between our Combined map and Miettinen2016 on the east coast of Sumatra ( Figure 5A) are likely coconut and sago plantations misclassified as oil palm in our Combined map. Our results support previous studies that found higher accuracies using supervised methods [61][62][63].
Austin2017 manually digitized oil palm by visually interpreting Landsat imagery, using rectilinear patterns and context indicators, such as road networks, to identify oil palm plantations [64]. Only largescale oil palm plantations that were organized in a grid pattern and associated with physical structures, such as roads and mills, were mapped. In addition, recently cleared extents proximal to existing plantations were included. The focus of Austin2017 on large-scale oil palm plantations and exclusion of small-scale plantations could be the reason why the Combined map detected a larger extent of oil palm area (22% higher than Austin2017). The focus on land-use by Austin2017's visual interpretation approach is an alternative to automated detection methods, such as our Combined map, which limits the interpretation of physical data to that of land-cover [15]. Hence, the differences in misclassifications of land-cover under Google Earth Pro, where Austin2017 identified a higher number of forest shrub mosaic and bare ground cover as oil palm, may not be so much a misclassification of oil palm land-cover under Austin2017's approach, but rather may reflect the different development phases of the plantation or the condition of the crop (e.g., immature, diseased) due to Austin2017's use of auxiliary information and knowledge of the area. This difference was also demonstrated using automated and visual interpretation approaches to oil palm mapping in Borneo by Miettinen et al. (2017) and Gaveau et al. (2016), respectively [15].

Conclusions
Our study shows the usefulness of the Google Earth Engine for processing large satellite datasets and producing classification results of land-cover change at regional scales. Specifically, our study highlights the advantages of combining radar and optical data to improve the accuracy of oil palm land-cover classification. The accessibility of optical and radar datasets on the Google Earth Engine, along with publicly available documentation of code, allows for collaboration and cooperation among different users to conduct land-use and land-cover change for environmental monitoring [65]. Our code for pre-processing Landsat and SAR data and conducting land-cover classification is available under the Supplementary Methods and can be used to improve semi-automated techniques to map and monitor oil palm expansion across the tropics [66]. However, our study also highlights limitations of automated mapping techniques for oil palm plantations and seconds the conclusion by Miettinen et al. (2018) about the need to combine use of visual and automated oil palm mapping approaches for monitoring oil palm expansion. Contextual factors, such as the presence of other palm plantations and the phase of development of oil palm plantations, are some examples of information that cannot be inferred using satellite imagery alone and revealed the differences between our combined optical and radar automated approach and mixed visual and automated mapping approaches from Miettinen2016 and Austin2017.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/7/1220/s1: Table S1: Description of 5 land-cover classes used in image classification; Table S2: Number of training and testing pixels for each image stack; Table S3: Error matrices for Landsat classified image; Table S4: Error matrices for SAR classified image; Table S5: Error matrices for Combined classified image; Table S6: Identification reference of the 40 Landsat images used in the classification; Table S7: Producer's (%) and user's (%) accuracy of Combined, Landsat, and SAR land-cover maps, classified with a Classification and Regression Tree (CART) algorithm. Code for data pre-processing and land-cover classification is available in Github repository (https://github.com/thuansarzynski/GEE_CombinedLandsatSAR_oilpalm).