Open ecosystems in Mexico are under increasing pressure, due particularly to the expansion of cities and agricultural activities. These developments occur without integrating biodiversity concerns in land use planning and result in extensive fragmentation and transformation of the landscapes. The semiarid region of Mesa Central was characterized using ten precipitation-based indices. Using multivariate statistical and geostatistical spatial analysis techniques, the influence of those indices on five land use strata was explored. Land use analysis indicated that the maximum values of the five significant precipitation-based indices were found in Grasslands, Agricultural Use, and Shrubs; minimum values were characteristic of substrates Secondary Desert Vegetation and Other Use. Our results suggest that the greatest number of extreme precipitation events is likely to occur in open ecosystems and consequently will have a strong influence on landscaping and land use. The semivariogram analysis and geostatistical layers demand attention from research institutions, policy makers, researchers, and food producers to take the appropriate and coordinated actions to propose scenarios to deal with climate change. Perhaps this study can stimulate thought concerning research endeavours aimed at promoting initiatives for biodiversity conservation and planning programs for climate change mitigation.

1. Introduction

Climate results from a combination of atmospheric factors and environmental conditions that operate at different levels [1]. Indices for climate variability and extremes have been used for a long time, often by assessing days with precipitation observations above or below specifically based thresholds [2]. Because vegetation covers most of the global land surface, it strongly affects the land-atmosphere exchanges of energy, momentum, and materials [3] through the distinctive combination of interacting elements that are repeated in similar form through the landscape; it also impacts runoff and erosion rates along with soil stability. Most semiarid regions suffer severe rainfall erosion [4]; since water has a high erosive capacity, more damage is likely to occur in these areas because of reduced vegetative protection.

Climate system warming is unequivocal, and since the 1950s many of the observed changes are unprecedented when seen over decades, even millennia [5]. According to Zeng et al. [3], arid and semiarid regions are dominated by shrub communities with many subtypes. In addition, within the spatial limits of these regions, an exceptional richness of biodiversity and an astonishing variety of biomes are found; exceptionally high species diversity and endemism occur. Plant community spatial distribution is strongly influenced by rainfall distribution, soil water retention rate, and the soil itself as substrate. Additionally, temperature and geomorphic factors influence plant community distribution in the landscape by influencing evaporation rate and the amount of sun insolation. According to Kuyler [6], the origin and development of landscapes are influenced by a combination of natural processes and human influences. Land use, particularly its intensity, is considered one of the determining factors of the abundance and richness in populations of soil organisms [7]. The diversity and heterogeneity of land use processes require detailed analysis because of their differential effects on the environment. Land use degradation occurs when the ecosystem is not able to regenerate its structures and ecological processes (resilience) in order to recover its stability. Most semiarid regions are characterized by a hot summer and extremely limited and spatially erratic annual precipitation (less than 300 mm). Changes in temperature and rain events have received extensive attention from researchers, since many regions worldwide have experienced significant variations in climate extremes during the past few decades [8, 9]. The number and frequency of rainfall events along with extended droughts raise the question as to whether extreme climate events are truly increasing and also whether they are associated with climate change. Knapp et al. [10] documented growing evidence at global, regional, and local levels that intra-annual precipitation regimes have become more extreme; global precipitation records during the 20th century show an average increase of only 9 mm over land areas, excluding Antarctica. Therefore, according to Pearson et al. [11], climate data are essential input variables for ecological modelling and play a significant role in flora and fauna distributions.

Climate change impacts add further complication to the already demanding water management challenges in arid and semiarid regions. Its effect on hydrological processes such as infiltration, percolation, runoff, and soil water storage is highly complex [12]. In Mexico an issue demanding attention is limited and often lack of serial climate data records on open ecosystems because of the absence of meteorological devices. Along with the high cost of maintenance, the spatial representation of ground-weather monitoring stations is a function of regional climatic variability, measurement site characteristics, observational practices [13], and frequency of instrument calibration to standards. Some stations may be strongly influenced by an unusual or changing microclimate within their immediate surroundings and therefore have less utility representing climatic dynamism over large spatial scales. The scientific community has extensively documented the robust influence of geomorphic factors (elevation, slope [north, south, east, west facing, or level]), and topographic setting (ridge top, valley, etc.) as biased sources for recorded data and trend analysis. A station near hills or hollows, in or out of nocturnal drainage channels, or influenced by other subtle features may result in errors in climate data records. Concerning this issue, the World Meteorological Organization (WMO) contends that for developing high quality datasets the database must consider these potential limitations: general measurement and sampling errors; lack of homogeneity (where external factors can influence the record, e.g., new trees or buildings near the observation site); and the practice of having to make statistical averages over large geographical regions when limited data are available to represent those areas. According to Jamaludin and Suhaimi [14], dealing with these issues implies a high cost and an overwhelming number of procedures, in fact way too many. Scientists, accordingly, have developed a number of interpolating methods to keep inherent data uncertainty at a minimum to ensure enough accuracy for its use in climate analyses.

Ordinary Kriging is an interpolator that works well on any scale because of its accuracy in improving the calculation of model error. In a geographic information system (GIS) environment it uses statistical models that allow a variety of thematic output maps that include predictions, standard errors of predictions, probability, and quantile (ArcGIS). With the use of interpolation, the idea is to have an estimate of the distance one would need to travel before data points separated by that distance are uncorrelated [15]. This information is usually presented in the form of the variogram, which is a function of the semivariance versus the distance lag. The basic difference between variance and semivariance data is that the first category is represented by a number (scalar) and the second one by a curve (vector).

The use of indices instead of data is a better way to make the influence of climate change on the landscape more evident. An index represents information derived from the data itself and is useful for supporting a variety of climate change analyses. The RClimDex is an R based package designed to provide a user-friendly interface to compute indices of climate extremes [16].

This study endeavours to provide new information which integrates a serial dataset of direct measurements of daily precipitation in order to derive precipitation-based climate indices, their spatial mapping interpretation, and trend analysis. The study area is the Mesa Central, a semiarid physiographic province of Mexico. This region has highly fragmented land use and a landscape of contrasting elevation and landforms. Daily precipitation data of a sufficient number of ground-weather stations were used to derive climate indices datasets as well as the slope of the estimated trend line equation. All datasets were fully incorporated into a GIS project to extend the analysis to regional scale, taking into consideration the land use categories. The goal is to improve the knowledge of how precipitation-based indices may help to clarify the way regional land use may be affected by climate change.

2. Materials and Methods

2.1. Description of the Study Area

The Mesa Central is a physiographic province that is part of the northern high plateau region in Mexico and covers two-thirds of the country’s arid and semiarid regions. It is located in central Mexico (Figure 1) and embodies ~8.6 Mha. Its spatial limits were originally delineated based on morphologic and geologic characteristics contrasting with those of neighbouring physiographic provinces [17]. It is a large area of folds and separated mountain failures with vast flood plains that become part of large and extended flat basins. Two physiographic features characterize this region: a semiarid area of volcanic cones separated by high alluvial basins adjacent to Sierra Madre Oriental (SMOr) and, secondly, desert lands and semiarid areas along with inner basins, located in the lower central zone and eastern region of the plateau [18]. The province combines rough topographic relief and flat areas. Over half of its surface is above 1,900 masl (Figure 1(a)) and the inside topographic elevations are moderate, generally forming ramps at 600 m or less. On the north and east it is bordered by the SMOr, on the west by the Sierra Madre Occidental (SMOcc), and on the south by the depression known as El Bajio. The spatial limits are in accord with the Instituto Nacional de Estadística y Geografía (INEGI) of Mexico. Regional land use incorporates five gross classes: the Agricultural Class, covering ~2.65 Mha (30.81% of the area); Shrubs (chaparral), comprising ~2.76 Mha (32.09%); Secondary Desert Vegetation (SDV), embodying ~1.89 Mha (22.98%); Grassland, totalling ~0.86 Mha (10.00%); and Other Use, accounting for the rest of the study area ~0.44 Mha (5.12%). Agricultural activity is under both irrigation and precipitation-dependent regimes; the shrub class includes crassicaule, rosetophyll, microphyll, and submountainous species; the SDV stratum has holm oaks, pine, mesquite, and their associations; the Grassland category incorporates halophytic, along with introduced and native species (Figure 1(b)); the Other Use stratum includes water bodies, human settlements, and bare soil.

Because of its considerable latitudinal size, the province was subdivided into three virtual regions, each corresponding to ~175 km of latitude distance. The Agricultural Class is regional and widely distributed. In the northern region the horizontal spatial distribution is dominated by a transition of the Shrubs, SDV, and Agricultural Classes. This region also displays an outstanding contrast in elevation with a transition from flat areas (northern) to peaks (east, central, and west). The southern region is a blend of all land use strata where the transition among classes is not specifically dominated by any particular class; there is, however, a distinct transition in elevation. Most of the Other Use class is located in this region, along with significant areas of the other land use strata. The central region exhibits a pattern of land use distribution comparable to the north region.

2.2. Data Description and Quality

Daily precipitation data of 32 ground-weather stations (Figure 1 and Table 1) were used to derive precipitation-based indices (Table 2).

The “” and “” columns refer to coordinate data: North American Lambert Conformal Conic Projection, Datum ITRF92. Bold numbers denote a simple consecutive list corresponding to the station numbers in Figure 1. The CNT column designates the number of years with data. Columns Yr STRT and Yr END indicate initial and final years of recorded data. Elevation data are provided in the ELEV column; data was acquired by point extraction from a digital elevation model (DEM) with ~30 m spatial resolution. The serial dataset range was 41 > years < 86. Elevation ranged from 1,594 to 2,310 masl.

Data quality was a priority because the indices are sensitive to changes in station location, sun exposure, equipment precision, and observer’s practice [19]. Before calculating the indices, the dataset preparation and processing instructions of the RClimDex user’s manual were followed, as proposed by Zhang and Yang [20] to insure data quality. Two phases are involved: the Quality Control (QC) and the Homogeneity Test (HT). The goal of QC is to identify errors in daily datasets that may potentially interfere with the correct assessment of the extremes; that is, all missing values (currently coded as −99.9) are replaced with an internal code that the software recognizes (i.e., NA: not available); and all allegedly aberrant data are substituted with the NA code. The quality test also identifies outliers in daily precipitation data. An outlier is an observation that falls outside the statistical limits of probability or else is negative. Currently these limits are defined as times standard deviation () of the value for the day; that is, (, ), where stands for the day and is an input from the user. Here we assigned a threshold value of four stdev’s. According to Vincent et al. [21], HT consists in the detection of shifts in climate time series, which are often due to station relocation, changes in instruments observation practices, and automation.

The RClimDex User’s Manual is the best source for a more detailed and extensive description of these indices; it is available for online consultation and downloading.

2.3. Statistical Analysis

All descriptive statistics were processed according to Chebyshev’s theorem; it formally states that “the proportion of observations falling within standard deviations of those numbers of the mean of those numbers is at least ” [22]. The advantage of this theorem is that it applies to any dataset regardless of the distribution shape of the data. For calculation, the serial maximum datum was indicated as the “upper limit” (1); then, datum was further added to (2): In (1), the calculated is the number of standard deviations that the serial maximum datum amounts to by dividing it by the standard deviation. Stated simply, the theorem gives the minimum proportion of data which must lie within a given number of standard deviations of the mean; true proportions found within the indicated regions could be greater than what the theorem guarantees. The correct interpretation is that at least Chebyshev’s value (percent) of each index is in the interval of and the rest of the proportion is outside this interval.

The geostatistical layers of surface map predictions were obtained by applying the Ordinary Kriging interpolation technique. It is considered superior to other commonly used interpolation techniques for precipitation estimation [23]. This technique provides an unbiased interpolation with minimum square mean estimation error and it is used for its ability to incorporate regional data and data indicating local trends. It has also been used frequently in soil science [24, 25], hydrology [26], and more recently in forestry [27], ecology, and climatology [2830]. Its calculation is based on a semivariogram, a geostatistical tool that analyses the spatial variability of one variable within the spatial limits of a specific area. A semivariogram is commonly used to model spatial structure in a single variable by measuring the strength of statistical correlation as a function of distance. According to Gringarten and Deutsch [31], the expected squared difference between two data values separated by a distance vector, h, is the variogram. The semivariogram (h) is one half of the variogram 2 (h). The variogram is a measure of variability; it increases as samples become more dissimilar. According to Webster and Oliver [32], the parameters that define the variogram are (1) the sill, which is the total variance and represents the variability in the absence of spatial correlation; (2) the range, which is the distance at which the variogram approaches the sill; and (3) the nugget effect, which is a combination of spatially unstructured variance (e.g., attribute error) and spatially structured variance at distances shorter than the minimum measurement separation. The sill minus nugget is known as partial sill or structural variance.

Because of the lack of the specific data for calibrating the semivariogram runs, we used fixed parameters for all runs; the primary variable was the precipitation-based index and the covariable was the elevation grid. Run-parameters were the model, ordinary co-Kriging; lag, 2,318.92; nugget, 0, 0; range, 204,168.71; and sill, 1.5631. According to Schabenberger and Pierce [33], the decision about whether or not to include a nugget effect is difficult. Its inclusion probably results in a better-fitting model [34], but considering that the geomorphology of the study area varies considerably and also that land use is disparate and delineated on a 1 : 250,000 scale, the decision to set the nugget value at 0, 0 seemed appropriate. This decision follows the recommendation of Schabenberger and Pierce [33] that nugget inclusion is not necessary.

There are several models of semivariance to choose from. According to Glover [15], the linear model seems to provide the best results because the data do not support evidence for a sill or a range; rather, they appear to have increasing semivariance as the lag increases. Gallardo [35] proposed a rule of thumb about the number of pairs to represent a point in the semivariogram; it must be greater than 30. This general rule means that the number of data on the space object would not be less than 50. The same author mentions that this rule serves as a guide and does not have to be followed absolutely. In this study the number of ground-weather stations is 32, which seems adequate to support the research.

For statistical analyses, the multivariate analysis of variance (MANOVA) along with the post hoc Bonferroni test for homogeneous groups was applied. To make the database for analysis suitable for the general MANOVA model, one which accepts categorical variables, the dataset was augmented with two categorical predictors: the three virtual regions (North, Central, and South) and the five land use classes. Resampling the geostatistical surface map layer with a randomly distributed dataset of 3,569 coordinate pairs provided the land use analysis. Each data point is an n-dimensional vector whose coordinates are , at a time () [36]. The required number of points, along with their distribution, was determined by applying a random distribution function and a restriction distance of 200 m between points. Because land use was originally delineated at a 1 : 250,000 scale and also divided into five appreciable strata for this study, it seemed imperative to select a minimum number of points by class. Calculating an appropriate sample size relies on the subjective choice of certain factors and sometimes crude estimates of others; as a result, this outcome may seem rather artificial. However, at worst it is a well educated guess, and it is considerably more useful than a completely arbitrary choice [37].

3. Results and Discussion

3.1. RClimdex

A variation in statistical significance of trend lines was evident. Of the ten precipitation-based indices, five were seen as significant (): PRCPTOT, R10 mm, R20 mm, R25 mm, and SDII. Hence, the discussion of the results and analysis will be limited to these indices.

To choose indices for mapping purposes and statistical analysis, along with a consideration of significance, we separated indices by applying two additional criteria: a minimum of 5 ground-weather stations had to result as significant for the index and their spatial distribution could not be concentrated in reduced areas. Of the calculated 320 trend lines (32 ground-weather stations × 10 indices), a positive trend was observed for 38 (~12%), whereas a negative trend was perceived for ~5% (15); the remaining trend line equations, 267 (~83%), were seen as nonsignificant (). These results appear to accentuate a regional spatial distribution for precipitation events which lacks homogeneity in frequency and intensity; they also seem to indicate that certain areas may be more susceptible than others to a strong influence of topography. Those regions may provide a good target for focusing research on the effects of land use and climate change for integrating plant communities and how this may influence spatial landscaping. As observed in Table 3, indices that met the proposed criteria and resulted in a positive trend were those related to heavy rain (R10 mm), total precipitation (PRCPTOT), days with rain ≥25 mm (R25 mm), and SDII; the only index displaying a negative trend was R20 mm. Hereafter the results and discussion will be limited to these five indices.

3.2. Surface Analysis

The variations and trend of climate indices as well as the observed semivariogram must be understood for reliable interpretation and modelling. The surface prediction maps on Figure 2 represent the individual spatial variability of each index along with the positive/negative trend of each ground-weather station. The legend (classes) for geostatistical surface layers is in accord with regional land use spatial distribution. Di Gregorio and Jansen [38] define land use spatial distribution as the spatial arrangement of life forms (e.g., trees, Shrubs, Grasslands, etc.) in a given region. Land use spatial distribution reflects an ecological or an evolutionary aspect of vegetation (e.g., scattered vegetation in arid areas, agricultural encroachment inside forest areas, and degradation due to overgrazing). In defining a particular classification design in a particular area two concepts is considered: the Minimum Mapable Area (MMA) and the Mixed Mapping Units (MMU). Di Gregorio and Jansen [38] mention that MMA is applied by cartographers when addressing the smallest area that can be shown on a map; it is scale-dependant and not related to classification. The same authors mention that the MMA is cartography related. The classifier is free to decide if land use implies an intricate mixture of classes where no one class is clearly dominant. Our study region was considered a spatially separate entity (e.g., patches of each stratum of land use inside all) in an intricate mixture (crops, Shrubs, and Grasslands). In the case of spatially separated entities of two or more classes, we followed the general criteria proposed by Di Gregorio and Jansen [38] that the cover of each class considered must be more than 20 percent (and consequently less than 80 percent) of the mapping unit. Based on these criteria, we defined four major classes for all significant indices.

The PRCPTOT index showed a pronounced effect on model performance because of the strong influence of the spatial autocorrelation range. Despite its designation as zero, a slight nugget effect was observed in two locations (Table 4). The nugget effect and the origin of slope are the most important features for fitting the semivariogram [39] and the nugget effect is also the most unpredictable because of the lack of closer samples in the dataset.

The discontinuity in the semivariogram’s origin can be viewed as a nonsense situation since landscape realities [40] and environmental conditions may be considered a continuous surface. According to Clark and Harper [41], the discontinuity of semivariogram origin should not occur in most spatial environmental processes because space is generally continuous; however, some individual processes in nature, such as gold mineralization, are exceptions to this pattern. Nevertheless, the semivariogram should be continuous at the origin. “One of our basic assumptions is physical continuity of the phenomena being measured” [41]. As evident in the PRCPTOT index, the nugget effect indicates erratic behaviour over very short distances along with considerable variability over distances less than specified lag spacing or sampling interval. This would imply that “perfect” sampling would eliminate the nugget effect entirely, but the methodology and the applied analysis techniques of this study are limited and circumscribed by the number of ground-weather stations; this limitation prompts us to consider the resulting value of the nugget effect as inconclusive. In addition, a key issue is the interpolation approach for a given set of input data [42]. This is especially true for areas such as those with contrasting topography (mountainous regions), where data collection is sparse and measurements for given variables may differ significantly even at relatively reduced spatial scales [43]. Burrough and McDonnell [42] state that when data are abundant most interpolation techniques give similar results. When data are sparse, the underlying assumptions about the variation among sampled points may differ and the choice of interpolation method and parameters may become critical. Results of spatial interpolation contain a certain degree of error, and this error is sometimes measurable [44]. As apparent in this study, land use in this semiarid region is highly fragmented and the number of stations seems to be locally adequate; this is notable in central region. According to Ahrens [45], the geographical representativeness of the stations diverges, but this is not systematically dependent on station elevation. Since this is not addressed by geographical distance weighting and by the knowledge that precipitation increases with elevation, it is expected that interpolated precipitation will be tendencially overestimated in those areas with rough terrain. One limitation that has been reported when using geostatistical spatial analysis is that the underlying stochastic concept relies on the distribution and on the stationarity of mean and spatial covariance [46]. Some studies have used prescribed data transformation: for example, Schuurmans et al. [47] use a square root transformation in application with situations of abundant daily rainfall; Verworn and Haberlandt [48] use a log transformation for selected flood cases. Nevertheless, in our study region precipitation events are erratic and even absent; to our knowledge the application of transformation is not a common practice, and the effect of transformation as well as the role of transformation choice has not been systematically examined so far. Nonetheless, this feature recommends further exploration for future research studies and may suggest incorporating additional datasets for remotely sensed net shortwave and longwave radiation, latent and sensible heat flux, and surface and subsurface runoff.

These results seem to support those reported by Wagner et al. [49] where they indicate that in an interpolation scheme the choice of a covariate had a significant impact on estimated precipitation and runoff amounts. Our results demonstrate that the incorporation of an elevation grid as a covariate should improve and enhance the spatial representativeness of ground-weather station datasets. Naimi et al. [50] mention that the value of the covariate of an interpolated location would be similar to that of an actual location. According to Gringarten and Deutsch [31] at some scales the processes are highly nonlinear and chaotic, leading to variations that have no spatial correlation structure. Typically, only a small amount of variability is explained by random behaviour. The results here seem to validate the erratic distribution and randomness of precipitation events in this semiarid region.

All semivariograms showed that spatial correlation to elevation decreases with separation distance. For all indices the size of the range of correlation scale depends on direction; that is, the vertical range of correlation is much less than the horizontal range due to the greater lateral distance between stations.

The interpretation of the geostatistical layers resulted in mean differences in all indices. Surfaces A and B resulted in an almost comparable spatial distribution with an evident latitude gradient that starts with the lowest area of the north region and extends to the highest area in the south region. In Figure 2(c), all of the north region and most of central region fundamentally had the lowest number of heavy precipitation days. For geostatistical layer D, the maximum number of days in a year with precipitation ≥25 mm is limited to a small area in the south region where the most of the Other Use class is present. For layer E, results indicated that the mean precipitation in wet days (6.8–8.3) was essentially homogeneous along the central north/south length of the study area.

3.3. Land Use

Based on the Chebyshev’s theorem it appears evident that the north region had the greatest dispersion among land use strata and significant indices indicating the proportion of observations from the mean (>41% proportion <82%). This was especially notable for the Other Use and Shrubs classes, but this was in contrast to the value which regionally showed the closest approach to the mean, as observed on the graphic with a >1.2 <2.4 (Figure 3). Consistently, the central region resulted in the highest proportion around the mean (72–91%), but in contrast with the highest range of (1.6–3.4); south region was influenced by the values of the SDII index resulting in a smaller proportion around the mean (47–66%), but the closest range from the mean (>1.4 < 1.7) (Figure 3). These results emphasize the variation among strata by the latitude gradient and land use.

Among significant indices the PRCPTOT indicated the lowest value for the Other Use stratum. This result may suggest evidence that anthropogenic activity is influencing the hydrologic cycle by affecting the total amount of precipitation within the region. Jacobson and Kaufman [51] reported a reduction in precipitation in California due to emissions of anthropogenic aerosol particles and precursor gases; that is, anthropogenic SOx, NOx, NH3, and speciated organic gases, but not CO2, CH4, or N2O. This feature deserves further research and extended study, but it is beyond the scope of this project.

According to Gringarten and Deutsch [31], the extent of spatial correlation decreases with separation distance until reaching a distant point where no spatial correlation exists. Virtually all climatic relationships impart trends which affect land use. By way of examples, water distribution upland or a systematic decrease in reservoir quality from proximal to distal locations may have an impact on an ecosystem’s integration processes, gas exchange, and functionality; likewise, geomorphology exerts a significant influence on the amount of solar radiation incident on a surface, since radiation varies with slope and aspect as a function of slope angle; elevation influences spatial distribution of wetness on soil, biotic composition of plant communities, pedogenic process, and runoff. This is especially important for terrain with highly complex geomorphology, where topography shading generally has important effects on energy balance. The cause for night time warming and variations in energy balance has been attributed to changes in atmospheric water vapour, cloud cover, jet contrails, and changes in surface characteristics such as land cover and land use [5254]. Local effects such as urban growth, irrigation, desertification, and variations in local land use can affect the wetness distribution, particularly in urban areas. Large-scale climatic effects on regional wetness distribution may include increases in cloud cover, surface evaporative cooling from precipitation, greenhouse gases, tropospheric aerosols [52], and extreme meteorological events, along with hydric and energy imbalance. This in turn may have consequences for runoff rate, sediment transport, groundwater recharge, and of course land use conformation.

3.4. Statistical Analysis

Data analysis was performed on the sample dataset points with added categorical variables for land use and regions. The MANOVA model provided significant outcomes () for both categorical predictors. As expected, most of the variation was observed in land use (Wilk’s lambda = 0.18, , ) and less by regions (Wilk’s lambda = 0.96, , ). The post hoc Bonferroni test showed significant differences () for regions and indices, except for PRCPTOT where no significance () was observed between the central and north regions (Table 5). The most outstanding results were found in land use.

The PRCPTOT index indicated significantly different results () when pairing the Shrubs, Grasslands and Agricultural Classes with the SDV stratum for all indices. For the remaining paired associations, no significance () was observed (Table 5). For the R10 mm, R20 mm, and SDII indices, results for a number of pairs were similar to PRCPTOT index; a notable exception was in the R25 mm index, which indicated two well defined groups. The first group included the SDV, Shrubs, and Agricultural strata; the second one encompassed the Other Use and Grasslands strata.

Regional variation is apparent by observing the latitude gradient. The south region consistently received the highest annual precipitation, the greatest number of extreme precipitation events, and the largest number of wet days per year (SDII index). Land use analysis indicated two well defined groups for all indices, except the most extreme index (R25 mm). Regionally the Shrubs, Agricultural Use, and Grasslands strata ranked highest in all indices, showing yearly more total precipitation, the greatest number of heavy precipitation days, and the greatest number of wet days. By contrast the SDV stratum ranked lowest in all indices.

In general, documentation showed that in the Mesa Central region precipitation-based indices indicated a definite influence on land use integration and the latitude gradient. Mahmood et al. [55] showed that land cover changes exert an important role in variations of atmospheric temperature, humidity, cloud cover, circulation, and precipitation. These changes have an impact ranging from local and regional scale to subcontinental and even global scales; they affect agriculture, deforestation/afforestation, desertification, and urbanization. From this perspective Betts et al. [56] have shown that land cover change can influence climate by modifying the physical properties of land surface.

4. Conclusions

This study, using multivariate statistical and geostatistical spatial analysis techniques, explored how different land use strata in the semiarid region of Mesa Central respond to ten precipitation-based indices. The very first conclusion of this study is that spatial autocorrelation, as observed through the analysis of the indices’ interpolated geostatistical layers, reduced the positional uncertainty in semivariograms and decreased with distance. The nugget value appearing in the PRCPTOT prompts the incorporation of other datasets, especially those that indirectly register sun insolation, soil humidity, and runoff. These datasets would help to measure the response of soil cover to precipitation and humid conditions.

The second conclusion is about the role of elevation as a dispersion factor affecting the proportion of observations from the mean; the dispersion was greater in the north region than in the central and south regions of the Mesa Central; but this dispersion is in contrast to the number of observations closest to the mean (>1.2 <2.4) and suggests further research initiatives.

The last conclusion is based on an interpretation of the semivariograms and predicted surface, along with the trend in indices. The spatial interpretation of resulted geostatistical layers deserves consideration as an important step in defining programs for biodiversity conservation and land use planning, in designing climate change mitigation programs, and in pursuing carbon sequestration and biomass production studies. This study, therefore, provides a foundation that can be confidently used in future research for exploring these important issues.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The authors wish to express their gratitude to the Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias (INIFAP) and Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE) for their financial and technical support. This gratitude is extended to Mario Salazar, Jaime Luévano, and Eulogio López for their assistance and technical comments. Special recognition goes to Wally Weerts for advising on grammar and English structure.