A comprehensive evaluation of disturbance agent classification approaches: Strengths of ensemble classification, multiple indices, spatio-temporal variables, and direct prediction

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

A comprehensive evaluation of disturbance agent classification approaches: Strengths of ensemble classification, multiple indices, spatio-temporal variables, and direct prediction

T

Katsuto Shimizua,b,1, Tetsuji Otac, , Nobuya Mizouea, Shigejiro Yoshidaa ⁎

a

Faculty of Agriculture, Kyushu University, 744 Motooka, Fukuoka 819-0395, Japan Research Fellow of the Japan Society for the Promotion of Science, Kojimachi, Chiyoda-ku, Tokyo 102-0083, Japan c Institute of Decision Science for a Sustainable Society, Kyushu University, 744 Motooka, Fukuoka 819-0395, Japan b

ARTICLE INFO

ABSTRACT

Keywords: Landsat Time series Tropical seasonal forest Causal agent

Landsat time series images are used for the detection of forest disturbance and the classification of causal agents. Various studies have classified disturbance agents with respect to forest disturbance detected using Landsat time series images. However, the accuracy of the finally classified disturbance agents in different approaches is rarely evaluated. In this study, we investigated the effectiveness of using ensemble classification, and multiple spectral and spatio-temporal information for the accuracy of the classification of disturbance agents in two-stage prediction (i.e., disturbance agents are classified with respect to the detected disturbance) and direct prediction (i.e., disturbance agents are directly classified from Landsat temporal information). Predictor variables were derived from the results of the trajectory-based temporal segmentation of five spectral indices using an annual Landsat time series (2000–2018). We compared six approaches of classifying disturbance agents. For two-stage prediction, we investigated four disturbance detection approaches: threshold-based detection with a single spectral index, random forest (RF) model with a single spectral index, RF model with multiple spectral indices, and RF model with spatio-temporal variables. The detected disturbance pixels were aggregated to disturbance patches and classified into disturbance agents. For direct prediction, two RF models one with only temporal variables and the other with spatio-temporal variables were constructed to classify pixel-based disturbance agents. The overall accuracy of the RF model using spatio-temporal variables for direct prediction was 92.4% and significantly higher than that of the RF model for two-stage prediction (90.9%). The use of an RF model based only on a single spectral index in disturbance detection was not effective for improving accuracy compared with threshold-based detection; however, the use of an RF model based on multiple spectral indices in disturbance detection improved the accuracy of the final classification of disturbance agents. Introducing spatial variables in RF models was effective for improving the overall classification accuracy in pixel-based direct prediction. However, it was not necessary in two-stage prediction because of spatial information contained in the patches. Although a spatially discontinuous appearance was observed for the RF model for directly classifying disturbance agents, this could be an alternative approach to two-stage prediction when considering the relative classification performance and simplicity of implementation.

1. Introduction The terrestrial carbon budget is significantly affected by forest disturbance and recovery, largely because of land use change, deforestation, and regrowth in tropical forests (Pan et al., 2011). In the last decade, tropical forests have been estimated as a net carbon source when considering carbon emissions from natural and anthropogenic

disturbances (Baccini et al., 2017). Quantifying the location and extent of these disturbances potentially provides an understanding of the proximate and underlying causes of forest disturbances. The identification of the proximate causes of the change (i.e., causal agents or disturbance agents), such as logging or road construction, could help to explain the underlying causes that could lead to effective policy actions to reduce carbon emissions from forest disturbances (Finer et al., 2018).

Corresponding author. E-mail address: [email protected] (T. Ota). 1 Present address: Forestry and Forest Products Research Institute, Matsunosato 1, Tsukuba, Ibaraki 305-8687, Japan. ⁎

https://doi.org/10.1016/j.isprsjprs.2019.10.004 Received 4 July 2019; Received in revised form 9 October 2019; Accepted 9 October 2019 0924-2716/ © 2019 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Thus, both accurate monitoring of forest disturbances and identification of causal agents of disturbances in the tropics have become increasingly important. Satellite imagery plays a vital role in investigating past and present forest dynamics in the tropics because historical ground survey data that cover a large area are rarely available. One approach to identify the causal agents of disturbances is a visual interpretation of forest disturbances. De Sy et al. (2015) quantified the drivers of deforestation by visually interpreting high spatial resolution satellite images. However, a visual interpretation highly depends on trained experts and thus is not suitable for wall-to-wall mapping. An alternative approach is the automated classification of disturbances detected by satellite imagery. Since the United States Geological Survey (USGS) made the Landsat archive freely available (Woodcock et al., 2008), many studies have developed automated change detection algorithms to detect forest disturbance using annual or sub annual Landsat time series (Zhu, 2017). A variety of algorithms using Landsat time series are available in the literature, such as Continuous Change Detection and Classification (Zhu and Woodcock, 2014), Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) (Kennedy et al., 2010), Composite2Change (Hermosilla et al., 2015a, 2015b), Breaks For Additive Seasonal and Trend (BFAST) (Verbesselt et al., 2012, 2010) and Vegetation Change Tracker (Huang et al., 2010). Several studies have classified disturbances detected by Landsat time series into causal agents using a machine learning algorithm after aggregating disturbance pixels as disturbance patches (e.g., Huo et al., 2019; Kennedy et al., 2015) (Fig. 1). These approaches that depend on automated disturbance detection and classification have significant benefits, particularly for mapping a large area. Regarding disturbance detection, most change detection algorithms enhance the pattern and variation of the temporal trajectories of a single spectral index to characterize the signs of disturbance (Fig. 1a). However, this approach has two major problems that prevent accurate disturbance detection. First, based on temporal changes in spectral indices, several algorithms such as LandTrendr use a threshold to label disturbance or no disturbance by enhancing the magnitude of the

difference from previous observations or predicted values (i.e., threshold-based disturbance detection). This means that only a change that has a magnitude greater than a predefined threshold is labeled as a disturbance. Because different disturbance agents have different magnitudes of changes in the spectral index, the optimal threshold to detect disturbance varies for each disturbance agent. For example, we have to set a moderate threshold to detect small-scale or subtle disturbances, such as those caused by selective logging, but this leads to false positive detection for observational noise in a time series. Conversely, a strict threshold can reduce false positive detection; however, it causes omissions for small disturbances. Thus, using a single threshold cannot properly discriminate real disturbances from false detection in areas with diverse disturbance agents. To overcome this, several studies have used machine learning algorithms such as ensemble classification by random forests (RFs) (Breiman, 2001) to detect disturbances based on predictor variables derived from a temporal trajectory (Liang et al., 2016; Nguyen et al., 2018). Predictor variables in previous studies included spectral magnitude (e.g., Cohen et al., 2018), change duration, slope (i.e., the ratio of magnitude and duration), spectral value, onset year (e.g., Nguyen et al., 2018), and max/min spectral value and slope of spatially adjacent pixels (e.g., Liang et al., 2016) based on the variety of spectral indices. Second, relying on a single spectral index can cause significant false detection. The signals of different spectral indices associated with disturbance are quite different from each other, particularly for subtle changes (Hislop et al., 2018; Kennedy et al., 2010; Schultz et al., 2016). In this context, the simultaneous use of several spectral indices for disturbance detection might provide a clear advantage for enhancing the ability to separate forest disturbance from noise compared with detection using only a single spectral index. Using a machine learning algorithm, multiple spectral indices can be incorporated for accurate disturbance detection and agent classification. However, the application of a machine learning algorithm with several spectral indices was rarely investigated in previous studies (e.g., Cohen et al., 2018; Schultz et al., 2016), particularly for tropical forests where diverse agents cause different patterns of disturbance events. In addition to the temporal domain of forest disturbance, the spatial

Fig. 1. Schematic flow of (A) two-stage prediction and (B) direct prediction for the classification of disturbance agents. Grid cells indicate disturbance pixels in Landsat images. Different colors represent different disturbance agents. 100

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

pattern is another important factor for characterizing the dynamics of disturbance and causal agents. Because forest disturbance occurs at different spatial scales from the single tree level to the stand level (e.g., fire and clear-cut), it is necessary to monitor disturbance processes as spatio-temporal events (Hughes et al., 2017). However, most algorithms only rely on temporal dynamics and consider each pixel independently to characterize disturbance events (Zhu, 2017) (Fig. 1b). Only a few approaches use spatio-temporal information to detect disturbance (e.g., Hamunyela et al., 2016a, 2016b; Hughes et al., 2017). When only the temporal patterns of the spectral index are used to detect disturbance in each pixel independently, a detected disturbance might exhibit an irrelevant change in the spatial context. Although the classification of causal agents is sometimes conducted based on a patch to incorporate spatial information, which is formed by grouping spatially adjacent disturbance pixels at the same date (Fig. 1c), the detection of change is merely at the pixel-level in this approach. The effectiveness of spatial information for detecting disturbance and classifying causal agents is not known. In previous studies, the location and timing of disturbances were first detected at the pixel level, and then after aggregating disturbance pixels as disturbance patches, disturbance agents were classified only for detected patch-based disturbances using a machine learning algorithm (Fig. 1d) (hereafter referred to as two-stage prediction, Fig. 1A). This approach is frequently applied for tropical, temperate, and boreal forests (e.g., Ahmed et al., 2017; Kennedy et al., 2015; MurilloSandoval et al., 2018; Oeser et al., 2017). Murillo-Sandoval et al. (2018) applied the BFAST algorithm to detect pixel-level disturbances from Landsat time series and classified disturbance agents of aggregated disturbance patches in the Colombian Andes. They demonstrated overall accuracies of 96% and 98% for pixel-level disturbance detection and patch-level agent classification, respectively. For two-stage prediction, accuracy assessments are conducted for both disturbance detection and agent classification. However, the agents of omitted disturbance in disturbance detection are often unknown because disturbance agents are classified only for detected disturbances. The omission level of disturbance may vary substantially depending on causal agents because of the different response of spectral indices caused by different disturbance agents. Consequently, some disturbances caused by specific agents, such as subtle disturbances, might be largely omitted; however, this is not reflected in the accuracy assessments. As Kennedy et al. (2015) discussed, an approach to minimize the omission of disturbance, such as applying a moderate threshold to detect a subtle change, is a possible approach if the following patchbased classification of disturbance agents can remove false positive detection. However, the omission of disturbance cannot be removed in any type of approach. In this context, the omission error of disturbance should be presented for each disturbance agent in the accuracy assessment. One possible approach to assess the performance of disturbance detection and agent classification is to use finally mapped disturbance agents in a pixel-based accuracy assessment because it reflects both the results of disturbance detection and agent classification. By doing this, the omission and commission of each disturbance agent can be assessed. In addition to two-stage prediction, the classification of disturbance agents without labeling disturbance, which means that disturbance agents are directly classified by the temporal patterns of spectral indices using a machine learning algorithm, might be an alternative approach (hereafter referred to as direct prediction, Fig. 1B). Different from twostage prediction, there is no binary labeling of disturbance or no disturbance in this approach. Instead, a machine learning algorithm is used to classify disturbance agents such as logging, dam construction, and shifting cultivation, and no disturbance (Fig. 1e). Because this approach skips disturbance detection, it provides clear insights for mapping disturbance agents. Although direct prediction is easy to implement and might achieve comparable accuracy for classifying disturbance agents, there are no previous studies that have investigated

the relative performance of the classification of disturbance agents using direct prediction. Ideally, disturbance agents should be identified using temporal, spatial, and spectral characteristics. Previous studies have well utilized these characteristics by forming spatial patches that represent primary unit of disturbance. This approach is reasonable for identifying disturbance agents; however, forming spatial patches after disturbance detection requires a great deal of processing time. Other approaches such as the classification of disturbance agents without forming disturbance patches are possible to reduce processing time. However, the role of different characteristics on the accuracy of classification has not been evaluated, and thus possible variables that should be included in other approaches, such as direct prediction, are not known. Investigations on direct and two-stage prediction with different types of variables can inform approaches for accurately classifying disturbance agents in a simple way. In this study, we compared six approaches for the classification of disturbance agents in tropical seasonal forests in Myanmar. Our specific objectives were to investigate the effectiveness on the accuracy of the classification of disturbance agents in terms of the following: (1) an ensemble classification (i.e., RF model) compared with threshold-based disturbance detection, (2) multiple spectral indices compared with a single spectral index, (3) spatio-temporal variables compared with only temporal variables, and (4) direct prediction compared with two-stage prediction. We used LandTrendr temporal segmentation for change characterization, because it has been widely applied for the detection of forest disturbance and has been shown to be capable of addressing many research objectives. We compared the accuracy of the classification of disturbance agents based on finally generated disturbance agent mapping to assess both the results of disturbance detection and agent classification. 2. Study area The study area is located in the Bago mountain range in the central part of Myanmar (Fig. 2). It comprises approximately 2.3 million ha, mainly dominated by mixed deciduous forests. The main species in the forests are teak (Tectona grandis) and pyinkado (Xylia xylocarpa). Mean annual rainfall and temperature in this region are approximately 2000 mm and 27 °C, respectively. The elevation of the study area ranges from 10 to 600 m above sea level with rugged terrain. Historically, there have been several disturbance agents in this region. The main disturbance agents in our study area are logging, dam construction, infrastructure development, establishment of timber plantations and agricultural expansion. Because of the rich environment for timber harvesting, the Myanmar government has been managing forests in the mountainous areas for selective logging operations. Selective logging has been conducted to harvest economically important hardwood species that are larger than the minimum stem size in Reserved Forests, in which access and harvest of forest products by local people has been restricted (Tun et al., 2016). From 2016, the government banned logging operations for 10 years in this region to preserve forest resources in response to forest cover loss. Other than disturbances caused by legal selective logging, forests have been subject to illegal logging (Khai et al., 2016; Mon et al., 2012; Win et al., 2018); however, the details of illegal logging activities across the entire region are unknown because there are no estimates of the extent of the activities. The construction of hydroelectric dams causes large scale forest loss (Prescott et al., 2017). The recent construction of dams in the study area has had a significant impact on the forest condition. Infrastructure development, mainly the construction of roads or settlements, is another important disturbance agent (Ling et al., 2017) especially in the region affected by the development of Nay Pyi Taw, which became the new capital city of Myanmar after 2006. Timber plantations, mainly for teak and other valuable species, have been established in Reserved Forests in the mountain region (Maung and Yamamoto, 2008). 101

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 2. Location of the study area in Myanmar. Natural color Landsat 8 Operational Land Imager images acquired in 2014 were used as the background. The World Borders Dataset provided by Thematic Mapping (http://thematicmapping.org/) was used as the country border.

Agricultural expansion in this study includes shifting cultivation, rubber plantations, and conversions to permanent agriculture land. Shifting cultivation has been practiced by local people mainly in lowland areas near villages. The average size of cleared fields ranges from 1.72 to 2.70 ha with a fallow period from 9 to 12 years in this region (Ei et al., 2017). Although the area affected by shifting cultivation has declined in Myanmar (Li et al., 2018), shifting cultivation is still the main agent of disturbance in the study area. Over past decades, rubber plantations have expanded in Myanmar (Chen et al., 2018; Fox et al., 2014). The majority of rubber plantations have been established in the northern and southern part of the country; however, the Bago region also had approximately 45,000 ha of rubber plantations in 2014 (Kenney-Lazar et al., 2018). Forests near the ridges of villages are sometimes subject to conversion to permanent agricultural land. These diverse agents of disturbance make this study area suitable for investigating the approach for the prediction of causal agents.

were converted from Landsat annual surface reflectance image composites. Then, spatial and temporal variables were calculated using the result of temporal segmentation for each pixel in each year. Based on the calculated variables, classification models were constructed using two-stage and direct prediction. For two-stage prediction, thresholdbased disturbance detection and three RF models were constructed to detect disturbances. Then, disturbance pixels were aggregated to patches for each model based on the spatial adjacency of disturbance pixels in the same year. Predictor variables were calculated for each patch, and RF agent classification models were constructed to classify patchbased disturbance agents. For direct prediction, two RF models to classify pixel-based disturbance agents were constructed with spatial and temporal variables derived from LandTrendr temporal segmentation. Finally, pixel-based accuracy assessments were conducted for each approach using reference samples collected by visual interpretation. 3.1. Data

3. Methods

We used Landsat Collection 1 Surface Reflectance Level-2 products of the Landsat Thematic Mapper (TM), Enhanced Thematic Mapper plus (ETM+), and Operational Land Imager (OLI) images acquired from 1999 to 2018 in the five path/rows of World Reference System 2 (i.e.,

To achieve our research objectives, a total of six approaches for classifying disturbance agents were investigated (Fig. 3). We first ran LandTrendr temporal segmentation for five spectral indices, which 102

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 3. Flowchart of the overall analysis.

resampled and co-registered to the pixels of the Landsat images. The elevation and slope were extracted from the DEM to be included in the RF models.

132/47, 132/48, 133/46, 133/47, and 133/48). Only images that were acquired from October 1st to January 31st with cloud cover of less than 70% were used to avoid leaf off conditions in deciduous forests. A total of 847 Landsat images satisfied the criteria. We downloaded Surface Reflectance products that were processed using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al., 2006) for TM/ETM+ images and with the Landsat Surface Reflectance Code (LaSRC) (Vermote et al., 2016) for OLI images from the USGS Landsat archive (https://earthexplorer.usgs.gov/). This data product contained a pixel quality assessment band created from an Fmask algorithm (Zhu and Woodcock, 2012), which was used for cloud and shadow removal in this study. All images were normalized to the base image in each path/row to minimize the variation of surface reflectance values mainly caused by vegetation seasonality using the iteratively reweighted multivariate alteration detection algorithm (Canty and Nielsen, 2008). After normalization, C-correction (Teillet et al., 1982) was applied to all the images to remove the effects of topography. Then, annual pixelbased image composites were developed based on the set of calculations using the day of year, sensor, opacity, and distance from clouds and shadows in each pixel for a given year, as indicated in White et al. (2014). We used the same criteria as White et al. (2014) to identify the best available pixels. Pixels that were missing after image compositing were filled with the average values of temporally adjacent pixels (i.e., pixels in −1 and +1 year of a given year) in annual image composites. By doing this, gap-free annual Landsat image composites were prepared for LandTrendr temporal segmentation. A digital elevation model (DEM) was used for the calculation of topographic variables in the topographic correction and RF models. We used 1 arc second data from the Shuttle Radar Topography Mission (SRTM) obtained from the USGS archives. The original DEM was

3.2. Temporal segmentation We used LandTrendr (Kennedy et al., 2010) to implement the temporal segmentation of the trajectory-based analysis. Using the temporal change of a spectral index in each pixel, LandTrendr finds vertices of a spectral trajectory and fits several straight-line segments to minimize noise for the detection of both abrupt and gradual disturbances. Vertices indicate the timing of potential vegetation change and temporal segments following these vertices generally represent change condition such as stable, disturbance, and recovery. A set of parameters to control the detection of vertices and the determination of straight-line fitting is required for the algorithm. We used parameters tuned in previous studies for this region (Shimizu et al., 2017). Tasseled Cap Brightness (TCB), Greenness (TCG), Wetness (TCW), and Angle (TCA) that were derived from Tasseled Cap transformation (Crist, 1985; Huang et al., 2002), and Normalized Burn Ratio (NBR) were used as spectral indices to characterize the temporal dynamics. This procedure resulted in temporal segments that represented the timing and duration of potential change on the spectral data for five spectral indices. The temporal segments and their vertices of different spectral indices might indicate the different timing of potential change because LandTrendr independently generated temporal segments for each spectral index. We did not force these vertices to comply with the same timing of one spectral index, such as implemented in (Kennedy et al., 2015). Rather, we preserved these differences to represent the characteristics of each spectral index for disturbance events in prediction. 103

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

Topographic elevation Topographic slope

Mean source spectral index value in pixel windows (3 × 3, 15 × 15, 31 × 31 and 51 × 51) for focal year Number of vertices in pixel windows (3 × 3, 15 × 15, 31 × 31 and 51 × 51) for focal year Mean spectral magnitude in pixel windows (3 × 3, 15 × 15, 31 × 31 and 51 × 51) for focal year TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR

Pixel-based reference samples were collected based on stratified random sampling. Because there was no available map for stratification, we first conducted purposive sampling. Then, we constructed a preliminary pixel-based RF model to map annual disturbance agents to be used for stratification. All extracted predictor variables in Table 1 were used to classify the disturbance agents and no disturbance at all pixels of study area across all years. We defined seven disturbance agents in this study; logging, dam construction, infrastructure development, timber plantations, shifting cultivation, rubber plantations and agricultural expansion (Table 2). Based on the disturbance agent map, we conducted stratified random sampling across the entire study area. The sample size was determined using the equation described in Olofsson et al. (2014). The mapped proportion of classes in the equation was derived from the result of the disturbance agent map. We finally collected a total of 7246 reference samples, which included at least 200 samples for each class. A visual interpretation was conducted for each sample to label the agent class and year using Landsat time series images, temporal trajectories of spectral indices and high spatial resolution images from Google Earth™. For the interpretation of logging activities, we also used the record of selective logging operations by the government, which indicated the year, total number of harvested trees, and rough location of logging activities. After visual interpretation, each sample had the agent class label and disturbance year if the pixel experienced disturbance during study period (i.e., from 2001 to 2017). We allocated 70% of the reference samples as training samples and the remainder as test samples. In each training sample, we only used one annual observation from entire study period for constructing models because we aimed at classifying annual disturbance and disturbance agents. When training samples had experienced disturbances,

Elevation Slope Topographic variables

1 1

Mean spectral value in window Number of vertices in window Mean spectral magnitude in window Spatial variables

4 4 4

Source spectral index value for pre, focal and post year Fitted spectral index value for pre, focal and post year Magnitude of source spectral index value for pre and post year Magnitude of fitted spectral index value for pre and post year Binary variable of the presence or absence of vertex for focal year Source spectral value (pre, focal, post) Fitted spectral value (pre, focal, post) Source spectral magnitude (pre, post) Fitted spectral magnitude (pre, post) Presence of vertex (focal) Temporal variables

3 3 2 2 1

TCB, TCB, TCB, TCB, TCB,

TCG, TCG, TCG, TCG, TCG,

TCW, TCW, TCW, TCW, TCW,

TCA, TCA, TCA, TCA, TCA,

NBR NBR NBR NBR NBR

Description SI

The predictor variables were calculated from the results of LandTrendr segmentation for each spectral index to represent the condition of a target pixel in the focal year (Table 1). These variables were categorized as temporal and spatial variables. For calculating the temporal variables, we set a three-year temporal window to extract pre and post information of the focal year (Fig. 4). The source and fitted spectral index values were derived from the temporal window. By differentiating the spectral index values of the pre, focal, and post years, the magnitude of the pre and post years were calculated both for the source and fitted spectral index. The binary presence or absence of the vertex derived from LandTrendr in the focal year was also extracted. The spatial variables were derived based on spatial moving windows ranging from 3 × 3 to 51 × 51 pixels. The spatial variables were calculated using the results of temporal segmentation in each pixel within the spatial moving windows. The total number of pixels that had a vertex in the focal year in the spatial moving windows was counted. Mean spectral values and the magnitude of spectral values in the pixels were also calculated. All temporal and spatial variables were generated for every pixel in each year except the first and last year because several temporal variables could not be calculated for those years. This indicates that predictor variables were generated for all years, regardless of the results of LandTrendr temporal segmentation. We extracted source spectral index value and magnitude with the intention of detecting subtle disturbances that were missed in LandTrendr temporal segmentation even though there was the risk of including noise. Topographic variables (i.e., elevation and slope) were also included as predictor variables. As a result, a total of 117 variables (i.e., 23 variables for each of five spectral indices and two topographic variables) were extracted. The specific combination of these variables was used to construct the disturbance detection models for two-stage prediction and RF classification models for direct prediction. 3.4. Reference sample collection

Variable

Variables for each SI

3.3. Pixel-based variable extraction

Class

Table 1 Predictor variables used in the random forest (RF) models. Spatial and temporal variables were calculated for each spectral index (SI).

K. Shimizu, et al.

104

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 4. Temporal variables derived from the results of LandTrendr temporal segmentation (black lines and dots) and the source spectral values (white dots). Within a three-year temporal window, pre, focal, and post fitted and source spectral values were extracted to calculate the variables.

Table 2 Disturbance agents and no disturbance defined in this study. Agent class

Description

Logging

Stand replacing forest disturbances caused by logging without transformation to other land use types. This also includes disturbance caused by logging activities, such as log decks and logging roads Construction of dam and reservoirs in forests Land cover change from forest to infrastructure caused by the construction of settlements and roads Establishment of timber plantations mainly for teak and other valuable species Small-scale disturbances caused by shifting cultivation. This excludes conversion to permanent agricultural land Stand replacing forest conversion followed by the development of rubber plantations Agricultural land expansion permanently converted from forests. This excludes disturbance caused by shifting cultivation and rubber plantations No forest disturbances, including other changes in other land covers, such as agricultural land, shrub and water

Dam construction Infrastructure development Timber plantations Shifting cultivation Rubber plantations Agricultural expansion No disturbance

observations in the year of disturbance were selected. When training samples had no disturbance, observations were randomly selected (i.e., one year from 2001 to 2017 were randomly selected). By doing this procedure, annual observations with reference labels were extracted for training models. For direct prediction, these training samples were used to construct RF models (Fig. 5a). For two-stage prediction, however, the training samples were subdivided into disturbance detection (Fig. 5b-1) and agent classification (Fig. 5b-2) because the training samples were required to construct models in both stages. Although it was a possible option to use all training samples both in disturbance detection and agent classification, we subdivided training samples in two-stage prediction to ensure enough training samples for no disturbance in agent classification, especially for the approaches that used RF in disturbance detection. The training samples were relabeled as disturbance or no disturbance in disturbance detection of two-stage prediction. Patchbased training samples were generated to correspond with patches in RF agent classification in two-stage prediction. Disturbance agents were labeled to patches that overlapped with training samples, and these patches were used as training samples to construct agent classification models. The size, shape and location of mapped patches in disturbance detection were different for each model, and thus those of the patchbased training samples were different for each model.

3.5. Two-stage prediction: pixel-based disturbance detection For two-stage prediction, disturbance detection was conducted using (1) a threshold to identify disturbances based on TCW temporal segmentation (TCW threshold), (2) an RF model using temporal variables derived from TCW temporal segmentation and topographic variables (TCW RF temporal), (3) an RF model using temporal variables derived from temporal segmentation of all spectral indices and topographic variables (RF temporal), and (4) an RF model using both temporal and spatial variables derived from the temporal segmentation of all spectral indices and topographic variables (RF spatio-temporal). The TCW was selected to detect disturbances in the threshold-based detection and the corresponding RF model because of its sensitivity to disturbance events. In the threshold-based disturbance detection (i.e., TCW threshold), pixels that had a vertex and a fitted spectral magnitude greater than a predefined threshold were labeled as disturbance. The threshold was determined by adjusting potential thresholds to maximize the overall accuracy of disturbance detection using a total of 2536 training samples. In the corresponding approach (i.e., TCW RF temporal), an RF model was used to detect disturbances as a replacement for the threshold. Disturbances were detected by this RF model using temporal variables derived from the results of temporal segmentation in Table 1; however, the timing of disturbances was controlled by the 105

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 5. Reference samples allocated to the training and test samples for direct prediction and two-stage prediction.

vertices of TCW temporal segmentation. Conversely, there was no restriction controlled by the vertices of temporal segmentation in the two latter approaches (i.e., RF temporal and RF spatio-temporal). All pixels in all years were subject to disturbance detection by these two RF models. In the RF (temporal) model, the temporal and topographic variables in Table 1 were used. Conversely, all predictor variables in Table 1 were used in the RF (spatio-temporal) model. The RF models in the three disturbance detection approaches were developed to detect disturbances using training samples. Each RF model was tuned using 10-fold cross validation with 500 trees using the package “randomForest” (Liaw and Wiener, 2002) and “caret” (Kuhn, 2008) in R, which is software for statistical computing (R Core Team, 2016). 3.6. Two-stage prediction: patch-based agent classification Fig. 6. Patch-based temporal variables for constructing the RF models.

We mapped and grouped disturbance pixels in each disturbance detection model. The spatially adjacent disturbance pixels in the same year were aggregated as patches and assumed to be caused by the same disturbance agent. Spatial adjacency was defined by an 8-neighbor connectivity rule. To remove erroneous detection, only disturbance patches that were greater than three pixels were retained. The disturbance agents were assigned to these patches by building RF classification models with predictor variables defined at patch level. After patch aggregation, predictor variables were extracted for all patches in the study area. The predictor variables were defined as temporal, shape, and topographic variables. Temporal variables were extracted from pre, focal, and post event segments to represent spectral conditions of pre and post disturbance (Fig. 6). We extracted variables based on segments rather than year following the approach in previous studies (e.g., Kennedy et al., 2015; Nguyen et al., 2018). Duration was defined as the average year of segments in each patch. The magnitude was calculated by averaging the difference of spectral index values of focal segments in each patch. Pre and post segment variables were calculated by the average and standard deviation of the spectral index values of a patch. These temporal variables were calculated for TCB,

TCG, TCW, TCW, and NBR. The topographic variables were characterized by the average elevation and slope in each patch. The shape of each patch was characterized using area, perimeter, and compactness. In total, 45 predictor variables were extracted for each patch of four disturbance detection approaches (Table 3). The classification of the patch-based disturbance agents was performed by constructing RF models for each approach. The predictor variables used in the RF agent classification model were the same for each model; the only difference was the disturbance patches that resulted from each disturbance detection approach. The disturbance agents and no disturbance were the same as those in Table 2. No disturbance class was also considered to remove false disturbance detection. The patch-based training samples were used to build the RF models. Each RF model was tuned using 10-fold cross validation with 500 trees in R. The constructed RF models were applied to all disturbance patches in the study area to map disturbance agents in each year. As a result, disturbance agent maps were generated for a total of four approaches. 106

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Table 3 Predictor variables for patch-based RF models to classify disturbance agents. Class

Variable

Spectral index

Description

Focal event segment

Duration Magnitude

TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR

Average duration year of the focal segment in a patch Average magnitude of the focal segment in a patch

Pre-event segment

Duration Average Standard deviation

TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR

Average duration year of the pre-event segment in a patch Average value of the pre-event segment in a patch Standard deviation of the value of the pre-event segment in a patch

Post-event segment

Duration Average Standard deviation

TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR TCB, TCG, TCW, TCA, NBR

Average duration year of the post-event segment in a patch Average value of the post-event segment in a patch Standard deviation of the value of the post-event segment in a patch

Shape

Area Perimeter Compactness

Area of a patch Perimeter of a patch Perimeter/area2

Topography

Elevation Slope

Topographic elevation Topographic slope

prediction. We compared the classification performance between twostage and direct prediction using the combination of (5) RF (spatiotemporal) in two-stage prediction and RF (spatio-temporal) in direct prediction. To investigate the effectiveness of LandTrendr variables, we also compared RF (temporal) and RF (spatio-temporal) with and without LandTrendr variables in direct prediction.

3.7. Direct prediction: pixel-based classification of disturbance agents For direct prediction, pixel-based disturbance agents were directly classified by (1) an RF model using temporal variables derived from temporal segmentation of all spectral indices and topographic variables (RF temporal) and (2) RF model using both temporal and spatial variables derived from temporal segmentation of all spectral indices and topographic variables (RF spatio-temporal). These RF models directly predict disturbance agents using pixel-based predictor variables. To construct the RF classification models, disturbance agents and no disturbance were defined as in Table 2. The RF models were tuned using 10-fold cross validation with 500 trees in R. A total of 5072 training samples were used to tune the RF models. The constructed RF models were applied to all pixels in the study area to map disturbance agents. To investigate the effectiveness of using LandTrendr temporal segmentation, we also constructed models without LandTrendr variables for both RF temporal and RF spatio-temporal. If accuracies of these models were comparable to models with LandTrendr variables, we can reduce time for implementing LandTrendr temporal segmentation.

4. Results The results of the accuracy assessment of the six approaches for classifying disturbance agents are shown in Table 4. For two-stage prediction, the OA of TCW threshold, TCW RF (temporal), RF (temporal), and RF (spatio-temporal) were 84.8% ( ± 1.5%), 86.4% ( ± 1.4%), 90.9% ( ± 1.1%), and 90.9% ( ± 1.1%), respectively, based on error-adjusted matrices. For direct prediction, the OA of RF (temporal) and RF (spatio-temporal) were 89.0% ( ± 1.2%) and 92.4% ( ± 1.1%), respectively. The PA and UA of each disturbance agent in each approach are shown in Fig. 7. The PA of disturbance agents of TCW threshold and TCW RF (temporal) for two-stage prediction were generally lower than those of other approaches (Fig. 7). Particularly for logging, dam construction, infrastructure development, and shifting cultivation, the differences of PA compared with other approaches were more than 17%, which indicates that there was a large omission in these two approaches. The PA of infrastructure development (77%), timber plantation (74%), and shifting cultivation (75%) for RF (spatio-temporal) for two-stage prediction were the highest among the other approaches. The PAs of RF (temporal) for two-stage prediction were lower than those of RF (spatio-temporal) for two-stage prediction, except for agricultural expansion. The PA of logging (65%) and rubber plantations (62%) of RF (spatio-temporal) for direct prediction were the highest among the other models. The UAs of RF (temporal) for direct prediction were lower than those of RF (spatio-temporal) for direct prediction. Among

3.8. Accuracy assessment The accuracies of classifying disturbance agents both in two-stage and direct prediction were evaluated using 2174 pixel-based test samples. Mapped disturbance agents and no disturbance in each approach were used for the accuracy assessment to account for both results of disturbance detection and agent classification. The accuracies were assessed in terms of overall accuracy (OA) based on the sample countbased confusion matrices and error-adjusted matrices proposed by Olofsson et al. (2014). The producer’s accuracy (PA) and the user’s accuracy (UA) of each class were also calculated based on the erroradjusted matrices. Furthermore, McNemar’s test was applied to investigate whether there were significant differences between classifications using different approaches. McNemar’s test is a non-parametric test based on a binary confusion matrix and is suitable for assessing the classification performance with the same samples (Foody, 2004). To assess the effectiveness of using RF models, multiple spectral indices and spatial variables on the accuracy of the classification of disturbance agents in two-stage and direct prediction, we compared five combinations of approaches. We compared (1) TCW threshold and TCW RF (temporal) in two-stage prediction to reveal the effectiveness of the RF model compared with threshold-based detection. The effectiveness of multiple spectral indices was investigated by comparing (2) TCW RF (temporal) and RF (temporal) in two-stage prediction. The gain of using spatial variables was investigated in (3) RF (temporal) and RF (spatio-temporal) in two-stage prediction, and (4) RF (temporal) and RF (spatio-temporal) in direct

Table 4 Results of the classification of disturbance agents for each approach. The overall accuracies (OA) calculated by a sample count-based confusion matrix and an error-adjusted confusion matrix are shown. The 95% confidence interval is also shown with error-adjusted OA.

107

Prediction

Models

Count-based OA (%)

Error-adjusted OA (%)

Two-stage Two-stage Two-stage Two-stage Direct Direct

TCW threshold TCW RF (temporal) RF (temporal) RF (spatio-temporal) RF (temporal) RF (spatio-temporal)

83.9 84.6 88.3 88.2 85.6 89.9

84.8 86.4 90.9 90.9 89.0 92.4

( ± 1.5) ( ± 1.4) ( ± 1.1) ( ± 1.1) ( ± 1.2) ( ± 1.1)

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 7. (a) Producer’s accuracy and (b) user’s accuracy of each class in the error-adjusted matrix of each approach.

magnitude of TCA and TCW were the most important for temporal and spatial variables (Fig. 9). Other spectral indices used to calculate the magnitude of change also had relatively high importance in the RF model. The mapped disturbance agents and year of the RF (spatiotemporal) are shown in Fig. 10.

disturbance agents, except no disturbance, the UA of dam construction was the highest for all approaches. The examples of disturbance detection and agent classification mapping for each approach are shown in Fig. 8. Using multiple spectral indices in the RF model, the omission of disturbance detection was visually reduced for two-stage prediction. By including spatial variables in the RF model, salt-and-pepper effects were reduced in disturbance detection for two-stage prediction, although the final OA of the classification of disturbance agents was almost the same. Mapped disturbance agents for direct prediction showed that some adjacent disturbance pixels were classified as different disturbance agents despite their spatial adjacency. Although the OA of RF (spatio-temporal) for two-stage prediction was lower than that of RF (spatio-temporal) for direct prediction, the mapped disturbance agents were more consistent in the spatial context. The comparison results for classification using McNemar’s test are shown in Table 5. There was a significant difference between the approaches in the following combinations (p = 0.05): TCW RF (temporal) and RF (temporal) for two-stage prediction, RF (temporal) and RF (spatio-temporal) for direct prediction, and RF (spatio-temporal) for two-stage prediction and RF (spatio-temporal) for direct prediction. There was no statistically significant difference between TCW threshold and TCW RF (temporal) for two-stage prediction, although the OAs were slightly different from each other (84.8% and 86.4%, respectively). The comparison of RF (temporal) and RF (spatio-temporal) for two-stage prediction showed no statistically significant difference between these approaches. The OAs of approaches without LandTrendr variables for RF (temporal) and RF (spatio-temporal) were 86.3% ( ± 1.3%) and 89.0% ( ± 1.3%), respectively. These OAs were lower than those of approaches with LandTrendr variables. McNemar’s test showed that there were significant differences between RF (temporal) with and without LandTrendr variables (p value = 0.01) and RF (spatio-temporal) with and without LandTrendr variables (p value < 0.001) in direct prediction. The highest OA was observed for the RF (spatio-temporal) for direct prediction. The RF model showed that the variables associated with the

5. Discussion The accurate detection of forest disturbance and classification of causal agents are vital for understanding the process of forest disturbance. Previous studies have mainly identified disturbance agents by detecting pixel-based disturbance and then classifying patch-based disturbance agents. The use of multiple spectral indices, spatial and temporal variables might improve the accuracy of identifying direct causes of disturbance, when considering the dynamics of disturbance. Furthermore, other prediction approaches that can reduce processing time with reasonable accuracies need to be investigated for operational monitoring. In this study, we investigated the effectiveness of using an ensemble classification, multiple spectral indices, and spatio-temporal variables, in the direct and two-stage classifying approaches for the classification of disturbance agents in tropical seasonal forests in Myanmar. We used spatial and temporal variables that were derived from the results of Landsat trajectory-based multiple spectral analysis. 5.1. Effectiveness of the ensemble approach, and multiple spectral and spatio-temporal information The relative performance of using ensemble classification (i.e., RF) compared with threshold-based disturbance detection was evaluated with the TCW threshold and TCW RF (temporal) approaches. TCW RF (temporal) had a slightly higher OA (86.4%) of the classification of disturbance agents than the TCW threshold (84.8%); however, there was no statistically significant difference in McNemar’s test (p = 0.05). Although RF can use multiple factors of change by introducing multiple variables to detect disturbance, it was not effective for the final classification results of disturbance agents in this study. In TCW RF (temporal), the timing of a disturbance occurrence was determined by TCW 108

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 9. Top 20 important variables of 117 predictor variables in RF (spatiotemporal) for direct prediction.

improving the accuracy of the final classification of disturbance agents, an RF model using multiple spectral indices can improve the accuracy of the classification of disturbance agents. Previous studies have also emphasized the utility of RF models in terms of detecting disturbance (e.g., Nguyen et al. (2018); Schultz et al. (2018)). The gain of multiple spectral indices in the RF model is supported by the different responses of spectral indices at disturbance events (Hislop et al., 2018; Schultz et al., 2016). Schultz et al. (2016) showed that the fusion of spectral indices improves disturbance detection by enhancing their complementary radiometric information, even though spectral indices individually perform poorly. The higher accuracy of RF (temporal) in two-stage prediction in this study revealed the importance of using multiple spectral indices for the classification of disturbance agents. The simultaneous use of spatial and temporal variables significantly enhanced the accuracy of the classification of disturbance agents in direct prediction but not in two-stage prediction. In direct prediction, spatial information was not included in the RF (temporal). The spatial and temporal variables can be complementary information for classifying causal agents, and this resulted in a higher OA in the RF (spatio-temporal) model. Conversely, in two-stage prediction, spatial information was included in the patch-based agent classification as shape-related variables (e.g., area, perimeter, and compactness), even if the spatial variables were not used in disturbance detection. The spatial variables in pixel-based disturbance detection might improve accuracy when considering visually reduced salt-and-pepper effects in the RF (spatio-temporal) approach (Fig. 8); however, even without spatial variables in disturbance detection, the use of shape-related variables in patch-based classification and the removal of small patches in the subsequent procedure made the accuracy of this approach comparable to that of the RF (spatio-temporal). These differences might affect the statistical difference of the accuracy of causal agent classification with spatial variables for two-stage and direct prediction. In this study, we only assessed the overall results of the classification of disturbance agents, which involved

Fig. 8. Examples of disturbance detection and agent classification in 2005 for (1) the two-stage TCW threshold, (2) two-stage TCW RF (temporal), (3) twostage RF (temporal), (4) two-stage RF (spatio-temporal), (5) direct RF (temporal), and (6) direct RF (spatio-temporal). There is no mapping of disturbance detection for direct prediction.

temporal segmentation in the same manner as the TCW threshold. This means that the vertex of the segment was required to identify the disturbance. Conversely, when the predictor variables derived from LandTrendr temporal segmentation of multiple spectral indices were used in the RF model and the timing of the disturbance was not restricted by the temporal segmentation of single spectral index, the OA of the RF model (i.e., RF temporal) was significantly higher than that of TCW RF (temporal). Although the use of the RF model as a replacement for the threshold in disturbance detection was not effective for

Table 5 Comparison of approaches of classifying disturbance agents using McNemar’s test. Models

Chi-squared

p value

Two-stage TCW threshold vs two-stage TCW RF (temporal) Two-stage TCW RF (temporal) vs two-stage RF (temporal) Two-stage RF (temporal) vs two-stage RF (spatio-temporal) Direct RF (temporal) vs direct RF (spatio-temporal) Two-stage RF (spatio-temporal) vs direct RF (spatio-temporal)

1.53 31.68 0.02 40.69 7.24

0.216 < 0.001 0.876 < 0.001 0.007

109

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Fig. 10. Disturbance agents map of RF (spatio-temporal) in direct prediction from 2001 to 2017 for the entire study area with the subsets of (a) disturbance agents and (b) disturbance year.

and infrastructure development). This might be because temporal variables largely explained difference from other disturbance agents even without spatial variables. Furthermore, even though the accuracies improved using both spatial and temporal variables in direct prediction, several disturbance agents such as agricultural expansion had relatively low accuracies in the RF (spatio-temporal). Because agricultural expansion in this study occurred on forest ridges, which generally had lower forest cover, omission errors might have occurred, thereby recording no disturbance pixels where disturbances had occurred. The variable importance of RF (spatio-temporal) in direct prediction showed that spatial, temporal, and topographic variables were useful to classify disturbance agents. Source magnitude variables were the most important to distinguish disturbance agents, indicating that source spectral variables as well as fitted spectral variables were useful for classification. Spatial variables, such as mean magnitude of adjacent pixels, were also important factors. The number of vertices in spatial moving indicates potential disturbance around the pixel because vertices are produced if LandTrendr detects disturbances. A small-scale disturbance generally has a few vertices within the window; however, a

both disturbance detection and agent classification. Thus, the effectiveness of spatial information in disturbance detection was not revealed in the results. Our results merely indicated that spatial information should be included both in the two-stage and direct prediction, either as pixelbased or patch-based variables in RF models. The spatial and temporal variables affected the accuracy of the overall classification of causal agents. Despite this, the gain of incorporating spatial or temporal information was different for each disturbance agent. For example, the RF (spatio-temporal) model both in two-stage and direct prediction had a higher PA for timber plantation compared with RF (temporal). Because each disturbance agent has a unique spatial extent (e.g., area) and temporal dynamics (e.g., pre and post condition of disturbance) (McDowell et al., 2015), the information required to discriminate one specific agent from another might be different. Overall, our study indicated that spatial variables improved the PA (i.e., reduced omission errors) of disturbance agents that had unique spatial characteristics (e.g., logging, timber plantations, shifting cultivation, and rubber plantations); however, this was not applicable for disturbance agents that caused land cover change (e.g., dam construction 110

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

large-scale disturbance has many vertices within the window. This difference might be one reason for the high importance of the number of vertices in spatial moving. Although variables related to spatial patches cannot be derived in direct prediction, spatial variables served as proxies for these variables. Topographic elevation was an important variable. Because several agents, such as logging, only occurred in mountainous regions, topographic variables were informative to distinguish disturbance agents. Variable importance may vary depending on study area and disturbance agents because the factor affecting the occurrence of disturbance highly depends on situation in each region. LandTrendr temporal segmentation was important to improve the accuracy of agent classification. This is not a surprising result because LandTrendr can characterize the temporal change of a spectral index and that usually represents disturbance events. Even though implementing LandTrendr requires time for processing, it is recommended to include variables derived from LandTrendr temporal segmentation for the better classification of disturbance agents.

prediction. Even though the mapping appearance is better, the time required to form disturbance patches, map disturbance, and construct RF models has to be considered in two-stage prediction. In this context, direct prediction using spatio-temporal variables in an RF model is a possible alternative approach for operational monitoring. One limitation of direct prediction is the lack of disturbance mapping. This deficiency might be unsuitable for some analyses that use the timing and location of disturbances. Despite this, if the final objective is the prediction of disturbance agents, then the direct prediction approach in this study reduces the processing time and improves the overall accuracy of mapping disturbance agents. Further studies should focus on the applicability of these approaches to other study areas, particularly for other climate regions that have different types of disturbance agents. 6. Conclusions This study investigated the effectiveness of ensemble classification, and multiple spectral and spatio-temporal information for the accuracy of the classification of disturbance agents in patch-based two-stage prediction and pixel-based direct prediction. Spatio-temporal variables were derived from the results of LandTrendr temporal segmentation and used for disturbance detection and agent classification. The results presented in this study indicated that the use of an RF model based only on TCW temporal segmentation in disturbance detection as a replacement of threshold-based detection was not effective for improving the overall accuracy of mapping disturbance agents. Conversely, the use of an RF model based on temporal segmentation of multiple spectral indices in disturbance detection improved the accuracy of the final classification of disturbance agents. Spatial variables were effective for improving the overall classification accuracy in pixel-based direct prediction; however, it was not necessary in patch-based two-stage prediction because the patch statistics already contained spatial information. Furthermore, the pixel-based RF model for directly classifying disturbance agents can be an alternative approach for operational uses because of the relative performance and reduced processing time compared with patch-based two-stage prediction, which is a frequently used approach in other studies.

5.2. Accuracy of the classification for direct and two-stage prediction The direct prediction of disturbance agents investigated in this study showed significantly higher OA compared with two-stage prediction. Despite the statistically significant difference in McNemar’s test, the difference of OA for RF (spatio-temporal) between two-stage and direct prediction was only 1.5% (90.9% and 92.4%, respectively). Although the accuracy of the classification for direct prediction was high, the mapped disturbance agents presented a spatially inconsistent appearance. The location of the disturbance has possibilities to across sub-pixel because disturbance agents were mapped at 0.09 ha (i.e., 30 × 30 m); however, the occurrence of several disturbance agents in a short distance, as in this study, does not reflect a real scenario on the ground. One possible solution to avoid inconsistent appearance in direct prediction is to form patches and apply a median filter to each patch after mapping disturbance agents. Although this requires additional post-processing, predicted disturbance agents can be represented in spatially consistent manner. Another possible solution might be the use of object-based image analysis before running LandTrendr. Although an object-based image analysis using Landsat time series needs additional investigations on variables and spatial segmentation for accurate prediction (e.g., Hughes et al., 2017), this approach has a possibility to cope with spatial discontinuity. The classification results produced by two-stage prediction provided a visually better appearance and contiguous mapping because the disturbance agents were finally mapped based on disturbance patches. Previous studies have indicated that object-based and pixel-based classification using machine learning algorithms can achieve similar overall accuracy (Duro et al., 2012). Although the classification approach of patch-based classification in this study is not the same as object-based classification in previous studies, our results also showed a similar accuracy of the classification of disturbance agents for two-stage and direct prediction. The accuracy of the classification of disturbance agents was assessed using pixel-based reference samples in this study. In previous studies, accuracy assessments were conducted both in pixel-based disturbance detection and patch-based agent classification. Even if the accuracy of a specific agent class is high in patch-based classification, the omission of disturbance in the pixel-based accuracy assessment has to be considered to understand the final mapping accuracy. Furthermore, the omission of disturbance might be different for each disturbance agent, although the omission of each disturbance agent is not typically revealed in the accuracy assessment. The accuracy assessment of the classification in this study evaluated overall classification results. The omission of disturbance in each disturbance agent was clearly revealed in the accuracy, although the accuracy of disturbance detection and agent classification in each step was not known. From an operational point of view, two-stage prediction requires a great deal of time to generate mapping results compared with direct

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This study was supported by the Japan Society for the Promotion of Science KAKENHI Grant Number JP18J12647. The work of Katsuto Shimizu was supported by Research Fellowships of Japan Society for the Promotion of Science for Young Scientists. We thank Leonie Seabrook, PhD, and Maxine Garcia, PhD, from Edanz Group (www. edanzediting.com/ac) for editing a draft of this manuscript. References Ahmed, O.S., Wulder, M.A., White, J.C., Hermosilla, T., Coops, N.C., Franklin, S.E., 2017. Classification of annual non-stand replacing boreal forest change in Canada using Landsat time series: a case study in northern Ontario. Remote Sens. Lett. 8, 29–37. https://doi.org/10.1080/2150704X.2016.1233371. Baccini, A., Walker, W., Carvalho, L., Farina, M., Sulla-Menashe, D., Houghton, R.A., 2017. Tropical forests are a net carbon source based on aboveground measurements of gain and loss. Science 358, 230–234. https://doi.org/10.1126/science.aam5962. Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. Canty, M.J., Nielsen, A.A., 2008. Automatic radiometric normalization of multitemporal satellite imagery with the iteratively re-weighted MAD transformation. Remote Sens. Environ. 112, 1025–1036. https://doi.org/10.1016/j.rse.2007.07.013. Chen, G., Thill, J.-C., Anantsuksomsri, S., Tontisirin, N., Tao, R., 2018. Stand age estimation of rubber (Hevea brasiliensis) plantations using an integrated pixel- and objectbased tree growth model and annual Landsat time series. ISPRS J. Photogramm.

111

ISPRS Journal of Photogrammetry and Remote Sensing 158 (2019) 99–112

K. Shimizu, et al.

Maung, T.M., Yamamoto, M., 2008. Exploring the socio-economic situation of plantation villagers: A case study in Myanmar Bago Yoma. Small-scale For. 7, 29–48. https:// doi.org/10.1007/s11842-008-9039-1. McDowell, N.G., Coops, N.C., Beck, P.S.A., Chambers, J.Q., Gangodagamage, C., Hicke, J.A., Huang, C., Kennedy, R., Krofcheck, D.J., Litvak, M., Meddens, A.J.H., Muss, J., Negrón-Juarez, R., Peng, C., Schwantes, A.M., Swenson, J.J., Vernon, L.J., Williams, A.P., Xu, C., Zhao, M., Running, S.W., Allen, C.D., 2015. Global satellite monitoring of climate-induced vegetation disturbances. Trends Plant Sci. 20, 114–123. https://doi. org/10.1016/j.tplants.2014.10.008. Mon, M.S., Mizoue, N., Htun, N.Z., Kajisa, T., Yoshida, S., 2012. Factors affecting deforestation and forest degradation in selectively logged production forest: a case study in Myanmar. For. Ecol. Manage. 267, 190–198. https://doi.org/10.1016/j. foreco.2011.11.036. Murillo-Sandoval, P.J., Hilker, T., Krawchuk, M.A., Van Den Hoek, J., 2018. Detecting and attributing drivers of forest disturbance in the Colombian Andes using Landsat time-series. Forests 9, 269. Nguyen, T.H., Jones, S.D., Soto-Berelov, M., Haywood, A., Hislop, S., 2018. A spatial and temporal analysis of forest dynamics using Landsat time-series. Remote Sens. Environ. 217, 461–475. https://doi.org/10.1016/J.RSE.2018.08.028. Oeser, J., Pflugmacher, D., Senf, C., Heurich, M., Hostert, P., 2017. Using intra-annual Landsat time series for attributing forest disturbance agents in Central Europe. Forests 8, 251. https://doi.org/10.3390/f8070251. Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E., Wulder, M.A., 2014. Good practices for estimating area and assessing accuracy of land change. Remote Sens. Environ. 148, 42–57. https://doi.org/10.1016/j.rse.2014.02.015. Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L., Shvidenko, A., Lewis, S.L., Canadell, J.G., Ciais, P., Jackson, R.B., Pacala, S.W., McGuire, A.D., Piao, S., Rautiainen, A., Sitch, S., Hayes, D., 2011. A large and persistent carbon sink in the world’s forests. Science 333, 988–993. https://doi.org/10. 1126/science.1201609. Prescott, G.W., Sutherland, W.J., Aguirre, D., Baird, M., Bowman, V., Brunner, J., Connette, G.M., Cosier, M., Dapice, D., Alban, J.D., Diment, A., Fogerite, J., Fox, J., Hlaing, W., Htun, S., Hurd, J., LaJeunesse Connette, K., Lasmana, F., Lim, C.L., Lynam, A., Maung, A.C., McCarron, B., McCarthy, J.F., McShea, W.J., Momberg, F., Mon, M.S., Myint, T., Oberndorf, R., Oo, T.N., Phelps, J., Rao, M., Schmidt-Vogt, D., Speechly, H., Springate-Baginski, O., Steinmetz, R., Talbott, K., Than, M.M., Thaung, T.L., Thawng, S.C., Thein, K.M., Thein, S., Tizard, R., Whitten, T., Williams, G., Wilson, T., Woods, K., Ziegler, A.D., Zrust, M., Webb, E.L., 2017. Political transition and emergent forest-conservation issues in Myanmar. Conserv. Biol. 31, 1257–1270. https://doi.org/10.1111/cobi.13021. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Found. Stat. Comput., Vienna. Schultz, M., Clevers, J.G.P.W., Carter, S., Verbesselt, J., Avitabile, V., Quang, H.V., Herold, M., 2016. Performance of vegetation indices from Landsat time series in deforestation monitoring. Int. J. Appl. Earth Obs. Geoinf. 52, 318–327. https://doi. org/10.1016/j.jag.2016.06.020. Schultz, M., Shapiro, A., Clevers, J., Beech, C., Herold, M., 2018. Forest cover and vegetation degradation detection in the Kavango Zambezi Transfrontier Conservation Area using BFAST monitor. Remote Sens. 10, 1850. https://doi.org/10.3390/ rs10111850. Shimizu, K., Ponce-Hernandez, R., Ahmed, O.S., Ota, T., Win, Z.C., Mizoue, N., Yoshida, S., 2017. Using Landsat time series imagery to detect forest disturbance in selectively logged tropical forests in Myanmar. Can. J. For. Res. 47, 289–296. https://doi.org/ 10.1139/cjfr-2016-0244. Teillet, P.M., Guindon, B., Goodenough, D.G., 1982. On the slope-aspect correction of multispectral scanner data. Can. J. Remote Sens. 8, 84–106. Tun, K.S.W., Stefano, J. Di, Volkova, L., 2016. Forest management influences aboveground carbon and tree species diversity in Myanmar’s mixed deciduous forests. Forests 7, 217. https://doi.org/10.3390/f7100217. Verbesselt, J., Hyndman, R., Zeileis, A., Culvenor, D., 2010. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 114, 2970–2980. https://doi.org/10.1016/j.rse.2010.08.003. Verbesselt, J., Zeileis, A., Herold, M., 2012. Near real-time disturbance detection using satellite image time series. Remote Sens. Environ. 123, 98–108. https://doi.org/10. 1016/j.rse.2012.02.022. Vermote, E., Justice, C., Claverie, M., Franch, B., 2016. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 185, 46–56. https://doi.org/10.1016/j.rse.2016.04.008. White, J.C., Wulder, M.A., Hobart, G.W., Luther, J.E., Hermosilla, T., Griffiths, P., Coops, N.C., Hall, R.J., Hostert, P., Dyk, A., Guindon, L., 2014. Pixel-based image compositing for large-area dense time series applications and science. Can. J. Remote Sens. 40, 192–212. https://doi.org/10.1080/07038992.2014.945827. Win, Z.C., Mizoue, N., Ota, T., Wang, G., Innes, J.L., Kajisa, T., Yoshida, S., 2018. Spatial and temporal patterns of illegal logging in selectively logged production forest: a case study in Yedashe. Myanmar. J. For. Plan. 23, 15–25. https://doi.org/10.20659/jfp. 23.2_15. Woodcock, C.E., Allen, R., Anderson, M., Belward, A.S., Bindschadler, R., Cohen, W.B., Gao, F., Goward, S.N., Helder, D., Helmer, E., Nemani, R., Oreopoulos, L., Schott, J., Thenkabail, P.S., Vermote, E.F., Vogelmann, J., Wulder, M.A., Wynne, R., 2008. Free access to Landsat imagery. Science 320, 1011. Zhu, Z., 2017. Change detection using landsat time series: A review of frequencies, preprocessing, algorithms, and applications. ISPRS J. Photogramm. Remote Sens. 130, 370–384. https://doi.org/10.1016/j.isprsjprs.2017.06.013. Zhu, Z., Woodcock, C.E., 2014. Continuous change detection and classification of land cover using all available Landsat data. Remote Sens. Environ. 144, 152–171. https:// doi.org/10.1016/j.rse.2014.01.011. Zhu, Z., Woodcock, C.E., 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 118, 83–94.

Remote Sens. 144, 94–104. https://doi.org/10.1016/J.ISPRSJPRS.2018.07.003. Cohen, W.B., Yang, Z., Healey, S.P., Kennedy, R.E., Gorelick, N., 2018. A LandTrendr multispectral ensemble for forest disturbance detection. Remote Sens. Environ. 205, 131–140. https://doi.org/10.1016/j.rse.2017.11.015. Crist, E.P., 1985. A TM Tasseled Cap equivalent transformation for reflectance factor data. Remote Sens. Environ. 17, 301–306. https://doi.org/10.1016/0034-4257(85) 90102-6. De Sy, V., Herold, M., Achard, F., Beuchle, R., Clevers, J.G.P.W., Lindquist, E., Verchot, L., 2015. Land use patterns and related carbon losses following deforestation in South America. Environ. Res. Lett. 10, 124004. https://doi.org/10.1088/1748-9326/10/ 12/124004. Duro, D.C., Franklin, S.E., Dubé, M.G., 2012. A comparison of pixel-based and objectbased image analysis with selected machine learning algorithms for the classification of agricultural landscapes using SPOT-5 HRG imagery. Remote Sens. Environ. 118, 259–272. https://doi.org/10.1016/J.RSE.2011.11.020. Ei, Kosaka, Y., Takeda, S., 2017. Underground biomass accumulation of two economically important non-timber forest products is influenced by ecological settings and swiddeners’ management in the Bago Mountains, Myanmar. For. Ecol. Manage. 404, 330–337. https://doi.org/10.1016/j.foreco.2017.09.015. Finer, M., Novoa, S., Weisse, J., Petersen, R., Souto, T., Stearns, F., Martinez, R.G., 2018. Combating deforestation: from satellite to intervention. Science 360, 1303–1305. https://doi.org/10.1126/science.aat1203. Foody, G.M., 2004. Thematic map comparison: evaluating the statistical significance of differences in classification accuracy. Photogramm. Eng. Remote Sens. 70, 627–633. https://doi.org/10.14358/PERS.70.5.627. Fox, J.M., Castella, J.C., Ziegler, A.D., Westley, S.B., 2014. Rubber plantations expand in mountainous Southeast Asia: what are the consequences for the environment? Asia Pac. Issues 114, 1–8. Hamunyela, E., Verbesselt, J., de Bruin, S., Herold, M., 2016a. Monitoring deforestation at sub-annual scales as extreme events in Landsat data cubes. Remote Sens. 8, 651. https://doi.org/10.3390/rs8080651. Hamunyela, E., Verbesselt, J., Herold, M., 2016b. Using spatial context to improve early detection of deforestation from Landsat time series. Remote Sens. Environ. 172, 126–138. https://doi.org/10.1016/j.rse.2015.11.006. Hermosilla, T., Wulder, M.A., White, J.C., Coops, N.C., Hobart, G.W., 2015a. Regional detection, characterization, and attribution of annual forest change from 1984 to 2012 using Landsat-derived time-series metrics. Remote Sens. Environ. 170, 121–132. https://doi.org/10.1016/j.rse.2015.09.004. Hermosilla, T., Wulder, M.A., White, J.C., Coops, N.C., Hobart, G.W., 2015b. An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites. Remote Sens. Environ. 158, 220–234. https://doi.org/10.1016/j.rse.2014.11.005. Hislop, S., Jones, S., Soto-Berelov, M., Skidmore, A., Haywood, A., Nguyen, T.H., 2018. Using landsat spectral indices in time-series to assess wildfire disturbance and recovery. Remote Sens. 10, 460. Huang, C., Goward, S.N., Masek, J.G., Thomas, N., Zhu, Z., Vogelmann, J.E., 2010. An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sens. Environ. 114, 183–198. https://doi.org/10. 1016/j.rse.2009.08.017. Huang, C., Wylie, B., Yang, L., Homer, C., Zylstra, G., 2002. Derivation of a tasselled cap transformation based on Landsat 7 at-satellite reflectance. Int. J. Remote Sens. 23, 1741–1748. https://doi.org/10.1080/01431160110106113. Hughes, M., Kaylor, S., Hayes, D., 2017. Patch-based forest change detection from Landsat time series. Forests 8, 166. https://doi.org/10.3390/f8050166. Huo, L.-Z., Boschetti, L., Sparks, A., 2019. Object-based classification of forest disturbance types in the conterminous United States. Remote Sens. 11, 477. https://doi.org/10. 3390/rs11050477. Kennedy, R.E., Yang, Z., Braaten, J., Copass, C., Antonova, N., Jordan, C., Nelson, P., 2015. Attribution of disturbance change agent from Landsat time-series in support of habitat monitoring in the Puget Sound region, USA. Remote Sens. Environ. 166, 271–285. https://doi.org/10.1016/j.rse.2015.05.005. Kennedy, R.E., Yang, Z., Cohen, W.B., 2010. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr – Temporal segmentation algorithms. Remote Sens. Environ. 114, 2897–2910. Kenney-Lazar, M., Wong, G., Baral, H., Russell, A.J.M., 2018. Greening rubber? Political ecologies of plantation sustainability in Laos and Myanmar. Geoforum 92, 96–105. https://doi.org/10.1016/J.GEOFORUM.2018.03.008. Khai, T.C., Mizoue, N., Kajisa, T., Ota, T., Yoshida, S., 2016. Stand structure, composition and illegal logging in selectively logged production forests of Myanmar: Comparison of two compartments subject to different cutting frequency. Glob. Ecol. Conserv. 7, 132–140. https://doi.org/10.1016/j.gecco.2016.06.001. Kuhn, M., 2008. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26. Li, P., Feng, Z., Xiao, C., Boudmyxay, K., Liu, Y., 2018. Detecting and mapping annual newly-burned plots (NBP) of swiddening using historical Landsat data in Montane Mainland Southeast Asia (MMSEA) during 1988–2016. J. Geogr. Sci. 28, 1307–1328. https://doi.org/10.1007/s11442-018-1527-4. Liang, L., Hawbaker, T.J., Zhu, Z., Li, X., Gong, P., 2016. Forest disturbance interactions and successional pathways in the Southern Rocky Mountains. For. Ecol. Manage. 375, 35–45. https://doi.org/10.1016/j.foreco.2016.05.010. Liaw, A., Wiener, M., 2002. Classification and regression by randomForest. R News 2, 18–22. Ling, L.C., Prescott, G.W., De Alban, J.D.T., Ziegler, A.D., Webb, E.L., 2017. Untangling the proximate causes and underlying drivers of deforestation and forest degradation in Myanmar. Conserv. Biol. 31, 1362–1372. https://doi.org/10.1111/cobi.12984. Masek, J.G., Vermote, E.F., Saleous, N.E., Wolfe, R., Hall, F.G., Huemmrich, K.F., Gao, F., Kutler, J., Lim, T.K., 2006. A Landsat surface reflectance dataset for North America, 1990–2000. IEEE Geosci. Remote Sens. Lett. 3, 68–72. https://doi.org/10.1109/ LGRS.2005.857030.

112