Directional Differences in Thematic Maps of Soil Chemical Attributes with Geometric Anisotropy

: In the study of the spatial variability of soil chemical attributes, the process is considered anisotropic when the spatial dependence structure differs in relation to the direction. Anisotropy is a characteristic that inﬂuences the accuracy of the thematic maps that represent the spatial variability of the phenomenon. Therefore, the linear anisotropic Gaussian spatial model is important for spatial data that present anisotropy, and incorporating this as an intrinsic characteristic of the process that describes the spatial dependence structure improves the accuracy of the spatial estimation of the values of a georeferenced variable in unsampled locations. This work aimed at quantifying the directional differences existing in the thematic map of georeferenced variables when incorporating or not incorporating anisotropy into the spatial dependence structure through directional spatial autocorrelation. For simulated data and soil chemical properties (carbon, calcium and potassium), the Moran directional index was calculated, considering the predicted values at unsampled locations, and taking into account estimated isotropic and anisotropic geostatistical models. The directional spatial autocorrelation was effective in evidencing the directional difference between thematic maps elaborated with estimated isotropic and anisotropic geostatistical models. This measure evidenced the existence of an elliptical format of the subregions presented by thematic maps in the direction of anisotropy that indicated a greater spatial continuity for greater distances between pairs of points.


Introduction
Geostatistical modeling analyzes the relationship between the values of a georeferenced variable observed at different locations and also considers that each sample unit carries information about its neighborhood.Thus, the representation of the spatial continuity of a georeferenced variable in the study area is performed with geostatistical techniques and information collected in a sample.One of the geostatistical methodologies consists in the construction of directional semivariograms, which allows for a preliminary analysis of the existence and type of anisotropy in the phenomenon studied [1].
In this context, the spatial dependence structure of a regionalized variable is said to be isotropic when the spatial continuity presents similar behaviors for all directions.Therefore, any function used to describe the spatial dependence structure, such as the semivariance (γ(h)) and covariance (C(h)) functions, depends only on the modulus (h) of the vector distance h between two locations, where h = h is the Euclidean distance [2][3][4].
Consequently, the estimated values of the parameters that define the spatial dependence structure (practical range, nugget effect and sill) are similar in all directions.However, if the spatial dependence structure differs in relation to the direction, then the process that describes the spatial dependence structure is said to be anisotropic [5].
If the presence of anisotropy is identified in the phenomenon studied, then this characteristic is incorporated in the covariance structure, since the correct specification of the spatial dependence structure can influence the spatial estimation of values of the georeferenced variable in non-sampled locations, as well as its accuracy.For example, when the spatial dependence structure is anisotropic, the spatial prediction assigns greater kriging weights to the samples in the direction of greater spatial continuity, and both the variance and the error of spatial prediction that are obtained in kriging are smaller when the anisotropic model in the kriging rather than the isotropic model is considered [6,7].
Studies developed by Facas et al. [8] and Guedes et al. [5,9], in simulated and real data sets (both with geometric anisotropy), showed relevant differences between the thematic maps elaborated considering or not considering the presence of the anisotropy.The overall difference between these thematic maps, ranging from pixel to pixel, varied from 5% to 35%, and this variation was influenced by the size and sampling design, as well as by the anisotropy factor.
These authors also concluded that the main difference presented in these thematic maps was the formation of elliptic subregions in the direction of anisotropy when an anisotropic model was considered in the analysis of spatial variability.However, this result was only observed through a visual analysis of these thematic maps as well as omnidirectional measures (omnidirectional accuracy indices) to quantify a global difference between the thematic maps.Thus, the objective of this work was to identify, in a quantitative way, the directional difference that exists in the thematic maps when the anisotropy is incorporated in the geostatistical model, using Moran's directional spatial autocorrelation index (I(d)).

Simulation Study
Four trials, each with 100 simulations, from the Monte Carlo experiment were performed, where each simulated data set represents a set of stochastic process realizations {Z(s i ), s i ∈ S}, where s i = (x i , y i ) T is a vector that represents the i-th location in the study area (i = 1, . .., n), such that S ⊂ R 2 and R 2 is the two-dimensional Euclidean space [10].This simulation study reproduces possible real data sets and improves the theoretical and practical knowledge about the anisotropy.
For each simulated data set, an irregular sampling design with 100 points (n = 100) was considered, with ordinates ranging from 0 to 1 on the X and Y and where each Z(s i ), com i = 1, . .., n, represents a georeferenced variable.
This georeferenced variable was expressed by the isotropic linear Gaussian spatial model Z(s i ) = µ(s i ) + (s i ) [11], where the deterministic term (µ(s i ) = µ), (s i ) is the stochastic term for i = 1,..., n, such that both depend on the spatial location in which is observed Z(s i ), such that E[ (s i )] = 0 and the variation between points in space is determined by a covariance function C(s i , s j ) = cov[ (s i ), (s j )].Also, we have assumed that Z = (Z(s 1 ), . .., Z(s n )) T has a Gaussian distribution n-variate, with vector mean equal to µ1, n × 1, and the covariance matrix ∑ = [(σ ij )], n × n, with σ ij = C(s i , s j ), i, j = 1, . .., n, with particular parametric form given by ∑ = ϕ 1 I n + ϕ 2 R, where I n is an identity matrix, n × n, R = R(ϕ 3 ) = [(r ij )] is a symmetric matrix, n × n, with diagonal elements r ii = 1 and r ij = ϕ 2 −1 σ ij for i = j = 1, . .., n.Thus, the covariance function is the function associated with semivariance by γ(h ij ) = C(0) -C(h ij ) for many isotropic and stationary Gaussian processes, where h ij = s is j is the value of Euclidian distance between the points s i and s j [4].Moreover, this covariance matrix has the following parameters: ϕ 3 is a function of the range (a > 0), ϕ 1 is the nugget effect and ϕ 2 is the partial sill.
In the first trial, we are assuming that µ = 200, and an exponential function correlation [12] composed of the following parameters that define the structure of spatial dependence: ϕ 1 = 0, ϕ 2 = 1, and a = 0.6.
In the order trials (2 • , 3 • , and 4 • ), this georeferenced variable was expressed by the anisotropic linear Gaussian spatial model [9], with the following technical characteristics: the same correlation function and one change in the linear Gaussian spatial model described above, where Ah ij expresses the Euclidean distance between pairs of points in n sampled locations in the semivariance function, i.e., γ( Ah ij ) = C( 0 ) − C( Ah ij ), considering a linear transformation at those locations expressed by the product matricial Ah ij , with h ij = s i − s j , where the transformation matrix is equal to A = cosβ − senβ /F a senβ cosβ /F a , with β as the highest spatial continuity angle in π radians (0 ≤ β ≤ π) defined in the azimuth system and F a = is the anisotropic ratio (F a > 1), where α β is the spatial dependency distance (range) in the direction of higher spatial continuity (β) and α β + π 2 is the spatial dependency distance (range) in the direction of lower spatial continuity (β + π 2 ) [5,7,9].In the order trials (2 • , 3 • e 4 • ), we are assuming that the angle of greater spatial continuity is equal to 90 • and the anisotropic ratio is (F a = 2, 3, and 4).Moreover, the first trial (isotropic) is a particular case of the anisotropic model (F a = 1) [4].
For each simulation, the estimation of the anisotropic Gaussian space linear model was carried out using the maximum likelihood method obtained by the following parameters: mean (µ), nugget effect (ϕ 1 ), partial sill (ϕ 2 ), practical range (a) and anisotropic ratio (F a ) [9,11].
Then, ordinary kriging was used in the spatial estimation of each simulated data set at unsampled locations.In order to compare the differences in the spatial estimation when considering the presence of geometric anisotropy in the application of the linear Gaussian spatial model, kriging was used for each simulated data set with an estimated isotropic linear Gaussian spatial model.
For each simulated data set, the directional comparison between the two maps using the predicted values (one considering the isotropic model and the other considering the anisotropic model) was made.The directional comparison for I(d) (Equation ( 1)), proposed by Rosenberg [13] was measured.For this calculation, five distance classes (0.15, 0.30, 0.45, 0.60, and 0.75) and two directions (90 • and 0 • ) were considered.Each distance was chosen to guarantee a relevant number of points pairs, greater than or equal to 40, for the calculation of I(d) [13].
where n is the number of unsampled points that were considered in the kriging; m l and m k are the predicted values in the point s l and s k (l, k = 1, . . ., n); m is the mean of predicted values for all points; W(d) = ∑ n l=1 l =k ∑ n k=1 k =l w lk (d) is the sum of the elements of a spatial weights matrix W (n × n) elaborated with the distance class d; and w lk (d) = w lk is the element at the l-th row and k-th column of the weights matrix W , expressed by: where for which d lk is the Euclidean distance between points s l and s k (l, k = 1, . . ., n and l = k); θ is the direction (in radians) of interest for the calculation of the Moran directional index; and α lk is the angular direction (in radians) between points s l and s k (i.e., between a line parallel to the X axis and the line formed by the points s l and s k ), such that 0 ≤ α lk < π and α lk = α kl , calculated by: α lk = atan where s l = (x l , y l ) T and s k = (x k , y k ) T are the coordinates of the points s l and s k (l, k = 1, . . ., n and l = k), respectively.Equation ( 2) demonstrates that Rosenberg [13] defined the Moran directional considering, in this index, the inclusion of a weight (or penalty) factor weights matrix W , expressed by cos 2 (α lk − θ), with l, k = 1, . . ., n and l = k.In this way, w lk = w lk , if the direction between points s l and s k is equal to the direction of interest (i.e., α lk = θ).However, the greater the distance between the direction of interest (θ) and the direction between points s l and s k (α lk ), the smaller the value of w lk (relative to the value of w lk ).
High values of I(d) indicate that there is a spatial autocorrelation of the georeferenced variable [14].Thus, the influence of distance class on the magnitude of the index can be explained by the basic geostatistical concept that closer data look more like more distant data [1].Thus, the analysis of the variation of I(d) as a function of the distance class (d) serves as an indication of the range of spatial dependence of the georeferenced variable for a specific direction.
Also, the significance of this index was evaluated by the pseudo-significance test, with a 5% probability [14].

Soil Chemical Property Study
The agricultural data were collected in a commercial area of grain production in Cascavel, Paraná, Brazil (Figure 1), with a total area of 167.35 ha.The area is located at approximately 24.95 • of South and 53.57• of West, with a mean altitude of 650 m above sea level.The soil is classified as typical Dystroferric Red Latosol, with a clayey texture.The region's climate is classified as mesothermal and super-humid temperate, climate type Cfa (Köeppen classification system), and the mean annual temperature is 21 • C. The sampling points were georeferenced and localized using a GNSS receiver (GeoExplorer, Trimble Navigation Limited, Sunnyvale, CA, USA) in a Datum WGS84 coordinate reference system, UTM (Universal Transverse Mercator) projection.The lattice plus close pairs sampling with 102 points was sampled [1,15].This sampling consists of a regular grid, with a minimum distance between points of 141 m; 19 points were randomly added in this regular grid, such that the smallest distance between the points added and the point of this regular grid is 75 and 50 m (Figure 1).The sample data were obtained through a routine chemical analysis in the soil analysis laboratory of COODETEC (Cooperativa Central de Pesquisa Agrícola) of representative samples of each plot weighing approximately 500 g and collected at each demarcated point (Figure 1).The following soil chemical properties with geometric anisotropy were considered: carbon content (C) (g dm −3 ), calcium content (Ca) (cmolc dm −3 ) and potassium content (K) (cmolc dm −3 ).The chemical analyses were performed using the Walkley-Black method [16].
Descriptive and geostatistical analyses were performed for each soil chemical property to verify the existence of directional trends and spatial dependence.Directional trends represent a linear association between the respective values of the soil chemical properties with the coordinates of the X or Y axis, and were assessed by Pearson's linear correlation coefficient (r(x), r(y)), in which values above 0.30 in a module indicate a directional trend [17].Spatial dependence was assessed by the spatial dependence index (SDI), being classified as weak when SDI ≤ 9%, moderate when 9% < SDI ≤ 20% and strong when SDI ≥ 20% [4,18].Still, the directional comparison analyses described above were also performed.For the calculation of (Equation ( 1)), five distance classes (150, 300, 450, 600 and 750 m) and two directions (the direction of greater spatial continuity (θ) and the direction orthogonal to it) were considered.
Moreover, the global comparison between the two maps using the predicted values (one considering the isotropic model and the other considering the anisotropic model) were made.The measurements for global comparison were performed by the Global Accuracy and the Kappa and Tau concordance indices [19,20].

Computational Resources
Simulated data sets, statistical and geostatistical analyses and the calculation of the Moran directional index were made in the geoR package [21] of software R version 3.5.1 [22].

Simulation Study
For the simulated data, the results obtained for the values of the Moran directional index (I(d)) calculated for the estimated values of the georeferenced variable in unsampled locations are presented in Table 1 and Figure 2.For the isotropic model (F a = 1) and for all distance classes, similar values of I(d) were observed in both directions considered (0 • and 90 • ) (Figure 2).In all directions, most of the simulations presented significant values of I(d) (Table 1).In addition, as the distance class increased, there was a similar trend of decreasing of the values of I(d) in both directions (Figure 2), as well as the percentage of significant values of this index (Table 1).In this way, under the presence of isotropy the results indicate that the pairs of unsampled locations considered in kriging presented the same degree of spatial similarity in both directions (0 • and 90 • ).Thus, in the 0 • and 90 • directions there was a similar description of the spatial continuity of the variable georeferenced in the thematic map.
For all other anisotropic ratios considered in the simulations (F a = {2, 3, 4}), the higher the distance class, the higher the decay rate and percentage of significant values of I(d) in the 0 • direction compared with the 90 • direction (Figure 2 and Table 1).
Even in the largest distance class (d = 0.75), the 90 • direction presented the higher values of I(d), as well as the higher percentage of I(d) significant values.These relevant differences occurred in almost all the simulations with the largest values of the anisotropic ratio (F a = {3, 4}).
These results suggest that the 90 • direction presented spatially similar predicted values within a radius of distance between their respective locations equal to 0.75, while in the 0 • direction, the distance radius separating pairs of predicted values that are considered similar is less than 0.60.

Soil Chemical Property Study
For the soybean crop, the difference between maximum and minimum samples (Table 2) indicated that each soil chemical property is classified from medium to very high level [23].According to Pimentel Gomes [24], the sample values of carbon content (C), calcium content (Ca) and potassium content (K) presented medium (10% ≤ CV ≤ 20%), high (20% ≤ CV ≤ 30%) and very high (CV > 30%) dispersion in relation to their corresponding average, respectively (Table 2).In addition, the soil chemical properties, when associated with the x and y coordinates, presented the estimated value of Pearson's linear correlation coefficient as lower than 0.30 (r(x), r(y)) (Table 2).This result suggests an absence of a linear trend of these properties with the X and Y axes [17].However, as for some properties the p-value is lower than 0.05 (Table 2), the linear correlation can be considered statistically significant, and thus the coordinate x has a significant linear association with the carbon content (C) and calcium content (Ca), while the y coordinate has a significant linear association with the potassium content (K) (regardless of the correlation coefficient value) (Table 2).According to the Postplot graph (Figure 3), the soil chemical attributes presented discrepant points that are located in the central region of the study area for carbon and calcium contents (Figure 3a,b), whereas for the potassium content, these are located in the southern region of the study area (Figure 3c).
In addition, for each soil chemical property, the Postplot graph showed the formation of clusters of sample points with similar values, mainly with a higher concentration of points in the 90 • direction (Figure 3).This result indicates that the soil chemical properties presented a directional tendency regarding spatial continuity, that is, geometric anisotropy [5].
Table 3 presents the results of the univariate geostatistical analysis with the linear anisotropic Gaussian spatial model estimates by maximum likelihood-ML, which describes the dependence structure of the soil chemical properties in the study area.For all soil chemical properties, the cross-validation criteria, Akaike information (AIC) and Bayesian information (BIC) showed that the covariance function described by the Gaussian model was the best fit [25,26].The estimated covariance matrix (Σ) is expressed according to Equations ( 3)-( 5), respectively, considering the attributes C, Ca and K.
where I is an identity matrix, n × n, with n = 102; Âh ij is the value of Euclidian distance between the points s i = (x i , y i ) T and s j = x j , y j T [20], considering a linear transformation The estimated values of the practical range and the spatial dependence index (SDI) indicated the existence of spatial dependence for all soil chemical properties, with a radius of spatial dependence from 149.57 to 371.89 m, and spatial dependence classified as weak for the soil carbon and calcium contents (SDI ≤ 9%) and moderate for the soil potassium content (9% < SDI ≤ 20%) [4,18].
For all the soil chemical properties, the presence of geometric anisotropy in the 90 • direction was observed, and the estimated value of the anisotropic ratio indicated a spatial dependence radius of 3.754 to 6.078 times greater in the 90 • direction than in its orthogonal direction (0 • ).
Figure 4 presents the thematic maps of soil chemical properties generated by kriging, also considering anisotropic and isotropic models.The thematic maps generated by the anisotropic model presented subregions with an elliptical format, whereas the same subregions presented a more circular format in the maps generated by the isotropic model.Comparing the thematic maps of soil chemical attributes considering the isotropic and anisotropic model, an estimated ÔA value of less than 85% was found, as well as K and T of less than 67%, which indicates that the maps are not similar, that is, that the maps prepared considering both models are not similar in terms of the distribution of attribute content in the area under study ( ÔA < 0.85 e K, T < 0.67; [19,20]) (Table 4).The similarity measures describe differences in spatial estimation globally, that is, in all directions.Conversely, Table 5 presents the values of the directional Moran I (I(d)), calculated for the predicted values by kriging, considering or not considering the presence of anisotropy in the carbon, calcium and potassium contents in the soil.Thus, for each soil chemical property and a direction of interest, I(d) describes whether or not there are directional differences between thematic maps elaborated with anisotropic and isotropic models.For all soil chemical attributes and geostatistical models (anisotropic and isotropic), almost all of the estimated values of I(d) are significant at 5% probability (Table 5).Also, for all distance classes, higher values of I(d) were obtained in the direction of greater spatial continuity (90 • ) compared to the estimated values of I(d) in the 0 • direction (Table 5).
Indeed, the estimated values of I(d) showed that, in the 90 • direction, there was spatial autocorrelation between points spaced at 750 m, and, in the 0 • direction, there was a weak spatial correlation from 450 m for the carbon and calcium content in the soil, with estimated values of I(d) closer to zero (Table 5).
These results indicate that the presence of anisotropy in the carbon, calcium and potassium contents influenced the shape of the subregions (showing a greater spatial continuity towards the anisotropic direction), even though the spatial prediction was performed with the isotropic model.
Moreover, there was a greater difference between the estimated values of I(d) for the 0 • and 90 • directions in the spatial prediction performed with the anisotropic model, with higher values of I(d) in the 90 • direction (Table 5).Thus, from the distance class of 450 m, the estimated values of I(d) in the 90 • direction were greater than 1.65 to 7.83 times at the estimated values of I(d) in the 0 • direction (Table 5).
Therefore, the thematic maps elaborated with the anisotropic model presented a greater spatial continuity in the 90 • direction compared with those elaborated with the isotropic model, evidencing the existence of directional differences between thematic maps elaborated when considering or not considering the anisotropy in the geostatistical model.

Discussion
The simulations with F a = 1 showed that the 90 • direction is the direction with the highest directional spatial autocorrelation, which indicates the higher similar estimated points density in this direction [13,28,29].Thus, under the presence of anisotropy, there was a significant difference between the direction that defines the presence of geometric anisotropy and its orthogonal direction, in relation to the spatial continuity of the estimated values of the georeferenced variable in unsampled locations.
There was an inverse trend between the values of the distance classes and I(d).However, the value of the coefficient in the largest distance class is often unreliable due to the few pairs of points in this class [13].In this way, the same problem that occurs with the directional semivariogram, when compared with the omnidirectional semivariogram, in relation to the existence of a smaller number of pairs of points at each distance [8], also occurs when we compare the calculation of directional and omnidirectional spatial autocorrelation [29].It may not be possible to analyze the directional spatial correlation for some data sets with a smaller sample size.The directional spatial autocorrelation could be calculated for various classes of distances/directions and the results are plotted in a diagram known as a Windrose, which is a complementary method for identifying the anisotropy [29].
Considering the soil chemical property study, the analyses indicated relevant differences in the thematic maps, when considering the incorporation or not of the geometric anisotropy, mainly in the direction of greater spatial continuity (90 • ).These differences and the shape of the subregions were also observed in thematic maps generated in studies conducted by Guedes et al. [5], even with data collected in another sampling design, and occur because it is assumed that the ellipsoid representing the anisotropy of the estimated property is centered on each node to be estimated [30].
A low value of the similarity measures ( ÔA < 0.85 e K, T < 0.67; [19,20] was obtained for each soil chemical property (Table 4).According to Richetti et al. [31], this low similarity represents a relevant difference between the generated maps, when the geometric anisotropy was or was not incorporated, in relation to the percentage of the area classified in each class.Thus, this result shows that the incorporation of the anisotropy, when it exists, alters the spatial distribution of these soil chemical properties in the study area.
As reported by Guedes et al. [9], these relevant differences in the spatial estimation of georeferenced variables presented by the similarity measures are justified by the high estimated values of the anisotropic ratio factor, since the values of the anisotropy factor were higher than 2.
These results corroborate other results obtained in the literature for simulated and real data, from the point of view of a traditional anisotropy analysis (modeling, visual analysis of the thematic map) [5,[7][8][9].But the main contribution of this article to research already carried out on anisotropy is the possibility of using directional spatial autocorrelation as a metric to highlight the existence of directional differences in the thematic maps regarding the subregions format, elaborated when considering or not considering the anisotropy in the geostatistical model.
Furthermore, the linear anisotropic Gaussian spatial model is important and suitable for the simulated and real data of this work, because when the data present anisotropy, it must be incorporated as an intrinsic characteristic of the process that describes the spatial dependency structure in order to improve the accuracy of the spatial estimation of the values of a georeferenced variable in unsampled locations [9], and this results in the generation of maps that describe with greater accuracy the spatial variability of the variable throughout the area under study [5,7].

Conclusions
Under the presence of geometric anisotropy, the Moran directional index evidenced the existence of directional differences in the thematic maps regarding the subregions format, elaborated when considering or not considering the anisotropy in the geostatistical model.In addition, the higher the anisotropy factor, the lower the decay rate of the Moran directional index values, indicating a greater grouping of similar values, that is, a greater spatial continuity of the georeferenced variable, in the direction that defines the geometric anisotropy.It was possible to use the Moran directional index as a measure of directional comparison between the thematic maps, with or without the anisotropy in the geostatistical analysis.
The results obtained allow us to conclude that improving the use of geostatistics applied to sustainable Precision Agriculture can provide important information for better use of the soil, as it allows for an increase in soybean productivity without affecting the environment with unnecessary applications of inputs in agricultural areas.

Figure 1 .
Figure 1.Area of study and locations of sampled points.

Figure 2 .
Figure 2. Boxplot for the values of the Moran directional index (I(d)), calculated according to the distance class (d) and the direction (θ), considering the predicted values of the simulations, for different values of the anisotropic ratio (F a ) ( • indicates discrepant points).

Figure 3 .
Figure 3. Postplot graph spatial representation of the locations sampled in the study area, classified into equal amplitude for the (a) carbon content (C), (b) calcium content (Ca) and (c) potassium content (K).

Figure 4 .
Figure 4. Thematic map of soil chemical property, with the isotropic ("bottom panel") and anisotropic ("top panel") model.

Table 1 .
Percentage of simulations with significant values of the Moran directional index (I(d)) (with 5% significance by the pseudo-significance test) calculated for the predicted values considering the estimated anisotropic model, in each distance class and direction; according to the anisotropic ratio considered in the simulations.

Table 2 .
Descriptive statistics of the soil chemical properties and Pearson's linear correlation coefficient between the soil chemical properties and the x and y coordinates (r(x) and r(y)).
r(x), r(y): Pearson's linear correlation coefficient between the soil chemical properties and the x and y coordinates; *: significant p-value; ns : p-value not significant.

Table 3 .
Estimated values of parameters obtained for anisotropic Gaussian spatial model and similarity measures in the comparison of spatial estimation between anisotropic and isotropic models.

Table 4 .
Estimated values of global comparison measures between thematic maps considering estimated anisotropic and isotropic models for the soil chemical properties.

Table 5 .
Estimated value of Moran's directional spatial autocorrelation index with the predicted values calculated by kriging, considering anisotropic and isotropic models, for the soil chemical properties.