RELATIONSHIP BETWEEN SAMPLE DESIGN AND GEOMETRIC ANISOTROPY IN THE PREPARATION OF THEMATIC MAPS OF CHEMICAL SOIL ATTRIBUTES

Spatial variability depends on the sampling configuration and characteristics associated with the georeferenced phenomenon, such as geometric anisotropy. This study aimed to determine the influence of the sampling design on parameter estimation in an anisotropic geostatistical model and the spatial estimation of a georeferenced variable at unsampled locations. Datasets were simulated with geometric anisotropy, considering five values for


INTRODUCTION
In a spatial variability study, the spatial dependence structure described by the semivariance function is considered anisotropic when it depends on the distance and direction separating the locations observed (Maity & Sherman, 2012).For transitive models (with sill) that describe the semivariance function, there are detailed reports in the literature that describe the concept of anisotropy and types, the geostatistical tools that permit the identification of anisotropy and the description of the anisotropic spatial models and their parameters (Boisvert et al., 2009;Facas et al., 2010;Cressie, 2015).
One of the most common forms of anisotropy is the geometric anisotropy.This type of anisotropy occurs when a transitive model of semivariance function presents constant values of nugget and sill in all directions but different ranges.When geometric anisotropy is identified, it must be incorporated as an intrinsic feature of the process that describes the spatial dependence structure to improve the accuracy of the spatial estimat ion of the values of a georeferenced variable in unsampled locations, which results in the generation of maps that more accurately describe the s patial variability of the variable throughout the area under study (Guedes et al., 2008;2013).
In an anisotropic spatial model, the spatial interpolation by kriging method assigns greater weight to sampled values located along the direction of greater spatial continuity.In addition, knowledge of the anisotropy of the spatial dependence structure facilitates planning of sampling configurations for further studies of the spatial variability of this variable in the same area.Mucha & Blaszczyk (2012) used s imu lated data to determine that the influence of geometric anisotropy on spatial estimation depends on several factors associated with the spatial dependence structure, such as the distance over which the samples exh ibit spatial autocorrelat ion (range), the degree of intensity of spatial dependence (relative nugget effect) and the distance between the sampling units.According to the studies conducted by Guedes et al. (2013), for simulated datasets, it follows that fro m the anisotropic ratio of 2, the spatial estimates in unsampled locations differ depending on the incorporation or o mission of geometric an isotropy.
However, Mucha & Blaszczyk (2012) and Guedes et al. (2013) only considered systematic design (lattices) in their studies of spatial dependence.There are no studies that show the influence that the interaction between different samplings and the presence of anisotropy have on the estimat ion of the spatial dependence structure and on the estimat ion of values of the georeferenced variable in non-sampled locations.
Thus, considering simulated data and real data (soil chemical properties) with the presence of geometric anisotropy, the objectives of this work are the following: (a) to determine the implications of changing the sampling design for parameter estimat ion of the anisotropic spatial model; (b) to determine the implications of changing the sampling design for the spatial estimation of the georeferenced variable at unsampled locations; and (c) to evaluate different sampling designs to determine for which set of anisotropic parameters the isotropy assumption is unreasonable.

MATERIAL AND MET HODS
Sampling designs with 100 samp le points arranged in a regular area with a maximu m limit of x and y coordinates equal to 100 were considered.Three sampling configurations were simulated: lattice 10 x 10, random and lattice 9 × 9 with the addition of 19 nearby points (lattice plus close pairs).The latter sampling design, the simu lation study and its sample size were chosen because of the type of design used in the agricultural experiment for the study of spatial variability of soil properties performed in this study.
Datasets for each of these designs representing embodiments of mu ltivariate stochastic processes were simu lated for each of these configurations, assuming Gaussian stationary variables , with a linear space model expressed by ( is the deterministic term of the model, represented by a constant average ; ) ( i Z s corresponds to the observed values of the variable under study at the i th known location, denoted by the vector with ( is the stochastic error with , and the change in space between points separated by Euclidean distance , such that is determined by a covariance function , with ( Kempen et al., 2012;De Bastiani et al., 2015).
Through the covariance function, we find the  Opazo et al., 2012).
However, in this study, we have assumed that the simu lated georeferenced variable is stationary and exhibits a stochastic process with geometric anisotropy.Then, the covariance function is the function associated with the semivariance by the [eq.( 1)], in wh ich expresses the Euclidean distance between pairs of points in n sampled locations, considering a linear transformation at those locations (Soares, 2014). (1) The semivariance function [eq.( 1)] can be rewritten as a function of the Euclidian d istance between two points and , as shown in [eq.( 2)].
(2) where, is the highest spatial continuity angle in radians ( ) defined in the azimuth system; , and is the anisotropic ratio, where is the spatial dependency distance (range) in the direction of higher spatial continuity ( ) and is the spatial dependency distance (range) in the direction of lower spatial continuity ( ), assuming that the directions of higher and lower spatial continuity are orthogonal.
The simu lated data have an anisotropic semivariance function [eq.( 1)] with an exponential model with parameter vector , and is equal to 0°, following the system of azimuth direction.For each samp ling design, five datasets were simu lated, in wh ich took on each of the values of .Thus, to combine the sample values of the setting and anisotropic ratio, 1000 sets of data were simu lated.
In each set of simu lated data and in each of sampling designs, we estimated the parameter vector by the maximu m likelihood method, considering the following situations: (1) anisotropic model with the parameter vector , with fixed; and (2) isotropic model with the parameter vector .Considering the anisotropic model, to estimate asymptotic standard errors of , we used the inverse of the expected information matrix (Uribe-Opazo et al., 2012).The expected informat ion matrix to an anisotropic process is given by [eq.( 3)], where and , with for .
(3) Then, the first derivatives of with respect to and are given, respectively, by the following: , , and , with and , for .
The first derivatives of with respect to and for the exponencial and Gaussian covariance functions are given, respectively, by [eq.( 4)] and [eq.( 5)]: and ( 4) and ( 5) for , and .
The first derivatives of with respect to and for the Matérn family covariance function (Matérn, 1986) are given, respectively, by [eq.( 6)] and [eq.( 7)]: , where is the gamma function, ; is the modified Bessel function of the third type of order , with fixed; and .
Considering the estimated geostatistical models, the values of the georeferenced variables at unsampled locations were estimated in the presence and absence of geometric anisotropy by the ordinary Kriging method.
The sets of estimated values of the georeferenced variables at unsampled locations (considering the anisotropic and isotropic models) were co mpared for the similarity measures: overall accuracy ( ) Kappa concordance index ( ) and Tau concordance index ( ), which aim to compare the themat ic maps generated by the two interpolations, with the same rating fro m the error matrix, considering 10 intervals of values or classes (De Bastiani & Uribe-Opazo, 2012).Thus, using bootstrap resampling, confidence intervals were constructed for the mean similarity measures using the normal approximat ion (Efron & Tibshirani, 1993), considering the initial sample to perform the bootstrap resampling with 1000 values of each measure similarity obtained in the 1000 simulations (Fang & Wang, 2012;Rossoni et al. 2014).
In each set of simulated data and in each sampling design, tests of isotropy from Guan et al. (2004) (GSC) and Maity & Sherman (2012) (MS) were applied.For both tests, the null hypothesis that the georeferenced variable is isotropic was elaborated by considering a set of linear contrasts.Maity & Sherman (2012) elaborated that the null hypothesis is equivalent to and the alternative hypothesis to .However, Guan et al. (2004) defined the null hypothesis as and the alternative hypothesis as , for , where is a set of chosen spatial lags, and is a matrix of contrasts, whose dimension depends on the definition of .These tests propose to compare the values of covariance (Guan et al., 2004) and semivariance (Maity & Sherman, 2012) for pairs of lags described by .
The tests can be used when the sampling region is any convex subset in .The MS test was developed for sampling locations that are irregularly spaced.The covariance function is calculated by the kernel estimator, and the grid-based block bootstrap is used to calculate the test s tatistic.The GSC test was developed for gridded and non-gridded sampling locations.The semivariance is calculated by Matheron's estimator, and the test statistic is evaluated using a subsampling technique.
For each simu lation, these tests were constructed based on the set of lags .The choice of was made by considering short lags and the direction of the weakest and strongest spatial correlat ion in each set of simulated data.In the MS test, the block length was equal to 10 and 200 resamples for each set of simu lated data.In the GSC test, to preserve enough pairs of sampling locations in each subblock, a block length of 2 was chosen.
The geostatistical analysis, identification of anisotropy and the analysis of the influence of the incorporation of anisotropy on the spatial estimation were also performed on a set of real data fro m a co mmercial area of grain production in Cascavel, Paraná, located at approximately 24.95º of South latitude and 53.Simu lated datasets and geostatistical analyses were prepared in the software R (R Development Co re Team, 2016) using the geoR package (Ribeiro Jr. & Diggle, 2016).Tests of isotropy were conducted using the sm package (Bo wman & Azzalin i, 2015).

Analysis of simul ated data
For the simulated data and considering the anisotropic model for the semivariance function, the nugget and anisotropic ratio parameters had a greater influence on the change in the sampling design.We obtained the worst results in terms of quality of the estimat ion of these parameters for the lattice design: the estimated values were more distant fro m the nominal values of the parameters and (Fig. 1-a  and 1-b) and the estimated values of the standard errors of and were highest (Fig. 1-c and 1-d) co mpared with the random design and lattice plus close pairs.By contrast, imp roved estimates of the parameter were obtained by random design: the estimated values of were closest to its nominal value (Fig. 1-a) and the lowest estimated values of the standard errors of (Fig. 1-c).
The estimation of the influences the determination of the weights assigned to the samples in the spatial estimat ion of a georeferenced variable at unsampled locations, performed by Krig ing (Soares, 2014).Thus, the results for the estimation of the parameter demonstrated that the use of random design in the study of a georeferenced variable with geometric anisotropy permits a more efficient description of the spatial variability, as expressed by the map.These graphs illustrate that, for all studied sampling designs, the similarity measures decreased as the anisotropic ratio increased.For random sampling and an anisotropic ratio factor of 2, 94% of the simulat ions did not exhibit high similarity (Krippendorff, 1980;De Bastiani & Uribe-Opazo, 2012), indicating a relevant difference in the spatial estimat ion based on the incorporation or omission of the geometric anisotropy, fro m the anisotropic ratio of 2. For systematic samp ling and lattice plus close pairs, this relevant difference occurred in mo re than 90% of the simulat ions for anisotropic ratio values ( ) equal to 3.
However, all graphs (Fig. 2) demonstrated that, particularly for systematic sampling, there is high variability of these similarity measures and many small values (outliers) in appro ximately 10% of the simulat ions.For the random and lattice plus close pairs, outliers were observed in approximately 3.5% of the simulat ions.
These outliers occurred in simulat ions that produced an overestimat ion of the anisotropic ratio, in agreement with results obtained for the systematic sampling by Guedes et al. (2013).However, these authors emphasize that this overestimat ion of the anisotropic ratio occurred only for small sample sizes.
These conclusions are also supported by the confidence intervals constructed for the average of similarity measures using the bootstrap method and 95% confidence level (Dalposso et al., 2016).These confidence intervals are shown in Table 1 for all simu lations grouped according to the anisotropic ratio, the samp ling configuration and type of comparison performed.These confidence intervals indicated that, for random sampling fro m the anisotropic ratio of 2, the confidence intervals for the mean accuracy measurements did not contain the minimu m amounts required to be classified as high similarity for spatial estimates fro m the geostatistical models incorporating or o mitting geo metric anisotropy.
However, for the lattice plus close pairs and systematic sampling, the confidence intervals for the mean accuracy measurements did not contain the minimu m amounts required to be classified as high similarity for anisotropic ratio ( ) value of 3.These results provide evidence that there are differences in the classifications of generated thematic maps when considering the incorporation of anisotropy or not for the lattice plus close pairs and systematic sampling.
These differences can be exp lained by the fact that when the thematic maps are elaborated considering the presence of geometric anisotropy, these maps present subregions with greater continuity in the direct ion in which the anisotropy was identified (Guedes et al., 2008;2013).
Table 1 also showed that the similarity measures decrease significantly as the anisotropic ratio increases at 5% significance, considering all studied sampling designs.Similar results were described by Rossoni et al. (2014), who considered the mean square error of the spatial estimation of a georeferenced variable for unsampled locations as a measure of the impact of the incorporation of geometric anisotropy in spatial estimation in simu lations generated for random samp ling.TABLE 1. Confidence intervals for the average Overall Accuracy and Kappa and Tau concordance indexes obtained by the bootstrap method, comparing the spatial estimation using anisotropic and isotropic models, grouped according to the anisotropic ratio and the sampling design.

Design
Anisotropic ratio ( 1 (a)  2 3 Overall Accuracy ( ) Co mparing the sampling configurations for the confidence intervals of these similarity measures, mainly fro m the anisotropic ratio ( ) equal to 2, the similarity measures of the studied sampling configurations differed significantly at 5% p robability.The rando m sampling had the lowest values of these measures, whereas the systematic samp ling had the greatest values.
Table 2 shows the percentage of rejections of the nonparametric tests of isotropy at the 5% nominal level fro m 100 simulations.The method from Maity & Sherman (2012) considering the random sampling and lattice plus close pairs was applied.The method fro m Guan et al. (2004) considering the systematic sampling was also applied.For all samp ling designs, an inflated type I error (case isotropic with factor anisotropy equal to 1) was observed.
The power of tests using random and systematic sampling is relatively larger than the lattice plus close pairs.In addit ion, for all configurations examined, the emp irical power increased as the anisotropic ratio increased.However, accord ing to Guan et al. (2004) and Maity & Sherman (2012), emp irical powers tend be low when the sample size is small and/or the anisotropy is weak.
The real data grid was also used in our simu lations.The results not presented in this article showed that for the irregular area, a higher inflat ing of type I error (cas e isotropic), compared with the rectangular area, was observed.Table 3 shows the results of univariate geostatistical analysis with the estimated spatial parameters for the variables C and H + Al +3 for isotropic model and geometric anisotropic model.According to the cross-validation criteria (Monego et al., 2015) (Table 4), the exponential model was the best fit for C and H + Al +3 by the maximu m likelihood method, and the best estimates for the measure cross-validation were obtained for the anisotropic models.

ANALYS IS OF SOIL CHEMICAL PROPERTIES
Higher spatial continuity directions were established using the conventional directions adopted by Guedes et al. (2013) considering the azimuth system.In most models, there is a moderate (Cambardella et al., 1994; 25% < RNE  75%) spatial dependence, where RNE = 100 /( + ).Only the isotropic model estimated for the variable C exh ibited a strong spatial dependence (Cambardella et al., 1994; RNE  25%).The ratio is near 2, and the estimated values of the range and its standard error were higher in the isotropic model.Moreover, applying the MS nonparametric test (Table 3), the null hypothesis of isotropy evidently can be rejected (p-value < 0.05).: Ani ratio; Iso: isotropic; Ani: anisotropic, and Expo.: exponential model; * significant at the 5 percent level, by M S nonparametric test fort isotropy.
The results of the comparison of spatial estimates obtained by the isotropic and anisotropic models (Table 4) revealed that the anisotropic model resulted in lower values of the average Krig ing variance for both C and H + Al +3 .According to Yamamoto (2000), the kriging variance is used as a quality indicator of spatial estimation.This result indicates that the anisotropic model produced the best efficiency and spatial estimation.
In addition, the values of the similarity measures indicate a relevant difference in the spatial estimates of the maps generated by the isotropic and anisotropic models (Kripendorff, 1980;De Bastiani & Uribe-Opazo, 2012).The results for the similarity measures in this work are similar to conclusions presented by Guedes et al. (2013) considering simulated and real data from a systematic sampling.
The differences in the spatial estimat ion of the variables under study upon the incorporation of geometric anisotropy in the spatial dependence structure can also be displayed in thematic maps (Fig. 3-(c) to 3-(f)).In thematic maps constructed from the anisotropic models (Fig. 3-d and 3-f), a greater spatial continuity of subregions was observed in the higher spatial continuity directions compared to the thematic maps generated by isotropic models (Fig. 3-c and 3-e).Similar results were described by Guedes et al. (2013), who estimated chemical soil properties for unsamp led locations, considering stratified systematic unaligned sampling.TABLE 4. Adjusted spatial models with their respective measurements obtained by cross -validation and associated measures the spatial estimation and the similarity in the co mparison between anisotropic and isotropic models , its spatial estimation. : average kriging variance; M E: mean error; S ME : standard deviation of the mean error; SM E: mean standardized error; S SME : standard deviation of the standardized error; AE: absolute error; Iso: isotropic; Ani: anisotropic.

CONCLUS IONS
The analyses of the simulat ions demonstrated that the systematic sampling exh ibited the worst performance in relation to the quality of estimat ion of the geostatistical anisotropic model co mpared with the random samp ling and lattice plus close pairs.This lower performance generally resulted in overestimated values of the nugget and anisotropic ratio, the highest estimated values of the nugget, range and anisotropic ratio, and the highest estimated values of the standard errors of the nugget, range and anisotropic ratio.
Moreover, the results of the simulations demonstrated that the random sampling and lattice p lus close pairs are more appropriate for the incorporation of geometric anisotropy with respect to spatial estimat ion for georeferenced variables at unsampled locations.The similarity measures and the nonparametric tests demonstrated that random samp ling, lattice p lus close pairs and lattice sampling have a relevant influence on the incorporation of geometric an isotropy in most simu lations fro m an isotropic ratio values of 2.
Most of the similarity measures differed significantly for the different sampling designs.These confidence intervals did not overlap at 95% confidence when the anisotropic ratio was changed.
The similarity measures for the chemical attributes of soil sampled in an irregular area with a lattice plus close pairs sampling design indicated a relevant difference in the spatial estimat ion depending on the incorporation or omission of geometric anisotropy in the geostatistical model.Moreover, the smallest values of the average kriging variance and improved detail of the subregions in the thematic maps were observed in the spatial estimation that incorporated geometric anisotropy.These results demonstrate that the incorporation of geometric anisotropy produces more reliable estimates that best indicate the spatial continuity of the studied attributes.
37º of West longitude and at an average height of 650 m above sea level.The dominant soil type is the Rhodic Hapludox, which has a clay-like texture.The data refer to the crop year 2010/ 2011 and an area of 167.35 ha of experiments conducted by researchers from the Spatial Statistics Laboratory of the State University of West Paraná -Brazil.The sampling configuration used in this study was a lattice p lus close pairs, with 102 sampling points.All samples were georeferenced and localized with the aid of a signal receiving apparatus with a global positioning system (GPS) Geoexp lore 3 set up for the Universal Transverse Mercator (UTM ) coordinate system.The area was planted with soybeans, and the data used in this study were related to soil chemical properties that exhibited geometric anisotropy in the geostatistical analysis: carbon (C -) and potential acidity (H + Al +3 -).The data sets were obtained by performing routine chemical analysis in the soil analysis laboratory of Cooperativa Central de Pesquisa Agrícola Ltda.(COODETEC) of representative samples of each plot of approximately 50 grams, which were obtained by mixing four rep licates of the parcel.H + Al +3 were obtained by the Shoemaker, Mac Lean and Pratt (SM P) buffer method.

FIGURE 1 .
FIGURE 1. Bo xplots for: estimated values of (a) nugget effect ( ) and (b) an isotropic ratio ( ); and estimated values of the standard error (SE) of (c) nugget effect () and (d) anisotropic ratio ( ) for an isotropic geostatistical model.

Figure 2 p
Figure2p resents boxplot graphs of the similarity measures that compare the spatial estimation performed using the anisotropic and isotropic exponential models, both with estimated parameters.The dashed lines in these graphs correspond to ranges of values for these accuracy measures in which the spatial estimates were of similar intensity.For overall accuracy, the minimu m level of (a) Isotropic cas;.M S: method from M aity & Sherman (2012); GSC: method fromGuan et al. (2004).

Fig. 3 -FIGURE 3 .
Fig. 3-(a) and Fig. 3-(b) illustrate the values of variables C and H + Al +3 for the sampled points in the study area, using quartile classificat ion.Groups of points with the same classificat ion extended in the direction of 90° (azimuth system), indicating the presence of geometric anisotropy in that direction.

TABLE 2 .
Percentage of rejections at 5% no minal level of the nonparametric tests of isotropy.

TABLE 3 .
Adjusted spatial models and estimated values of parameters by ML.The standard deviations of the estimates are in parentheses.