Integrated Landscape Change Analysis of Protected Areas and their Surrounding Landscapes: Application in the Brazilian Cerrado

Remote sensing tools have been long used to monitor landscape dynamics inside and around protected areas. Hereto, scientists have largely relied on land use and land cover (LULC) data to derive indicators for monitoring these dynamics, but these metrics do not capture changes in the state of vegetation surfaces that may compromise the ecological integrity of conservation areas’ landscapes. Here, we introduce a methodology that combines LULC change estimates with three Normalized Difference Vegetation Index-based proxy indicators of vegetation productivity, phenology, and structural change. We illustrate the utility of this methodology through a regional and local analysis of the landscape dynamics in the Cerrado Biome in Brazil in 2001 and 2016. Despite relatively little natural vegetation loss inside core protected areas and their legal buffer zones, the different indicators revealed significant LULC conversions from natural vegetation to farming land, general productivity loss, homogenization of natural forests, significant agricultural expansion, and a general increase in productivity. These results suggest an overall degradation of habitats and intensification of land use in the studied conservation area network, highlighting serious conservation inefficiencies in this region and stressing the importance of integrated landscape change analyses to provide complementary indicators of ecologically-relevant dynamics in these key conservation areas.


Introduction
Protected areas (PAs) currently cover around 15% of the world's land and inland waters and remain the cornerstone of global biodiversity conservation strategies [1]. However, the capacity of PAs to preserve biodiversity greatly relies on the ability to maintain the ecological processes that span beyond their boundaries [2,3].
The anthropogenic-driven landscape dynamics in the PAs and surroundings (hereafter referred to as 'interface areas') highly influence the ecological processes that extend across these areas, such as those related to wildlife movement (e.g., population dynamics [4] or gene flows [5]), the ecosystem Since 41% of the selected 61 CUs do not have, to date, an approved management plan (Appendix B), and a clear delimitation of the BZ is missing from the majority of management plans, we decided to use the 10 km radius as defined by the earlier, more restrictive, resolution for the delimitation of the BZs. Our resultant study areas extend over 90,747 km 2 , spanning across twelve Brazilian states. Figure 1 shows the location of the study areas, and Figure 2 shows the location and landscape in 2001 and 2016 of the local case study.    Since 41% of the selected 61 CUs do not have, to date, an approved management plan (Appendix B), and a clear delimitation of the BZ is missing from the majority of management plans, we decided to use the 10 km radius as defined by the earlier, more restrictive, resolution for the delimitation of the BZs. Our resultant study areas extend over 90,747 km 2 , spanning across twelve Brazilian states. Figure 1 shows the location of the study areas, and Figure 2 shows the location and landscape in 2001 and 2016 of the local case study.

Land Use and Land Cover Data
The land use and land cover (LULC) maps of 2001 and 2016 of the Brazilian Land Cover and Use Collection 4.0 Map Series of the MapBiomas project [46] were acquired for analyzing land use and land cover changes (LULCC) in the study areas (presented in Section 3.1). The potential of these remote Remote Sens. 2020, 12, 1413 5 of 36 sensing-based maps to locate and quantify large-scale LULC changes remains largely unexplored by the scientific community, and only a few recent studies have attempted to do so (e.g., [47][48][49][50]).
These maps, produced at a 30 m spatial resolution over the entire Brazilian territory, are obtained through annual automatic processing and classification of Landsat imagery. For our study, we used the levels 2 and 3 of the 2001 and 2016 classifications (a detailed explanation of the processing chain and validation methods used can be found in [51], and the results of the accuracy assessment are presented in Table 1). Table 1. Results of the accuracy assessment of the levels 2 and 3 of the 2001 and 2016 classifications of Collection 4.0 MapBiomas land use and land cover maps (for information on quantity and allocation disagreement statistics see [52] We extracted the data corresponding to our study areas from the original raster maps, and, to provide a synthetic overview of the trajectories related to the main LULC types, we reclassified the original maps to create a simplified version containing ten categories of the original classification ( Table 2). Due to pixels with missing values in the 2001 and/or 2016 layers, we were unable to analyze a total of 164 km 2 , which represents 0.18% of the study area.

Satellite Data
Part of our study aims to examine the vegetation dynamics of various vegetation cover types, including natural and anthropogenic vegetation cover. SRS vegetation indices have long been used for this purpose, the Normalized Difference Vegetation Index (NDVI) being the most widely used due to its strong relationship with vegetation biophysical parameters such as leaf area index, green leaf biomass, and leaf photosynthetic activity [53]. Although NDVI is affected by shadows, and is more affected by noise (e.g., aerosol effects) and by radiometric saturation over high biomass vegetation than other vegetation indices [54], it allows scientists to efficiently capture the vegetation dynamics of a wide range of vegetation cover types at short time scales. This is due to its high sensitivity to slight variations in the quantity and state of photosynthetic vegetation [55].
We used two different datasets of NDVI data: one dataset composed of MODIS time series, with a high temporal resolution to capture the vegetation dynamics at a fine temporal scale, and a complementary dataset consisting of high spatial resolution Landsat data. All the images were projected to the common Projected Coordinate System 'SIRGAS 2000/Brazil Polyconic (EPSG: 5880)'.

MODIS Data
Two annual time series of the MODIS Terra Vegetation Indices 16-Day Global 250 m product (MOD13Q1 v.5) [56] were obtained from the collection 'MODIS/MOD13Q1' of the Earth Engine Data Catalog through the Google Earth Engine (GEE) platform [57]. Although the MOD13Q1 version 6.0 collection extends the series of data beyond 2017 (when MOD13Q1 version 5.0 stopped being produced by the NASA Land Processes Distributed Active Archive Center), it currently presents issues related to "unexpected missing data, incorrect instances of No Data and spikes in the NDVI values and Usefulness band values not correctly assigned" (https://lpdaac.usgs.gov/products/mod13q1v006). We therefore downloaded the first and the last annual time series of the MOD13Q1 Version 5.0 product with a full-year coverage (23 composite images per year) for 2001 and 2016. This product is processed as 16-day syntheses of the daily surface reflectance series 'MOD09', and provides an NDVI band at 250 m resolution with minimized atmospheric and directional effects [56,58].
To reduce any potential remaining atmospheric effects, we applied the noise-reduction algorithm presented in [59] adapted from the iterative Savitzky-Golay smoothing algorithm developed by Chen et al. [60]. To ensure the noise reduction is also applied to the first and last images of the time series, we downloaded five extra images at both extremes of the two annual time series.
Finally, to linearly interpolate the smoothed time series pixel-wise, and place their values at regular intervals of 16 days, we used the Composite Day of Year (CDOY) band of the MOD13Q1 product, that contains the acquisition date of the pixels retained to produce the original composite images [56]. This step ensured a uniform sampling of the NDVI for each of the evaluated years. The different pre-processing steps applied to the MODIS data are shown in Figure 3.
To attain a valid NDVI range from 0.0 to 1.0, a scale factor of 0.0001 was applied to all the bands in the 2001 and 2016 annual NDVI time series.

Landsat Data
We used Landsat data to obtain high spatial resolution images representing the landscapes of each CU and their BZs in 2001 and 2016. The choice of Landsat scenes was made according to these four criteria: 1. Scenes had to be cloud-free over the studied areas; 2. The number of scenes per CU and its corresponding BZ was minimized; 3. If multiple scenes were needed to cover a CU and its BZ, the dates of acquisition between scenes needed to be the closest possible in time; 4. For the same CU and BZ, the acquisition date of the 2001 scene(s) needed to be the closest possible in time to the acquisition date of the 2016 scene(s).
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing of the time series, we downloaded five extra images at both extremes of the two annual time series. Finally, to linearly interpolate the smoothed time series pixel-wise, and place their values at regular intervals of 16 days, we used the Composite Day of Year (CDOY) band of the MOD13Q1 product, that contains the acquisition date of the pixels retained to produce the original composite images [56]. This step ensured a uniform sampling of the NDVI for each of the evaluated years. The different pre-processing steps applied to the MODIS data are shown in Figure 3. The original data derives from 5 MODIS tiles h12v10, h12v11, h13v09, h13v10, h13v11, which have been cropped to the extent of the study areas and mosaicked to a single image mosaic per date (i.e., 23 mosaics per year). 2 The final annual time series start and end on the day-of-year 1 and 353, respectively.
To attain a valid NDVI range from 0.0 to 1.0, a scale factor of 0.0001 was applied to all the bands in the 2001 and 2016 annual NDVI time series.

Landsat Data
We used Landsat data to obtain high spatial resolution images representing the landscapes of each CU and their BZs in 2001 and 2016. The choice of Landsat scenes was made according to these four criteria: 1. Scenes had to be cloud-free over the studied areas; 2. The number of scenes per CU and its corresponding BZ was minimized; 3. If multiple scenes were needed to cover a CU and its BZ, the dates of acquisition between scenes needed to be the closest possible in time; 4. For the same CU and BZ, the acquisition date of the 2001 scene(s) needed to be the closest possible in time to the acquisition date of the 2016 scene(s).
A total of 113 scenes were selected (the acquisition details and identification of the Landsat scenes are presented in Appendix C). Since no single Landsat mission was operational in both 2001 and 2016 (except for the Landsat 7 mission in which, however, the Enhanced Thematic Mapper Plus (ETM+) sensor is failing since 2003 [61]), the scenes were selected from three missions and sensors of the Landsat Data Continuity Mission: Landsat 5 Thematic Mapper (TM) and Landsat 7 ETM+ scenes for the year 2001 and Landsat 8 Operational Land Imager (OLI) scenes for 2016. Despite the inherent differences between the three sensors due to the design improvements and adjustments between successive missions, the latest Tier-1 Surface Reflectance collections are considered suitable for time-series analysis since they include the Landsat scenes with the highest available data quality, which are consistent and inter-calibrated across the different Landsat instruments [62][63][64].
The Near-Infrared (NIR) and Red (R) bands of the selected scenes were retrieved through the GEE platform from the following collections: In the GEE cloud-computing platform, the NDVI was subsequently computed from the NIR and Red bands of the 113 scenes that share a common spatial resolution of 30 m, and were mosaicked to a single 2001 Landsat TM/ETM+ mosaic and a 2016 Landsat OLI mosaic. Finally, the mosaics were cropped to a 300 m buffer from the study areas' limits (to avoid edge effects when calculating the landscape metric related to the spatial variability of NDVI, presented in Section 3.2).

Land Use and Land Cover Change Analysis
The reclassified 2001 and 2016 MapBiomas LULC maps (Section 2.2) were used to determine the LULC transitions in the study areas (presented in Section 4.1.1). The analysis was carried out from pixel-area estimates and separately for the CUs (total area evaluated = 2,755,413 ha) and BZs (total area evaluated = 6,302,902 ha), by applying two widely used tools:

•
Transition matrices (also known as cross-tabulation matrices) [65,66], in which each row represents the LULC category in 2001, each column the LULC category in 2016, each entry represents the surface area experiencing LULCC or maintaining the same LULC type, and the last column and row represent the total area covered by each LULC category in 2001 and 2016, respectively. • Sankey diagrams [67,68], derived from the transition matrices, in which the total area covered by each LULC category in 2001 and 2016 is represented as stacked vertical bars (the height of each bar representing the relative proportion of each LULC type in the study areas), and in which transition lines, connecting each sequential pair of bars, represent the surface area experiencing LULCC or maintaining the same LULC type (the thickness of each line is proportional to the total surface area experiencing the corresponding LULC transition or remaining unchanged). The 'sankeyNetwork' function of the 'networkD3' R package v.0.4. [69,70] was used to derive the Sankey diagrams.

Vegetation Change Analysis
To analyze the vegetation dynamics of the major vegetation cover types (i.e., 'Natural Forest', 'Non-Forest Natural Formation', 'Pasture', and 'Agriculture'), we first extracted the pixels belonging to these vegetation cover types in both years from the 2001 and 2016 MapBiomas LULC maps (i.e., the areas where the vegetation cover type remained unchanged).
To assess the dynamics of each vegetation cover type separately (i.e., measuring the "within-class" dynamics) and considering that the landscape metrics used for this purpose are calculated at the MODIS spatial resolution of 250 m, we narrowed down the analysis to exclusively 'pure pixels' (i.e., 250 m resolution pixels with a uniform vegetation cover type).
To this end, the mask of "unchanged vegetation cover type pixels" at the original MapBiomas 30 m resolution was overlapped by a regular square-grid fitting the 250 m resolution MODIS cells. The grid cells overlapping different vegetation cover types' pixels or any missing values were used as a mask, applied to all of the satellite imagery to discard 'mixels' and missing values for subsequent analysis. By doing this, we avoided the mixed signals from cells covering heterogeneous vegetation types and ensured that the observed vegetation dynamics corresponded exclusively to the targeted vegetation types under evaluation. A total of 438,348 cells, equivalent to approximately 27,397 km 2 were retained for the vegetation change analysis.
The masked satellite imagery was used to calculate the following three NDVI-based landscape metrics that provide complementary information on the spatiotemporal dynamics of the four major vegetation cover types.
A first landscape metric (NDVI inter ) was derived from the 2001 and 2016 annual MODIS NDVI time series to evaluate the magnitude and overall trend of the inter-annual change of the four vegetation cover types. This first landscape metric was calculated pixel-wise as the difference between the annual mean NDVI in 2016 and 2001 (Figure 4a shows a schematic representation of the calculation of NDVI inter ): where X 0 and X t are the images in the first and last annual time series, respectively, and n is the number of images in the annual time series. The valid range of NDVI inter is from -2 to 2.
Considering that the NDVI is a well-established and widely-used spectral-based proxy for above-ground vegetation productivity (e.g., [71][72][73][74][75][76]), NDVI inter was therefore considered here as an indicator of vegetation productivity change: negative values represent areas where the mean annual productivity has decreased, and conversely, positive values indicate a general increase in the mean annual productivity. (1) where 0 and are the images in the first and last annual time series, respectively, and is the number of images in the annual time series. The valid range of is from -2 to 2.  A second landscape metric (NDVI intra ) was derived from the 2001 and 2016 annual MODIS NDVI time series to evaluate the magnitude and overall trend of the intra-annual variability change of the vegetation cover types.
We first calculated pixel-wise the cumulative intra-annual variations of NDVI in 2016 (i.e., the sum of the absolute differences between each pair of consecutive images in the annual time series), and divided the result by the number of sequential image pairs in the annual time series (i.e., 22) to obtain an average estimate of the intra-annual variability in 2016. We then repeated the same operation for the 2001 series and calculated the difference between the 2016 and 2001 results.
This second landscape metric was calculated pixel-wise as follows (Figure 4b shows a schematic representation of the calculation of NDVI intra ): where X 0 and X t are the images in the first and last annual time series, respectively, and n is the number of images in the annual time series. The valid range of NDVI intra is from -2 to 2.
The NDVI intra provides an estimate of the mean intra-annual variability change of NDVI, which allows quantifying variations in the seasonal patterns of vegetation cover. Taking the definition of 'Land Surface Phenology', "the seasonal pattern of variation in vegetated land surfaces observed from remote sensing" [77], as a reference, the NDVI intra was considered here as an indicator of vegetation phenology change.
Unlike threshold-based or change-detection-based phenology metrics [78], NDVI intra does not aim to identify the changes in different phenological phases (e.g., start-of-season, green-up, peak of growing season, etc.) but rather to quantify the changes in the overall oscillatory pattern. For instance, negative values of NDVI intra indicate a decrease in the mean intra-annual variability of NDVI observations, due to either a decrease in the amplitude or the number of seasonal oscillations. Conversely, positive values indicate that seasonal oscillations have become more pronounced or numerous, therefore increasing the mean intra-annual variability.
A third landscape metric (NDVI spatial ) was derived from the 2001 Landsat TM/ETM+ and 2016 Landsat OLI mosaics to evaluate the magnitude and overall trend of the spatial variability change of the vegetation cover types.
We first spatially overlapped the Landsat images with the 250 m resolution grid and calculated the standard deviation of each set of Landsat pixels falling within each grid cell of 6.25 ha (i.e., approximately 64 Landsat pixels) for 2016 and 2001, separately. We then calculated pixel-wise the difference between the standard deviation values in 2016 and 2001. This allowed us to summarize the spatial NDVI variability change of the high-resolution Landsat data at the MODIS pixel size and bring this metric to a common spatial sampling scale with the other two metrics.
This third landscape metric was calculated as follows ( Figure 4c shows a schematic representation of the calculation of NDVI spatial ): where Y 0 and Y t are the set of high-resolution pixels within a coarse-resolution cell in the first and last year images, respectively, and n is the number of high-resolution pixels that fit within a coarse-resolution cell. The valid range of NDVI spatial is from ≈-1 to ≈1. The landscape metric NDVI spatial , which estimates the spatial variability change of NDVI, was used in this study as a proxy for vegetation structure change. Vegetation heterogeneity change (i.e., vegetation horizontal structure change) can be quantified through NDVI texture analysis techniques (e.g., [79][80][81]), which measure the variability of pixel values in a given area [82]. A highly heterogeneous vegetation cover, either in terms of the spatial distribution of plants or plant species composition therefore presents a high variability in NDVI, whereas a homogeneous vegetation cover will present a low variability [82].
The standard deviation of the NDVI in a given area provides a simple estimate of this spatial variability and has thus been used in previous studies to estimate the heterogeneity of different vegetation cover types (e.g., [76,83,84]) including landscapes in the Cerrado Biome from Landsat data [85,86]. Negative values in NDVI spatial therefore indicate a homogenization of the vegetation cover, and conversely, positive values depict areas that have become more heterogeneous.
The correlation between the three metrics for each vegetation cover type was calculated to test if the three metrics contribute with complementary and non-redundant information on the vegetation dynamics. The general trends of vegetation productivity, phenology, and structure change (as indexed by the landscape metrics) of each vegetation cover type at the regional level are presented in Section 4.1.2 and Appendix E. To facilitate their interpretation and determine any differences in the trends according to particular locations, their respective statistical distributions are represented in boxplots for the CUs and BZs, separately, and per Brazilian state. In addition, the local vegetation dynamics trends of the local case study are presented in Section 4.2.

Land Use and Land Cover Change Analysis Results
The LULC transition matrices of the 61 studied CUs and their BZs are presented in Appendix D and summarized in Table 3 and Figure 5.  In total, the absolute surface area changes (considering loss and gain) in the CUs and BZs summed up to 75,468 ha and 777,544 ha, respectively, but a total of 348,463 ha in the CUs and 1,522,656 ha in the BZs experienced some type of LULCC from 2001 to 2016. These LULCC represent 12% of the CUs' surface area and 24% of the BZs' surface area.
As shown in Figure 5a, the main LULCC that took place in CUs were changes between natural forests (i.e., natural savannah and forest formations) and non-forest natural formations (i.e., grassland and wetland vegetation formations). Accordingly, the natural vegetation (i.e., 'Natural Forest' and 'Non-Forest Natural Formation') inside CUs has experienced little absolute change in surface area between 2001 and 2016 (Table 3). Moreover, in CUs, a large amount of pastures evolved to natural vegetation in CUs (46,883 ha), and a large amount of natural vegetation in the CUs (51,787 ha) was transformed to farming land (i.e., 'Pasture', 'Agriculture', and 'Mosaic of Agriculture and Pasture'), mostly pastures (32,010 ha) but also agriculture (19,741 ha).
In BZs, the changes between natural forests and non-forest natural formations were also widespread ( Figure 5b), but the most important LULCC (in number of hectares) was the conversion of natural forests to pastures (202,690 ha). This decline of natural forests in BZs has been considerably counterbalanced by a large amount of pastures being converted to natural forest (167,162 ha). These relatively balanced and opposed conversions, which also took place between non-forest natural formations and pastures (Figure 5b), result in relatively little net loss in natural vegetation in the BZs (Table 3). Agriculture, on the other hand, experienced a considerable surface gain (Table 3). This important expansion of seasonal and perennial crops was mostly due to conversion of pastures, but there was also substantial conversion of natural vegetation to agriculture (Figure 5b).
In total, agriculture experienced the greatest expansion, and pastures experienced the greatest surface loss in both the CUs and the BZs, due, to a great extent, to their conversion to agriculture but mostly to natural vegetation (Table 3). Forest Plantations and Urban Infrastructures also expanded considerably in BZs (Table 3), however, since these, and other LULCC were relatively small in comparison with the aforementioned major transitions, further discussion will be focused on the results of the latter.   Despite the changes, the great majority of the surface in the study areas did not experience LULCC ( Figure 5). The four major vegetation types (i.e., Natural Forest, Non-Forest Natural Formation, Pasture, and Agriculture) remained unchanged in 75% of the surface area (6,821,747 ha) between 2001 and 2016. For example, 90.6% of the natural forests in the CUs and 81.5% in the BZs did not experience LULCC, and 87.7% of the non-forest natural formations in the CUs and 71.2% in the BZs also remained unchanged. Figure 6 shows the distribution of the three landscape metrics for the 438,348 cells retained for the vegetation change analysis (see Section 3.2). As expected, in all three cases, the actual empirical range of the metrics was considerably lower than their theoretical range (see Section 3.2), the latter being calculated considering maximal theoretical variability in NDVI, which is virtually impossible to reach in real conditions. Out of the analyzed cells, 57.3% belonged to natural forest areas, 27.3% to non-forest natural formations, 8.3% to pastures, and 7.1% to agricultural land. For the four vegetation types, no important relationships were found between the three landscape metrics: nearly all the landscape metrics pair-wise Pearson correlation coefficients (r) were between -0.1 and 0.1, except for the and correlation for the pastures, the and correlation for the natural forests (r = 0.2), and the and correlation for the natural forests (r = 0.4). All correlations were statistically significant with a p-value < 0.001. This suggests that each landscape metric provides complementary information on their vegetation dynamics.

Vegetation Change Analysis Results
To provide a concise overview of the regional-level vegetation dynamics' results, we only exhaustively present and discuss the results for the Natural Forest category. The results for the other three vegetation types are presented in Appendix E. The distributions of the landscape metrics' values for the Natural Forest are presented in Figure 7, for the whole study area (CUs and BZs, separately) and for each of the twelve Brazilian states (CUs and BZs combined). Out of the analyzed cells, 57.3% belonged to natural forest areas, 27.3% to non-forest natural formations, 8.3% to pastures, and 7.1% to agricultural land. For the four vegetation types, no important relationships were found between the three landscape metrics: nearly all the landscape metrics pair-wise Pearson correlation coefficients (r) were between -0.1 and 0.1, except for the NDVI intra and NDVI spatial correlation for the pastures, the NDVI inter and NDVI spatial correlation for the natural forests (r = 0.2), and the NDVI inter and NDVI intra correlation for the natural forests (r = 0.4). All correlations were statistically significant with a p-value < 0.001. This suggests that each landscape metric provides complementary information on their vegetation dynamics.
To provide a concise overview of the regional-level vegetation dynamics' results, we only exhaustively present and discuss the results for the Natural Forest category. The results for the other three vegetation types are presented in Appendix E. The distributions of the landscape metrics' values for the Natural Forest are presented in Figure 7, for the whole study area (CUs and BZs, separately) and for each of the twelve Brazilian states (CUs and BZs combined).   The NDVI inter metric, showed a general productivity loss in both the CUs and the BZs (72.4% of the analyzed Natural Forest cells in the BZs experienced a decrease in the annual average of NDVI and 71.9% in the CUs) (Figure 7a). This trend was widespread across the CUs and BZs of nearly all the states, except for Distrito Federal and Bahia (where nearly 70% of their Natural Forests experienced an increase in productivity), and was especially pronounced in the CUs and BZs of Piauí and Tocantins, where 96.3% and 89.3% of the natural forests have loss productivity (Figure 7b).
The NDVI intra metric showed a slight trend towards a more stable NDVI throughout the year for the natural forest areas (59.7% and 54.7% of the analyzed area in the BZs and CUs, respectively, experienced a decrease in the intra-annual variability of NDVI) (Figure 7c). The natural forests located in the CUs and BZs of Bahia were the most representative of this trend (with 73.8% of the forests obtaining negative NDVI intra values), and on the other hand, the forests in Piauí followed a clear opposite trend (83% of the forests experienced an increased intra-annual variability of NDVI) (Figure 7d).
Furthermore, the general negative trend in the spatial variability of NDVI measured by NDVI spatial (Figure 7e) suggested that the canopy of natural forests was, in general, experiencing homogenization. This trend was followed by 96.8% and 92.4% of the forests in Paraná and São Paulo, respectively (Figure 7f). A relatively smaller proportion of the forests in the other states follow this trend, except for those in Bahia, which experienced the opposite trend at a larger proportion (i.e., the forest cover is becoming more heterogeneous) (Figure 7f).

Local Case Study Results
The Parque das Nascentes do Rio Taquari and its BZ are dominated by natural forest, pasture, and agriculture areas (Figure 8a). The great majority of these areas have not experienced LULCC between 2001 and 2016 (Figure 8b), except for some forested areas having been cleared and transformed into pasture (these areas, characterized by a loss in productivity, are highlighted by low NDVI inter values e.g., top zoom of Figure 8c). Nevertheless, the results of the three landscape metrics over the apparently stable areas show that different vegetation dynamics have taken place (Figure 8c-e). The three landscape metrics are represented in Figure 8 as a continuous gradient over all types of LULC (depicting both vegetation dynamics related to LULCC processes and other processes) but only the masked satellite imagery (Section 3.2) is analyzed.
The NDVI inter results (Figure 8c) reveal that pastures have experienced different productivity changes locally, 40.4% have increased their average annual NDVI while 59.6% have experienced loss in productivity (the bottom-left zoom shows an example of loss in average annual NDVI of nearly 0.1). Regarding forests, productivity trends also vary locally. As for agriculture, NDVI inter results (Figure 8c) show that the majority of areas (67.6%) have increased their productivity (as in the bottom-left zoom, where the mean annual NDVI has increased by nearly 0.2).
The NDVI intra results (Figure 8d) show that the magnitude of phenology change was the highest in the agricultural land, in both positive and negative directions. This result reflects a change in the number of annual crop cycles (shown in the zoom images and annual temporal profiles of NDVI in Figure 8d). High positive NDVI intra values (in red) thus reflect an increase in the number of cropping cycles and high negative values (in blue) show the opposite (Figure 8d).
The NDVI spatial results (Figure 8d) show an overall increase in the spatial variability of NDVI over the pasture land (68% of the pastures have positive NDVI spatial values, indicating a general increase in heterogeneity). They also show an overall decrease in the spatial variability over the agricultural areas, where 60.4% is becoming more homogeneous, despite the contrasted changes over the edges of the field plots (related to changes in their size and shape). However, the clearest trend reflected by the NDVI spatial metric concerns the homogenization of the forest cover (represented by negative values). This trend, which also represents the majority of forest areas at the regional level (Section 4.1.2), represents 92.5% of the forest in this local case study, and suggests general densification of the forest canopy (zoom image in Figure 8d).

Discussion
Our study illustrates the potential of combining LULC change-detection with VI-based change-detection methods to derive spatially explicit and complementary information on the landscape changes taking place in PAs and their surrounding zones. Additionally, our combined approach has supplemented current knowledge of the Brazilian Cerrado's interface areas' landscape changes (e.g., [33,34]) with new information on the types of LULC and vegetation dynamics. We discuss these in Sections 5.1 and 5.2, and in Section 5.3, we address some operational and research perspectives of our work.

Land Use and Land Cover Change Analysis in the Cerrado Biome's Protected Areas and Buffer Zones
The regional-level LULCC analysis showed relatively little net loss in native vegetation, (especially in the CUs) but underlying high levels of conversion, which raises concerns about the general conservation effectiveness of these core PAs and their legally managed BZs. The major LULC conversion between 2001 and 2016 in these areas, other than the conversion between different natural vegetation types, was related to the transformation of natural vegetation to farming land: 51,787 ha in the CUs, despite this conversion being forbidden in strictly protected CUs, and 446,155 ha in the BZs.
These natural habitats were mainly cleared to install new pastures for cattle ranching; however, the conversion to crops was also responsible for a non-negligible portion of the deforestation. Even if a nearly proportional amount of natural vegetation has increased over former pastures, especially in CUs, secondary savannas recovered from abandoned pastures hardly return to the biodiversity levels of old-growth savannas without active restoration interventions [87]. The proposed methodology allows us to locate these areas, but complementary in situ monitoring would be relevant to investigate how these areas are being managed and to assess their levels of biodiversity in comparison with old-growth savannas.
The LULCC results have also highlighted the fast expansion between 2001 and 2016 of agricultural land in the BZs (increase rate of 94%) and to a lesser extent within the CUs limits (around 22,000 ha in 15 years). This result echoes the alarming rate of agricultural expansion in the Cerrado Biome as a result of increasing international demand for commodity crops [29,88,89]. Given that most areas in the BZs correspond to private landholdings, that a considerable proportion of the CUs has not been yet expropriated [90], and that many CUs still do not have approved management plans, this result may suggest that the regulation of the landowners' decisions to secure conservation targets in these areas has been greatly limited. Yet, considering the rate of habitat loss in the sustainable use CUs reported by Françoso et al. [34], the situation for the whole PA-network may be even worse than for the 61 strictly protected CUs investigated here.
Complementary assessments of these landscape dynamics through comparative analysis of the LULCC within interface areas and the rest of the Cerrado would provide further insight into the effectiveness of these areas to reduce landscape composition transformation as compared to unprotected land. Despite misclassifications, which are common in LULC data (especially between classes with similar spectral properties), the regular and open remote sensing-derived LULC datasets provided by the MapBiomas (https://mapbiomas.org/) and the TerraClass Cerrado [91] initiatives offer promising tools for this type of large-scale analyses.

Contribution of the Vegetation Change Analysis
The results of the vegetation change analysis offered complementary indications on environmental changes related to processes other than LULCC. Even if ground-surveys would ultimately be needed to explain the local biophysical processes behind the observed trends, some general hypotheses may be drawn from the main results.
The general Natural Forest productivity loss (as measured by the NDVI inter ) in both in the CUs and BZs suggests a general decline in above-ground biomass, which could reflect a potential degradation of the savannah and forest ecosystems functions. Bustamante et al. [92] indicated that a decrease in the primary productivity of the Cerrado was expected due to the projected climate and fire-regime changes. Besides, the homogenization of the canopy cover, shown by the NDVI spatial results may reflect the natural progression of the Cerrado vegetation from open to more closed formations with a higher density of woody plants [93]. This pattern was observed in the local case study (Figure 8e) and may suggest a general natural evolution towards late successional stages and forest encroachment (e.g., [94]).
Furthermore, the local case study results have illustrated the potential of the landscape metrics to monitor subtle changes in farming areas in a spatially explicit way. For example, the metrics highlighted some dynamics that were also predominant at the regional level, such as a global increase in crop productivity (Figures 8c and A4b), and a prevalent homogenization of the agricultural land (Figures 8e and A4f), which may ultimately reflect a general intensification of the agricultural practices. Moreover, NDVI intra allowed to identify changes in the number of crop cycles (Figure 8d). This metric may, therefore, be relevant to identify areas with an increase in the number of annual cropping cycles, as a complementary indicator of intensification of cropping systems (e.g., [95,96]).
Finally, the results showed a particularly high dispersion of the landscape metrics in the BZ agricultural areas, as compared to the other evaluated vegetation types ( Figure A4). This may indirectly reflect a wide variety of agricultural practices, which likely results in highly differentiated protected-area-agriculture interface landscapes. Additionally, the contrasting distributions per state of some landscape metrics (e.g., Figure A2b), show strong regional variability in the dynamics of certain vegetation types. To assess potential differences in the management of these areas (e.g., [34]), it is thus important to conduct comparative landscape change analyses at the state level, as well as by type of CU category or jurisdictional level.

Perspectives for Future Research in Remote Sensing and Interface Areas
Most studies on interface areas are still mostly using LULCC analysis to monitor landscape dynamics (e.g., [8][9][10][11][12][13][14][15][16][17]). By combining LULCC analysis with VI-based change detection, we were able to characterize LULCC and explore vegetation changes in areas where no LULCC was observed, thus offering complementary information. This complementarity was particularly evident in our studied areas, since 75% of the vegetation cover had not experienced LULCC during the study period, but had experienced changes in phenology, productivity, and structural changes. It was further evident in the local case study, where relatively little LULCC obscured substantial change in pastures, forests, and cropland.
Our methodology, therefore, offers a practical and promising avenue for producing more detailed assessments of landscape dynamics in interface areas. Integrated methodologies, such as the one presented, can provide a more complete picture of functional landscape changes, which can, in turn, improve our understanding about the links between landscape dynamics and ecological processes.
Beyond expanding the change detection capacity by allowing us to detect conditional change in addition to transitional change [97], the presented NDVI-based metrics also offer the advantages of pixel-based analyses. These advantages include temporal and spatial robustness, fine-scale observations and a gradient-based representation of landscapes (especially relevant for fine structure change detection [98]). However, unlike LULCC analysis, which provides simple and intuitive information, field surveys are required to obtain precise interpretations of the local processes behind observed vegetation dynamics. The NDVI-based metrics results can, nevertheless, help orientate field surveys to sample areas where particular vegetation dynamics are observed or over a gradient of dynamics to gain insight on the underlying social-ecological processes and their spatial distribution.
Through this procedure, NDVI-based metrics can offer valuable indications on some processes that may be relevant for ecological assessments such as habitat simplification and degradation, or intensification of agricultural practices. They can, therefore, help locate zones in the interface areas that could be potentially experiencing loss in ecological functionality or could be prone to future land conversion or degradation. This study, therefore, builds on the recent research exploring the use of SRS-derived products as indicators of ecosystem function [99][100][101] and Essential Biodiversity Variables (EBVs) [102][103][104][105] and could be potentially integrated into ecosystem services assessments [106], ecosystem models [107], species distribution models (e.g., [100,108,109]), movement ecology studies (e.g., [110]), and even social-ecological systems analyses (e.g., [111]).
From an operational point of view in the Brazilian context, this study shows the potential of using existing SRS-derived LULC maps from the MapBiomas project (https://mapbiomas.org/) to detect illegal deforestation within 'strictly protected' CUs. Alert systems such as the one proposed by the MapBiomas initiative could include analyses at the CUs level and BZs to offer an operational tool to alert the official managing agencies of these PAs. Moreover, the proposed methodology may represent a significant contribution to current information systems both in the interface areas in Brazil (where it could be implemented in monitoring systems of the Chico Mendes Institute for Biodiversity Conservation to complement in situ monitoring [112]), as well as in other world regions, since existing systems are still mostly focused on monitoring habitat loss in forest areas [13].
The new cloud-computing platform GEE [57] offers unprecedented computing possibilities that facilitate the reproduction of the proposed methodology. The applicability of the methodology and its broad-scale operationalization is, however, limited by the availability of appropriate input data. Especially, timely and reliable LULC maps must be available for the analysis, and while these products are increasingly available, they are still missing in many regions and require a high-investment to produce them. Additionally, the surface area over which the vegetation change analysis can be carried out is limited by the spatial resolution of the landscape metrics and the LULC data (e.g., in our study, a mask had to be applied to discard 'mixels' for the vegetation change analysis, due to the spatial resolution mismatch of landscape metrics at 250 m and LULC data at 30 m).
In particular, the data from the Sentinel-2 constellation of the Copernicus program, which are available at 10 m spatial resolution and 5 days temporal resolution, from 2017 onwards, alone or in combination with other high spatial resolution data (e.g., [113]), are key products to overcome the limitations of 'mixel' issues and will allow the monitoring of landscape dynamics over a larger extent, including small, fragmented landscapes. Additionally, the new Landsat Tier-1 Collections [63,114] offer an unprecedented opportunity for high-resolution time-series analyses with the improved usability and consistency of the Landsat archives and could facilitate the upgrading of bi-temporal change detection to continuous monitoring through time series or trajectory analysis [97]. The proposed landscape metrics could provide valuable information on the exact timing of major landscape changes and on subtler changes in the ecosystems' conditions if they were calculated over successive years for multiple-year periods. However, to reveal any additional processing that may be required for time series analysis, studies evaluating the radiometric consistency of the vegetation indices and their spatial variability from one instrument to another are still needed (e.g., [115][116][117][118]).

Conclusions
This study illustrates the potential of combining SRS-derived LULC data and NDVI-based landscape metrics to obtain complementary, spatially explicit information on the landscape changes in protected areas and their surrounding landscapes (interface areas). The presented methodology allows the localization of major long-term land conversions in a network of interface areas in the Cerrado Biome and the detection of trends in the productivity, phenology, and structure changes of the major natural and anthropogenic vegetation cover types. This study allows us to fill a knowledge gap in this geographic region in particular and opens new perspectives of the analysis of interface areas in general. Ultimately, this work contributes to the development of satellite remote sensing tools to improve the monitoring capacity and guide and support conservation planning in key conservation areas.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A.
The rainfall anomalies' analysis was performed with the historical data from the Meteorological Database for Teaching and Research (BDMEP) of the Brazilian National Institute of Meteorology (INMET) [119]. This database contains the data collected by the 'conventional meteorological stations' from the INMET stations network. The rainfall anomalies' analysis was carried out in 5 steps:

1.
We collected the monthly total precipitation data (mm) of the BDMEP meteorological stations that had a nearly complete 30 year time series of data from 1986 to 2016 (with no more than 1 consecutive month of missing records in the time series). 2.
The data was transformed to spatial point data, using the meteorological stations' geographic coordinates, and the data from the meteorological stations closest to the 61 conservation units (CUs) analyzed in this study was retained: the data from a total of 36 meteorological stations, spanning from 60.0 • W to 35.0 • W and 0.1 • S to 23.5 • S, was selected. 3.
The missing data gaps in the time series were temporally interpolated with the average value of the previous and next months and the monthly data was aggregated to retrieve the annual total precipitation (mm) for each year. 4.
The Inverse Distance Weighting (IDW) algorithm implemented in the Interpolation Plugin of the Quantum GIS (QGIS) software [120] was then applied to the annual total precipitation data to obtain spatially continuous annual spatial interpolations over the whole study areas at a 5 km spatial resolution. The Zonal Statistics Plugin [120] was then used to extract the total annual precipitation values for the 30 year period for each CU centroid.

5.
To evaluate the rainfall conditions in the two analyzed years (2001 and 2016), we used the climatological normal (i.e., the average value of a meteorological element over 30 years) as a reference base for the calculation of the Modified Rainfall Anomaly Index (mRAI) [121]. This index, based on the Rainfall Anomaly Index (RAI) developed by van Rooy (1965), is well adapted to short time periods of precipitation data [121], and thus to the length of our time series, and allows to classify the rainfall anomalies into 9 classes (Table A1). The mRAI was calculated as follows: where, P i is the total precipitation sum for the year i, P is the median of all annual precipitation sums of the 30 year period, SF is the scaling factor of 1.7 for positive anomalies P i ≥ P and −1.7 for negative anomalies (P i ≤ P), and E is the mean of the 10% most extreme precipitation annual sums. For positive anomalies (P i ≥ P), E is the mean of the annual precipitation sums of the 3 most rainy years of the 30 year period; for negative anomalies (P i ≤ P), E is the mean of the annual precipitation sums of the 3 least rainy years of the 30 year period.
The results show a significant positive correlation between the total annual rainfall in 2001 and 2016 for the 61 conservation units (Figure A1a), which suggests that both years received similar rainfall amounts. In addition, the two evaluated years present similar mRAI results ( Figure A1b), both years being classed normal to moderately dry with regards to the 30 year period and thus do not present any extreme rainfall anomalies.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing The mRAI was calculated as follows: where, is the total precipitation sum for the year , ̅ is the median of all annual precipitation sums of the 30 year period, is the scaling factor of 1.7 for positive anomalies ( ≥ ̅ ) and -1.7 for negative anomalies ( ≤ ̅ ), and ̅ is the mean of the 10% most extreme precipitation annual sums. For positive anomalies ( ≥ ̅ ), ̅ is the mean of the annual precipitation sums of the 3 most rainy years of the 30 year period; for negative anomalies ( ≤ ̅ ), ̅ is the mean of the annual precipitation sums of the 3 least rainy years of the 30 year period.
The results show a significant positive correlation between the total annual rainfall in 2001 and 2016 for the 61 conservation units (Figure A1a), which suggests that both years received similar rainfall amounts. In addition, the two evaluated years present similar mRAI results ( Figure A1b), both years being classed normal to moderately dry with regards to the 30 year period and thus do not present any extreme rainfall anomalies.

Appendix B.
List of the 61 conservation units (CUs) in the Cerrado Biome, designated before 2001 and with a 'strictly protected' status (Table A2).      The non-forest natural formations experienced divergent dynamics in the BZs and CUs in terms of productivity change, as shown by the opposed distribution of their values ( Figure E1a). In the BZs, 58.9% of the grassland and wetland areas were characterized by a decrease in the annual average of NDVI, whereas in the CUs, 60.4% of these areas experienced an increase in NDVI. However, in Figure E1b, we can observe a clear disparity of trends depending on the location of these areas. For example, the near totality (90.5%) of this vegetation type in Piauí has experienced a decrease in productivity, while the opposite trend can be observed in Bahia, where the non-forest natural vegetation productivity has increased in 92.8% of the sampled cells.
Similarly, the results for the metric were diverse according to the location of the non-forest natural formation areas ( Figure E1d). Hence, at the Cerrado level, these areas did not present a clear trend inside the BZs and slightly tended to an increase in the intra-annual variability of NDVI in the CUs (59.9% of the analyzed cells in the CUs followed this trend) ( Figure E1c).
The non-forest natural vegetation obtained similar results to the natural forest areas for the metric (Figure E1e), suggesting a general homogenization of its surface. This trend is representative of a larger proportion of surface area in the CUs (65.3%) than in the BZs (53.1%). However, the magnitude of change is relatively low ( Figure E1e; Figure E1f). The non-forest natural formations experienced divergent dynamics in the BZs and CUs in terms of productivity change, as shown by the opposed distribution of their NDVI inter values ( Figure A2a). In the BZs, 58.9% of the grassland and wetland areas were characterized by a decrease in the annual average of NDVI, whereas in the CUs, 60.4% of these areas experienced an increase in NDVI. However, in Figure A2b, we can observe a clear disparity of trends depending on the location of these areas. For example, the near totality (90.5%) of this vegetation type in Piauí has experienced a decrease in productivity, while the opposite trend can be observed in Bahia, where the non-forest natural vegetation productivity has increased in 92.8% of the sampled cells.
Similarly, the results for the NDVI intra metric were diverse according to the location of the non-forest natural formation areas ( Figure A2d). Hence, at the Cerrado level, these areas did not present a clear trend inside the BZs and slightly tended to an increase in the intra-annual variability of NDVI in the CUs (59.9% of the analyzed cells in the CUs followed this trend) ( Figure A2c).
The non-forest natural vegetation obtained similar results to the natural forest areas for the NDVI spatial metric (Figure A2e), suggesting a general homogenization of its surface. This trend is representative of a larger proportion of surface area in the CUs (65.3%) than in the BZs (53.1%). However, the magnitude of change is relatively low ( Figure A2e; Figure A2f).  Regarding the intra-annual variability of the pastures (as measured by the , presented in Figure E2c), no clear trend is noticeable except locally, in the CUs and BZs of Bahia, where the intra-annual variability of NDVI decreased in 82.3% of the analyzed pastureland ( Figure E2d).
The spatial variability of NDVI in pastures did not follow any clear trend ( Figure E2e). A slight trend towards a general increase in the spatial variability of NDVI is noted within the pasture areas in the BZs (Figure E2e), which is influenced by the dynamic of the pastures in the states of Bahia, Tocantins, and Mato Grosso do Sul ( Figure E2f). Pastures in the BZs generally experienced a loss in productivity (64.8% of the analyzed cells in the BZs experienced a decrease in the annual average of NDVI) ( Figure A3a), the most affected areas being located in the state of Bahia (where 93% of the pasture cells have a negative NDVI inter value) ( Figure A3b). In the CUs, the number of pasture cells analyzed was limited to 1319 cells (approximately 8244 ha), 57.2% of which have experienced a gain in the annual average of NDVI.
Regarding the intra-annual variability of the pastures (as measured by the NDVI intra , presented in Figure A3c), no clear trend is noticeable except locally, in the CUs and BZs of Bahia, where the intra-annual variability of NDVI decreased in 82.3% of the analyzed pastureland ( Figure A3d).
The spatial variability of NDVI in pastures did not follow any clear trend ( Figure A3e). A slight trend towards a general increase in the spatial variability of NDVI is noted within the pasture areas in the BZs (Figure A3e), which is influenced by the dynamic of the pastures in the states of Bahia, Tocantins, and Mato Grosso do Sul ( Figure A3f). The results of the monitored agricultural areas in the BZs ( Figure E3a) show a general trend of productivity increase (64.6% of the analyzed Agriculture cells in the BZs experienced an increase in the annual average of NDVI, thus a positive value). This trend is especially clear in the states of Goiás and Distrito Federal, where the proportion of agricultural areas with an increased productivity accounted for 82.4% and 84.7%, respectively ( Figure E3b).
However, the agricultural land in the BZs showed the highest dispersion of values ( Figure E3a) of all the four vegetation types, reflecting the rich diversity of cropping systems. This variability in the results is illustrated by the different trends followed by the agricultural land in the different states, and by the high dispersion of values within each state ( Figure E3b). The few agricultural areas monitored within the CUs' limits (41 cells or 256 ha) did not show a clear trend ( Figure E3a).
Just as for , the agricultural land in the BZs exhibited the highest dispersion of values of all vegetation types ( Figure E3c). However, in this case, the high diversity of results at the state level ( Figure E3d), resulted in a lack of a common trend at the biome level. Nevertheless, the agricultural land of the study areas in Distrito Federal and Piauí, exhibited a clear, increased intra-annual variability of NDVI trend (accounting for 88.5% and 79.7% of their total agricultural area, respectively), even if the sampled surface area was relatively small for both states (3918 ha for Distrito Federal and 956.25 ha for Piauí). Conversely, most of the limited sample of agricultural land inside the CUs (73.2%) experienced a decrease in the intra-annual variability of NDVI ( Figure E3c).
As for the other landscape metrics, the agricultural land in the BZs exhibited a greater dispersion of the values than the rest of vegetation cover types ( Figure E3e), and, in this case, a slight The results of the monitored agricultural areas in the BZs ( Figure A4a) show a general trend of productivity increase (64.6% of the analyzed Agriculture cells in the BZs experienced an increase in the annual average of NDVI, thus a positive NDVI inter value). This trend is especially clear in the states of Goiás and Distrito Federal, where the proportion of agricultural areas with an increased productivity accounted for 82.4% and 84.7%, respectively ( Figure A4b).
However, the agricultural land in the BZs showed the highest dispersion of NDVI inter values ( Figure A4a) of all the four vegetation types, reflecting the rich diversity of cropping systems. This variability in the NDVI inter results is illustrated by the different trends followed by the agricultural land in the different states, and by the high dispersion of values within each state ( Figure A4b). The few agricultural areas monitored within the CUs' limits (41 cells or 256 ha) did not show a clear trend ( Figure A4a).
Just as for NDVI inter , the agricultural land in the BZs exhibited the highest dispersion of NDVI intra values of all vegetation types ( Figure A4c). However, in this case, the high diversity of results at the state level ( Figure A4d), resulted in a lack of a common trend at the biome level. Nevertheless, the agricultural land of the study areas in Distrito Federal and Piauí, exhibited a clear, increased intra-annual variability of NDVI trend (accounting for 88.5% and 79.7% of their total agricultural area, respectively), even if the sampled surface area was relatively small for both states (3918 ha for Distrito Federal and 956.25 ha for Piauí). Conversely, most of the limited sample of agricultural land inside the CUs (73.2%) experienced a decrease in the intra-annual variability of NDVI ( Figure A4c).
As for the other landscape metrics, the agricultural land in the BZs exhibited a greater dispersion of the NDVI spatial values than the rest of vegetation cover types ( Figure A4e), and, in this case, a slight trend towards a general homogenization of the agricultural landscapes is observed. This trend represented 58.7% of the monitored agricultural land in the BZs and was representative of the general