Detecting & Modeling Geometric Anisotropy in Python

TL;DR: Build several DirectionalVariogram(coords, vals, azimuth=a, tolerance=22.5, bandwidth=b) objects across azimuths, read each .parameters[0] range, and find the longest (major axis) and its perpendicular (minor axis). The anisotropy ratio is major_range / minor_range. Correct it by rotating coordinates to the major axis and rescaling the minor axis by the ratio, then fit one isotropic Variogram.

Why This Matters #

Most real spatial processes are directional: contamination follows groundwater flow, ore follows structural trends, temperature follows prevailing wind. Assuming an isotropic variogram in those settings smears the surface across the grain of the process and biases every prediction. Detecting and correcting direction is the whole purpose of studying anisotropy and directional variograms, a specialization within the broader variogram modeling and semivariance analysis workflow. scikit-gstat makes the detection tractable through its DirectionalVariogram, and a simple coordinate transformation turns the finding into a model kriging can use.

Under geometric anisotropy the sill is constant but the range varies with direction, tracing an ellipse. If the major-axis range is $a_{max}$ at azimuth $\theta$ and the minor-axis range is $a_{min}$, the anisotropy ratio and angle are:

$$R = \frac{a_{max}}{a_{min}} \ge 1, \qquad \theta = \arg\max_{\phi} \; a(\phi)$$

Rescaling the minor axis by $R$ maps the ellipse of ranges back to a circle, restoring isotropy in the transformed coordinates.

Environment #

pip install "scikit-gstat>=1.0" "numpy>=1.23" "matplotlib>=3.7"

import numpy as np
import matplotlib.pyplot as plt
from skgstat import Variogram, DirectionalVariogram

Step-by-Step Implementation #

Step 1 — Build a field with directional structure #

rng = np.random.default_rng(29)
coords = rng.uniform(0, 1000, size=(260, 2))

# Correlation reaches farther along the east-west (x) axis than north-south (y):
# stretching x by 1/3 before the smooth kernel lengthens the x-range.
xs = coords[:, 0] / 3.0
ys = coords[:, 1]
vals = 10 + 4 * np.sin(np.sqrt(xs ** 2 + ys ** 2) / 40.0) + rng.normal(0, 0.7, 260)

This synthetic field has a longer correlation range roughly east-west, so we expect the major axis near an azimuth of 0 degrees.

Step 2 — Compute directional variograms across azimuths #

azimuths = [0, 45, 90, 135]     # degrees, counter-clockwise from the +x axis (east)
dir_ranges = {}

for az in azimuths:
    DV = DirectionalVariogram(
        coords, vals,
        azimuth=az,
        tolerance=22.5,     # half-width of the directional search cone
        bandwidth=80.0,     # max perpendicular offset of a pair from the search line
        n_lags=12,
        maxlag="median",
        model="spherical",
        use_nugget=True,
    )
    dir_ranges[az] = DV.parameters[0]
    print(f"azimuth {az:3d} deg  ->  range {DV.parameters[0]:7.1f}")

Each DirectionalVariogram only counts point pairs whose orientation falls within tolerance of the azimuth and within bandwidth of the search line. The fitted range at each azimuth samples the ellipse of ranges.

Step 3 — Identify the major and minor axes #

major_az = max(dir_ranges, key=dir_ranges.get)
minor_az = (major_az + 90) % 180        # perpendicular direction

# Ensure we have the perpendicular measured; recompute if it was not scanned
if minor_az not in dir_ranges:
    DVp = DirectionalVariogram(coords, vals, azimuth=minor_az,
                               tolerance=22.5, bandwidth=80.0,
                               n_lags=12, maxlag="median",
                               model="spherical", use_nugget=True)
    dir_ranges[minor_az] = DVp.parameters[0]

major_range = dir_ranges[major_az]
minor_range = dir_ranges[minor_az]
print(f"major axis {major_az} deg  range {major_range:.1f}")
print(f"minor axis {minor_az} deg  range {minor_range:.1f}")

The azimuth with the longest range is the major axis; the direction 90 degrees from it is the minor axis.

Step 4 — Quantify the anisotropy ratio and angle #

anisotropy_ratio = major_range / minor_range      # >= 1
anisotropy_angle = major_az
print(f"anisotropy angle: {anisotropy_angle} deg")
print(f"anisotropy ratio: {anisotropy_ratio:.2f}")

A ratio near 1 means the structure is effectively isotropic; a ratio of 2 or more is strong geometric anisotropy worth modelling explicitly.

Step 5 — Correct anisotropy by transforming coordinates #

theta = np.deg2rad(anisotropy_angle)
rot = np.array([[ np.cos(theta), np.sin(theta)],
                [-np.sin(theta), np.cos(theta)]])

rotated = coords @ rot.T          # align major axis with the new x-axis
rotated[:, 1] *= anisotropy_ratio # stretch the minor axis so ranges match

V_iso = Variogram(rotated, vals, model="spherical", n_lags=12,
                  maxlag="median", use_nugget=True)
print("isotropic-equivalent range:", round(V_iso.parameters[0], 1))
print("isotropic-equivalent sill :", round(V_iso.parameters[1], 2))

After rotating to the major axis and rescaling the minor axis by the ratio, the ellipse of ranges becomes a circle, so a single isotropic Variogram now describes the structure in every direction. This transformed model is what you feed to an anisotropic kriging routine.

Interpreting the Output #

The set of directional ranges traces the anisotropy ellipse. Two numbers summarize it: the anisotropy angle, the azimuth of longest correlation, and the anisotropy ratio, how many times farther correlation reaches along the major axis than the minor. A constant sill across directions with only the range changing confirms geometric anisotropy; if the sill itself changes with direction, the process shows zonal anisotropy, which the coordinate transformation alone does not fix.

Confirming a constant sill also depends on the variable being stationary in the first place, which ties back to stationarity and trend analysis.

Critical Best Practices #

Scan enough azimuths to resolve the ellipse #

Four azimuths (0, 45, 90, 135 degrees) are the minimum to detect anisotropy; use finer steps such as every 22.5 degrees when the angle matters for downstream kriging. Too few directions can miss a major axis lying between your sampled azimuths.

Keep tolerance and bandwidth balanced #

A narrow tolerance sharpens directional resolution but starves each direction of pairs. Widen tolerance or bandwidth until every directional variogram clears the roughly 30-pairs-per-bin threshold, a trade-off explored in choosing lag bins and bandwidth for variograms.

Distinguish geometric from zonal anisotropy #

The coordinate transformation only corrects geometric anisotropy, where the sill is constant and only the range changes. If the sill itself varies with direction, you have zonal anisotropy and need a nested model with a directional component rather than a simple rescaling.

Verify the ratio is real, not sampling noise #

Directional ranges wobble when pairs are sparse. Before acting on an anisotropy ratio, confirm the pattern is stable across nearby azimuths and that each direction had enough pairs, so you are modelling the process and not the sampling geometry.

Troubleshooting #

Next Steps #

With the anisotropy angle and ratio in hand, revisit how the underlying models are fitted in fitting spherical, exponential and Gaussian variogram models, and confirm your directional binning is well populated using choosing lag bins and bandwidth for variograms.

Related

• Anisotropy & Directional Variograms — the directional modelling framework
• Choosing Lag Bins and Bandwidth for Variograms — tuning tolerance and bandwidth for directional cones
• Fitting Spherical, Exponential & Gaussian Variogram Models — the models fitted along each direction

← Back to Anisotropy & Directional Variograms

Frequently Asked Questions #

What is geometric anisotropy in a variogram?

Geometric anisotropy means the range of spatial correlation depends on direction while the sill stays the same. Correlation reaches farther along one axis, the major axis, than along the perpendicular minor axis, tracing an ellipse of ranges rather than a circle. It is common in wind-driven, fluvial and structurally controlled processes.

How does scikit-gstat define the azimuth for a DirectionalVariogram?

The azimuth is given in degrees and measured counter-clockwise from the positive x-axis, so 0 degrees points east and 90 degrees points north. The tolerance sets the angular half-width of the search cone and the bandwidth caps how far a pair may sit off the search line.

How do I correct geometric anisotropy before kriging?

Rotate the coordinates so the major axis aligns with an axis, then rescale the minor axis by the anisotropy ratio so ranges become equal in every direction. Fitting a single isotropic variogram on the transformed coordinates then models the anisotropic structure correctly.