Causal Analysis of Accuracy Obtained Using High-Resolution Global Forest Change Data to Identify Forest Loss in Small Forest Plots

: Identifying areas of forest loss is a fundamental aspect of sustainable forest management. Global Forest Change (GFC) datasets developed by Hansen et al. (in Science 342:850–853, 2013) are publicly available, but the accuracy of these datasets for small forest plots has not been assessed. We used a forest-wide polygon-based approach to assess the accuracy of using GFC data to identify areas of forest loss in an area containing numerous small forest plots. We evaluated the accuracy of detection of individual forest-loss polygons in the GFC dataset in terms of a “recall ratio”, the ratio of the spatial overlap of a forest-loss polygon determined from the GFC dataset to the area of a corresponding reference forest-loss polygon, which we determined by visual interpretation of aerial photographs. We analyzed the structural relationships of recall ratio with area of forest loss, tree species, and slope of the forest terrain by using linear non-Gaussian acyclic modelling. We showed that only 11.1% of forest-loss polygons in the reference dataset were successfully identiﬁed in the GFC dataset. The inferred structure indicated that recall ratio had the strongest relationships with area of forest loss, forest tree species, and height of the forest canopy. Our results indicate the need for careful consideration of structural relationships when using GFC datasets to identify areas of forest loss in regions where there are small forest plots. Moreover, further studies are required to examine the structural relationships for accuracy of land-use classiﬁcation in forested areas in various regions and with di ﬀ erent forest characteristics.


Introduction
Because forest losses greatly affect the benefits derived from the preservation of ecosystems [1][2][3], understanding their spatial distribution is important for effective forest management. Many research projects have identified the ecological importance of the spatial patterns of forest and areas of forest degradation or loss [4,5]. For example, forest loss in riparian areas often results in deterioration of stream water quality [6][7][8]. However, appropriate application of a mosaic approach to early stage reforestation can improve biodiversity by providing habitats for specific species [9,10]. Understanding the regional distribution of areas of forest, forest loss, and forest degradation is also useful for taking stock of the current situation, predicting future conditions, and undertaking sustainable forest planning [11][12][13]. Thus, detecting forest losses is a fundamental requirement for sustainable forest management.
Monitoring forest disturbances is a major interest in the remote-sensing community and many researchers have produced forest disturbance maps based on time-series satellite imagery at regional or

Study Area
Our study area is in the Bungo-Ono municipality of Oita Prefecture, south-west Japan (centered on 32 • 58 N, 131 • 35 E; Figure 1a,b). Its total area is 60,314 ha and it includes 44,659 ha of forest. Apart from areas of urban development and cropland, almost all of the study area is mountainous, with a mean slope of 28.8 • in forested areas. The region is dominated by plantations of Japanese cedar (Cryptomeria japonica), Japanese cypress (Chamaecyparis obtusa), and red pine (Pinus densiflora), and deciduous broad-leaved forest including sawtooth oak (Quercus acutissima) (Figure 1c). Within the study area, 84% of the forested area is privately owned and 69% of the private forest plots are of less than 1 ha in area. Forestry is a major industry in the region and timber supply in the region has been required to meet increasing demand in recent years. Because there have been no large-scale natural disasters in the region for at least a decade, human activities have been the main cause of forest loss. Most of the harvesting of timber in the study area during the recent decades has been by clear-cutting; there have been episodes of moderate forest thinning.

Reference Data
We generated a reference map of areas of forest loss by using aerial photographs acquired in 2014 and 2016 (Figure 2a,b). In 2014, photos were taken of the northern part of the study area on 1 June and those of the southern part on 9 September. In 2016, photos were taken of the entire study area on 8 December. We generated three-dimensional point clouds from both vintages of aerial photographs by using the Structure from Motion approach in PhotoScan Professional software (Agisoft LLC, Saint-Petersburg, Russian Federation). Both point clouds were then converted to digital surface models (DSMs) and ortho-images were produced at 0.5-m resolution. Then, the 2016 DSM was subtracted from the 2014 DSM and locations where the difference exceeded 5 m were extracted to represent areas of forest loss between 2014 and 2016. The forest loss map for the entire study area was then visually checked for credibility. We used the 2014 ortho-images to visually identify four categories of forest species present at that time, Japanese cedar, cypress, red pine, broad-leaved trees, as well as non-forested areas. Given the high resolution of the aerial photographs and the careful

Reference Data
We generated a reference map of areas of forest loss by using aerial photographs acquired in 2014 and 2016 (Figure 2a,b). In 2014, photos were taken of the northern part of the study area on 1 June and those of the southern part on 9 September. In 2016, photos were taken of the entire study area on 8 December. We generated three-dimensional point clouds from both vintages of aerial photographs by using the Structure from Motion approach in PhotoScan Professional software (Agisoft LLC, Saint-Petersburg, Russian Federation). Both point clouds were then converted to digital surface models (DSMs) and ortho-images were produced at 0.5-m resolution. Then, the 2016 DSM was subtracted from the 2014 DSM and locations where the difference exceeded 5 m were extracted to represent areas of forest loss between 2014 and 2016. The forest loss map for the entire study area was then visually checked for credibility. We used the 2014 ortho-images to visually identify four categories of forest species present at that time, Japanese cedar, cypress, red pine, broad-leaved trees, as well as non-forested areas. Given the high resolution of the aerial photographs and the careful interpretation of the aerial photographs, the uncertainty in the reference dataset was considered to be sufficiently small.

Global Forest Cover Change Datasets
We used the 10° lat × 10° long annual forest-loss map data (top left corner at 40°N, 130°E) from the GFC dataset ver. 1.6 [17] and used SAGA GIS v.2.3.2 open-source software to vectorize the data. To cover the period of our reference forest-loss map, we extracted polygons for forest losses from June 2014 to December 2016 and masked the data outside our study area (Figure 2c,d).

Accuracy Assessment
To assess the accuracy of the GFC dataset we compared the areal extent, local slope, and tree species of the forest-loss polygons obtained from the GFC dataset with those from our reference dataset. Because forests of the study area are too dense to generate digital terrain model (DTM) from

Global Forest Cover Change Datasets
We used the 10 • lat × 10 • long annual forest-loss map data (top left corner at 40 • N, 130 • E) from the GFC dataset ver. 1.6 [17] and used SAGA GIS v.2.3.2 open-source software to vectorize the data. To cover the period of our reference forest-loss map, we extracted polygons for forest losses from June 2014 to December 2016 and masked the data outside our study area (Figure 2c,d).

Accuracy Assessment
To assess the accuracy of the GFC dataset we compared the areal extent, local slope, and tree species of the forest-loss polygons obtained from the GFC dataset with those from our reference dataset. Because forests of the study area are too dense to generate digital terrain model (DTM) from aerial photographs, slope was calculated at 5-m resolution based on a 5-m DTM published by the Geospatial Information Authority of Japan. The raster ortho-image showing the initial (2014) tree species was Remote Sens. 2020, 12, 2489 5 of 13 resampled to 5-m resolution to reduce computation time and then masked with the forest-loss polygons derived from the 2016 GFC forest-loss or the reference dataset. Areas classified as non-forest were included in forest-loss polygons because the raster image that we used for land classification was resampled to a lower resolution; consequently, the edges of some of the polygons categorized as forest loss may actually represent non-forest areas ( Figure 3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 14 aerial photographs, slope was calculated at 5-m resolution based on a 5-m DTM published by the Geospatial Information Authority of Japan. The raster ortho-image showing the initial (2014) tree species was resampled to 5-m resolution to reduce computation time and then masked with the forest-loss polygons derived from the 2016 GFC forest-loss or the reference dataset. Areas classified as non-forest were included in forest-loss polygons because the raster image that we used for land classification was resampled to a lower resolution; consequently, the edges of some of the polygons categorized as forest loss may actually represent non-forest areas ( Figure 3). To quantify the accuracy of the GFC dataset in our study area, we used a method developed by Linke et al. [28]. We overlaid the forest-loss polygons of the reference data and GFC data and calculated the areas of overlap. We defined the ratio of the area of overlap of each GFC polygon to the entire area of the corresponding reference polygon as the "recall ratio", and the ratio of the area of overlap of each GFC polygon with the corresponding reference polygon to the entire area of the GFC polygon as the "precision ratio" (Figure 4). We considered reference polygons with high recall ratios to have been detected in the GFC dataset, and GFC polygons with high precision ratios to have detected the forest loss of the reference data.
To investigate the influence of polygon area on the accuracy of the GFC dataset, we estimated the distributions of recall and precision ratios for five ranges of reference and GFC polygon area (ha) determined on the basis of the number of GFC raster cells they enclosed: ≤36, >0. 36 Reference forest-loss polygons with recall ratios >50% were considered to have been successfully detected by from the GFC dataset. Likewise, GFC forest-loss polygons with precision ratios >50% were considered to have successfully detected the forest loss of the reference dataset. We also compared the proportions of the areas of these polygons for the five ranges of reference polygon area defined above.
We made all of the above calculations by using the Pandas v. 1.0.3 and Geopandas v. 0.7.0 modules of the Python v. 3.8.1 programing language. The confidence intervals of the accuracies were not required, because of the forest-wide assessment approach we used in this study. To quantify the accuracy of the GFC dataset in our study area, we used a method developed by Linke et al. [28]. We overlaid the forest-loss polygons of the reference data and GFC data and calculated the areas of overlap. We defined the ratio of the area of overlap of each GFC polygon to the entire area of the corresponding reference polygon as the "recall ratio", and the ratio of the area of overlap of each GFC polygon with the corresponding reference polygon to the entire area of the GFC polygon as the "precision ratio" (Figure 4). We considered reference polygons with high recall ratios to have been detected in the GFC dataset, and GFC polygons with high precision ratios to have detected the forest loss of the reference data.

Inferences on a Structural Relationship
For local forest management, overlooking the occurrence of forest loss should be avoided in order to reduce disaster risk. Determining whether factors such as the size of the areas of forest loss or the local slope directly influence the accuracy of forest loss determined from the GFC dataset (and if so, by how much) should further indicate the utility of the GFC data for this purpose. Since identifying where omission occurs is important in local forest management, relationships between recall ratio and forest losses were analyzed. We used a linear non-Gaussian acyclic model (LiNGAM) to infer a structure of the direct relationships. LiNGAM uniquely and empirically identifies a model structural equation with a semiparametric approach that supposes linearity for relationship functions and non-Gaussian distributions for error variables [29]. Using the assumption of linearity for causal Reference forest-loss polygons with recall ratios >50% were considered to have been successfully detected by from the GFC dataset. Likewise, GFC forest-loss polygons with precision ratios >50% were considered to have successfully detected the forest loss of the reference dataset. We also compared the proportions of the areas of these polygons for the five ranges of reference polygon area defined above.
We made all of the above calculations by using the Pandas v. 1.0.3 and Geopandas v. 0.7.0 modules of the Python v. 3.8.1 programing language. The confidence intervals of the accuracies were not required, because of the forest-wide assessment approach we used in this study.

Inferences on a Structural Relationship
For local forest management, overlooking the occurrence of forest loss should be avoided in order to reduce disaster risk. Determining whether factors such as the size of the areas of forest loss or the local slope directly influence the accuracy of forest loss determined from the GFC dataset (and if so, by how much) should further indicate the utility of the GFC data for this purpose. Since identifying where omission occurs is important in local forest management, relationships between recall ratio and forest losses were analyzed. We used a linear non-Gaussian acyclic model (LiNGAM) to infer a structure of the direct relationships. LiNGAM uniquely and empirically identifies a model structural equation with a semiparametric approach that supposes linearity for relationship functions and non-Gaussian distributions for error variables [29]. Using the assumption of linearity for causal functions, causality is inferred in LiNGAM according to the equation where x i and e i are continuous observed and exogenous variables, b ij is the connection strength from x j to x i [30]. Equation (1) can be converted into matrix form, where x = (x i | i ∈ 1, · · · , n) and e = (e i | i ∈ 1, · · · , n) are vectors of observed and exogenous variables, respectively, B = b ij i, j ∈ 1, · · · , n is an n × n matrix of coefficients, where n is the number of observed variables. If b ij is not zero, x i has a direct influence on x j and larger values of b ij represent stronger connections. Diagonal components of B (i.e., b ii ) are zero because the model is assumed to be acyclic. LiNGAM infers B from observed data using either an independent component analysis [29] or an approach using regression analysis and examination of independence (DirectLiNGAM algorithm) [31]. We used the DirectLiNGAM algorithm in this study. The parameters we input to LiNGAM as observed data for each reference forest-loss polygon were area of forest loss, mean slope, mean value of the digital canopy model (DCM), dominant tree species, and recall ratio. DCM was estimated by subtracting the 2014 DEM from the 2014 DSM. Although tree species is a categorical variable, we assigned four real numbers to represent an order of ease of the detection [32][33][34]

Descriptions of Forest Loss
The GFC dataset within our study area contained 1141 forest-loss polygons with a total area of 586.7 ha. The reference data contained 1480 polygons with a total area of 780.0 ha. The mean areas of Remote Sens. 2020, 12, 2489 7 of 13 forest loss were 0.51 ha for the GFC dataset and 0.52 ha for the reference data. Although these means are almost the same, their distributions differ. The distribution in the GFC data was biased towards polygons of smaller area than those of the reference data ( Figure 5). This difference might reflect the different intervals of data collection for the two sources: the reference data identified the forest loss over the two-year interval between the summers of 2014 and 2016, whereas the GFC data identified annual forest loss. Consequently, an apparently single period of forest loss in the reference data might be represented by more than one polygon in the GFC dataset.

Descriptions of Forest Loss
The GFC dataset within our study area contained 1141 forest-loss polygons with a total area of 586.7 ha. The reference data contained 1480 polygons with a total area of 780.0 ha. The mean areas of forest loss were 0.51 ha for the GFC dataset and 0.52 ha for the reference data. Although these means are almost the same, their distributions differ. The distribution in the GFC data was biased towards polygons of smaller area than those of the reference data ( Figure 5). This difference might reflect the different intervals of data collection for the two sources: the reference data identified the forest loss over the two-year interval between the summers of 2014 and 2016, whereas the GFC data identified annual forest loss. Consequently, an apparently single period of forest loss in the reference data might be represented by more than one polygon in the GFC dataset.
. Forest losses identified from the GFC dataset were on slightly steeper slopes than those in the reference data ( Figure 6). The mean slopes of these areas in the GFC and reference datasets were 25.9° and 24.0°, respectively. The mean slope of the entire forested part of the study area was 28.8°. Forest losses recorded in both datasets were from forest on slopes gentler than the mean slope for the entire forested area, likely reflecting preferential selection of gentler slopes for forest harvesting. Forest losses identified from the GFC dataset were on slightly steeper slopes than those in the reference data ( Figure 6). The mean slopes of these areas in the GFC and reference datasets were 25.9 • and 24.0 • , respectively. The mean slope of the entire forested part of the study area was 28.8 • . Forest losses recorded in both datasets were from forest on slopes gentler than the mean slope for the entire forested area, likely reflecting preferential selection of gentler slopes for forest harvesting. According to both the reference and GFC datasets, the proportions of the losses of Japanese cedar and cypress forest from 2014 to 2016 were considerably higher than their proportions in the entire forested area in 2014 (Figure 7). This difference is probably because these species are planted for timber production, so they are harvested more often than other species. The amount of broad-leaved forest loss was greater in the reference dataset than in the GFC dataset, which suggests that pre- According to both the reference and GFC datasets, the proportions of the losses of Japanese cedar and cypress forest from 2014 to 2016 were considerably higher than their proportions in the entire forested area in 2014 (Figure 7). This difference is probably because these species are planted for timber production, so they are harvested more often than other species. The amount of broad-leaved forest loss was greater in the reference dataset than in the GFC dataset, which suggests that pre-existing forest species might influence the performance of the GFC dataset for detection of forest loss. According to both the reference and GFC datasets, the proportions of the losses of Japanese cedar and cypress forest from 2014 to 2016 were considerably higher than their proportions in the entire forested area in 2014 (Figure 7). This difference is probably because these species are planted for timber production, so they are harvested more often than other species. The amount of broad-leaved forest loss was greater in the reference dataset than in the GFC dataset, which suggests that preexisting forest species might influence the performance of the GFC dataset for detection of forest loss.

Recall and Precision Ratios
Only 11.1% of the reference forest-loss polygons had recall ratios >50%, indicating that forest loss in our study area would be under-reported on the basis of the GFC dataset. The polygons of the reference dataset with recall ratios <50% amounted to 488.3 ha, whereas those with recall ratios >50% amounted to only 291.8 ha. However, forest losses were more successfully detected in the larger polygons ( Figure 8). Of the forest losses with areas >2.25 ha, 60.4% had recall ratios higher than 50%. These results suggest that many of the small forest-loss polygons were not detected in the GFC dataset. The predominance of recall ratios <50% for all reference polygons of <2.25 ha area indicates that other factors might have affected the recall ratios obtained.
Our results for precision ratio presented higher accuracies. Of the forest-loss polygons from the GFC dataset, 57.3% had precision ratios >50%. Precision ratio improved with increasing polygon area, in particular in the uppermost two ranges of polygon size (Figure 9). In the lowest size range (<0.36 ha), 49.5% of GFC polygons had precision ratios >50%, but in the highest size range (>2.25 ha), 96.6% had precision ratios >50%. As was the case for recall ratio, the distributions of precision ratio demonstrated that larger GFC polygons returned better results. Success in the identification of forest loss from the GFC dataset on the basis of precision ratio was clearly superior to that measured by recall ratio, particularly for larger areas of forest loss. However, the greater predominance of smaller polygons in the GFC dataset compared to the reference dataset ( Figure 5) may have compromised the overall results of detection of forest loss from the GFC data. had precision ratios >50%. As was the case for recall ratio, the distributions of precision ratio demonstrated that larger GFC polygons returned better results. Success in the identification of forest loss from the GFC dataset on the basis of precision ratio was clearly superior to that measured by recall ratio, particularly for larger areas of forest loss. However, the greater predominance of smaller polygons in the GFC dataset compared to the reference dataset ( Figure 5) may have compromised the overall results of detection of forest loss from the GFC data.
.  For each size range, the black dots mark the proportion of forest-loss polygons successfully detected (precision ratio >50%) in the GFC dataset.

Inferences on a Structural Relationship
The relationship structure determined by our application of the DirectLiNGAM algorithm indicates that area of forest loss, DCM, and tree species directly affected the recall ratios we obtained ( Figure 10). As noted in the preceding section, we showed that larger areas of forest loss were better identified by the GFC dataset. Higher mean DCM values and sharper crown forms also had positive effects on recall ratio. These results indicate that forest losses that provided large amounts of timber (i.e., large areas with high DCM) from major forestry species (Japanese cedar and Japanese cypress)

Inferences on a Structural Relationship
The relationship structure determined by our application of the DirectLiNGAM algorithm indicates that area of forest loss, DCM, and tree species directly affected the recall ratios we obtained ( Figure 10). As noted in the preceding section, we showed that larger areas of forest loss were better identified by the GFC dataset. Higher mean DCM values and sharper crown forms also had positive effects on recall ratio. These results indicate that forest losses that provided large amounts of timber (i.e., large areas with high DCM) from major forestry species (Japanese cedar and Japanese cypress) were more easily detected in the GFC dataset. Slope had no immediate relationship with recall ratio in our study area, although both area of forest loss and DCM showed positive associations with slope. These associations can possibly be explained in terms of the profitability of forestry operations. If both the area of harvestable timber and the height of the canopy (DCM) were large, harvesting of forest stands might be economically viable even on relatively steep slopes.

Discussion
Our aim was to examine the accuracy of the use of the GFC dataset [17] to identify areas of forest loss. Because many past analyses of the accuracy of land-cover maps have been done using a probability sampling approach [35,36], assessment methods based only on selected samples have been extensively developed and applied [27]. In contrast, we used a polygon-based forest-wide assessment approach. In particular, we aimed to analyze the effectiveness of the use of the GFC data to identify forest loss in a forest region with a substantial number of small forest plots, where we could assess the resolution of small polygonal areas of forest loss in the GFC dataset. Such a polygonbased analysis has practical advantages for users of the GFC dataset who require forest management at the forest-stand scale, such as the owners of the small forest plots in our study area. The availability of maps that satisfy specific users' requirements will contribute to sustainable forest management [27,37].
The accuracies of our results were considerably lower than those of previous analyses (e.g., accuracies of 83.8% in [25] and 57.6% in [26] using pixel-based sampling points); we successfully detected forest loss (>50% recall ratio) for only 11.1% of our reference polygons. However, we achieved greater accuracy for the larger reference polygons. Thus, we infer the presence of a considerable number of small forest-loss polygons to be the reason for the poorer overall performance of our approach. Linke et al. [28] also reported that the accuracy of the GFC dataset was low for forestloss polygons of <2 ha. These observations are consistent with the generally accepted view that patch area has a strong influence on the accuracy of interpretations of satellite imagery (e.g., [38,39]).
The results of our inference modelling on a structural relationship suggest that area of forest loss,

Discussion
Our aim was to examine the accuracy of the use of the GFC dataset [17] to identify areas of forest loss. Because many past analyses of the accuracy of land-cover maps have been done using a probability sampling approach [35,36], assessment methods based only on selected samples have been extensively developed and applied [27]. In contrast, we used a polygon-based forest-wide assessment approach. In particular, we aimed to analyze the effectiveness of the use of the GFC data to identify forest loss in a forest region with a substantial number of small forest plots, where we could assess the resolution of small polygonal areas of forest loss in the GFC dataset. Such a polygon-based analysis has practical advantages for users of the GFC dataset who require forest management at the forest-stand scale, such as the owners of the small forest plots in our study area. The availability of maps that satisfy specific users' requirements will contribute to sustainable forest management [27,37].
The accuracies of our results were considerably lower than those of previous analyses (e.g., accuracies of 83.8% in [25] and 57.6% in [26] using pixel-based sampling points); we successfully detected forest loss (>50% recall ratio) for only 11.1% of our reference polygons. However, we achieved greater accuracy for the larger reference polygons. Thus, we infer the presence of a considerable number of small forest-loss polygons to be the reason for the poorer overall performance of our approach. Linke et al. [28] also reported that the accuracy of the GFC dataset was low for forest-loss polygons of <2 ha. These observations are consistent with the generally accepted view that patch area has a strong influence on the accuracy of interpretations of satellite imagery (e.g., [38,39]).
The results of our inference modelling on a structural relationship suggest that area of forest loss, tree species, and the height of the forest canopy are the major influences on accuracy of detection of forest loss. Previous studies have investigated the importance of geographic [40][41][42] and other forest characteristics [43] on the accuracy of land classification based on satellite imagery. However, there has been little research on the causal relations between interpretation accuracy and forest conditions. Gislason et al. [41] reported that geographic features are an important influence on interpretation accuracy, but demonstrated no immediate relationship between terrain slope and accuracy in this study. Although at first glance the difference of the distributions of slope in the reference and GFC datasets ( Figure 4) seems to suggest that slope would affect accuracy, we found no immediate relationship. The relationships that we identified of slope with both polygon area and DCM might have led to this pseudo correlation ( Figure 8). Clearly, understanding a structure of relationships is important for meaningful assessment of the utility of our method for application in other regions. Inferences on structural relationships are also needed for the design of sampling to avoid biases that might lead to erroneous estimations of forest loss [44]. Indeed, our results showed the benefits of modelling by revealing the unsatisfactory performance of the GFC dataset in areas of small forest plots.

Conclusions
We examined the accuracy of results obtained by using the GFC dataset to identify forest loss in a region with numerous small privately owned forest plots. Generally in such regions, areas of forest loss are also small. Our results indicate that use of the GFC dataset did not lead to accurate detection of areas of forest loss in our study area, and that the small size of those areas was the major factor of the inaccuracy. These results imply that forest conditions in a particular region inevitably affect the accuracy of the use of satellite imagery to detect forest loss. If satellite images are used for this purpose, the causal relationships between forest loss and the particular characteristics of that region must be considered. To support accurate application of satellite imagery for future regional-scale sustainable forest management, further studies are required to examine the causal relationships between accuracy of land-use classification and different forest conditions in various regions.