Detecting and Attributing Drivers of Forest Disturbance in the Colombian Andes Using Landsat Time-Series

The Colombian Andes foothills have seen an expansion of forest disturbance since the 1950s. While understanding the drivers of disturbance is important for quantifying the implications of land use change on regional biodiversity, methods for attributing disturbance to specific drivers of change at a high temporal and spatial resolution are still lacking in the Andes region, in part due to persistent cloud cover. Using 20 years of Landsat images (1996–2015) covering Picachos National Park in the Colombian Andes, we detected sub-annual forest cover disturbances using the Breaks For Additive Season and Trend (BFAST) Monitor algorithm; characterized different types of disturbance using spectral, spatial, and topographic indicators; and attributed causes of forest disturbance such as conversion to pasture, conversion to agriculture, and non-stand replacing disturbance (i.e., thinning) using a Random Forest (RF) classifier. Conversion to pasture has been the main driver of forest disturbance in Picachos, responsible for 11,395 ± 72 ha (17%) of forest cover loss, followed by non-stand replacing disturbance and conversion to agriculture. Disturbance detection had 96% overall agreement with validation data, although we had a high omission error of 21% primarily associated with forest to agriculture conversion. Other change drivers had a much more reliable attribution with forest to pasture conversion or non-stand-replacing disturbance, showing only 1–5% commission and 2–14% omission errors. Our results provide spatially-explicit information on sub-annual disturbances and associated drivers of change that are necessary for evaluating and improving domestic conservation efforts and establishing systematic ecological observations, which is currently absent from Colombia. While effective at revealing forest change dynamics in a geographically remote and socio-politically complex region like Picachos, our approach is highly automated and it can be easily extended to the rest of Andes-Amazon transition belt where low availability of remote sensing data and high cloud cover impede efforts at consistent monitoring of forest cover change dynamics and drivers.


Introduction
Colombia is a "megadiverse" country given its high level of endemism and species richness [1].However, over the last 50 years, the spread of informal settlements into forested regions, internal social conflict, and the growing demand for agricultural land have increased the pressure on Colombia's forest ecosystems even in remote or protected regions [2].As a consequence, between 1990 and 2015, Colombia lost more than six million hectares of forest, and more significantly, witnessed a 44% increase in deforestation in 2016 compared to 2015 [3].Despite the duration and geographic breadth of these changes, there remains limited understanding of the relative influence and timing of different drivers on forest cover change across much of Colombia.The Colombian Andes foothills were rather inaccessible three decades ago but have seen recent forest disturbance [4][5][6] due to favorable incentives for cattle ranching in the 1980s [7], the enlargement of coca cultivation in the 1990s through 2002 [2], and the more recent expansion of cattle ranching and other agricultural land uses [8,9].Though National Protected Areas make up 12% of Colombia's land mass [10], protected areas are also threatened by migration and displacement, economic development, and illicit cultivation of coca (Erythroxylum coca Lam.) [11], with attendant effects on socio-ecological relationships [4,7,[12][13][14].
Despite many methods having been recently developed for change detection [15], there remains a lack of studies that identify drivers of disturbances with high spatial and temporal detail in rich biodiversity regions, and the Colombian Andes-Amazon area is no exception.Previous studies have identified drivers of Colombian forest cover change at the sub-national level through comparisons of two or three Landsat images collected five or more years apart [2,10].However, a bi-or tri-temporal change assessment limits our ability to monitor more temporally complex or spectrally subtle land cover change processes that can only be detected by annual or sub-annual time-series analysis [16].Other studies have used more regularly collected Moderate Resolution Imaging Spectroradiometer images for measuring the long-term effectiveness of protected areas at limiting deforestation [17,18] or detecting the main drivers of forest cover change at the municipality-level [12].However, those analyses have incorporated coarse spatial resolution data (>250 m) that are unable to capture small-scale processes of change common in mountains and inaccessible places [19,20].Because of these limitations, such studies have not been able to gauge the effect of drivers of forest cover change with adequate spatial or temporal detail to support regional forest conservation efforts.
Remote sensing time-series analysis that exploits the full archive of available imagery offers the means of capturing disturbance processes over their entire duration, whether short-or long-term.A diversity of time-series algorithms has been developed to monitor forest disturbance including LandTrendr [21], the vegetation change tracker (VCT) [22], and trend detection approaches based on "best-available-pixel" composites [23][24][25].However, these approaches produce annual disturbance maps, which may miss short-lived or otherwise sub-annual changes [26] such as post-disturbance recovery or changes in condition (low impact selective logging, small induced fires, and foliage condition) that are better captured using dense time-series [19,27,28].Further, while such algorithms have been reliable in regions with relatively low cloud cover and a relatively uninterrupted Landsat archive of >30 years (e.g., United States and Canada; [21,23,24,29], these disturbance detection algorithms are broadly untested in regions with regular cloud cover or satellite image scarcity. The development of algorithms such as the Noise Insensitive Trajectory Algorithm [30], the Detecting Breakpoints and Estimating Segments in Trend [31], and the Breaks For Additive Season and Trend (BFAST)-Monitor [28,32] have improved the capacity for disturbance detection despite persistent cloud cover by building a dense time-series stack of all possible dates of satellite observation.Of these, BFAST-Monitor has had the widest adoption for identifying disturbance events, having been used to model seasonal and trend components in tropical forested regions [33] with a single vegetation index (e.g., NDVI) [34,35], multiple vegetation indices [36], or by taking into account spatial context [37,38].
The detection of forest disturbance events and processes has been highly operationalized at the global scale [39], yet the identification of proximate drivers of forest disturbance remains a challenge.Most driver detection studies have focused on North American forests and used spectral-temporal metrics to map common drivers such as fires, clear-cut harvests, insect outbreaks, and road expansion [21,23,24].While these approaches have been shown to be effective in data rich In this study, our objective is the detection and attribution of drivers of forest disturbance in an area with abrupt and subtle gradual changes using Landsat time-series (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).Using the foothills of Picachos National Park as a case of study, we (1) detected forest disturbances using BFAST-Monitor and (2) attributed the disturbance driver using a Random forest classifier and various time-series metrics.We focused on the most common drivers of disturbance, including conversion to pasture, conversion to agriculture, and non-stand replacing disturbance.These drivers were selected given their recognized negative impact on landscape fragmentation, forest structure, and biodiversity [14,40].Drivers of forest disturbance were characterized using metrics derived from time-series trajectories such as spectral metrics (e.g., average spectral values before and after disturbance event), trend metrics related to magnitude of change, spatial pattern of disturbance (e.g., patch area and perimeter), and topographic metrics (e.g., elevation).Using previous metrics and training data, we identified the most probable driver of disturbance.The results contribute to the adoption of conservation monitoring strategies (e.g., zonification processes) in tropical forested areas that are threatened by unregulated and often undocumented human encroachment, which are difficult to monitor due to regular cloud cover and smaller, although progressive, disturbance patches.

Study Area
The study covers 67,284 ha in the foothills of Colombia's Picachos National Park and includes the basins of the Templado, Platanillo, and Yulo rivers (Figure 1).Picachos is quite rugged with an elevational range of 300-2500 masl, has a temperature range of 5-25 • C, and annual accumulated rainfall varying from 1050-3250 mm [14].The Park has great hydrological significance, as river headwaters here feed the basins of the Orinoco and Amazon rivers.A majority (62%) of the Picachos' foothills are part of the Andes-Orinoquia biome characterized by broad-leaf evergreen forest, with additional land covers including pasture and mosaic cropland (shrub, herbaceous cover) [14].
Forests 2018, 9, x FOR PEER REVIEW 3 of 17 [21,23,24].While these approaches have been shown to be effective in data rich contexts, they have not been broadly tested in other regions complicated by fewer satellite observations, persistent cloud cover, or smaller disturbance patches.In this study, our objective is the detection and attribution of drivers of forest disturbance in an area with abrupt and subtle gradual changes using Landsat time-series (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).Using the foothills of Picachos National Park as a case of study, we (1) detected forest disturbances using BFAST-Monitor and (2) attributed the disturbance driver using a Random forest classifier and various time-series metrics.We focused on the most common drivers of disturbance, including conversion to pasture, conversion to agriculture, and non-stand replacing disturbance.These drivers were selected given their recognized negative impact on landscape fragmentation, forest structure, and biodiversity [14,40].Drivers of forest disturbance were characterized using metrics derived from time-series trajectories such as spectral metrics (e.g., average spectral values before and after disturbance event), trend metrics related to magnitude of change, spatial pattern of disturbance (e.g., patch area and perimeter), and topographic metrics (e.g., elevation).Using previous metrics and training data, we identified the most probable driver of disturbance.The results contribute to the adoption of conservation monitoring strategies (e.g., zonification processes) in tropical forested areas that are threatened by unregulated and often undocumented human encroachment, which are difficult to monitor due to regular cloud cover and smaller, although progressive, disturbance patches.

Study Area
The study covers 67,284 ha in the foothills of Colombia's Picachos National Park and includes the basins of the Templado, Platanillo, and Yulo rivers (Figure 1).Picachos is quite rugged with an elevational range of 300-2500 masl, has a temperature range of 5-25 °C, and annual accumulated rainfall varying from 1050-3250 mm [14].The Park has great hydrological significance, as river headwaters here feed the basins of the Orinoco and Amazon rivers.A majority (62%) of the Picachos' foothills are part of the Andes-Orinoquia biome characterized by broad-leaf evergreen forest, with additional land covers including pasture and mosaic cropland (shrub, herbaceous cover) [14].Picachos has one of the largest, relatively intact regions of forest cover within Colombia, connecting the Orinoquia, Amazon, and Andes ecosystems [14], and was designated as a Category II Park by the International Union for Conservation of Nature (IUCN) in 1977 to conserve large-scale ecological processes, species, and ecosystems for current and future generations [41].Although people are not legally allowed to reside within the Park, more than 200 families and 119 small houses and plots of 50-250 ha were documented in Picachos as of 2016 [14].The settlement of peasant communities results from historical processes of colonization and forced displacement and has expanded present-day grazing land and subsistence agricultural production [8,14,42,43].

Image Processing
Imagery from the Landsat archive was selected given its long record of regular observation and intra-annual observation frequency of nearly eight images per year over Picachos, which supports sub-annual forest disturbance [6]; and its 30 m (0.09 ha) spatial resolution, which is suitable for detecting small disturbance patches common in Picachos (mean =2 ha).We used a time-series of 149 Level-1 Terrain-Corrected (L1T) Landsat 5, 7, and 8 surface reflectance images (path/row 8/58) from 1996 to 2015 with cloud cover less than or equal to 50%; 225 additional images with more than 50% cloud cover were removed prior to analysis.Clouds and cloud shadows, water, and snow were masked using the CFmask cloud-shadow mask product [44].Cloud confidence and pixel description from pixel Quality Assessment (QA) band's quality were used to remove medium and high cloud confidence, but also pixels surrounding high cloud confidence.
For each image in the time-series, we calculated NDMI (Normalized Difference Moisture Index; Equation ( 1)) [45], which combines short-wave infrared (SWIR) and near infrared (NIR) to measure canopy moisture content [46].For the NDMI calculation, spectral bands B4 and B5 for Landsat 5 and 7 were used, whereas B5 and B6 were employed for the Landsat 8 case.
NDMI has been successfully used to monitor disturbance and regrowth in tropical forests [47], has shown better sensitivity for disturbance detection than greenness indices such as Enhanced Vegetation Index (EVI) and the Normalized Difference Vegetation Index (NDVI) [48], and has provided a similarly high disturbance detection accuracy in tropical forests compared to Tasseled Cap Wetness (TCW) [49] or the Normalized Difference Fraction Index (NDFI) [33,50].Additionally, NDMI is effective even with the level of intra-annual variability in image availability common in the Andes [33,34].With NDMI, we produced a binary forest/non-forest mask using a Support Vector Machine algorithm in ENVI 5.1 on a cloud-free image on 22 February 1998.This mask represents all forested pixels at the beginning of the study period and was applied to all subsequent image dates to support the identification of forest disturbances, excluding other land cover dynamics.

Disturbance Detection Using BFAST-Monitor
Given its proven robustness for monitoring changes across various tropical forests [34,35,38], we selected the BFAST-Monitor algorithm to map disturbances in Picachos.BFAST-Monitor iteratively decomposes time-series data into trend, seasonal, and noise components, and detects significant breakpoints where the difference between the actual and expected signal deviates beyond a 95% significance boundary [32,51].
We applied BFAST-Monitor in one-year sequential periods starting in 1999 using a first-order harmonic curve over the NDMI time-series; the 1996-1998 period was used as a historical baseline.Due to the irregularity of clear-sky Landsat observations, the trend component was omitted to avoid erroneous NDMI values in the monitoring period that may cause the detection of false breakpoints [34].Although the breakpoint detection approach is robust to noise and other variability in the time-series [32], a threshold based on NDMI change magnitude was applied to avoid the detection of spurious breakpoints due to outliers in the monitoring period.To establish this threshold, we visually compared potential breakpoints with a WorldView-1 true color (0.5 m resolution) image from 12 December 2009, and selected an NDMI change magnitude threshold of −0.01.Therefore, all breakpoints with an NDMI decline lower in magnitude than this threshold were classified as disturbed.Disturbances identified by BFAST-Monitor involve three components: location, change date, and change magnitude.Individual disturbance pixels were grouped into patches based on Rook's case adjacency; patches with an area less than 1 ha (~9 pixels) were discarded to avoid isolated pixels and the salt-pepper effect.Forest pixels that did not show any disturbance over our study period were considered stable forest (non-disturbed).

Disturbance Validation Using TimeSync
An accuracy assessment identifies the errors during the classification of disturbance, and the sample data can be used for estimating both accuracy and area along with the uncertainty of these estimates [52].We used TimeSync [34,53], an interactive tool that supports the visual interpretation of forest cover change across the time-series, to validate the disturbance timing and location.We sampled 614 sites from across Picachos using a stratified sampling procedure following Cochran's recommendations [54]; twenty-two sites were removed because of excessive cloud cover in either Landsat or high-resolution Google Maps-hosted reference.Within TimeSync and referencing Google Earth and RapidEye high-resolution images, we interpreted annual mean Landsat NDVI and NDMI composite time-series to identify disturbance dates.Of our 592 sample sites, we interpreted 211 disturbed sites and 381 non-disturbed sites.The initial area estimates were derived directly from a disturbed/non-disturbed map (i.e., pixel counting), and the area-adjusted error was measured using a stratified estimator for every class to remove the biased effect of unequal omission and commission errors (see Supplementary Material).For disturbances, a confusion matrix was created by cross-tabulating the BFAST Monitor-derived disturbance information with TimeSync-derived validation data.For drivers of disturbance, a confusion matrix from the Random Forest output was used.For both matrices, total, producer's, and user's accuracies were calculated.

Characterizing Drivers of Change
Disturbance detection offers fundamental information about changes in landscape form and function that result from drivers of change, which differ in their spectral-temporal response function [19,55].Long-term, high magnitude disturbances are commonly related to land cover conversion while short-term, low magnitude disturbances are more commonly a result of land cover alteration (Figure 2).For example, conversion to agriculture might present an abrupt change followed by an annual cycling of spectral values indicative of vegetation harvest and growth.The spectral trend of non-stand replacing disturbance, on the other hand, has a subtle but prolonged decrease in condition, and conversion to pasture commonly presents an abrupt break followed by a sustained period of relatively low vegetative condition (e.g., Figure 3).
To effectively characterize each driver, we measured patch-level metrics based on spectral indicators, trend, spatial pattern, and topographic condition (Table 1).First, spectral metrics consisted of average and standard deviation of near infrared (NIR), shortwave infrared 1 (SWIR1), shortwave infrared 2 (SWIR2), and NDMI reflectance, before and after the date of disturbance [23,56].Second, trend metrics were calculated using the mean NDVI trend at each disturbance pixel for each year using ordinary least squares (OLS) regression [26]; we also included the kurtosis and skewness for both the magnitude of change and NDVI trend [24].Pixel-level trends and magnitudes were then averaged across each patch.Third, for the spatial pattern, we measured the patch area, perimeter [57], shape index [58], and fractal dimension [59].Finally, for topographic condition, we measured the elevation, slope, and topographic position index defined as the difference between the elevation of a cell and the mean value of its eight surrounding cells [21].

Attributing Drivers of Change
From the total population of 3450 disturbance patches, we visually interpreted and categorized 576 (17%) patches as follows: conversion to agriculture (50 patches, 146 ha total); non-stand replacing disturbance (147 patches, 361 ha total); and conversion to pasture (379 patches, 1171 ha total).This categorization was carried out using land cover maps based on Corine Land Cover (CLC) data generated by Colombian government agencies for 2002, 2007, and 2012 [8,14], which provide the most reliable reference data for interpreting land cover conversion in Colombia.Additionally, we used high resolution (0. Patch-level metrics (Table 1) were used for the training and validation (60% and 40% of sampled patches, respectively) of a Random Forests (RF) classification of drivers for each visually interpreted disturbed patch [21,23,36,60].With 36 input metrics, we limited the potential for multicollinearity by removing nine metrics with Pearson's r values greater than 0.70.Drivers attributed to each disturbed patch using majority votes from the RF were compared with the TimeSync validation sample to measure the attribution accuracy.We also calculated the distance to the second most voted class as a measure of classification confidence, calculated as the ratio d2/d1, where d2 corresponds to the second most voted class and d1 is the assigned class [61].Resulting values are between 0 (when d2 < d1) and 1 (when d2 = d1), and patches with classification confidence values higher than 0.6 were re-labeled as unclassified given the unreliable driver attribution [23].

Disturbance Detection Agreement
The total disturbed forest area within the Picachos foothills from 1998-2015 was 14,213 ± 1095 ha, while stable forest was recorded across 53,076 ± 1094 ha.Our estimate of disturbed forest within Picachos is 14% more than the 12,128 ha of forest cover loss estimated by official reports based on CLC maps during the similar study period of 1999-2012.Using visual interpretation results at 592 validation pixels, we estimated a very high total accuracy (TA) of 96% ± 1 (see pixel-level confusion matrix and full area-weight matrix Tables S1 and S2 in supplementary material).The disturbance class had user's (UA) and producer's accuracies (PA) of 96% ± 2 and 79% ± 8, respectively; the stable forest class had UA of 94% ± 2 and PA of 99% ± 0.4.Disturbed pixels had low commission (4%) and moderate omission errors (21%), whereas stable forest pixels were more affected by commission (6%) and less so by omission errors (1%) (Table 2).The higher rate of omission for the disturbance class is a consequence of low image availability (e.g., average of 51 cloud-free Landsat observations in the period 1999-2015), which is much lower compared to other tropical forest disturbance studies that have used BFAST-Monitor [33,47].

Characterization of Drivers'
Overall, spectral and trend metrics had the most predictive power, followed by topographic and spatial pattern indicators.The mean value post-change NDMI was the most important metric in the RF classification, followed by post-change standard deviation SWIR2, NDVI trend, and mean pre-change NDMI (Figure 4).If these metrics were excluded from the RF model, the accuracy would decrease by 33%, 21%, 15%, and 13%, respectively.Other metrics such as mean post-change NIR (11%), standard deviation of NDVI trend (8%), elevation (7.5%), and median magnitude of NDMI (4.6%) had much lower respective importance.and less so by omission errors (1%) (Table 2).The higher rate of omission for the disturbance class is a consequence of low image availability (e.g., average of 51 cloud-free Landsat observations in the period 1999-2015), which is much lower compared to other tropical forest disturbance studies that have used BFAST-Monitor [33,47].

Characterization of Drivers'
Overall, spectral and trend metrics had the most predictive power, followed by topographic and spatial pattern indicators.The mean value post-change NDMI was the most important metric in the RF classification, followed by post-change standard deviation SWIR2, NDVI trend, and mean prechange NDMI (Figure 4).If these metrics were excluded from the RF model, the accuracy would decrease by 33%, 21%, 15%, and 13%, respectively.Other metrics such as mean post-change NIR (11%), standard deviation of NDVI trend (8%), elevation (7.5%), and median magnitude of NDMI (4.6%) had much lower respective importance.The most important metrics for characterizing drivers were those related to post-disturbance NDMI and SWIR2, which reflects the variation in post-disturbance spectral trends across drivers.Spatial pattern metrics such as patch area or perimeter displayed little contribution to the The most important metrics for characterizing drivers were those related to post-disturbance NDMI and SWIR2, which reflects the variation in post-disturbance spectral trends across drivers.Spatial pattern metrics such as patch area or perimeter displayed little contribution to the classification, likely because the area of a disturbed patch was similar between drivers, and larger patches affected by other drivers agents (e.g., wildfires) are not common in the study area [10].Drivers had varying degrees of separability across metrics.Considering the six most important metrics, each driver showed considerable variation in median and interquartile ranges (Figure 5).For less important metrics, there was greater overlap in the interquartile range.For example, non-stand replacing and conversion to agriculture are distributed in a wider elevation range than conversion to pasture, which is mostly located at lower elevations.Conversion to pasture was well distinguished from the other two drivers across metrics, while conversion to agriculture and non-stand replacing disturbance showed less spectral separability.Of note, non-stand replacing disturbance had the smallest standard deviation in NDVI, indicating the subtle spectral changes generated by forest thinning over time.Conversion to agriculture exhibited narrowed values, especially for NDVI trend and pre-change NDMI, with the latter indicating stable values before disturbance detection.
Forests 2018, 9, x FOR PEER REVIEW 9 of 17 classification, likely because the area of a disturbed patch was similar between drivers, and larger patches affected by other drivers agents (e.g., wildfires) are not common in the study area [10].Drivers had varying degrees of separability across metrics.Considering the six most important metrics, each driver showed considerable variation in median and interquartile ranges (Figure 5).For less important metrics, there was greater overlap in the interquartile range.For example, non-stand replacing and conversion to agriculture are distributed in a wider elevation range than conversion to pasture, which is mostly located at lower elevations.Conversion to pasture was well distinguished from the other two drivers across metrics, while conversion to agriculture and non-stand replacing disturbance showed less spectral separability.Of note, non-stand replacing disturbance had the smallest standard deviation in NDVI, indicating the subtle spectral changes generated by forest thinning over time.Conversion to agriculture exhibited narrowed values, especially for NDVI trend and pre-change NDMI, with the latter indicating stable values before disturbance detection.

Driver Attribution Agreement
The total agreement of driver attribution at the weight-area level was 98%±0.008using a 95% level of confidence (Table 3) (patch -level confusion matrix and full area-weight matrix are in Supplementary Tables S3 and S4, respectively).Based on area-weighted methods, conversion to pasture exhibited the lowest omission and commission (1%) errors and was most readily discriminated from non-stand replacing disturbance and conversion to agriculture.Non-stand replacing disturbance had 5% and 12% omission and commission errors, respectively, and the lowest agreement was found in the conversion to agriculture class, which exhibited high omission (30%) and commission (4%) errors.

Driver Attribution Agreement
The total agreement of driver attribution at the weight-area level was 98%±0.008using a 95% level of confidence (Table 3) (patch -level confusion matrix and full area-weight matrix are in Supplementary Tables S3 and S4, respectively).Based on area-weighted methods, conversion to pasture exhibited the lowest omission and commission (1%) errors and was most readily discriminated from non-stand replacing disturbance and conversion to agriculture.Non-stand replacing disturbance had 5% and 12% omission and commission errors, respectively, and the lowest agreement was found in the conversion to agriculture class, which exhibited high omission (30%) and commission (4%) errors.Conversion to agriculture was commonly confused with non-stand replacing disturbance, likely due to spectral similarity in agricultural cycling and subtle short-and long-term changes in the forest (Figure 3).Making a clear distinction between both classes is difficult since agriculture has a relatively low coverage area and is generally surrounded by other land covers, which can introduce Forests 2018, 9, 269 10 of 16 a mixed pixel effect, reducing the potential for accurate classification.Moreover, despite non-stand replacing disturbance often preceding deforestation driven by subsistence agriculture [16,36] or pasture expansion [62,63], we were unable to identify causal linkages between specific drivers.

Driver Dynamics
The annual area of disturbance increased from 1999 through 2007, except for 2005-2006 when the rate of disturbance was relatively low and 2003 when imagery was not available (Figure 6).This increase in disturbed area was mainly the result of pasture expansion for cattle grazing by people living within Picachos.Much of this conversion took place in favorable flat and productive soil along the Templado, Platanillo, and Yulo rivers, and was likely stimulated by economic incentives conferred by the Colombian government to cattle ranchers to expand their production in the 1980s and 1990s [14].After 2007, the rate of annual disturbance had a more inconsistent trajectory.Of note, the increase in disturbed area in 2007 is not a product of increased Landsat image availability, since other years within our temporal study window had a similar or higher number of images.
Forests 2018, 9, x FOR PEER REVIEW 10 of 17 Conversion to agriculture was commonly confused with non-stand replacing disturbance, likely due to spectral similarity in agricultural cycling and subtle short-and long-term changes in the forest (Figure 3).Making a clear distinction between both classes is difficult since agriculture has a relatively low coverage area and is generally surrounded by other land covers, which can introduce a mixed pixel effect, reducing the potential for accurate classification.Moreover, despite non-stand replacing disturbance often preceding deforestation driven by subsistence agriculture [16,36] or pasture expansion [62,63], we were unable to identify causal linkages between specific drivers.

Driver Dynamics
The annual area of disturbance increased from 1999 through 2007, except for 2005-2006 when the rate of disturbance was relatively low and 2003 when imagery was not available (Figure 6).This increase in disturbed area was mainly the result of pasture expansion for cattle grazing by people living within Picachos.Much of this conversion took place in favorable flat and productive soil along the Templado, Platanillo, and Yulo rivers, and was likely stimulated by economic incentives conferred by the Colombian government to cattle ranchers to expand their production in the 1980s and 1990s [14].After 2007, the rate of annual disturbance had a more inconsistent trajectory.Of note, the increase in disturbed area in 2007 is not a product of increased Landsat image availability, since other years within our temporal study window had a similar or higher number of images.Conversion to pasture led to the greatest area of disturbed land (17% of Picachos foothills; 11,425 ± 59 ha), which is slightly more than the CLC-calculated area of 10,827 ha during the similar period of 1999-2012.The second most dominant driver was non-stand replacing disturbance with an area of 1085 ± 76 ha; there is no documented precedent for this driver in CLC assessments.Conversion to agriculture contributed 290 ± 50 ha, which is also in good agreement with the CLC estimate of 380 ha.Finally, 1175 ha, corresponding to 8% of the total disturbed area, could not be confidently attributed and was left as unclassified.
Disturbance was most prominent along the margins of rivers with the Platanillo watershed experiencing the most disturbance; the earliest disturbances took place close to these rivers and have since expanded outwards (Figure 7).Conversion to pasture was found along the Templado, Plantanillo, and Yulo rivers, while conversion to agriculture and unclassified patches were mainly found near Yulo River.As Yulo River forms the southern border of Picachos, disturbance here likely reflects the influence of agricultural expansion and grazing from nearby peasant reserves.On the other hand, non-stand replacing disturbance was more dispersed over the entire study area and even evident at higher elevations.Notably, conversion to pasture and non-stand replacing disturbance are commonly in close proximity to each other, suggesting a gradual land conversion process.Though Picachos is a strategic regional corridor between the Andes and Amazonia, there remains limited biophysical and sociocultural information to support effective decision-making processes.Detecting patterns of pasture intensification with remote sensing-derived spatio-temporal information is a critical tool in designing and evaluating land management policies that support desired, sustainable land use strategies [67].The results of this study record in great detail the short-

Discussion
This study is the most comprehensive analysis of disturbances and drivers in the Andes, and highlights the value of using the full time-series of Landsat observations within an automated disturbance detection and driver attribution framework.Given the high spatio-temporal resolution, this approach was able to map land cover dynamics that are broadly absent from Colombia's Corine Land Cover mapping framework, which records changes on a five-year timescale.Though the current study was constrained to the Picachos foothills, the disturbances and drivers in Picachos are typical across the rest of the Andes-Amazonia transition belt region [7,64], especially pasture expansion and intensification [65].In Picachos, as well as the rest of the Andes-Amazon region, conversion to pasture is highly related to cattle expansion [40] and has detrimental effects on landslide susceptibility in high relief regions [7,66], as well as ecosystem connectivity between the Andes and Orinoco and Amazon ecosystems [2].Given the common challenges to disturbance detection and driver attribution across the Andes-Amazonia region, this study offers a proof-of-concept for a robust ecological monitoring system that has applicability well beyond the borders of Picachos, including the surrounding protected areas of Tinigua and Macarena.
Though Picachos is a strategic regional corridor between the Andes and Amazonia, there remains limited biophysical and sociocultural information to support effective decision-making processes.Detecting patterns of pasture intensification with remote sensing-derived spatio-temporal information is a critical tool in designing and evaluating land management policies that support desired, sustainable land use strategies [67].The results of this study record in great detail the short-and long-term effects of human influence on regional forest cover dynamics, facilitating data-driven dialogue with the community and developing land-use planning such as better management zonification in the foothills.
We detected, characterized, and attributed disturbances with relatively high agreement, but limitations of this approach may be summarized as four aspects.First, we were unable to identify multiple, sequential disturbances at a given pixel over the time-series, nor identify causal linkages between specific drivers, such as non-stand replacing disturbance preceding agricultural conversion [16,36] or pasture expansion [62,63].In future research, engaging local community knowledge could be relevant for establishing linkages between localized drivers and disturbances and identifying the most relevant metric for disturbance detection or driver attribution [68,69].Second, the training data used in building the RF classifier were not randomly sampled but rather sampled based on the availability of reference datasets (e.g., CLC, RapidEye, and high-resolution images).Despite limitations in our sampling approach, the out-of-bag sample provided by the RF model gave us an unbiased estimate of the classification error, which achieved a value of 13%.Third, the discrimination and attribution of drivers was more accurate for conversion to pasture and non-stand replacing disturbance than for conversion to agriculture, and distinguishing between non-stand replacing disturbance and conversion to agriculture remains a challenge.Fourth, although overall agreement was accurate for both disturbance and drivers' attribution detection, omission errors are significantly high, especially for disturbances (21%) and conversion to pasture (30%), which reflects how low image availability after cloud removal may affect area estimates.
Future work should thus include additional information on disturbance detection and driver attribution, for instance, the previous number of breaks before a given disturbance event [70], the incorporation of change duration [21] such as greening periods (positive trends) or browning periods (negative trends), and spatial context (pixel neighborhood), which has shown improvement in the mapping of changes [38].Finally, ensemble methods which account for the contributions of many different change detection algorithms have enhanced forest change mapping [71].

Conclusions
This study integrated spectral, trend, pattern, and topographic metrics derived from BFAST-Monitor and the Landsat archive to detect disturbances and identify related drivers in the Picachos foothills.We have been able to consistently differentiate the major drivers of forest disturbance within Picachos National Park.We confirmed that conversion pasture is the main cause of change, especially from 2007-2014, followed by non-stand replacing disturbance and conversion to agriculture.The results of this study should encourage the remote sensing community to further develop robust frameworks using dense time-series metrics to analyze the drivers and consequences of human-induced changes, especially in regions with low image availability or regular cloud cover, in support of sustainable resource use strategies that support people and forests alike.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/9/5/269/s1,Table S1: Confusion matrix at pixel-level for disturbances and stable forest, Table S2: Full area-weight confusion matrix for disturbances and stable forest, Table S3: Confusion matrix at object-level for drivers of disturbance using 40% data training from Random Forest, Table S4: Full area-weight confusion matrix for drivers of forest disturbance.
Author Contributions: P.J.M.-S., T.H. and J.V.D.H. conceived and designed the study.P.J.M.-S.performed the data processing.P.J.M.-S., J.V.D.H., and M.A.K. analyzed and interpreted the data, and wrote the paper.

Forests
have not been broadly tested in other regions complicated by fewer satellite observations, persistent cloud cover, or smaller disturbance patches.

Figure 1 .
Figure 1.Picachos National Park foothills study area, outlined in red.

Figure 1 .
Figure 1.Picachos National Park foothills study area, outlined in red.

Figure 3 .
Figure 3. Example of trajectories and disturbance detection using all available observations and associated drivers of change.This example illustrates the presence of abrupt and gradual changes: (a) Forest-to-agriculture conversion with a cyclical pattern after the breakpoint; (b) Non-stand replacing disturbance where NDMI trend exhibits a subtle decrease over time; and (c) Forest to pasture conversion where an abrupt breakpoint is followed by low, stable NDMI values.Vertical red lines indicate a breakpoint, blue lines indicate a seasonal + trend fitted model, and M denotes the magnitude of change.

Figure 2 .
Figure 2. Conceptual flowchart for disturbance detection and driver attribution.

Figure 3 .
Figure 3. Example of trajectories and disturbance detection using all available observations and associated drivers of change.This example illustrates the presence of abrupt and gradual changes: (a) Forest-to-agriculture conversion with a cyclical pattern after the breakpoint; (b) Non-stand replacing disturbance where NDMI trend exhibits a subtle decrease over time; and (c) Forest to pasture conversion where an abrupt breakpoint is followed by low, stable NDMI values.Vertical red lines indicate a breakpoint, blue lines indicate a seasonal + trend fitted model, and M denotes the magnitude of change.

Figure 3 .
Figure 3. Example of trajectories and disturbance detection using all available observations and associated drivers of change.This example illustrates the presence of abrupt and gradual changes: (a) Forest-to-agriculture conversion with a cyclical pattern after the breakpoint; (b) Non-stand replacing disturbance where NDMI trend exhibits a subtle decrease over time; and (c) Forest to pasture conversion where an abrupt breakpoint is followed by low, stable NDMI values.Vertical red lines indicate a breakpoint, blue lines indicate a seasonal + trend fitted model, and M denotes the magnitude of change.

Figure 4 .
Figure 4. Metric importance ranking for discriminating drivers based on the mean decrease accuracy of the RF classification.

Figure 4 .
Figure 4. Metric importance ranking for discriminating drivers based on the mean decrease accuracy of the RF classification.

Figure 5 .
Figure 5. Boxplots of the relationships between the eight most important metrics for classifying drivers of change showing median (black line) and interquartile range.

Figure 5 .
Figure 5. Boxplots of the relationships between the eight most important metrics for classifying drivers of change showing median (black line) and interquartile range.

Figure 6 .
Figure 6.Temporal distribution of drivers of change at patch-level.Year 2003 is not included due to the lack of Landsat imagery.

Figure 6 .
Figure 6.Temporal distribution of drivers of change at patch-level.Year 2003 is not included due to the lack of Landsat imagery.

Figure 7 .
Figure 7. Spatial distribution of drivers from 1999-2015 (upper).Disturbances are concentrated in close proximity to rivers.Temporal distribution of disturbance (below).Inset map represents disturbance locations of in Picachos' highest elevations (A) and Platanillo settlement (B).

Table 1 .
Summary of patch-level metrics used to characterize drivers of change.SRTM DEM = Shuttle Radar Topography Mission Digital Elevation Model (90 m).

Table 1 .
Summary of patch-level metrics used to characterize drivers of change.SRTM DEM = Shuttle Radar Topography Mission Digital Elevation Model (90 m).

Table 1 .
Summary of patch-level metrics used to characterize drivers of change.SRTM DEM = Shuttle Radar Topography Mission Digital Elevation Model (90 m).

Table 2 .
Area-weighted error matrix comparing BFAST-Monitor-and TimeSync-derived classification at validation sites.

Table 2 .
Area-weighted error matrix comparing BFAST-Monitor-and TimeSync-derived classification at validation sites.

Table 3 .
Area-weighted error matrix for disturbance driver attribution based on Random Forest output.

Table 3 .
Area-weighted error matrix for disturbance driver attribution based on Random Forest output.