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:

Figure — from a covariance matrix to a drawable ellipse, in four exact steps.

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:

Figure — scale factors for two-dimensional confidence levels; the one-dimensional numbers do not apply.

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"

Reading an Ellipse Field #

A single ellipse describes one point; a field of them describes a network, and it is read differently. Uniform, near-circular ellipses across the whole area mean the geometry is even and no direction is privileged — the ordinary result of a well-distributed network. Ellipses that elongate consistently along one bearing mean the observations constrain that direction poorly, usually because the control lies along a line rather than around a perimeter. Ellipses that grow steadily with distance from a single monument mean the network is effectively cantilevered off that mark, and its error propagates into everything downstream of it.

The practical value of plotting the field rather than tabulating the semi-axes is that these three patterns are obvious at a glance and nearly invisible in a table of numbers. Plot the ellipses at the same exaggeration factor, state that factor on the figure, and keep the axes equally scaled — an ellipse field on unequal axes reproduces the elongation pattern of the plot rather than of the network.

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.

Frequently Asked Questions #

Why does my ellipse look like a circle when the covariance is clearly not isotropic? #

Nine times out of ten the axes are not equally scaled. Matplotlib will happily stretch a plot to fill its box, which turns an elongated ellipse into a circle and hides exactly the property you drew it for. Set the aspect ratio to equal on any plot showing a spatial error region — an error ellipse plotted on unequal axes is not just imprecise, it is misleading about the direction of the weakness.

Which way is the orientation angle measured? #

The eigenvector gives an angle measured counter-clockwise from the positive x axis, which in an east-north plot is measured from east. Surveyors normally quote the ellipse orientation as an azimuth, measured clockwise from north. The two differ by a reflection and a 90 degree offset, so an ellipse that is correct on the plot can be quoted incorrectly in the accompanying table. Convert explicitly and state which convention the number is in.

Should the plot show one-sigma or 95 per cent ellipses? #

Whichever the specification asks for — but label it on the figure itself, not only in the caption. Figures get extracted from reports and pasted into other documents, and an unlabelled ellipse is read as whatever the reader expects. Where both are useful, draw them as nested contours from the same covariance so the relationship is visible.

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