Study area

The watershed area that makes up Tomatlán municipality is equivalent to 552 680 ha (Figure 1). According to the Comisión Nacional del Agua (National Water Commission) (Conagua), it belongs to the VIII^th Administrative Hydrological Region, Lerma-Santiago-Pacific, and to the 16^th Hydrological Region, Armería-Coahuayana, and includes the watersheds of Ipala, María García, Tomatlán A, Tomatlán B, and San Nicolás A rivers (^{CEA, 2015}). It is located in western Jalisco, at 19°56’3" N and 105°14’8" W (^{Inegi, 2009}). The study area is included in a decree on water reserves for the environment, issued on June 6, 2018 (^{DOF, 2018}).

Figure 1 Watersheds that make up Tomatlán municipality, Jalisco. 

Vegetation covers and land use maps of 1997 and 2017 were based on satellite images taken by the Landsat platform during the dry season, on May 4^th by the TM sensor, and on May 27^th by the OLI-TIRS sensor. Subsequently, atmospheric corrections, corrections for the brightness temperature, and a false-color image with bands 3, 4 and 5 were made. Training areas were determined for each cover, and a supervised classification was carried out with the adjustment of the maximum-likelihood algorithm. At the same time, the sample size was calculated in order to proceed to the field verification and the thematic evaluation; likewise, the confidence interval and the Kappa coefficient were estimated for each forest cover and land use map. Change rate was estimated for the 1997-2017 period in order to obtain the change matrix proposed by ^{Pontius et al. (2004)}, which indicates the tendency of the changes.

The watersheds were compared in terms of the forest covers and land use; the result was the change rate and matrix.

The images for this study were taken by various satellites: Landsat 8 OLI-TIRS and Landsat 5 TM. The scene of the take corresponds to PATH=030 and ROW=046 of the WRS2 system of Landsat 5, which has seven reflective bands with a 30 m spatial resolution.

Preprocessing or preclassification

This methodological phase included the application of processes related to the DOS1 atmospheric and radiometric correction proposed by ^{Chávez (1996)}. Based on the combination of TM3 (visible spectrum), TM4 (near infrared) and TM5 bands (visible red), it was possible to identify and interpret the analysis categories, and, therefore, 140 training fields or areas (number of representative or prototypic pixels per class or category) were digitalized. The location and spatial distribution of the training fields were carried out in a random, homogeneous way (20 fields minimum for each one of the covers) and through in situ rounds, along with the use and management of a global geopositioning system. The delimitation of the training areas (total number of pixels) was carried out based on ^{Jensen’s general rule (1996)}, which considers at least ten times the number of bands utilized for developing the statistics or classification parameters. The classifier was calibrated using the prototypical pixels of each category (training areas), which were graphically represented by means of histograms in order to obtain the spectral signatures for each category. These were subjected to a separability analysis (^{García and Mas, 2008}) based on the transformed divergence (TD) method, which uses the maximum-likelihood algorithm proposed by ^{Maussel et al. (1990)}.

Classification

Classification was made by the maximum-likelihood algorithm, which considers the previously estimated classification parameters (calibrated classifier) for labeling each pixel of the image in one of the defined classes (^{Lira, 2010}). The land use and vegetation maps of the watersheds for 1997 and 2017 were generated based on this algorithm. The methods, techniques and procedures described above (pre-classification and classification) were implemented with the ArcGIS 10.5 and ENVI 5.3 softwares.

Validation of the 1997 and 2017 forest cover and land use maps

In order to determine the degree of certainty of the data of the forest cover and land use maps, these were subjected to a process of evaluation of their thematic reliability, which was performed at the following methodological phases:

Sampling design. This consisted in determining the sampling units (pixel), the selection method or simple random sampling, and the number of sampling units or sample size. The latter was determined according to ^{Congalton (1988}; ¹⁹⁹¹), who suggests verifying at least 50 points per category. This principle was applied to the classified images, and seven categories or analysis classes were determined through the verification of 350 sampling points.

Evaluation of the verification sites. The sampling points were registered in situ; i.e. their spatial location (UTM coordinates, area 13N, datum WGS84) was determined by means of field rounds (direct observation) and with the use and management of a Global Geopositioning System; inaccessible points were verified using Google Earth 7.3 images. The sampling points were overlayed in the land use and vegetation cover maps (1997 and 2017) resulting from the visual interpretation of the satellite images, in order to associate them with the categories or classes defined in each of the classified maps.

Data analysis. This was based on the development of a confusion matrix or error matrix, which is bi-dimensional and on which the rows usually represent the reference classes (truth-ground), and the columns, the categories of the map (ground truth-image). The diagonal of the matrix expresses the number of verification points, for which a correspondence exists between the categories of the classified map and the reference data, while the marginal points indicate assignation errors. Changes were made by the methodology proposed by ^{Pontius et al. (2004)}, which consisted in developing a cross-tabulation or change matrix.

Kappa index

The multiple relationships between the various categories were analyzed with multivariant categorical analysis techniques, which are ideal for the global study of confusion matrixes (^{Campbell, 2007}). One of the most commonly used indexes is the Kappa statistic (k), which measures the difference between the map-observed reality agreement and the agreement that is expected to occur at random. It attempts to delimit the degree of fit due only to the accuracy of the classification, dispensing with the one that is caused by random factors. The purpose of the test is to evaluate whether the classification discriminates between the categories of interest with a significantly greater precision than that which would have been obtained with a random assignation (^{Chuvieco, 1990}; ^{Morales, 2014}). The following formula was used:

Where:

K = Kappa coefficient

N = Equal to the total sampled data

X¡¡ = Indicates the observed agreement

(X¡* X+¡) = Expected agreement (product of the marginals)

The obtained value represents the percentage in which the classification is better than the randomly expected one; if the obtained value is 0.80, the classification of the map is 80 % better. Three classes of results are considered: a value of K<0.4 represents a poor agreement; between 0.4 and 0.8, a moderate agreement, and a value above 0.8, a strong agreement.

Change rate in the 1997-2017 period

Based on the data derived from the change matrix (surface area of each category), change rates were estimated by forest cover and land use for the analysis period (1997-2017). The best way to quantify them is by means of the analysis of the dynamics of the vegetal cover, relating it to space and time (^{Berry et al., 1996}). ^{FAO (1996)} proposed Equation 2, which expresses the change in the surface area of the initial year as a percentage:

Where:

CR= Change rate (%)

S1= Surface area at the initial date

S2= Surface area at the final date

n = Difference of years between the initial and the final dates

Transition matrix

This was obtained through the overlaying of two forest cover and land use maps with different dates and the resulting data analysis arranged in a matrix, which allowed efficient identification of changes between categories. The matrix is arranged in such a way, that the cover classes of the initial date (T₁) are placed on the rows, and the categories of the final date (𝑇₂) are placed in the columns, whereby only the diagonal of the matrix shows the total amount of stable landscape between the two dates, while the rest corresponds to all possible combinations of exchange between categories (Table 1).

Gain index

The proportion of the landscape that experiences an increase between the initial date (𝑇₁) and the final date (𝑇₂) corresponds to 𝐺 𝑖𝑗 and is calculated by estimating the difference between the column of total time 2 𝑃 +𝑗 and the persistence 𝑃 𝑗𝑗 :

Losses index

Loss Lij, represents the proportion of the landscape that dwindles between dates; it is estimated by means of the difference between the row of total time 1 Pj+ and persistence Pjj:

Forest cover and land use maps

Two forest cover and land use maps were generated for the years 1997 and 2017 (Figures 2 and 3) by applying processes related to the radiometric and atmospheric correction, the interpretation of satellite images (Landsat TM and OLI-TIRS), and the supervised maximum-likelihood method. Seven types of vegetal covers and land use were determined: 1. Urban area (𝑍𝑈); 2. Temperate conifer and broadleaf forest (𝐵𝑇𝐶𝐿); 3. Agricultural (𝐴𝐺); 4. Grasslands (𝑃𝐴); 5. Low deciduous tropical forest (𝑆𝐵𝐶); 6. Medium growth semi-deciduous tropical forest (𝑆𝑀𝑆𝐶); and 7. Water bodies (𝐶𝐴).

Figure 2 Forest cover and land use map for the year 1997. 

Figure 3 Forest cover and land use map for the year 2017. 

The images were divided into seven general classes that combined more specific covers (Table 3), for if the study area were to be divided into a larger number of classes, the data contained in the training area would exhibit variability and often polymodal distributions, and, therefore, would detract from the efficiency of the maximum-likelihood algorithm (^{Campbell, 2007}).

Kappa index

Results from the process of evaluation of the thematic reliability show that the forest cover and land use maps for 1997 and 2017 were 79.54 %, 0.71 and 83.40 %, 0.77, respectively. The confidence interval and the Kappa index turned out to be acceptable in all cases (Fm >75 %) and (k<0.4) (^{Lesschen et al., 2005}). It must be mentioned that the validation of these maps was satisfactory, a product of the application of the simple random sampling method (^{Congalton, 1988}), as well as of the field collection and spatial location of the verification sites (sample size) (^{Valdez et al., 2011}).

Change rate

Table 4 lists the forest covers and land use, plus the surface area and the percentage for 1997 and 2017 in each of the assigned categories.

The change field refers to the hectares that increased and diminished in each category (Table 4), and annual rate is the increase or reduction in hectares per year, during the analyzed period. As for the cover change rate, the expansion of agricultural areas and grasslands is the main deforestation factor with an impact on the low deciduous forest, with 0.8 %. These values are similar to those cited in other areas of the country.

Deforestation rates of 4.2 % have been estimated in the forest of Los Tuxtlas, Veracruz (^{Dirzo and García, 1992}); of 1.58 to 2.3 % for Huistán and 0.46 to 3.42 % for Chanal, in Chiapas (^{Ochoa-Gaona and González-Espinosa, 2000}). Averages of 3.31 to 2.14 % have been documented for the Lacandona rainforest (^{Mendoza and Dirzo, 1999}), and of 1.4 % for the state of Morelos (^{Trejo and Dirzo, 2000}). The urban areas were the covers with the slightest change, i. e. 0.6 %; this is considered to be the result of the low land availability for the construction of infrastructure in the lower parts of the watersheds, as about 179 487 ha of the surface area are covered with grasslands and agricultural lands, and 8.0 %, by lower deciduous tropical forests.

Change matrix for the 1997-2017 period

The columns of the matrix represent the (aerial) portion of the categories of date 1 (𝐹₁); the rows contain the proportion (surface area) of the categories of date 2 (𝐹₂) (Table 5). The main diagonal corresponds to persistence (in red fonts) of the categories preserved between the dates 𝐹₁ and 𝐹₂ for each period. In the temporal analysis of the dynamic of the changes occurred in a period of twenty years, the grasslands have been observed to maintain a cover of 81 517 ha; the cover that most contributed to this increase was the medium growth semi-deciduous tropical forest, with a surface area of 28 530 ha; while the greatest loss was for the medium growth semi-deciduous tropical forest, which preserved a surface area of 92 906 ha.

Forest cover and land use change maps for the 1997-2017 period

The end product was a forest cover and land use change map indicating the percentage of change and the percentage to which the analyzed covers were preserved (Figure 4).

Figure 4 Changes in the forest covers and land use. 

In general, Figure 4 shows all the changes observed in the vegetation cover and land use, as well as the structures that were maintained. The total of covers that exhibited changes was 23 %, while 70 % remained unchanged; this is similar to the estimations cited by ^{Camacho-Sanabria et al. (2017)} in a study performed in Progreso, Hidalgo State, where 95 % of the cover was preserved and land use changes affected a mere 5 %.