PyQGIS Engineering & Automation

Troubleshooting CRS Mismatches in PyQGIS Scripts

Diagnose and fix coordinate reference system mismatches in headless PyQGIS scripts: invalid layer CRS, missing transform context, PROJ grid failures, and…

TL;DR: Verify layer.sourceCrs().isValid(), explicitly inject a QgsCoordinateTransformContext from QgsProject.instance(), and always validate the assembled QgsCoordinateTransform against a known control point — never let headless scripts inherit silent GUI defaults.

This page is part of the coordinate transformations and CRS handling guide, which sits inside the broader PyQGIS Core Architecture & Data Handling reference.

Complete Runnable Diagnostic Script

Drop this template into any headless script or CI runner. It covers every common failure mode: invalid layer CRS, empty transform context, broken transformation chains, and missing PROJ grid data. Replace the EPSG codes and control-point coordinates to match your target projection.

python
"""
crs_mismatch_diagnostic.py

Full diagnostic and transform-validation utility for PyQGIS scripts.
Compatible with QGIS 3.28+ / PROJ 8+.

Usage (standalone, no QGIS Desktop):
    python crs_mismatch_diagnostic.py
"""
from __future__ import annotations

import os
import sys
from typing import Optional

from qgis.core import (
    QgsApplication,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsCoordinateTransformContext,
    QgsCoordinateTransformException,
    QgsPointXY,
    QgsProject,
    QgsVectorLayer,
)


def build_crs(epsg: int) -> QgsCoordinateReferenceSystem:
    """
    Construct and validate a CRS from an EPSG code.

    Raises:
        ValueError: if the code produces an invalid CRS object.
    """
    crs = QgsCoordinateReferenceSystem(f"EPSG:{epsg}")
    if not crs.isValid():
        raise ValueError(
            f"EPSG:{epsg} did not produce a valid CRS. "
            "Check the code and ensure PROJ data is available."
        )
    return crs


def diagnose_layer_crs(layer: QgsVectorLayer) -> None:
    """
    Print a structured CRS diagnosis for a single layer.

    Compares sourceCrs() (raw data projection) against the active
    project CRS and reports whether a valid transform can be built.
    """
    src: QgsCoordinateReferenceSystem = layer.sourceCrs()
    proj_crs: QgsCoordinateReferenceSystem = QgsProject.instance().crs()
    ctx: QgsCoordinateTransformContext = QgsProject.instance().transformContext()

    print(f"[Layer]   {layer.name()}")
    print(f"  sourceCrs : {src.authid()} | valid={src.isValid()}")
    print(f"  project   : {proj_crs.authid()} | valid={proj_crs.isValid()}")

    if not src.isValid():
        print("  DIAGNOSIS: Layer has no valid CRS — assign one explicitly.")
        return

    if not proj_crs.isValid():
        print("  DIAGNOSIS: Project CRS is unset — call QgsProject.instance().setCrs().")
        return

    transform = QgsCoordinateTransform(src, proj_crs, ctx)
    print(f"  transform : isValid={transform.isValid()}")
    if not transform.isValid():
        print("  DIAGNOSIS: Transformation chain invalid — check PROJ grids and EPSG support.")


def transform_with_validation(
    x: float,
    y: float,
    src_epsg: int,
    dst_epsg: int,
    *,
    control: Optional[tuple[float, float]] = None,
    tolerance_m: float = 0.001,
) -> QgsPointXY:
    """
    Transform a point and optionally assert it matches a known control value.

    Args:
        x:            Source X (easting or longitude).
        y:            Source Y (northing or latitude).
        src_epsg:     Source CRS EPSG code.
        dst_epsg:     Destination CRS EPSG code.
        control:      Expected (X, Y) for the result, used for accuracy assertion.
        tolerance_m:  Max allowable deviation from the control point (metres).

    Returns:
        QgsPointXY in the destination CRS.

    Raises:
        ValueError: CRS construction failed.
        RuntimeError: Transform invalid or coordinate value out of bounds.
    """
    src_crs = build_crs(src_epsg)
    dst_crs = build_crs(dst_epsg)

    # Always pull the populated context from the active project (even in
    # headless mode — QgsProject.instance() returns a valid singleton after
    # QgsApplication.initQgis()).
    ctx = QgsProject.instance().transformContext()
    transform = QgsCoordinateTransform(src_crs, dst_crs, ctx)

    if not transform.isValid():
        raise RuntimeError(
            f"QgsCoordinateTransform(EPSG:{src_epsg} → EPSG:{dst_epsg}) is invalid. "
            "Verify PROJ_LIB contains required grid shift files."
        )

    try:
        result: QgsPointXY = transform.transform(QgsPointXY(x, y))
    except QgsCoordinateTransformException as exc:
        raise RuntimeError(f"Transform raised an exception: {exc}") from exc

    if result.x() in (0.0, float("inf")) or result.y() in (0.0, float("inf")):
        raise RuntimeError(
            "Transform returned a degenerate result (0 or inf). "
            "This usually means swapped axis order or an unsupported datum."
        )

    if control is not None:
        dx = abs(result.x() - control[0])
        dy = abs(result.y() - control[1])
        if dx > tolerance_m or dy > tolerance_m:
            raise AssertionError(
                f"Control-point assertion failed: got ({result.x():.4f}, {result.y():.4f}), "
                f"expected ({control[0]:.4f}, {control[1]:.4f}), "
                f"delta=({dx:.4f}, {dy:.4f}) > tolerance={tolerance_m}."
            )

    return result

Diagnostic Flow

The decision tree below captures the complete triage path from “something is wrong with my coordinates” to the correct corrective action. Follow each branch in order — fixing an upstream problem (invalid layer CRS) usually resolves symptoms that look like downstream failures (broken transform chains).

CRS Mismatch Diagnostic Flow Flowchart showing how to triage a CRS mismatch: start with layer CRS validity, then transform context, then transform validity, then PROJ grids, and finally control-point assertion. CRS mismatch suspected layer.sourceCrs() .isValid()? no Fix .prj or assign CRS yes transformContext populated? no Inject project transform ctx yes QgsCoordinateTransform .isValid()? no Check PROJ grids/EPSG yes control-point within tolerance? no Wrong grid or axis order issue yes Transform validated

Architecture Breakdown

QgsCoordinateReferenceSystem — contract and gotchas

QgsCoordinateReferenceSystem wraps PROJ’s internal CRS representation. Construction accepts EPSG strings ("EPSG:4326"), WKT2 text, or PROJ strings, but does not raise an exception on failure — you must call isValid() after construction. An object built from a malformed WKT string will silently produce an invalid CRS and every downstream transform will also silently fail.

In headless scripts, PROJ resolves CRS definitions against its data directory. If PROJ_LIB is not set or points to an incomplete installation, even standard EPSG codes may return invalid objects. Always verify the environment variable at script startup.

python
import os
from qgis.core import QgsCoordinateReferenceSystem

# Validate PROJ data directory before constructing any CRS
proj_lib: str = os.environ.get("PROJ_LIB", "")
if not proj_lib or not os.path.isdir(proj_lib):
    raise EnvironmentError(
        f"PROJ_LIB is unset or invalid: '{proj_lib}'. "
        "QGIS requires a complete PROJ data directory for CRS resolution."
    )

crs = QgsCoordinateReferenceSystem("EPSG:25832")
assert crs.isValid(), f"CRS construction failed for EPSG:25832"
print(crs.toWkt())  # inspect the resolved WKT2 definition

sourceCrs() vs crs() on a layer is a common point of confusion. crs() may reflect on-the-fly reprojection applied during the desktop session, while sourceCrs() always returns the CRS embedded in the data file. In automated scripts, always use sourceCrs() to keep transformations deterministic and independent of any inherited project state.

QgsCoordinateTransformContext — why it matters headlessly

The transform context stores user-defined datum shift preferences — for instance, preferring the NTv2 grid shift over the approximate 3-parameter Helmert when both are available for a given CRS pair. In QGIS Desktop, this context is populated from the project file. Headless scripts that call QgsCoordinateTransformContext() directly get an empty context, causing PROJ to fall back to whatever transformation it deems most appropriate, which may be a low-accuracy approximation.

The correct pattern is to pull the context from QgsProject.instance(). In standalone scripts where no project file is loaded, build a minimal context and document which transformation pipeline you are relying on:

python
from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransformContext,
    QgsProject,
)


def get_transform_context() -> QgsCoordinateTransformContext:
    """
    Return the best available transform context for headless use.

    Prefers the project context (populated when a .qgz file is loaded).
    Falls back to a default context for standalone scripts — callers
    should document which transformation accuracy this implies.
    """
    project_ctx: QgsCoordinateTransformContext = (
        QgsProject.instance().transformContext()
    )
    # A non-empty context implies the project carries explicit datum
    # shift preferences; reuse it directly.
    return project_ctx

QgsCoordinateTransform — validity and exception behaviour

A QgsCoordinateTransform object can be constructed without error even when the underlying PROJ pipeline is broken. Call isValid() before any transform attempt. Additionally, transform() raises QgsCoordinateTransformException for points outside the CRS domain of validity — for example, transforming a latitude-longitude pair in the Southern Ocean into a CRS defined only for northern Europe. Always wrap transform calls in a try/except block in production code.

The transformBoundingBox() method applies the same validity rules but also handles axis-order differences transparently, making it safer than manually transforming corner points when dealing with geographic CRS (where X/Y ordering conventions vary between PROJ versions).

Wiring the Diagnostic into a Standalone Script

For a typical headless batch job — reading shapefiles, reprojecting features, and writing output — wire the diagnostic utility as a pre-flight check before entering the main processing loop:

python
"""
batch_reproject.py

Standalone PyQGIS script demonstrating pre-flight CRS validation
before batch reprojection. No QGIS Desktop required.
"""
from __future__ import annotations

import os
import sys

from qgis.core import (
    QgsApplication,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsCoordinateTransformException,
    QgsFeature,
    QgsProject,
    QgsVectorFileWriter,
    QgsVectorLayer,
)

# 1. Bootstrap QGIS application (required for all standalone scripts)
QgsApplication.setPrefixPath("/usr", True)
qgs = QgsApplication([], False)
qgs.initQgis()

SRC_PATH = "/data/input/survey_points.shp"
DST_EPSG = 32632  # UTM zone 32N
OUT_PATH = "/data/output/survey_utm32.gpkg"

# 2. Load layer and run pre-flight CRS check
layer = QgsVectorLayer(SRC_PATH, "survey", "ogr")
if not layer.isValid():
    sys.exit(f"Failed to load layer: {SRC_PATH}")

src_crs: QgsCoordinateReferenceSystem = layer.sourceCrs()
if not src_crs.isValid():
    sys.exit(
        f"Layer '{layer.name()}' has no valid CRS. "
        "Add a .prj file or set the CRS explicitly before running this script."
    )

dst_crs = QgsCoordinateReferenceSystem(f"EPSG:{DST_EPSG}")
if not dst_crs.isValid():
    sys.exit(f"Destination EPSG:{DST_EPSG} is not a recognised CRS.")

ctx = QgsProject.instance().transformContext()
transform = QgsCoordinateTransform(src_crs, dst_crs, ctx)
if not transform.isValid():
    sys.exit(
        "Transformation pipeline is invalid. "
        "Verify PROJ_LIB and the presence of required grid shift files."
    )

# 3. Control-point pre-flight assertion
# Replace with a known survey benchmark in your area of interest.
CONTROL_SRC_X, CONTROL_SRC_Y = 13.4050, 52.5200   # Berlin, WGS 84
EXPECTED_DST_X, EXPECTED_DST_Y = 392969.0, 5820045.0
TOLERANCE_M = 1.0  # 1-metre tolerance for this quick sanity check

try:
    from qgis.core import QgsPointXY
    result = transform.transform(QgsPointXY(CONTROL_SRC_X, CONTROL_SRC_Y))
    dx = abs(result.x() - EXPECTED_DST_X)
    dy = abs(result.y() - EXPECTED_DST_Y)
    if dx > TOLERANCE_M or dy > TOLERANCE_M:
        sys.exit(
            f"Control-point assertion failed: "
            f"got ({result.x():.1f}, {result.y():.1f}), "
            f"expected ({EXPECTED_DST_X}, {EXPECTED_DST_Y}), "
            f"delta=({dx:.1f}, {dy:.1f}) m."
        )
    print(f"Pre-flight OK: control point within {max(dx, dy):.3f} m of expected.")
except QgsCoordinateTransformException as exc:
    sys.exit(f"Control-point transform raised an exception: {exc}")

# 4. Main processing loop (features are now safe to transform)
# ... feature iteration, geometry transformation, and output writing go here

qgs.exitQgis()

See optimizing feature iteration with QgsVectorLayer.getFeatures() for efficient iteration patterns to combine with the CRS validation above.

Production Best Practices

  • Always use sourceCrs(), not crs(). The latter reflects desktop on-the-fly reprojection state and is non-deterministic in headless contexts.
  • Validate isValid() on every QgsCoordinateReferenceSystem object immediately after construction — before building any transform.
  • Verify PROJ_LIB at startup. Missing grid files cause silent accuracy degradation, not hard failures. Log the resolved path with QgsMessageLog so it appears in QGIS log panels.
  • Wrap transform.transform() in try/except QgsCoordinateTransformException. Points outside a CRS domain of validity raise this exception rather than returning a sentinel value.
  • Assert against a control point in every batch job. A ±1 m tolerance is suitable for most projected CRS; tighten to ±0.001 m for surveying or cadastral workflows.
  • Document the PROJ version and grid files your script depends on. Grid availability differs between QGIS packaging channels (OSGeo4W, conda-forge, system packages). Pin these in your project’s requirements.txt or Docker image.
  • Test with QgsCoordinateTransformContext() deliberately empty to confirm your script does not silently fall back to low-accuracy transforms when the project context is missing.
  • Check transform.sourceDatumTransformId() and transform.destinationDatumTransformId() to confirm PROJ selected the datum shift you intended, not a fallback.