Assessing Accuracy of Land Cover Change Maps Derived from Automated Digital Processing and Visual Interpretation in Tropical Forests in Indonesia

Abstract

This study assessed the accuracy of land cover change (2000–2018) maps compiled from Landsat images with either automated digital processing or with visual interpretation for a tropical forest area in Indonesia. The accuracy assessment used a two-stage stratified random sampling involving a confusion matrix for assessing map accuracy and by estimating areas of land cover change classes and associated uncertainty. The reference data were high-resolution images from SPOT 6/7 and high-resolution images finer than 5 m obtained from Open Foris Collect Earth. Results showed that the map derived from automated digital processing had lower accuracy (overall accuracy 73–77%) compared to the map based on visual interpretation (overall accuracy 80–84%). The automated digital processing map error was in differentiating between native forest and plantation areas. While the visual interpretation map had a higher accuracy, it did not consistently differentiate between native forest and shrub areas. Future improvement of the digital map requires more accurate differentiation between forest and plantation to better support national forest monitoring systems for sustainable forest management.

1. Introduction

The fastest rate of forest degradation is occurring in the tropics, where more than 7 million hectares per year are deforested annually [1]. Addressing the drivers of forest loss requires accurate maps to enable a time-series analysis of land-use change to support the development of effective mitigation actions. While remotely sensed images can be automatically processed to predict global forest cover [2], the value of the maps for conservation management at local and regional scales relies on the accurate identification of land cover classes [3]. Furthermore, in many tropical areas, the quality of remotely sensed information is degraded by cloud cover or impacted by the nature of the terrain, posing technical challenges for accurate digital mapping.

This study addresses the development of land cover change maps for Indonesia, which has lost more than 68 million hectares to deforestation over recent decades (1950–2015) [4]. The rate of deforestation in Indonesia has varied over time, with a peak of around 1.9 million hectares per year in the 1990s, decreasing to around 0.68 million hectares per year between 2010 and 2015 [5,6]. The main drivers of this deforestation have also varied over the last few decades, from timber extraction (1966–1980) to conversion to shrublands, agriculture, oil palm, and timber plantations [7].

In 1990 the Indonesian Government established a national forest monitoring system (NFMS) that developed and applied a visual method to interpret Landsat images to map land cover according to 23 classes [8] that included 7 forest cover classes (primary dryland forest, secondary dryland forest, primary swamp forest, secondary swamp forest, primary mangrove forest, secondary mangrove forest, and plantation forest). Indonesia’s Ministry of Environment and Forestry (MoEF), with considerable input of local knowledge and ground-truthing, applied stratified random sampling to determine the accuracy of 88% to 90% for visual interpretation (VI) maps but did not provide fine details of the method for accuracy assessment and its uncertainty [8,9]. The VI mapping is subjective in the delineation of the land classes and is considered an overly time- and labor-intensive method for producing maps [8,10]. Although the VI map is currently Indonesia’s official map to support forest management [11,12], it is recognized that automated methods offer the ability to more rapidly and consistently map land cover classes in a measurable, repeatable, and verifiable way [8,13]. Further impetus for the development of automated mapping methods stems from the requirement for repeatable, accurate, and consistent land cover change mapping according to the Intergovernmental Panel on Climate Change (IPCC) principles of transparency, comparability, consistency, completeness, and accuracy [14].

Indonesia’s National Institute of Aeronautics and Space (LAPAN) has used automated digital processing (ADP) to map forest and non-forest land cover classes, and this map has been used by the MoEF for checking its VI maps [8,15]. The LAPAN ADP map was based on temporal Landsat images such as the Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) to produce two land cover classes of forest and non-forest that is mapped using different indices for each region of Indonesia [16]. However, a major drawback of these ADP maps is their inability to differentiate natural forest from tree plantations [17], thereby limiting their application in tracking the loss of native forest areas through time.

A statistically rigorous sampling-based accuracy assessment of land cover change maps is essential to quantify data quality thus that map users, such as policymakers, can evaluate the utility of maps to build trust and ensure confidence in public reporting, policy, and investment decisions. High-resolution imagery, such as images provided by the Open Foris Collect Earth or SPOT 6/7, in combination with a rigorous sampling strategy, can be used to verify and test the accuracy of land cover change maps. For example, a stratified random sampling with reference images from high resolution of the Google Earth images time-series was applied to determine the accuracy of the U.S. National Land Cover Database (NLCD) [18].

The primary aim of this study was to assess the map accuracy and the uncertainty of land cover classes of Indonesia’s ADP and VI maps against high-resolution reference images using statistically rigorous methods. This study will be used to identify the limitation of the ADP map for future improvement of the map accuracy and to reduce its uncertainty for building a robust national forest monitoring system, to support government policies that rely on the information of land cover and land cover change.

2. Materials and Methods

2.1. Experimental Approach and Study Area

To evaluate the accuracy of both ADP and VI maps, a sample study area of more than 1.9 million hectares was selected in east Kalimantan on the island of Borneo (0°50′20.32″ S–2°31′56.61″ S and 114°54′49.33″ E–116°37′22.46″ E), within the 4 regencies of Barito Timur, Tabalong, Balangan, and Paser [19,20,21]. The area varies topographically from swampy plains to hilly lowlands with elevations up to about 1550 m and with rainfall between about 2000 and 4000 mm per year. Notably, oil palm and rubber plantations have increased in areas in this region of Indonesia since the 1990s [22,23], thus that by 2017 over 350,000 ha of plantations had been established in the study area [19,20,21].

The steps in the accuracy assessment of both the Landsat based VI and ADP land cover change maps are set out below. To assess forest conversion to non-forest for the VI map, the following 6 land cover classes, from the total of 23 land classes mapped, were considered native forest (primary dryland forest, secondary dryland forest, primary mangrove forest, secondary mangrove forest, primary swamp forest, secondary swamp forest). The VI map includes a plantation forest class that was considered as non-forest for the purposes of compiling the land cover change maps. The land cover classes have also been used for determining forest cover in the National Forest Reference Emission Level for Deforestation and Forest Degradation that differentiates between natural forest and plantation forest [15].

The ADP map was classified into forest and non-forest. The forest cover class was defined as an area of trees with >5 m height, >30% canopy cover, and excludes plantations, such as oil palm. The ADP technique for mapping the annual forest extent and change maps is described in Section 2.2, broadly following the methodology outlined in [16]. Once the ADP forest extent and change maps were developed, their accuracy was assessed using a time-series of high-resolution images alongside the VI maps, as described in Section 2.3.

2.2. Developing Automated Digital Land Cover Maps for 2000–2018

The ADP maps for 2000–2018 used the Landsat time-series from the Land Cover Change Analysis (LCCA) project [16], based on images provided by the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/) data acquired in 2019. The methods are outlined in [16] and consist of image pre-processing (manual scene selection, orthorectification, radiometric terrain correction, cloud masking, and mosaic) and multitemporal land cover classification. Because manual scene selection was used to produce the cloud-free mosaic, a cloud masking process was required. Previous studies have shown limitations of optical images due to cloud cover, in which the addition of more temporal images (based on a higher cloud cover threshold) produces lower land cover classification accuracy compared to datasets composed with less cloud cover [24].

The ADP land cover classification detects forest and non-forest cover from spectral indices that are different for each region due to their unique combination of geography and vegetation type. The multitemporal land cover classification for the ADP maps was processed using the Canonical Probability Network [16] to produce annual forest and non-forest classes. The multitemporal processing reduces the uncertainty in the land cover classification using probabilistic rules in the iteration process and predicts missing data such as gaps caused by persistent cloud cover [25]. The target accuracy of the ADP map was 85%.

2.3. Assessment of the Land Cover Change Map Accuracy

The accuracy of ADP and VI land cover change maps was assessed by 2 methods: The confusion matrix, and by estimating the areas of land cover classes and associated uncertainty [26,27]. The land cover change maps were assessed in 3 periods of 6 years (2000–2006, 2006–2012, and 2012–2018). The VI maps were produced every 3 years, 2000, 2003, 2006, and 2009, and afterward annually between 2011 and 2018.

The assessments of the accuracy of the land cover change maps were based on a comparison of the map labels to the reference labels of high-resolution images.

The high-resolution images for assessing the accuracy of the map were obtained from:

• Commercial high-resolution SPOT 6/7 images (1.5 m resolution) from 2014 to 2019. Nearly a full coverage of SPOT 6/7 images was obtained for the year 2018. Whilst for the years before 2018, the availability of the SPOT 6/7 images was limited due to cloud cover. The SPOT 6/7 images were provided by LAPAN (https://inderaja-catalog.lapan.go.id/dd4/) data acquired in 2019.
• High-resolution images for the years 2000–2018 from the open-source platform Open Foris Collect Earth, developed by the Food and Agriculture Organization (FAO) of the United Nations [28]. Images with spatial resolution finer than 5 m were obtained from Google Earth and Microsoft Bing Maps, such as Quickbird, GeoEye-1, and Worldview-1 and -2 imagery [29,30].

2.3.1. Validation Data Selection

This study used the 2-stage stratified random sampling method due to the limitation of the high-spatial-resolution images used as the reference data, particularly for years before 2006. The stratified random sampling was used to maintain spatial randomness and to distribute the points proportionally across land cover classes [31,32]. The classes were formed from the difference in annual land cover map classes (forest and non-forest) of the ADP map. The land cover change classes used in this study were stable non-forest, stable forest, forest loss, and forest regrowth.

The 2-stage stratified random sampling used hierarchical grids for the primary sampling unit (PSU) and the secondary sampling unit (SSU). A comparable previous study used 4 pixels (30 m) for the SSU within each PSU of 6 × 6 km [33]. The current study used the PSU grids of 2 × 2 km resolution with the SSU block of 2 × 2 Landsat sized pixels (50 × 50 m) from the land cover change map. This SSU resolution was aligned with the minimum forested area (0.25 hectares) mapped in the VI map [11]. The selection of both PSU and SSU was based on the proportion of the land cover change classes within the stratification zone.

The overall sample size (n = 392 sample) was derived from 2 units of SSU, which were randomly selected in each of 196 PSU (4.08%) from a total of 4808 PSU. The overall sample size was comparable to a previous study that used 4.19% (n = 419 sample) from 10,000 simple size [34]. The overall sample size aimed to be distributed proportional to the probability of the land cover classes of stable non-forest, stable forest, forest loss, and forest regrowth within each stratification zone. To assess the accuracy of the map based on the distribution of the sample units per class and zones, refer to

Table 1. Distribution of the sample units per classes and zones (%) for automated digital processing (ADP) map in 2012–2018.

	Distribution Per Class	Distribution Per Zones *)
		Zone 3	Zone 9	Zone 2
Stable non-forest	26.9	81.4	7.8	10.8
Stable forest	51.2	89.0	4.5	6.5
Forest loss	12.6	70.8	18.0	11.2
Forest regrowth	9.3	89.4	6.8	3.8

*⁾ Percentage of the sample units per class distributed in each zone.

and Figure 1 for the ADP map in 2012–2018. A similar process was applied to the 2000–2006 and 2006–2012 time intervals, including for the VI maps, with sample distribution based on the availability of reference images.

The Open Foris Collect Earth and SPOT 6/7 based validation data were used to assess the accuracy of ADP and VI maps, with agreement assessed at the SSU level.

2.3.2. Confusion Matrix Accuracy Assessment

The map accuracy was assessed using the confusion matrix, the standard error, and the 95% confidence interval to measure the uncertainty of the map classes as outlined in [26,27,35]. The confusion matrix provided information on the producer’s accuracy, user’s accuracy, and overall accuracy. The overall accuracy was the sum of correctly classified pixels divided by the total sample size pixels. The user’s accuracy corresponded to errors of commission and was used to inform the user (consumer) on the reliability of the classifier that the areas classified as a specific land cover type or a certain percentage were correct. The producer’s accuracy corresponded to errors of omission. It was used to inform the analyst (producer) of the proportion of correctly classified areas in a particular category [35]. The overall producer’s and user’s accuracy in the confusion matrix of the population data and total population where $p_{ab}=n_{ab}$/n was determined from

Table 2. Population of data (n) and the map proportion area (p) confusion matrix with map as rows (m) and reference as columns (r).

		Reference (r)						Under/Over Estimation, SSU	Estimated Area, Hectares
		Class a	Class b	Class c	Class d	Total	${\mathit{A}}_{\mathit{m}}$
Map (m)	class a	$p_{aa}$	$p_{ab}$	$p_{ac}$	$p_{ad}$	$p_{a+}$	$A_{ma}$	$E_a$	$H_a$
	class b	$p_{ba}$	$p_{bb}$	$p_{bc}$	$p_{bd}$	$p_{b+}$	$A_{mb}$	$E_b$	$H_b$
	class c	$p_{ca}$	$p_{cb}$	$p_{cc}$	$p_{cd}$	$p_{c+}$	$A_{mc}$	$E_c$	$H_c$
	class d	$p_{da}$	$p_{db}$	$p_{dc}$	$p_{dd}$	$p_{d+}$	$A_{md}$	$E_d$	$H_d$
	Total	$p_{+a}$	$p_{+b}$	$p_{+c}$	$p_{+d}$	n
	$A_r$	$A_{ra}$	$A_{rb}$	$A_{rc}$	$A_{rd}$		A		H

where: Class a = stable non-forest, b = stable forest, c = forest loss, d = forest regrowth; the total map proportion area is 1; A is the numbers of sample units in SSU; $H_a$ is the estimate area of each land cover class in hectares.

.

The confusion matrix used the total sample size per class in map (m) (class a), referred to as $n_{a+}$, and the total sample size in reference data (r) (class a) referred to as $n_{+a}$. The overall accuracy was obtained by the diagonal sum of m and r divided by the total sample size n. The user’s accuracy for map class a was determined by dividing the sample number of class a by the total sample number of class a in m. The producer’s accuracy for reference class a was determined by dividing the sample number of class a by the total sample number of class a in r.

The variance of user’s accuracy, producer’s accuracy, and overall accuracy were based on the stratified random sampling equation described in [27]. The formulas assumed proportional allocation and ignored the cluster sampling feature of the 2 SSUs selected within each PSU. The estimated standard errors were, therefore, likely to be slight underestimates. The standard error with 95% confidence interval (z = 1.96) for the user’s accuracy for map class a, was estimated as follows:

$$\hat{\mathrm{S}}\hat{\mathrm{E}}\left({\widehat{U}}_a\right)=\mathrm{z}\sqrt{\frac{\frac{n_{aa}}{n_{a+}} \left(1-\frac{n_{aa}}{n_{a+}}\right)}{\left(n_{a+}-1\right)}}$$

(1)

Variance for producer’s accuracy for map class a, is estimated as follows:

$$\widehat{V}\left({\widehat{P}}_a\right)=\frac{1}{N_{+a}^2}\left[\frac{N_{a+}^2{(1-{\widehat{P}}_a)}^2{\widehat{U}}_a\left(1-{\widehat{U}}_a\right)}{n_{a+}-1}+{\widehat{P}}_a^2{{\displaystyle \sum }}_{i\ne a}^d{N}_{i+}^2\frac{n_{ia}}{n_{i+}}\left(1-\frac{n_{ia}}{n_{i+}}\right)/\left(n_{i+}-1\right)\right]$$

(2)

where $N_{+a}$ is the estimated marginal total number of SSU of reference class a, $N_{a+}$ is the estimated marginal total number of SSU of map class a.

Variance for overall accuracy, is estimated as follows:

$$\widehat{V}\left(\widehat{O}\right)={{\displaystyle \sum }}_{a=1}^d{W}_a^2{\widehat{U}}_a\left(1-{\widehat{U}}_a\right)/\left(n_{a+}-1\right)$$

(3)

where $W_a$ is area proportion of map class a.

The standard error with 95% confidence interval for producer’s accuracy is estimated as $\pm 1.96\sqrt{\widehat{V}\left({\widehat{P}}_a\right)}$ and for overall accuracy $\pm 1.96\sqrt{\widehat{V}\left(\widehat{O}\right)}$.

2.3.3. Land Cover Class Area and Uncertainty

The numbers of sample units for map in SSU for class a referred as $\left(A_{ma}\right)$ is estimated as the numbers of sample units (A) multiplied by $\left(p_{a+}\right).$ The numbers of sample units for reference in SSU for class a referred as $\left(A_{ra}\right)$ was estimated as the numbers of sample units (A) multiplied by $\left(p_{+a}\right)$. Map uncertainty (under or over-estimation based on the SSU) for class a is expressed as $E_a=A_{ma}-A_{ra}$.

The area of each land cover in hectares was estimated by $H_a=\left(A_{ma}\right)\ast k$, where $k$ = 0.25. The uncertainty in the area estimate was based on the variance, the standard error, and the 95% confidence interval as follows:

$$\hat{\mathrm{S}}\hat{\mathrm{E}}\left({\widehat{P}}_{+a}\right)=\sqrt{{{\displaystyle \sum }}_i{W}_i^2\frac{\frac{n_{ia}}{n_{i+}} \left(1-\frac{n_{ia}}{n_{i+}}\right)}{\left(n_{i+}-1\right)}}=\sqrt{{{\displaystyle \sum }}_i\frac{W_i{\widehat{P}}_{ia}-{\widehat{P}}_{ia}^2}{n_{i+}-1}}$$

(4)

where $W_i$ is the area proportion of map class i. The uncertainty area for class a (in hectares) is determined as $\hat{\mathrm{S}}\hat{\mathrm{E}}(H_a)=(H_a)\ast \hat{\mathrm{S}}\hat{\mathrm{E}}\left({\widehat{P}}_{+a}\right).$ The 95% confidence interval obtained as $H_a\pm 1.96\ast \hat{\mathrm{S}}\hat{\mathrm{E}}(H_a)$.

3. Results

3.1. Landcover Change 2000–2018

The ADP land cover change map (Figure 2) showed that forest loss in the initial period (2000–2006) largely occurred in the east and south-east of the region, dominated by dryland forests, and in the central region, dominated by lowland dipterocarp forests with elevation less than 300 m [36,37]. This initial deforestation was generally scattered within the existing stable forest (see Figure 2 zoom box A). Deforestation continued in the second period (2006–2012) and was detected in the south-east and the north-east of the region, dominated by lowland dipterocarp forests such as Shorea and Dipterocarpus with elevation less than 300 m [36,37,38]. The rapid forest loss was caused by the development of plantations (oil palm and rubber). During the third period (2012–2018), forest loss largely occurred in the western part of the region on the lowland wet-peat forests.

The forest regrowth was detected throughout the whole ADP mapping period and largely occurred in the area of oil palm and rubber plantation development (see in Figure 2 zoom box B and C). Forest regrowth was observed in the eastern, middle, and south-west of the region. The stable forested area identified from the ADP maps was predominantly in the dryland forest in the northern and middle of the region with elevation largely <500 m.

3.2. Map Accuracy

The 392 SSUs were used to assess map accuracy for each period. The map percent area of stable forest in the ADP map showed a gradual decrease from 65% in the period of 2000–2006 to 45% in 2012–2018 (Figure 3). The map percent area of stable non-forest increased steadily from 20% in 2000–2006 to 33% in 2012–2018, respectively. The map percent area of forest regrowth also increased from 1% in 2000–2006 to 8% in 2012–2018. However, the map percent area of forest loss was relatively stable 13% in 2000–2006, 11% in 2006–2012, and 14% in 2012–2018 (Figure 3).

The map percent area for the visual interpretation SSU showed a similar increase for stable non-forest and a decrease for the stable forest. However, for the VI map, there was a very small percent area of forest regrowth (ca. 0.4%) for 2006–2012 (Figure 3).

The overall accuracy of the VI map was higher compared to the ADP map for all the years of observations (

Table 3. Overall map accuracy of the ADP and VI maps for 2000–2018, standard error (SE) is given in parentheses.

Maps	2000–2006	2006–2012	2012–2018
ADP	76% (±5.6%)	77% (±5.8%)	73% (±4.4%)
VI	82% (±5.5%)	80% (±6%)	84% (±5.4%)

).

In terms of user’s accuracy, the ADP method was more accurate than the VI method in mapping non-forest but not as accurate in mapping stable forest areas (Figure 4). Assessment of stable forest in the ADP map for 2000–2006 showed a low user’s accuracy of 76% (±5.2%) due to the commission error of 24%. A similar commission error caused low user’s accuracy for the forest loss class of 52% (±13.7%). More importantly, the ADP method classified plantation as regrowth forest with a proportional area of 1.3% in the initial period (Figure 3). An example of ADP misclassification was the commission error in the stable forest class, caused by plantation (oil palm and rubber) being labeled as forest cover. The VI method had higher user’s accuracy in detecting stable forest and forest loss classes compared to the ADP method for all years (Figure 4). The map accuracy for each land cover class, based on the confusion matrix from 2000 to 2018, is shown in

Table A1. Description of population data and map proportion area confusion matrix for the automated digital processing of land cover change map in 2000–2006.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	75	2	2	0	79	0.95 (±0.049)	0.52(±0.021)	1,550,344	2,825,944	−1,275,600	387,586 (±3609)
	b) Map. pr	0.191	0.005	0.005	0	0.202
Stable forest	P. data	47	195	14	0	256	0.76 (±0.052)	0.95 (±0.069)	5,023,901	4,023,046	1,000,855	1,255,975 (±33,469)
	Map. pr	0.120	0.497	0.036	0	0.653
Forest loss	P. data	19	6	27	0	52	0.52 (±0.137)	0.63 (±0.107)	1,020,480	843,858	176,622	255,119 (±2932)
	Map. pr	0.048	0.015	0.069	0	0.133
Forest regrowth	P. data	3	2	0	0	5	0	0	98,123	0	98,123	24,530 (±108)
	Map. pr	0.008	0.005	0	0	0.013
Total	P. data	144	205	43	0	392	Overall 0.76 (±0.056)		7,692,848	7,692,848		1,923,212
	Map. pr	0.367	0.523	0.110	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

,

Table A2. Description of population data and map proportion area confusion matrix for the visual interpretation of land cover change map in 2000–2006.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	138	37	23	0	198	0.70 (±0.064)	0.96 (±0.088)	3,885,673	2,825,944	1,059,729	971,418 (±22,651)
	b) Map. pr	0.352	0.094	0.059	0	0.505
Stable forest	P. data	3	165	0	0	168	0.98 (±0.020)	0.80 (±0.019)	3,296,935	4,023,046	−726,111	824,234 (±6970)
	Map. pr	0.008	0.421	0	0	0.429
Forest loss	P. data	3	3	20	0	26	0.77 (±0.166)	0.47 (±0.067)	510,240	843,858	−333,618	127,560 (±1089)
	Map. pr	0.008	0.008	0.051	0	0.066
Forest regrowth	P. data	0	0	0	0	0	0	0	0	0	0	0
	Map. pr	0	0	0	0	0
Total	P. data	144	205	43	0	392	Overall 0.82 (±0.055)		7,692,848	7,692,848		1,923,212
	Map. pr	0.368	0.523	0.110	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

,

Table A3. Description of population data and map proportion area confusion matrix for the automated digital processing of land cover change map in 2006–2012.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	101	5	5	0	111	0.91 (±0.053)	0.60 (±0.029)	2,168,520	3,306,747	−1,138,228	542,130 (±7343)
	b) Map. pr	0.258	0.013	0.011	0	0.282
Stable forest	P. data	30	172	7	0	209	0.82 (±0.052)	0.93 (±0.065)	4,101,544	3,630,553	470,991	1,025,386 (±23,640)
	Map. pr	0.077	0.439	0.018	0	0.533
Forest loss	P. data	16	2	27	0	45	0.60 (±0.145)	0.70 (±0.117)	883,108	755,548	127,560	220,777 (±2545)
	Map. pr	0.041	0.005	0.069	0	0.115
Forest regrowth	P. data	22	6	0	0	28	0	0	539,677	0	539,677	134,919 (±1208)
	Map. pr	0.055	0.015	0	0	0.070
Total	P. data	169	185	39	0	392	Overall 0.77 (±0.058)		7,692,848	7,692,848		1,923,212
	Map. pr	0.430	0.472	0.098	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

,

Table A4. Description of population data and map proportion area confusion matrix for the visual interpretation of land cover change map in 2006–2012.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	157	45	11	0	213	0.74 (±0.059)	0.94 (±0.076)	4,170,230	3,277,310	892,920	1,042,558 (±25,591)
	b) Map. pr	0.399	0.115	0.028	0	0.542
Stable forest	P. data	8	135	8	0	151	0.90 (±0.049)	0.73 (±0.038)	2,953,504	3,630,553	−677,049	738,376 (±12,411)
	Map. pr	0.019	0.344	0.020	0	0.384
Forest loss	P. data	3	4	21	0	28	0.76 (±0.162)	0.53 (±0.079)	539,677	784,984	−245,308	134,919 (±1185)
	Map. pr	0.008	0.009	0.054	0	0.070
Forest regrowth	P. data	0	2	0	0	2	0	0	29,437	0	29,437	7,359 (±28)
	Map. pr	0	0.004	0	0	0.004
Total	P. data	167	185	40	0	392	Overall 0.80 (±0.060)		7,692,848	7,692,848		1,923,212
	Map. pr	0.426	0.472	0.102	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

,

Table A5. Description of population data and map proportion area confusion matrix for the automated digital processing of land cover change map in 2012–2018.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	124	4	1	0	129	0.96 (±0.033)	0.61 (±0.019)	2,531,575	4,003,421	−1,471,846	632,894 (±6697)
	b) Map. pr	0.316	0.010	0.003	0	0.329
Stable forest	P. data	26	149	2	0	177	0.84 (±0.054)	0.87 (±0.059)	3,473,556	3,355,809	117,748	868,389 (±18,049)
	Map. pr	0.066	0.380	0.005	0	0.452
Forest loss	P. data	35	7	14	0	56	0.25 (±0.114)	0.82 (±0.180)	1,098,978	333,618	765,360	274,745 (±3361)
	Map. pr	0.089	0.018	0.036	0	0.143
Forest regrowth	P. data	19	11	0	0	30	0	0	588,738	0	588,738	147,185 (±1446)
	Map. pr	0.048	0.028	0	0	0.077
Total	P. data	204	171	17	0	392	Overall 0.73 (±0.044)		7,692,848	7,692,848		1,923,212
	Map. pr	0.520	0.436	0.043	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

and

Table A6. Description of population data and map proportion area confusion matrix for the visual interpretation of land cover change map in 2012–2018.

Map		Reference					c) User’s Accuracy with SE in Parentheses	c) Producer’s Accuracy with SE in Parentheses	Number of SSU for Map d)	Number of SSU for Reference d)	Over/Under Estimate d) (in SSU)	c) Estimation Area with SE in Parentheses (Hectares)
		Stable Non-Forest	Stable Forest	Forest Loss	Forest Regrowth	Total
Stable Non-forest	a) P. data	184	40	1	0	225	0.82 (±0.051)	0.92 (±0.061)	4,415,538	3,924,922	490,615	1,103,884 (±26,795)
	b) Map. pr	0.469	0.102	0.003	0	0.574
Stable forest	P. data	11	137	2	0	150	0.91 (±0.045)	0.75 (±0.037)	2,943,692	3,571,679	−627,988	735,923 (±11,651)
	Map. pr	0.028	0.349	0.005	0	0.383
Forest loss	P. data	5	5	7	0	17	0.41 (±0.241)	0.70 (±0.255)	333,618	196,246	137,372	83,405 (±492)
	Map. pr	0.013	0.013	0.018	0	0.043
Forest regrowth	P. data	0	0	0	0	0	0	0	0	0	0	0
	Map. pr	0	0	0	0	0
Total	P. data	200	182	10	0	392	Overall 0.84 (±0.054)		7,692,848	7,692,848		1,923,212
	Map. pr	0.510	0.464	0.026	0	1

^a) P. data = population data, ^b) Map. pr. = map proportion area, ^c) SE with 95% confidence interval (CI), ^d) SSU is a block of 2 × 2 Landsat sized pixels (50 × 50 m).

.

3.3. Estimating the Area and Uncertainty of Mapped Land Cover Change Classes

The ADP method underestimated stable non-forest areas and overestimated other land classes for all years (Appendix A, Table A1, Table A2, Table A3, Table A4, Table A5 and Table A6). For example, in 2000–2006, the stable non-forest class was 706,486 ha, but the ADP estimated 387,586 ha. Therefore, the ADP map underestimated stable non-forest by 45% or 318,900 ha of the total non-forest areas. The ADP map indicated the extent of the stable forest at 1,255,975 ha, while the reference estimated 1,005,761 ha—an overestimate of 25% (250,213 ha) of the total stable forest areas. This trend was reversed for the VI method, where stable non-forest was overestimated by 38% (264,932 ha), stable forest underestimated by 18% (181,528 ha), and the forest loss class underestimated by 40% (83,404 ha). Moreover, the VI method showed some inconsistency in the attribution of the forest regrowth class due to interpreter subjectivity. For example, some instances of the interpreter attributing dry shrub in 2006 (non-forest class) and secondary dryland forest in 2009, while validation using the Collect Earth and the SPOT 7 images showed the areas were natural forest from 2000 to 2018 (Figure 2).

The ADP method detected a gradual increase of stable non-forest and forest regrowth, while stable forest decreased through the observation period. The stable non-forest class of the ADP map increased by 63% from 387,586 ha (±3609) in 2000–2006 to 632,894 ha (±6697) in 2012–2018. Nevertheless, the stable non-forest of the ADP map was smaller by 53% on average compared to the VI maps for all years. The decrease of the stable forest was detected from both the ADP and VI maps between 2000 and 2018. The stable forest estimated using the ADP map changed from 31% from 1,255,975 ha (±33,469) in 2000–2006 to 868,389 ha (±18,049) in 2012–2018. The stable forest estimated using the VI maps changed from 824,234 ha (±6970) in 2000–2006 to 735,923 ha (±11,651) in 2012–2018. The ADP map detected significant forest loss at 750,641 ha from 2000–2018, or 54% greater compared to the VI map. Yet, the ADP map detected larger forest regrowth area compared to the VI map for all periods.

4. Discussion

This study presents a comprehensive assessment of the accuracy of forest cover maps derived from automated digital processing and visual interpretation from 2000–2018. Based on the overall accuracy analysis, the current ADP map had lower accuracy compared to the VI map for all periods. The lower accuracy of the ADP map was due to commission error of labeling plantation (oil palm and rubber) as forest cover. Low user’s accuracy of the ADP map in the first period (2000–2006) led to an overestimate of stable forest area compared to the VI map. Moreover, the ADP method incorrectly mapped the dynamic growth of plantations in the later periods as forest regrowth (Figure 2). This low accuracy of the ADP method may be partially attributed to the reliance on optical images. Previous studies have shown that the use of optical images makes it difficult to differentiate between forest and dense canopy crops such as mature oil palm [17,39] as well as perennial agricultural crops and plantations of less than 5 m height such as pineapple, tea, and soybean as tree cover [3,40] due to similar spectral signal.

Another challenge of using optical images for monitoring forest extent and change in tropical regions is persistent cloud cover [41]. For example, only 5% of Landsat image scenes for Kalimantan and Papua in 2010 were free from cloud cover [16]. This limits the application of optical images to produce continuity of temporal data. Therefore, the accuracy of the current map may be improved by using multi-source data such as radar and optical systems. Radar systems such as Sentinel-1 can penetrate clouds, smoke, haze and operate in both day and nighttime modes [42,43].

In agreement with other studies, the multi-temporal classification process was useful in reducing uncertainty [25,44,45]. The other advantage of multi-temporal processing is to estimate land cover classes in areas of missing data due to cloud cover, using interpolation and extrapolation of the maps spatially and temporally [16].

Consistent with previous studies [46], the VI process in this study was superior for differentiating between plantation and natural forest. The VI method uses a comprehensive and contextual assessment of land cover based on an interpreter’s knowledge. The interpreter uses visual pattern recognition to detect an object based on shape, size, pattern, texture, shadow, site, and association. Moreover, one advantage of VI is in detecting heterogeneity of land cover that affects spectral response [47]. However, the VI map was inconsistent in detecting forest regrowth because accuracy varies with interpreter training, experience, and skill. Because VI is labor-intensive and time-consuming, the VI map consistency is also subject to operator fatigue [48,49]. The overall accuracy of the VI map was 80–84%, which was in contrast to previous studies claiming above 90% accuracy [8,9]. Such a discrepancy can be attributed to a more comprehensive assessment method used in our study as well as reliance on high resolution satellite such as SPOT 6/7 imagery rather than Landsat used for the accuracy assessment in previous studies [27,50].

In accordance with the IPCC principles, the ADP method produced land cover maps in a repeatable and operationally efficient way. However, improvement of map accuracy is required.

The current study used the two-stage stratified random sampling method for accurate assessment of the maps because the method is cost-effective in assessing large areas with limited reference data [33,51]. The method stratified the land cover classes based on geographical stratification zones, with the flexibility to allocate the sample based on the proportion of each land cover class. The two-stage stratified random sampling method is suitable for assessing maps with a spatial unit larger than a pixel (blocks) and can incorporate reference data error such as imprecise position [52]. The stratified sampling method also provides a precise estimate for both the user’s accuracy and the producers’ accuracy, and it estimates a variance. An increase in precision may be obtained by selecting fewer sample units within each PSU, and re-allocating these sample units to additional PSUs [52].

The limitation of the current study was the small sample size that requires further steps of cross-validation for multiple processing times to increase the accuracy of the results [34]. However, advanced technologies provide the availability of various high-resolution images of sufficient temporal representation allowing the selection of large sample sizes to increase accuracy that is suitable for assessing large geographic areas with heterogeneous class [27,34].

5. Conclusions

The Indonesian Government uses the VI maps supported by the ADP maps for national forest monitoring using Landsat images as the main input data. Accuracy assessment of the land cover change maps is essential to quantify data quality for future improvement of Indonesia’s national forest monitoring system. The improvement of ADP maps is required to develop a robust NFMS and for the accurate updating and quality control of the VI maps to monitor the forested area and its changes. We applied the two-stage stratified random sampling process to assess the accuracy of the ADP and VI maps over 2000–2018 for Kalimantan. The two-stage stratified random sampling provides a cost-effective method that is based on limited reference data. The stratified sampling provides precision for user’s accuracy, producers’ accuracy, and its estimate variances are important for reporting precision of estimate area and understanding sources of errors in maps.

The accuracy of the ADP and the VI maps was assessed using two methods: (1) A confusion matrix with overall accuracy, user’s accuracy, and producer’s accuracy that was quantified using reference sample units based on manual interpretation of SPOT 6/7 images (1.5-m resolution), and the temporal high-resolution images from the Open Foris Collect Earth, and (2) by estimating the uncertainty associated with areas of land cover change classes. The result showed that the ADP method produced lower overall accuracy compared to the VI maps. This is primarily because of the commission error in mapping plantation as natural forest cover due to their similar spectral signatures that caused high variances in the forest loss class.

The ADP map can be improved by refining the classification method to differentiate plantation from natural forest. Improvement of the ADP map accuracy is important to estimate forest area, carbon stocks, and changes as a part of a rigorous national forest monitoring system. A repeatable, accurate, and consistent land cover change map would comply with the IPCC principles of transparency, comparability, consistency, completeness, and accuracy and would allow policymakers to track and manage land cover change for sustainable forest management and investment decisions.