remote sensing Land Cover Change in the Andes of Southern Ecuador—Patterns and Drivers

: In the megadiverse tropical mountain forest in the Andes of southern Ecuador, a global biodiversity hotspot, the use of fire to clear land for cattle ranching is leading to the invasion strongly fragmented due to the rapid expansion of the arable frontier and the even more rapid invasion by bracken. Unexpectedly, more bracken-infested areas were converted to pastures than vice versa , a practice that could alleviate pressure on forests if promoted. Road proximity was the most important spatial element determining forest loss, while for bracken the altitudinal range conditioned the degree of invasion in deforested areas. The annual deforestation rate changed notably between periods: ~1.5% from 1975 to 1987, ~0.8% from 1987 to 2000, and finally a very high rate of ~7.5% between 2000 and 2001. We explained these inconstant rates through some specific interrelated local and national political and socioeconomic drivers, namely land use policies, credit and tenure incentives, demography, and in particular, a severe national economic and bank crisis.


Introduction
Human-induced land cover change is a threat to biodiversity and ecosystem services [1], particularly in biodiversity hotspots of the tropics [2][3][4]. Here, land cover change dynamics are extremely diverse and intricate because they depend strongly on complex regional interactions between people and the specific environment [5], crosscutting local to global scales [6,7]. Therefore, the patterns and pace of local land cover change are difficult to deduce from broader scales or from other local studies, even in the same country [8].
In Ecuador, most of the studies have focused on the Amazonian lowlands [9][10][11][12][13][14][15][16][17][18][19][20][21][22] and coastal [5,8,23,24] regions, while land cover change in the Andes has been studied less, due, among other reasons, to the inaccessibility of the region. It is obvious that in remote areas like the Ecuadorian Andes, wide analysis of land cover change would not be possible without the use of remotely sensed data. However, the topographical complexity in such locations and particularly the high cloud frequency of the Andes makes the use of remote sensing techniques a challenge [25][26][27]. Accurate atmospheric and topographic corrections of the satellite scenes are a prerequisite to obtain satisfactory classification results [25,[28][29][30][31][32], but are not always performed due to time and resource limitations [33].
In the San Francisco Valley breaching the south-eastern Ecuadorian Andes and its surroundings (refer to Section 2.1), large areas of tropical mountain forest have been cleared, mostly for agricultural expansion of smallholder cattle-ranching [34][35][36][37][38][39]. Since 2002, this location has been the subject of two consecutive multidisciplinary projects that have been investigating the ecosystem functioning and services for sustainable land use management [38,40]. However, previous research only focused on particular subjects related to land cover and land cover change, with some of them relying only on analysis of one or a few satellite images covering short time periods. These few studies revealed a high degree of deforestation related to topography and road infrastructure [35,41,42]. Other socioeconomic studies from this area found a strong relationship between land cover change and private land tenure systems, cattle ranching, and decrease of local traditional ecological knowledge on plant use [34,36,43]. Moreover, a recurrent problem recognized in all the studies is the invasion of bracken fern (Pteridium arachnoideum and Pteridium caudatum), a fire-resistant and very competitive weed [38,[44][45][46] that colonizes unused burnt or inactive pasture areas, or even active pastures [44][45][46]. Bracken invasion frequently leads to pasture abandonment, which takes the area out of the production process, degrades ecosystem functions and services, and threatens the remaining natural forest [47]. Possible solutions for avoiding deforestation in the study area have been explored by means of payment schemes that include productive land use concepts [48].
Despite the knowledge hitherto gained, a comprehensive approach to study land cover change in the SE-Ecuadorian Andes over an extended period considering fragmentation, fire-pasture-bracken dynamics, and the study of proximate and local and non-local underlying drivers is still lacking.
Landsat data, with its spatial resolution of 30-60 m, longer term (since 1972), and free availability, is especially suitable for land cover change detection, as already demonstrated elsewhere in a variety of studies [49][50][51][52][53]. Due to the complex terrain and the difficult weather conditions in the study area, an accurate change detection analysis is only possible if an adequate atmospheric and topographic correction is applied to the satellite images. Additionally, other pre-processing steps like (i) image co-registration; (ii) cloud/cloud-shadow detection; and (iii) intercalibration are also necessary to produce accurate results. Hence, a framework to conduct the required processing steps with a minimum of manual work is needed to fully exploit this valuable source of information efficiently.
Consequently, the main aims of this paper are (i) the development and implementation of a semi-automated processing framework to allow an accurate classification of Landsat scenes in complex terrain; and, based on this, (ii) to present an integrated view of the land cover change process in the Upper Zamora and San Francisco Valleys from 1975 until 2001. The results of this study are fundamental for the ongoing research in this region, particularly to identify other potential alternatives of sustainable land management adapted to the geographical, ecological, cultural, and historical reality [54][55][56][57]. Moreover, the development of the pre-processing tool for the Landsat scenes is of benefit not only for future monitoring of the land cover change in this area but also for research on land cover change in tropical mountain forests in general.

Geographical Setting
The study area (~25 by 25 km) is located at the eastern slopes of the South Ecuadorian Andes in the Huancabamba depression, a region where the Andes are comparatively lower in elevation but with a higher topographical complexity ( Figure 1). This complexity and the strong altitudinal gradient from ~900 to ~3500 m above sea level (a.s.l.) result in a high heterogeneity of environmental conditions, among others, leading to an outstanding diversity of vascular plants [58] and other organisms, mostly harbored by the natural mountain rain forest [38,40]. The region has a tropical humid climate with an extremely wet period from April to July and a relatively dry period from September to December [59]. According to Homeier et al. [60], the forest upper boundary reaches an altitude of ~2700 m a.s.l., where it is replaced by subpáramo and páramo vegetation. Since 1982, the southwestern part of the study area has belonged to the Podocarpus National Park, and since 2000 the northern part has belonged to the Corazón de Oro Protected Forest. Politically, the area under investigation comprises the rural parishes (parroquias) of Imbana, Sabanilla, and Zamora in the Zamora-Chinchipe province, where most of the inhabitants are Mestizo settlers (~80%, [61]). The northern part of the study area has been mainly colonized by the Saraguros, a Quechua indigenous group from the highlands of Loja province. The main road crossing the area connects Loja (capital of the Loja province) and Zamora (capital of the Zamora-Chinchipe province).

Land Use
The cultural landscape in the study area is principally composed of a mosaic of forest patches, pastures, and bracken-infested areas (Figure 2), the latter being a less diverse and probably long-lasting system [39,44].
New pastures are established for three main reasons: (i) to counterbalance the fertility loss of old paddocks with low stock rates; (ii) to increase production [8]; and (iii) because of the invasion of bracken fern (Pteridium arachnoideum and Pteridium caudatum) [43][44][45].
As in similar areas of Latin America [62], slash and burn is the traditional way of establishing and rejuvenating pastures. After choosing an area of about one hectare, the inhabitants remove high value forest timber species [37]. During the beginning of the dry period, they cut parts of the understory and the remaining trees, let the area dry out, and burn the plot [37]. Normally, fire stops at the plot edge, as forest moisture inhibits its passage. However, strong winds, dry periods, the presence of neighboring dry bracken-infested areas [39], or/and the use of fire accelerators can lead to burning adjacent forest areas [37]. The clearing of more land than farmers can hold under use (due to labor scarcity or other household factors) leads to the establishment of bracken-dominated vegetation [37], mostly found in the steepest areas [34,36]. Once the soil cools down, farmers quickly plant pasture grasses, given that about three weeks after fire bracken fern sprouts vigorously [45]. If the emerging bracken sprouts are not eliminated, they will expand, making their eradication very time-consuming and labor-intensive [63].

Methods
As stressed in the introduction, atmospheric and topographic corrections and other pre-processing steps are a fundamental prerequisite to obtain satisfactory classification results for land cover change studies in regions of complex terrain. Because of the high amount of data and time and resource limitations, a semi-automated processing chain with a minimum of manual intervention was developed. In the following sections, we present a description of the pre-processing steps (3.1), the image classification (3.2), and the change detection analysis (3.3 and 3.4).

Pre-Processing of Landsat Scenes
We chose all available Landsat satellite images for the study area from the United States Geological Survey (USGS) website [64] which had less than 50% cloud cover and no data gaps originating from the failure of the system's scan line corrector (SLC) of Landsat ETM+. Our intention was to cover the longest attainable time period. This resulted in five images covering the period from 1975 to 2001 (Table 1). For spatial accuracy of change detection, the satellite images were co-registered to one reference image (2001) using the first-order polynomial function. Root mean square errors of less than half a pixel were achieved for all images.

Masking Clouds and Cloud-Shadows
Cloud and cloud-shadow masks were developed for the four scenes with cloud cover ( Table 1). The cloud and cloud-shadow mask of the scene from 1980 was digitized directly from a false color RGB composite due to the presence of only a few isolated clouds. For the other three cloud-contaminated scenes, we adapted the method of Martinuzzi et al. [27]. First clouds were identified by classifying band 1 (or 4 in the case of Landsat MSS) and 6 (thermal band only available in Landsat TM and ETM+) using digital number (DN) brightness ranges ( Table 2). Each classification was buffered to assure a complete inclusion of the clouded areas. In the case of Landsat TM and ETM+, the cloud mask resulted from the intersection of both classifications, since the cloud classification of band 1 included bare soil and some urban regions and the cloud classification of band 6 included dense forested areas. Due to the lack of the thermal band in the MSS scene, we visually corrected the cloud classification of band 4. Additionally, a buffer was applied to the cloud masks. Afterwards, cloud-shadow masks were created by shifting the cloud masks over the cloud-shadowed surfaces using the nearest neighbor algorithm. For this, cloud-shadows were identified in band 4 (or 7 for Landsat MSS) and distinctive points both in the cloud mask and cloud shadows were selected. Then, the cloud and cloud-shadow masks were combined and a further buffer was applied to remove remaining cloud-shadow-contaminated pixels originating from topographic distortions. The final cloud and cloud-shadow masks were manually inspected and corrected.

Atmospheric and Topographic Correction
To compensate for the varying atmospheric effects within the study area, we developed a fully automatic module (AtToCor) combining an extended version of the Second Simulation of the Satellite Signal in the Solar Spectrum radiative transfer code (6S code [65]) and a topographic correction module.
In order to apply the 6S code for area-wide and automatic operation, two modifications are implemented in AtToCor: (i) an area-wide adaptation to highly variable terrain altitudes; and (ii) an automatic scheme to consider an as-exact-as-possible representation of the atmosphere over the area of interest for the time of image acquisition.
(i) In the original version of the 6S code, the atmospheric absorption, reflection, and scattering coefficients are calculated for one single altitude. In the modified version, the parameters are calculated for 27 altitude levels spanning the entire range of the terrain in the respective satellite image (Figure 3). The 27 parameters are then used to estimate the effect of the altitude on each parameter by fitting second-order polynomial regression models between atmospheric absorption, reflectance, and scattering coefficients and altitude. The number of altitude levels was set to 27 to ensure a sufficiently high number of samples to train the polynomial regression models. This approach allows for the estimation of atmospheric effects in relation to the pixel's altitude from a co-located digital elevation model. (ii) The atmospheric pressure, water vapor, and temperature profiles in the original 6S version either rely on values of standard atmosphere or require data from radiosondes, which are mostly unavailable outside the industrialized countries, and not in mainland Ecuador. We implemented the option to extract and download the profiles of air temperature and specific humidity at different pressure levels for the over-flight day automatically from reanalysis data [66] at the National Center for Environmental Prediction/National Center of Atmospheric Research (NCEP/NCAR) website [67] ( Table 3). Using the center location of the image and the time of acquisition, the profiles are temporarily and spatially interpolated. The 6S code requires the proper definition of several other parameters for the conversion of top-of-atmosphere radiance values into at-surface reflectance values. We relied on the information of the sun-sensor-soil geometry (extracted automatically from the metadata of the scenes), on the atmospheric ozone content data and horizontal visibility, and on topographic variables. The total column ozone content value at the time of image acquisition is automatically extracted and downloaded to AtToCor from the Total Ozone Mapping Spectrometer (TOMS) satellite data [68]. Horizontal visibility was manually determined using data from a scatterometer installed in 2002 at a meteorological station in the study area [59]. For imagery acquired before 2002, we used monthly averages of the observation period 2002-2007. Topographic variables were derived from a digital elevation model (DEM) from the Shuttle Radar Topography Mission 3 (90 m resolution), which was bilinearly resampled to the resolution of the satellite imagery.
To test the importance of the factor terrain altitude for the atmospheric correction, we ran the original and the modified version of the 6S code with the bands of the Landsat ETM+ scene of 2001. We then compared the reflectance values of the bands corrected with these two modes.
To compensate for illumination effects arising from the rugged topography in our study area after the atmospheric correction, a topographic correction based on correlations between illumination and reflectance values was applied [69]. We used a semi-empirical approach, which first establishes a linear regression model between reflectance values and illumination. The corrected reflectance values are then calculated by subtracting the illumination effect estimated by the linear model from the uncorrected reflectance values derived by the 6S code.

Radiometric Intercalibration
To use the same spectral signatures for the selected land cover classes (see Table 4) in different years and for different sensors, the satellite scenes were radiometrically normalized (or intercalibrated). A straightforward regression method was applied [70], which uses the least squares of all pixels of one band to calculate normalization coefficients. It is based on the assumption that stable land cover types should have the same spectral signals in different images [71]. The accuracy of this method depends on the accuracy of the images' co-registration, as well as on the absence of significant reflectance changes between the scenes [72], such as changes in vegetation phenology [73]. Since the image co-registration performed well (i.e., the root mean square error was smaller than half a pixel) and changes in phenology do not play an important role in tropical evergreen mountain forests, we considered the regression method suitable for our case. We used the scene of 2001 as a reference for the correction of the TM and ETM+ imagery, and the scene of 1980 to correct the MSS scene of 1975. To evaluate the accuracy of this method we selected a subset of forest that was assumed to have not significantly changed during the observation period and we compared the deviation of the normalized reflectance values from the reference scenes before and after the intercalibration (Figure 4). In general, the intercalibration reduced the difference of reflectance between the scenes, improving the comparability among them.

Land Cover Classification
Prior to the classification, the area above 2700 m a.s.l. was masked out. This altitude is covered by subpáramo [60], which presents low human impact and a spectral signal similar to bracken-infested areas. Then, spectral signatures were extracted from training sites in the reference scenes of 2001 and 1980 for the land cover classes presented in Table 4. Due to the coarser spatial and spectral resolution of MSS scenes and the unavailability of ancillary data for image interpretation, only four classes were separated in the two oldest scenes: "forest", "pasture and bracken", "non-vegetated", and "water".
Training sites in the Landsat ETM+ image of 2001 were selected based on expert knowledge and visual interpretation of an ortho-photograph mosaic (1 m resolution), derived from several flights in 2001 and processed by aero-triangulation techniques (aerial photograph: [74]; training sites: [35]). Generally, a sample size 10 times higher than the number of bands is recommended to account for the spectral variability within each of the classes [75]; given the use of six bands, we collected between 70 and 409 training pixels for each class. Training sites for the scene of 1980 were collected only by satellite image interpretation because no auxiliary data taken at the time of image acquisition were available. False RGB composites were generated to enhance the contrast between different land covers and to facilitate visual interpretation. Additionally, we overlaid the RGB composite images of the other years to compare possible land cover changes as a further interpretation aid. Given that Landsat MSS has four bands, we collected between 43 and 341 training pixels. The spectral signatures of the land cover classes generated from the training sites of the 2001 scene were also used for classifying the 2000 and 1987 scenes. Likewise, the signatures extracted from the 1980 scene were also used for the classification of the 1975 scene.
The maximum likelihood classification (MLC) approach was applied because it was already successfully used for other Landsat classifications in the study area [35,41,42]. Moreover, for classification approaches in the Brazilian Amazon based on Landsat data, MLC has been found to outperform new methods such as the artificial neural networks and classification tree [76]. After classification, the class "water" was merged with the class "non-vegetated." The accuracy of the classification results was assessed by compiling a contingency matrix and calculating the kappa index of agreement [77,78]. Validation sites in the image of 2001 were selected based on expert knowledge and visual interpretation of the ortho-photograph mosaic (as mentioned above for the training sites collection). For all the other scenes (2000, 1987, 1980, and 1975), validation sites were selected based on visual interpretation of the RGB composites of each year. None of the validation sites were used to train the MLC.

Quantification of Land Cover Change and Habitat Fragmentation
We quantified land cover change by means of the post-classification inter-comparison technique, nowadays one of the most used change detection approaches [79]. It implies a pixel-by-pixel comparison of separately classified images where the same cell size is required. We resampled the two coarse resolution maps of 60 m to 30 m using the nearest neighbor algorithm to avoid information loss in the three fine spatial resolution maps. This procedure has been found to perform similarly to the resampling method, which coarsens the fine resolution map [79].
The proportion of transition between the periods 1975-1980, 1980-1987, 1987-2000, and 2000-2001 was calculated using the following formula: Pij = (Sij(t1)/Si(t0)) × 100 (1) where Pij is the proportion of transition from class "i" to class "j", Si(t0) is the surface of the "i" element of the transition matrix for the initial year, and Sij(t1) is the surface of the matrix element "ij" in the following time step. Annual rates of deforestation (r) were calculated as well, using the following formula [14,16]: where A1 and A2 are the forest areas at the beginning and end of a given time period, and t is the number of years in that time period. From an ecological perspective, the investigation of change rates and their reasons is not sufficient. Knowledge of habitat fragmentation in the natural system is fundamental for a better understanding of the ecological consequences, such as biodiversity loss, threat by e.g., invasive species, habitat decline, and climate change [80][81][82]. To detect habitat fragmentation in the different scenes, we used the program Fragstats [83]. We calculated the number of patches and the mean patch size of the main land cover classes (forest, pasture, and bracken), as both parameters are considered excellent indicators of landscape fragmentation [84,85].
The use of multiresolution imagery for land cover change detection can result in a bias. For each of the four abovementioned periods (1975-1980, 1980-1987, 1987-2000, and 2000-2001), this bias was estimated by consecutively reducing the spatial resolution to 60, 90, 120, and 240 m and to 90, 120, and 240 m for the imagery with the higher and the lower spatial resolution, respectively. The nearest neighbor resampling method was applied for resolution reduction. At each of the resolutions, the size of the deforested area was counted.

Analysis of Land Cover Change in Relation to Distance from Roads, Altitude, and Slope
Land cover change does not occur homogeneously in space. Certain site conditions or adjacency relationships-also called attractors of landscape change-result in some areas changing more rapidly than others [55].
Following the findings of Goerner et al. [41], Göttlicher et al. [35], and Thies et al. [42], we considered the distance from roads, altitude, and slope as the most important attractors of landscape change in our study area. Different studies have shown that the more easily reachable a forest is (close to roads and/or at lower altitude and/or slope), the higher the possibility of it being disturbed [86][87][88].
We quantified the annual land cover changes of forest, pasture, and bracken-infested areas in relation to these three attractors. Altitude and slope were obtained from the DEM and a roads vector file of the province of Zamora-Chinchipe was acquired from the Instituto Geográfico Militar del Ecuador (IGM) and adapted to historical state using the Landsat scenes.

AtToCor
RGB color composites of the Landsat ETM+ scene of 2001 before and after atmospheric and topographic corrections are presented in Figure 5. Top-of-atmosphere reflectance values exhibit a large proportion of haze (Figure 5a). After atmospheric correction the haze was completely removed both in the valleys (at low altitudes) as well as at the mountaintops (Figure 5b). After atmospheric and topographic correction, the differences in illumination and shadowing effects were clearly eliminated ( Figure 5c). However, the semi-empirical topographic correction produces an error on mountain crests. Given that this area is found in the subpáramo region, which was excluded from our classification, the abovementioned inaccuracy does not represent an error source in further processing steps. For future work where mountain crests are important, other topographic correction algorithms should be tested and implemented. The incorporation of 27 altitudes in the 6S code, to improve the atmospheric correction, had the clearest effect in the blue and NIR bands of the Landsat data (Figure 6a,d-e). The distinct outcomes for each band reflect the general pattern of atmospheric effects on electromagnetic radiation [89]. The results indicate that differences between the original 6S code (using only one fixed altitude) and the modified version were more pronounced at higher altitudes. Thus, we conclude that the data processed with the modified version allow for a better discrimination of vegetation types. For our purposes, this is of particular relevance for the differentiation between pasture lands and bracken-infested areas.

Accuracy Assessment of Land Cover Classification
High accuracy was obtained for all the classifications, with overall accuracies ranging from 94.5% to 98% and Kappa values from 0.75 to 0.98 (Table 5). These results suggest that the methodology followed to pre-process the Landsat images, including the atmospheric and topographic corrections and the further intercalibration, was successful.

Land Cover Change and Habitat Fragmentation
In 1975 ~28,500 ha (88.5% of the study area) were covered by forest. By 2001, ~8500 ha (~20% of the study area) were deforested and mainly converted into pastures or invaded by bracken fern (Table 6; Figure 7). The annual deforestation rate between 1980 and 2000 was 1%, which is slightly lower compared to the rate of 1.3% calculated by the Food and Agriculture Organization (FAO) [90] for the whole of Ecuador in the same period. In the first two periods (1975-1980 and 1980-1987) the annual deforestation rate was ~1.5%. In the third period (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) it slowed down to approximately half (0.8%), while in the fourth period (2000)(2001) an extremely high deforestation rate of 7.5% was observed. A big fire event (see Figure 7e, upper middle of the scene), which accounted for 3.5% of the deforestation rate in that year, contributed to this high rate.  Table 7). This practice of converting bracken-infested areas into pastures or reforestation areas, if strengthened, could be a measure to alleviate, at least in the short to medium term, the pressure on the natural forest. As described by Knoke et al. [48], reforestation of bracken-infested areas could reduce deforestation by up to 45% and, at the same time, may increase farmers' profits by up to 65%.
From 1975 to 2001, an accelerated forest fragmentation tendency was observed. Forest patches increased from 66 to 581 and the mean forest patch size decreased from 433 to 34 ha (Figure 8a). Forest patches are characterized by dry fire-prone edges that neighbor frequently burned pastures [91] or dry bracken fern-dominated vegetation [39]. Moreover, they are often degraded by logging, which increases the process of desiccation and fuel loading [92]. Hence, as fire is frequently used in our study area, fragmentation may further increase the ecosystem vulnerability to fires [93] by positive feedback loops, leading to further forest degradation and loss [18,92,94].  Table 7. Transition matrices of land cover change in the periods 1975-1980, 1980-1987, 1987-2000, and 2000-2001. 1980 Forest   Between 1987 and 2000, the number of bracken patches remained almost constant but the mean patch size increased from ~1.6 to ~2.6 ha (Figure 8b). Burning for pasture maintenance or expansion affected adjacent dry bracken-infested areas and the vulnerable borders of disturbed forest. Burnt forest edges were invaded by bracken and hence bracken patches expanded. This mechanism confirms the observations on forest-to-bracken cover change patterns observed in the field [37].
Between 2000 and 2001, the bracken spatial pattern was considerably affected by a big fire event. The number of bracken patches rapidly increased from 1153 to 1364 and the mean patch size decreased to 2.2 ha (Figure 8b). As shown in Figure 7c (in the upper middle sector), bracken started to rapidly spread in a patchy way after the fire event. However, at the time of the image acquisition in 2001 it had not yet colonized the complete burned area.
The use of two different spatial resolutions in the change detection analysis may result in a bias, because small-scale changes are only detected in the higher resolution data. By quantifying this deviation, we found a potential average bias of ~9% for the deforestation change trajectories between the Landsat MSS scenes (60 m of spatial resolution) and the Landsat TM and ETM+ scenes (30 m of spatial resolution), because only ~6% to ~12% of the deforested area was not detected if the 30 m spatial resolution data of the periods 1987-2000 and 2000-2001 was resampled to 60 m resolution ( Figure 9). This implies that the use of multiresolution maps in this study only introduces small errors since ~90% of the deforested patches were detected when the resolution was coarsened to 60 m. Even when the deforestation maps were further coarsened to 240 m, ~50% of deforested patches were still detected.  1975-1980, 1980-1987, 1987-2000, and 2000-2001 at different spatial resolutions.

Land Cover Change in Relation to Distance from Roads, Altitude, and Slope
As in other studies in montane forest areas [95] and in the Amazon lowlands [88,[96][97][98][99], we identified roads as the most significant attractor of landscape change and proximate driver of deforestation (Figures 10a and 11).
Between 1987 and 2001, new pastures were mostly created along the road Loja-Zamora, from Zamora city upstream until the confluence with the San Francisco River (Figure 10b) at altitudes between 1000 to 1500 m a.s.l. (Figure 11c). On the other hand, bracken-infested areas occurred mostly in the Upper Zamora River Basin and not always close to roads (Figure 10c). New bracken-infested areas mainly formed at higher altitudes than pastures (between 1800 and 2100 m a.s.l.; Figure 11d), possibly due to the higher competitiveness of bracken under colder and moister conditions [100]. Moreover, this altitudinal range coincides with the altitudinal range where both bracken fern species coexist (1500 and 2400 m a.s.l.) [101]. In a recent study, Laurance et al. [102] pinpointed the Andes of Ecuador as a "conflict area" for road construction, given its high agricultural economic potential, high biodiversity, and provision of valuable ecosystem services. Considering the history of land cover change in the southeastern Ecuadorian Andes, with roads as major attractors of deforestation, any future expansion of road infrastructure as foreseen by the National Strategic Mobility and Transport Plan (PEM) in Ecuador should be addressed with extreme care to avoid serious environmental impacts with irreversible consequences.
We would have expected to find more new bracken-infested areas on steeper slopes where pastures are less wanted [103]. However, new bracken-infested areas were found in flatter slopes compared to pastures (Figure 11c,d). The reason why pastures were found on steeper slopes can be explained by the expansion of old pastures into steeper neighboring lands, and by the depletion of flat areas near the roads. Bracken fern expands rapidly where burnt areas or poorly maintained pastures exist and where it has competitive advantages. Its lower invasion rate in steeper areas corresponds to the relatively open canopy of bracken's horizontal-oriented fronds, which facilitate the penetration of wind-dispersed seeds from competitor species and light [44,45,101].  (second row, b), pasture lands (third row, c), and bracken-infested areas (fourth row, d) with respect to the same attractors of land cover change. The grey tones of the bars indicate the total available area covered by forest (in hectares).

Potential Deforestation Drivers
The main potential drivers supposed to affect the different deforestation rates found in the studied periods are summarized in Figure 12. The link between potential drivers, actors, and deforestation will be discussed in the following sections.

Precipitation during the Dry Season
The Mestizo and Saraguro's land use system is highly dependent on weather. During the dry spells of the pre-humid Andes, farmers are used to burning land. Particularly, fires during drier periods favor fire dispersion and hence deforestation [93,105,106]. Therefore, we expected higher deforestation particularly if precipitation amounts during the dry season exhibit a clear negative anomaly.
Results, however, showed no relationship between the frequency of negative precipitation (or drought) anomalies (rainfall less than average rain −1 standard deviation) during the dry seasons and higher deforestation rates in the studied time frames (Figure 12c): The period 1975-1987 was characterized by a high deforestation rate of ~1.5%, but only one dry season presented a negative precipitation anomaly; for 1987-2000, with a deforestation rate of ~0.8%, we identified two negative precipitation anomalies. Finally, the extremely high deforestation rate observed in the period 2000-2001 does not coincide with any drought anomaly.

Rural Population Dynamics
From 1974 until 2001, the rural population in the study area increased from 1908 to 2427 inhabitants (Figure 12b). Annual growth rates varied in the observation periods: ~4.5% between 1974 and 1982, a decrease of −2.5% between 1982 and 1990, and a very small increase of ~1.5% between 1990 and 2001.
Consistently with many other studies in other regions [20,[107][108][109], we found a correspondence between population increment and the decrease of forest coverage in the first observation period (1975)(1976)(1977)(1978)(1979)(1980). Likewise, the lower rate of deforestation between 1987 and 2000 is probably a result of a decrease in rural population between 1982 and 1990 as well. However, the increase in population between 1990 and 2001 does not have a repercussion in forest loss dynamics. The latter suggests that other drivers or a combination of these (including population change) may have a stronger relevance in explaining forest cover change, as will be discussed in the next section.

Socioeconomic and Political Factors
Despite the national economic growth in Ecuador at the end of the 1960s and 1970s based on the discovery of oil in the north east and its export [110], a severe drought in Loja province in 1968 and the finalization of the Loja-Zamora road in 1962 resulted in the migration of thousands of smallholders to other regions of the country, including the study area [111].
Moreover, with the Agrarian Reform and Colonization Laws of 1964 and 1973, land clearing became the main condition for access to land titles, [9] and these to credits [112]. Land titles were issued at the regional office of the National Institute of Agrarian Reform and Colonization (IERAC) in Zamora city, which was opened in 1975. In accordance with the findings of Southgate et al. [9], Wunder [113], and Alvarez and Naughton-Treves [96] elsewhere, the high deforestation rate in the period 1975-1987 can be explained by the high number of newcomers conducting subsistence agriculture, their need for land, and the policies facilitating land titles and credits in non-forested areas (Figure 12a).
During the 1980s and 1990s, the Ecuadorian economy faced severe budget crises, high fluctuations of the export prices of oil, currency devaluations, earthquakes that heavily affected oil production, and a frontier war with Peru in 1995. The period 1997-1998 brought a very destructive El Niño event and the Asian economic crisis. In this period the government reduced credits offered to the farmers and the Agrarian Development Law in 1994 eliminated the condition of forest clearing for adjudications [43]. The IERAC office was closed and land titling processes were suspended from 1996 to 2002.
Parallel to these events, since the mid-1990s, as the Ministry of the Environment was created and the Law of Environment Administration was issued, the influence of environmental policy regarding protected areas became an issue in the country. In the study area, the declaration of the Podocarpus National Park (in 1982) allocated land to conservation purposes, despite lack of information among the local population. Formally, this declaration settled legal barriers to the conversion of forests to agricultural land and to the acquisition of property titles.
Between the 1970s and the end of the 1990s, there arose a rush for romerillo timber (Podocarpus oleifolius and Prumnopitys montana) in the study area, which brought a noticeably large gross income that was in part redeemed by the high costs of extraction [113]. Settlers were involved in the selective extraction of the few wild growing high value timber species and thus had less interest and time to invest in clearing forest for new pastureland. All these facts could explain the lower deforestation rate in the period 1987-2000 ( Figure 12a).
By the end of the 1990s, romerillo timber became scarce in the study area. Due to the extractive activity, the forest had already suffered disturbance through microclimatic alterations and an increase in the understory litter, resulting in a higher susceptibility to fire [114]. Moreover, since 1997 the authorities had begun to confiscate timber transported without permission [37] and extensive cattle farming turned into the unique relevant production activity.
In the same period, Ecuador's economy was in a tailspin of hyperinflation, which provoked the default of the foreign debt, the collapse of the Ecuadorian bank sector, high unemployment rates, and finally, the dollarization of the national currency at the beginning of 2000 [115]. The crises led to an extremely high migration to the USA and the EU, especially during the end of the 1990s [116,117]. Wunder and Sunderlin [118] expected lower deforestation in the period after the dollarization due to the loss of competitiveness, which harms the land-using sectors. Messina and Cochrane [18] relate the dollarization with land abandonment, forest transition, and a possible migration to other forested areas in the north Ecuadorian lowlands.
As in other regions of the country, Gerique [37] reports that several settlers migrated to Spain due to the economic crisis. Nonetheless, as previously found by Jokisch [119], Carr [120], Gray [121], and Gray and Bilsborrow [122], our results suggest that emigration did not lead to agricultural abandonment. The dollarization could have led people directly affected by the crisis and with more possibilities to move outside the country. However, poor people in rural areas, which did not have US money at their disposal, probably opted to increase subsistence agriculture or to expand their property as a security measure. In other tropical countries like Cameroon and Indonesia, it was found that economic crises motivated farmers to clear land [123] to, among other reasons, increase income security [124], particularly due to the rising prices of milk.
Furthermore, the profound mistrust of the banking sector after the crisis and the credit and financing difficulties in the agrarian sector could have provoked a relevant part of the remittances sent by the migrants, which were the second most important source of income for the country in 2000 [125], being invested in a way that resulted in the expansion of pastures. Cattle ranching is preferred to the production of crops because of many reasons: (i) the prices of milk and meat in regional markets are more stable than those of other products; (ii) the income is higher; (iii) it is less risky; and (iv) it awards a prestigious social status because it enables an accumulation of capital in regions with unsure loan and pension systems [37,62,126,127]. Stoorvogel et al. [128] mention that in this period the production of milk was preferred to other activities because it did not require expensive external inputs. Guerrero Cazar and Ospina Peralta [125] show that in 1999 even the farm gate prices had grown by 30.8% for cattle and 35.9% for milk. However, as in other areas of Ecuador [119], remittances were probably not dedicated to technological improvements of agricultural production, nor to diversification, but to the expansion of arable land.
The extremely high deforestation rate between 2000 and 2001 is then most probably an indirect effect of the wood scarcity and strengthened controls on wood extraction, and in particular a consequence of the economic crisis in the livelihood of the local population (Figure 12a). There were only two options for lumberjacks in particular and for local livelihoods in general, to abandon the area or to invest in cattle ranching as a mean of subsistence. In addition, even if it has been proven that protected areas [88], particularly in the Amazon, can act as effective barriers against forest fires and deforestation [129,130], the big fire (see Figure 7e) revealed by the satellite image classification within the perimeter of the Corazón de Oro Protected Forest (created in June 2000) shows that bad fire management and weak governance, especially during times of economic crisis, are important issues to be addressed with regard to forest protection.

Conclusions
The described semi-automated pre-processing framework for Landsat satellite data provides a solid basis for land cover change analysis in the tropical Andes of southeastern Ecuador. Although remote sensing applications in this area are challenging because the complex topography distorts the satellite signal and the high cloud frequency complicates the acquisition of images, we achieved good results in the correction of five Landsat scenes from 1975 to 2001, an essential step towards accurately classifying the land cover. To our knowledge, this is the longest period of satellite-based land cover change research in the region and it is the first time that an analysis of the fire-bracken-pasture dynamics has been included.
Between 1975 and 2001, in the Upper Zamora and San Francisco Valleys, strong forest loss and habitat fragmentation originated in the need for pasture expansion. However, non-sustainable fire practices accounted for almost one third of forest loss and degradation of usable land resources. The issue of burning more land than required supported the infestation of the land by bracken, whose spatial expansion was two times faster than pastureland between 1987 and 2001. The increase in forest patches and the decrease in their mean size, which made the remaining forest more vulnerable to fire, led to the observed rapid expansion of bracken fern surface. To our surprise, a considerable percentage of bracken-infested areas was also converted to pastureland, a practice that could alleviate pressure on forests if further supported by incentives.
The most important spatial element conditioning deforestation was the proximity to roads, given that pastures are easily reachable. Pasture expansion resulted in the use of neighboring areas at higher slopes. Bracken colonized all unused burnt areas and poorly managed pastures, with a higher incidence on flatter slopes. A higher increase in new pasture areas was found at lower elevations, while the altitudinal range where most of the new bracken-infested areas appeared coincided with the altitudes where both bracken fern species coexist.
Land cover change detection revealed that the deforestation rate in the study area was highly variable between 1975 and 2001. Precipitation during the dry season was not a relevant driver in any of the studied periods. The high deforestation from 1975 to 1987 is explained by the recently arrived settlers, who were motivated by land use policies favoring land titles and credit. The deforestation slowdown between 1987 and 2000 seems to be related to a decrease in the rural population, the cessation of the aforementioned land use policies, and more selective timber extraction activity. The creation of the Podocarpus National Park and the new environmental laws apparently helped forest protection as well. The extremely high deforestation between 2000 and 2001 is partially explained by wood scarcity and strengthening of controls over wood since the end of the 1990s, which made cattle ranching the only profitable productive activity. Still, we strongly believe that the economic crisis in Ecuador, which led to the dollarization in 2000 and to a national emigration wave to USA and EU, motivated farmers to expand their properties and clear more forest as a security measure, while emigrants contributed to this expansion through remittances. This is a subject that would merit more attention in future works.