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

: 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 stratiﬁed 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 ﬁner 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.


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

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.

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

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:

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 and Figure 1  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.

Confusion Matrix Accuracy Assessment
The map accuracy was assessed using the confusion matrix, the standard error, and  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.

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. 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 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: Variance for producer's accuracy for map class a, is estimated as follows: 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: where W a is area proportion of map class a.
The standard error with 95% confidence interval for producer's accuracy is estimated as ±1.96 V P a and for overall accuracy ±1.96 V Ô .

Land Cover Class Area and Uncertainty
The numbers of sample units for map in SSU for class a referred as (A ma ) is estimated as the numbers of sample units (A) multiplied by (p a+ ). The numbers of sample units for reference in SSU for class a referred as (A ra ) was estimated as the numbers of sample units (A) multiplied by (p +a ). 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 = (A ma ) * 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: where W i is the area proportion of map class i. The uncertainty area for class a (in hectares) is determined asŜÊ(H a ) = (H a ) * ŜÊ P +a . The 95% confidence interval obtained as H a ± 1.96 * ŜÊ(H a ).

Landcover Change 2000-2018
The ADP land cover change map ( Figure 2) showed that forest loss in the initial period (2000)(2001)(2002)(2003)(2004)(2005)(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)(2007)(2008)(2009)(2010)(2011)(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.

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

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 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 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). 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 Tables A1-A6. The overall accuracy of the VI map was higher compared to the ADP map for all the years of observations (Table 3). 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 Tables A1-A6.

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, Tables A1-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.

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)(2001)(2002)(2003)(2004)(2005)(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 comprehen-sive 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 highresolution 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].

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 costeffective 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.