Anisotropy & Directional Variograms

Spatial continuity is rarely the same in every direction. Contaminant plumes stretch downstream, ore bodies follow structural trends, and rainfall correlates further along prevailing winds than across them. When correlation depends on direction as well as distance, the field is anisotropic, and an isotropic variogram will systematically misrepresent it. This page, part of Variogram Modeling & Semivariance Analysis, shows how to detect directional dependence with directional variograms, quantify it with an anisotropy ratio and angle, and reduce it to isotropy through a rotation and rescaling transform that feeds directly into kriging.

Prerequisites #

• Python 3.10+
• numpy>=1.24, scipy>=1.10, scikit-gstat>=1.0, geopandas>=0.14, matplotlib>=3.7
• A fitted isotropic baseline from theoretical variogram models for comparison
• Sample coordinates in a projected (metric) CRS so azimuths and ranges are geometrically meaningful
• Enough samples per direction: roughly 30+ pairs per lag per azimuth, which usually means 100+ observations

Mathematical Core #

An isotropic variogram assumes semivariance depends only on the scalar separation $h = \|\mathbf{h}\|$. Under anisotropy it depends on the full lag vector $\mathbf{h}$, including its direction $\phi$:

$$\gamma(\mathbf{h}) = \gamma(h, \phi), \qquad h = \|\mathbf{h}\|, \quad \phi = \operatorname{atan2}(h_y, h_x).$$

A directional variogram estimates $\gamma(h, \phi)$ by restricting each lag class to pairs whose separation vector points within an angular tolerance $\Delta\phi$ of a target azimuth $\phi$, optionally capped by a lateral bandwidth $b$ so distant pairs do not leak across directions.

Geometric anisotropy #

Under geometric anisotropy the sill is constant but the range varies with direction, tracing an anisotropy ellipse. The longest range $a_{\max}$ lies along the major axis at azimuth $\theta$ (the anisotropy angle); the shortest range $a_{\min}$ lies along the perpendicular minor axis. The anisotropy ratio is

$$\lambda = \frac{a_{\min}}{a_{\max}} \in (0, 1],$$

where $\lambda = 1$ is isotropy and small $\lambda$ is strong directional continuity. Geometric anisotropy is removable: rotate coordinates so the major axis aligns with a reference axis, then rescale the minor axis by $1/\lambda$. In the transformed space the field is isotropic and a single ordinary model can be fitted. The transform is

$$\begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1/\lambda \end{bmatrix} \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix},$$

a rotation by $-\theta$ followed by a stretch of the minor axis. Equivalently, the effective isotropic distance between two points is the reduced distance

$$h_{\text{eff}} = \sqrt{\left(\frac{h_{\parallel}}{a_{\max}}\right)^2 + \left(\frac{h_{\perp}}{a_{\min}}\right)^2},$$

where $h_{\parallel}$ and $h_{\perp}$ are the lag components along the major and minor axes.

Read physically, $h_{\text{eff}}$ is distance measured in units of range rather than in metres. Two points 400 m apart along the major axis, with $a_{\max} = 800$ m, sit at $h_{\text{eff}} = 0.5$ — exactly as correlated as two points only 160 m apart across it when $a_{\min} = 320$ m, since $\lambda = 0.4$. The contour $h_{\text{eff}} = 1$ is the anisotropy ellipse: the locus of separations at which correlation has died out. Every intermediate contour is a scaled copy of it, so the correlation structure is a family of concentric, similarly oriented ellipses. That is why the transform works — rescaling the minor axis by $1/\lambda$ turns those concentric ellipses into concentric circles, and a circle is precisely what an isotropic model assumes. It also states plainly what the transform cannot do: because one $\theta$ and one $\lambda$ apply to the entire domain, the ellipse is assumed identical everywhere in it.

Zonal anisotropy #

Under zonal anisotropy the sill itself differs by direction, usually because a directional trend injects extra variance along one axis so the variogram never reaches the same plateau. No coordinate transform can remove it. It is handled by nesting structures: an isotropic component plus a directional structure with a very large range in the low-variance direction, so it contributes variance only along the high-variance axis.

The anisotropy ellipse #

Every sample “sees” further along the major axis than across it. Kriging that ignores this ellipse treats all directions as the minor axis (over-smoothing) or the major axis (over-reaching); supplying $\theta$ and $\lambda$ makes the search neighbourhood and weights honour the true geometry.

Annotated Implementation #

scikit-gstat exposes directional variograms through the DirectionalVariogram class, which adds an azimuth, tolerance, and bandwidth to the standard estimator. The workflow is: estimate along several azimuths, fit each, read the ranges, and derive the ratio and angle.

import numpy as np
import geopandas as gpd
from skgstat import DirectionalVariogram

# ── 1. Load projected samples (metres) ──────────────────────────────────────
gdf = gpd.read_file("plume_samples.gpkg").to_crs("EPSG:32633")
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
values = gdf["concentration"].to_numpy(dtype=float)

# ── 2. Estimate directional variograms at 0, 45, 90, 135 degrees ────────────
azimuths = [0, 45, 90, 135]     # degrees, measured from the x (east) axis
tolerance = 22.5                # +/- half-sector so the four cover the circle
directional = {}

for az in azimuths:
    dv = DirectionalVariogram(
        coordinates=coords,
        values=values,
        azimuth=az,
        tolerance=tolerance,
        bandwidth="q33",     # lateral cap so far pairs don't leak across sectors
        model="spherical",
        n_lags=12,
        maxlag="median",
        fit_method="trf",
        fit_sigma="linear",
    )
    rng, sill, nugget = dv.parameters
    directional[az] = {"vario": dv, "range": rng, "sill": sill, "nugget": nugget}
    print(f"azimuth {az:3d} deg -> range = {rng:8.1f} m,  sill = {sill:.3f}")

The bandwidth argument is what stops a nominally 90-degree lag from including a pair that is far to the side but technically within the angular wedge at long distance; capping it keeps directional variograms crisp. With sparse data you can widen tolerance toward 45 degrees to gather enough pairs, at the cost of directional resolution.

Identifying axes and deriving ratio and angle #

def anisotropy_from_directional(directional):
    """Derive anisotropy angle (major-axis azimuth) and ratio a_min / a_max
    from a dict of {azimuth_deg: {'range': ...}} directional fits."""
    az_range = {az: d["range"] for az, d in directional.items()}

    major_az = max(az_range, key=az_range.get)   # longest range = major axis
    a_max = az_range[major_az]

    # minor axis is the azimuth closest to 90 deg from the major axis
    minor_az = min(az_range, key=lambda az: abs(((az - major_az) % 180) - 90))
    a_min = az_range[minor_az]

    ratio = a_min / a_max
    return {"angle_deg": float(major_az), "ratio": float(ratio),
            "a_max": float(a_max), "a_min": float(a_min),
            "major_az": major_az, "minor_az": minor_az}

aniso = anisotropy_from_directional(directional)
print(f"anisotropy angle  = {aniso['angle_deg']:.1f} deg (major axis)")
print(f"anisotropy ratio  = {aniso['ratio']:.3f}")
print(f"major range a_max = {aniso['a_max']:.1f} m")
print(f"minor range a_min = {aniso['a_min']:.1f} m")

A ratio near 1.0 across all four azimuths means the field is effectively isotropic and you can fall back to a single isotropic model. A ratio well below 1.0 with a stable sill across directions is the signature of geometric anisotropy that the transform below can remove.

The rotation and rescaling transform #

def reduce_to_isotropy(coords, angle_deg, ratio):
    """Rotate to the anisotropy angle and rescale the minor axis so the field
    becomes isotropic. Returns transformed coordinates for isotropic fitting."""
    theta = np.deg2rad(angle_deg)
    c, s = np.cos(theta), np.sin(theta)

    rotation = np.array([[c, s], [-s, c]])        # rotate by -theta
    scaling = np.array([[1.0, 0.0], [0.0, 1.0 / ratio]])  # stretch minor axis
    transform = scaling @ rotation

    return coords @ transform.T                    # shape (n, 2)

coords_iso = reduce_to_isotropy(coords, aniso["angle_deg"], aniso["ratio"])

# Refit a single isotropic model in the transformed space
from skgstat import Variogram
V_iso = Variogram(coords_iso, values, model="spherical", n_lags=12,
                  maxlag="median", fit_method="trf", fit_sigma="linear")
print("isotropic-equivalent range:", V_iso.parameters[0])

In the transformed coordinates the directional differences collapse, so an isotropic variogram fits cleanly. That single model, together with the recorded angle and ratio, is exactly what an anisotropic kriging call needs.

Passing the anisotropic model to kriging #

# pykrige expects the major-axis range plus an anisotropy scaling and angle.
# anisotropy_scaling is a_max / a_min = 1 / ratio; angle is measured clockwise.
from pykrige.ok import OrdinaryKriging

OK = OrdinaryKriging(
    coords[:, 0], coords[:, 1], values,
    variogram_model="spherical",
    variogram_parameters={
        "sill": aniso["a_max"] and directional[aniso["major_az"]]["sill"],
        "range": aniso["a_max"],
        "nugget": directional[aniso["major_az"]]["nugget"],
    },
    anisotropy_scaling=1.0 / aniso["ratio"],
    anisotropy_angle=aniso["angle_deg"],
)

The kriging engine now stretches its neighbourhood along the major axis so aligned samples receive more weight, which is the whole point of characterising anisotropy in the first place.

Output Interpretation #

Read the four directional fits as a set, not individually:

• Same sill, different ranges is geometric anisotropy. The ranges trace the ellipse; take the longest as the major axis and derive the ratio and angle. This is the correctable case.
• Different sills by direction signals zonal anisotropy, usually an unmodelled directional trend. A coordinate transform will not fix it; detrend first (see theoretical variogram models for handling trend-inflated variograms) or add a nested zonal structure.
• Ratio interpretation. $\lambda > 0.8$ is mild anisotropy often not worth modelling; $0.3$ to $0.8$ is moderate and worth incorporating; $< 0.3$ is strong directional continuity that materially changes kriging results.
• Angle stability. Re-estimate with a slightly different tolerance; if the major-axis azimuth jumps around, the anisotropy is weak or the directional variograms are pair-starved. A robust anisotropy angle is stable across reasonable tolerance choices.

Warning signs include directional variograms with fewer than ~30 pairs in the near lags (unreliable ranges), a major axis that aligns suspiciously with the sampling grid orientation (a sampling artefact, not real continuity), and a ratio that changes drastically when you widen the bandwidth.

Production Considerations #

Pair counts drive everything. Splitting the data into four directions divides the available pairs roughly fourfold per direction, and the bandwidth cap trims further. Inspect dv.bin_count per azimuth; if near lags are thin, widen the tolerance, drop to two directions (major and minor guesses), or accept an isotropic model. Directional estimation is still $O(n^2)$ in pairs, so subsample large datasets before scanning many azimuths.

Scan before you commit. Rather than four fixed azimuths, sweep azimuths in 15-degree steps and plot range versus azimuth (a rose diagram of ranges). The maximum locates the major axis more precisely than a coarse four-direction grid and reveals whether the ellipse is well defined or noisy.

When the anisotropy is not stationary. Everything above assumes one ellipse for the whole domain, and a directional variogram averages over every pair in the study area to estimate it. A field whose strike rotates gradually — a river valley that bends, a fold whose orientation changes across the property — therefore returns a major-axis range shorter than either local range and a ratio closer to 1 than the truth. The anisotropy is diluted rather than detected, and the failure is silent: the fit still converges and the numbers still look plausible. Two symptoms give it away. The range-versus-azimuth rose comes out broad and flat instead of peaked, and splitting the samples into spatial halves yields major axes that disagree by more than the angular tolerance. The remedies are locally varying anisotropy (kriging with a per-node angle field derived from a smoothed trend gradient or from mapped structural data) or simply partitioning the domain into zones and modelling each separately. Both cost sample density, so partitioning only pays when each zone still clears the 30-pairs-per-lag bar on its own.

Consistency with kriging conventions. Angle conventions differ: scikit-gstat measures azimuth from the x-axis counter-clockwise, while some kriging libraries measure clockwise from north. Convert carefully and validate by kriging a known transect; a sign error in the angle produces a surface elongated in exactly the wrong direction.

Reproducibility. Store the anisotropy angle, ratio, major range, sill, and nugget together with the CRS and the tolerance/bandwidth used, so the anisotropic model can be reconstructed exactly for the kriging pipeline.

Troubleshooting #

Next Steps #

For a focused, reproducible walk-through on real data, see detecting and modeling geometric anisotropy in Python. Because directional ranges are simply the isotropic ranges resolved by direction, revisit theoretical variogram models to fit each directional curve, and empirical variogram estimation for the binning that underlies every directional estimate. The anisotropic model you build here is consumed by the kriging, interpolation and surface generation techniques workflow, where the ratio and angle stretch the prediction neighbourhood along the true axis of continuity.

Related

• Detecting & Modeling Geometric Anisotropy in Python — end-to-end worked example on real samples
• Theoretical Variogram Models — fitting spherical, exponential and Gaussian curves to each directional variogram
• Empirical Variogram Estimation — the lag binning that directional estimation extends
• Kriging, Interpolation & Surface Generation Techniques — where the anisotropy ratio and angle shape prediction weights

← Back to Variogram Modeling & Semivariance Analysis