Spatial near future modeling of land use and land cover changes in the temperate forests of Mexico

Study area

The study area is located in the western part of the state of Chihuahua, Mexico. It is part of the ‘Sierra Tarahumara’ and have a surface area of 497,159 ha. Its extreme coordinates are 108°00′W, 29°00′N and 107°10′W, 27°30′N (Fig. 1). It is one of the regions of temperate forests, which has experimented the greatest disturbances in the past years in the state of Chihuahua (Herrera, 2002). It belongs to the most extensive forest areas in North America. It is immersed within a complex orography composed of large canyons and deep canyons, which results in a mixture of temperate and tropical ecosystems. It is characterized by its high biodiversity and number of endemic species, estimating the presence of around 4,000 species of plants. Also, it is recognized by the International Union for the Conservation of Nature as one of the megacenters of plant diversity (Felger & Wilson, 1995). The main land uses in the area include: pine forests, oak forests, pine-oak and oak-pine forest associations, agriculture and grassland communities. The economic activities in the region are forestry, extensive livestock and rainfed agriculture (INEGI, 2003).

Data source

For the analysis of the LULCC, three scenes of the Landsat sensor (Path 33, Row 41), with a spatial resolution of 30 m, were used. The scenes corresponded to the years 1990, 2005 and 2017 and they were acquired from clear sky days and each of them taken during the same month to reduce the temporal variation. The scenes were downloaded from the United States Geological Survey (USGS, 2018). The characteristics of each scene can be seen in Table 1.

The scenes were radiometrically corrected. The radiometric correction was carried out with the QGis software v.2.8 through the SemiAutomatic Classification plugin (Congedo, 2017).

Integration and composition of bands

Once the scenes were corrected, they were integrated into a layer stack. False color composites for the Landsat TM5 were then generated, with a combination of the bands 5, 4 and 3. Band 5 corresponds to the infrared channel (1.55–1.75 µm), band 4 to the near infrared (0.76–0.90 µm) and the band 3 to the red channel (0.63–0.69 µm). This combination was applied to the scenes of 1990 and 2005. Regarding the scene of 2017, the combination for Landsat OLI8 was applied and corresponded to the bands 6, 5 and 4, where band 6 corresponds to the medium infrared channel (1.55–1.65 µm), band 5 to the near infrared channel (0.85–0.88 µm) and band 4 to the red channel (0.64–0.67 mn) (Lillesand & Kiefer, 2000).

Land use and land cover classification

The Suport Vector Machine (SVM) classification was applied to the 1990, 2005 and 2017 images through the software R (R Core Team, 2016) with the R package “caret” (Kuhn, Wing & Weston, 2015) to obtain LULC information. The SVM classifier is a supervised technique of nonparametric statistical methods (Mountrakis, Im & Ogole, 2011). The SVM classification has been used in several research studies in the past (Kavzoglu & Colkesen, 2009; Otukei & Blaschke, 2010; Shao & Lunetta, 2012). For the supervised classification, five classes of land use were defined; (1) primary forest, (2) secondary forest, (3) human settlements, (4) areas without vegetation and (5) water bodies (Table 2).

Modeling and spatial simulation with Dinamica-EGO

The LULCC scenarios were made based on the historical trends of change in forest cover during 1990–2017 of the supervised classifications using Dinamica-EGO (Soares-Filho, Cerqueira & Pennachin, 2002). The historical trends of LULCC is based on the transition matrix (Monteiro Junior, Silva E De Amorim Reis & Mesquita Souza Santos, 2018). Dinamica-EGO uses the algorithm of cellular automata, and the method Weights of Evidence (Olmedo et al., 2018). For the simulation of deforestation, the following steps were undertaken: (1) selection of change drivers as well as transitions, (2) exploratory analysis of the drivers of deforestation, (3) simulation and (4) validation. These four steps are described in the following sections.

Selection of variables and transitions

The selection of the set of exploratory variables to simulate the LULCC is essential for the modeling success (Miranda-Aragón et al., 2012; Pérez-Vega, Rocha Álvarez & Regil García, 2016). In this study, 19 variables were used; 17 static and two dynamic variables. Static variables remain constant during model execution. Dynamic variables change during the execution of the model and they are continuously updated in each iteration (Olmedo et al., 2018). The set of variables used is shown in Table 3.

The transition refers to the total amount of LULCC that occurred in the simulation period. In this study, the transitions of interest were: (a) primary forest to secondary forest, (b) primary forest to areas without apparent vegetation, (c) primary forest to urban areas and (d) secondary forest to areas without apparent vegetation (Table 4).

Exploratory analysis of the data

When we modeled LULCC dynamics, Weights of Evidence (WoE) were applied to project transition probabilities. Regarding deforestation, degradation or any other type of change, we previously know about the location of favorable conditions for LULCC. The influence of static and dynamic variables and the elaboration of the LULC maps was performed with WoE in the Dinamica-EGO software (Soares-Filho et al., 2010).

The positive values of WoE represent an attraction between a transition of land use and a specific variable. The greater the value of W⁺, the greater the probability of transition. Negative values of W⁻ indicate low probabilities of transition instead (Maeda et al., 2010a; Maeda et al., 2010b). By using the WoE values of the variables used in the analysis of LULCC, the Dinamica-EGO model calculates the transition probability of each pixel to change. Thus, the pixels are assigned with a probability value for a given transition and probability maps are generated for the transitions of interest (Soares-Filho, Rodrigues & Costa, 2009; Soares-Filho et al., 2010; Mas & Flamenco, 2011).

Given that the basic hypothesis of the WoE technique is that the driving variables must be independent, for this study the correlation between the variables was tested through the Cramer Coefficient (V), represented by Eq. (1). (1)${\displaystyle V=\sqrt{\frac{{\chi }^2}{\mathit{\Gamma }\dots M}}}$

Where: χ² = is the chi-square statistic of the contingency between two variables, Γ = denotes the sum of the values of contingency, M = is the minimum of n − 1 or m − 1, where n denotes the number of rows and m the number of columns. Bonham-Carter (1994) mentioned that values lower than 0.5 for the Cramer Coefficient (V) suggest independence, while values higher than 0.5 involve a greater association (Almeida et al., 2003; Teixeira et al., 2009).

Simulation of land use and land cover changes

Three types of scenarios were used for 2050; they were called pessimistic, optimistic and stationary. For the three scenarios, the modeling base was the period 1990–2017. The transition matrix of 1990 and 2017 were used to estimate the possible change in forestry coverage in the future, taking 2017 as the beginning year and 2050 as the final year. In the pessimistic scenario, the transition probability matrix and the change function (patcher and expander) were modified, increasing the deforestation and fragmentation rates between 1990 and 2017. This was done based on the hypothesis that the development of road infrastructure, urban expansion, fires, uncontrolled exploitation, among others, will produce strong spatial changes of land use. For the optimistic scenario, the state and national forest development plans were considered. Such plans promote the protection and conservation of forest resources (CONAFOR, 2001). For this scenario, the conservation and promotion of strategies to protect forests were represented by reducing the transition matrix value, as well as the patcher and expander change functions. Regarding the stationary scenario, transitions or change functions were not modified. In this case, it is assumed that the trend will be the same as the one between 1990 and 2017.

Validation

To evaluate the model performance, we used a Fuzzy Similarity Index (FSI), where the representation of a pixel is influenced by itself and its neighborhood (Ximenes et al., 2011; Yanai et al., 2012; Chadid et al., 2015). The FSI employed in this study was developed by Hagen (2003), modified by Soares-Filho et al. (2017) and implemented in Dinamica-EGO. The FSI verifies the agreement between the observed and the simulated land use and land cover datasets by obtaining the number of coincident cells within increasing window sizes of a neighborhood (Costanza, 1989; Soares-Filho et al., 2017). The validation process was carried out by comparing a simulated map and a reference map. The simulation of the 2017 LULCC map was generated. To generate the simulation of 2017, the transition matrix was used between 1990 and 2005. The comparison through the FSI allowed to evaluate the areas of coincidence of change and no change between the real and simulated map of 2017. Finally, the general procedure used in this study is outlined in the flowchart depicted in Fig. 2.

Detection of land use/land cover changes

Results from the analysis of LULCC show a considerable gain for secondary forest. The forest cover of the primary forest was reduced from 55.8% of the study area in 1990 to 37.7% in 2017. The areas without vegetation increased their area from 4.11% to 4.87% during 1990–2017 (Table 5). Regarding human settlements and water bodies, they showed a positive trend with an increase from 0.03% and 0.01 in 1990 to 0.1% and 0.03 in 2017, respectively. In general, the primary forest was the land use that experimented a negative trend. The rest of the land uses showed surface gains. The rate of change obtained indicate that the secondary forest, the human settlements and the water bodies were the land uses with the greatest transformation rates, with 8.03, 12.58 and 27.48, respectively, for the period of 1990–2017 and with 10.68, 15.96 and 12.3, respectively, from 2005 to 2017. Figure 3 shows the area occupied by the land uses studied. Likewise, it shows the rate of change of these land use/land cover for the periods 1990–2005 and 2005–2017. The calculated global precision, based on the Kappa Index, presented values of 80%, 85% and 84% for 1990, 2005 and 2017, respectively.

Table 6 shows the land use/land cover change dynamics. The primary forest lost the greatest surface area (28,406 ha) during 1990–2005, increasing the surface lost to 63,546 ha during 2005–2017. In contrast, the secondary forest showed the largest increases in area with 87,800 ha in the period 1990–2017.

Transition matrix

The transition probabilities of LULCC for the periods 1990–2005 and 2005–2017 are shown in Table 7. The diagonal of the matrix represents the permanence probability, i.e., the probability of a LULC type to remain unchanged. The areas without vegetation showed a 90% probability of transition from 1990 to 2005, lowering it to 62% from 2005 to 2017. The areas of primary forest presented a negative trend with a 71% probability of permanence in the period 1990 to 2005, and changing it to 61% for the period 2005–2017.

Weights of evidence (WoE) analysis

The WoE of the 19 variables were analyzed to eliminate those values that were above 0.5, based on the Cramer Coefficient (V). The distance to urban locations showed positive values of WoE from 1,000 to 9,000 m distance and from 42,000 to 47,000 m indicating an influence for cover change from secondary forest to area without vegetation. The distance to rural localities showed positive values of WoE in distances from 0 to 700 m. The topographic position index showed positive values in the ranges of −150 to −60 and 120 to 240. The distance to sawmills indicates that deforestation appears from 0 to 16,000 m with respect to the process of change between secondary forest to areas without vegetation. The transition from primary forest to area without vegetation is likely to occur in distances to the main roads between 13,000 and 21,000 m. The density of main streams such as rivers and creeks had an influence in densities from 0.039 to 0.079 m²/km².

In the transition from primary forest to secondary forest, the variable altitude showed positive values of WoE in the range of 1,200–1,300 m, suggesting that most of the changes occur in this range. The slope showed that the process of change between primary forest and secondary forest is located on slopes of 45–60 and 60–75 degrees. The transition from primary forest to human settlements was influenced by the distance to secondary streams from 500 to 1,000 m. The distance to sawmills presented an influence from 0 to 6,000 m. The distance to mines showed that the attraction to change occurs between 2,000 and 10,000 m.

Model validation

The model validation was carried with the simulated and the true land use classification of 2017. The FSI was applied for neighborhoods from 1 ×1 to 7 ×7 pixels. The minimum value reported for FSI was 49% in 1×1 pixels, while in 7×7 pixels the value of FSI was 91%. These results indicate that the real and simulated land use changes agree from 49% to 91%. Simulation starts with 49% and adjusts to 91%, reaching a similarity adjustment value at a distance of 210 m. These results agree with that obtained by Ximenes et al. (2011). According to Soares-Filho et al. (2017), and similar studies (Carlson et al., 2012; De Rezende et al., 2015; Elz et al., 2015), for the resolution and the number of transitions considered in the model, the values obtained for the FSI suggest that the models are good and can be used in the simulation of LULCC scenarios. Figure 4 represents the FSI in relation to the size of the window.

Scenarios

The LULCC based on the transitions between 1990 and 2017 for the stationary, optimistic and pessimistic scenarios are presented in Table 8.

Figure 5 shows the LULC classification of 2017 and the stationary, optimistic and pessimistic scenarios for 2050, after the model calibration.

In the stationary scenario the area without vegetation would increase from 4.8% in 2017 to 5.27%. Likewise, the secondary forest would increase from 57.7% (2017) to 73%. For this scenario, the changes in human settlement and water bodies would not increase or reduce their area. Conversely, the rate of change of primary forest and secondary forest were the greatest between 2017 and 2050. Regarding the optimistic scenario, it showed reductions in areas of primary forest; however, in lower magnitudes than for the stationary and pessimistic scenarios. For the pessimistic scenario, the Markov matrix was modified considering a greater pressure on the forest ecosystem. The area without vegetation showed a positive trend, with 4.8% in 2017 and an increase to almost 8% in 2050. The secondary forest would go from 57.7% to 85.6% in 2050. Finally, the primary forest would reduce its area to a 8% and isolated forest areas would appear. The rate of change for this scenario were the ones that showed the highest values. The LULCC dynamics projected for 2050 for the three scenarios (stationary, optimistic, pessimistic) is presented in Table 9.