Lessons Learned While Implementing a Time-Series Approach to Forest Canopy Disturbance Detection in Nepal

: While deforestation has traditionally been the focus for forest canopy disturbance detection, forest degradation must not be overlooked. Both deforestation and forest degradation inﬂuence carbon loss and greenhouse gas emissions and thus must be included in activity data reporting estimates, such as for the Reduced Emissions from Deforestation and Degradation (REDD+) program. Here, we report on efforts to develop forest degradation mapping capacity in Nepal based on a pilot project in the country’s Terai region, an ecologically complex physiographic area. To strengthen Nepal’s estimates of deforestation and forest degradation, we applied the Continuous Degradation Detection (CODED) algorithm, which uses a time series of the Normalized Degradation Fraction Index (NDFI) to monitor forest canopy disturbances. CODED can detect low-grade degradation events and provides an easy-to-use graphical user interface in Google Earth Engine (GEE). Using an iterative process, we were able to create a model that provided acceptable accuracy and area estimates of forest degradation and deforestation in Terai that can be applied to the whole country. We found that between 2010 and 2020, the area affected by disturbance was substantially larger than the deforested area, over 105,650 hectares compared to 2753 hectares, respectively. Iterating across multiple parameters using the CODED algorithm in the Terai region has provided a wealth of insights not only for detecting forest degradation and deforestation in Nepal in support of activity data estimation but also for the process of using tools like CODED in applied settings. We found that model performance, measured using producer’s and user’s accuracy, varied dramatically based on the model parameters speciﬁed. We determined which parameters most altered the results through an iterative process; those parameters are described here in depth. Once CODED is combined with the description of each parameter and how it affects disturbance monitoring in a complex environment, this degradation-sensitive detection process has the potential to be highly attractive to other developing countries in the REDD+ program seeking to accurately monitor their forests.


Introduction
Forest loss due to deforestation and forest degradation is an important concern in combating climate change because it reduces the capacity of forests to act as a carbon sink, becoming instead a source of greenhouse gas emissions. While significant improvements in estimating emissions from deforestation have been made globally, tracking and monitoring forest degradation remains a challenge for many tropical forested countries. The Conference of the Parties (COP 13) to the United Nations Framework Convention on Climate Change (UNFCCC) acknowledged the impact of forest degradation in 2007, including it as a proposed mechanism for reducing emissions within the Reduced Emissions from Deforestation and Degradation (REDD+) program [1]. As many countries progress towards accessing and implementing result-based payments under the REDD+ mechanism, forest-related emissions from both deforestation and degradation must be assessed within a reasonable degree of uncertainty. To estimate these forest-related emissions, countries use activity data, which describes the extent of human activity causing emissions including deforestation and degradation and is often measured using estimates of the area of land cover change, coupled with carbon factors to estimate the loss of carbon pools [1,2].
Forest degradation has no universally accepted definition, which is one of the reasons it has been more difficult to quantify in a standardized procedure. The United Nations Food and Agriculture Organization (FAO) defines degradation as the "reduction of the capacity of a forest to provide goods and services" [3] while the Global Forest Observation Initiative (GFOI) has described it as a disturbance caused by human intervention in a forested landscape that results in carbon emissions but not a land cover change [2]. Satellite-based remote sensing observations are often the tool used to quantify forest loss of any sort, but the types of degradation that can be observed from space may not be congruent with policy definitions. Despite these challenges, carbon-focused policies require action by countries, many of which are drafting and adopting their own national definitions. For a forest assessment to be useful, the definition of degradation should be explicitly stated and linked with the tools used to measure it.
Although measuring degradation with remote sensing is challenging, it has become increasingly tenable. Detecting change with remote sensing can be as simple as comparing images before and after an event, but this coarse method is often insufficient for degradation processes observed in forests. Degradation, which includes fuel wood collection, smallscale timber harvesting, biomass loss from insects, and subcanopy fires, can be reliably detected only when the temporal density of the images is frequent enough to observe both long-duration and ephemeral changes in the forest structure [4]. Forests often regrow following a degradation event, and the emissions from the temporarily reduced forest would be missed if repeated observations were not made throughout the event's time period [5,6]. Additionally, degradation events are often spatially fine-scaled. Disturbances can be smaller than the resolution of commonly used Earth observations, detectable only through subpixel analysis [7][8][9]. Freely available satellite imagery may overcome some of these challenges. For example, Landsat 7 and 8 data have repeated observations every 16 days, 8 days when used in combination, at 30 m resolution, making degradation detection more feasible for countries struggling with the prohibitive cost of commercial remote sensing observations [6,9,10]. Both of these Landsat data sets are integrated into Google Earth Engine (GEE), a free cloud-computing platform that links vast image data sets with algorithms and large computational capacity [11]. As these advances make measurement of degradation more accessible in theory, it becomes important to test and document utility in practice.
Like many other countries, Nepal has established an emissions reduction payment agreement with the World Bank's Carbon Fund to manage forests sustainably. The Carbon Fund provides a mechanism for developing countries to receive monetary compensation for reducing forest-related greenhouse gas (GHG) emissions. Carbon Fund payments are designed to help countries and their stakeholders achieve long-term sustainability in financing forest conservation (https://www.forestcarbonpartnership.org/carbon-fund (accessed on 1 May 2021)). Countries in the program must establish a measurement, reporting, and verification (MRV) system with a national forest monitoring system (NFMS) to regularly report estimates of changes in forest emissions and carbon stocks [6,12,13].
If degradation accounts for more than 10% of these emissions, it must be included in the country's reporting [14,15]. Pearson et al. [14] performed an estimation of forest degradation emission for 74 developing countries between 2005 and 2010, finding that only 11 of the 74 countries had degradation that fell below this threshold. Nepal was reported to have around 90% of its emissions from degradation [14].
A key component of any forest monitoring program is identification of the processes that cause forest loss. Nepal's REDD+ Readiness Program Proposal identified the major deforestation and degradation drivers in the country [16][17][18]. Estimates from 2005 to 2010 indicated the degradation-specific drivers producing the most emissions in Nepal are, in order of decreasing magnitude, fuel wood collection, timber harvesting, and fire [14]. Because Nepal began promoting scientific forest management (SciFM) in 2012, the distribution of emission drivers may have shifted since then [19]. Alternatively, Nepal's own Strategy and Action Plan for 2018-2022 defined the leading drivers as unsustainable harvesting, overgrazing, fuelwood collection, forest fires, and conversion of forests to other land uses [20,21]. The relative impact of each driver varies over time according to sociopolitical and environmental conditions [17].
Through its documentation process, Nepal has made progress characterizing forest degradation. A key step was the development of the Forest Reference Level (FRL)/Forest Reference Emission Level (FREL) and Emission Reduction Project Document [18,21]. During the preparation of the FRL/FREL, only degradation due to fuelwood collection was assessed. Forest degradation due to grazing, timber extraction, and fire were not included because of a lack of reliable data. Nepal initially quantified its forest loss and degradation, or activity data, by analyzing the canopy change maps produced by Hansen et al. [22]. Later, the country worked to provide further detail to its assessment by using the Morphological Spatial Pattern Analysis Tool (MSPA) [23][24][25] for assessment of deforestation and forest degradation. This assessment categorized Nepal's forests as Inner Core or Edge Forests and mapped the changes between these categories from 2004 to 2014. The transition from Inner Core to Edge was considered a proxy for forest degradation. This strategy was supported by subsequent analysis of above-ground biomass estimates completed using LiDAR, which showed a distinct difference between the Inner Core and Edge Forest classes [21]. LiDAR data was available only in the Terai Arc Landscape portion of Nepal, which encompasses the northwest portion of Terai and extends into India. The data set was not sufficient for forest degradation analysis based on spatial point density (the data set has 0.8 points per square meter on average) and temporal density (only one time period is available). The FRL estimation process did not use LiDAR data because of expected difficulties with financing imaging again in the future [20]. Further, this approach uses the distance from the forest edge as a proxy for degradation activities, but in reality, these activities are not evenly distributed across forest edges across the landscape. Because of these limitations, an improved methodology was needed for estimating emissions due to forest degradation.
Here, by explicitly mapping locations of forest degradation activities, we report on efforts to continue developing forest degradation mapping capacity in Nepal. These maps will support Nepal's efforts to more accurately estimate activity data across the country within the context of REDD+. Our goal was to create a process to obtain more precise unbiased area estimates for activity data, particularly forest degradation. To achieve this goal, we created a map that represents forest degradation hotspots. This map is used to allocate samples using the map strata to target collection of data in a more efficient manner than a simple random sample. Meeting the desired precision goals using a stratified random sample requires interpretation of fewer points than a simple random sample. The maps cannot be used alone to estimate activity data because all maps have biased errors, for example due to class confusion [2]. We demonstrate the usefulness of employing freely available resources, which means this work could easily serve as a model for other developing countries working to establish reliable forest degradation monitoring methods using remote sensing. This project aims to assess degradation in Nepal using mapping and sample-based area estimation methods that are relevant to many other countries working with the Carbon Fund [2]. Countries including Nepal use these maps and sample interpretation to estimate activity data-the area of land undergoing land cover/use changes-with associated uncertainty intervals. We used the Continuous Degradation Detection (CODED) algorithm to detect changes in forests. This algorithm evaluates a time series of the Normalized Degradation Fraction Index (NDFI) to monitor forest canopy disturbances [26,27]. A key benefit of CODED is its use of linear spectral mixture analysis (SMA) to represent grades of degradation using the NDFI metric, by capturing subpixel and low-grade degradation events [27]. Another benefit is its application through a graphical user interface in Google Earth Engine (GEE), which capitalizes on cloud computing and integrated free Earth observation data sets including Landsat. CODED also differentiates between forest disturbances of deforestation and degradation, which is ideal for GHG emission estimations in REDD+ reporting.
CODED's many benefits make it a promising tool for estimating activity data; however, some challenges arise when implementing the tool in an applied setting. CODED requires specifying values of over a dozen parameters, which should be selected based on the forest dynamics and spectral conditions of the study region. Additionally, CODED was originally developed and applied in the Amazon Forest [9,26], so we needed to calibrate the parameter settings to a new geography. To address these challenges, we used in-depth parameter testing and results validation on an ecologically diverse region of Nepal. As a result of this process, we present an approach to parameter selection that can be applied to all of Nepal and to other study regions to inform sample allocation in degradation hotspots to estimate activity data. We also discuss the lessons we learned during this analysis and describe how individual parameter adjustments affect the analysis results. Our comprehensive description of the algorithm adjustment procedure will be a learning tool for future users and could save them a great deal of time and effort when tuning CODED to new locations. Consequently, this disturbance detection process has the potential to be highly attractive to other developing countries in the REDD+ program seeking to estimate activity data due to forest degradation and deforestation.

Study Region
Nepal is a diverse country-culturally, ecologically, and geographically. Within its borders there are 118 ecosystem types and 35 distinct forest types. Nepal is composed of five physiographic regions: the High Himal (part of the Himalayas including Mt. Everest; up to 8848 m elevation above mean sea level), the High Mountains (543-4951 m elevation), the Middle Mountains (110-3300 m elevation), the Churia (93-1955 m elevation), and Terai (the lowland plains; 63-330 m elevation) [28]. In this study, we focus on the Terai region ( Figure 1), which makes up 13.717% of Nepal's geographic area. The region, which is 20.88% forest and other wooded land area, contains a great amount of forest biodiversity, varied types of forest management practices, and a disproportionately high amount of forest disturbance activities. The four most prevalent forest types in Terai are Khair Sissoo, Sal, Tropical Mixed Hardwood, and Subtropical Mixed Hardwood. The climate of the Terai region is subtropical and characterized by monsoon seasons (June-September) and dry winters (December-February) [29]; in general, trees undergo leaf bud burst in February-April as temperatures rise and the rainy season begins, with greatest leaf cover around August and minimum leaf cover around January [30]. Average above-ground biomass is approximately 196 t/ha [29]. The Terai has been experiencing more urbanization and population densification due to migration from other physiographic regions, which has led to a higher demand for forest products [31]. Terai is the only region in Nepal where demand consistently outpaces supply for both timber and fuelwood. This pattern is projected to continue well into the future [17,21]. We studied these forest changes within the time period of 2010-2020. has led to a higher demand for forest products [31]. Terai is the only region in Nepal where demand consistently outpaces supply for both timber and fuelwood. This pattern is projected to continue well into the future [17,21]. We studied these forest changes within the time period of 2010-2020. The direct drivers of forest disturbance in the Terai physiographic region have variable impacts on the landscape that result in different patterns of forest disturbance. For example, urbanization and infrastructure development result in deforestation, permanently converting forest to another land cover type, while selective/illegal logging and livestock grazing result in forest degradation, the reduction in forest cover/quality without a full land cover conversion. Forest degradation events vary both in degree (partial to entire loss of forest cover) and temporal magnitude (short-term to long-term reduction in forest cover). Of particular concern in the Terai region is selective/illegal logging, which results in small-scale and widespread canopy loss. Recent studies suggest that forest degradation is an underreported source of carbon emissions [14,[32][33][34], making accurate area estimates of both deforestation and forest degradation critical to accurately estimating carbon impacts.

CODED
To detect both deforestation events and the more subtle small-scale forest degradation events caused by forest disturbance drivers of interest in the Terai ecoregion, we used the Continuous Degradation Detection (CODED) algorithm. CODED is a free and opensource tool that runs on GEE. By default, CODED uses freely available optical Earth observation imagery provided by Landsat satellites. The combination of a free and opensource algorithm, freely available imagery, and comprehensive documentation (https://coded.readthedocs.io/en/latest/background.html (accessed on 20 April 2021)) makes CODED an attractive tool for groups with limited budgets and capacity for developing custom-built solutions for detecting forest disturbance. CODED was recently updated to include new functionality offered in GEE. CODED produces maps with four strata classes: stable forest, stable nonforest, forest degradation, and deforestation.
All Landsat Collection 1 surface reflectance images available within the study period were cloud masked and used in the algorithm's analysis. CODED uses subpixel spectral mixture analysis (SMA) to calculate and evaluate time series changes in NDFI [32,35]. See Equations (4) and (5) in Souza et al. for how NDFI is calculated [27]. CODED uses linear spectral unmixing to represent subpixel fractions of key spectral endmembers, including green vegetation (GV), nonphotosynthetic vegetation (NPV), soil, shade, and clouds [9,27]. CODED then calculates NDFI, which is designed to analyze the endmember The direct drivers of forest disturbance in the Terai physiographic region have variable impacts on the landscape that result in different patterns of forest disturbance. For example, urbanization and infrastructure development result in deforestation, permanently converting forest to another land cover type, while selective/illegal logging and livestock grazing result in forest degradation, the reduction in forest cover/quality without a full land cover conversion. Forest degradation events vary both in degree (partial to entire loss of forest cover) and temporal magnitude (short-term to long-term reduction in forest cover). Of particular concern in the Terai region is selective/illegal logging, which results in small-scale and widespread canopy loss. Recent studies suggest that forest degradation is an underreported source of carbon emissions [14,[32][33][34], making accurate area estimates of both deforestation and forest degradation critical to accurately estimating carbon impacts.

CODED
To detect both deforestation events and the more subtle small-scale forest degradation events caused by forest disturbance drivers of interest in the Terai ecoregion, we used the Continuous Degradation Detection (CODED) algorithm. CODED is a free and open-source tool that runs on GEE. By default, CODED uses freely available optical Earth observation imagery provided by Landsat satellites. The combination of a free and open-source algorithm, freely available imagery, and comprehensive documentation (https://coded.readthedocs.io/ en/latest/background.html (accessed on 20 April 2021)) makes CODED an attractive tool for groups with limited budgets and capacity for developing custom-built solutions for detecting forest disturbance. CODED was recently updated to include new functionality offered in GEE. CODED produces maps with four strata classes: stable forest, stable nonforest, forest degradation, and deforestation.
All Landsat Collection 1 surface reflectance images available within the study period were cloud masked and used in the algorithm's analysis. CODED uses subpixel spectral mixture analysis (SMA) to calculate and evaluate time series changes in NDFI [32,35]. See Equations (4) and (5) in Souza et al. for how NDFI is calculated [27]. CODED uses linear spectral unmixing to represent subpixel fractions of key spectral endmembers, including green vegetation (GV), nonphotosynthetic vegetation (NPV), soil, shade, and clouds [9,27]. CODED then calculates NDFI, which is designed to analyze the endmember fractions within each pixel in a way that emphasizes the disparities between undisturbed and disturbed forest canopies. When there are greater proportions of soil and NPV, the NDFI value is low and indicative of degraded forest conditions or naturally sparse canopy cover [27,32]. As a result, CODED is able to detect even small-scale disturbances in the forest canopy, which is typical of forest degradation. CODED uses time-series analysis of NDFI to detect ephemeral forest degradation while reducing false positives due to seasonal vegetation patterns. A baseline NDFI pattern, which accounts for seasonal variation using harmonic regression, is established from a defined reference period (here pre-2010). To reduce processing requirements, CODED uses a forest mask, which may be generated using random forest models or provided as a separate layer.
Training data, comprising different land cover types in the study region, are required by the algorithm to distinguish between forest disturbances that cause a change in land use (deforestation) and those that do not (forest degradation). In this study, 3631 training data points were collected by the first author (RRA). These training data points were collected in Collect Earth (CE) desktop. Initially, a two-by-two kilometer point grid was laid over the country, and those points were classified into the six IPCC classes. For this research, we created a subset of points in the Terai physiographic region and converted these points to forest and nonforest points.
Changes in NDFI for each pixel were evaluated during the detection period (here 2010-2020) by comparing the observed level of NDFI with a prediction of the expected range of values from the baseline NDFI pattern. Forest disturbance events were then identified when a pixel fell outside of the expected range of NDFI a defined number of times. Multiple observations of change were required to lower the chance of false positives caused by noise.

Parameter Testing
Multiple parameters control how sensitive CODED is to detecting changes. These are particularly important for the model's sensitivity in detecting forest degradation, as forest degradation involves subtle changes. We iteratively tested several of these parameters to find the best performing combination for our study area in the Terai. The parameters we adjusted include two different versions of CODED ("Original" and "Updated") [36], different training data sets, the forest mask (including generating a forest mask within CODED), the tree cover threshold, the number of consecutive observations required to flag a change event, and the chi-squared probability, which controls the width of the statistical change window (Table 1). Parameters not shown in this table were left as the default value. The "Updated" version of CODED differs slightly from the original in the naming of the parameters and makes use of the GEE Continuous Change Detection and Classification (CCDC) implementation to efficiently run over larger areas [36]. Our original training data set contained points from multiple years, but the vast majority were from 2017. We therefore tried reducing the training data set to only these 2017 values to better agree with the requested single training year in the CODED interface. The forest mask can be imported or defined as a percent of tree canopy cover, the tree cover threshold parameter, from the Global Forest Watch (GFW) data set. The number of consecutive observations is the number of abnormal NDFI values required for the algorithm to flag a potential forest disturbance. Chi-squared controls the size of the moving time window used for detecting statistical changes. The change magnitude threshold is a post-processing parameter that determines what severity of NDFI change events will be counted as a disturbance in the final output.
We created 11 map iterations in CODED during our testing process (Table 1). Our first map (iteration 1), created in the original version of CODED, produced moderate results, so we decided to try using the updated version of CODED because we believed the parameters were easier to understand and adjust in that version. In contrast to the original version of CODED, the updated CODED does not create a forest mask on the fly, so we began with testing the cutoff threshold for defining the forest mask from the default GFW input [37]. Iteration 2 used a threshold of 80% canopy cover, the default for CODED, and resulted in high overestimation of forest disturbances. We reduced this default threshold with GFW to 50% in iteration 3 to examine the impact of this threshold on the results; this change reduced the overall forested area and the mapped degradations only slightly. The definition of forest used by Nepalese authorities requires only 10% tree canopy cover, but preliminary studies by Nepal's Department of Forest Research and Survey (DFRS) indicate that areas with at least 30% tree canopy cover best match typical Nepalese forest characteristics [21]. This discrepancy and the fact that canopy density varies across space and ecosystem type influenced our decision to test a wide range of tree canopy cover thresholds in our parameter testing.
At this point we hypothesized that using only one year (2017) of training data points might improve the results, as our original training data set contained some data points from 2016 and 2018 as well, and we learned that the updated CODED version assumes only one year of training data is being used. We did see visible improvement with this adjustment in map iteration 4, so decided to continue using this training data set from 2017 for future iterations.
In iteration 5, we used a GFW forest mask threshold of 20% to include a broader definition of forest, but this produced little change in the amount of forest degradation detected in the final results. We began testing alternative forest masks that we imported into CODED because we wanted to test whether the forest mask layer, which defines what CODED considers forest and controls the starting forest area, was the reason we were not detecting deforestation and degradation in areas where it was expected. Iteration 6 used the Global Land Analysis & Discovery group at University of Maryland's Total Canopy Cover layer with a 10% threshold; we chose this mask as an easily accessible free alternative to the GFW mask. Iteration 7 used the NLCMS forest class as a mask; we chose to test this mask because it was produced with the input of local experts on Nepal's forests instead of using global models. However, we saw little change with either option. We inferred that once a reasonable forest mask is found, it does not markedly improve the results to continue to search for a more accurate forest mask because these iterations had poorer accuracies than earlier attempts.
For iteration 8, we switched back to the original CODED version to see if our one-year training data improved the accuracy with that algorithm as well. Through qualitative inspection, we thought this was our most accurate result so far, and quantitative accuracy assessment confirmed this (overall accuracy 92.2%). The user's and producer's accuracies were also higher, especially with deforestation-type disturbances ( Figure 2). This is the map we stratified in order to produce our verification data set.   Table 1. Figure 2. Comparison of overall map accuracy and user's and producer's accuracy for the map iterations of CODED described in Table 1.
After further consideration, we decided to also test the parameter for the number of consecutive observations required to label a disturbance event. We increased the required consecutive observations from 4 to 6 (iteration 9) and to 10 (iteration 10). However, these both decreased the number of detected disturbances, even in areas where we knew from on-the-ground observations that degradation had occurred.
At this point, we had multiple variables that were known to influence the qualitative accuracy of the map, and others that we suspected may be important, specifically the chi-squared parameter. Based on our accuracy assessments (Figure 2), we noted that the forest mask choice (iterations 5 and 7) and the number of consecutive observations (iterations 9 and 10) decreased accuracy notably by changing the amount of deforestation and degradation detected. We decided to conduct methodical tests of these parameters in a small region of the Terai that we had identified as having frequent classification errors. We also examined variations of the chi-squared and change magnitude threshold parameters in this manner. Based on the outcome of these methodical tests, we created map iteration 11 using the best combination of parameter values and the forest mask we generated in map iteration 8 (Table 1).

Accuracy Assessments and Unbiased Area
To assess map and CODED model accuracy, we used both qualitative and quantitative measures. Qualitatively, we visually compared the areas mapped as degradation for each map iteration. We had access to local information of known regions that had experienced degradation within our study period and were able to use these areas for qualitative checks and comparison of each map's performance. In addition, we noted locations where there was disagreement between the interpreter and the CODED map.
For each map, we calculated quantitative measures of users' and producers' accuracies using confusion matrices. Validation sample locations were chosen based on the stratified map output of iteration 8. Using GEE's built-in stratified random sample algorithm, we created a random stratified sample of the map with 250 sample points in each of the four map classes (stable forest, stable nonforest, forest degradation and deforestation) for a total of 1000 points [38,39]. We used Collect Earth Online (CEO), a free, open-source, web-based tool that facilitates data collection, to collect validation data based on this random stratified sample [40]. We leveraged CEO's degradation tool, which calculates a time-series graph of the normalized difference fraction index (NDFI) for each validation point based on a time series of Landsat imagery, allowing the interpreter to detect forest degradation [27]. For consistency, the author most familiar with the Terai landscape (RRA) labeled all of the validation data plots and created an interpretation key with examples of deforestation and forest degradation. Bias introduced by interpretation errors is a well-known issue for studies like this one using interpretation of remotely sensed imagery for validation. We recognize this is a reality in most practical applications of forest disturbance monitoring and assessment, and we designed our method of addressing bias in a way that could be easily implemented by most organizations. This bias cannot be fully eliminated, but it can be mitigated by increasing the number of trained interpreters used or through continued monitoring of interpretation accuracy and consistent feedback and discussion to improve interpretations [41,42]. Often only one skilled interpreter with local forestry knowledge is available; this was true for our team, so we utilized the latter approach. Quality assurance was conducted by reviewing a subset of the validation data points twice and discussing interpretation labels of samples with imagery that was challenging to interpret with the broader team. The interpretation key used for guiding the process to label samples was improved through collaborative review and discussion, providing a standardized interpretation process (Interpretation Supplementary Materials Key S1).
In line with the CODED algorithm, for the validation effort, forest disturbance was defined as a temporary reduction and subsequent recovery of NDFI and vegetation cover without a land cover change from forest, while deforestation was defined as a reduction in NDFI and conversion to another land use, such as agriculture or water (e.g., a riverbed shift). Of note, these differed from the Nepalese working definitions, which require a minimum contiguous area of canopy cover to define an area as forest and then use changes in canopy cover to determine whether forest degradation or deforestation have occurred [21,27]. As CODED is a pixel-based algorithm, there is no minimum size requirement for forest land cover.
The resulting validation data were used to calculate an error matrix and accuracy metrics for each map, using an area weighting approach when sampling strata differed from those of a given map iteration. We used the error matrix in an unbiased ratio estimator described in [38] to calculate unbiased area estimates and confidence intervals for each class [41]. These unbiased area estimates and associated confidence intervals are the key measures of activity data in the REDD+ context and are critical for meeting contractual agreements [2].

Parameter Testing
We created 11 iterations in CODED during our testing process, which were qualitatively and quantitatively evaluated and compared with prior iterations (Table 1; Figure 2). Through our iterative process we were able to considerably increase the producer's accuracies for both deforestation and degradation ( Figure 2). The most influential variables we adjusted included the forest mask used, chi-squared, and the post-processing magnitude threshold.

Final Map
Using our final map, we calculated unbiased area estimates indicating 105,650 hectares (+/− 37,264 hectares) of forest degradation and 2753 hectares (+/− 3525 hectares) of deforestation occurred during the period of 2010-2020 in the Terai region of Nepal (Figure 3). This translates to 4.63% and 0.12% of the total area of the Terai and 18.73% and 0.69% of the starting area of forest, respectively.
3). This translates to 4.63% and 0.12% of the total area of the Terai and 18.73% and 0.69% of the starting area of forest, respectively.
Overall map accuracy was 91.7%. Producer's accuracy was 16% for degradation events and 84% for the deforestation events. User's accuracy was 34% for degradation and 22% for deforestation ( Table 2). While some of the disagreement between the interpreter and CODED map results occurred sporadically throughout the landscape, there was consistent disagreement in areas of sparse forest and in contiguous forests near rivers.   Overall map accuracy was 91.7%. Producer's accuracy was 16% for degradation events and 84% for the deforestation events. User's accuracy was 34% for degradation and 22% for deforestation (Table 2). While some of the disagreement between the interpreter and CODED map results occurred sporadically throughout the landscape, there was consistent disagreement in areas of sparse forest and in contiguous forests near rivers. Table 2. Error matrix with sample counts for the Map 11 iteration. Each grid cell is the correspondence between the final map (rows) and the reference data (columns). Note that the row (map values) totals do not correspond to the 250 samples originally derived from each class because the map is different from the stratification used to create the sample (see Section 2.4). Stable forest  208  25  92  26  351  Stable nonforest  1  249  37  1  288  Degradation  28  13  116  21  178  Deforestation  0  35  88  60  183  Total  237  322 333 108

Discussion
Iterating across multiple parameters using the CODED algorithm in the Terai region of Nepal and generating maps and unbiased area estimates of activity data has provided a wealth of insights not only for detecting forest degradation and deforestation in Nepal but also for the process of using tools like CODED in applied settings. Using this process, we were able to create unbiased area estimates for forest degradation activity data (105,650 +/− 37,264 hectares) with a margin of error (precision) of 35%, improving on previous approaches and providing pathways to both fulfill contractual agreements within the REDD+ context [2] and allow stratified samples derived from the outputs for subsequent fieldwork to be small enough to be financially and logistically feasible for Nepal. Because the Terai region is the most diverse and challenging area to classify in Nepal, we believe that this process will achieve similar or better results when applied to the entire country. Our results also highlight the importance of degradation as a source of carbon emissions [14,[32][33][34]. This work can serve as a model for other developing countries working to establish reliable forest degradation area estimation methods using remote sensing.
We found that CODED model performance, measured quantitatively and qualitatively using producer's and user's accuracy, varied dramatically based on the model parameters specified, with certain parameters particularly influencing accuracy (Figure 2). For example, the effect of the magnitude threshold parameter was clearly seen spatially as overestimation or underestimation when known areas of forest degradation were examined. The reduction of chi-squared led to markedly higher producer's accuracy for deforestation and degradation. The choice of the tree cover threshold seemed to have a large impact on the user's accuracy for deforestation detection.
While the majority of the work done for this project was focused on generating the best-quality map of activity data-specifically forest degradation hot spots-ultimately, our approach to obtain unbiased area estimates relies on the interpretation of images across a probabilistic sample allocated using a stratified random sample. The activity data maps have biased errors due to class confusion or edge effects, so we use these maps with a stratified sample for area estimation [38,39,43]. The samples are then labeled as being impacted by forest degradation activities or not. While ground observations are often considered more accurate sources for assessing activities at these sample locations, these can be cost-prohibitive, and it can be difficult to assess presence of and timing of past degrad ation activities in the field, depending on how fast the forest has recovered. An alternative approach is the interpretation of a time series of available high and moderate resolution imagery, a chronosequence approach. Remote sensing teams led by McRoberts [41] and Pengra [42] demonstrated that this approach too introduces bias because interpreters can mislabel conditions, especially for more-difficult-to-interpret land cover classes, and depending on the quality of the imagery. Their reccomendation to mitigate for these errors and bias is to increase the quality control procedures involved with photo-interpretation such as by having a greater number of interpreters label each sample and training protocols to enhance interpreter agreement. Recommended quality assurance procedures will be used during the operational process for estimating and reporting emissions estimates [41,42].
We observed that most disagreements between the CODED maps and the interpreter occurred in areas of sparse forest and in contiguous forests near rivers. Because of this, we believe our lower producer's accuracy for degradation was likely due to classifications in our reference sample rather than the CODED parameter choices. We will test this hypothesis in future planned work through a rigorous examination of samples, especially in the difficult-to-evaluate sparse forest regions. Additionally, the fact that CODED's spectral end-members were trained in a tropical environment may also have depressed our accuracies for degradation events. These same end-members have been used successfully many times in environments outside the Amazon [27,32,44,45]. While development of new end-members is outside the scope of this work, we hope to examine this in future work. However, because degradation is harder to detect remotely than deforestation, and given that our achieved precision met our reporting goals, we believe that our process is successful despite this.

Lessons Learned from Using CODED in an Applied REDD+ Context
While exploring the different parameters in CODED, we found two key things directly applicable to using tools like CODED in an applied REDD+ context where mapping and estimating the area of degradation are the objective. First, model accuracy was highly sensitive to parameter choice (Figure 2), and second, while pre-built tools drastically reduce the time needed to create classifications, significant effort and domain knowledge are still needed to create a good classification result.
The CODED algorithm uses chi-squared to control the size of the window for detecting statistical change. A change magnitude threshold variable that filters out low-magnitude change detections that may be erroneous or insignificant is applied in post-processing. Changing either of these parameters notably increased or decreased the amount of change that was detected (Table 3, Figure 4). The data pre-processing parameters that we tested, including creating temporal composites of the input data, did not have as strong an impact. To better prepare future users for what to expect when iterating through some of the most influential variables within CODED, we completed model runs changing only one parameter at a time and displayed the resulting map effect on a region from our study area known to have degradation from timber harvesting and deforestation (Figure 4).   Further, as we were reviewing model iterations in CEO, we found patterns of omission error that provide important insight into NDFI, SMA, and CODED. One pattern of omission error is the result of NDFI saturating at 1 and not capturing variation in Terai ( Figure 5; note comparison with SWIR2, which was not used for analysis but is shown here to demonstrate a nonsaturating signal). Specifically, under certain ecological conditions, NDFI reaches a maximum value prematurely such that it does not capture all of the variation. The ecological conditions noted here are high leaf density during the height of the monsoon season around August [30]; further exploration of the saturation phenomenon may reveal additional conditions. As the CODED model is based on SMA, this observation suggests that the spectral end member approach may need to be recalibrated for our study area, and potentially other study areas as well. The spectral signature definitions for the end members used in CODED and NDFI were generated for the Amazon [26,27]. Although both CODED [32] and NDFI have been successfully used in other environments (e.g., West Africa; [44,45]), the assumption that the same end members work for Nepal and areas outside of the Amazon must be further examined. Other researchers working in Southeast Asia have also made note of issues with the SMA approach. For example, in Kalimantan researchers found that spectral end-member selection needed to be adjusted to account for differences in terrain and atmospheric conditions causing variation in spectral reflectance [46][47][48].
A second pattern of forest degradation omission arose from the approach to forest masking used by CODED. In the Terai region, there were some areas that were classified as various nonforest land covers in 2010 when the forest mask was applied. The verification data collected identified these areas as forest degradation or deforestation based on second-growth forest that had grown in the intermediate time period. For CODED to detect these areas as degradation, however, they must be identified in the forest mask as forest because CODED does not process the nonforest areas under the often-correct assumption that nonforested areas like agriculture, water, and urban areas are not likely to become forested [26]. However, this assumption may not hold in Southeast Asia, where cyclical forest-agriculture patterns are common (e.g., shifting agriculture), or in Nepal, where rivers may change course. This finding highlights that use of a forest mask limits the detection of disturbances in secondary forest, and that careful consideration of the appropriate time period of the forest mask is critical if disturbance to secondary forest is a concern.
Other sources of disagreement between the map results and the interpreter in CEO may be due to complicated environmental scenarios. For example, one issue was mixed  Further, as we were reviewing model iterations in CEO, we found patterns of omission error that provide important insight into NDFI, SMA, and CODED. One pattern of omission error is the result of NDFI saturating at 1 and not capturing variation in Terai ( Figure 5; note comparison with SWIR2, which was not used for analysis but is shown here to demonstrate a nonsaturating signal). Specifically, under certain ecological conditions, NDFI reaches a maximum value prematurely such that it does not capture all of the variation. The ecological conditions noted here are high leaf density during the height of the monsoon season around August [30]; further exploration of the saturation phenomenon may reveal additional conditions. As the CODED model is based on SMA, this observation suggests that the spectral end member approach may need to be recalibrated for our study area, and potentially other study areas as well. The spectral signature definitions for the end members used in CODED and NDFI were generated for the Amazon [26,27]. Although both CODED [32] and NDFI have been successfully used in other environments (e.g., West Africa; [44,45]), the assumption that the same end members work for Nepal and areas outside of the Amazon must be further examined. Other researchers working in Southeast Asia have also made note of issues with the SMA approach. For example, in Kalimantan researchers found that spectral end-member selection needed to be adjusted to account for differences in terrain and atmospheric conditions causing variation in spectral reflectance [46][47][48].
documented. Another challenge is that we ran 11 different iterations to arrive at a satisfactory model. Each iteration involved model specification, model execution, and model evaluation phases. Model evaluation required both local knowledge, to be able to qualitatively identify areas of success and failure in the model, and statistical knowledge, in order to avoid re-stratification and revalidation for each model iteration [38].
However, for many use cases of CODED in applied settings, including this project, the main goal is not map accuracy but instead creating unbiased estimates of the area of forest degradation and deforestation with associated uncertainty. These estimates are critical in the REDD+ context, where estimates of activity data to required levels of precision are contractually obligated. Thus, spending substantial effort and time on the improvement of map accuracy is not as worthwhile as improving precision in area estimates. This precision is achieved through generation of more reference samples. For other applied teams working to achieve satisfactory classification results using CODED for area estimation purposes, some of our lessons learned may be insightful. We learned the importance of creating definitions for forest degradation early in a project's lifecycle. At the beginning, the author collecting validation data (RRA) used the Nepalese definition of degradation (partial forest loss while still above 10% tree canopy cover; [21]), but our map was generated assuming degradation meant any forest disturbance that does not result in a land cover change from forest. This definition is in accordance with how CODED defines degradation as multiple anomalously low NDFI readings (specifically five consecutive readings for our final map iteration) not followed by a detected land cover change from forest. Once this miscommunication was realized and resolved, our estimates of overall map accuracy improved by 13%. A second pattern of forest degradation omission arose from the approach to forest masking used by CODED. In the Terai region, there were some areas that were classified as various nonforest land covers in 2010 when the forest mask was applied. The verification data collected identified these areas as forest degradation or deforestation based on secondgrowth forest that had grown in the intermediate time period. For CODED to detect these areas as degradation, however, they must be identified in the forest mask as forest because CODED does not process the nonforest areas under the often-correct assumption that nonforested areas like agriculture, water, and urban areas are not likely to become forested [26]. However, this assumption may not hold in Southeast Asia, where cyclical forest-agriculture patterns are common (e.g., shifting agriculture), or in Nepal, where rivers may change course. This finding highlights that use of a forest mask limits the detection of disturbances in secondary forest, and that careful consideration of the appropriate time period of the forest mask is critical if disturbance to secondary forest is a concern.
Other sources of disagreement between the map results and the interpreter in CEO may be due to complicated environmental scenarios. For example, one issue was mixed pixels, where a pixel was included in the forest mask but could be seen through highresolution imagery to be partially composed of another land cover class. Further discussion is needed regarding how interpreters should handle labeling these types of pixels. How CODED handles mixed pixels is also still somewhat unclear. Regions along land use class borders, which vary slightly over time and with imagery horizontal accuracy, are likely more prone to labeling confusion by both the algorithm and human-led validation efforts. It would be advantageous to have repeated validation efforts by different interpreters on all of these points of confusion, but this is not always possible when resources are limited.
Second, achieving good classification results using CODED required significant effort and domain knowledge. Understanding the challenges faced in applied settings will help to build better tools. One challenge was that the connection between parameter and model output was not always clear. Similarly, the parameter names changed between the previous iteration of CODED and the current version, and not all of these changes were documented. Another challenge is that we ran 11 different iterations to arrive at a satisfactory model. Each iteration involved model specification, model execution, and model evaluation phases. Model evaluation required both local knowledge, to be able to qualitatively identify areas of success and failure in the model, and statistical knowledge, in order to avoid re-stratification and revalidation for each model iteration [38].
However, for many use cases of CODED in applied settings, including this project, the main goal is not map accuracy but instead creating unbiased estimates of the area of forest degradation and deforestation with associated uncertainty. These estimates are critical in the REDD+ context, where estimates of activity data to required levels of precision are contractually obligated. Thus, spending substantial effort and time on the improvement of map accuracy is not as worthwhile as improving precision in area estimates. This precision is achieved through generation of more reference samples.
For other applied teams working to achieve satisfactory classification results using CODED for area estimation purposes, some of our lessons learned may be insightful. We learned the importance of creating definitions for forest degradation early in a project's lifecycle. At the beginning, the author collecting validation data (RRA) used the Nepalese definition of degradation (partial forest loss while still above 10% tree canopy cover; [21]), but our map was generated assuming degradation meant any forest disturbance that does not result in a land cover change from forest. This definition is in accordance with how CODED defines degradation as multiple anomalously low NDFI readings (specifically five consecutive readings for our final map iteration) not followed by a detected land cover change from forest. Once this miscommunication was realized and resolved, our estimates of overall map accuracy improved by 13%.
Additionally, we learned that model parameter selection is best done by adjusting one parameter at a time in a methodical way. A full sensitivity analysis may produce marginally higher map accuracies, but an analysis of this complexity is impractical for most real-world applications. We recommend thoroughly testing one parameter at a time so the users can clearly observe the impacts of adjusting that parameter in their particular environment. Our general explanation and examples of how each CODED parameter impacts the map (Table 3 and Figure 4) should provide support for this process, but each environment will vary slightly in the magnitude and distribution of these impacts, so careful parameter testing for each use case is recommended. A more complete analysis would involve retroactively altering parameters after making alterations to others, analyzing how the parameter changes interact, and seeing if a change to one alters the previously ideal setting for another. However, a full sensitivity analysis or this comparatively complex level of testing would not be feasible for many users. In this study, we seek to explain a realistic approach for practitioners, detailing an achievable methodology that results in an acceptable level of accuracy for applied uses.
The graphical user interface makes this process of repeated parameter adjustment simple for users at any skill level. It is also important for users to understand the relative magnitude of impact from each parameter and how it interacts with changes in the landscape. For example, we found the chi-squared probability parameter, which determines the size of the statistical boundary surrounding the model residuals that is used to identify change, to lead to the most obvious map changes. Adding map examples of relative magnitude and spatial impacts on forest degradation and deforestation detection from each parameter (e.g., Figure 4) would be a useful addition to the documentation of the CODED algorithm as it becomes more widely used. Further, our parameter selection occurred in a subtropical forest; users in other ecosystems, such as those with periodic or permanent snow cover, will have different concerns and may have different key model parameters. For example, in regions with significant snow cover, NDFI would need to be adjusted to account for the fact that green vegetation fractions are probably not visible year-round. Ecosystems with lower tree cover will have correspondingly lower NDFI due to exposure of the understory and understory vegetation. These regions will not have the same issues with NDFI saturation as observed in the Terai. Additionally, while spectral mixture analysis is specifically designed to account for the problem of mixed pixels, in regions with less canopy cover a mixed forest pixel will look similar to a degraded forest, which is why the time series component is so important. For the Terai and other ecosystems, a critical first step is to model the "typical" state of the forest and look for changes from that baseline.
Moving forward, the usability, documentation, and knowledge transfer to applied teams examining local deforestation and forest degradation patterns will be critical. Future creation and modifications to tools like CODED meant for use within the broader community should include not only comparisons of accuracy to existing products but also comparisons of usability, ease of use for applied users, limits to the generalizability of the model or metrics used, and advice for effectively tuning the model. Future research in Nepal should include additional examination of whether different drivers of forest degradation are differently detectible and comparison of approaches for detecting degradation and deforestation, including using different SMA and modelling approaches. More broadly, future research in the area of forest degradation detection must address how SMA approaches NDFI saturation in certain environments like the Terai in Nepal.

Conclusions
This project was a first attempt at applying a process for tuning CODED, a method for detecting forest disturbance, in the Terai region of Nepal to create unbiased area estimates of activity data with sufficient precision to support contractual agreements [2] and allow stratified samples derived from the outputs for subsequent fieldwork to be small enough to be financially and logistically feasible. Success in this ecologically and geographically diverse study area with multiple forest deforestation and degradation drivers suggests that the process will be successful when applied to all of Nepal. We focused much of our attention on testing parameters within the CODED algorithm-first to improve the accuracy and precision of this analysis, but also to contribute to the literature a more thorough understanding of the relative implications of each parameter adjustment.
We found the chi-squared parameter, the number of consecutive observations required to define a disturbance event, and the NDFI change magnitude threshold to be especially influential to the results, with the quality of the forest mask and the training data also important aspects to consider. It took us 11 iterations to establish an acceptable set of parameters for our analysis. Lessons learned through the lengthy process highlighted the positive qualities and areas for potential improvement within the CODED algorithm. CODED's use of NDFI and the NDFI prediction algorithm makes it useful for detecting forest degradation, which is often small-scale and ephemeral but nonetheless a large portion of GHG emissions [14,15]. The graphical user interface on a free and open-source platform with cloud-computing also gives it wide appeal. However, CODED requires a high level of domain knowledge to achieve accurate classifications, as well as local knowledge of the forest conditions for useful validation. We have provided an in-depth examination and documentation of the parameters with examples to contribute to the future usefulness of CODED for other countries looking to accurately measure deforestation and forest degradation activity data for implementing result-based payments and for forest management. Others looking to implement similar frameworks will need to work through these same steps to tune the CODED algorithm to their specific ecological context; the most important parameters in other forest ecologies will likely differ from those we found in the Terai.
Nepal is in the process of implementing the Emission Reduction Program and following the required MRV process. In the past, Nepal used proxy-based approaches to estimate forest degradation, which allocated a 15% conservativeness factor as an uncertainty buffer to total estimated emission reductions associated with forest degradation [49]. Uncertainty buffers are utilized by the Emission Reduction Program to encourage countries to reduce the uncertainty in their estimates and to mitigate the risk of overestimation [50]. Our efforts using CODED are an important step in order to reduce the uncertainty buffer so that Nepal can realize more of the monetary benefits from their emission reduction efforts. This study provides the activity data for forest degradation due to selective logging. This type of local data source reduces the uncertainty buffer required for emission reduction payment mechanisms. This will reduce the buffer estimating for carbon payment. At a later date, the process described here will be applied to the entire country and additional drivers of forest degradation and deforestation in Nepal in this ongoing collaboration.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13142666/s1, The forest degradation interpretation key used for this analysis is provided as S1: Nepal Forest Degradation Interpretation Key.  Data Availability Statement: For this analysis we used publicly available tools: CODED, available at https://github.com/bullocke/coded (accessed on 28 October 2020). CODED draws on publicly available LANDSAT 5, 7, and 8 data. The data presented in this study are available on request from the corresponding author. The data are not publicly available as the work is still ongoing.