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

qgis-impl-georeferencing — Anti-Patterns

Anti-Pattern 1: Calling warpLayer() on QgsVectorWarper

Wrong

warper = QgsVectorWarper(method, gcps, dest_crs)
warper.warpLayer(source_layer, output_path)  # AttributeError -- method does NOT exist

Why It Fails

QgsVectorWarper has NO warpLayer() method. This method does not exist in the C++ API or the Python bindings. The ONLY method for transforming features is transformFeatures().

Correct

warper = QgsVectorWarper(method, gcps, dest_crs)
sink = QgsFeatureStore()
success = warper.transformFeatures(
    source_layer.getFeatures(),
    sink,
    QgsProject.instance().transformContext(),
)

---

Anti-Pattern 2: Using QgsImageWarper from Python

Wrong

from qgis.analysis import QgsImageWarper  # ImportError -- no Python bindings
warper = QgsImageWarper(method, gcps)
warper.warp(input_path, output_path)

Why It Fails

QgsImageWarper exists in the C++ API but has NO SIP bindings. It is NOT available in the qgis.analysis Python module. Attempting to import it raises ImportError.

Correct

Use GDAL Python bindings for raster georeferencing:

from osgeo import gdal, osr

gcps = [gdal.GCP(map_x, map_y, 0, pixel_x, pixel_y)]
src_ds = gdal.OpenShared(input_path, gdal.GA_ReadOnly)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
src_ds.SetGCPs(gcps, srs.ExportToWkt())
vrt_ds = gdal.AutoCreateWarpedVRT(src_ds, None, srs.ExportToWkt(), gdal.GRA_Bilinear, 0.125)
# ... write output

---

Anti-Pattern 3: Using Too Few GCPs for the Transformation Type

Wrong

# Only 2 GCPs but requesting PolynomialOrder1 (needs minimum 3)
transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    [QgsPointXY(0, 0), QgsPointXY(100, 0)],
    [QgsPointXY(10.0, 50.0), QgsPointXY(11.0, 50.0)],
)
# transformer is None -- fit failed silently
ok, tx, ty = transformer.transform(50, 50)  # AttributeError: NoneType

Why It Fails

Each transformation type has a strict minimum GCP requirement. createFromParameters() returns None when the requirement is not met. Calling methods on None crashes.

Correct

method = QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1
transformer = QgsGcpTransformerInterface.createFromParameters(
    method, source_pts, dest_pts
)
if transformer is None:
    raise RuntimeError(
        f"Fit failed: need at least "
        f"{QgsGcpTransformerInterface.create(method).minimumGcpCount()} GCPs"
    )

---

Anti-Pattern 4: Ignoring the Return Value of transform()

Wrong

# Assuming transform() modifies variables in place (C++ behavior)
x, y = 50.0, 50.0
transformer.transform(x, y, False)
print(f"Transformed: ({x}, {y})")  # Still prints original values!

Why It Fails

In Python, transform() returns a tuple (bool, float, float) due to SIP bindings. The original x and y variables are NOT modified. This is different from the C++ API where output parameters are modified by reference.

Correct

success, tx, ty = transformer.transform(50.0, 50.0, False)
if success:
    print(f"Transformed: ({tx}, {ty})")

---

Anti-Pattern 5: Not Checking Residuals Before Applying Transformation

Wrong

# Create transformer and immediately apply to data without validation
transformer = QgsGcpTransformerInterface.createFromParameters(method, src, dst)
warper = QgsVectorWarper(method, gcps, crs)
warper.transformFeatures(layer.getFeatures(), sink, context)
# No idea if the transformation is accurate

Why It Fails

A successful fit does NOT guarantee accuracy. GCPs may contain errors (wrong coordinate pairs, typos, wrong units). Without residual checks, bad transformations silently corrupt data.

Correct

import math

transformer = QgsGcpTransformerInterface.createFromParameters(method, src_pts, dst_pts)
if transformer is None:
    raise RuntimeError("Fit failed")

# Compute RMSE before applying
sum_sq = 0.0
for s, d in zip(src_pts, dst_pts):
    ok, tx, ty = transformer.transform(s.x(), s.y(), False)
    if ok:
        sum_sq += (tx - d.x()) ** 2 + (ty - d.y()) ** 2

rmse = math.sqrt(sum_sq / len(src_pts))
if rmse > acceptable_threshold:
    raise RuntimeError(f"RMSE {rmse:.6f} exceeds threshold -- check GCPs")

# Only proceed if accuracy is acceptable
warper = QgsVectorWarper(method, gcps, crs)
warper.transformFeatures(layer.getFeatures(), sink, context)

---

Anti-Pattern 6: Mixing GCP Destination CRS Values

Wrong

# Different CRS per GCP without realizing the warper expects consistency
gcps = [
    QgsGcpPoint(src1, dst1, QgsCoordinateReferenceSystem("EPSG:4326"), True),
    QgsGcpPoint(src2, dst2, QgsCoordinateReferenceSystem("EPSG:3857"), True),  # Different CRS!
    QgsGcpPoint(src3, dst3, QgsCoordinateReferenceSystem("EPSG:4326"), True),
]

Why It Fails

Each QgsGcpPoint stores its own destination CRS. When creating a QgsVectorWarper, the class uses transformedDestinationPoint() internally to reproject all destination coordinates to the warper's destination CRS. If CRS values are wrong or inconsistent, the internal reprojection produces incorrect coordinates.

Correct

ALWAYS use a consistent CRS for all destination points, or ensure each GCP's CRS accurately reflects its coordinate values:

dest_crs = QgsCoordinateReferenceSystem("EPSG:4326")
gcps = [
    QgsGcpPoint(src1, dst1, dest_crs, True),
    QgsGcpPoint(src2, dst2, dest_crs, True),
    QgsGcpPoint(src3, dst3, dest_crs, True),
]

---

Anti-Pattern 7: Using Thin Plate Spline for Extrapolation

Wrong

# GCPs only cover a small area, but transforming points far outside
transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.ThinPlateSpline,
    source_pts,  # All clustered in one corner
    dest_pts,
)
# Transform a point far from any GCP
ok, tx, ty = transformer.transform(9999, 9999, False)  # Wildly inaccurate

Why It Fails

Thin Plate Spline interpolates exactly through GCP locations but can produce extreme distortion when extrapolating outside the convex hull of GCPs. The further from GCPs, the more unpredictable the results.

Correct

For areas outside GCP coverage, use polynomial transformations (which extrapolate more smoothly) or add GCPs that bracket the full extent of the data:

# Use PolynomialOrder1 for smoother extrapolation behavior
transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts, dest_pts,
)

---

Anti-Pattern 8: World File with Wrong Pixel Height Sign

Wrong

# Positive pixel height -- image appears flipped vertically
write_world_file("output.tfw", 1.0, 0.0, 0.0, 1.0, 500000.0, 5200000.0)

Why It Fails

In world files, the pixel height (line 4) MUST be negative for north-up images because pixel rows increase downward while map Y coordinates increase upward.

Correct

write_world_file("output.tfw", 1.0, 0.0, 0.0, -1.0, 500000.0, 5200000.0)
#                                                ^^^^  negative for north-up

qgis-impl-georeferencing — Examples

Example 1: Complete Vector Georeferencing Pipeline

Georeference a vector layer with 4 GCPs using polynomial order 1, including residual computation.

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

# 1. Create source layer with test features
source_layer = QgsVectorLayer(
    "Point?field=name:string&field=value:integer", "source", "memory"
)
provider = source_layer.dataProvider()

f1 = QgsFeature()
f1.setAttributes(["station_A", 1])
f1.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(150, 250)))

f2 = QgsFeature()
f2.setAttributes(["station_B", 2])
f2.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(300, 150)))

provider.addFeatures([f1, f2])

# 2. Define GCPs (source local coords -> destination map coords)
dest_crs = QgsCoordinateReferenceSystem("EPSG:4326")
gcps = [
    QgsGcpPoint(QgsPointXY(0, 0), QgsPointXY(10.0, 50.0), dest_crs, True),
    QgsGcpPoint(QgsPointXY(500, 0), QgsPointXY(11.0, 50.0), dest_crs, True),
    QgsGcpPoint(QgsPointXY(500, 500), QgsPointXY(11.0, 51.0), dest_crs, True),
    QgsGcpPoint(QgsPointXY(0, 500), QgsPointXY(10.0, 51.0), dest_crs, True),
]

# 3. Compute residuals BEFORE applying to real data
source_pts = [gcp.sourcePoint() for gcp in gcps]
dest_pts = [gcp.destinationPoint() for gcp in gcps]
method = QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1

transformer = QgsGcpTransformerInterface.createFromParameters(
    method, source_pts, dest_pts
)
if transformer is None:
    raise RuntimeError("Fit failed")

for i, (src, dst) in enumerate(zip(source_pts, dest_pts)):
    ok, tx, ty = transformer.transform(src.x(), src.y(), False)
    if ok:
        residual = math.sqrt((tx - dst.x()) ** 2 + (ty - dst.y()) ** 2)
        print(f"GCP {i}: residual = {residual:.8f}")

# 4. Warp features
warper = QgsVectorWarper(method, gcps, dest_crs)
sink = QgsFeatureStore()
success = warper.transformFeatures(
    source_layer.getFeatures(),
    sink,
    QgsProject.instance().transformContext(),
)

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

# 5. Output results
for feature in sink.features():
    geom = feature.geometry()
    attrs = feature.attributes()
    print(f"{attrs[0]}: {geom.asWkt()}")

---

Example 2: Geometry-Level Transformation

Transform individual geometries without creating a full warper pipeline.

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

# Define corresponding point pairs
source_pts = [
    QgsPointXY(0, 0),
    QgsPointXY(1000, 0),
    QgsPointXY(1000, 1000),
    QgsPointXY(0, 1000),
]
dest_pts = [
    QgsPointXY(500000, 5200000),
    QgsPointXY(501000, 5200000),
    QgsPointXY(501000, 5201000),
    QgsPointXY(500000, 5201000),
]

# Create transformer
geo_tf = QgsGcpGeometryTransformer(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts,
    dest_pts,
)

# Transform a polygon
polygon = QgsGeometry.fromWkt(
    "POLYGON((100 100, 200 100, 200 200, 100 200, 100 100))"
)
result, ok = geo_tf.transform(polygon)
if ok:
    print(f"Transformed polygon: {result.asWkt()}")

# Transform a line
line = QgsGeometry.fromWkt("LINESTRING(50 50, 150 300, 400 200)")
result, ok = geo_tf.transform(line)
if ok:
    print(f"Transformed line: {result.asWkt()}")

---

Example 3: Raster Georeferencing with GDAL

Georeference an unreferenced TIFF using GDAL Python bindings.

from osgeo import gdal, osr

input_path = "/path/to/scan.tif"
output_path = "/path/to/georeferenced.tif"

# 1. Define GCPs: gdal.GCP(map_x, map_y, map_z, pixel_col, pixel_row)
gcps = [
    gdal.GCP(500000.0, 5200000.0, 0, 0, 0),           # top-left
    gdal.GCP(501000.0, 5200000.0, 0, 1000, 0),         # top-right
    gdal.GCP(501000.0, 5199000.0, 0, 1000, 1000),      # bottom-right
    gdal.GCP(500000.0, 5199000.0, 0, 0, 1000),         # bottom-left
]

# 2. Open and assign GCPs
src_ds = gdal.OpenShared(input_path, gdal.GA_ReadOnly)
if src_ds is None:
    raise RuntimeError(f"Cannot open {input_path}")

srs = osr.SpatialReference()
srs.ImportFromEPSG(32633)  # UTM zone 33N
src_ds.SetGCPs(gcps, srs.ExportToWkt())

# 3. Create warped VRT
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(32633)
vrt_ds = gdal.AutoCreateWarpedVRT(
    src_ds, None, dst_srs.ExportToWkt(), gdal.GRA_Bilinear, 0.125
)

# 4. Write output GeoTIFF
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(
    output_path,
    vrt_ds.RasterXSize,
    vrt_ds.RasterYSize,
    src_ds.RasterCount,
    gdal.GDT_Byte,
)
dst_ds.SetProjection(dst_srs.ExportToWkt())
dst_ds.SetGeoTransform(vrt_ds.GetGeoTransform())
gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Bilinear)

# 5. Cleanup
dst_ds = None
vrt_ds = None
src_ds = None
print(f"Georeferenced raster written to {output_path}")

---

Example 4: Point-by-Point Transformation

Transform individual coordinates without geometry objects.

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

source_pts = [
    QgsPointXY(100, 100),
    QgsPointXY(900, 100),
    QgsPointXY(900, 700),
]
dest_pts = [
    QgsPointXY(5.0, 52.0),
    QgsPointXY(6.0, 52.0),
    QgsPointXY(6.0, 53.0),
]

transformer = QgsGcpTransformerInterface.createFromParameters(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    source_pts,
    dest_pts,
)

if transformer is None:
    raise RuntimeError("Fit failed")

# Forward transform: source -> destination
coordinates_to_transform = [(500, 400), (200, 600), (750, 150)]
for x, y in coordinates_to_transform:
    ok, tx, ty = transformer.transform(x, y, False)
    if ok:
        print(f"({x}, {y}) -> ({tx:.6f}, {ty:.6f})")

# Inverse transform: destination -> source
ok, sx, sy = transformer.transform(5.5, 52.5, True)  # True = inverse
if ok:
    print(f"Inverse: (5.5, 52.5) -> ({sx:.2f}, {sy:.2f})")

---

Example 5: Comparing Transformation Methods

Evaluate multiple transformation types on the same GCP set.

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

TM = QgsGcpTransformerInterface.TransformMethod

source_pts = [
    QgsPointXY(50, 50), QgsPointXY(450, 50), QgsPointXY(450, 350),
    QgsPointXY(50, 350), QgsPointXY(250, 200), QgsPointXY(150, 300),
]
dest_pts = [
    QgsPointXY(10.0, 50.0), QgsPointXY(11.0, 50.0), QgsPointXY(11.0, 51.0),
    QgsPointXY(10.0, 51.0), QgsPointXY(10.5, 50.5), QgsPointXY(10.2, 50.8),
]

methods_to_test = [
    TM.Linear,
    TM.Helmert,
    TM.PolynomialOrder1,
    TM.PolynomialOrder2,
    TM.Projective,
    TM.ThinPlateSpline,
]

for method in methods_to_test:
    transformer = QgsGcpTransformerInterface.createFromParameters(
        method, source_pts, dest_pts
    )
    if transformer is None:
        name = QgsGcpTransformerInterface.methodToString(method)
        print(f"{name}: fit FAILED (need {transformer.minimumGcpCount()} GCPs)")
        continue

    # Compute RMSE
    sum_sq = 0.0
    for src, dst in zip(source_pts, dest_pts):
        ok, tx, ty = transformer.transform(src.x(), src.y(), False)
        if ok:
            sum_sq += (tx - dst.x()) ** 2 + (ty - dst.y()) ** 2
    rmse = math.sqrt(sum_sq / len(source_pts))

    name = QgsGcpTransformerInterface.methodToString(method)
    print(f"{name}: RMSE = {rmse:.8f}")

---

Example 6: Write Georeferenced Output to File with Processing

Save warped vector features to a GeoPackage using QgsVectorFileWriter.

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

# Setup source and GCPs (abbreviated)
source_layer = QgsVectorLayer("path/to/input.shp", "input", "ogr")
dest_crs = QgsCoordinateReferenceSystem("EPSG:4326")
gcps = [
    QgsGcpPoint(QgsPointXY(0, 0), QgsPointXY(10.0, 50.0), dest_crs, True),
    QgsGcpPoint(QgsPointXY(100, 0), QgsPointXY(11.0, 50.0), dest_crs, True),
    QgsGcpPoint(QgsPointXY(100, 100), QgsPointXY(11.0, 51.0), dest_crs, True),
]

# Warp
warper = QgsVectorWarper(
    QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1,
    gcps, dest_crs,
)
sink = QgsFeatureStore()
success = warper.transformFeatures(
    source_layer.getFeatures(),
    sink,
    QgsProject.instance().transformContext(),
)

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

# Write to GeoPackage
output_path = "/path/to/output.gpkg"
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "GPKG"

writer = QgsVectorFileWriter.create(
    output_path,
    source_layer.fields(),
    source_layer.wkbType(),
    dest_crs,
    QgsCoordinateTransformContext(),
    options,
)

for feature in sink.features():
    writer.addFeature(feature)

del writer  # Flush and close
print(f"Written to {output_path}")

qgis-impl-georeferencing — Methods Reference

QgsGcpPoint (since QGIS 3.26)

from qgis.analysis import QgsGcpPoint

Constructors

QgsGcpPoint(
    sourcePoint: QgsPointXY,
    destinationPoint: QgsPointXY,
    destinationPointCrs: QgsCoordinateReferenceSystem,
    enabled: bool = True
)

QgsGcpPoint(other: QgsGcpPoint)  # Copy constructor

PointType Enum

QgsGcpPoint.PointType.Source       # 0
QgsGcpPoint.PointType.Destination  # 1

Methods

---

QgsGcpTransformerInterface (since QGIS 3.20)

from qgis.analysis import QgsGcpTransformerInterface

TransformMethod Enum

QgsGcpTransformerInterface.TransformMethod.Linear            # 0
QgsGcpTransformerInterface.TransformMethod.Helmert           # 1
QgsGcpTransformerInterface.TransformMethod.PolynomialOrder1  # 2
QgsGcpTransformerInterface.TransformMethod.PolynomialOrder2  # 3
QgsGcpTransformerInterface.TransformMethod.PolynomialOrder3  # 4
QgsGcpTransformerInterface.TransformMethod.ThinPlateSpline   # 5
QgsGcpTransformerInterface.TransformMethod.Projective        # 6
QgsGcpTransformerInterface.TransformMethod.InvalidTransform  # 65535

Factory Methods

Instance Methods

CRITICAL: In Python, transform() returns a 3-tuple (bool, float, float) due to SIP bindings converting C++ output parameters into return values. The C++ signature shows bool return only, but Python receives (success, x, y).

---

QgsGcpGeometryTransformer (since QGIS 3.18)

from qgis.analysis import QgsGcpGeometryTransformer

Constructors

# From existing transformer
QgsGcpGeometryTransformer(gcpTransformer: QgsGcpTransformerInterface | None)

# Direct initialization with coordinates
QgsGcpGeometryTransformer(
    method: QgsGcpTransformerInterface.TransformMethod,
    sourceCoordinates: Iterable[QgsPointXY],
    destinationCoordinates: Iterable[QgsPointXY]
)

Methods

---

QgsVectorWarper (since QGIS 3.26)

from qgis.analysis import QgsVectorWarper

Constructor

QgsVectorWarper(
    method: QgsGcpTransformerInterface.TransformMethod,
    points: Iterable[QgsGcpPoint],
    destinationCrs: QgsCoordinateReferenceSystem
)

Methods

CRITICAL: There is NO warpLayer() method on QgsVectorWarper. The ONLY method for warping is transformFeatures().

---

GDAL Python Bindings (for Raster Georeferencing)

from osgeo import gdal, osr

Key Functions

Resampling Methods