Study area

The research was carried out in the area located between Maderas del Carmen and Ocampo Natural Protected Areas, both included in the Wildlife Protection Area (WPA) category (^{Conanp, 2016}). The study area encompasses part of the Acuña, Múzquiz and Ocampo municipalities in the northeast of the state of Coahuila (Figure 1).

Figure 1 Location of the study area in the WPAs (APFF) Maderas del Carmen and Ocampo (APFF), Coahuila. 

Acquisition and processing of satellite images

Satellite images available from the Landsat programs for the years 1985, 2001 and 2019 were located and downloaded from the Global Visualization platform of the U. S. Geological Survey (^{GloVis, 2005}). Other images, for the 1986, 2002 and 2020 years, were also obtained, all of them in a 30 X 30 m spatial resolution. Images of the rainy season were also incorporated including the red and infrared bands: b3 and b4 from the 1980´s Landsat 5 and from Landsat 7 of the beginning of 2000, and b4 and b5 from the 2019 Landsat 8, in which the following characteristics converge: the lowest possible percentage of cloudiness and no default remote sensing scanning to confirm the graphic tendencies in the calculated index, or identify inconsistencies or lack of data for the rainy season, between August and October.

The season of the images was selected with the expectation of the maximum spectral response as a result of the annual growth of the vegetation. The selection of the season of the year for the images was made expecting the maximum spectral response due to the annual growth of the vegetation. In addition, the WorldClim 2.1 database of historical monthly climate data of 1970-2000 was consulted (^{Fick and Hijmans, 2017}). The variables utilized were the monthly mean temperature (°C) and the monthly precipitation (mm).

The dates September 24^th, 1985, October 14^th, 2001, and October 8^th, 2019 were selected based on this analysis. An atmospheric correction was applied to all bands. Compound images were created to include the red and the near infrared bands according to the platform and sensor of each image. The plots located in the WFA Maderas del Carmen and Ocampo were selected, based on the property files of the region in order to generate a range polygon of the area of interest, in which a 1 000 m buffer was generated in order to in order to crop the three satellite images (Figure 1).

Supervised classifications

In terms of the different times elapsed between the three images and in order to generate three training areas in the supervised classifications, the thematic vectorial covers of the National Institute of Statistics and Geography (Instituto Nacional de Estadística y Geografía, Inegi) regarding the Soil Use and Vegetation (1:250 000) were utilized as a base and reference for selecting sample polygons as spectral signatures of the various types of vegetation or ecosystems; series II for the 1985 image, series IV for the 2001 one, and series VI for the training areas of the 2019 image. The soil use and vegetation types corresponded to the Inegi series: forests (all their variants and combinations, including chaparral), microphytic and rosetophyle scrubs, grasslands (natural, induced, and halophytic), rocky surface area (bedrock and alluvium), and surface with water (permanent rivers and small dams).

Recodings were performed at the ecosystem level in order to obtain forest, scrubland, and grassland areas. The coding of sample polygons used as a spectral signature was utilized in the digital classification processes. The ArcGIS 10.8.1 image classification module of ESRI^© 2020 was utilized to classify the three images with the Maximum Likelihood Supervised Classification method. For this purpose, a unique set of spectral signatures was edited and used for each processed date.

Post-classification processes

Once the classification and the soil use and vegetation map were generated for each date, with the initially assigned class names from the utilized training areas preserved in each map, the original classes were visually validated with Google Earth Pro and with two visits to the study area during the summer of 2020, in addition to a previous knowledge of the site since the year 2000. The results were vectorized in order to generate shapefiles that preserved a field with the name of the class used as training. Additional fields with the name and the ecosystem made a revision and a reclassification possible, in which the generated spectral signatures were recoded to the types of soil use and vegetation created according to the assigned classes.

Once the vegetation types and ecosystems were validated for each date, the three dominant ecosystems in terms of landscape level coverage in the area were analyzed, and the areas were calculated in square kilometers for the three dates, in order to analyze the changes between 1985, 2001 and 2019.

In order to ensure the recording of the vegetation index values calculated for all ecosystems over the time considered, the three vegetation maps were intersected through a geoprocess in order to locate those areas with permanence of the three most widely distributed vegetation types ―forests, scrublands and grasslands―, and thus the calculated values of the SAVI index between 1985 and 2019 were recorded. This made it possible to register the change in values during the period and to interpret their cover.

Use of the SAVI index

In each original image cropped to the boundaries of the study area, the R (red) and NIR (near infrared) bands of each date were used in the generation of the SAVI index with the Map Algebra module of Quantum GIS^© 3.16 (^{QGIS Development Team, 2020}); vegetation indexes were generated for each satellite image (Landsat 1, 2, 5 or 8, depending on the date) and band number (b) of the electromagnetic range corresponding to the red (R) and near infrared (NIR) bands to estimate the Soil Adjusted Vegetation Index (SAVI). Formula (1) was applied in order to calculate the index for each date.

Where:

R = Red band.

NIR = Near infrared band.

L = Constant value of 0.5 for soils in vegetation considering a semi-arid and temperate climate in the region (^{White and Swint, 2014}).

Classified vegetation types (ecosystem) and plots with different management intensities in the study area

The total area analyzed is 2 800 km²; 50.8 % of it are subjected to diverse management activities, while the remaining 49.2 % have been under conservation since 2000. There is a detailed background of vegetation classification (^{Muldavin et. al., 2004}) considering up to 33 plant associations grouped into five communities. However, in order to simplify this study regarding the analysis of vegetation types, and given that the difference in spectral response is considered to be greater as the cover increases at the ecosystem level, the results of each classification were based on the location of training areas used as a sample of the main ecosystems existing in the area, according to the Land Use and Vegetation series II and VI of Inegi (^{INEGI, 2001}; ^{INEGI, 2016}), in which six land uses and vegetation were identified, some in secondary condition ―verified through field visits and Google Earth―, by direct observation of the polygons of the training areas.

As a post-processing step, all polygons smaller than 10 000 m² (1 ha) were eliminated in order to avoid the "salt and pepper" effect. Subsequently, the polygons were reclassified into three major generalized ecosystems. In the reclassified and vectorized maps resulting from these classifications, the three most extensive and widespread ecosystems ―forests, scrublands, and grasslands― were identified.

Finally, based on a vector file with the zone's property boundaries and on the consultation of the property regime, plots with ejido community management, private plots with livestock management, and private lands without livestock and with conservation vocation since 2000 were recognized. A field was added to the database for the type of management, divided by intensity of use. The type of management was captured in the property file and intersected with the ecosystem polygons for each of the three dates analyzed, so that maps indicating ecosystems and management type were generated for 1985, 2001 and 2019.

Capture of SAVI index values by ecosystem and management regime

Within the distribution area of the three classified ecosystems, a point type shapefile was made on a grid pattern, systematically divided every 500 m using the Repeating Shapes module for ESRI (^{Jenness, 2012}). The SAVI index values corresponding to the three dates analyzed were captured in the shapefile thus created, and the ecosystem and management regime records were kept using the Extract Multi Values to Point module of ArcGIS 10.8 (^{ESRI, 2020}). Thus, a database of 22 250 records was generated with the UTM coordinate of each specific record in the study area. The base arrangement was organized using the SAVI index values for three dates as the main variable, and the ecosystem factors (forest, scrubland and grassland areas), and a management factor (conservation areas, private areas with traditional management and communal ejido areas with extensive exploitation). The methodological scheme is described in Figures 2a and 2b.

Figure 2 (a) Methodological scheme of preprocessing and digital classification of Landsat images of the study area, acquisition, atmospheric correction, cropping, composite images, training areas and supervised classification; (b) Verification and recoding of classes by ecosystem and generation of a network of points to capture these values and ecosystem and management factors. 

Statistical data analysis

The resulting database was converted to a DBF vector file format from the shapefile to Excel. Each of the columns was identified with a representative name, and the data was separated by ecosystem for analysis through filtering, and the intensity of management by ecosystem in databases separated by date was left as a factor for comparison analysis, in order to test for normality. The mean values by ecosystem and date were plotted to observe the trend over time. Normality tests were performed using Pearson's Chi-square test with R Commander of R^© (^{Fox, 2020}) for SAVI ecosystem and management data on all three dates. When normality was not observed in the distribution of the data for the analysis of variance, the Kruskal Wallis H-test was applied to determine differences between the means of each management intensity in the three ecosystems and the three dates, using the Statistix 8.1 package (^{Statistix, 2003}).

The results were concentrated by ecosystem and date to analyze differences between the means of each management intensity by ecosystem and to determine in which ecosystem there was greater coverage, based on the mean value of the index.

Spatial analysis of ecosystems and tendencies between 1985 and 2019

According to the results of the supervised and recoded classification at the ecosystem level, of the three main ecosystems in the area, forests have had a spatial reduction of just over 300 km² between 1985 and 2019. The areas of scrubland have exhibited only minor changes, decreasing by approximately 4 km² between 1985 and 2019. On the other hand, grasslands have increased by almost 3 600 km² over a 34-year period (Figure 3a).

Figure 3 (a) Change analysis in forests, scrublands, and grasslands; average SAVI vegetation index values between 1985 and 2020 under three different management intensities; (b) Forests; (c) Scrublands; (d) Grasslands. 

Tendencies in SAVI index values

The average values of the SAVI index calculated for forest areas showed a gradual decline from 1985 until 2002, when the value began to recover, until 2020; the values were always higher for the managed private areas dedicated to conservation, while privately lands under management exhibited intermediate values. Finally, the ejidos registered the lowest index values on all dates (Figure 3b).

In scrublands, a gradually decreasing trend was observed from 1985 to 2002, with an increase in 2019; the highest values corresponded to the land under conservation management at the beginning of the period and towards the end, with similar figures in the years 2002 and 2019; private farms with management registered the lowest values during the entire period analyzed, while community-managed lands had intermediate values between those of private lands and areas with conservation status (Figure 3c).

The grasslands exhibited a similar tendency to that of the forests and scrublands, with a gradual decrease in values between the 1980´s and the beginning of the 21^st century, with an increase after 2002. In this ecosystem, the lands dedicated to conservation did not have higher index values. The highest values were obtained from 2002 to 2020. Among the grasslands, managed private plots registered the highest SAVI indexes but have exhibited intermediate values in the last 20 years. Finally, the areas with communal management had a general constant of similarly low vegetation indexes from the beginning of 2001 to 2002, and in the final years of the study they exhibited the lowest SAVI index among the different property regimes (Figure 3d).

Data normality tests

In the analyzed years, between 1985, 2001 and 2019, all the mean SAVI values for the three main ecosystems yielded lower probability values than the Pearson Chi-square p-values calculated for forests (41.6-105.3), scrublands (431.9-2 078), and grasslands (231.9-2 342); i.e., there was no normality in data distribution (with p-values of 0.032 to 2.2e^-16).

Management factor in the forest ecosystem

According to the non-parametric analysis of variance in the Kruskal-Wallis H-test, in 1985 there were significant differences (α=0.05; critical value Z=2.394) between private plots with conservation management, private plots with management, and ejidos. The first had the highest value, followed by the private plots with management, and the ejidos had the lowest average index value. In 2001, there were differences between private conservation and managed plots and land under communal management (α=0.05, critical value Z=2.394); communal management presented the lowest index value. In 2019, both managed and conservation private plots differed from ejidos according to the H-test (α=0.05, critical value Z=2.394); the lowest value corresponded to the plots under community management (Table 1).

Management factor in the scrubland ecosystem

In scrublands, the non-parametric analysis of variance in the 1985 Kruskal-Wallis H-test indicated significant differences (α=0.05; critical value Z=2.394) between private land with conservation management, private land with management and ejidos; from these, the first exhibited the highest value, followed by the private plots with management, and the lowest average index value was estimated in the ejidos (Table 1). Significant differences persisted in 2001 (α=0.05, critical value Z=2.394) between the three types of management; the lands with the highest values were those under community management, followed by private lands under conservation management, while private lands under management had the lowest average value. In 2019, there was a change, and no significant differences were detected between management intensities for the scrublands (α=0.05, critical value Z=2.394) (Table 1).

Management factor in the grassland ecosystem

The non-parametric analysis of variance of the Kruskal-Wallis H-test showed significant differences for 1985 in the grasslands (α=0.05; critical value Z=2.394) between private lands under conservation management and private lands under management and ejidos; the former had the highest value, followed by the private plots with management, and the ejidos had the lowest average index (Table 1). In 2001, there was a significant difference between the three management intensities according to the Kruskal Wallis H-test (α=0.05; critical value Z=2.394); in this case, private plots with management exhibited the highest value, followed by the ejidos, while the plots under conservation management had the lowest rate (Table 1). In 2019, the situation changed because significant differences were detected by the H-test (α=0.05; critical value Z=2.394); the areas with conservation management registered the highest value, followed by privately managed land and community-managed plots, which were statistically equal, below the other two management types (Table 1).