Reconstructing Long Term High Andean Forest Dynamics Using Historical Aerial Imagery: A Case Study in Colombia

: High Andean forests are biodiversity hotspots that also play key roles in the provisioning of vital ecosystem services for neighboring cities. In past centuries, the hinterland of Andean fast-growing cities often experienced a dramatic decline in forested areas, but there are reports that forest cover has been recovering recently. We analyzed aerial imagery spanning the years 1940 to 2007 from nine administrative localities in the Eastern Andean Cordillera of Colombia in order to elucidate precise patterns of forest vegetation change. To this aim, we performed image object-based classiﬁcation by means of texture analysis and image segmentation. We then derived connectivity metrics to investigate whether forest cover trajectories showed an increase or decrease in fragmentation and landscape degradation. We observed a forest cover recovery in all the examined localities, except one. In general, forest recovery was accompanied by an increase in core habitat areas. The time scale of the positive trends identiﬁed partially coincides with the creation of protected areas in the region, which very likely furthered the recovery of forest patches. This study unveils the long-term dynamics of peri-urban high Andean forest cover, providing valuable information on historical vegetation changes in a highly dynamic landscape.


Introduction
Landscape change is the result of a broad variety of human activities with various effects on ecosystems, which depends on the nature and intensity of these activities [1,2]. Land cover change does not only interfere with the structure and functioning of ecosystems, but also affects the interactions with the atmosphere, ultimately contributing to local and global climate change [3,4]. In this context, habitat fragmentation modifies the availability of resources by creating artificial patch edges and by altering the degree of isolation of patches [5,6]. This, in turn, alters species distribution, abundance, and diversity [7,8], as well as carbon stocks [9], and ecosystems nutrient fluxes [10]. Particularly, landscape fragmentation poses severe threats to the native tree flora [11] and induces higher individual mortality [12] and extinction rates either at the species level [13] or at the population level [14,15]. Ultimately, the resulting biotic homogenization affects the resilience of forest remnants to current and future changes in environmental conditions [16,17]. the permanency of forest fragments, as well as the clear identification of regrown forest patches with corresponding age estimates (secondary succession) and their spatial configuration through time [44].
This study aims to investigate long-term changes in the extension and structure of high Andean forests by gathering detailed information on landscape history and patterns of landscape change through several decades. Specific objectives are: (i) to assess the forest cover change in the surroundings of Bogotá by means of aerial imagery covering a range of 60 years; and (ii) to evaluate the forest cover structure as described by spatial pattern indicators in relation to forest cover trends in the region.

Study Area
The study area encompassed nine localities in the hinterland of Bogotá, the capital city of Colombia ( Figure 1). Six of the study localities are situated in an arc north of the city in the department of Cundinamarca (Tabio, Soacha, San Francisco, El Rosal, Guasca, and Guatavita) and three are located within the capital district Bogotá D.C. (Torca, Pasquilla, and Sumapaz). These localities were selected because they match a network of permanent plots that was established as a basis for various botanical and ecological studies between 2013 and 2019 [39][40][41]44]. The high Andean forests close to Bogotá reach a modest height (15-25 m) and are typically dominated by trees and tall shrubs of the genera Weinmannia, Miconia, Cavendisha, Myrsine, Myrcianthes, Xylosma, and Daphnopsis [38,41,44,56]. Around the city, the forest remnants appear as scattered and fragmented patches. Some of these forested areas form part of a natural reserves network, created in the 1970s, and some have received their protection status only recently (as in the case of Los Encenillales in Pasquilla and El Encenillo in Guasca). The remaining forest patches are reported to be threatened by urban expansion, road construction [27,30], and by the opening of new areas to create pastures, especially in more rural localities [31,57,58]. Detailed information related to the study localities is presented in Supplementary Table S1. information on the permanency of forest fragments, as well as the clear identification of regrown forest patches with corresponding age estimates (secondary succession) and their spatial configuration through time [44]. This study aims to investigate long-term changes in the extension and structure of high Andean forests by gathering detailed information on landscape history and patterns of landscape change through several decades. Specific objectives are: (i) to assess the forest cover change in the surroundings of Bogotá by means of aerial imagery covering a range of 60 years; and (ii) to evaluate the forest cover structure as described by spatial pattern indicators in relation to forest cover trends in the region.

Study Area
The study area encompassed nine localities in the hinterland of Bogotá, the capital city of Colombia ( Figure 1). Six of the study localities are situated in an arc north of the city in the department of Cundinamarca (Tabio, Soacha, San Francisco, El Rosal, Guasca, and Guatavita) and three are located within the capital district Bogotá D.C. (Torca, Pasquilla, and Sumapaz). These localities were selected because they match a network of permanent plots that was established as a basis for various botanical and ecological studies between 2013 and 2019 [39][40][41]44]. The high Andean forests close to Bogotá reach a modest height (15-25 m) and are typically dominated by trees and tall shrubs of the genera Weinmannia, Miconia, Cavendisha, Myrsine, Myrcianthes, Xylosma, and Daphnopsis [38,41,44,56]. Around the city, the forest remnants appear as scattered and fragmented patches. Some of these forested areas form part of a natural reserves network, created in the 1970s, and some have received their protection status only recently (as in the case of Los Encenillales in Pasquilla and El Encenillo in Guasca). The remaining forest patches are reported to be threatened by urban expansion, road construction [27,30], and by the opening of new areas to create pastures, especially in more rural localities [31,57,58]. Detailed information related to the study localities is presented in Supplementary Table S1.

Acquisition and Processing of Aerial Photos
A total of 71 aerial greyscale photos of the selected localities ( Figure 1) was obtained from the Colombian Geographical Institute Agustín Codazzi (IGAC), covering a period from the 1940s up to the year 2007 (i.e., roughly one photo per decade and locality). The photos were scanned in high resolution and georeferenced in two consecutive steps using QGIS 2.18 [59]. To this aim, we used the thin plate spline as a transformation method for data interpolation, since it led to the highest precision among all the methods tested (first and second order polynomials). After a first georeferencing step, the overlapping areas were clipped to obtain a series of aerial scenes for each of the localities, which is equivalent in shape, size, and position ( Figure 2a). A subsequent, second georeferencing step was performed to achieve the highest precision possible. All photos were resampled at 1 m pixel size. From the initial 71 scenes, we selected 41 that covered the very same area in each of the study localities through the time series (Table 1).

Acquisition and Processing of Aerial Photos
A total of 71 aerial greyscale photos of the selected localities ( Figure 1) was obtained from the Colombian Geographical Institute Agustín Codazzi (IGAC), covering a period from the 1940s up to the year 2007 (i.e., roughly one photo per decade and locality). The photos were scanned in high resolution and georeferenced in two consecutive steps using QGIS 2.18 [59]. To this aim, we used the thin plate spline as a transformation method for data interpolation, since it led to the highest precision among all the methods tested (first and second order polynomials). After a first georeferencing step, the overlapping areas were clipped to obtain a series of aerial scenes for each of the localities, which is equivalent in shape, size, and position ( Figure 2a). A subsequent, second georeferencing step was performed to achieve the highest precision possible. All photos were resampled at 1 m pixel size. From the initial 71 scenes, we selected 41 that covered the very same area in each of the study localities through the time series (Table 1).   Traditionally, aerial photographs were classified through manual interpretation into land-cover/land-use categories. This required well trained specialists, and imposed numerous challenges [45]. In recent years, automated classification techniques have been used in interpreting aerial imagery [49,[60][61][62]. Object-based classification is used to efficiently delineate homogeneous polygons of land-cover/land-use based on panchromatic aerial imagery [63]. Object-based approaches rely on grouping neighboring pixels of similar properties to form objects, a process known as segmentation, before carrying out further image-processing analysis [64]. Resulting objects can be classified using properties such as tone and color, size, shape, texture, and contextual relationships [45]. Texture analysis was performed to enhance the pixel classification of the greyscale pictures using the glcm package [65] in RStudio 1.0.153 [66]. A moving window size of 7 × 7 pixels was adopted as an average value in order to avoid the drawbacks of either too large or too small window sizes [67,68]. Texture metrics were calculated for the four directions and the generated raster layers were examined. Based on the most informative metrics to visually distinguish our land cover classes, we retained three metrics for classification purposes: entropy, mean and second moment [69]. Segmentation for the object-based image analysis (OBIA) was carried out using the 'Seeded Region Growing Algorithm' [70] on a layer stack of the selected three texture metrics and the greyscale original aerial picture with the OBIA module in SAGA GIS 2.1.4 [71] to enhance the image classification at the pixel level by reducing the within-class variability [72,73]. We specified 5 pixels as the minimum object size in accordance with the smallest recognizable entity, i.e., the approximate average width of a tree canopy in the high Andean forests under study.

Supervised Classification and Accuracy Assessment
In order to perform the supervised classification, we manually selected our training polygons among the ones generated by the segmentation algorithm. We distinguished three classes: B (bare ground/built-up including main and secondary roads, houses, and bare fields), G (grass, cultivated fields, and low shrubs), and T (trees and tall shrubs). We did not distinguish forest plantations and natural forests because the quality of the aerial imagery did not allow for it. Supervised classification was carried out with all the selected aerial photos using the maximum likelihood criterion in SAGA GIS 2.1.4. [71] (Figure 2b) due to its efficacy in performing a classification when there are well separated land cover classes, as was the case in our study [74].
An accuracy assessment was carried out to derive the overall accuracy (OA), class-specific user accuracy (UA) and the producer accuracy (PA). We used a balanced, random sampling response design to obtain a more accurate estimate of the user's accuracy for the proportionally smallest class B and allocated an adequate sample size to it [75]. The three land cover classes (B, G, T) were defined as strata and 100 validation points were randomly generated for each class. When object-based segmentation is used as a means to improve classification results, it is recommendable to adopt a pixel-based response sampling design, due to its simplicity and lack of specific boundary or structural quality control [76]. Thus, our response design was built of 300 blocks [77] of 3 × 3 pixels around the randomly generated points. Within the blocks, we applied a majority criterion. When a precise block-based class attribution was not possible due to the low quality of the greyscale picture, an additional visual interpretation of the land cover classes based on the neighboring pixels was carried out. Good practice recommendations stipulate that the reference classification should be of higher quality than the map classification [78]. For this reason, the visual validation was carried out on the original higher-resolution greyscale aerial picture. Stratified area estimators as well as the overall and class-specific accuracy were obtained following the equations proposed by Oloffson et al. [79]. The error matrices resulting from this response sampling design, all accuracy values, and class area statistics are presented in the Supplementary Material (Tables S2 and S3).

Morphological Spatial Pattern Analysis (MSPA)
In order to get some insights into the quality of the underlined forest recovery, we complemented the land cover classification with the morphological spatial pattern analysis (MSPA). In principle, a classification of forest fragments would also need data on the floristic composition and structure of the vegetation over time. However, by deriving MSPA classes from geometric patterns within the forest cover, we can infer the information on the structural features of the forest cover trends. We are well aware that other fragmentation metrics and frameworks exist and are available [80,81]. However, we chose to derive a limited set of robust indicators, which provides a simple interpretation of features in the landscape and which relates to the needs of wildlife populations and plant communities in the study area, while maintaining considerable sensitivity to pattern changes over time [82].
The classified images of all localities were analyzed for geometry and connectivity using the MSPA QGIS plugin [83]. The MSPA uses a series of mathematical morphological operators to describe the geometrical properties of raster binary maps, e.g., forest (or foreground) = 1, non-forest (or background) = 0. The MSPA returns a raster layer with a minimum of seven mutually exclusive feature classes (core, islet, loop, bridge, perforation, edge, and branch) [83] that trace the geometry and connectivity of the binary layer components' spatial arrangement [84]. We performed the MSPA focusing on the forest class (T) and thus defined the remaining two classes (B and G) as background.
The definition of each MSPA category is provided in Table 2. According to their morphology, these categories can be subdivided into four main groups: (i) core, (ii) islet, (iii) connectors and branches, and (iv) edge [80]. Table 2. Definition of the class groups and categories for the morphological spatial pattern analysis (MSPA) adapted from [85].

MSPA Class Group MSPA Class Category Definition
Core Core Internal forest pixels whose distance to the background is higher than a specified threshold value (here: 5 m) We produced forest/non-forest binary layers from the land cover imagery, removed the raster polygons smaller than the selected threshold size (50 pixels), and replaced them with the pixel value of the largest neighboring polygon (Figure 2c). Subsequently, we performed the MSPA by using the following parameterization: Foreground Connectivity = 4 pixels, Edge Width (s) = 5, Transition = 1, and Intext = 0. Foreground Connectivity was set to four pixels to constrain to cardinal directions only [84], but also to take into account the maps' high spatial resolution and to avert an overestimation of the connectivity in the connector's classes. Edge width was set to five meters, which we considered to be an adequate threshold in our particular study system, to distinguish between the fragment edge and core areas. Transition was set to one to illustrate all the detected connections [84]. The Intext value was set to zero, as opposed to the default Intext = 1, to avoid the subdivision of the internal background into the two categories core-opening and border-opening [84] (Figure 2d). The processing workflow is exemplified in Figure 2 and all the processing steps of the aerial imagery are summarized in the Supplementary Figure S4. Finally, an analysis of the correlation was performed to investigate the relationships between the class T area and MSPA groups defined in Table 2.

Results
The supervised land cover classification performed satisfactorily, returning good OA values (see Supplementary Tables S2 and S3 for the complete error matrices and accuracy indices of all the classes and localities through the time series). Regarding the class-specific accuracy values, the class T and G returned good values of UA and PA. Finally, the class B returned satisfactory values for UA but had low PA values in few cases indicating omission errors.
In the period between 1940 and 2007, we found an overall increasing proportion of forested surface in eight of the nine study localities (Figure 3 and Supplementary Table S5). The only exception with a net loss of forest cover was found in Pasquilla (−6.77 ha), while the highest increase was observed in Torca (+289.40 ha, Figures 3 and 4a-e). A consistent increase in forest cover was detected in Guasca and Guatavita (+141.19 ha and 115.61 ha, respectively), and also in Soacha (+106.88 ha) and Tabio (+99.64 ha; Figure 3 and Supplementary Table S5).
Based on the forest cover trends identified, the study localities can be divided into three groups: (i) Group 1 includes Guatavita, Soacha, Guasca, and Tabio. This group is characterized by a well defined increase in forest cover with no or very little decrease during the study period; (ii) Group 2 includes Torca, San Francisco, and Sumapaz, and is characterized by an initial increase in forest cover, followed by a decrease occurring around the 1990s and then, again, a final increase; (iii) Group 3 includes Pasquilla and El Rosal, which are characterized by an initial decrease of forest cover and a later increase that started in the 80ies.
The results of the Morphological Spatial Pattern Analysis (MSPA) are presented in Figure 5 and in the Supplementary Table S6. Overall, due to the peri-urban nature of the sampled area, the class background (i.e., non-forested areas corresponding to classes B and G) was the most abundant one in terms of percentage through time, except for El Rosal in 2007, Guatavita in 2007, and Torca from 1962 until 2004.
The share of the connector group (loops, bridges, and branches), which includes the potential corridors for the species dispersal that connects the core areas [86], decreased noticeably through time only in Guasca and Pasquilla. The share of the edge group (perforations and edge MSPA classes), which is directly related to fragmentation [5,86,87], increased constantly through time in Guatavita and Tabio, and also slightly in Sumapaz. (1) well defined increase in forest cover, with no or very little decrease during the study period; (2) initial increase in forest cover, followed by a decrease around the 1990s and then, again, a final increase; (3) initial decrease in forest cover and then, an increase around 1980. B = bare ground/built up (including main and secondary roads, houses, and bare fields); G = grass/cultivated fields (including areas with low shrubs), and T = trees and tall shrubs (forest cover). (1) well defined increase in forest cover, with no or very little decrease during the study period; (2) initial increase in forest cover, followed by a decrease around the 1990s and then, again, a final increase; (3) initial decrease in forest cover and then, an increase around 1980. B = bare ground/built up (including main and secondary roads, houses, and bare fields); G = grass/cultivated fields (including areas with low shrubs), and T = trees and tall shrubs (forest cover). Based on the forest cover trends identified, the study localities can be divided into three groups: (i) Group 1 includes Guatavita, Soacha, Guasca, and Tabio. This group is characterized by a well defined increase in forest cover with no or very little decrease during the study period; (ii) Group 2 includes Torca, San Francisco, and Sumapaz, and is characterized by an initial increase in forest cover, followed by a decrease occurring around the 1990s and then, again, a final increase; (iii) Group 3 includes Pasquilla and El Rosal, which are characterized by an initial decrease of forest cover and a later increase that started in the 80ies.
The results of the Morphological Spatial Pattern Analysis (MSPA) are presented in Figure 5 and in the Supplementary Table S6. Overall, due to the peri-urban nature of the sampled area, the class background (i.e., non-forested areas corresponding to classes B and G) was the most abundant one in terms of  Considering all the localities, the correlation analysis carried out between the area of class T (trees or high shrubs) and MSPA groups, i.e., core, islet, connectors and edge (see Supplementary Table S7 for correlation coefficients and relative p-values), showed a significant positive correlation between the forest cover and the group core (r = 0.97, p < 0.001) and a slightly positive correlation with the group edge (r = 0.59, p = 0.001), which was often was stronger at an individual locality level (e.g., El Rosal: r = 0.987, p = 0.02; San Francisco: r = 0.913, p = 0.031; Tabio r = 0.969, p = 0.031). Only two localities that showed a positive correlation between the forest cover and core (Guasca: r = 0.868, p = 0.057; Torca: r = 0.987, p = 0.002), also showed a negative correlation between the forest cover and the group connectors (Guasca: r = -0.858, p = 0.063; Torca: r = −0.904, p = 0.035). Interestingly, the localities belonging to the first and second groups defined on general forest cover trends, had overall positive values of correlation with the group core (e.g., Guatavita: r = 0.953, p = 0.012; Torca r = −0.987, p = 0.002; Sumapaz: r = 0.95, p = 0.046) and negative values with the islet (e.g., Guatavita: r = −0.924, p = 0.025; Torca: r = −0.983, p = 0.003). The localities belonging to the third group did not show relevant correlation with the core, islet or connectors, but showed a positive correlation tendency with the group edge (e.g., El Rosal).  Table 2).  Table 2).

Spatial Resolution Potential of Historical Aerial Imagery and Time Coverage
Our approach carries the advantage of analyzing a longer time series than normally allowed by mere satellite data at a very high spatial resolution. It allows tracking forest cover changes at a forest fragment scale, which becomes particularly useful for discussing the actual forest species composition and structural parameters, when related to vegetation assessment in the same sites.
On the other hand, the land cover classification generated by aerial images lacks multispectral information and thus limits classification performance. We found that this is particularly true for the accuracy of the rare class B. Even though OA, UA and PA values for classes T, G and B generally fulfilled a minimum threshold of goodness, some omission errors are noticeable for the PA values for class B in some of the aerial scenes. User's and producer's accuracy are sensitive to a number of factors, including the number of random validation points generated and the sampling method employed [88]. However, in our case the low values are most likely due to the relatively small proportion of class B (cf. Figure 2b). In some cases, the area's proportion of the mapped class (W), as incorporated in equations 1-5 of the accuracy calculations found in Oloffson et al. [79], was several times smaller in comparison to the other classes, which led to overall lower accuracy values for class B. Additionally, the lack of ground-truth data for past aerial scenes may have also influenced the goodness of selected training polygons. However, classes B and G were not in the focus of the present study.
Due to the unavailability of more recent, aerial photographs for our study area, our analysis stops in 2007. Our limited timeframe does not allow us to assume that general land cover trends have not changed in the more recent past (i.e., after 2007), considering that there are several urbanization and infrastructure development projects currently being planned or already implemented that may interfere with the recovery of the forests under study. For instance, there is evidence that forested areas close to El Rosal have been substantially decreasing between 2000 and 2015 [89]. Today, unfortunately, only a very little proportion of the study area in El Rosal (and contiguous Subachoque) is under some regime of protection, i.e., the nature reserve Hacienda la Laja, Subachoque, created in 2015 [90]. To the contrary, in Guasca, the suspension of mining activities in the 1990s culminating with the creation of the protected area "El Encenillo" in 2007, leading to a continuous recovery of the forest cover (Figure 3). In a few cases, the limitations to the interpretation of forest cover trends are caused by the incomplete time series of the aerial images, as in the case of Soacha, for which only three aerial scenes were available.

Individual History of Localities and Forest Cover Trends
The observed locality-specific variation in the forest cover trends appears to be connected to the individual history of each of the localities. In those rural localities that did not experience the worst effects of the country's five decade-long armed conflict between the government, the paramilitary groups and the guerrilla [91,92], and which possess a rather 'smooth' topography, at least within the analyzed area, as Guasca, Tabio, Soacha, San Francisco, El Rosal, and Guatavita, the general forest cover recovery identified here probably is related to a mix of abandonment of rural practices by those seeking better economic opportunities in urban conglomerates [22] and an increase in low-density settlements of secondary homes of people living in Bogotá, who allow forest regrowth within their properties [27], jointly with the widespread implementation of protection efforts since the 1970s [24]. Within this general outline, Torca is exceptional since it is contiguous to the main urbanization of Bogotá and hosts various highly diverse and partially well conserved forest fragments in its steeper parts [41,44]. During the past century, Torca experienced a considerable decrease in agriculture and animal husbandry in favor of an urban expansion, which attracted middle-and upper-class residents that do not engage in rural activities [93]. The latter, together with both the implementation of the protected area "Bosque Oriental de Bogotá" in the 1970s [94] and the harsh topography of a large portion of this area, played a key role in regulating land use conversion. The major part of forest recovery took place in areas of low elevation close to the main motorways, which used to be cleared regularly. However, recent plans to expand the motorways [95] and to build up a massive urbanization project [96] pose new threats to recovered forest fragments, especially at low elevations, where the topography is less pronounced.
The case of Pasquilla, as the only locality in our study with an overall net forest cover loss, can be illustrative of a well populated agricultural landscape in the close vicinity of Bogotá. During the 1940s, Pasquilla became heavily populated and agriculturally exploited, due to the increasing violence in the country and the migration towards the cities [57]. This is consistent with our observations of an initial decrease in forest cover in Pasquilla, right from the 1940s, followed by a phase of slow recovery from the 1980s. Between 1995 and 2009, the main land use shifted from agriculture to animal husbandry, especially affecting the forest cover and its ecological integrity [58]. This is visible in our results by the increased area of size class G (Figure 3). The land use shifts led to an expansion of human-induced azonal páramos and the opening of forest areas for cattle grazing, which affect both the margins and the interior of forest patches, leading to the perforation and size reduction of the fragments, altering the forest structure and its regenerative potential [31,97]. Finally, the considerable recovery of forest core areas in 2006 may be directly linked to the increased protection efforts, together with an emerging environmental commitment of the local population, culminating with the creation of the protected area "Los Encenillales" in 2004 [58]. The impact of the 1940s migrations on forest cover changes is also evident in Sumapaz, which is far from the city of Bogotá. In the mountainous region of Sumapaz, and the eastern portion of the Tolima department, the violence and the armed conflict started already in 1920 and caused the displacement of the population [98]. Consistently, we found an initial increase in forested areas during the 1950s and 1960s in the Sumapaz region, which may result from the rural depopulation that brought farmers from remote areas such as Sumapaz to Pasquilla and other localities considered safer. It may be noted that our study locality in Sumapaz is outside the limits of the National Park (PNN Sumapaz) and thus is not under any protection regime.

Insights from MSPA on Forest Cover Trends
Besides the core areas, a particularly relevant MSPA class for the dynamic of secondary forest fragments is islets. These correspond to patches of forest vegetation that are potentially vulnerable to disappearance due to their shape and size (i.e., they are small and/or elongated, thin, and isolated). Islets may also act as nuclei for forest growth by serving as stepping stones for pollination and plant dispersal [99,100]. A decrease in islets alone can be either interpreted as a sign of a re-gained connectivity between fragments and the inclusion of former stepping stones into a more extended forest fragment network or as a proof of a more stringent delimitation of cultivated areas, and therefore a more extensive exploitation of the land by the cutting of standing trees within pastures. In our sample, all the localities that recorded a decrease in the class islet for the last period of the time series were also characterized by an increase in the class core. This indicates that the observed islet decrease is mostly due to their inclusion into the forest fragment network and a general restoration of connectivity. Overall, we found a general positive correlation tendency between the forest cover and the core group, and a negative correlation tendency between the forest cover and islets. This was particularly significant for Soacha and Guatavita, where the islet decrease was constant through time, and for Torca, where the decrease trend was only interrupted in the year 1993, as underlined by the positive correlation of forested area with the group core (Guatavita: r = 0.953, p = 0.012; Torca: r = 0.987, p = 0.002) and its negative correlation with the group islet (Guatavita: r = −0.924, p = 0.025; Torca: r = −0.983, p = 0.003).
In the localities with increased islet values for the last period (Sumapaz, San Francisco, and Tabio) the interpretation is more complex: in Sumapaz, the final increase in islets was accompanied by a total gain in the forest class T (+18% in 1996, Supplementary Tables S5 and S6), a slight increase in core (+6%), a stationary edge (+0%), and a slight decrease in the background (−2%) and connectors (−5%). This indicates that new forest nuclei have formed, as there is no major change in any of the other classes. In San Francisco, the increase in islets was particularly high and it was accompanied by a total gain of forest (+27% in 2007), a slight decrease in background (−13%), a little increase in core (+4%), and a marked increase in connectors (+113%) and edges (+43%). This indicates that new tree nuclei have formed and then connected to other fragments, however, the high edge value implies that cores might be broken or changed in shape and so an increased fragmentation of the main forested area is assumed here and is supported by the positive correlation between the forested area and the edge group (r = 0.913, p = 0.031). In Tabio, the increase in islet is linked to an overall forest gain (+17% in 2007, Supplementary Table S6), a decrease in core (−6%) and background (−9%), and a marked increase in connectors (+64%) and edges (+4%). This indicates the formation of new tree nuclei in previous background areas, but also a certain degree of fragmentation at the borders of the major forest fragments. This interpretation is confirmed by an increase in edges and connectors, and a decrease in core and background, which is also underlined by the strong positive correlation of forested area and the edge group (r = 0.969, p = 0.031). The mountains of Tabio are included together with other portions of our study area in the special management area "Reserva Forestal Protectora Productora de la Cuenca Alta del Río Bogotá" [101] since the 1970s. This light protection status still allows for agricultural activities and infrastructure development, inevitably resulting in forest clearing to create pastures, and thus the fragmentation of the previously continuous forest cover. However, except for Guatavita and Tabio, the other localities covered by this special management area are also under other protection figures, and the limits of the Reserva overlap with protected areas with stricter regulations. For Guatavita, the establishment of the Reserva seems to have been beneficial in promoting forest cover and islet inclusion in the main core areas, even though the positive correlation between the forest cover and the edge group for this locality may indicate fragmentation. Conversely, our findings for Tabio underscore that a mild protection regime, extensive agricultural exploitation, and infrastructure development together can hamper the protection and regeneration of forests considerably. Thus, the hazards associated with urbanization and road construction, already reported as key drivers of deforestation in the area [27,30,31,57,58], should be taken into consideration when formulating forest management strategies.

Possible Factors Influencing Forest Recovery
Our results show a significant increase in forest cover and a general recovery of forest connectors in the examined localities between the 1940s and the early 2000s. Forest regeneration on abandoned agricultural lands, by enhancing landscape connectivity, can potentially contain the ecological damage of deforestation and land degradation [53,102,103]. Restored connectivity, together with the increase in size forest fragments is essential to guarantee the freedom of movement and availability of resources to local fauna populations [104][105][106][107][108].
The overall forest recovery, first described by Etter [47] for two localities in the Bogotá high plain for the period 1940-1996, was recently confirmed for the region north of Bogotá for the years 1985-2015 [27], and for the Colombian Andes in general, for the years 2001-2010 [52]. Our study, which incorporates localities in the southern part of the Bogotá high plain and evaluates a large time scale, further supports and extends these regional evidences.
However, the identified positive trend diverges from the high deforestation rates outlined for Andean lowlands [24] or the Amazon region [109,110]. It also differs from another locality in the Colombian high Andean department of Quindío, where a reduction in forest cover was reported for the period 1954-2009 [48]. The trend of forest recovery in the hinterland of Bogotá also differs from the results from other high Andean regions, such as southern Ecuador [111,112] and southeastern Bolivia [113], which experienced continuous forest loss. Nevertheless, forest recovery was also identified in other similar high Andean systems, such as in the San Martín department, in northern Peru [114].
As observed by Hecht [21], the patterns and drivers of deforestation are likely to vary between localities due to their individual history and accessibility. Dramatic deforestation takes place mainly in economic frontier zones [21] such as the Andean lowlands, the Amazon [24,109,110] or Southern Ecuador [111,112], where population growth, the expansion of road networks, and the lack of protection favors timber extraction [112]. On the other hand, the hinterland of heavily populated cities, such as our study region, is a complex system, where the observed forest regrowth is mainly driven by rural-urban migration and population decline, eventually leading to the abandonment of rural activities [22,27,47,53]. These factors considerably contribute to forest regrowth together with the inaccessibility of steep and remote forested areas [27,47].
Even though forest cover dynamics are likely to be higher in peri-urban areas than in more remote forest fragments [27], the positive impacts of protected areas are evident in our study area, despite the proximity of the city. Conservation efforts through the establishment of protected areas appear to be effective at limiting deforestation and land-use change, as previously reported for Colombian montane forests [24,47], even without the rigorous and efficient implementation of conservation strategies in place and despite the transformation of areas inside national parks to some degree through anthropogenic activities [30].

Conclusions
The unique time coverage of the here employed, generally under-utilized, aerial imagery allowed us to provide a detailed reconstruction of vegetation dynamics at the landscape level, which in several cases clearly mirrors the historical changes that occurred during the occupation and exploitation of these rural areas. Forest recovery alone is expected to improve the systems' ecological conditions, but both the land cover types prior to a transition, and the land cover established by the transition, impact the ecosystem structure and functioning (see [115,116]). Accordingly, in the continuously intervened systems this study is focused on, it is highly likely that the legacy of centuries of land use is still relevant for todays' ecosystem remnants [117,118]. Similarly, when it comes to forest tree diversity and resilience, the effects from interventions in past centuries are likely to be relevant and noticeable even in a not-so-close future [119,120]. For this reason, generating locality-scale, longer-term information on land cover changes is key, considering that the majority of the studies focuses on broader geographical scales and employs satellite data with a shorter time-coverage. Moreover, elucidating not only the extent but also the nature of forest recovery as to its spatial configuration provides deep insights into forest fragment-related cover changes, and helps to understand the drivers and the processes shaping plant communities and the associated provision of ecosystem services [121].
Our findings of forest-recovery in the peri-urban region around Bogotá are encouraging. Nevertheless, they only apply to a limited sampling area. The aerial pictures that we used are available for all Colombia and are currently being georeferenced by the IGAC. Thus, an upscaling of the approach described in this study is possible and would be important to model the developments in adjacent administrative departments and high Andean forests ecosystems in general. This is very relevant, considering that the understanding of regional patterns of land use and land cover conversion in heavily populated high Andean systems is necessary to formulate sound land management plans that can potentially prevent broader scale, irreversible landscape degradation [113].
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/8/788/s1, Table S1: "Detailed information on study localities", Table S2: "Accuracy indices values for each class in each locality through time series", Table S3: "Error matrices for each locality through time series", Figure S4: "Processing steps of the aerial imagery", Table S5: "Forest cover extent and its changes in the study localities", Table S6: "MSPA classes groups percent change through time series in each locality", Table S7: "Correlation coefficients of forest cover and MSPA groups".