Monitoring long-term forest dynamics with scarce data: a multi-date classification implementation in the Ecuadorian Amazon

ABSTRACT Monitoring long-term forest dynamics is essential for assessing human-induced land-cover changes, and related studies are often based on the multi-decadal Landsat archive. However, in areas such as the Tropical Andes, scarce data and the resulting poor signal-to-noise ratio in time series data render the implementation of automated time-series analysis algorithms difficult. The aim of this research was to investigate a novel approach that combines image compositing, multi-sensor data fusion, and postclassification change detection that is applicable in data-scarce regions of the Tropical Andes, exemplified for a case study in Ecuador. We derived biennial deforestation and reforestation patterns for the period from 1992 to 2014, achieving accuracies of 82 ± 3% for deforestation and 71 ± 3% for reforestation mapping. Our research demonstrated that an adapted methodology allowed us to derive the forest dynamics from the Landsat time series, despite the abundant regional data gaps in the archive, namely across the Tropical Andes. This study, therefore, presented a novel methodology in support of monitoring long-term forest dynamics in areas with limited historical data availability.


Introduction
The Amazon rainforest constitutes one of the biologically most diverse, structurally complex, and carbon-rich bioregions of the world (Asner et al., 2014). It performs essential global-scale functions and provides a multitude of ecosystem services (Paula et al., 2014). Tropical deforestation is a major threat to the region and a driver of climate change with potentially critical impacts on the biosphere (Fearnside, 2005). Large-area deforestation assessments indicate that the Amazon Basin lost 13.3% of its forest from 2000 to 2013, where the headwater basins suffered most of the pressure (RAISG, 2015). This is also particularly alarming for the region itself, as the highland Amazon (or Tropical Andes) is highly susceptible to global warming (Karmalkar, Bradley, & Diaz, 2008), while being under-researched in deforestation studies (Armenteras, Rodríguez, Retana, & Morales, 2011). Moreover, only little is known about forest succession (Barbosa, Broadbent, & Bitencourt, 2014) or land-cover intensities (Kuemmerle et al., 2013) in this subregion, which are also important components for understanding the impacts on ecological services (e.g. on biodiversity, carbon sequestration, or nutrient sinks; Brown & Lugo, 1990;Edwards, Massam, Haugaasen, & Gilroy, 2017;Poorter et al., 2016). Monitoring forest dynamics in the Tropical Andes, therefore, plays a key role for informing policymakers and resource managers in their decisionmaking processes over the next few years (Angelsen & Wertz-Kanounnikoff, 2008;De Koning et al., 2011). While lowland tropical forests have been well researched, closing the remaining knowledge gaps on forest dynamics in Andean tropical forests is of prime importance (Armenteras, María, Rodríguez, & Retana, 2017;Da Ponte et al., 2015;Oliveira, Eller, Bittencourt, & Mulligan, 2014;Spracklen & Righelato, 2014).
Comprehensive forest dynamics monitoring has traditionally implied decades of field observations (Fragal, Silva, & Novo, 2016). Such observations are costly, and it might even be impossible to collect the necessary data in the field. Remote sensing offers a unique alternative for supplying this information for large areas and in different spectral, spatial, and temporal resolutions. Terra and Aqua on board of the Moderate Resolution Imager Spectroradiometer (MODIS) have regularly been used for broad-scale mapping of forest dynamics (Hansen et al., 2008). However, the spatial resolution of MODIS data is limited when fine-scale disturbance regimes prevail such as from selective logging, skid trails, or disturbances related to landslides or local windthrow. The Landsat sensor family with its 30-m spatial resolution and its 45-year observation record is better suited for capturing long-term and fine-scale processes. It is the most widely used observation system for land-cover and land-use change (LCLUC) assessments and forest dynamics monitoring programs (Camara, 2013;Hansen & Loveland, 2012). Since the launch of Landsat-1 in 1972, the Landsat Program has continuously collected data across the globe andsince the launch of the Landsat Thematic Mapper (TM) in 1982in six spectral bands covering the optical, near-infrared, and shortwave infrared wavelength regions. For these reasons, it is the most long-term medium-resolution Earth observation satellite archive available. Due to the open data policy since 2008, Landsat data are available free of charge as a standard high-level product for long-term LCLUC analysis . This development allowed major improvements in automated time-series analysis, leading to a range of novel algorithms such as TIMESAT, LandTrendr, BFAST, and CCDC (Eklundh & Jönsson, 2015;Kennedy, Yang, & Cohen, 2010;Verbesselt, Hyndman, Newnham, & Culvenor, 2010;Zhu & Woodcock, 2014) that allow for extracting forest dynamics information. However, in some regions around the globe, the archive data density is considerably lower, mostly due to persistent cloudiness (Arvidson, Gasch, & Goward, 2001;Chance, Hermosilla, Coops, Wulder, & White, 2016) reduce data quantity and quality. This results in time series with poor signal-to-noise ratio, as useful information about forest status is weak and not significant to differentiate from random noise. This is a limitation for transferring these novel algorithms to datascarce regions, such as the Tropical Andes (Santos, Dubovyk, & Menz, 2017), as time series analyses require rather dense data stacks over time.
Conceptual approaches including those based on image compositing, multi-sensor fusion, and postclassification change detection have demonstrated their potential to overcome these limitations (Griffiths et al., 2014;Hansen et al., 2013;Potapov et al., 2012). Multi-date classification (Zhu, 2017) has the advantage to cope with noise-prone observations, especially when based on cloud-free composites that not only allow identifying deforestation or reforestation processes but also to characterize land-use intensities from postdeforestation dynamics (Rufin, Müller, Pflugmacher, & Hostert, 2015). For these reasons, multi-date classification is the methodology of choice for monitoring longterm forest dynamics in areas such as the Tropical Andes. However, implementing a multi-date classification scheme under specific regional conditions can still be challenging, for example, due to poor data availability in the past depending on historical data receiving strategies or increased cloudiness due to topography. Consequently, we applied this methodology to a study case located in the Amazon region of Ecuador, the Upper Napo Watershed (UNW), where heavy rainfall regimes (Espinoza et al., 2015) and complex landscapes (Asner et al., 2014) are major impediments. Our overarching objective was to monitor long-term forest dynamics and identify deforestation/ reforestation for the period between 1992 and 2014. We pose the following research questions related to these objectives: • Which processing steps and techniques are needed to implement multi-date classification in a data sparse region, as exemplified in the UNW? • How well does a multi-date classification approach perform when monitoring long-term forest dynamics in the environments of the Tropical Andes?

Study area
The UNW is located between 78°25′W and 76°25′W longitude and 0°10′N and 1°30′S latitude (Figure 1 (a)). It covers an area of about 12,500 km 2 in the Ecuadorian Amazon, spreading across the three provinces Napo (63% of the watershed), Orellana (26%), and Pastaza (9%). The altitudinal gradient of the Andes covers~260-5600 m above sea level (a.s.l.). Mean temperatures varying from −0°C to 26°C and annual precipitation from 1100 to 5300 mm (MAE, 2013). The core rainy season extends from December to May, but fog and clouds are abundant throughout the year, especially at higher elevations (Ramírez, Teuling, Ganzeveld, Hegger, & Leemans, 2017). The complex geology creates diverse edaphic conditions that in combination with topographic gradients and climatological impacts results in a multitude of extremely species-rich ecosystems (Hoorn et al., 2010). The percentage of natural vegetation in the UNW was estimated to be 79% in 2014, with 76% of forest, 2% of shrub-dominated landscapes, and 1 % of grasslands (MAE, 2017). Forests are generally evergreen, but forest ecosystems vary substantially in tree composition, flood regimes, and topographic and bioclimatic boundary conditions (MAE, 2013). According to Guariguata and Ostertag (2001), reforestation occurs as quickly as in 5 years after a disturbance in the evergreen forest ecosystems of the Ecuadorian Amazon, with variations depending on past land use practices and abiotic site conditions. This was verified during fieldwork in May 2017 at three reforestation sites in the UNW, where 46 hemispherical photographs were acquired. We followed Pueschel, Buddenbaum, and Hill (2012) to binarize these photographs and derive the canopy closure index. We found that canopy closure after 5 years can be greater than of a 20-year-old forest (Figure 1(b)). We accordingly used a time threshold of 5 years as a reference for mapping reforestation (compare "Materials and methods" section). Four National Protected Areas, mainly created during the 1990s, are present in the UNW and cover 25% of the total area. 78,300 ha of native forest were reported to be converted to pastures and croplands between 2000 and 2014, resulting in an average annual deforestation rate of~6300 ha or 0.5% of the forested land (MAE, 2017). While these numbers vary, deforestation rates do not exceed 1.2-1.6% yr −1 (Sierra, 2000).

Materials and methods
We organized the processing in five main steps ( Figure 2) when implementing our multi-date classification. All procedures were developed in the R language (R Development Core Team, 2017), and different strategies were applied to improve the processing (e.g. parallelization, vectorization, and c-code libraries; Bengtsson, 2016;Clayden, 2016;GDAL Development Team, 2017;Hijmans, 2016;Weston, 2015) in complex computations.

Data and preprocessing
For this study, we downloaded 1350 images for four Landsat footprints (09-60, 09-61, 10-60, and 10-61) and for the period 1989-2016. They were processed to surface reflectance and acquired from the United States Geological Survey (USGS) Global Archive, sourced through the Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA; USGS, 2014). This data-set included Landsat TM, Enhanced Thematic Mapper plus (ETM+), and the Operational Imager (OLI) data. This ready-to-use data-set is radiometrically calibrated by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS; Masek et al., 2012) and orthorectified using a digital elevation model (DEM) and ground control points (NASA, 2011). The percentage of masked areas (or no-data pixels) for each image was calculated using C code based on the Function of Mask (Fmask) algorithm (CFmask; Zhu & Woodcock, 2012) and LEDAPS by merging cloud, cloud shadow, glacier, and water areas into a unique class (Figure 3(a)). Images with percentages above 90% of no-data pixels and images without orthorectification according to their metadata were omitted, reducing the time series of 27 years to effectively 23 years since 1991 (Section "Postclassification change detection"). In total, 288 images were used (i.e. 68, 79, 74, and 67 images per Landsat footprint, respectively) to complete the multitemporal composites. Data were mostly available during the dry season (67% of images, Figure 3(b)), while the fewer images acquired during the rainy season avoid additional data gaps in single years for our forest dynamics analyses (Kimes, Nelson, Skole, & Salas, 1998;Lunetta, Johnson, Lyon, & Crotwell, 2004). This is in contrast to other studies, which selected images from specific periods within a year (Müller, Griffiths, & Hostert, 2016). However, we preferred to maintain all images as a strategy to reduce data loss in this data-sparse environment. The average time interval between consecutive images for each footprint was 141, 122, 130, and 154 days, respectively. Nevertheless, in all footprints, a data gap from August 1992 to July 1996 in the Landsat archive, introduced an interruption of 3.9 years in our time series. Finally, a set of vegetation indices, band ratios, and Tasseled Cap transformation derivatives were calculated from the Landsat images (Table 1) following recommendations from similar forest dynamics studies (Kennedy et al., 2010;Müller et al., 2016;Potapov et al., 2012). To overcome the topographic effects, a c-correction algorithm (Riaño, Chuvieco, Salas, & Aguado, 2003) was applied to Landsat bands and its derivatives to evaluate whether it contributes to improving classification results (Section "Accuracy metrics results").
We collected 872 image chips from different sources acquired between April 2000 and August 2016 from high-and very high-resolution data for validating our medium-resolution remote sensing outputs (Olofsson et al., 2014): aerial photography (1-m spatial resolution), pan-sharpened images from the Advanced Land Imager (ALI, 10 m), Sentinel-2a (10 m), and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery (15 m). The latter is available for free (National Institute of Advanced Industrial Science and Technology & Geological Survey Japan, 2017), and we downloaded the whole archive to densify our high-resolution validation dataset. All images were preprocessed including manual coregistration (using Landsat imagery as a reference) and cloud masking including shadows (applying simple thresholds and manual screening in some cases). Finally, multispectral imagery were stacked, that is, ASTER, ALI, and Sentinel, for display them as falsecolor composites during the construction of validation samples plots (Section "Accuracy assessment").
To guide our implementation (Section "Postclassification change detection"), the informat ion used for the establishment of Ecuador's Forest Reference (EFR) Emission Level (MAE, 2017) was collected. This data-set constitutes a series of land-cover and vegetation maps based on Landsat, ASTER, and Rapid Eye imagery that were visually interpreted with accuracies around 70% for different periods between 1990 and 2014.
For the elevation source, we used the three arcsecond (90 m) DEM from the Shuttle Radar Topography Mission (SRTM; CGIAR -CSI, 2008). This data-set corrected for data gaps and its quality is, especially for mountainous regions of Ecuador, higher than the one arc-second (30 m) resolution product. We derived elevation, slope, aspect, roughness, the topographic position index (TPI), and the terrain ruggedness index (TRI; Wilson, O'Connell, Brown, Guinan, & Grehan, 2007) to evaluate their contribution to classification performance.

Standardized biennial compositing
For compositing, we discarded Landsat bands 1 and 2 as they are known to be more sensitive to atmospheric effects (Zhang, Carder, Muller-Karger, Lee, & Goldgof, 1999), while Landsat bands 3, 4, 5, and 7 and the calculated derivatives described in Table 1 Crist and Cicone (1984) were grouped in biennials according to the acquisition date of the image used in their calculation. The biannual time step was chosen as it resulted in 5 ± 2 images being available for most composite cases. This arrangement resulted in a no-data percentage average of 34 ± 13% (Figure 4(a)). All biannual input pixels were z-transformed for the compositing: with ρ ij being the pixel vector at position i and for date j, and μ and σ being the pixel's mean and standard deviation at date j. We then calculated the median Z i value for each vector Z i as this metric is known to be less affected by atmospheric contamination or phenological variation in image compositing (Potapov et al., 2012). Other metrics (e.g. quantiles, maximum, minimum, and variance) that are regularly applied in similar studies (De Fries, Hansen, Townshend, & Sohlberg, 1998) were also tested. However, the data-scarce situation required a conservative approach using the median. Finally, we normalized Z i and stored it as the output value β for a given biennial composite according to As residual radiometric offsets occur in the overlap areas between footprints, we required a pixel-level radiometric alignment (Pflugmacher, Cohen, & Kennedy, 2012). We selected the 2002 composite from path-row 09-60 and 09-61 as reference composites, as they had few no-data values and a low atmospheric aerosol load. Values in the overlapping footprints (10-60 and 10-61) were aligned based on histogram matching across the overlap areas to neighboring footprints. The same procedure was then applied across the time series within each footprint. In total, 52 biennial composites where aligned with references footprints, reducing, for example, the variance of the tasseled cap brightness (TCB) from 0.1 to 0.01 after histogram matching (Figure 4(b)).

Model training
We defined four classes to map permanent forest cover and deforested/reforested areas. Permanent forests included, on one hand, evergreen forests (encompassing montane, foothill, lowland, and flooded forests) and, on the other hand, Guadua spp. forests with their spectrally distinct patterns due to their different species composition, canopy height, and overall lower biomass (Silman, Ancaya, & Brinson, 2003). Other nonforest"" vegetation above 3300 m a. s.l. such as grasslands or shrubs were not considered in this research. Conversely, change classes included human land use for agricultural production, that is, pastures and croplands (including early revegetation, commercial, and subsistence plantations), and nonvegetated areas, that is, bare soils and urban areas.
We gained initial knowledge of the approximate class distribution by running a first and nonvalidated image classification with a few training samples, which built the basis for distributing the training samples for training in a guided fashion (Table 2). We then interpreted the almost 1000 training samples on-screen based on a mosaicked, biennial color composite from 2002 and on a recent vegetation map (MAE, 2017). Since the four classes represented a complex spectral feature space and their visual separation was challenging ( Figure 5(a)), we tested different classifiers using the caret package (Kuhn, 2016). This software applies a parameter tuning of classifiers and bootstraps training samples to determine their effect on performance and decided on an optimal model. As most classifiers provided similar out-of-bag error ( Figure 5(b)), we decided to use a least squares support vector machine (SVM) with polynomial kernel (svmPoly; Karatzoglou, Smola, Hornik, & Zeileis, 2004), as it achieved the highest overall correlation (0.697) with the land-cover reference maps (MAE, 2017). Other classifiers including Random Forest (rf, 0.684), Stochastic Gradient Boosting (gbm, 0.682) or Neural Networks (pcaNNet, 0.680) showed good correlations, but they were not higher than svmPoly.

Postclassification change detection
By classifying biennial composites, 13 land-cover maps were obtained for each footprint, covering the period 1989-2016. A 3 × 3 median filter was applied to eliminate spurious pixels within a land-cover map, but for data gaps, values were input calculating the per-pixel time-series mode from all land-cover maps. While random noise and discontinuities were eliminated, artifacts from clouds and cloud shadow remnants, sensor noise, or simply misclassified pixels were still present in the data. Therefore, we further applied a temporal filter with transition rules (Clark, Aide, Grau, A class-wise standard deviation of 0.9 and a confidence interval of 0.9 were considered to obtain a reasonable class-wise sample size. & Riner, 2010) to identify illogical land cover and land-cover change patterns and reclassified errors according to a set of rules based on contextual knowledge. For example, it is impossible that bare soil becomes an evergreen forest in 1 year and returns to the bare soil class the next year again. Instead, this may represent either cropped land (bare soil or cropped land in 1 year) in the case of agricultural land, or it simply may be a misclassification. In any case, it will not represent land change associated with forest cover. We accordingly implemented a temporal filter based on a moving window of three consecutive observations and a set of allowed transition rules (Table 3).
Since the first and the last land-cover maps could not be temporally filtered according to this scheme, we omitted those years (i.e. 1989-1991 and 2015-2016 periods), thereby reducing our time series to the period 1992-2014, that is, 11 observations for each footprint.
We then derived deforestation and reforestation dates from the series of land-cover maps. We flagged the first year of a forest pixel being mapped as one of the nonforest classes as the deforestation year. Reforestation, though, is a continuous process that can only be mapped from satellite data once a certain threshold of "forestedness" has transgressed, that is, a previously nonforested pixel spectrally resembles a forest class for a minimum period of time (defined in this article in 5 years as described in Section "Study area").
Finally, the outputs of the postclassification change detection were mosaicked and the UNW area extracted, considering the treeline (3300 m a.s.l.) to exclude nonforested areas. As our aim was to evaluate different algorithms for multi-date forest change classification, we iterated all these steps for each preprocessing approach (i.e. surface reflectance or topographic correction) and omitted/applied temporal filtering to evaluate their individual contribution to the overall accuracies of these maps.

Accuracy assessment
We calculated confusion error matrices for deforestation and reforestation maps (Cohen et al., 2017;Olofsson et al., 2014;Thomas et al., 2011). We employed a stratified random sampling based on the deforestation and reforestation classes, arranging samples of 5-3 pixels in a cross shape (Figure 6(a)). We sampled a minimum of 50 samples per class and 100 for the larger classes of stable forest and stable nonforest. Since reforestation stable classes could be affected during postclassification change detection due to its assigned regrowth-time threshold, we sampled these classes independently to ensure their performance. Finally, as spatial autocorrelation can bias the accuracy assessment (Congalton, 1991), a minimum threshold distance between samples of the same class was applied (Table 4).
Each sample was interpreted on-screen based on high-resolution imagery (Figure 6(b)) and Landsat color composites (Figure 6(c,d)). We calculated the overall accuracy, the kappa index, and class-wise commission and omission errors.

Variable importance and svmPoly optimization report
We decomposed the svmPoly model to observe which variables contributed to the respective results. Figure 7(a) reveals that Landsat bands 4-7 and TCB, TCW, and tasseled cap greenness (TCG) yielded an importance above 90% in all classes. Conversely, vegetation indices and band ratios as well as terrain parameters were overall less significant but contributed to separating pastures/croplands from the Guadua spp. forests. These overlapped spectrally, but separability improved when integrating terrain derivatives in the classification process. Regarding optimization, the caret software explored three parameters of the svmPoly classifier (degree, scale, and cost) to maximize its classification accuracy (Figure 7(b)). Its final calibration yielded 274 support vectors with a cost of constraints (C) of 1, and the hyperparameter values were set to a degree of 3, scale of 0.001, and a default offset of 1.

Accuracy metrics results
Deforestation and reforestation maps based on different filtering techniques varied substantially (Table 5).
Overall accuracies were significantly better when applying surface reflectance and temporal filtering to land-cover classifications, achieving 82 ± 3% and 71 ± 3% (calculated with a 95% confidence interval) for deforestation and reforestation maps, respectively. Topographic correction led to poorer overall accuracies by 13 ± 2% in both maps when compared to the     best result. Temporal filtering improved accuracies substantially by 21 ± 1% for both deforestation and reforestation maps. Commission and omission errors for stable and change classes are shown for the filtered product ( Table 6). The overall commission and omission errors were lower for the deforestation map (mean of 19% and 5%, respectively) than reforestation map (mean 32% and 11%, respectively). Moreover, stable classes were less prone to commission and omission errors (mean 1-11%) compared to change classes (mean 0-46%).

Deforestation and reforestation maps
Maps of deforestation and reforestation years are shown in Figures 8(a) and 9(a). In general, the patterns follow the description of Wasserstrom and Southgate (2013) for the Ecuadorian Amazon during its oilrelated colonization . For instance, deforested areas along the E45 highway (built in 1975) and the banks of the Napo River relate to settlements that already existed before the period we analyzed. The age of the deforestation patches along the E20 highway (built in 1983) decreased with increasing distance from the highway (Figure 8(b-1)). In the mountainous areas, the detection of landslide scars (Figure 8(b-2)) was accurate, and topographic shadows did not apparently inhibit the change detection. However, falsepositive errors were observed not only in areas with many mixed pixels where mostly evergreen forests and pasture/cropland occurred (Figure 8(b-3)) but also in areas with no-data values due to the occurrence of sparse observations or an inaccurate water mask.
Areas of reforestation seemed to be more prominent along the E45 highway, where deforestation was less intense. Known areas of reforestation since the 1990s were well represented (Figure 9(b-1)), as was forest succession after landslides in mountainous areas (Figure 9(b-2)). Overall, the reforestation year map was more affected by mixed pixel problems and  mask errors than the map of deforestation year (Figure 9(b-3)). Following Rudel, Bates, and Machinguiashi (2002), we calculated overall deforestation and reforestation by applying a buffer distance of 3 km along the two main highways E45 and E20 in the UNW to corroborate our observations. Accumulated deforestation along highway E45 summed up to an area of 3320 ha and 11,403 ha along highway E20. In contrast, reforestation along highway E45 accumulated to 7458 ha and 5415 ha along highway E20.

Comparison with other sources
We compared our implementation with two different sources: Forest Loss Year (FLY) according to Hansen et al. (2013) and EFR Emission Level information (MAE, 2017). Both sources were cropped with the UNW and relabeled to match our classes ( Figure 10).
Results are similar across the three classifications. However, differences specifically exist with FLY for specific time periods such as 2000-2002 and 2010-2012 (Figure 10(a)) or EFR reforestation between 2008 and 2014 ( Figure 10(b)). On average, deforestation was 2757 ha year −1 for the period from 2000 to 2014, in FLY data and 4394 ha year −1 in EFR. Our estimates are comparably conservative with 2319 ha year −1 . According to FAO (Puyravaud, 2003), these values represented annual deforestation rates of −0.35%, −0.57%, and −0.31%, respectively.
Furthermore, reforestation summed up to 574 ha year −1 for the 2000-2014 period in FLY, indicated 2277 ha year −1 in EFR, and 1504 ha year −1 in our analysis, representing annual reforestation rates of  0. 07%,0.28%,and 0.19%,respectively,in EFR,796,982 ha stayed unchanged, while our analysis yielded 748,688 ha of stable forests.

Biennial image compositing and preprocessing effects on results
Our image compositing technique was based on the standardization and median calculation, which is an effective strategy to maximize information extraction when the number of observations is limited. We chose a biennial classification scheme (Griffiths et al., in review). We thereby improved the signalto-noise ratio, as 1-year composites may be inferior in data-scarce conditions (Potapov, Turubanova, & Hansen, 2011). Composites from longer time periods, though, may not be adequate to monitor subtle processes such as reforestation (Bustamante et al., 2016).
The histogram matching algorithm that we used for radiometrically aligning the composites enabled a regional-scale classification and at the same time created consistency across the time series. This was supported by the high consistency between the land-cover classifications from different years and the results after applying our postclassification change detection algorithm. However, cloud-free composites as references and sufficient spatial overlap between the target and reference footprints are mandatory for the proper functioning of the histogram matching (Benjamin & Leutner, 2017).
We also accommodated for correcting the radiometric distortions due to topography. In our case, the c-correction algorithm principally improved the homogeneity of the imagery across sunlit and shaded slopes. However, commission errors increased after applying the topography correction. This is in line with the findings in Chance et al. (2016), which reported negative effects of a topographic correction on change detection analysis. Others found that the application of topography correction generally had a smaller influence on the overall accuracy of a classification when compared to the selection of a classifier (Vanonckelen, Lhermitte, & Rompaey, 2015). Future work should further improve the results from topographic correction by employing the best DEMs available (Chance et al., 2016;Pimple et al., 2017).

Postclassification change detection performance
Our postclassification change detection strategy was based on land-cover maps (MAE, 2017) as reference to validate model training and classifier outputs. This allowed the selection of the most precise classifier based on the correlation between the classification and the land-cover maps. While the limited size of our sampling set may not be representative for some classifiers (Zhu et al., 2016), SVMs apparently performed well, as SVMs support small training samples (Wieland, Torres, Pittore, & Benito, 2016).
The original Landsat bands and derived tasseled cap components had a considerable predictive power (Figure 7). This was specifically true for bands 3-7, which are known to be important predictor variables not only in forest/nonforest classifications in the tropics (Potapov et al., 2012) but also in dry regions of the world (Mellor, Haywood, Stone, & Jones, 2013). Spectral mixtures and spectral similarity of landcover types (see Figure 5(a)) limit the separability at 30-m Landsat spatial resolution. In this regard, elevation and terrain derivatives from the DEM (slope, aspect) contributed to class separation, despite their predictive power not being as high as that of the Landsat bands or tasseled cap components.
While some gaps related to cloud and cloud shadow remnants remained after the classification, we were able to demonstrate that temporal filtering is a powerful technique for removing these artifacts and considerably improving the results (comparison in Table 5). The set of transition rules allowed us to filter most of the illogical class transitions; however, some highly dynamic events were still missed due to the data-scarce setting and accordingly introduced omission errors.

Performance of multi-date classification in the UNW
Despite the remaining limitations of our multi-date classification implementation, the spatially explicit forest dynamics patterns at the UNW allow for novel insights beyond what was already known from previous satellite data analyses relying on only two or just a few points in time (Sierra, 2000) or only spectral information (Walsh, Shao, Mena, & McCleary, 2008) ideally in the Tropical Andes. Dynamics along the E45 highway after 1992 mostly related to reforestation on peripheral lands, while deforestation rates were comparably low in that region (Figure 11(a)). This may be explained by the population census, where the population in the urban centers of Tena and Archidona increased by 233% between 1990233% between and 2010233% between (INEC, 2010, suggesting a rearrangement in the population distribution between the rural and the urban areas. This assumption is supported by similar findings by Rudel et al. (2002) in the Southern Ecuadorian Amazon. In contrast, deforestation was principally identified along the E20 highway. Since this highway was constructed more recently, new settlements and commercial activities linked to oil extraction have triggered deforestation (Wasserstrom & Southgate, 2013; Figure 11(b)).
A comparison with other sources revealed further details. For instance, deforestation showed to be more similar to FLY than to EFR, most likely because the FLY data-set is also based on a multi-date classification, while the EFR was based on an object-based classification that generalized deforestation patches. In the case of reforestation, all results differed markedly. Different conceptualizations of reforestation (Hansen et al., 2013;MAE, 2017) and confusion with secondary forests are likely to be the main reasons. According to Cohen et al. (2017), it is not surprising that forest disturbance maps differ due to semantic and methods differences. Different accuracies and results accordingly relate to multiple factors such as differences in the change detection algorithm, in the quality of specific satellite imagery used, metrics and training information, the time-series density, or the thresholds applied to identify change.
Overall, our multi-date classification implementation was demonstrated to be far less sensitive to data scarcity and atmospheric contamination than other approaches using automated time-series analysis algorithms (Santos et al., 2017).

Conclusions and outlook
Forest dynamics in the complex and vulnerable regions of the Tropical Andes are still underresearched considering remote sensing data analyses (Da Ponte et al., 2015;Oliveras, Anderson, & Malhi, 2014). To the best of our knowledge, this was the first study of its kind that specifically focused on the challenges related to scarce data and the poor signal-to-noise ratio in a long time series for automated forest change analyses in the Tropical Andes. We demonstrated that an adapted implementation of multi-date classification based on image compositing, multi-sensor fusion, and postclassification change detection could mitigate most of these limitations. Our findings add to the expanding body of literature on such approaches with a focus on data-scarce situations and highlight the importance of the Landsat archive for monitoring decadal land-cover change even in cloudy regions of the world.
Future research should focus on diversifying data sources and predictors as our findings provide further evidence that classification results, specifically when using machine learners, will improve in data-rich environments. Moreover, the increasing web-based availability of high-and very high-resolution data will in the future allow to further improve sample quantity and quality, while semi-automatic approaches (Huang, Weng, Lu, Feng, & Zhang, 2015) and temporal-spectral profiles sampling (Senf, Pflugmacher, Wulder, & Hostert, 2015) are also promising alternatives. Furthermore, since our methodology requires a reforestation time threshold, it would be beneficial considering specific thresholds for different forest communities. This should ideally be based on forest growth models such as FORMIND (Paulick, Dislich, Homeier, Fischer, & Huth, 2017) that support specifying distributions for reforestation time threshold. Additionally, improvements may well be possible with further refined transitioning rules in postclassification filtering or automated solutions with more complex transition rules with an increasing number of land-cover classes (Ahlqvist, 2008;Abercrombie & Friedl, 2016).
As other areas may experience similar or even more severe data scarcity than the UNW, image compositing might be limited to lower observation frequencies. In this regard, regions such as north-central Africa and northern Russia, which have the sparsest Landsat coverage compared to Ecuador , may constrain multi-date classification usability to frequencies greater than biennials. Such limitations may, for example, constrain its usability in the frame of reducing emissions from deforestation and forest degradation (REDD+), which requires biennial updates reports for Forest Reference Level Information especially in developing countries (UN- REDD Programme, 2015). Finally, Figure 11. Deforestation/reforestation area for the period 1992-2014 (a) along the highway E45 and (b) along the highway E20.
implementing such a multi-date classification for larger study areas requires cloud-based or high-performance computing (HPC) environments, as the processing is demanding, and it is more effective to "bring the algorithm to the data" than to download massive data-sets. Currently, some alternatives are available (e.g. EODC, 2018;Gorelick et al., 2017;Open Foris, 2015), which allow implementing similar methodologies for large areas. Cloud-based or HPC environments also provide novel opportunities to develop monitoring systems based on sensor constellations, such as Landsat and Sentinel-2 .
As optical remote sensing of the core tropics regularly suffers from high cloud cover, integrating newly available imagery will increase change map reliabilities. Linking the vast Landsat archive with the quickly expanding Sentinel-2 archive is, therefore, one of the cornerstones for future improvements (Drusch et al., 2012;Wulder et al., 2016). Such a strategy will also allow to extend the applicability of our approach to larger regions such as the entire Tropical Andes and to ecosystems with more diverse land cover.