Remote Sensing Based Spatial Statistics to Document Tropical Rainforest Transition Pathways

In this paper, grid cell based spatial statistics were used to quantify the drivers of land-cover and land-use change (LCLUC) and habitat degradation in a tropical rainforest in Madagascar. First, a spectral database of various land-cover and land-use information was compiled using multi-year field campaign data and photointerpretation of satellite images. Next, residential areas were extracted from IKONOS-2 and GeoEye-1 images using object oriented feature extraction (OBIA). Then, Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) data were used to generate land-cover and land-use maps from 1990 to 2011, and LCLUC maps were developed with decadal intervals and converted to 100 m vector grid cells. Finally, the causal associations between LCLUC were quantified using ordinary least square regression analysis and Moran’s I, and a forest disturbance index derived from the time series Landsat data were used to further confirm LCLUC drivers. The results showed that (1) local spatial statistical approaches were most effective at quantifying the OPEN ACCESS Remote Sens. 2015, 7 6258 drivers of LCLUC, and (2) the combined threats of habitat degradation in and around the reserve and increasing encroachment of invasive plant species lead to the expansion of shrubland and mixed forest within the former primary forest, which was echoed by the forest disturbance index derived from the Landsat data.


Introduction
Madagascar is well known for the richness of its flora and fauna and is considered one of the top biodiversity hotspots worldwide [1,2].In combination with its long isolation, Madagascar's continent-like diverse landscapes and regionally varied climate has given rise to its extraordinarily high levels of biodiversity and endemism [3,4].With about 80% of its biota found nowhere else on Earth [5,6], conserving the country's biodiversity is a high national and international priority [4,7].However, Madagascar is also seeing the unprecedented declines in animal and plant populations that are occurring throughout the world [8].Land-cover and land use change (LCLUC) is the leading driver of biodiversity loss [7,8] and can also significantly alter local, regional, and global climates [9,10].Because tropical forests harbor the highest concentration of the world's terrestrial biodiversity and are critical to sustaining a range of essential ecosystem services, they are a focus of many conservation efforts [11].
Madagascar has experienced widespread deforestation [12][13][14][15], which has resulted in the loss of approximately 90% of its original forest cover [14].A growing body of research indicates that most of the loss is from anthropogenic activities, including deforestation for human settlements and agriculture, overexploitation of timber for cooking fuel and construction, logging operations for precious hardwoods and clearance for mining [7].The impact of deforestation extends well beyond habitat loss through the unleashing of secondary factors that, alone or in combination, further erode the viability of already vulnerable populations.For example, shrinking and degraded habitats increase opportunities for hunters, pathogens and invasive species to prey upon, infect or out-compete resident populations while species unable to disperse between fragmented habitats can suffer from inbreeding depression [16][17][18][19].Together, these factors may have dramatically outpaced the effects of climate change.
Madagascar's original network of Réserve Naturelles Intègrales was established in 1927 although their boundaries were not set by Decree 66-242 unit 1966 [20].At 2,228 ha, Réserve Naturelle Intègrale de Betampona (BNR) is an isolated remnant of lowland rainforest that once stretched across the entire east coast and is representative of the type of pressures that threaten Malagasy forests.Madagascar's population, growing at 2.8% [21], remains steeped in poverty which is most deeply felt in the rural areas [7].Most families living around BNR are subsistence farmers who practice slash and burn agriculture (tavy), hunt wildlife for consumption and depend on wood they collect for cooking fuel and construction [18].The tavy system begins with burning a patch of forest, using the land for multiple years to grow rice and other crops after which the land is left to lay fallow for several years.This process is repeated, but as the time needed to restore soil fertility increases with each cycle, the duration of the fallow period is correspondingly lengthened [22].However, the traditional tavy system breaks down when the land available is inadequate to meet the needs of growing families and, as a result, many farmers have been forced to decrease fallow times and accept lowered yields [22].
Tavy fires escaped into BNR after its formation and account for some sections of existing secondary forest that were colonized by invasive exotics and the pioneer native species Madagascar cardamom (Aframomum angustifolium) [20].Infractions continued through the 1980's, including the conversion of the purported 200 meter "zone of protection" into tavy fields, due, in part, to insufficient numbers of trained staff overseeing the reserve [14,20,23].In the early 1990s, Madagascar National Parks developed a collaboration with the conservation non-governmental organization (NGO), Madagascar Fauna and Flora Group (MFG), to serve as their research partner.The presence of a research team in the forest and the development of positive relationships and programs with community members have greatly reduced infractions, especially when compared to other sites [18].Although the MFG's programs have lessened the likelihood of fires entering BNR, the environmental repercussions of previous fires now represent the greatest threat to its biodiversity.At least 15 invasive plants are established in the reserve, and some of them have the potential to irreversibly alter forest structure and cause the extirpation of countless endemic animal and plant species [24].Other existing pressures can serve to hasten those losses.For example, as the density of humans living near the Reserve increases, so does the density of their domestic animals and associated pest species.The negative impact of highly invasive exotics such as rats (Rattus rattus and R. norvegicus), domestic cats (Felis catus) and dogs (Canis lupus familiaris) on native species has been documented [25,26] and is the focus of current research in BNR.Invasive exotic and domestic animals, whether due to predation, competition or disease transfer, can interact with other factors to reduce the population size of a species to a point from which it cannot recover.The challenge of developing strategies to conserve Madagascar's increasingly threatened biodiversity while also improving the health and welfare of the Malagasy people is one shared throughout the island.Interdisciplinary teams are essential to meeting this complex challenge.
Remote sensing is a powerful tool that enables monitoring temporal dynamics of forest changes over large areas [27], as well as the drivers of biodiversity loss in tropical rainforests [6].Multiple studies report successful applications of remote sensing in monitoring forest degradation and loss in Madagascar [14,28,29], however, few studies focus on quantitative estimation of forest degradation, transition pathways, and their anthropogenic drivers.Efforts should be made to capture not only information on deforestation but also drivers of forest change, pathways, and where and to what extent forest degradation occurs to aid in developing targeted conservation actions, monitoring and evaluating executed conservation management plans or, for example, exploring the feasibility, potential benefits and challenges of implementing a REDD+ program.The goals of this paper are to quantify the drivers of LCLUC dynamics, forest disturbance, and change trajectories in the study area over the last two decades.

Study Site
The study area covers an overall surface of 9.8 km × 9.8 km, centered over and including the 2,228 ha BNR.It is located approximately 40 km northwest of Toamasina and 45 km to the west of the Indian Ocean coast (Figure 1).The topography of the study area is characterized by rolling-to-hilly terrain with an elevation ranging from 2 m to 580 m.Mean elevation from sea level is 270 meters.The climate in the region is hot and humid with annual rainfall amounts of over 2000 mm and annual average humidity between 81% and 91%.The region has two seasons: a hot, rainy season from November to April and a cooler, drier season from May to October [6].The average annual temperature is 24 °C with annual lows of around 16 °C between June and August and annual highs of 32 °C possible between December and February.The soil is a red-yellow ferralitic.The dominant type of crop in the area is rice, which includes wet paddy rice being planted mostly at the bottom of valleys and dry hillside rice.Maize is grown in small patches along with sugar cane, chili and a few other crops, but these are far less common.

Data
Several field campaigns were conducted between 2010 and 2014.Data collected included land cover training samples, climate data, and location surveys of major invasive plant species.Thirty by thirty meter plots consisting predominantly of invasive plants were identified in the field, and mapped with hand-held GPS units.Additional training samples were photo-interpreted using high-resolution satellite images, i.e., GeoEye-1, IKONOS-2.
Various remote sensing data were collected from sensors including Landsat, Ikonos-2, Geoeye-1 and Hyperion (Table 1).All of these data were pre-processed, which included radiometric and atmospheric correction of passive optical images.Gap filling on Landsat ETM+ images collected after 2003 was conducted using local linear histogram matching.Prior to gap filling, cloud and cloud shadow were extracted using reflectance threshold values, and replaced using other images from the respective season.Please refer to the details of both field data collection and satellite data processing in authors' recent work [24].

Methodology
Figure 2 shows the workflow of the study, which is detailed in the following sections.

Land Cover Classification
Reference spectra of LCLU classes were developed by extracting from Hyperion image over pure patches identified during the field campaigns.All samples were verified against the 0.5 m orthorectified and pan sharpened Geoeye-1 image and screened relative to the field photos taken in the field.Additional training samples of agriculture, water and other LCLU classes were photo-interpreted from the image.Developed reference spectra were re-sampled to match sensor wavelengths, and were used as training endmembers for subsequent analysis.Figure 3 shows selected training samples using field, high resolution and middle resolution satellite imagery.It is possible to discriminate these land cover types in greater detail at higher spatial resolution; however, spectral confusions may occur between certain land cover types at the 30 m spatial resolution of Landsat and Hyperion.For example, residential areas with fallow, and agriculture with grassland appear similar in their spectra (Figure 4).The spatial extent of grassland land cover was very small in the study area [24], so we did not consider grassland as a separate land cover class in our classification.The classification scheme included eight classes: agriculture, evergreen forest, fallow, mixed forest, Molucca raspberry, residential, shrubland and water.Molucca raspberry, one of the 15 invasive plant species found in the reserve, can be easily identified at different resolutions due to its distinct spectral appearance.However, smaller isolated patches of Molucca raspberry were also found within the mixed forest stand that were confused with mixed forest to some extent in Landsat images classification.Therefore, Molucca raspberry was merged into mixed forest after classification.It is also worth noting that residential areas were extracted using object oriented feature extraction (OBIA) using high spatial resolution Geoeye-1 and IKONOS-2 imagery.Then, it was masked out during time series Landsat classification, thus enabling a better delineation of fallow land cover.Training samples were selected of prominent land cover types which were consistent on the Landsat images from 1990 to 2010.The use of the same training samples on multi-temporal image classification improves consistency and comparability of classification results.Fallow land cover is the area used for the production of crops that do not exhibit visible vegetation as a result of being tilled in a management practice that incorporates prescribed alternation between cropping and tillage.It can be identified on a single date image as a harvested bare soil, cleared land for slash-and-burn, or landslide-slope area triggered by heavy rain; however, it is subject to change on multi-temporal images.Therefore, fallow training samples were determined on each classification image, respectively.Fallow land cover class was added to agriculture during the post-classification process.The Terrain Categorization (TERCAT) tool available with ENVI 5.0 (ENVI® image processing and analysis software, from Exelis Visual Information Solutions) was employed to create an output product in which pixels with similar spectral properties were clumped into classes.The TERCAT tool provides most of the standard classification algorithms in addition to an additional algorithm called Winner Takes All (WTA).WTA is a voting method that classifies pixels based on the majority compiled from all of the other classification methods that were conducted [24].We used the following supervised classification methods in our TERCAT process: Maximum Likelihood, Spectral Angle Mapper (SAM), Minimum Distance, Mahalanobis Distance, and Support Vector Machine (SVM).To determine the best possible combination of input parameters specific to each classification method, inputs were determined in an iterative process that previews the classification results and compares them with training samples.
The classification accuracy was assessed with producer's and user's accuracies and kappa coefficient.The producer's and user's accuracies are measures of omission and commission errors [30].The kappa coefficient is a measure of agreement that accounts for the rate of correct classification occurring by chance.

Grid Cell Processes
LCLUC analysis often focuses on the overall trend and qualitative description of spatial change over a period of time.Quantitative spatial analysis of LCLUC characteristics over time and space is of paramount significance for understanding the impacts of human activity and natural environmental change [31].Grid cell method provides a quantitative measure of spatial-temporal dynamics of LCLUC [32,33].It allows the researcher(s) to analyze different spatial characteristics of LCLUC within and outside BNR to further reveal anthropogenic disturbance on the ecosystems.Our study area is approximately 100 square kilometers.Based on multiple experiments, we decided to use a 100 m square grid cell.The reasoning for this is: the smaller the grid cell size, the finer the spatial details of change observed from the grids; however, the cell size should be greater than the size of 30 m Landsat pixels since the goal here is to quantify percentage change of land cover types within a landscape unit.The attribute of a 30 m Landsat pixel is a single land cover type.Aggregation of surrounding pixels to a larger grid allows the estimation of percent cover of a land cover type within the grid (e.g., 100 m grid cell).
Vector shapefiles are best suited for the depiction and analysis of discrete objects in space while raster grids are most effective for mapping the qualities of space itself [34].Equally sized vector grid cells combines the advantages of both vector and raster representation of space, providing a robust approach for spatial statistics.For example, a classification based change detection using raster data results in a binary map which shows the trajectory of change between the land cover types; with the grid cell approach, the degree of change of a specific land cover can be demonstrated (e.g., in percent) in addition to the trajectory of the change.We created a blank vector polygon square 100 meter grid cells using fishnet tool of ArcGIS 10.1 software.Land cover classes were assigned to empty grid cells using spatial join.The percentage of a specific land cover within each cell was calculated and added to a new attribute table.It is worth noting that converting the land cover data into a square grid is not straightforward, especially at the grid cell boundary where different land cover types intersect.For pixels located on a boundary grid line or spanning multiple grid lines, it is difficult to judge which grid these pixels should be assigned to.In most cases, these pixels were assigned to the cell that covered the largest part of the pixel [32].The problem with this method is that the sum of all land cover areas within one grid can exceed the area of the grid cell itself [33].To avoid this problem, we converted each raster LCLU map to a vector map, then intersected each land cover type of LCLU vector map with 100 square meter empty grid cell map, separately.This process divides the land cover patch in LCLU vector map into pieces so that it fits into the grid cells proportionally.Then, we sum the area of each of the patch pieces (Ai) within the grid cell and calculated the percentage (Pi) of each land cover in a grid following the equation below.
where, i is the piece of the patch that falls in a grid cell, Ai is the area of the land cover within the piece of the patch, Ag is the total area of the grid cell (100 square meters in our case).

Moran' I
Spatial autocorrelation is an indicator of the degree of clustering, randomness, or fragmentation of a pattern.Moran' I is designed to detect spatial autocorrelation [35].The value of Moran' I ranging from −1 to 1, the values −1, 0, and 1 indicating negative spatial autocorrelation, randomness and positive spatial autocorrelation, respectively.It is positive when the observed value of locations within a certain distance or their contiguous locations tends to be similar (clustered), negative when they tend to be dissimilar (dispersal), and approximately zero when the observed values are distributed randomly and independently over space.

Forest Disturbance Mapping
The forest disturbance index (DI) [36] is a metric indicative of possible drivers of forest change, which provides another perspective of forest degradation in addition to the spatial statistics of land use change.It can be used to map windfall disturbance and timber collection infractions.We calculated DI using Landsat data from 1990 to 2011.The DI is computed using linear combination of normalized Tasseled Cap Transformation bands, i.e., brightness (Br), greenness (Gr) and wetness (Wr) [37] as shown in the equation below.
(2) Br, Gr, Wr is proven capable of measuring variations in soil background, vegetation vigor and senescence, respectively [37].The DI index is based on the fact that disturbed forest patches generally exhibit higher brightness, lower greenness and wetness due to the exposure of undercanopy soil signal from disturbance compared to undisturbed areas.Positive values indicate forest disturbance.The DI has been successfully tested in forest disturbance mapping [36,38,39], as well as differentiating windfall disturbance from selective logging [40] in different study regions across the globe.

Temporal Dynamics LCLU from 1990 to 2011
Figure 5 shows the land cover changes from 1990 to 2011 derived from Landsat TERCAT classification.Most land cover types displayed a coherent trend of increase or decrease, but some land cover shows irregular trend profiles.These anomalies can be attributed to the seasonal differences in rainfall, vegetation growth and agricultural activities during the time the data for each year was collected.As mentioned in Section 2.1, rain can be expected for more than 300 days per year, but is most intense between December and February during which time the average annual temperature is the highest around 32 °C.Average monthly rainfall in December can reach 370 mm, whereas average monthly rainfall in July can drop to 290 mm.In addition, there was minor cloud contamination in at least one of the Landsat images collected on 10 June 2007.Considering these factors, we selected images of 19 June 1990, 16 May 2004, and 1 May 2010 for final analysis (Figure 6).As shown in Figure 6, agriculture clearly increased at the cost of evergreen forest and mixed forest.Shrubland significantly increased probably owing to environmental degradation from human activities, as well as from emerging drought regimes in the region because of significant increase in temperature and a steady decline of rainfall [29].Mixed forest class also exhibited a declining trend, which might have been converted to shrubland and/or agricultural lands.While the temporal analysis provided an overall outlook of LCLUC in the study area, the results of the grid cell approach (see Section 4.2) revealed the causal relationship between land use changes and spatial details of change in space between different land cover classes.Overall accuracy, Kappa coefficient, producer's accuracy, and user's accuracy were calculated using the training samples and field survey plots as the "ground truth" data.As shown in Table 2, classification results show overall very good agreement.Overall accuracy of the classification maps of 1990, 2004 and 2010 were 75.5%, 82.6%, 86.8%, and kappa = 0.71, 0.79, 0.84, respectively.As mentioned in Section 3.1, fallow land cover is actually harvested cropland or agricultural fields left unplanted during a specific time.Its spectra are very similar to bare soil spectra, which is drastically different than vegetation spectra.For this reason, fallow training samples were determined on each classification image separately, and analyzed as an individual land cover for accuracy assessment.It should be noted that the fallow class was merged into agriculture land cover for spatial statistics of decadal LCLUC analysis.Both producer's and user's accuracy to classify agriculture, mixed forest and shrubland was lower when compared to other categories of land cover, particularly for the1990 and 2004 images (Table 2).Potential explanations for why the three categories were difficult to distinguish include: (1)

Spatial-Temporal Change of LCLU
In order to understand the spatial distribution of land cover changes, we created grid cell maps by overlying land cover maps on the empty grid cells and calculated the percentage for each cell of the five land cover types within it.Figure 7 shows the grid cell based spatial changes of agriculture from 1990 to 2010.The value of each grid cell of the change map was calculated by subtracting the agricultural area of 1990 from that of 2010 within the same grid, and then dividing the changed area by the cell area.As illustrated in Figure 7, agriculture increased in most of the study areas except for within the limits of BNR, where it is illegal.The increase of agriculture was most obvious in the northwest part of the study area (replacing the evergreen forest), and in the close proximity of the reserve.Evergreen forest decreased mostly over the outside of the protected zone, especially in the areas where the agriculture increased (Figure 8).A mix of decrease and increase in evergreen forest was observed within the reserve.The decrease could be attributed to the degradation of primary forest from windfalls, illegal logging, encroachment of invasives [29], and/or cyclone damage while the increase was explained by forest growth over time.Figure 9 shows the spatial change of mixed forest during 1990-2010.The obvious trend was that the most increase took place within the reserve, and significant reductions occurred in areas surrounding human settlements indicating strong anthropogenic influence on mixed forest.The spatial and temporal dynamics of shrubland from 1990 to 2010 were mostly manifested outside the reserve (Figure 10).The shrubland land cover class increased especially in the north and northeast of the study region where dramatic decrease in mixed forest was observed.This might indicate the degradation of environment due to slash-and-burn agriculture, increased firewood collection by a growing human population and abiotic factors, e.g., changes in precipitation [29].

Quantitative Relationship between Land Cover Categories
In order to further explore spatial characteristics of LCLUC, we calculated Pearson's correlation coefficient (r) among four land cover categories based on the grid cells.Table 3 presents a summary of the coefficient matrix between changes of land cover categories from 1990 to 2010.Agriculture change from 1990 to 2010 was negatively correlated with all other land cover classes with relatively strong significance.The correlation coefficient between the agriculture and mixed forest and shrubland were −0.47 and −0.45, respectively, which indicated agriculture increased at the cost of shrubland and mixed forest to some extent.Evergreen forest was negatively correlated with other land cover types as well, and the relationship was relatively strong between evergreen forest and mixed forest (r = −0.30),which was followed by the next strongest correlation between evergreen forest and shrubland.The implication was that evergreen forest was converted to mixed forest and shrubland.Relatively strong and negative correlation was also found between mixed forest and shrubland.The correlation coefficient between mixed forest and shrubland was −0.31, which indicated mixed forest was converted to shrubland to some degree.The correlations mentioned above were all statistically significant at the 0.01 level.The possible change trajectories between these land cover types is illustrated in Figure 11.It should be noted that Figure 11 employs two examples illustrating potential forest degradation pathways and associated declines in biodiversity, ecosystem functions and services (see [41,42] for other examples).Intact forests damaged by fires, cyclones or landslides may be replaced by shrubland or grasslands without first passing through a mixed forest phase.Using multi-temporal and high resolution satellite images of our study area, we observed examples of primary or mixed forests cleared and converted to croplands and then unmanaged fallows (a common practice in the tavy system) which transitioned into grasslands.Eckert et al. [28] noted that most pristine forests in their study site were directly converted to non-forest formations with low carbon stocks as a result of the widespread practice of slash and burn agriculture in eastern Madagascar.Invasive species can contribute to habitat degradation through outcompeting resident species leading to reductions in biodiversity and negatively altering effected ecosystem functions [41,43].Ghulam [29] documented a significant increase in the penetration of invasive plant species along the edge and into BNR from 2005 to 2012.

Spatial Autocorrelation Analysis for LCLU Categories
We also detected the degree of clustering and dispersion of each LCLU category by calculating Moran' I.As shown in Table 4, there were several obvious trends found from the index.First, regardless of 1990 or 2010, Moran' I for each LCLU class was relatively high and positive, which indicated that each land cover was spatially clustered, particularly the agriculture class.This was because the agriculture land cover type was directly related to human activity, and, therefore, it was concentrated in space within the proximity of human settlements.Second, the value of Moran' I were increased from 1990 to 2010 for all of the land cover classes.The increase of Moran' I for agriculture implied that the increase of agricultural land mainly happened near or around the agricultural land in 1990 with a sprawl increasing pattern.As shown in LCLU maps, there was evergreen forest outside BNR, which disappeared over the past 20 years resulting in a limited and more concentrated distribution only within the reserve limit.This increased clustering in space, explained the increased value of Moran' I for evergreen forest.A similar process were found for mixed forest over the past 20 years.Overall, the area of shrubland was relatively small and the spatial distribution was relatively dispersed in 1990, but it significantly increased during 1990-2010, resulting in increased value of Moran' I for shrubland.

Forest Disturbance
Figure 12 shows the areal mean forest disturbance index (DI) derived from Landsat data (gray bars) from 1990 to 2011.Except for 1996, DI values were positive for all of the Landsat images used in this study, and both frequency and magnitude increased since 2004.The implication is that forest degradation has been intensified over the last several decades probably due to damage from windfall as much as from selective logging.This was consistent with recent publications that the region has been experiencing emerging drought characterized by increased temperature and decreasing rainfall [29].We did not attempt to distinguish windfall disturbance from selective logging as the focus here was documenting forest transition pathways and its potential impacts on biodiversity.It is worth noting that the reserve has been hit by a number of severe cyclones coming from the Indian Ocean since 2007.Additionally, periodic illegal timber collection inside the reserve has been observed in recent years.

Conclusions
This paper examined the drivers of land-cover/land-use dynamics in Réserve Naturelle Intègrale de Betampona (BNR) and surrounding areas from 1990 to 2011.Detailed land cover maps were developed using multi-source and multi-resolution images, and a grid cell based spatial statistical approach was used to quantify the drivers of land use changes.Forest transition pathways were depicted using Moran's I at 100 m grid cell.The results showed that grid cell based Moran's I approach is a powerful method for quantifying spatial and temporal dynamics of forest change.It was observed that agriculture

Forest index
significantly increased over the last two decades mostly over the patches of evergreen forest and mixed forest outside the protected area.Shrubland followed a similar trend to agriculture.A significant increase in mixed forest was observed within the limits of the reserve, indicating transition from evergreen to mixed forest.On the other hand, the negative correlation between mixed forest and shrubland indicates that mixed forest was converted to shrubland.These findings were echoed by the forest distance index calculated from Landsat data, which showed consistent and increased forest disturbance over the last two decades.

Figure 1 .
Figure 1.Location map of the study area.Land-cover/land-use (LCLU) map shown on the right is derived from IKONOS-2 multispectral imagery collected in 1 May 2010.The LCLU map is overlain on a 2 m resolution shaded relief map.

Figure 2 .
Figure 2. Overall workflow of the methodology.

Figure 3 .
Figure 3. Major land cover types at different spatial resolution.

Figure 5 .
Figure 5. Temporal dynamics of land use changes in the study area from 1990 to 2011.

Figure 6 .
Figure 6.Land cover classes derived over the same season in different years (adapted from [29]).

Figure 7 .
Figure 7. Percentage change of agriculture in each grid cell from 1990 to 2010.

Figure 8 .
Figure 8. Percentage change of evergreen in each grid cell from 1990 to 2010.

Figure 9 .
Figure 9. Percentage change of mixed forest in each grid cell from 1990 to 2010.

Figure 10 .
Figure 10.Percentage change of shrubland in each grid cell from 1990 to 2010.

Figure 11 .
Figure 11.A conceptual framework of ecosystem degradation in the study area showing transition pathways due to natural forces and human activities (adapted from [41]).Invasive plants occupy areas cleared from landslides and abandoned croplands.Evergreen forests are converted to mixed forest, agricultural lands, grassland and shrubland, and/or eventually to degraded shrublands.

Figure 12 .
Figure 12.Temporal dynamics of forest disturbance (DI).Negative DI means no disturbance and optimal forest growth, and positive values mean forest disturbance due to windfalls, logging or other abiotic stress.Mean DI value calculated over the BNR was negative in 1996 and positive for the rest of the sampling years.An increase in frequency and magnitude of the positive values since 2004 indicates significant forest disturbance.

Table 1 .
Remote sensing data used in this study.

Table 2 .
at Landsat 30 m spatial resolution, mixed forest with shrubland and shrubland with mature taller crops (e.g., maize) Accuracy assessment of classification results based on confusion matrix for 1990, 2004 and 2010 (pixels).

Table 4 .
Moran' I for LCLU categories in 1990 and 2010.