The KRIGE2D Procedure

Geometric Anisotropy

Geometric anisotropy is the simplest type of anisotropy. It occurs when the same sill (or scale) parameter $c_0$ is present in all directions but the range $a_0$ changes with direction. In geometric anisotropy the covariance model uses the same forms in all directions.

Therefore, geometric anisotropy features one single sill value, and depending on the direction the semivariogram reaches the sill within a different distance. This is illustrated in Figure 55.12, where an anisotropic exponential semivariogram is plotted. Assume that the two curves displayed in this figure have the same sill $c_0 = 1.5$ and are generated using the ranges $a_{0,1}=3$ in the direction $\theta _{1}=30^{\circ }$ (effective range is $r_{\epsilon ,1}=9$) and $a_{0,2}=1$ in the direction $\theta _{2}=120^{\circ }$ (effective range is $r_{\epsilon ,2}=3$).

As you can see from the figure, the ratio of the shorter to longer range is $R = 1/3$. The anisotropy factor R is the value to use in the RATIO= parameter in the MODEL statement in PROC KRIGE2D. When you model geometric anisotropy $R \le 1$. In fact, isotropy is a partial case of geometric anisotropy for which ${a_{0}}^{\mathit{min}} = {a_{0}}^{\mathit{max}}$ and $R = 1$.

Figure 55.12: Geometric Anisotropy with Major Axis in the Direction $\theta _{1}=30^{\circ }$

Geometric Anisotropy with Major Axis in the Direction θ1=30○


The values of the RANGE= and ANGLE= parameters in the MODEL statement in PROC KRIGE2D are set based on the major anisotropy axis characteristics. Specifically, the RANGE= parameter is the value of the major axis range ${a_{0}}^{\mathit{max}} = a_{0,1}$, and the ANGLE= parameter is the angle $\theta _{1}$ of the major axis measured clockwise from north (angles measured in this way are also known as azimuths). You can then specify the following MODEL statement in PROC KRIGE2D to approximate the covariance structure: 

model form=exp range=3 scale=1.5 angle=30 ratio=0.3333;

If you use a nested model, provide the type for each one of the nested structures with the FORM= option, and assign the individual SCALE= parameters so that they add up to the total sill (include in the sum the nugget effect, if present). In the typical case, all of your nested structures have the same anisotropy axes. This means that you specify the same ANGLE= parameter value for all structures. Each structure likely has its own values for the RANGE= and RATIO= parameters depending on the degree of its contribution to the nested model.

The terminology associated with geometric anisotropy is that of ellipses. To see how this comes about, consider the following hypothetical set of calculations. Let {$Z(\bm {s}), \bm {s} \in D \subset \mathcal{R}^2$} be a geometrically anisotropic process, and assume sufficient data points are present to calculate an experimental semivariogram at a large number of angle classes $\theta \in \{ 0,\delta \theta , 2\delta \theta , \cdots ,180^{\circ }$}. At each of these angles $\theta $, the experimental semivariogram is plotted and the range $a_{0}$ is recorded. A diagram in polar coordinates $(a_{0},\theta )$ yields an ellipse with the major axis ${a_{0}}^{\mathit{max}}$ in the direction of the largest $a_{0}$ and the minor axis ${a_{0}}^{\mathit{min}}$ perpendicular to it. For the example in Figure 55.12, the ellipse is shown in Figure 55.13(a). Its major axis has size ${a_{0}}^{\mathit{max}}$ situated at angle $\theta _{1}$ clockwise from north, and the minor axis has size ${a_{0}}^{\mathit{min}}$ oriented at angle $\theta _{2}$ clockwise from north.

The KRIGE2D procedure handles geometric anisotropy by applying a reversible transformation in two steps that converts geometric anisotropy into isotropic conditions.

The first step is to align your coordinates axes with the anisotropy ellipse axes. Specifically, you choose to rotate by an angle $\varphi $ the standard Cartesian orientation of the $(x,y)$ coordinates system shown in Figure 55.13(a) so that the Y axis coincides with the ellipse minor axis. The rotation result is illustrated in Figure 55.13(b). The second step is to elongate the minor axis so its length equals that of the major axis of the ellipse. You can see the result in Figure 55.13(c). The computational details are shown in the following.

Figure 55.13: Transformation Applied to Geometric Anisotropy

Transformation Applied to Geometric Anisotropy


The transformation angle $\varphi $ is measured in standard Cartesian orientation counterclockwise from the X axis (east). If the major axis azimuth is $\theta _{1}$, then the Cartesian system of $(x,y)$ needs to be rotated by $\varphi =90^{\circ }-\theta _{1}$ so that the Y axis can coincide with the ellipse minor axis; see Figure 55.13(a).

Let us call the ellipse major axis X’ and the minor axis Y’. The transformation that converts any coordinates in the $(x,y)$ system into $(x’,y’)$ coordinates in terms of $\varphi $ is given by the matrix:

\[ \mb{H}= \left( \begin{array}{cc} \cos (\varphi ) & \sin (\varphi ) \\ -\sin (\varphi ) & \cos (\varphi ) \end{array} \right) \]

The elongation of the minor axis in the second step is performed with the matrix:

\[ \mb{D}_{R}= \left( \begin{array}{cc} 1 & 0 \\ 0 & 1/R \end{array} \right) \]

Note: These two steps are sequential and their order cannot be reversed. For any point pair $P_1$ and $P_2$ with respective coordinates $\bm {s}_{1}=(x_{1},y_{1})$ and $\bm {s}_{2}=(x_{2},y_{2})$ in the $(x,y)$ axes, their distance is given by

\[ \mid P_ iP_ j \mid _{(x,y)}\ = h = \sqrt {(\delta x)^2+(\delta y)^2} \]

where the distance components $\delta x = x_{2}-x_{1}$ and $\delta y=y_{2}-y_{1}$. Based on the previous, the corresponding distances $\delta x’$ and $\delta y’$ in the $(x’,y’)$ coordinates system are given by the vector:

\[ \left( \begin{array}{c} \delta x’ \\ \delta y’ \end{array} \right) = \mb{D}_{R} \ \mb{H} \left( \begin{array}{c} \delta x \\ \delta y \end{array} \right) = \left( \begin{array}{cc} \cos (\varphi ) & \sin (\varphi ) \\ -\sin (\varphi )/R & \cos (\varphi )/R \end{array} \right) \left( \begin{array}{c} \delta x \\ \delta y \end{array} \right) \]

The transformed interpair distance is then:

\[ \mid P_ iP_ j \mid _{(x',y')}\ = h’ = \sqrt {(\delta x’)^2+(\delta y’)^2} \]

As a result, the original anisotropic semivariogram in Figure 55.12 that was a function $\gamma (\bm {h}) = \gamma (h,\theta )$ of both h and $\theta $ is then transformed to an equivalent function $\hat{\gamma }(h’)$ only of $h’$:

\[ \hat{\gamma }(h’) = \gamma (\bm {h}) \]

This single isotropic semivariogram $\hat{\gamma }(h’)$ is then used for kriging purposes.

The two steps used by PROC KRIGE2D in the previous analysis can also be performed in a different manner. For instance, you might equivalently choose to rotate the $(x,y)$ Cartesian coordinates so that the Y axis coincides with the ellipse major axis, rather than with the minor axis as was shown earlier. Also, you might prefer to compress the major axis rather than elongating the short one. In any case, you need to perform the appropriate computations for the transformation of your choice.