*2.3. Methods*

#### 2.3.1. Data Preprocessing

Primary data were obtained and processed according to the methods presented in Table 2, while a flow chart illustrating the methodology employed in this study is presented as Figure 2:


**Table 2.** Data and methods of potential factors.

Notes that LUC, P, GDP, A, KI, DI, MAP, MAT, and S represent land-use change, population, gross domestic product, accessibility, karstification intensity, drought index, mean annual precipitation, mean annual temperature, and slope, respectively.

**Figure 2.** Flow chart of methodology employed.

Forest change: Based on a Markov transition matrix [55] and spatial analysis function of ArcGIS, a transfer matrix of different land types was obtained, and on this basis, land-use changes, including forest cover, were estimated. Chord diagrams of land transformation were constructed for Guizhou province and each of its nine major municipalities using R.

Land-use transition involves the changes in regional land-use patterns, and they are significant components of land-use-transition studies [56]. Markov modeling is commonly used to consider the processes and mechanisms of landscape-dynamics changes over the longer term [55,57]. We adopted a transition matrix as the core part of the Markov model (see Formula (1)), which is generally applied in estimations of land-cover changes [58,59]. While different types of conversion may occur, more attention was paid to those that account for most of the forest change. In addition, we calculated the conversion ratio of the main converted types through Formula (2). All analyses were conducted by ArcGIS and R.

$$\mathbf{T} = \begin{bmatrix} D\_{11} & D\_{12} & \dots & D\_{1n} \\ D\_{21} & D\_{22} & \dots & D\_{2n} \\ \dots & \dots & \dots & \dots \\ D\_{n1} & D\_{n2} & \dots & D\_{nn} \end{bmatrix} \tag{1}$$

where T refers to the conversion matrix of different types of land-use change from 1980 to 2018, *Dnn* refers to the change in the area (unit: km2) from one land-use type to another during the study period, and n refers to the area of a certain type of land that was involved in the computation.

$$R\_{ij} = A\_{ij} / B\_i \tag{2}$$

where *Rij* refers to the rate of the *i* type of land use converted to *j* type, *Aij* represents the area of *i* type of land use converted to *j* grade, and *Bi* is the total area of *i* type of land in 1980 (unit: km2).

Population and Gross Domestic Product (GDP): These data were analyzed spatially in ArcGIS 10.3. According to REDCP, population data were processed as follows:

$$\text{POP}\_{\text{ij}} = \text{POP} \times \left( \text{Q}\_{\text{ij}} / \mathcal{Q} \right) \tag{3}$$

where POP*ij* is the population-spatial-distribution data in a 1 km × 1 km grid, Q*ij* is the total weighting of land-use type, night light, and settlement density in a grid, POP is the population of the county-level administrative unit in which the grid is located, and *Q* is the total weighting of land-use type, night light, and settlement density for the county-level administrative unit in which the grid is located.

According to REDCP, GDP data were processed as follows:

$$\text{GDP}\_{i\bar{j}} = \text{GDP} \times \left(\text{Q}\_{i\bar{j}}/Q\right) \tag{4}$$

where GDP*ij* is the GDP-spatial-distribution data in a 1 km × 1 km grid, Q*ij* is the total weighting of land-use type, night light, and settlement density in a grid, GDP is the GDP of the county-level administrative unit in which the grid is located, and *Q* is the total weighting of land-use type, night light, and settlement density for the county-level administrative unit in which the grid is located.

Accessibility (A): This parameter refers to the accessibility of a location in terms of transportation, including national, provincial, county, and township roads. The Euclidean distance was calculated for each of the different levels of road, weighted according to ranking of their importance (national > provincial > county > township), and then analyzed spatially.

Karstification intensity (KI): The degree of karstification was classified as follows, according to the purity of carbonate: detrital carbonate-reservoir rock (DCRR), carbonate-reservoir rock with detrital reservoir rock (CRR-DRR), and non-carbonate rock (NCR). Weights were assigned according to karstification intensity, whereby DCRR > CRR-DRR > NCR.

Drought index and mean annual precipitation/temperature: A drought index [60] was obtained from Formulae (5)–(7). Next, based on Kriging interpolation in ArcGIS10.3, the map of drought index was made.

$$\mathbf{K} = \mathbf{E}'/\mathbf{P}'\tag{5}$$

$$\mathbf{P'} = \mathbf{P/P\_P} \tag{6}$$

$$\mathbf{E}' = \mathbf{E}/\mathbf{E}\_P \tag{7}$$

where K denotes the drought index of any period, P describes the relative rate of change in precipitation during the period (1980–2015), P represents annual total precipitation for 2015, P*<sup>P</sup>* is the mean annual precipitation during the period (1980–2015), E describes the relative rate of change in evaporation during the period (1980–2015), E represents the evaporation in 2015, and E*<sup>P</sup>* is the mean annual evaporation during the period (1980–2015).

Slope: Following image cutting and splicing, DEM was used to classify slopes, as follows: 0–6◦, 6–15◦, 15–25◦, 25–35◦, and >35◦.
