Qgis Impl Georeferencing

qgis-impl-georeferencing

Quick Reference

Core Classes

Transformation Types

Enum Access

from qgis.analysis import QgsGcpTransformerInterface
TM = QgsGcpTransformerInterface.TransformMethod

TM.Linear            # 0
TM.Helmert           # 1
TM.PolynomialOrder1  # 2
TM.PolynomialOrder2  # 3
TM.PolynomialOrder3  # 4
TM.ThinPlateSpline   # 5
TM.Projective        # 6

Critical Warnings

NEVER call warpLayer() on QgsVectorWarper -- this method does NOT exist. ALWAYS use transformFeatures() with a feature iterator and feature sink.

NEVER use QgsImageWarper from Python -- it has NO Python bindings. ALWAYS use GDAL Python bindings (osgeo.gdal) for raster georeferencing.

NEVER use fewer GCPs than the minimum required for the chosen transformation type. The createFromParameters() factory returns None when GCP count is insufficient.

NEVER ignore residuals after fitting a transformation. ALWAYS compute and check residuals to verify accuracy before applying the transformation to data.

NEVER assume all GCPs share the same destination CRS -- each QgsGcpPoint stores its own destinationPointCrs. ALWAYS set the CRS explicitly for every GCP.

ALWAYS check the return value of transformFeatures() and createFromParameters() -- both return a success indicator. On failure, call warper.error() for diagnostics.

Decision Tree: Choosing a Transformation Type

START: How many GCPs do you have?
|
+-- 1-2 GCPs
|   +-- Need rotation? --> Helmert (min 2)
|   +-- No rotation?   --> Linear (min 2)
|   +-- Only 1 GCP?    --> ThinPlateSpline (exact fit, 1 point = translation only)
|
+-- 3-5 GCPs
|   +-- General purpose         --> PolynomialOrder1 (min 3)
|   +-- Perspective correction  --> Projective (min 4)
|
+-- 6-9 GCPs
|   +-- Moderate distortion --> PolynomialOrder2 (min 6)
|   +-- Simple distortion   --> PolynomialOrder1 (extra GCPs improve accuracy)
|
+-- 10+ GCPs
|   +-- Heavy distortion    --> PolynomialOrder3 (min 10)
|   +-- Exact GCP fit needed --> ThinPlateSpline (interpolates through all points)
|   +-- General purpose      --> PolynomialOrder1 (overdetermined = more robust)
|
ACCURACY PRIORITY:
  - Highest global accuracy: PolynomialOrder1 with many well-distributed GCPs
  - Exact fit at GCP locations: ThinPlateSpline (but may oscillate between points)
  - Rigid transformation (no distortion): Helmert

Essential Patterns

Pattern 1: Create Ground Control Points

from qgis.analysis import QgsGcpPoint
from qgis.core import QgsPointXY, QgsCoordinateReferenceSystem

# Each GCP maps a source coordinate to a destination coordinate
gcp = QgsGcpPoint(
    QgsPointXY(100, 200),                           # source (pixel/local coords)
    QgsPointXY(15.5, 47.1),                         # destination (map coords)
    QgsCoordinateReferenceSystem("EPSG:4326"),       # destination CRS
    True                                              # enabled
)

# Access and modify
src = gcp.sourcePoint()          # QgsPointXY
dst = gcp.destinationPoint()     # QgsPointXY
crs = gcp.destinationPointCrs()  # QgsCoordinateReferenceSystem

gcp.setEnabled(False)  # Disable without deleting

Pattern 2: Create a Transformer and Transform Points

from qgis.analysis import QgsGcpTransformerInterface
from qgis.core import QgsPointXY

source_pts = [QgsPointXY(0, 0), QgsPointXY(100, 0), QgsPointXY(100, 100)]
dest_pts = [QgsPointXY(15.0, 47.0), QgsPointXY(16.0, 47.0), QgsPointXY(16.0, 48.0)]

# Create and fit transformer in one step
transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts,
    dest_pts
)

if transformer is None:
    raise RuntimeError("Transformation fit failed -- check GCP count and distribution")

# Transform a single point -- returns (success, x, y) tuple in Python
success, tx, ty = transformer.transform(50.0, 50.0, False)  # False = forward
if success:
    print(f"Transformed: ({tx}, {ty})")

Pattern 3: Vector Georeferencing with QgsVectorWarper

from qgis.analysis import (
    QgsGcpPoint, QgsGcpTransformerInterface, QgsVectorWarper
)
from qgis.core import (
    QgsCoordinateReferenceSystem, QgsFeatureStore,
    QgsPointXY, QgsProject, QgsVectorLayer
)

# Define GCPs
dest_crs = QgsCoordinateReferenceSystem("EPSG:4283")
gcps = [
    QgsGcpPoint(QgsPointXY(90, 210), QgsPointXY(8, 20), dest_crs, True),
    QgsGcpPoint(QgsPointXY(210, 190), QgsPointXY(20.5, 20), dest_crs, True),
    QgsGcpPoint(QgsPointXY(350, 220), QgsPointXY(30, 21), dest_crs, True),
    QgsGcpPoint(QgsPointXY(390, 290), QgsPointXY(39, 28), dest_crs, True),
]

# Create warper
warper = QgsVectorWarper(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    gcps,
    dest_crs
)

# Transform features into a feature store
source_layer = QgsVectorLayer("Point?field=name:string", "source", "memory")
sink = QgsFeatureStore()
success = warper.transformFeatures(
    source_layer.getFeatures(),
    sink,
    QgsProject.instance().transformContext()
)

if not success:
    raise RuntimeError(f"Warp failed: {warper.error()}")

for feature in sink.features():
    print(feature.geometry().asWkt(), feature.attributes())

Pattern 4: Geometry Transformation with QgsGcpGeometryTransformer

from qgis.analysis import QgsGcpGeometryTransformer, QgsGcpTransformerInterface
from qgis.core import QgsGeometry, QgsPointXY

source_pts = [QgsPointXY(0, 0), QgsPointXY(100, 0), QgsPointXY(100, 100)]
dest_pts = [QgsPointXY(15.0, 47.0), QgsPointXY(16.0, 47.0), QgsPointXY(16.0, 48.0)]

# Create geometry transformer directly from coordinates
geo_transformer = QgsGcpGeometryTransformer(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts,
    dest_pts
)

# Transform any QgsGeometry
geom = QgsGeometry.fromPointXY(QgsPointXY(50, 50))
transformed_geom, ok = geo_transformer.transform(geom)

if ok:
    print(f"Transformed: {transformed_geom.asWkt()}")

Pattern 5: Raster Georeferencing with GDAL

from osgeo import gdal, osr

# Step 1: Create GDAL GCPs (target_x, target_y, z, pixel_x, pixel_y)
gcps = [
    gdal.GCP(15.0, 47.0, 0, 0, 0),
    gdal.GCP(16.0, 47.0, 0, 100, 0),
    gdal.GCP(16.0, 48.0, 0, 100, 100),
]

# Step 2: Open source raster and assign GCPs
src_ds = gdal.OpenShared("/path/to/unreferenced.tif", gdal.GA_ReadOnly)
gcp_srs = osr.SpatialReference()
gcp_srs.ImportFromEPSG(4326)
src_ds.SetGCPs(gcps, gcp_srs.ExportToWkt())

# Step 3: Create warped VRT (auto-calculates output dimensions)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(4326)
tmp_ds = gdal.AutoCreateWarpedVRT(
    src_ds, None, dst_srs.ExportToWkt(), gdal.GRA_Bilinear, 0.125
)

# Step 4: Write to GeoTIFF
dst_ds = gdal.GetDriverByName("GTiff").Create(
    "/path/to/georeferenced.tif",
    tmp_ds.RasterXSize,
    tmp_ds.RasterYSize,
    src_ds.RasterCount,
)
dst_ds.SetProjection(dst_srs.ExportToWkt())
dst_ds.SetGeoTransform(tmp_ds.GetGeoTransform())
gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Bilinear)

# Cleanup
dst_ds = None
src_ds = None

Common Operations

Compute Residuals for Accuracy Assessment

from qgis.analysis import QgsGcpTransformerInterface
from qgis.core import QgsPointXY
import math

def compute_residuals(transformer, source_pts, dest_pts):
    """Compute per-GCP residuals and total RMSE."""
    residuals = []
    for src, dst in zip(source_pts, dest_pts):
        success, tx, ty = transformer.transform(src.x(), src.y(), False)
        if not success:
            residuals.append(float("inf"))
            continue
        dx = tx - dst.x()
        dy = ty - dst.y()
        residuals.append(math.sqrt(dx * dx + dy * dy))

    rmse = math.sqrt(sum(r * r for r in residuals) / len(residuals))
    return residuals, rmse

# Usage
source_pts = [QgsPointXY(0, 0), QgsPointXY(100, 0), QgsPointXY(100, 100)]
dest_pts = [QgsPointXY(15.0, 47.0), QgsPointXY(16.0, 47.0), QgsPointXY(16.0, 48.0)]

transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts, dest_pts
)
residuals, rmse = compute_residuals(transformer, source_pts, dest_pts)
print(f"Per-GCP residuals: {residuals}")
print(f"RMSE: {rmse:.6f}")

Write a World File Manually

def write_world_file(path, pixel_width, rotation_x, rotation_y, pixel_height,
                     upper_left_x, upper_left_y):
    """Write a 6-line world file (.tfw, .jgw, .pgw)."""
    with open(path, "w") as f:
        f.write(f"{pixel_width}\n")    # Line 1: pixel width (x scale)
        f.write(f"{rotation_x}\n")     # Line 2: rotation about y axis
        f.write(f"{rotation_y}\n")     # Line 3: rotation about x axis
        f.write(f"{pixel_height}\n")   # Line 4: pixel height (y scale, negative)
        f.write(f"{upper_left_x}\n")   # Line 5: x coordinate of upper-left center
        f.write(f"{upper_left_y}\n")   # Line 6: y coordinate of upper-left center

# Example: 1m resolution, no rotation, origin at (500000, 5200000)
write_world_file(
    "/path/to/image.tfw",
    1.0, 0.0, 0.0, -1.0, 500000.0, 5200000.0
)

World File Extensions by Image Format

Generate World File via GDAL

from osgeo import gdal

# Create GeoTIFF with world file sidecar
ds = gdal.Open("/path/to/georeferenced.tif", gdal.GA_ReadOnly)
gdal.GetDriverByName("GTiff").CreateCopy(
    "/path/to/output.tif", ds, options=["TFW=YES"]
)
ds = None

Disable Individual GCPs for Leave-One-Out Validation

def leave_one_out_validation(gcps, source_pts, dest_pts, method):
    """Disable each GCP in turn and check its residual."""
    from qgis.analysis import QgsGcpTransformerInterface

    for i in range(len(source_pts)):
        # Build lists without the i-th point
        src_subset = [p for j, p in enumerate(source_pts) if j != i]
        dst_subset = [p for j, p in enumerate(dest_pts) if j != i]

        transformer = QgsGcpTransformerInterface.createFromParameters(
            method, src_subset, dst_subset
        )
        if transformer is None:
            print(f"GCP {i}: fit failed without this point")
            continue

        success, tx, ty = transformer.transform(
            source_pts[i].x(), source_pts[i].y(), False
        )
        if success:
            import math
            dx = tx - dest_pts[i].x()
            dy = ty - dest_pts[i].y()
            residual = math.sqrt(dx * dx + dy * dy)
            print(f"GCP {i}: leave-one-out residual = {residual:.6f}")

Reference Links

• references/methods.md -- API signatures for QgsGcpPoint, QgsGcpTransformerInterface, QgsGcpGeometryTransformer, QgsVectorWarper
• references/examples.md -- Complete working examples for vector and raster georeferencing
• references/anti-patterns.md -- What NOT to do when georeferencing

Official Sources

• https://qgis.org/pyqgis/master/analysis/QgsGcpPoint.html
• https://qgis.org/pyqgis/master/analysis/QgsGcpTransformerInterface.html
• https://qgis.org/pyqgis/master/analysis/QgsGcpGeometryTransformer.html
• https://qgis.org/pyqgis/master/analysis/QgsVectorWarper.html
• https://qgis.org/pyqgis/3.40/analysis/index.html