Comprehensive Framework for Analysing the Intensity of Land Use and Land Cover Change in Continental Ecuadorian Biosphere Reserves

Abstract

The conventional methods of land use and land cover (LULC) analysis are frequently based on crosstabulation matrices. However, this analysis becomes complex when including sites with multiple management zones and categories at different time points. This is challenging regarding data processing and the presentation of numerous results. We transformed the graphical representation of conventional Intensity Analysis to assess and compare the intensity and magnitude of LULC changes in six Continental Ecuadorian Biosphere Reserves at three levels: interval, category, and transition, and at three time intervals. A dimension reduction strategy was used to convert multiple bar charts into three composite heat maps. The results confirm the global conservation effectiveness in the core zones among the biosphere reserves with less than 10% change, but worrisome dynamics were detected in buffers and transitions with a change of up to 10% for all periods. Deforestation and agriculture were detected as the most relevant land processes. This work highlights the dimension reduction by summarizing 17, 51, and 312 bar charts from conventional Intensity Analysis in three composite heat maps, one for each level of analysis. Systematic suspicious transitions occurred in the water body category because of its dynamics and misclassification in LULC maps.

1. Introduction

Over time, human interactions with nature have been reflected in changes in land use and land cover (LULC). These have impacted ecosystem services [1], biodiversity and climate [2], productivity [3,4], etc., thus affecting the sustainability of land use. Therefore, the interest in promoting natural resource management and ecosystem monitoring at different scales requires a deep understanding and study of LULC dynamics, its spatial distribution, and the identification of patterns and changing processes [5]. In Agenda 2030, LULC plays an important role in SDG 15 due to its interlinkages with other SDGs [6]. Therefore, the detailed study of LULC dynamics is the main basis for supporting several targets.

Landscape transformation occurs when an existing LULC at an initial time point is converted to another LULC at a later time point; this is called a transition [7,8]. LULC transitions are of great interest since they are the territorial expression of human actions or natural processes [9] that affect land conservation or transformation. Anthropic transitions are the result of decisions made by multiple actors and interests [10]; therefore, their control is essential for the integral sustainability of the territory. Understanding transitions in a specific context of global change allows for identifying whether they are due to tendential, constant, or emerging spatial processes [11,12].

The analysis of LULC transitions is essential to assess the magnitude of the change and its effect relative to baseline conditions for policy development and to avoid adverse and undesirable side effects. In addition to LULC transitions, the analysis of LULC persistence is fundamental for conservation purposes or the maintenance of protected areas and for educational, political, economic, and social decisions in the area. Both transitions and persistence have been widely studied and modeled worldwide using a variety of approaches, from simple spatial to agent-based models. However, these are rarely used by decision makers [13] because of the complexity of implementing the methods and tools. This creates a disconnect between science and policy makers.

Overall, conventional methods of LULC analysis are based on crosstabulation matrices and cross-classification maps to quantify gains and losses and to represent spatially explicit changes over a time interval, respectively. In the last two decades, crosstabulation matrix analysis has been extended to three levels: interval, category, and transition, to provide additional information on the intensity and magnitude of change in LULC maps through Intensity Analysis (IA) [14]. However, these methods apply to single site maps, the same set of categories, and for more than two or more time points. Furthermore, the results are presented in individual graphs, limiting comprehensive analysis of the study area.

LULC analysis for integrated territorial management requires the analysis of more than one site and the inclusion of multiple management zones [15], which is useful for decision making, fund investment, prioritization, and monitoring areas of interest. Terrestrial protected areas are an example of a LULC analysis unit where all the above elements converge, hindering the understanding of territorial dynamics and, therefore, the achievement of conservation goals [16]. Furthermore, within these territorial units, several LULC transitions interact, such as agricultural and livestock activities, land abandonment, deforestation, reforestation, and the expansion of roads and urban areas, which in turn increases the complexity of these transition areas. The effective representation and communication of these complex dynamics is currently an open question for land change science.

In this sense, the LULC change analysis of different sites with multiple management zones [17,18] and multiple categories [19,20] at different time points [21,22] is challenging in terms of data processing and the presentation of numerous results [23]. Currently, the representation of both the intensity and the magnitude of the change in conventional IA [14] is through bar graphs, which works and is easily understandable for low-dimension LULC analysis. However, when performing more complex territorial analysis, it is necessary to design a graphical strategy to communicate the results effectively [24], i.e., feature-expression heat maps [25]. This can bridge the gap between knowledge based on complex models [26,27] and decision making [28,29] at different decision levels of land management and stakeholder engagement events.

Our work aims to enhance the advancements in IA presented in [14,30], specifically focusing on the graphical representation of the results for analyzing the LULC evolution within the three management zones (core, buffer, and transition) across six biosphere reserves in continental Ecuador. The study spans intermediate intervals between 1990 and 2018. We used a dimension reduction strategy to transform the graphical representation utilized in both [14,30]. This transformation resulted in a comprehensive framework that summarizes the LULC dynamics in three composite heat maps, one for each analysis level (interval, category, and transition). These composite heat maps include a ranking system for annual change intensity to facilitate the identification of the largest changes at each level of analysis, including LULC change intensity and change magnitude. Furthermore, the composite heat map at the transition level of analysis allows for the identification of specific land processes and potentially suspicious transitions. Our goal is to provide a more insightful and condensed overview of LULC dynamics in biosphere reserves.

2. Materials and Methods

2.1. Study Area

Our study includes the six Continental Ecuador Biosphere Reserves (CEBRs) designated under the intergovernmental UNESCO’s Man and the Biosphere (MAB) Programme: Yasuní, Sumaco, Podocarpus–El Cóndor, Macizo del Cajas, Bosque Seco, and Chocó Andino. These are distributed throughout the country and cover approximately 21% of the territory. Each CEBR includes different ecosystems from the Amazon, Andean, and coastal regions (Figure 1). Complete details about the ecological, socio-economic, and hydroclimatic characteristics of the Ecuadorian biosphere reserves can be found at https://en.unesco.org/biosphere/lac#ecuador (accessed on 20 September 2023). Yasuní was the first CEBR to be designated in 1989. Primary humid tropical rainforest is the dominant ecosystem in this area, which is part of the upper Amazon basin with predominant Ultisols and Inceptisols soil types. This ecosystem shelters a wide diversity of species per square meter and is considered to be the most diverse ornithological site in the world. The population in this biosphere reserve (approximately 20,000 inhabitants) is composed of ethnic groups (Waorani, Kichwa, and Shuar) and people living in voluntary isolation. Coffee and cocoa production destined for agro-industry and export is one of the economic activities in this area. Petroleum activities in the Ishpingo Tambococha–Tiputini (ITT) block within the biosphere reserve generate significant income for the country. The core zone of this biosphere reserve is the same as the boundary of the Yasuní National Park. Sumaco was designated in 2000 and is located near Yasuní and the metropolis of Quito. It protects different ecosystems that have developed around the Sumaco volcano (3900 m.a.s.l.), from the Andean Mountain range to the Amazonian Plain, with Inseptisols and Andisols as representative soil types. This biosphere reserve is located near the metropolis of Quito. The local economy depends on agricultural production. The Sumaco and Napo Galeras National Parks are defined as the core zones of this biosphere reserve.

Podocarpus–El Condor, designated in 2007, is a highly biodiverse area in the south of the country, with an ecosystem gradient and contrasting climatic conditions. This area harbors 25% of the recorded flora in the country and is home to half of Ecuador’s bird species. There are approximately 1,729,000 inhabitants concentrated in the transition zone. The local economy is mainly based on subsistence farming and livestock raising. The Podocarpus and Yacuri National Parks and Cerro Plateado Biological Reserve are the core zones in this area. Macizo del Cajas designated in 2013, is the first biosphere reserve located in the western foothills of the Andes with diverse soil types from Entisols and Inceptisols on the coast to Andisols and Mollisols on the Andean side. Its natural ecosystems (i.e., páramo, cloud forest, montane forest, ocean, and mangrove) play a fundamental role in the provision of water from the Atlantic and Pacific slopes of the Andes Mountains range. The city of Cuenca, with a population of around 640,000 people, is located in the transition zone of this biosphere reserve. Both the El Cajas National Park and the Quimsacocha National Recreation Area are core zones. Although Macizo del Cajas also includes marine areas, this study has only focused on an analysis of the inland areas.

Bosque Seco, designated in 2014, is located within the Bosques de Paz transboundary biosphere reserve. This CEBR is characterized by a seasonally dry forest, which is the habitat of endemic species of vascular plants and birds. The population living in the biosphere reserve is predominantly rural (106,000 inhabitants). The main economic activities include fishing, agriculture, and tourism. This last activity is associated with the massive flowering of the guayacán, which occurs after the first rains at the beginning of each year. The core zone of the Bosque Seco is made up of municipal reserves. Its predominant soil types are Inceptisols, Alfisols, Entisols, and Ultisols. Chocó Andino was the last CEBR to be designated in 2018. It is located northwest of Quito in the Andean region and is dominated by soil types such as Andisols and Mollisols. This biosphere reserve includes two critical ecoregions: the humid moist forest of the Chocó–Darien region, which extends from Panama to western Ecuador, and northern Andean Mountain forests. Both are highly endemic. This area has a population of around 880,000 inhabitants. The main economic activities are retail trade and industrial manufacturing, and the secondary activities are agriculture, fish farming, and livestock raising. The core zone in this area consists of nine protected forests, three conservation and sustainable use areas, several private reserves, and the Pululahua National Park, which is also a geobotanical reserve.

The transformation of land cover in Ecuador due to anthropogenic activities has also affected the CEBRs, modifying land cover patterns and the landscape. The CEBRs are areas of interest since they represent a wide range of ecosystem characteristics, which in turn serve as a sample to assess the country’s LULC dynamics.

2.2. Input Data and Zoning

We used the official dataset of the Ministry of the Environment, Water, and Ecological Transition of Ecuador–MAATE (by its acronym in Spanish), Figure 2. This dataset includes four raster maps of land use and land cover (LULC) categories for 1990 [31], 2000 [32], 2008 [33], and 2018 [34], as well as the official zoning of the CEBRs, as shown in Figure 2a. The maps were derived from Landsat 4–5 TM, 7 ETM, and 8 OLI and from Aster and RapidEye. The satellite images used to generate the maps had a difference of 2 years before and after the reference date to select cloudless satellite images. The overall accuracy of map classifications is 69% (1990), 73% (2000), and 76% (2008) [35], and for 2018, this information is not available. The methods used in terms of satellite image correction, map classification, and validation for the 1990, 2000, and 2008 maps can be consulted in [36].

This study considered the following categories at level 1: forest (F), agricultural land (AL), built-up (BU), shrub–herbaceous vegetation (SH), water bodies (WB), and other (Ot). The spatial resolution of the raster maps is 30 m. All maps have the same categories and spatial extent. The first level of land use categories from the Ministry of Environment Water and Ecological Transition of Ecuador MAATEE is recommended since this level of analysis guarantees the temporal stability of categories. Level 2 in the maps does not preserve the same categories, causing inconsistencies.

The official zoning of the CEBRs corresponds to the core, buffer, and transition zones and has been used in this study, except for Yasuní, as there is no official zoning for this biosphere reserve. However, based on a zoning proposal in 2007 [37], we used the following rules to define Yasuní’s management zones. The core is the union of the Yasuní National Park boundary and the intangible zone of Executive Decree 751 of 21 May 2019. The buffer zone is the symmetrical difference between the boundary of the biosphere reserve of MAATE and the core area. The buffer zone of the Yasuní biosphere reserve agrees with the buffer zone defined in Executive Decree 751 of 21 May 2019, which established a buffer of 10 km around the intangible zone. The transition zone has not yet been officially defined and was not taken into account in this work. Each raster was clipped with the zoning boundaries of each CEBR (Figure 1(1)).

2.3. Change Analysis for Conservation Assessment of Biosphere Reserves

The analysis of change was based on the Intensity Analysis (IA) method proposed by [13], which uses crosstabulation matrices and bar charts to analyze LULC changes between two maps at three levels: interval, category, and transition. At each level, the uniformity value is calculated, which is the cut-off for determining the intensity of the change. IA also considers the magnitude of the change at three analysis levels and works for a single site, the same categories, and for more than two points in time. Our work extends the conventional approach of IA by considering maps of multiple time intervals, sites and zones, which means that these maps correspond to the zoning of each CEBR from 1990 to 2018. We also transformed the bar charts into heat maps and added graphical features to consolidate all IA information into three heat maps, one for each level of analysis. This improvement is described in detail in the next section. In the approach used in [13], the first level of analysis is Interval Level Intensity Analysis (IIA). This involves two equations, the first calculating the annual change intensity for a time interval (S_t), and the second calculating the interval uniformity value (IUV), which is a uniform rate that would exist if the annual changes were uniformly distributed over the entire time extent. If the value of S_t (Equation (1)) exceeds the IUV (Equation (2)), the change is fast; otherwise, it is slow, Figure 2b1.

$$S_t=\frac{\left\{{\sum }_{j=1}^J\left[\left({\sum }_{i=1}^J{C}_{tij}\right)-C_{tjj}\right]\right\}/\left[{\sum }_{j=1}^J\left({\sum }_{i=1}^J{C}_{tij}\right)\right]}{Y_{t+1}-Y_t}·100\%$$

(1)

$$IUV=\frac{{\sum }_{t=1}^{T-1}\left\{{\sum }_{j=1}^J\left[\left({\sum }_{i=1}^J{C}_{tij}\right)-C_{tjj}\right]\right\}/\left[{\sum }_{j=1}^J\left({\sum }_{i=1}^J{C}_{tij}\right)\right]}{Y_T-Y_1}·100\%$$

(2)

where S_t represents the annual intensity of change for time interval [Y_t, Y_t+1]; IUV represents the value of the uniform line for time intensity analysis; J represents the number of categories; I is the index for a category at the initial time point for a particular time interval; j is the index for a category at the final time point for a particular time interval; T represents the number of time points; t is the index for the initial time point of interval [Y_t, Y_t+1], where t ranges from 1 to T − 1; Y_t represents the year at time point t; and C_ti,j represents the number of pixels that transition from category i at time Y_t to category j at time Y_t+1.

This notation is also used in the following category and transition analyses.

The second, Category Level Intensity Analysis (CIA), examines each category to measure the gross losses and gross gains across space. For this study, we only considered the gross gains because of the clarity of interpretation of LULC changes. The CIA compared the observed gain category intensities (Equation (3)) and the Category Uniformity Value (CUV), which is equivalent to the S_t of (Equation (1)) that would exist if the change within each interval were uniformly distributed over the entire spatial extent. An observed gain intensity value greater than the CUV indicates an active gain category or a relatively active change; otherwise, it is dormant (Figure 2b2). Equation (1) gives the uniform intensity for the time interval t for this category level analysis; thus, Equation (1) links the interval level analysis with the category level analysis.

$$G_{tj}=\frac{\left[\left({\sum }_{i=1}^J{C}_{tij}\right)-C_{tjj}\right]/\left(Y_{t+1}-Y_t\right)}{{\sum }_{i=1}^J{C}_{tij}}·100\%$$

(3)

where G_tj is the annual intensity of gross gain of category j for time interval [Y_t, Y_t+1].

The third, Transition Level Intensity Analysis (TIA), analyzes the transition in terms of the gain for a given category (Equation (4)). The transition reaches a uniform transition intensity rate of annual change (Transition Uniformity Value, TUV) (Equation (5)) if it is evenly distributed across the categories available for the transition. A transition intensity rate that exceeds the annual TUV rate is considered to be a targeted transition, while the opposite is considered to be an avoided transition for a given category (Figure 2b3).

$$R_{tin}=\frac{C_{tin}/\left(Y_{t+1}-Y_t\right)}{{\sum }_{j=1}^J{C}_{tij}}·100\%$$

(4)

$${TUV}_{tn}=\frac{\left[\left({\sum }_{i=1}^J{C}_{tin}\right)-C_{tnn}\right]/\left(Y_{t+1}-Y_t\right)}{{\sum }_{j=1}^J\left[\left({\sum }_{i=1}^J{C}_{tij}\right)-C_{tnj}\right]}·100\%$$

(5)

where C_tin represents the area of transition from i to n during [Y_t, Y_t+1]; R_tin represents the annual intensity of transition from category I to category n during time interval [Y_t, Y_t+1], where i ≠ n; TUV_tn represents the value of uniform intensity of transition to category n from all non-n categories at time Y_t during time interval [Y_t, Y_t+1].

Detailed calculations for IIA, CIA, and TIA can be consulted in [14]. The R package is available in https://github.com/dondealban/learning-intensity-analysis (accessed on 15 July 2023) and also available in a Microsoft Excel Macro.

The graphical representation of the three-level conventional approach to IA is based on bars (Figure 2b1–b3). The left bars represent the observed magnitudes of the annual transitions, and the right bars represent the intensities of these transitions. Uniformity is a vertical red dotted line located on the right side. The interpretation of these intensities when a bar crosses or does not cross this line varies from level to level, as described above. At all levels of analysis, the observed sizes of the annual transitions and intensities of the changes were compared ([14], Figure 3, Figure 4 and Figure 5).

2.4. IA from Bars Charts to Composite Heat Maps

Using the IA uniformity values calculated at the three analysis levels, we reduced the graphical dimension of conventional IA, i.e., the number of figures. This consisted of calculating the difference (d) between the highest value of each bar and IUV, CUV, and TUV, respectively. Using the calculated d values, we generated heat maps to compare different sites and zoning and their temporal evolution in three summary charts, one for each analysis level (Figure 2c1–c3).

In the Supplementary Material, the calculations and the transformation of the bar charts into the composite heat map are exemplified step by step.

In the heat maps, larger d values can be identified from the high saturation color ramp and smaller d values with the low saturation color ramp. High saturation color represents a fast time interval change, an active gain category, and a target transition in IIA, CIA, and TIA, respectively. On the other hand, low saturation color represents the opposite.

Additionally, we rank the uniformity values of each analysis level (the last column of each heat map) from high to low uniformity values and color from red to green to compare the magnitude of the change by time interval and zoning. High IUV means a higher uniform rate of change comparing all time periods, and high CUV means a higher uniform rate of change in the categories at different time intervals. For IUV and CUV, the color and rank are globally associated, but TUV shows rank and color at each time interval based on the transition from a particular category C_i to another C_j. TUV rank value is calculated from a particular C_i to a different C_j for the three zones and for each time interval. TUV color ramps were applied to each of the ranks formed. Since TUV rank values are defined from high to low uniformity values for each category, there is no correspondence of the change magnitude between C_n+i ranks (i.e., rank 1 of C_n compared to rank 1 of C_n+1). High-rank values for IUV, CUV, and TUV saturate high in the color ramp. Low-rank values for IUV, CUV, and TUV saturate in the opposite direction in the color ramp.

Transition magnitudes were also plotted over the heat maps using pie chart symbols (Figure 2c). These summarize the left side of the bar charts in Figure 2b. The filling of the pie charts is proportional to the annual transition sizes and is expressed according to each level of analysis. The coupled inspection of change intensity (heatmaps) and change size (pie charts) allows the user to easily follow the most striking changes for each level of analysis: (i) between time intervals in IIA, (ii) category changes within a time interval in CIA, and (iii) transitions by zoning and time intervals in TIA. In addition, we represented land processes in the pie chart in TIA by adding a color code, as shown at the bottom of Figure 2c3. The colors indicate land processes such as agricultural and livestock rotation, deforestation, reforestation, urbanization, and transport infrastructure. Within these processes, we also included the category of suspicious transitions for transitions that we suspect were errors in the classification of the maps rather than real change on the ground.

Our approach included all possible categories in the map and did not consider any aggregation strategy to simplify the analysis. Another aspect that allows a better understanding of the influence of dormant categories in larger change sizes can be easily achieved using the heat map approach, where all information can be displayed: intensity change, change size, uniform value, and rank. This allows the analysis of land management zones and transition dynamics to be improved and the level of detail to be increased compared to the conventional IA approach.

Considering the rank and color coding of IUV, CUV, and TUV is the first step in identifying major changes. It is then important to evaluate the size and the intensity of the change simultaneously. This analysis sequence allows us to assess the real change to be assessed at all analysis levels.

3. Results

3.1. Interval Intensity Analysis–IIA

Figure 3 shows a high saturation in the IUV column, mainly in the buffer and transition zones for all CEBRs. The buffer and transition zones of Chocó Andino rank first, with a high uniformity rate of change in all periods. The transition zones of Bosque Seco and Macizo del Cajas follow the high IUV rank. Conversely, the core zones show low IUV values in all time intervals.

In the first time interval, the intensity of LULC change was relatively fast in the transition zones of Chocó Andino and Sumaco, the buffer zone of Macizo del Cajas, and the core zone of Yasuní. The change intensity in the Chocó Andino transition zone was higher during this time interval. In the buffer and transition zones of all CEBRs, except Yasuní, the magnitude of change is greater than 5%. The core zone of Yasuní had a small magnitude of change, close to 0%.

The second time interval showed a fast rate of change intensity in the core zone of Bosque Seco and the buffer zone of Sumaco. In the other CEBR zones, the intensity of change was relatively slow. Overall, the annual transition magnitude was between 5% and 10%. In most cases, a fast rate of change was observed during the last time interval. The annual transition magnitude varies between 5% and 15%. None of the CEBR zones showed a uniform change intensity over the three time intervals.

Intuitively, a user can easily identify core zones with a greener IUV, indicating a low land change impact, as expected in biosphere reserve core zones across all time intervals (Figure 3). A more dynamic state of land change was detected in buffer and transition zones along the three time intervals. The change magnitudes in the core zones were smaller than in the buffer and transition zones, although fast change intensities were registered in the core zones of Yasuní (1990–2000), Macizo del Cajas, and Sumaco (2008–2018). Slow intensities with small change sizes were also identified in the second time interval for the same biosphere reserves in the core zone. The most striking findings were the IUV rank, the intensity and change in size in the buffer and transition zones of Chocó Andino and Bosque Seco over all time intervals, and the size of change in the core zone of Bosque Seco (Figure 3).

3.2. Category Intensity Analysis–CIA

The CUV rank in Figure 4 shows higher values, mainly in the buffer and transition zones in all CEBRs along the three time intervals. Chocó Andino ranks first in these zones from 1990 to 2018, followed by the transition zones of Bosque Seco in the first and third time intervals and the Podocarpus–El Cóndor’s transition zone in the second time interval. Macizo del Cajas ranks third from 1990 to 2000. The core zones of Bosque Seco and Chocó Andino have medium CUV values.

The CUV rank indicates a relative uniformity pattern in the core zones of Yasuní (16,17,17), Macizo del Cajas (13,13,12), Sumaco (15,15,15), Bosque Seco (9,7,9), Podocarpus–El Condor (12,12,13), and Chocó Andino (10,11,10), showing a global constant rank position from 1990 to 2018. This means that the pattern of change in the core zone is constant according to CUV. A similar finding is registered for the Yasuní buffer zone (17,16,16) and the Sumaco transition zone (14,14,14).

The intensity gains of all categories were mainly active during the third time interval rather than during the first and second time intervals in all CEBRs zones. The categories AL, BU, and Ot show an active change in the three time intervals. The SH category has a mostly dormant intensity of change, and the WB category generally has no transitions, with some exceptions in the transition zone in the first and second periods and in the core zone of the second period. The opposite is true for the third period, where the SH and WB categories show an active intensity of change in almost all cases. None of the categories show a uniform intensity in terms of change.

Overall, the AL category was the most active across all time intervals and zones with larger magnitudes of change, and the BU and Ot categories were the most active categories with the smallest magnitudes of change. The F category was mostly dormant, with change sizes greater than 5%, except for the core zone in Macizo del Cajas.

3.3. Transition Intensity Analysis—TIA

Across the three time intervals, in the AL gains, the TUV values are mainly high in the buffer and transition zones in all CEBRs. AL targets SH with different magnitudes of change and Ot mainly with the smallest magnitudes of change. AL mostly avoids F with a magnitude of change greater than 20%, except in the core zone of Macizo del Cajas, where AL targets F with a magnitude of change between 5% and 10% in the first time interval and greater than 20% from 2000 to 2018. AL also targets WB in the last time interval, which was considered a suspicious transition.

In the SH gains, the high TUV values are distributed over all time intervals and CEBR zones, highlighting the core zones of Macizo del Cajas from 2000 to 2018. The intensity and magnitude of change reach high values, especially when SH targets AL. SH avoids F with a tendency to a change in magnitude, mostly less than 5%. SH targets WB in the last time interval, which was again considered a suspicious transition.

WB category gains target other LULC categories in the last time interval, but these were mainly suspicious transitions.

BU’s gains have high TUV in the transition and buffer zones in all time intervals. Chocó Andino ranks first in the transition zone at all time intervals. BU targets AL intensively with a change magnitude greater than 20% in most cases, followed by SH and Ot. BU avoids F, except in the core zones of Podocarpus–El Cóndor from 1990 to 2018. BU targets F in the transition zone of Bosque Seco from 2008 to 2018 with a magnitude of change greater than 20%. Overall, BU avoids WB from 1990 to 2018, but in the last time interval, BU also targets WB in the buffer and transition zones with a change in magnitude of less than 15%. This was considered a suspicious transition.

Ot gains show the high TUV distributed over all time intervals and zones of the CEBRs, even in the core zones, which are protected areas. Ot targets AL and SH with different magnitudes of change. Ot avoids F, except in the buffer zones of Bosque Seco and Chocó Andino and in the transition zones of Sumaco (1990–2000) and Bosque Seco (2008–2018). The Ot category also targets WB from 2008 to 2018, which was considered a suspicious transition.

F gains also show a high TUV distributed across all time intervals and zones of the CEBRs. F targets AL mostly in the core zones and some buffer zones with a magnitude of change greater than 20%. In the transition zones, F avoids AL, except in the Sumaco transition zone from 1990 to 2000.

Crop and livestock rotation occurs when AL shifts to SH or vice versa, SH shifts to Ot or vice versa, and Ot changes to AL. This process is common to all time intervals and zones of the CEBRs and is associated with a common agricultural practice, the fallow period. It does not occur exclusively in the transition and buffer zones, as might be expected, but also in the core zones. The magnitude and intensity of change associated with this process are mostly high and targeted. The crop and livestock rotation is a striking process identified by TIA.

Deforestation occurs when F shifts to AL, SH, BU, or Ot. The intensity of change related to this process is mostly avoided, and the magnitude of the change is relatively high, especially when AL gains from F. There are some exceptions in the AL gains from F, such as the core, buffer, and transition zones of Macizo del Cajas, which show a targeting intensity of change in the three time intervals, with a high magnitude of change. Another exception is in the Ot gains, where the buffer zones of Bosque Seco and Chocó Andino had a targeting intensity of change with the magnitudes of change greater than 20% from 1990 to 2000; this was also observed in the transition zone of Bosque Seco from 2008 to 2018. The deforestation process is also a prominent process in the CEBRs of TIA.

Changes in urbanization and transport infrastructure occur when BU shifts to AL, SH, or Ot. The intensity of change in this process is targeted, especially when BU prefers AL. The magnitude of change associated with urbanization is mostly greater than 5%. This process also occurs in the core zones in all time intervals and zoning of the CEBRs.

The reforestation process occurs when F shifts to AL, SH, or Ot. This process has developed in all zones and time intervals with a magnitude of change greater than 20%.

Some transitions are presented as inconsistent, which is considered a suspicious transition by conventional IA. In Figure 5, these cases have been identified with a pink color for the pie charts, suggesting that this may be part of the classification error of the maps or the effect of climate change because of the dynamics between dry and wet years. The time interval 2008–2018 shows more suspicious transitions, especially in the WB category. This category shows high shifting dynamics between the categories.

Although Figure 5 summarizes 312 bar charts from conventional IA, we have described the most striking transitions. The former order of presenting the results serves as an example to show and discuss all the transitions associated with each land process. In addition, the suspicious transitions can also be traced in terms of intensity and magnitude of change. Specific examples, step-by-step of IIA, CIA, and TIA, are shown in the Supplementary Material.

4. Discussion

4.1. Interpretation of LULC Dynamic Changes in CEBRs

The CEBRs are living scenarios that exhibit (i) persistence and (ii) dynamic shifts in LULC. At the first analysis level, Figure 3 shows the IUV rank, which allows for the easy detection of major interval change intensities to prioritize future actions. A high persistence in the core zones is reflected by lower IUV, suggesting an effective conservation in these areas. Furthermore, the change in magnitude in the Bosque Seco core zone from the 1990s to 2018 is worrying. Other works report similar findings [38,39] but do not compare these across management zones, nor do they use similar management area comparisons with other biosphere reserves.

The year of designation of CEBRs is also a fact that affects the interpretation of LULC changes. This is the case of Chocó Andino, which ranks first in the IUV and has large magnitudes of change from 2000 to 2018 and fast interval change. However, Choco Andino was officially designated as a biosphere reserve in 2018. This biosphere reserve needs attention in the next UNESCO periodic review (2028) to monitor LULC changes. This finding reinforces the idea of reconciling land use conflicts in Choco Andino and supports immediate action, as the situation may be problematic. On the other hand, a slowdown in the intensity of interval change in the coming years could also serve as an example of enhanced land management measures that could be useful for extrapolation to other biosphere reserves. In line with this, the work of Wiegant [40] on the planning of land use realities of the forest and landscape restoration (FLR) objectives of the Bonn Challenge in Choco Andino and Bosque Seco could benefit from evaluating the effectiveness of local measures on land use relevant processes through IIA using IUV as a comparative metric.

In Figure 3, the Bosque Seco transition zone is also a notable site that ranks 3rd in IUV and shows a large size and fast LULC change from 2000 to 2008. This was driven by the expansion of extensive cattle rearing and maize production. The changes occurred in the wet forest remnants between 1000 and 2300 m.a.s.l [41]. This points to the prevention of similar conditions, especially in mountainous landscapes in other biosphere reserves such as Macizo del Cajas-ranked 4th, Podocarpus–El Cóndor-ranked 5th, and Sumaco –ranked 8th.

Composite heat mapping at the IIA level can contribute to LULC modeling during the early stages by displaying overall LULC results in participatory processes. A recent study in northern Ecuador by [42] further discussed the importance of understanding the complexity of LULC change to support LULC modeling and stakeholder engagement in sustainability actions at the local level. In addition, the study highlights the mismatch between national policies and local land use realities. This is an example of how both approaches, IIA composite heat mapping and LULC modeling, could improve land use policy making.

In the second level of analysis (Figure 4), the CUV rank confirms intensive changes among LULC categories in buffer and transition zones. Low overall CUV values indicate persistence in the core zone. In addition, lower values were also found in the Yasuní buffer zone, indicating a low intensive change between categories. This is partly because of the large persistence of the forest category and the smaller changes compared to larger categories, such as agriculture and forests, according to the uniform line of IA. This does not mean that intensive changes in these categories are neglected in LULC change analysis; on the contrary, further analysis with TIA is needed for these areas. One notable shift at the CIA level is between forest and agriculture, and some authors suggest that this shift is not only driven by a deforestation process. For example, in the Sumaco transition zone, part of the shift is due to the expansion of cocoa production using traditional agroforestry systems called chakra. This land process is driven by a change in area of between 0.5 and 4 ha [43]. Despite the impact of this land change in the area, other authors suggest that this process is an essential element for sustainable transition [44] compared to the monoculture agricultural system. On the other hand, the Yasuni buffer zone had a consistent CUV during all of the time intervals because of low anthropogenic intervention. These areas are managed by local communities [45], and low land transformation was confirmed from 1990 to 2018. Inhabitants are committed to the protection of the area since the main livelihoods are dominated by subsistence production and the area is strongly dependent on forests for developing socio-economic activities.

As shown in Figure 4, the CIA increases the level of detail in the analysis by introducing the LULC categories. Most of the changes occurred in the AL, BU, and Ot categories. The CUV in CIA and IUV in IIA for the core zones are similar, corroborating the protection of the core zone. Other approaches [46] also summarize the dynamics of LULC among categories and the issues of current land development patterns in fragile areas in the Three Gorges Reservoir area in China, which is especially associated with reforestation and urban sprawl. This work used the heat map over the transition matrix in percentage and used the concept of the dynamic degree of LULC to quantify the rates of change in different LULC categories within a given time span. However, a limitation in this work concerned the idea of a graphical abridged product being used to analyze LULC dynamics in the time span that includes all categories, as we propose with the composite heat mapping of IA.

LULC changes are the territorial expression of complex interactions that are nontransferable on a spatial or temporal scale between case studies because of the unique properties of territorial development. A challenging objective in monitoring LULC change in biosphere reserves is to provide evidence of financial support to ensure the actual implementation of various policies. Interval (IIA) and category (CIA) assessments will support the implementation of some general policies. However, land processes and ecosystem states will respond to a different set of specific policies, where more detailed information is required to address these policies with greater certainty. In Figure 5, the TIA provides insights into whether a LULC category avoids or targets another category, taking into account the interval and zoning at the same time, allowing comparison across all biosphere reserves. Although this figure reflects the evolution of current land use policies, it does not provide information on the degree of implementation of LULC policies. Nevertheless, the TUV rank between intervals can be compared to assess the progress of policy implementation and its evolution over time.

In Figure 5, one of the most striking and useful results is the color coding of the pie charts to identify land use processes and dynamics. However, in terms of LULC conservation, no transitions are desirable, as shown in light blue cells in the heat map. According to each transition column in Figure 5, WB is the category with the highest number of suspicious transitions, which, in most cases, is due to water logging shifting to vegetated areas, such as humid forests, mangroves, and wetlands during the image classification process. This classification problem has also been reported in other cases in China [47], Mexico [48], and Ecuador [49]. This can be attributed to the misclassification of shadows in images, leading to confusion with other land cover [50]. The WB category in the 2018 LULC map shows potential classification errors because of the suspicious transitions between water bodies and other LULC categories. This can be seen in Figure 5 in the magenta pie charts. Decisions on these transitions need further investigation, especially when WB targets AL and SH with large change magnitudes. In part, this can be attributed to basin irrigation systems and water logging due to heavy rains recorded along the Ecuadorian coast in 2018 [51]. The dynamics over Andean paramo lakes also represent shifts between vegetation and water bodies, which can be attributed to the misclassification of imagery in 2018. River sinuosity also affects inundation areas that shift with vegetation, and meandric river dynamics imply a potential classification error between vegetation and water bodies. Two other classification errors may occur in the coastal environment between water bodies and vegetation, corresponding to the aquaculture of shrimp farms and mangroves. Some authors dealing with these special cases in aquaculture [52] and mangroves [53,54] suggest the careful selection of pure training pixels when defining the classifiers. In the last period of analysis, WB transitions were systematically misclassified, indicating that the classifier or the image source was changed during the implementation of the image classification protocol. It is advisable to users implementing this approach that any interpretation of suspicious transition needs to be deeply assessed, especially for policy making.

In terms of conservation, if F targets AL with large change magnitudes in core areas, this denotes an important result with respect to the effectiveness of reforestation efforts. Similarly, when AL avoids F with high change magnitudes, these are considered prominent cases in almost all CEBR management zones because the change intensity is so low that F is very close to becoming a target of the AL category, implying exacerbated deforestation rates. Another dynamic LULC shift is from SH to AL and vice versa, denoting an active rotation between cropping and livestock rearing. These categories may change intermittently as livestock and cropping are parallel activities of land use in which the rotation may also serve as a soil conservation measure but also as an intensive one. There is no evidence in the literature about the impact of this rotation in Ecuador, but a dynamic between the years indicates an active use of land for cropping and livestock rearing. The effect of these LULC changes in terms of dynamics and intensity warrants further research. Figure 5 is useful for comparison between interval and zoning using TUV. Persistence or low-rank change in TUV indicates continuity in land change or stable management for each land category.

This is the first time that all the information transitions derived from IA are displayed to evaluate the LULC changes, taking into consideration the intensity, magnitude of change, no transitions, suspicious transitions, and land dynamics between different intervals and zoning. The use of heat maps serves to comprehensively display all the information in a condensed form to obtain a ranking to easily track relevant changes compared to changes in other intervals and zones.

4.2. Application of Dimension Reduction and IA in LULC Analysis

Conventional IA and its subsequent extensions have been widely used around the world to study the LULC dynamics of natural and urban landscapes in the context of climate change [55,56], biodiversity [57,58], deforestation [59,60,61,62], urban growth [63,64,65,66,67], environmental risks [68,69], ecosystems [70,71,72], error map quantification [73,74,75], and gender land use preferences [76], among others. These studies share the following common factors in their analyses: a single study area, a set of maps allowing analysis of up to four time intervals, and a set of categories associated with the map layer. Even areas designated as biosphere reserves, where there is zoning, are analyzed in the same way [77]. The results obtained in these studies report the bar charts associated with each level of analysis, which are manageable and understandable when the factors mentioned above are considered. However, if, for example, the zoning of biosphere reserves were included in the analyses, the number of bar graphs would increase, and so would the difficulty of managing and understanding the results.

We propose not only integrating information on the magnitude and intensity of the change in the graph of the third level of analysis, as suggested by [30], but also adding this information to the three levels of analysis. Our approach reveals the magnitudes and intensities of the changes in three heat maps and a classification of the uniform annual rate of change. In addition, the composite heat map of the third level of analysis identifies land processes and suspicious transitions.

This approach allows the inclusion of multiple scales: case studies, zoning, categories, and time intervals. This is useful for decision making in different contexts and management levels. This approach can even be used to report significant LULC changes in the national periodic review report of Biosphere Reserves required by UNESCO every ten years.

4.3. Advantages and Shortcomings

The reporting of LULC change analysis and its results is nowadays an important practice in terms of decision making. The method developed in this work is presented to overcome the problem of summarizing complex sets of IA results. Composite heat mapping is an intuitive tool for displaying multiple results to highlight and easily track changes. Comparisons between zoning areas are also of interest in evaluating evolution or persistence in protected areas. This method offers a straightforward framework to assess multiple aspects simultaneously using graphical summarization when reporting LULC change using the IA method. Four relevant advantages are presented in composite heat maps: (i) the ranking and color coding of uniformity values to quickly identify higher and lower annual changes, (ii) the intensity of change in land transformation through the heat map color ramp, (iii) the magnitude of the change presented in pie charts, and (iv) the color coding included in the pie charts to identify land processes. The intention of the framework proposed does not imply validation of data sets or official maps since the approach is a complementary analysis that highlights suspicious transition detections.

Most policy makers and decision makers are not interested in the processes behind models or methods, but they are interested in concrete results from complex analyses. On the other hand, researchers are constantly confronted with decisions and policy making with a common barrier of balancing between sophisticated methods to overcome complex issues and understanding the knowledge levels of all stakeholders. To overcome this barrier and effectively communicate LULC dynamics, it is necessary to summarize, prioritize, and simplify the results. The framework developed in this work provides the possibility to convey complex LULC dynamics in an inclusive manner, as this framework is a tool that will support dialogue and agreement between stakeholders related to deforestation, forest conservation, urban sprawl, and agricultural expansion.

However, the balance between complex LULC change, decision making, and policy making needs to be handled carefully, as simplification also carries social risks. In order to move beyond risky decisions and to ensure that the results are understood by non-experts, specific assessments need to be made to match the scientific results. This type of assessment, which was not developed in this work, needs to be addressed in the next steps. Advantageously, the approach presented in this work to promote IA in the study of LULC change could be further developed to improve land change science through feedback from participatory stakeholder workshops.

Several works propose different landscape indexes considering territorial processes such as urban sprawl [78,79,80,81,82,83,84,85], cropland expansion [86,87], or deforestation [88,89,90,91,92], in which expansion is usually measured via the buffering of core patches or cumulative maps. The index approach is compatible with composite heat maps with the common aim of characterizing territorial dynamics to quantify multi-temporal changes. Despite the index approach being based on geo-referenced data, it requires technical stakeholders to understand the dynamics of individual patches. Our approach is useful in terms of usability since it summarizes all patches within territorial zoning at three levels of analysis, increasing the detail in the multi-temporal understanding of changes by non-technical stakeholders. Thus, indexing and composite heat maps support various stakeholders to enhance landscape planning issues.

One of the most highlighted challenges is effectively contributing to SDG 15, which among diverse land processes, focuses on reversing land degradation. To halt the land process towards ecosystem degradation, target 15.b, which addresses financial resources for conservation effectiveness, requires comprehensive methods to facilitate financial support assessment using indirect methods as geospatial alternatives. This target is associated with forest management, which is central to most biosphere reserves. The changes in dynamics at three levels of analysis could benefit our understanding of the effectiveness of conservation measures to evaluate the degree of conservation status in temporal and spatial scales through the rankings of IUV, CUV, and TUV. Comparing all core areas of biosphere reserves at intervals, categories, and transitions allows us to confirm that land degradation processes in these areas are scarce, and management could be positively evaluated in the achievement of target 15.1 from a mapping analysis perspective. However, the social science analysis conducted by Krauss in 2022 [93] suggests that despite the SDGs being indivisible, some aspects, such as social and economic matters, are dealt with by other SDGs, and there is no need to address them in SDG 15: Life on Land. Thus, the composite heat mapping of IA is presented as being a valuable tool contributing to SDG 15, but the effect on social and economic interlinkages needs to be cautiously addressed since contextual data are required to be adequately treated in relation to other SDGs.

The Ecuadorian biosphere reserves were used in this work to generate composite heat maps of IA, but this approach is suitable for any case study in which the analysis includes several time intervals, zoning, and multiple categories. In the assessment of ecosystem change, it is recommended to establish common border rules to avoid ecosystem delimitation discrepancies. Image classification with a high accuracy of approximately 70% [94] is advisable when analyzing changes in sensitive ecosystems or highly detailed changes. Biosphere reserve maps derived from supervised classification at 30 m resolution are enough to perform the three levels of analysis, as most areas of the vegetation structural classes occurred at scales ≥ 30 m [95].

The comprehensive framework can support different fields of study related to LULC change to ensure a good understanding of land use changes at different user levels. In future research, various calculations may also be possible, such as composite heat maps to retrieve LULC indices, land process dynamics, and sensitivity analysis of suspicious transitions.

5. Conclusions

An important aspect of LULC change analysis is the effectiveness of the representation and interpretation of results in dynamic areas such as biosphere reserves. The location of specific land processes is also a prime concern for monitoring actions in the short term. Factors relating to time intervals, zoning, and administrative boundaries are needed to propose and design integrated remedial actions. However, the inclusion of all these factors increases the complexity of representing, understanding, and effectively communicating LULC dynamics to decision makers and stakeholders. The implementation of both aggregated and straightforward alternatives is needed to configure effective decisions once this complex LULC change analysis, such as the IA, has been carried out. The framework developed in this work is presented as an alternative to improve the understanding of complex territorial processes, specifically intended to analyze contrasting objectives or cases as the persistence and the dynamic shifts in a comprehensive manner.

The following aspects note the advantages of using composite heat maps for IA:

• Our work offers an alternative framework for IA to visually identify and rank LULC dynamics in three composite heat maps, one for each level of analysis: interval, category, and transition.
• The composite heat maps were created based on factors that are commonly considered in the LULC change analysis for decision making. These are multiple areas of interest, zoning, more than three map layer categories, and different time intervals.
• Each composite heat map integrates information derived from the IA, such as the uniform annual rate classification, the magnitude, and intensity of the change, and at the final level of analysis, it is possible to identify land use dynamics and suspicious transitions of maps through color coding.
• The simultaneous evaluation of the magnitude and intensity of the change allows an integrated assessment of LULC change.
• The ranking of uniformity values and color coding is used to identify and prioritize uniform intensity changes at each level of analysis.
• The composite heat maps provide evidence of the conservation effectiveness of core zones in all CEBRs. In addition, it warns of LULC changes in buffer and transition zones.