Visualizing Geodetic Error Ellipses with Matplotlib: Deterministic Confidence Regions

Rendering a 95 % positional confidence region from a station’s 2×2 covariance submatrix must reproduce the same semi-axes and orientation byte-for-byte across CI runs and stay within ±1e-9 m of the analytic ellipse, or the plotted uncertainty no longer matches the adjustment it claims to represent. This operation is the visualization endpoint of error distribution modeling in Python, the parent topic under Algorithmic Math & Geodetic Workflows: the covariance produced by a least-squares pass is decomposed here into a drawable Matplotlib Ellipse patch whose geometry a reviewer can trust for cadastral sign-off.

Concept: From Covariance to Confidence Region #

The positional uncertainty of a planimetric station is captured by the 2×2 submatrix extracted from the inverted normal equations of an adjustment. Its eigenvectors give the principal axes of the uncertainty, and its eigenvalues give the variance along each. To turn variances into an axis length at a chosen confidence level we scale by the chi-square quantile with two degrees of freedom — the correct distribution for a bivariate region, not the 1.96σ factor that applies only to a one-dimensional interval:

where $\lambda_1 \ge \lambda_2$ are the eigenvalues, $v_1$ is the eigenvector of $\lambda_1$, and $\chi^2_{2}(0.95) \approx 5.991$. The semi-major axis $a$, semi-minor axis $b$, and rotation $\theta$ fully determine the patch.

Two numerical hazards make this harder than a textbook plot. First, when affine and polynomial stages are chained, covariance matrices accumulate IEEE-754 drift and a near-singular matrix can present a tiny negative eigenvalue, sending sqrt into a domain error or — worse — silently producing an inverted ellipse. Using the symmetric/Hermitian solver numpy.linalg.eigh rather than the general numpy.linalg.eig guarantees real eigenvalues and orthonormal eigenvectors, and a clamp at machine-epsilon scale absorbs the residual drift. The same numerical-stability discipline governs least-squares adjustment for control networks, where a single non-finite entry voids the normal-equation system. Second, the eigenvector sign returned by LAPACK is platform-dependent, so $\theta$ is reduced modulo 180° to make the rendered orientation deterministic.

Complete Runnable Implementation #

The function below is self-contained and type-hinted. It validates symmetry, decomposes with eigh, sorts eigenvalues descending so a ≥ b, clamps near-zero negatives, scales by the chi-square quantile, and normalizes the rotation. It raises a typed error rather than emitting a corrupted patch.

from __future__ import annotations

from dataclasses import dataclass

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import chi2


class GeodeticPrecisionError(ValueError):
    """Raised when a covariance matrix violates survey-grade tolerance bounds."""


@dataclass(frozen=True)
class ErrorEllipse:
    a: float          # semi-major axis (metres)
    b: float          # semi-minor axis (metres)
    theta_deg: float  # rotation of major axis, degrees in [0, 180)
    chi2_scale: float


def error_ellipse_params(
    cov: np.ndarray,
    confidence: float = 0.95,
    tol: float = 1e-10,
) -> ErrorEllipse:
    """Decompose a 2x2 covariance matrix into confidence-ellipse geometry.

    Implements the bivariate chi-square confidence region (df=2) used for
    positional uncertainty in least-squares geodetic adjustments.

    Args:
        cov: 2x2 symmetric, positive-semidefinite covariance matrix (m^2).
        confidence: Target confidence level, e.g. 0.95 for the 95% region.
        tol: Absolute tolerance for the symmetry and definiteness checks.

    Returns:
        ErrorEllipse with semi-axes in metres and rotation in [0, 180).

    Raises:
        GeodeticPrecisionError: matrix is non-symmetric or not PSD within tol.
    """
    m = np.asarray(cov, dtype=np.float64)
    if m.shape != (2, 2):
        raise ValueError("Covariance matrix must be 2x2.")
    if not np.allclose(m, m.T, atol=tol):
        raise GeodeticPrecisionError("Non-symmetric covariance matrix.")

    # Symmetric solver: real eigenvalues, orthonormal eigenvectors.
    eigvals, eigvecs = np.linalg.eigh(m)
    order = np.argsort(eigvals)[::-1]          # descending -> a >= b
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]

    if eigvals[1] < -tol:
        raise GeodeticPrecisionError(
            f"Negative eigenvalue {eigvals[1]:.3e}: matrix not PSD."
        )
    eigvals = np.maximum(eigvals, 0.0)         # absorb sub-epsilon drift

    scale = float(chi2.ppf(confidence, df=2))  # chi^2(0.95, 2) ~= 5.991
    a = float(np.sqrt(eigvals[0] * scale))
    b = float(np.sqrt(eigvals[1] * scale))

    # Reduce mod 180 so LAPACK eigenvector-sign flips do not change output.
    theta = float(np.degrees(np.arctan2(eigvecs[1, 0], eigvecs[0, 0])) % 180.0)
    return ErrorEllipse(a=a, b=b, theta_deg=theta, chi2_scale=scale)


def plot_error_ellipse(
    ax: plt.Axes,
    x: float,
    y: float,
    ell: ErrorEllipse,
    edgecolor: str = "#1f77b4",
    facecolor: str = "none",
    linewidth: float = 1.5,
) -> Ellipse:
    """Add a confidence ellipse patch centred at (x, y) to a Matplotlib axis."""
    patch = Ellipse(
        xy=(x, y),
        width=2.0 * ell.a,    # Ellipse takes full axis lengths, not semi-axes
        height=2.0 * ell.b,
        angle=ell.theta_deg,
        edgecolor=edgecolor,
        facecolor=facecolor,
        linewidth=linewidth,
    )
    ax.add_patch(patch)
    return patch

Parameter and Return Reference #

Minimal Worked Example #

A station whose easting/northing covariance is [[4.0, 1.5], [1.5, 2.0]] m² yields a 95 % region with a major axis tilted toward the correlated direction:

cov = np.array([[4.0, 1.5], [1.5, 2.0]])
ell = error_ellipse_params(cov, confidence=0.95)
print(round(ell.a, 4), round(ell.b, 4), round(ell.theta_deg, 3))
# -> 5.3643 2.6783 28.155

fig, ax = plt.subplots()
plot_error_ellipse(ax, x=0.0, y=0.0, ell=ell)
ax.set_aspect("equal")
ax.set_xlim(-7, 7)
ax.set_ylim(-7, 7)

The semi-major axis of 5.3643 m and orientation of 28.155° are reproducible regardless of the LAPACK eigenvector sign convention, because the rotation is reduced modulo 180°.

Validation Check #

Gate the geometry before it reaches a renderer or an audit export. The area of the 95 % region equals π·a·b = χ²·√det(cov), an identity independent of the eigen-decomposition that catches any silent corruption of the axes:

import math

ell = error_ellipse_params(np.array([[4.0, 1.5], [1.5, 2.0]]))
expected_area = ell.chi2_scale * math.sqrt(np.linalg.det([[4.0, 1.5], [1.5, 2.0]]))
assert math.isclose(math.pi * ell.a * ell.b, expected_area, rel_tol=1e-9), \
    "ellipse area drifted from chi2 * sqrt(det(cov)) — geometry is corrupt"
assert ell.a >= ell.b >= 0.0, "axes not sorted / negative semi-axis"

Common Mistakes #

• Scaling by 1.96 instead of √χ²₂. The 1.96σ factor is the one-dimensional 95 % interval; a two-dimensional region needs √5.991 ≈ 2.448. Using 1.96 understates the ellipse by about 20 % and produces non-conservative confidence regions that fail cadastral review.
• Passing semi-axes to Ellipse. matplotlib.patches.Ellipse consumes full width/height, so the semi-axes must be doubled. Forgetting the factor of two draws a region half the correct size — a subtle error because the shape still looks plausible.
• Using eig and ignoring the sign / negative-eigenvalue drift. numpy.linalg.eig can return complex values and unsorted, sign-arbitrary eigenvectors. After chained affine transformations for local grids and polynomial shift algorithms for regional adjustments, a covariance can carry a sub-epsilon negative eigenvalue; without the eigh solver plus the PSD clamp this becomes a sqrt domain error or an inverted patch.

Related #

• Error Distribution Modeling in Python — the parent topic that produces the covariance matrices rendered here
• Python Script for Geodetic Inverse Problem Solving — generates the distance/azimuth observations whose residuals feed the covariance
• Least-squares Adjustment for Control Networks — the adjustment whose inverted normal equations supply the 2×2 submatrix
• Implementing Affine Transformations for Local Grids — a linear stage whose drift can perturb the covariance before plotting
• Algorithmic Math & Geodetic Workflows — the full deterministic pipeline reference