Developing Multi-Source Indices to Discriminate between Native Tropical Forests, Oil Palm and Rubber Plantations in Indonesia

: Over the last 18 years, Indonesia has experienced signiﬁcant deforestation due to the expansion of oil palm and rubber plantations. Accurate land cover maps are essential for policymakers to track and manage land change to support sustainable forest management and investment decisions. An automatic digital processing (ADP) method is currently used to develop land cover change maps for Indonesia, based on optical imaging (Landsat). Such maps produce only forest and non-forest classes, and often oil palm and rubber plantations are misclassiﬁed as native forests. To improve accuracy of these land cover maps, this study developed oil palm and rubber plantation discrimination indices using the integration of Landsat-8 and synthetic aperture radar Sentinel-1 images. Sentinel-1 VH and VV difference (>7.5 dB) and VH backscatter intensity were used to discriminate oil palm plantations. A combination of Landsat-8 NDVI, NDMI with Sentinel-1 VV and VH were used to discriminate rubber plantations. The improved map produced four land cover classes: native forest, oil palm plantation, rubber plantation, and non-forest. High-resolution SPOT 6/7 imagery and ground truth data were used for validation of the new classiﬁed maps. The map had an overall accuracy of 92%; producer’s accuracy for all classes was higher than 90%, except for rubber (65%), and user’s accuracy was over 80% for all classes. These results demonstrate that indices developed from a combination of optical and radar images can improve our ability to discriminate between native forest and oil palm and rubber plantations in the tropics. The new mapping method will help to support Indonesia’s national forest monitoring system and inform monitoring of plantation expansion.


Introduction
The increasing global demand for agricultural products has driven the rapid expansion of oil palm and rubber plantations. Globally, between 2000 and 2018, oil palm plantations increased by >220% from 7 to 22.5 Mha, while rubber plantations increased by 45% from 5.5 to 8 Mha [1]. As a consequence, considerable areas of native forest have been converted to agricultural land use. In Indonesia, between 2001 and 2016, 43% (3.96 Mha) of the total 9.2 Mha of native forest loss was due to plantation expansion [2]. This native forest clearing has adversely impacted ecosystems, biodiversity, and the terrestrial carbon cycle [3,4].
Remote sensing is widely used to derive information on land cover and land cover change for natural resources management and environmental applications [5]. Among the remote sensing data used, optical images such as those recorded by the Landsat satellites have been particularly important in supporting long-term forest monitoring of land cover change (>40 years) [6][7][8]. Currently, Indonesia's National Forest Monitoring The study area is situated across four regencies in Kalimantan, namely Barito Timur, Tabalong, Balangan and Paser (0 • 50 20.32" S-2 • 31 56.61" S and 114 • 54 49.33" E-116 • 37 22.46" E) covering an area of over 1.9 million ha [33][34][35]. The terrain in these regions is variable, from swampy plains to hilly lowlands ranging from sea level up to 1550 m. The land cover consists of 41% native forests, 14% oil palm plantations, and 6% rubber plantations, with the remaining 39% accounted for by non-forested areas such as settlements, paddy fields, mixed agriculture, shrubs, grasslands, and mining [33][34][35].
The native forest is dominated by species of the Dipterocarpaceae family, mainly the genus of Shorea, such as Shorea johorensis and Shorea leprosula, and Dipterocarpus, such as Dipterocarpus gracilis, at an elevation below 700 m, and the high-quality timber species Agathis endertii is mostly found at an altitude of 1450-1600 m [36][37][38]. These endemic trees usually have large trunks (up to 122 cm in diameter) [36]. Due to past logging history, these forests are classified as secondary forests [9].
The oil palm and rubber plantations in the tropics are evergreen trees [39,40]. Oil palm trees (Elaeis guineensis) have single branchless trunks of 12-17 m height with eight fronds forming a rank in succession. Each frond is 8 m long and consists of 250-300 leaves that are 1.3 m long with 6 cm broad length. Oil palm grows best at a temperature range of 24-28 • C and mean annual rainfall of 2000 mm per year or greater [39]. Topographically, oil palm is mainly planted below 200 m elevation and on flatter areas with slope <20%.
Rubber trees (Hevea brasiliensis) have cylindrical trunks that branch out at the top, forming a dense canopy. The leaves are in spirals and have three leaflets [41]. The leaves are 3-20 cm long for the main stalk and 3-10 cm for the leaflets. Each leaflet is 5-35 cm long and 2.5-12.5 cm broad. Similar to oil palm, rubber plantations are mostly grown in lowlands at an elevation of less than 600 m, where the temperature ranges between 20 and 30 • C, annual rainfall rate is >1500 mm per year, and the slope is below 30% [40].
Shrub and grassland areas that are included within the non-forest class are dominated by a mix of shrub species (up to 1.5 m tall), such as Hyptis capitate, Ageratum conyzoides L., and Eleocharis dulcis, and are usually dense in both drylands or wetlands [42,43]. Sample photos of the native forest, oil palm and rubber plantations, shrubs and grasslands in the study area are shown in Figure 1.

Satellite Data and Pre-Processing
The datasets used to develop plantation discrimination indices were derived from Landsat-8 and Sentinel-1 images. The Landsat-8 images were downloaded from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/, accessed on 3 July 2019), and Sentinel-1 images were downloaded from the Alaska Satellite Facility (ASF) (https://www.asf.alaska.edu/, accessed on 3 July 2019). The sensor details, preprocessing, and image details are as follows.

Satellite Data and Pre-Processing
The datasets used to develop plantation discrimination indices were derived from Landsat-8 and Sentinel-1 images. The Landsat-8 images were downloaded from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/, accessed on 3 July 2019), and Sentinel-1 images were downloaded from the Alaska Satellite Facility (ASF) (https://www.asf.alaska.edu/, accessed on 3 July 2019). The sensor details, pre-processing, and image details are as follows. • Landsat-8 The Landsat-8 satellite has a revisit time of 16 days and consists of two sensors: the Operational Land Imager (OLI) sensor with nine bands ranging from visible to short-wave infrared (SWIR) recorded at 30 m resolution and the Thermal Infrared Sensor (TIRS) with two thermal bands recorded at 100 m resolution. For this study, Landsat-8 OLI images were used. The pre-processing steps to produce mosaiced Landsat-8 images at 25 m resolution with six multi-spectral image bands (blue, green, red, near infrared, short-wave infrared 1, and short-wave infrared 2) were based on the LCCA method outlined in [10] and consist of scene selection; geometric, radiometric, and terrain correction; cloud masking; and mosaicking of the individual Landsat-8 images. The current study used 15 Landsat-8 scenes with date acquisitions provided in Table 1. Only scenes with less than 50% cloud cover were selected. •

Sentinel-1
The Sentinel-1A and Sentinel-1B satellites both include a C-band SAR sensor with frequency ~5.4 GHz (wavelength ~5.6 cm) to produce multi-polarized datasets of single (VV or HH) and dual polarization (HH + HV or VV + VH). The four different polarization

•
Landsat-8 The Landsat-8 satellite has a revisit time of 16 days and consists of two sensors: the Operational Land Imager (OLI) sensor with nine bands ranging from visible to short-wave infrared (SWIR) recorded at 30 m resolution and the Thermal Infrared Sensor (TIRS) with two thermal bands recorded at 100 m resolution. For this study, Landsat-8 OLI images were used. The pre-processing steps to produce mosaiced Landsat-8 images at 25 m resolution with six multi-spectral image bands (blue, green, red, near infrared, short-wave infrared 1, and short-wave infrared 2) were based on the LCCA method outlined in [10] and consist of scene selection; geometric, radiometric, and terrain correction; cloud masking; and mosaicking of the individual Landsat-8 images. The current study used 15 Landsat-8 scenes with date acquisitions provided in Table 1. Only scenes with less than 50% cloud cover were selected.
The Sentinel-1A and Sentinel-1B satellites both include a C-band SAR sensor with frequency~5.4 GHz (wavelength~5.6 cm) to produce multi-polarized datasets of single (VV or HH) and dual polarization (HH + HV or VV + VH). The four different polarization modes of Sentinel-1 satellites are VV (vertical transmit-vertical receive), HH (horizontal transmit-horizontal receive), HV (horizontal transmit-vertical receive), and VH (vertical transmit-horizontal receive). The dual-polarization of VV + VH bands was used in this study to analyze the backscatter intensity variations for all land cover classes. Each of the Sentinel-1 satellites can acquire data for the same ground location every 12 days or every six days when combining data from the two satellites. The current study used three scenes of Sentinel-1A Interferometric Wide (IW) swath mode provided in Level-1 Ground Range Detected (GRD) format in descending mode, with a pixel size of 10 m. The acquisition dates for the Sentinel-1 images are shown in Table 1.
The Sentinel-1 IW GRD images were multi-look (averaged) and ground range detected. However, these images required pre-processing to more accurately match ground coordinates and derive the backscatter intensity referred to as sigma nought (σ • ), with units of decibel (dB). The pre-processing steps were conducted using the SNAP open-source software; the steps consisted of applying an update of precise orbit acquisition information, radiometric calibration, and terrain correction to remove terrain effects using the digital elevation model developed from the SRTM (Shuttle Radar Topography Mission). This was followed by mosaicking the images to produce a 25 m resampled image.

Field Data and Training Sample Collection
Samples from field data and high-resolution images are required for accurate interpretation of land cover. Ground-truth information for determining the best SAR or optical index and relevant index threshold for discrimination of plantations were obtained from SPOT 6/7 images (resolution 1.5 m) and field measurements. The training samples were collected from a group of image regions of interest (ROI); 22 ROIs (12,994 pixels) for native forest, 33 ROIs (11,124 pixels) for rubber, 25 ROIs (13,008 pixels) for oil palm, and 26 ROIs (7648 pixels) for non-forest such as shrub in dryland and wetland regions.
A total of 32 field plots were measured in September 2019, 9 native forest plots, 8 oil palm plantation plots, 6 rubber plantation plots and 9 non-forest plots (shrubs and grass, 7 plots in dryland, and 2 plots in wetland). The plots were selected using a purposive sampling based on a minimum of 16 (4 × 4 pixels) visibly similar pixels of Landsat-8 images to represent one hectare of homogeneous land surface. A belt-shaped transect was used to collect information on vegetation characteristics, such as canopy cover, tree density, and ground cover [44], with a sample area of 500 m 2 (5 × 100 m), which was multiplied by 20 to represent an area of 1 hectare [45].
The canopy cover percentage was measured using a template developed by the USDA [46] and was assessed visually at two sampling locations for each plot. The ground vegetation cover was assessed visually using the step point method to estimate percentage (%) of live plant and bare soil cover [47]. The diameter at breast height (1.3 m) of all live trees with a diameter ≥5 cm was measured using a diameter tape. The tree height was measured on 3 to 4 representative trees using a Vertex III [48].

Methodological Framework
Decision trees have been widely used in land cover classification in remote sensing applications [49] and they are often used for integration of multi-source data such as optical and SAR images [23,50,51]. In this study, a decision tree was developed to discriminate Remote Sens. 2022, 14, 3 6 of 18 between non-forest, native forest, oil palm and rubber plantations to produce a land cover map for 2018.
The decision tree was applied by integrating forest and non-forest maps produced by Indonesia's LCCA program with indices derived from optical (Landsat-8) images and SAR (Sentinel-1) images to assist in identifying oil palm and rubber plantation areas. The decision tree rules were based on existing optical and SAR indices, specifically: (1) the difference of the two Sentinel-1 bands (VH and VV), (2) the Grey Level Co-occurrence Matrix (GLCM) texture, (3) spectral indices of Landsat-8 images, and (4) thresholds from a training sample validated with high-spatial resolution images of SPOT 6/7 and field measurements.
The mapping process to produce the land cover map for 2018 is shown in Figure 2. The decision tree classification was based on a series of binary decisions determined from thresholds applied to the SAR and optical indices. The LCCA forest/non-forest map was used as a root node of the decision tree, while oil palm and rubber plantations were discriminated at subsequent nodes. The steps to determine the indices and threshold for oil palm and rubber plantations discrimination are outlined in Sections 2.6 and 2.7. The decision tree used to produce the land cover map and the subsequent validation of this map using high resolution images are discussed in Section 3.4.

Methodological Framework
Decision trees have been widely used in land cover classification in remote sensing applications [49] and they are often used for integration of multi-source data such as optical and SAR images [23,50,51]. In this study, a decision tree was developed to discriminate between non-forest, native forest, oil palm and rubber plantations to produce a land cover map for 2018.
The decision tree was applied by integrating forest and non-forest maps produced by Indonesia's LCCA program with indices derived from optical (Landsat-8) images and SAR (Sentinel-1) images to assist in identifying oil palm and rubber plantation areas. The decision tree rules were based on existing optical and SAR indices, specifically: (1) the difference of the two Sentinel-1 bands (VH and VV), (2) the Grey Level Co-occurrence Matrix (GLCM) texture, (3) spectral indices of Landsat-8 images, and (4) thresholds from a training sample validated with high-spatial resolution images of SPOT 6/7 and field measurements.
The mapping process to produce the land cover map for 2018 is shown in Figure 2. The decision tree classification was based on a series of binary decisions determined from thresholds applied to the SAR and optical indices. The LCCA forest/non-forest map was used as a root node of the decision tree, while oil palm and rubber plantations were discriminated at subsequent nodes. The steps to determine the indices and threshold for oil palm and rubber plantations discrimination are outlined in Sections 2.6 and 2.7. The decision tree used to produce the land cover map and the subsequent validation of this map using high resolution images are discussed in Section 3.4.

Forest/Non-Forest Map from Landsat
The LCCA program defines a forest as a group of trees taller than 5 m with a canopy cover of more than 30%, as specified by Indonesia's Ministry of Forestry [10,16]. The LCCA's forest/non-forest map does not explicitly identify plantations, such as oil palm and rubber. These plantations can occur in both the forest and the non-forest classes [10].
The annual LCCA maps with forest/non-forest classes are used to derive a change product of forest loss and forest regrowth for each annual interval [10]. These maps use the succession of Landsat instruments as the primary data, which includes the Thematic Mapper (TM), the Enhanced Thematic Mapper Plus (ETM+), and the OLI sensors.
The production of the LCCA maps involves two main steps, as outlined in [10]. First, image pre-processing consists of manual scene selection, orthorectification, radiometric terrain correction, cloud masking, and mosaicking. Second, multitemporal land cover classification steps consist of developing a forest base from the training sample to determine forest indices and thresholds and application of a multitemporal processing step using a Canonical Probability Network [7,10]. The multitemporal process aims to combine and refine the entire time series of individual yearly classification maps and to fill the

Forest/Non-Forest Map from Landsat
The LCCA program defines a forest as a group of trees taller than 5 m with a canopy cover of more than 30%, as specified by Indonesia's Ministry of Forestry [10,16]. The LCCA's forest/non-forest map does not explicitly identify plantations, such as oil palm and rubber. These plantations can occur in both the forest and the non-forest classes [10].
The annual LCCA maps with forest/non-forest classes are used to derive a change product of forest loss and forest regrowth for each annual interval [10]. These maps use the succession of Landsat instruments as the primary data, which includes the Thematic Mapper (TM), the Enhanced Thematic Mapper Plus (ETM+), and the OLI sensors.
The production of the LCCA maps involves two main steps, as outlined in [10]. First, image pre-processing consists of manual scene selection, orthorectification, radiometric terrain correction, cloud masking, and mosaicking. Second, multitemporal land cover classification steps consist of developing a forest base from the training sample to determine forest indices and thresholds and application of a multitemporal processing step using a Canonical Probability Network [7,10]. The multitemporal process aims to combine and refine the entire time series of individual yearly classification maps and to fill the missing data due to cloud cover. The result of applying these steps is a forest/non-forest map for the entire country each year. The current study used a final forest/non-forest map for 2018 based on the Landsat-8 image data listed in Table 1.

Discrimination of Oil Palm Plantation from Native Forest
Previous studies showed that it is difficult to differentiate between native forest and oil palm plantations using Landsat images in the tropics [12,13]. Thus, the oil palm Remote Sens. 2022, 14, 3 7 of 18 discrimination indices and thresholds were developed from the single-date mosaic of Sentinel-1A SAR images using the training samples derived from ground truth plots and SPOT 6/7 images.
The best performing indices and thresholds were determined by assessing individual Sentinel-1 VH and VV bands, as well as several band combinations including the ratio of VV and VH (VV/VH), the difference of VH and VV (VH − VV), a normalized difference index (NDI) derived from (VV − VH)/(VV + VH), and a gray level co-occurrence matrix (GLCM) texture using reference data and feature separability [23,24].
The GLCM algorithm provides an ability to quantify texture based on mean, variance, homogeneity, and contrast within a pre-defined pixel window [52]. A previous study showed that GLCM variance reduced uncertainty and improved the ability to detect small regions of forest/plantation fragmentation and/or urban development [53,54]. To improve the accuracy of oil palm plantations discrimination, a GLCM variance based on standard deviation [55] with the input band of VH-VV and a window of 7 × 7 pixels was used.
The identification of relevant oil palm discrimination indices was based on the Jeffries-Matusita (JM) value estimated from the deviation between class means and its distance from the training samples. The JM value rates separability between two classes with output ranging from zero to two (two indicates the best feature separability) [56]. JM separability was computed using QGIS [57]. The plantation threshold was determined based on a 95% confidence interval of feature frequency histogram of oil palm samples, which was then validated visually.

Discrimination of Rubber Plantation from Native Forest
The vegetation indices often utilized for discrimination of rubber plantations are the Normalized Difference Vegetation Index (NDVI), the Normalized Difference Moisture Index (NDMI), and the Enhanced Vegetation Index (EVI) [58][59][60].
Mapping rubber plantations in the tropics is challenging because a rubber plantation has a similar structure to that of a native forest, and phenological phases may not be captured in optical images due to frequent cloud cover. Previous studies that have attempted to map the extent of rubber plantations using the EVI have shown misclassification with native forest and non-forest classes [61].
Here, we integrated Landsat-8 and Sentinel-1 imagery to improve the discrimination between rubber plantations and native forests. Several vegetation indices were evaluated to discriminate between rubber plantations and native forest. The C-band was used to discriminate between rubber and shrubs due to C-band penetration being more sensitive to sparse and low vegetation. The selected indices and thresholds for discrimination of rubber plantations were developed based on similar steps as used for oil palm indices. The vegetation indices that provided the greatest discrimination performance based on the JM values were NDVI and NDMI, as defined below: where RED is the red spectral band (0.64-0.67 µm) of the OLI, NIR is the near infrared band (0.85-0.88 µm) and SWIR 1 is the first OLI short-wave infrared band (1.57-1.65 µm).

Accuracy Assessment
A confusion matrix showing overall accuracy, producer's accuracy, and user's accuracy was used to assess the land cover map based on 784 validation points. The sample size was determined based on stratified random sampling, as outlined in [62], distributed proportionally to native forest, oil palm and rubber plantations, and non-forest classes. The minimum sampling unit was created from a block of 2 × 2 pixels in the land cover change map (50 × 50 m), which represents a minimum forested area of 0.25 hectares [16]. The land cover at each of the validation points was assessed with high-resolution images from SPOT 6/7 recorded in 2018 and accessed using Google Earth.

Field Results
Field measurement showed great variability of vegetation structure between land cover types. Native forest had the highest tree density and highest canopy cover out of the four classes (Figure 3a,b). Trees in native forests had an average diameter of 19 cm (Figure 3d), yet some trees had a diameter of up to 80 cm and almost no bare soil present (Figure 3e).
Oil palm plantations were generally planted in a symmetrical pattern on relatively flat ground with tree density of 160 to 180 trees/ha (Figure 3b), heights of 12 to 17 m (Figure 3c) and diameters ranging from 50 to 80 cm (Figure 3d). Rubber plantations had tree densities around 410 trees/ha (Figure 3b), diameters ranging from 10 to 26 cm (Figure 3d), and heights of 13 to 23 m (Figure 3c). The average ground cover on rubber plantations was 60% (Figure 3e).
Shrub and grassland vegetation, on average, were less than 0.84 m tall, with some shrubs reaching 2 m. Few trees were present in this land class (Figure 1). Field observations suggested that oil palm and rubber plantations may be discriminated based on characteristics of semi-open canopy, gaps between the trees, height, and ground cover vegetation (Figure 3a-d).

Algorithm for Mapping Oil Palm Plantation
The separability of oil palm plantations, native forests and non-forest (shrub/grass) in Sentinel-1 images was analyzed based on the frequency distribution histogram of training samples. Oil palm plantations had a VH backscatter with values between −16 and −12 dB (Figure 4a), while the VV ranged from −8 to −4 dB (Figure 4b). Oil palm plantations also produced the largest VH and VV difference (>7.5 dB) compared to rubber plantations and native forest. However, this index did not clearly discriminate between native forests and rubber plantations (Figure 4e). Oil palm plantations were not easily discriminated using the ratio VV/VH and NDI, as shown in Figure 4c,d.
The GLCM variance was applied to improve the ability to detect oil palm plantations and remove noise in the SAR images. The GLCM variance further increased separability between oil palm plantations and the other classes (Figure 4f). This was combined with the oil palm plantation threshold of VH backscatter to separate oil palm plantations from native forest and non-forest (Figure 4a). where RED is the red spectral band (0.64-0.67 µm) of the OLI, NIR is the near infrared band (0.85-0.88 µm) and SWIR 1 is the first OLI short-wave infrared band (1.57-1.65 µm).

Accuracy Assessment
A confusion matrix showing overall accuracy, producer's accuracy, and user's accuracy was used to assess the land cover map based on 784 validation points. The sample size was determined based on stratified random sampling, as outlined in [62], distributed proportionally to native forest, oil palm and rubber plantations, and non-forest classes. The minimum sampling unit was created from a block of 2 × 2 pixels in the land cover change map (50 × 50 m), which represents a minimum forested area of 0.25 hectares [16]. The land cover at each of the validation points was assessed with high-resolution images from SPOT 6/7 recorded in 2018 and accessed using Google Earth.

Field Results
Field measurement showed great variability of vegetation structure between land cover types. Native forest had the highest tree density and highest canopy cover out of the four classes (Figure 3a,b). Trees in native forests had an average diameter of 19 cm (Figure 3d), yet some trees had a diameter of up to 80 cm and almost no bare soil present (Figure 3e).
Oil palm plantations were generally planted in a symmetrical pattern on relatively flat ground with tree density of 160 to 180 trees/ha (Figure 3b

Algorithm for Mapping Oil Palm Plantation
The separability of oil palm plantations, native forests and non-forest (shrub/grass) in Sentinel-1 images was analyzed based on the frequency distribution histogram of training samples. Oil palm plantations had a VH backscatter with values between −16 and −12 dB (Figure 4a), while the VV ranged from −8 to −4 dB (Figure 4b). Oil palm plantations also produced the largest VH and VV difference (>7.5 dB) compared to rubber plantations and native forest. However, this index did not clearly discriminate between native forests and rubber plantations (Figure 4e). Oil palm plantations were not easily discriminated using the ratio VV/VH and NDI, as shown in Figure 4c,d.
The GLCM variance was applied to improve the ability to detect oil palm plantations and remove noise in the SAR images. The GLCM variance further increased separability between oil palm plantations and the other classes (Figure 4f). This was combined with the oil palm plantation threshold of VH backscatter to separate oil palm plantations from native forest and non-forest (Figure 4a). Evaluation of class separability was conducted using the JM statistic (Table 2). This showed that, from the training sample, the combination of the GLCM variance and the VH backscatter resulted in close to complete separability between oil palm plantations and native forest (1.978) and complete separability between oil palm and rubber plantations (2.0). However, these combined indices resulted in low separability between native forest and rubber plantations. Thus, the final indices and thresholds used to detect oil palm plantations were a GLCM variance of less than 40 and a VH backscatter intensity between −15 and −13 dB.

Algorithm for Mapping Rubber Plantation
Rubber plantations were discriminated from other classes using NDVI and NDMI. Rubber plantations had relatively high NDVI with values of 0.60-0.80. However, this NDVI range overlaps with values for native forest (0.65-0.85) and shrubs within the non-forest class (0.60-0.80) (Figure 5a). The discrimination between rubber plantations and native forest was clear when using NDMI (Figure 5b), although there was a significant NDMI overlap between rubber and shrub/grass in the dryland and wetland regions of the study area. Thus, the combination of VV and VH backscatter of rubber plantations was applied to discriminate the overlap between native forest, rubber plantation and shrub/grassland.  Evaluation of class separability was conducted using the JM statistic (Table 2). This showed that, from the training sample, the combination of the GLCM variance and the VH backscatter resulted in close to complete separability between oil palm plantations and native forest (1.978) and complete separability between oil palm and rubber plantations (2.0). However, these combined indices resulted in low separability between native forest and rubber plantations. Thus, the final indices and thresholds used to detect oil shrub/grassland. Rubber plantations had relatively higher VH (Figure 4a) and lower VV (Figure 4b) compared to oil palm plantations, while the rubber plantation VH and VV backscatter values overlapped with native forests. However, the rubber plantation is relatively well distinguished from shrub and grasslands within the non-forest class based on higher VH and VV backscatter (Figure 4a,b). Evaluation of rubber plantation indices based on the JM value showed that the use of NDVI and NDMI alone provided less accurate discrimination than when combined with VV and VH backscatter ( Table 3). The JM value for discrimination of native forest and rubber plantations was 1.993, while the JM values for discrimination of rubber and non-forest in dryland and wetland were 1.977 and 1.721, respectively. Thus, the final indices and thresholds to discriminate the rubber plantations were an NDVI between 0.67 and 0.77, an NDMI of less than 0.2, a VV backscatter between −8.5 and −6 dB, and a VH backscatter from −14 to −12 dB. Rubber plantations had relatively higher VH (Figure 4a) and lower VV (Figure 4b) compared to oil palm plantations, while the rubber plantation VH and VV backscatter values overlapped with native forests. However, the rubber plantation is relatively well distinguished from shrub and grasslands within the non-forest class based on higher VH and VV backscatter (Figure 4a,b).
Evaluation of rubber plantation indices based on the JM value showed that the use of NDVI and NDMI alone provided less accurate discrimination than when combined with VV and VH backscatter ( Table 3). The JM value for discrimination of native forest and rubber plantations was 1.993, while the JM values for discrimination of rubber and non-forest in dryland and wetland were 1.977 and 1.721, respectively. Thus, the final indices and thresholds to discriminate the rubber plantations were an NDVI between 0.67 and 0.77, an NDMI of less than 0.2, a VV backscatter between −8.5 and −6 dB, and a VH backscatter from −14 to −12 dB. Table 3. JM value for rubber plantation discrimination from native forest and shrubs.

Land Cover Map and Validation
The decision tree classification for discrimination of oil palm and rubber plantations to produce land cover map for 2018 is shown in Figure 6. This final land cover map comprised four classes: oil palm plantation, rubber plantation, native forest, and non-forest, as shown in Figure 7.
Spatial application of the decision tree showed that oil palm plantations were mainly located in the east of the study area (the region of Paser). Large-scale oil palm plantations were also detected in the central region of Tabalong (shown in the zoom box of Figure 7). This was validated with high spatial resolution SPOT 6 images acquired on 16 August 2018. Small areas of oil palm plantation were also scattered in the west of the region of Barito Timur.
Rubber plantations were mainly observed in the middle of the Tabalong and Balangan regions in relatively flat areas. Rubber plantations are also suitable to grow in relatively higher elevations and have been detected in the northwest of the Paser regency. The map also shows that rubber plantations have been established on a relatively small scale and are scattered within remnants of native forests.

Land Cover Map and Validation
The decision tree classification for discrimination of oil palm and rubber plantations to produce land cover map for 2018 is shown in Figure 6. This final land cover map comprised four classes: oil palm plantation, rubber plantation, native forest, and non-forest, as shown in Figure 7. Spatial application of the decision tree showed that oil palm plantations were mainly located in the east of the study area (the region of Paser). Large-scale oil palm plantations were also detected in the central region of Tabalong (shown in the zoom box of Figure 7). This was validated with high spatial resolution SPOT 6 images acquired on 16 August 2018. Small areas of oil palm plantation were also scattered in the west of the region of Barito Timur.
Rubber plantations were mainly observed in the middle of the Tabalong and Balangan regions in relatively flat areas. Rubber plantations are also suitable to grow in relatively higher elevations and have been detected in the northwest of the Paser regency. The map also shows that rubber plantations have been established on a relatively small scale and are scattered within remnants of native forests.
Validation of the land cover map was assessed using a confusion matrix ( Table 4). The overall accuracy was 92%, with producer's accuracy for all classes higher than 90%, except for rubber, which was 65%. The native forest class had a producer's accuracy of 97% and a user's accuracy of 93% and showed good discrimination of native forest from other classes, particularly oil palm and rubber plantations. The commission error of 7% was largely attributable to misclassification of rubber plantations.
The oil palm plantation had a higher producer's accuracy of 91% and a user's accuracy of 84% compared to the rubber plantation. The rubber plantation class had the lowest accuracy, indicated by a 65% producer's accuracy and an 80% user's accuracy. This was Figure 6. The decision tree to discriminate between native forest, oil palm plantation, rubber plantation, and non-forest. identified from the omission error of 35% for the rubber plantation and is likely due to difficulty discriminating between native forests with lower tree density and rubber plantations.   Validation of the land cover map was assessed using a confusion matrix ( Table 4). The overall accuracy was 92%, with producer's accuracy for all classes higher than 90%, except for rubber, which was 65%. The native forest class had a producer's accuracy of 97% and a user's accuracy of 93% and showed good discrimination of native forest from other classes, particularly oil palm and rubber plantations. The commission error of 7% was largely attributable to misclassification of rubber plantations. The oil palm plantation had a higher producer's accuracy of 91% and a user's accuracy of 84% compared to the rubber plantation. The rubber plantation class had the lowest accuracy, indicated by a 65% producer's accuracy and an 80% user's accuracy. This was identified from the omission error of 35% for the rubber plantation and is likely due to difficulty discriminating between native forests with lower tree density and rubber plantations.

Discussion
This study demonstrated the use of multi-source indices derived from optical and SAR images for land cover mapping. The study resulted in improvement in land cover classification accuracy relative to existing Indonesian forest maps (LCCA). Specifically, the method developed improved the ability to discriminate between native forest and oil palm and rubber plantations. Previous studies showed that the LCCA forest and nonforest map, which uses optical Landsat images, had lower overall accuracy (73-77%) due to misclassification of plantations [31]. The integration of Landsat-8 OLI and Sentinel-1 C-band increased land cover classification accuracy, resulting in an overall accuracy of 92%, with an increase in accuracy for the four classes (native forest, oil palm plantation, rubber plantation, and non-forest). Moreover, the multi-source approach resulted in producer's and user's accuracy for most land cover classes of more than 85% [63] (Table 4), which fulfills Indonesia's desired criteria for map accuracy [64]. The exception was the accuracy for the rubber plantation class, which had a producer's accuracy of 65% and a user's accuracy of 80%.
The native forest class had the highest producer's accuracy (97%) and user's accuracy (93%), while the non-forest class was the next most accurate with a producer's accuracy of 94% and a user's accuracy of 96%. This high classification accuracy of native forest and non-forest classes was a direct result of explicitly separating oil palm and rubber plantation classes. Although oil palm plantations had an omission error of 9% and a commission error of 16%, the oil palm class had a high producer's accuracy of 91% and user's accuracy of 84%.
The high classification accuracy for oil palm plantations showed the benefit of integrating C-band data from the Sentinel-1 images, particularly in discriminating other closedcanopy forests (cover > 50%) such as native forests and rubber plantations (Figure 3a). This is in contrast to optical images, which have shown limitations in discriminating vegetation with dense canopy cover.
In agreement with Miettinen et al. [23], the Sentinel-1 C-band could discriminate oil palm plantations based on the differences between VV and VH backscatter. This is likely due to the vegetation structure of oil palm plantations, which differs from that of rubber plantations and native forests. Oil palm plantations have a high VV backscatter due to canopy scattering from a large crown and fronds, while VH backscatter was low compared to other vegetation classes due to the low density of stems (Figure 3a,b).
The lower accuracy for mapping rubber plantations could be due to dependency on discriminating rubber plantations based on C-band SAR. The current study used optical Landsat-8 images of less than 50% cloud cover in combination with the relatively short wavelength of the Sentinel-1 C-band. This wavelength is less likely to penetrate through the canopy and differentiate between the distinct mid-and understory of native forest and rubber plantations. Use of longer wavelengths with greater penetration, such as the L-band (wavelength~23.5 cm) from ALOS PALSAR, would result in greater penetration through the canopy and more distinct scattering from trunks and branches, twigs, and ground cover vegetation [65].
The other causes of the low accuracy of rubber plantation discrimination may be the calibration of the SAR digital numbers using sigma nought. The sigma nought may produce significant variation of the backscatter coefficients depending on the incidence angle, wavelength, polarization, topography and land surface [66]. Therefore, previous studies showed that image calibration using gamma nought results in improvements in classification, particularly for mapping forest types based on variations in structure across varying topography [19,67].
The current study indicates that rubber plantations have an NDVI ranging from 0.67 to 0.77, which is lower than the NDVI expected for native forests (<0.80) (Figure 5a). This is in agreement with previous studies, where time-series images of rubber plantation NDVI showed values of less than 0.80, while native forests, including temperate and tropical forests, had an NDVI of up to 0.9 [68,69].
NDMI is sensitive to the moisture content of leaves. An NDMI of greater than 0.28 indicates healthy vegetation with dense canopy cover, while bare soil has low and negative values of the index [60,70]. For example, native evergreen forest has been shown to produce an NDMI of more than 0.40 all year [71]. The low value of NDMI for rubber plantations may be caused by low stem density and semi-open canopy, leading to drying of the soil [72]. For this reason, rubber plantations were discriminated from native forest using a threshold of less than 0.2. However, the NDMI threshold works as a function of vegetation area fraction that does not significantly differ for wet or dry soil areas [73].
In tropical regions, shrubs are often misclassified as plantation and forest regrowth, particularly in their early years [74][75][76]. Thus, shrubs and grassland had an overlap with rubber plantations with respect to their NDVI and NDMI distributions ( Figure 5). The shorter C-band wavelength is sensitive to shrub, grass, and low vegetation structure and showed a relatively low backscatter value for both VV and VH polarization compared to rubber plantations [66]. Therefore, backscatter of VV and VH was applied to discriminate shrubs from rubber plantations, which increased the rubber class separability quantified by the JM value (Table 3).
Our study resulted in the generation of a land cover map in 2018. However, this should not necessarily be seen as a final product. It is important to realize that spatio-temporal processing steps are required to improve the mapping accuracy and temporal consistency of the land cover classifications by reducing unlikely transitions between classes. For example, implementation of the multitemporal processing method using a Bayesian Hidden Markov model can be used to combine continuous and discrete variables to improve classification accuracy [77,78].

Conclusions
Accurate and timely information of land cover and land cover change is required to help countries in the tropics to comply with the IPCC principles of transparency, comparability, consistency, completeness, and accuracy for sustainable forest management and investment decisions. Indonesia's national forest monitoring system currently relies on mapping two land cover classes of forest and non-forest. This study developed multi-source indices to discriminate between native forests and tree plantations, such as oil palm and rubber, by integrating optical Landsat-8 and Sentinel-1 SAR images. The use of SAR images helped to improve the accuracy of forest mapping, which is currently strongly impacted by frequent cloud cover. Integration of multi-source images provides more comprehensive information on which to base land cover classification. Landsat-8 images can differentiate vegetation types such as rubber plantations from spectral indices such as NDVI and NDMI, while Sentinel-1 C-band SAR can be used to differentiate between native forest and tree plantations based on backscatter characteristics.
The current study shows integration of optical and SAR images improved classification with an overall accuracy of 92%. All four of the land cover classes fulfilled the minimum target accuracy, except for rubber plantations. Therefore, improvement of the current method using applications of multitemporal processing will be required in order to develop a final product. Once this has been achieved, this method can be used to extend Indonesia's national mapping to four land cover classes of native forest, oil palm plantation, rubber plantation, and non-forest.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.