Exploratory Analysis

The quality control applied to the data consisted of evaluating the following points:

Spatial Consistency - The coordinates of eight stations were updated due to inconsistency in their location: five in Coahuila (5016, 5032, 5035, 5048 and 5049) and three in Nuevo León (19012, 19054 and 19069).

Format - Improbable dates and repeated records were excluded, as an example in NOAA records, every month of every year has a 31-day format.

Completeness - For the annual data integration, 95 stations were selected that contained at least 25 years of records, considering the criterion of excluding the months with more than five missing daily records or more than three continuous daily records.

The basic statistics for the variables PMA Spring and Summer are described in Table 1. The three precipitation variables have a positive asymmetric distribution, so it is necessary to evaluate whether the "proportional effect" is presented in these cases, which is a particular form of heteroscedasticity (the variability of the data changes throughout the study area), particularly for positive asymmetric distributions, the local variance increases as its local mean increases (^{Goovaerts, 1997}). This proportional effect makes it impossible to interpret experimental variograms (^{Grimes & Pardo-Igúzquiza, 2010}).

To identify how changes in local variance are related to changes in the local average, it is common to use scatter plots, calculated from mobile windows (^{Goovaerts, 1997}). Statistics for mobile windows were calculated with an area of 100 km², with an overlap of 25 km², and 18 were selected, including at least 10 data. The overlap of the windows implies that some data were used repeatedly for the calculation of local means and variances, which is not relevant at the moment in which we only want to establish whether or not the proportional effect exists (^{Isaaks & Srivastava, 1989}). The result obtained for PMA and Summer shows that variances and local averages have low correlation ( ρrank=0.31), so that homoscedasticity can be assumed, while Spring has a proportional effect due to a higher correlation (ρrank=0.62). Figure 3 shows the scatter plots.

Figure 3 Scatter plots for 18 mobile windows with at least 10 data per window. 

As is known, the weather stations are located in places where certain characteristics are available in terms of operation and registration, so there are areas with higher concentrations (clustering), particularly when combining the proportional effect and the grouping of the sample data. This led to conflicts with the interpretation of the experimental variograms. One way to know the grouping and the proportional effect is by graphing the local means and variances as a function of distance (^{Goovaerts, 1997}).

For the case of the three variables, fluctuation around a unit value is observed, concluding that the variances and the means are not significantly affected by distance, so there is no need to perform any process to consider the effect of grouping, for example, using only regularly spaced data (^{Goovaerts, 1997}). Figure 4 presents the results obtained for PMA.

Figure 4 Means and variances as a function of distance (lags). Both statistics are normalized by their corresponding mean and global variance. 

Structural analysis

Considering that asymmetry and atypical data directly affect the modeling of variograms, a logarithmic transformation was applied to reduce the degree of asymmetry and analyze the behavior of atypical data. The logarithmic transformation used for Spring did not improve the distribution, so its application was not considered. Figure 5 shows the experimental variograms for the variables with moderate asymmetry ( PMA and Summer ). Their logarithmic transformation were included in the same graph after re-scaling the variance. The logarithmic variograms seem less erratic, with less nugget for PMA ; therefore, this variable is modeled without transformation, while Summer justifies the transformation due to the notable difference between the variograms.

The anisotropy analysis was performed by inspecting the experimental variograms in different directions: 0° (NS), 45° (NE-SW), 90° (EW) and 135° (NW-SE) with an angular tolerance of ± 22.5° (^{Wackernagel, 2003}). Since a portion of the Sierra Madre Oriental is located in the study area, and it influences precipitation patterns, variograms were constructed in the NW20° direction, in association with the regional structure of that mountain range, as well as the perpendicular direction to it (NE70°). Additionally the omnidirectional variogram was constructed (0° direction and angular tolerance of 90°). Gstat software (^{Pebesma & Wesseling, 1998}) and the RGEOESTAD (^{Díaz, Hernández, & Mendez, 2012}) were used to perform this analysis.

Figure 5 Omnidirectional experimental variograms for PMA and Summer. The logarithms variogram was re-scaled to the variance from the original data. 

According to Figure 6, given the inflection points shown in the experimental variograms, two spatial variation scales in the NE70° direction (perpendicular to the structure of the Sierra Madre Oriental) are seen: one with a range of 75 km and another of 125 km approximately. At a distance greater than 125 km, the structure of the variograms becomes erratic in that direction due to the decrease in the number of pairs that contribute to the values of the variances. For the NW20° direction, the variograms show greater spatial continuity (less variance) and follow a pattern similar to the omnidirectional variogram.

Figure 6 PMA, Spring and Log Summer variograms in NE70°, NW20° directions with angular tolerance of 22.5 (dashed lines), and the omnidirectional variogram (0°) with angular tolerance of 90° (continuous line). The numeric labels show the number of pairs. 

Another way to evaluate anisotropy is through the elaboration of a variographic map (^{Isaaks & Srivastava, 1989}), in which the values of the experimental variogram are plotted and the center of the map corresponds to the origin of the variogram. Figure 7 shows a preferential direction towards the NW-SE for Log Summer and PMA, roughly aligned with that presented by the Sierra Madre Oriental. For the case of Spring, this preferential direction appears to be less.

Figure 7 Variographic maps, PMA, Log Summer and Spring, values standardized by the sample variance. 

Taking into account the results of the variographic maps, the anisotropic models were used for the variables PMA and Log Summer .

Once the directions of the variograms were defined, spherical theoretical models were adjusted using least squares technique. This procedure was carried out with the RGEOESTAD (^{Díaz et al., 2012}) and gstat (^{Pebesma & Wesseling, 1998}) software, from which variographic parameters were obtained (nugget, sill and range). Table 2 shows the summary and Figure 8 illustrates the adjusted spherical models, including the omnidirectional. The models adjusted for the NE70° direction are not shown in this paper.

Figure 8 Directional variograms adjusted for PMA , Spring and Log Summer : a) NW20° direction, b) omnidirectional. 

All variograms for the NW20° direction exhibit a linear behavior at the origin, as well as small nugget effects, with nugget/sill proportions from 16% to 23%, and represent an origin discontinuity that is attributable to measurement errors and variation at lower distances in the sampling interval (28.4 km).

The variograms in the NE70° direction were modeled to obtain the anisotropy factor λ, defined as the ratio of the lower range (NE70°) to the larger range (NW20°).

Cross-validation

The cross-validation technique consists of taking an observation and estimating its value with the remaining observations. This is done repeatedly until all observations are concluded.

The errors are defined as the difference between the observed and estimated values, and are calculated using the following expression:

Where ε are the errors, Z the observed value and Z* the estimated value.

There are different indicators to measure the quality of the estimate globally. In this work, the following were applied:

Table 3 shows the validation results.

According to RMSE from cross-validation, the anisotropic model improved the PMA and Spring estimations. For the Summer estimation, the isotropic model was better, according to the results.

Estimation and Simulation Methods

The estimated rainfall and error maps are shown in Figures 9, 10 and 11. These maps were generated using the ArcMap software with geostatistical analyst extension (^{ESRI, 2016}), using previously adjusted variographic models and applying the ordinary kriging and Gaussian sequential simulations techniques. These techniques are discussed and widely analyzed in the works of some authors such as: ^{Chilès and Delfiner (1999)}, ^{Goovaerts (1997)}, ^{Wackernagel (2003)} and ^{Isaaks and Srivastava (1989)}.

Isotropic variographic models were used for the Gaussian sequential simulation technique, applying an anamorphosis transformation to the data. The open source program gstat (^{Pebesma & Wesseling, 1998}) was used to perform 100 simulations for each of the variables and to generate the maps. An average of these simulations was considered, using the ArcMap geostatistical analyst extension (^{ESRI, 2016}).

The anisotropic model better reproduced the spatial distribution of PMA rainfall data, emphasizing the preferential direction to the Sierra Madre Oriental. In the northern part of Nuevo León, there are very similar patterns among the three models (isotropic, anisotropic and simulations), and as is to be expected, the simulations reproduced the exact same extreme data values.

Figure 9 Nuevo León maps: estimation, simulation and estimate error for PMA , in the period 1930-2014. 

Figure 10 Nuevo León maps: estimation, simulation and estimate error for Spring , in the period 1930-2014. 

Figure 11 Nuevo León maps: estimation, simulation and estimate error for Summer , in the period 1930-2014.