Using Bayesian multitemporal classification to monitor tropical forest cover changes in Kalimantan, Indonesia

ABSTRACT

Significant areas of native forest in Kalimantan, on the island of Borneo, have been cleared for the expansion of plantations of oil palm and rubber. In this study multisource remote sensing was used to develop a time series of land cover maps that distinguish native forest from plantations. Using a study area in east Kalimantan, Landsat images were combined with either ALOS PALSAR or Sentinel-1 images to map four land cover classes (native forest, oil palm plantation, rubber plantation, non-forest). Bayesian multitemporal classification was applied to increase map accuracy and maps were validated using a confusion matrix; final map overall accuracy was >90%. Over 18 years from 2000 to 2018 nearly half the native forests in the study area were converted to either non-forest or plantations of either rubber or oil palm, with the highest losses between 2015 and 2016. Trending upwards from 2008 large areas of degraded or cleared forests, mapped as non-forest, were converted to oil palm plantation. Conversion of native forests to plantation mainly occurred in lowland and wetland forest, while significant forest regrowth was detected in degraded peatland. These maps will help Indonesia with strategies and policies for balancing economic growth and conservation.

KEYWORDS:

• Multitemporal classification
• Bayesian networks
• Landsat
• ALOS PALSAR
• sentinel-1
• Indonesia

1. Introduction

Tropical forests cover approximately 45% of world forested area, account for about 60% of the global terrestrial carbon pool and play a vital role in biodiversity conservation (FAO 2020a; Poker and MacDicken 2016). Between 1990 and 2015, the area of global tropical forests declined by 10%, from 1966 million hectares to 1770 million hectares (Keenan et al. 2015). The countries with the largest decreases of native tropical forest area were Brazil, Indonesia, and the Democratic Republic of the Congo (Turubanova et al. 2018). In Indonesia, from 2000 to 2012, native forest area decreased by approximately 6.02 million hectares (Margono et al. 2014).

Remote sensing has been widely used to monitor global forest changes (Achard and Hansen 2012). The Landsat series of satellites have been used for global forest change monitoring using optical imagery for more than 40 years (Hansen et al. 2013). However, forest monitoring programs that have largely relied on single satellite sensors are limited in their ability to distinguish between different land cover types and to track changes in their extent. For example under certain conditions, land cover mapped using optical Landsat has limited ability to distinguish native forests, plantations and perennial crops (Curtis et al. 2018; Tropek et al. 2014). These limitations are generally due to similar spectral characteristics of these land cover types (Torbick et al. 2016).

Integration of multisource remote sensing data can provide additional information on surface reflectance characteristics and structural properties of land cover types. Moreover, integration of multisource remote sensing data increases temporal data density for capturing land cover in cloudy regions, characterizing seasonal variations, and in reducing uncertainty caused by variations in illumination conditions. This includes incorporating radar data, which can be used to distinguish variations in canopy structure and moisture characteristics (Bannari et al. 1995; Hoekman 1985). Furthermore, the use of radar images is essential to increase the temporal frequency of available imagery in cloudy tropical regions, further improving classification accuracy. Satellite-based synthetic aperture radar (SAR) images, such as Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR), have been used to map global forest/non-forest areas (Shimada et al. 2014).

Free access to Global 25 m resolution PALSAR-2/PALSAR mosaic and forest/non-forest map and Sentinel-1 increases the data available for the system to support monitoring land cover changes (Altunel, Akturk, and Altunel 2020). Moreover, high-resolution SAR images such as the C band of Environmental Research Satellite (ERS) and Sentinel-1 are superior in agricultural classification (Blaes, Vanhalle, and Defourny 2005; Ban 2003; Miettinen, Liew, and Kwoh 2015). At the same time, free access to Sentinel-1 (10 m spatial resolution) since 2014 has also helped monitor global tree plantation development (Fagan et al. 2022).

Sarzynski et al. (2020) integrated Landsat and ALOS-2/PALSAR-2 images for the tropical forests of Sumatra to improve discrimination between forests and oil palm plantations from 74% to 84%. A time series of ALOS-1/2 PALSAR-1/2 and Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) was used to map oil palm plantations in Indonesia and Malaysia from 2001 to 2016, where a disturbance and recovery algorithm (Breaks For Additive Season and Trend – BFAST) of MODIS was applied to detect land cover change time points for the years without PALSAR data (2011–2014 and 2001–2006) (Xu et al. 2020; Verbesselt et al. 2010). ALOS PALSAR and MODIS images were also integrated to map rubber plantations in China for 2000–2010 based on a phenology-based method (Kou et al. 2018). However, oil palm and rubber plantations in the equatorial regions present as evergreen trees and do not exhibit a distinct seasonal variation in solar spectral reflectance characteristics (Cracknell et al. 2013).

A combination of multisource data such as L- and C-bands SAR images, Landsat and MODIS optical images increases near-real-time data for forest monitoring and plantation dynamics (Reiche et al. 2018; Dong et al. 2015). The integration of Landsat and Sentinel-1 imagery from 2018 was used to map the two largest in area plantation types (oil palm and rubber) in Kalimantan (Sari et al. 2022; Santosa, Sunkar, and Kwatrina 2020). Thus, integrating multisource and multitemporal remote sensing images can benefit the mapping of land cover dynamics in tropical regions for monitoring forest changes, including oil palm and rubber plantation expansion.

Indonesia’s current national forest monitoring process is based on visual interpretation (VI) of Landsat imagery, supplemented by automated digital processing (ADP) (Margono, Usman, and Sugardiman 2016; LAPAN 2014). However, technical improvements of the VI map are required due to the time-consuming and subjective nature of the process (Margono, Usman, and Sugardiman 2016). While improvements in the accuracy of the ADP maps is necessary, extension of this automated process to discriminate more land cover types (particularly identifying where native forest is being converted to oil palm and rubber plantations) is also critical (LAPAN 2014).

This study aims to improve the accuracy of Indonesia’s ADP mapping for better monitoring of native forest loss, particularly due to plantation expansion. The current study applied multitemporal processing using a Bayesian network to integrate Landsat, ALOS PALSAR, and Sentinel-1 to improve land cover maps that differentiate between native forest, non-forest, oil palm and rubber plantations for 2000–2018.

2. Materials and methods

2.1. Study area

The study area of more than 1.9 million hectares is in the east of Kalimantan on the island of Borneo (Figure 1). The terrain in this region is variable, with elevation ranging from sea level up to about 1550 m. The rainfall is generally more than 2000 mm. Native forests, oil palm plantations, rubber plantations, and non-forested areas such as settlements, rice fields, mixed agriculture, shrubs and grasslands, and mining, comprise the land cover in the region. Oil palm and rubber plantations have increased in this region of Indonesia since the 1990s (Ministry of Agriculture 2017a, 2017b). For example, from 2009 to 2018, oil palm and rubber plantations increased by 60% from 244,899 hectares to 391,286 hectares (Badan Pusat Statistik 2010a, 2010b, 2010c, 2010d, 2019a, 2019b, 2019c).

Figure 1. Location of the study area.

2.2. Data and image pre-processing

The current study used a time series multisource of optical and SAR images as follows: (1) from 2007 to 2010 using Landsat 5, Landsat 7, and ALOS PALSAR, (2) from 2015 to 2018 using Landsat 8 and Sentinel-1. These data were used to produce individual year land cover maps. All Landsat images from 2000 to 2018 were obtained from the Land Cover Change Analysis (LCCA) project (LAPAN 2014), and provided by the United States Geological Survey (USGS) (2019). Sentinel-1 images from 2015 to 2018 were downloaded from the Alaska Satellite Facility (ASF) (2019), and ALOS PALSAR images from 2007 to 2010 were downloaded from the Japan Aerospace Exploration Agency (JAXA) (2019). The sensor characteristics for each instrument, pre-processing, and image details are described in the following sections.

2.2.1. Landsat

The Landsat series of satellites have imagery from several earth-observing instruments available. Specifically, we utilized Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) sensor imageries. These sensors all include solar reflectance spectral bands ranging from visible, thermal, and shortwave infrared (SWIR), including panchromatic bands. The pre-processing steps to produce annual mosaicked Landsat images at 25 m resolution resampled with six multispectral image bands (blue, green, red, near infrared, and shortwave infrared 1 and 2) were based on the LCCA method outlined in (LAPAN 2014), and include scene selection; geometric, radiometric, and terrain correction; cloud masking; and mosaicking of the individual Landsat images.

2.2.2. Sentinel-1

The study used Sentinel-1A C-band (wavelength ∼5.6 cm) SAR imagery with the dual-polarization of VV (vertical transmit, vertical receive) and VH (vertical transmit, horizontal receive). Sentinel-1A satellites can acquire data for the same ground location every 12 days. The SAR mosaic images in the current study developed from annual mosaic images, including annual mosaic images of Global 25 m ALOS PALSAR and Sentinel-1 to differentiate between dense cover trees of native forest, oil palm and rubber plantation in the tropics based on their vegetation structure. Three scenes of Sentinel-1A Interferometric Wide (IW) swath mode were acquired between January and February (wet season) in Level-1 Ground Range Detected (GRD) format in descending mode to produce annual mosaiced images. The pre-processing steps used the Sentinel Application Platform (SNAP) open-source software version 7.0.4 to produce a 25 m resampled image (to match Landsat images resampled in the LCCA) and to derive the backscatter intensity to gamma nought, in units of decibel (dB). The steps consisted of applying precise orbit acquisition information, radiometric calibration, and terrain correction.

2.2.3. Alos PALSAR

Global 25 m resolution of annual mosaics of ALOS PALSAR L-band (wavelength ∼23.5 cm) provided by JAXA were used for the study. The mosaic images were SAR Analysis Ready Data (ARD), meaning that the images had geometric, terrain and radiometric corrections applied. The ALOS PALSAR data included dual polarizations of HH (horizontal transmit, horizontal receive) and HV (horizontal transmit, vertical receive). The backscatter intensity for each of these bands was provided in gamma nought as outlined in (JAXA 2016).

2.3. Methodological overview

Incorporating multitemporal information into land cover classification can improve accuracy compared to maps produced using only a single-date (Symeonakis et al. 2012). This is because multiple images can be used to build additional confidence in class attribution. Multitemporal classification can also be used to reduce misclassification issues associated with missing data due to cloud cover, variable atmospheric conditions, and sensor and illumination angle effects (LAPAN 2014; Caccetta et al. 2016).

Bayesian networks use conditional probability to determine the likely outcome given a set of independent variables (Pollino and Henderson 2010; Ghahramani 2001). Previous studies have demonstrated Bayesian multitemporal classification to integrate Landsat and ALOS PALSAR from 2005 to 2010 (Reiche et al. 2015). These Bayesian networks can estimate and track land cover change, with probability conditioned on the basis of previous land cover state (Lee, Cardille, and Coe 2020).

The current study applies a multitemporal classification using a Bayesian network for integrating Landsat, ALOS PALSAR, and Sentinel-1 data. The classification probability for each pixel at each date in the time series is based on the Bayesian network calculation and reflects each pixel’s relative certainty for the final land cover class designation.

2.4. Individual year classification using decision trees

A method to produce land cover maps for Kalimantan for individual years, with no multi-temporal dependency, using a decision tree based on indices and image thresholds, has been outlined in our previous study (Sari et al. 2022). In brief, a decision tree was used to combine the forest and non-forest maps produced by Indonesia’s LCCA program with Landsat, ALOS PALSAR, and Sentinel-1 images, to produce more accurate land cover classifications that incorporate four classes (native forest, oil palm plantation, rubber plantation, and non-forest).

However, during 2000–2006 and 2011–2014, there were no ALOS PALSAR or Sentinel-1 data available covering the study area (CEOS 2014). To fill this gap in SAR data, we developed a method to use the shortwave infrared (SWIR 1: $\lambda \cong 1.55\mu \mathrm{m}\mathrm{t}\mathrm{o}1.75\mu \mathrm{m}$) band from Landsat images to discriminate between native forest and other land cover classes. Previous studies have shown that SWIR 1 can be used to discriminate between vegetation types based on plant canopy water status (Olsen et al. 2013; Tucker 1980). In this study, we adapted the previously published method, by replacing decision tree nodes that rely on SAR data, with indices and thresholds based on SWIR 1, NDVI and Normalized Difference Moisture Index (NDMI) (Tucker 1980; Razak et al. 2017; Gao 1996; Sari et al. 2022).

The decision tree applied to Landsat images from 2000 to 2006, and 2011 to 2014, to discriminate between native forest and tree plantations, is shown in Figure 2.

Figure 2. The individual year classification decision tree to discriminate between native forest, oil palm, rubber plantations, and non-forest developed from Landsat images, for years when no SAR data were available; modified from Sari et al. (2022).

The SWIR 1 threshold was determined based on a training sample of Landsat images in 2018 with four land cover classes of oil palm (7500 pixels), rubber plantations (8000 pixels), native forest (7100 pixels), and non-forest (12,500 pixels distributed across shrub and grass in dryland and wetland). The SWIR 1 threshold was based on a 95% confidence interval of feature frequency histogram of training samples, which was then validated visually. The final SWIR 1 thresholds for forest cover were between 0.13 and 0.23 (Figure 3).

Figure 3. Mean reflectance of training sample for native forest, oil palm and rubber plantations, and non-forest (dryland and wetland) in Landsat bands with band 1 (Blue), band 2 (Red), band 3 (Green), band 4 (Near Infrared), band 5 (SWIR 1), and band 6 (SWIR 2).

The thresholds for NDVI and NDMI used for discrimination of native forest, oil palm and rubber plantation, were derived in a similar way to the SWIR 1 thresholds. The methodology used is outlined in (Sari et al. 2022).

The plantation occurrence images were employed in the multitemporal classification and were used to create boundaries for potential plantation areas. The plantation occurrence was developed from a linear combination of NDVI and SWIR 1 for oil palm, and NDMI and SWIR 1 for rubber, and processed with a median filter of three-by-three pixels.

2.6. Multitemporal classification

Multitemporal classification was used to refine the yearly four-class maps to produce the final time series of annual land cover maps. A Bayesian network was applied in reverse time order, where the classification for 2018 was used as an initial base classification image. Preceding classifications were then processed using Bayesian probability, based on the dependence of prior classification and the individual year probability of oil palm and rubber plantation occurrence.

The multitemporal algorithm is shown in Equations (1) and (2) and is adapted from Pollino and Henderson (2010). The graphical representation of the processing is shown in Figure 4. The multitemporal classification steps were:

1. determine the probability for each pixel to change class based on the two years of map classifications (e.g. 2017 and 2018) – designated as M. This probability P(M), was based on a transition matrix of 16 possible class changes, derived from the two individual land cover maps;

2. use Equation (1) to calculate pixel transition probability between forest and non-forest classes;

3. if the probability $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}}\right)$ of no-change is greater than the probability of change then the class will be stable. If the probability of change is greater, then the class for that pixel will change;

4. use Equation (2) for pixels that have transitions that involve oil palm and rubber plantation classes as either the initial or final classes to determine if the transition should occur, i.e. $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}},{\mathrm{R}}_{\mathrm{t}-1}\right)>$0.5. (1) $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}}\right)={\displaystyle \frac{\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}}|{\mathrm{M}}_{\mathrm{t}-1}\right)\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}\right)}{\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}}\right)}}$(1)

where $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}}\right)$ is the probability of the class determined by the decision tree in the preceding year (referred as ${\mathrm{M}}_{\mathrm{t}-1}$) given the posterior state from the decision tree in the following year $\left({\mathrm{M}}_{\mathrm{t}}\right)$. (2) $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}},{\mathrm{R}}_{\mathrm{t}-1}\right)={\displaystyle \frac{\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}},{\mathrm{R}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}-1}\right)\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}\right)}{\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}},{\mathrm{R}}_{\mathrm{t}-1}\right)}}$(2) where $\mathrm{P}\left({\mathrm{M}}_{\mathrm{t}-1}|{\mathrm{M}}_{\mathrm{t}},{\mathrm{R}}_{\mathrm{t}-1}\right)$ is the probability of the new class determined by the decision tree (referred as ${\mathrm{M}}_{\mathrm{t}-1}$) given the posterior state from the decision tree in 2018$\left({\mathrm{M}}_{\mathrm{t}}\right)$, and the individual year rubber plantation probability images for 2017, referred to as $\left({\mathrm{R}}_{\mathrm{t}-1}\right)$ Note that if the class assigned by the decision tree for 2017 is oil palm, then the probability of rubber plantation occurrence (${\mathrm{R}}_{\mathrm{t}-1}$) would be replaced by the probability for oil palm occurrence (${\mathrm{O}}_{\mathrm{t}-1}$).

Figure 4. Sample of graphical concept for the multitemporal process for 2018 and 2017, where M = individual year classification, O = approximate oil palm plantation occurrence, and R = approximate rubber plantation occurrence.

After the time series of classified maps is processed back to 2000, the next step is a second iteration of the same process in the forward direction from 2000 to 2018. These two iterations were used to reduce uncertainty of the final land cover classes (Guoping and Ren 2018). The final step was to reduce the occurrence of single pixels in the oil palm and rubber plantation classes based on a median filter of three-by-three pixels.

2.7. Map validation

To assess the accuracy of the temporally adjusted land cover maps over the time series, three land cover maps nine years apart were assessed: 2000, 2009, and 2018. These years represented the start, middle and end of the time series and for which high resolution images were available. We acknowledge that our approach contrasts with studies that evaluate time series classification accuracy in each year maps, such as from 2001 to 2013 (Tsutsumida and Comber 2015). However, there are also several studies on time series land cover classification that have used similar approach to temporally discrete validation (Chen et al. 2019; Yan et al. 2019). The maps were assessed using a confusion matrix. Overall accuracy (OA), producer’s accuracy (PA), and user’s accuracy (UA) were calculated based on 784 validation points for each map. The samples were selected based on stratified random sampling, as outlined in (Olofsson et al. 2014), and distributed proportionally across native forest, oil palm, 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 (Ministry of Environment and Forestry 2018). The land cover from the three maps at each of the validation datasets was assessed based on high-resolution images from SPOT 6/7 recorded in 2018 and temporal high-resolution images from the open-source platform for 2000 and 2009 (FAO 2020b).

3. Results

3.1. Classification results and accuracy assessment

Three land cover maps (2000, 2009, and 2018) were used for quantifying the accuracy of the native forest, oil palm, rubber plantation, and non-forest classes, relative to SPOT 6/7 and high-resolution validation data. Maps from these three years are shown in Figure 5.

Figure 5. Maps of land cover in the study area for 2000, 2009 and 2018.

Accuracy assessment of the three representative maps can be seen in Figure 6. Overall accuracy of the maps was above 90%. The native forest class had user’s and producer’s accuracy for all years of greater than 90%. The non-forest class had user’s and producer’s accuracy of more than 85% for all years. However, producer’s accuracy for rubber plantations was as low as 77%, and as low as 61% for oil palm in 2000.

Figure 6. Map accuracy assessment: (a) overall accuracy for the three validation years selected, (b) user’s accuracy and producer’s accuracy for the 2000 map, (c) user’s accuracy and producer’s accuracy for the 2009 map, and (d) user’s accuracy and producer’s accuracy for the 2018 map.

Figure 7 shows the estimated areas of each land cover class over the time series from 2000 to 2018. In 2000 the native forest covered 72.58% (1,390,080.69 hectares) of the total area, while non-forests covered 19.57% (376,286.25 hectares). Rubber plantations covered 7.21% (138,637.81 hectares) and oil palm plantations covered 0.95% (18,207.25 hectares) of the total area. Over the 18 years of the time series, native forest area gradually decreased, mainly due to an increase in non-forest area, followed by rubber and oil palm plantation expansion. By 2018, native forest area accounted for only 39.52% (760,055.57 hectares) of the total area. At the same time, non-forest areas increased to 37.35% (718,387.56 hectares) of the total area. Rubber was the dominant plantation type in the region, covering 17.20% (330,696.63 hectares) of the total area in 2018. However, oil palm plantation area increased significantly from 2011, reaching nearly half that of rubber plantations (114,072.25 hectares, 5.93%) (Figure 7).

Figure 7. Land cover changes from 2000 to 2018 across the study area.

3.2. Land-cover change results from 2000 to 2018

Figure 8 shows the dynamics of native forest change from 2000 to 2018. Over 18 years, native forest conversion to non-forest largely occurred in lowland and wetland regions. Native forest conversion to rubber plantations tended to occur in the hilly areas close to native forest in the north, and in the lowland forest in the centre and west of the region. Native forest conversion to oil palm plantations mainly occurred in the lowland areas in the east. However, a significant area of native forest regrowth was detected in peatlands in the southwest of the region (Figure 8).

Figure 8. Native forest conversion to non-forest, oil palm and rubber plantations, and native forest regrowth from 2000 to 2018 in Kalimantan.

Estimated areas of land cover conversion during the study period were derived from pixel-by-pixel analysis of the land cover map between 2000 and 2018, as shown in Table 1.

Table 1. Areas of land cover changed from 2000 to 2018 (Ha).

Land cover n 2000	Total in 2000	% of total land area in 2000	Land cover in 2018
			Native forest	Oil palm plantation	Rubber plantation	Non-forest
Native forest	1,390,080.69	72.28	753,047.32	25,290.88	153,790.13	457,952.37
Oil palm plantation	18,207.25	0.95	350.38	17,501.19	12.06	343.63
Rubber plantation	138,637.81	7.21	109.00	3321.63	127,718.06	7489.13
Non-forest	376,286.25	19.57	6548.88	67,958.56	49,176.38	252,602.44
Total in 2018	1,923,212.00		760,055.57	114,072.25	330,696.63	718,387.56
% of total land area in 2018			39.52	5.93	17.20	37.35

From 2000 to 2018 the study area in Kalimantan lost 630,025.12 hectares (45.32%) of native forest. More than 70% of this conversion was to non-forest (71.89%), followed by rubber plantation (24.14%) and oil palm plantation (3.97%). Non-forest areas increased by 342,101.31 hectares (90.92%), from 376,286.25 hectares in 2000 to 718,387.56 hectares in 2018 (Table 1). However, as previously shown in Figure 8, a relatively large area (6,548.88 hectares) of native forest regrowth from non-forest was detected during the study period (Table 1).

During the study period, rubber plantations increased by 192,058.81 hectares (138.53%) with native forest contributing 47% and non-forest contributing 15% of the new areas established (Table 1). These new rubber plantations that were established between 2000 and 2018 were largely distributed in the central region of the study area (Figure 5).

Oil palm plantations increased by 95,865 hectares (526.52%) from 18,207.25 hectares in 2000 to 114,072.25 hectares in 2018. Oil palm plantations mainly developed in non-forest areas (60%), followed by native forest conversion (22%). A small proportion of rubber plantations were converted to oil palm plantations (3%). Figure 5 shows the spatial distribution of the oil palm plantation expansion. Over 18 years, oil palm plantations largely developed in lowland areas in the east of the study area. A recent trend of oil palm plantation expansion was also detected in the west of the region (Figure 5).

The trend of land cover changes from 2000 to 2018 can be seen in Figure 9. Significant native forest loss occurred between 2015 and 2016. Rubber plantations increased significantly in the period from 2000 to 2006. The oil palm plantation area increased markedly between 2010 and 2011 and continued to increase until 2018.

Figure 9. Trend of land cover changes across the east Kalimantan study area from 2000 to 2018.

4. Discussion

This study developed a method to integrate a time series of Landsat, ALOS PALSAR, and Sentinel-1 satellite imagery, processed using Bayesian classification probability, to produce multitemporal land cover maps from 2000 to 2018 that differentiate between native forest, oil palm and rubber plantations, and non-forest. Over the 18-year time series, nearly half of the east Kalimantan study area native forests were converted to either non-forest, rubber plantations, or oil palm plantations.

The current land cover maps had an overall accuracy of 91–96%. This compares favourably to previous studies that showed an overall map accuracy of 73–77% using only optical images (Sari et al. 2021). Incorporating SAR images improved classification accuracy in discriminating closed-canopy cover (>50%), such as between native forests, oil palm and rubber plantations, due to SAR ability to penetrate the canopy and distinguish vegetation structural differences (Sari et al. 2022; Miettinen and Liew 2011). The time series Bayesian method was used to refine the individual year land cover classification to discriminate between native forest, oil palm and rubber plantation and non-forest.

ALOS PALSAR L-band images increased classification accuracy because of the greater penetration and sensitivity to vegetation structure to distinguish rubber plantations. L-band discriminated rubber plantations through the canopy with more distinct scattering from trunks and branches, twigs, and ground cover vegetation (Rosenqvist 1996). Whilst Sentinel-1 C-band was prominent in discriminating oil palm plantations based on marked differences between VV (high backscatter value) and VH (low backscatter value), due to the low density of stems (Sari et al. 2022; Miettinen and Liew 2011). Map validation in 2018 showed oil palm plantation had a relatively high user’s accuracy (>90%), and producer’s accuracy (>75%).

The highest accuracy was in mapping native forest and non-forest classes, with 88–99% for both user’s accuracy and producer’s accuracy. The oil palm class had user’s accuracy in the range of 92–98%, while producer’s accuracy varied from 61% to 78%. The rubber plantation class had user’s accuracy from 77% to 91% and producer’s accuracy from 77% to 87%. The method demonstrated lower map accuracy for oil palm plantations, particularly for areas in the earlier stages of plantation development. Similarly, there were higher errors in differentiation between young rubber and non-forests. For example, the final maps showed many small areas with mixes of oil palm plantation, non-forest, and young rubber plantation, particularly in the southeast of the region (Figure 5).

Future improvement of map accuracy can be achieved by using a higher resolution digital elevation model (DEM) for terrain correction, to reduce the uncertainty. SAR images are sensitive to topography that causes a bias, such as low certainty of forest biomass estimation, due to variation of the backscatter coefficients that are dependent on the incidence angle, wavelength, polarization, topography and land surface (CEOS 2018; Suzuki, Kim, and Ishii 2013). To improve DEM accuracy in spatial pixel spacing and elevation, FABDEM (Forest And Buildings removed Copernicus DEM) of TanDEM-X can be used that reduces a mean absolute vertical error of forests from 5.15 to 2.88 m (Hawker et al. 2022). Moreover, future improvements for oil palm plantation classification accuracy may be achieved by adding information of oil palm in the multitemporal Bayesian network developed using object-based segmentation or random forest from ALOS PALSAR (Ahmed et al. 2018; Dong et al. 2020; Sarzynski et al. 2020).

In agreement with previous studies, we showed that native forest loss in the east Kalimantan study area largely occurred in lowland and wetland areas, due to forest conversion, and possibly from forest fires between 2015 and 2016 (Figure 9), that covered approximately 849,703 hectares of the total area in east Kalimantan (Margono et al. 2014; Field et al. 2016; Decree of the Ministry of Environment and Forestry 2017). Fire data from MODIS and Suomi NPP Visible Infrared Imaging Radiometer Suite (VIIRS S-NPP), with an accuracy of >85%, showed that several fires occurred in the region between 2015 and 2016 (Figure 8) (NASA 2022).

A significant increase in oil palm plantation areas started in 2008. This is likely due to an increasing global demand for this product, with global palm oil production increasing by 242% from 120 Mt in 2000 to 410 Mt in 2018 (FAO 2021). Oil palm plantations largely grow below 150 m elevation and are located on mineral soils, in peatlands, coastal plains and swamp ecosystems (Ruysschaert et al. 2011). Oil palm plantation expansion from 2000 to 2018 was mainly in near-coastal areas of east Kalimantan. Recently oil palm plantation expansion was detected in lowlands, close to peatland in the west of the study region (scattered small areas) (Figure 8). Rubber plantations can grow at elevations up to 700 m and are mainly small-scale plantations (Phuc and Nghi 2014). We observed a significant rubber plantation expansion into native forest areas (153,790.13 hectares over 18 years) as well as in non-forest areas in lowland areas in the south of the region. This study also detected conversion of rubber plantations to oil palm plantations in the central region, which may have been driven by a relatively more profitable palm oil, and fuelled by the decline in rubber prices in 2011–2013 (Feintrenie, Chong, and Levang 2010; Krishna et al. 2017; ResearchInChina 2014).

Overall, over the 18 years, expansion of rubber and oil palm plantations contributed to native forest loss of approximately 179,081 hectares. Vast areas of plantation were also developed in previously non-forested land (117,134.94 hectares), such as degraded or cleared forests (Austin et al. 2017). Dramatic increase in rubber and oil palm plantations in areas of Kalimantan observed in this study is in line with a corresponding increase in global demand for rubber and palm oil products over the past two decades (FAO 2021; Wicke et al. 2011).

However, we also found a significant native forest regrowth in the degraded area of peatland in the study area, which can be a direct consequence of a recent Indonesian policy on peatland restoration (Presidential Instruction No. 5/2019). The results shown in this study can be an objective measure of the benefits of such a policy in terms of land restoration.

The time series of land cover maps developed in this study have demonstrated that they fulfil the Indonesian national standard of map accuracy of >85% and support sustainable plantation development in the Kalimantan region (Badan Standardisasi Nasional 2010; Badan Informasi Geospasial 2014). This information on land cover changes can be used to support Indonesia’s national forest monitoring system and sustainable forest management and investment decisions.

5. Conclusions

The current study used remote sensing data from Landsat, ALOS PALSAR, and Sentinel-1 satellite sensors over the period from 2000 to 2018 to improve automated mapping of spatial dynamics of native forest and plantation areas in east region of Kalimantan. Incorporating multisource SAR images of L- and C-bands and Landsat increased temporal data to improve land cover map accuracy, including ability to distinguish among vegetation types. The current study showed that nearly half of native forests in the east Kalimantan study area were converted to non-forest, or plantations of either rubber or oil palm between 2000 and 2018. Plantation expansion occurred mainly in lowland and wetland areas. However, significant native forest regrowth was also detected in peatland in the southwestern region of the study area over this period.

A multitemporal classification using a Bayesian network was applied to accurately map change and estimate areas of native forest conversion over the 18-year period. The maps had an overall accuracy above 90%. This fulfils the national map accuracy benchmark to comply with the IPCC requirement for transparent, repeatable, accurate, consistent, and complete land cover change mapping for sustainable forest management, including sustainable plantation development (Penman et al. 2003).

Accurate time series maps are critical to support Indonesia’s national forest monitoring system and to conserve the remaining native forests. Responsible management of this valuable resource will help to reduce emissions from forest degradation and deforestation, which is one of Indonesia’s strategies to reduce global greenhouse gas emissions (Indonesia LTS-LCCR 2050 2021). Moreover, land cover maps developed in this study provide insights into forest cover dynamics, including the rapid expansion of rubber and oil palm plantations. The ongoing generation of these land cover maps will help Indonesia to plan management strategies and policies that can balance the need for economic growth, and the need to conserve the valuable natural ecosystems and habitats across the region.

Acknowledgements

The authors thank the Research Center for Remote Sensing, National Research and Innovation Agency (BRIN), the LCCA teams, and the Indonesia–Australia Forest Carbon Partnership for sharing the satellite images and methods. I.L.S acknowledges the Indonesian LPDP for Ph.D scholarship.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author, ILS, upon reasonable request.