Batch Reproject Vector Layers in PyQGIS

Data collected over years arrives in a patchwork of coordinate systems: a national grid from the mapping agency, UTM from a GPS survey, WGS 84 from a web export, a legacy local system from an old CAD drawing. QGIS reprojects on the fly for display, but analysis, exchange and storage are far simpler when everything shares one CRS. Reprojecting dozens or hundreds of layers by hand is slow and error-prone; a script does it consistently — with the same datum transformation every time, which matters more than most people expect.

This recipe belongs to Coordinate Reference Systems. It inventories layers and their CRSs, reprojects them with an explicit transformation, writes the results into one GeoPackage while keeping names and styles, skips what does not need work, verifies positions, and logs every step. Rasters are covered separately in batch reprojecting raster datasets.

Inventory, transform, write, verifyEach source layer is inventoried for its CRS. Layers already in the target CRS are copied unchanged. Others are reprojected with native:reprojectlayer using one explicitly chosen datum transformation, written as tables in a single target GeoPackage with their original names and styles. A check point with known coordinates in both systems verifies the result, and a log records every layer.One target CRS, one transformation, one logsourcesshp, gpkg,geojsoninventoryCRS per layerskip if targetreprojectexplicittransformationGeoPackagenames, styles+ log

Prerequisites

Inventory the sources

Before transforming anything, list what is there and in which CRS. The inventory decides what needs work and exposes layers with missing or suspicious CRS definitions.

from pathlib import Path
from collections import Counter
from qgis.core import QgsVectorLayer, QgsProviderRegistry

def vector_sources(folder):
    for path in sorted(Path(folder).rglob("*")):
        if path.suffix.lower() not in {".shp", ".gpkg", ".geojson", ".fgb", ".gml"}:
            continue
        md = QgsProviderRegistry.instance().providerMetadata("ogr")
        for sub in md.querySublayers(str(path)) if path.suffix.lower() == ".gpkg" else [None]:
            uri = sub.uri() if sub else str(path)
            name = sub.name() if sub else path.stem
            layer = QgsVectorLayer(uri, name, "ogr")
            if layer.isValid() and layer.isSpatial():
                yield name, uri, layer

inventory = [(name, uri, lyr.crs()) for name, uri, lyr in vector_sources("/data/inbox")]
print(Counter(crs.authid() or "(no authid)" for _, _, crs in inventory).most_common())
for name, _, crs in inventory:
    if not crs.isValid():
        print("NO CRS:", name)

Breakdown: querySublayers lists the tables inside each GeoPackage so every layer is found, not just the first; plain files are treated as one layer each. Counting CRSs by authority id shows the spread at a glance — usually a handful of systems. A layer with no valid CRS cannot be reprojected safely; resolve it before continuing, as described in handling missing CRS. A layer with a custom CRS and no authority id is worth a look too: it is often a known system with a slightly different definition.

Choose the transformation explicitly

Between two CRSs on different datums there can be several transformations with very different accuracy — from a few metres to a few centimetres. Letting each machine pick its own default makes results depend on which grids happen to be installed.

Several ways from A to BFrom a national grid on an old datum to ETRS89 UTM there are typically several transformations: a three-parameter Helmert shift accurate to a few metres, a seven-parameter Helmert accurate to about a metre, and a grid-based transformation accurate to centimetres if the grid file is installed. A batch should pick one, record it, and use it for every layer.Pick one, record it, use it everywhere3-param Helmert2–5 malways available7-param Helmert≈ 1 musually availablegrid-basedcm levelneeds grid file

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

target = QgsCoordinateReferenceSystem("EPSG:25832")
source = QgsCoordinateReferenceSystem("EPSG:31467")           # DHDN / Gauss-Krüger zone 3

ops = QgsDatumTransform.operations(source, target)
for i, op in enumerate(ops):
    print(i, f"±{op.accuracy} m" if op.accuracy >= 0 else "accuracy unknown",
          "available" if op.isAvailable else "MISSING GRID", op.name)

best = next(op for op in ops if op.isAvailable)
context = QgsCoordinateTransformContext()
context.addCoordinateOperation(source, target, best.proj)
print("using:", best.name)

Breakdown: QgsDatumTransform.operations lists every PROJ operation between the two systems with its accuracy and whether its grid files are present. Choosing the first available operation — PROJ orders them by accuracy — and pinning it in a transform context makes every layer use exactly that pipeline. If the best operation is marked as missing its grid, install the grid first, as covered in installing PROJ datum grids, rather than silently falling back to a metre-level shift.

Reproject and write to one GeoPackage

With the transformation fixed, each layer is reprojected and appended as a table to a single target GeoPackage, keeping its name.

import processing
from qgis.core import QgsProcessingContext, QgsProcessingFeedback

out_gpkg = "/data/harmonised/all_layers_25832.gpkg"
pctx = QgsProcessingContext()
pctx.setTransformContext(context)
first, log = True, []

for name, uri, crs in inventory:
    if crs == target:
        layer = QgsVectorLayer(uri, name, "ogr")
        action = "copied"
    else:
        layer = processing.run("native:reprojectlayer", {
            "INPUT": uri, "TARGET_CRS": target, "OPERATION": best.proj,
            "OUTPUT": "memory:"}, context=pctx)["OUTPUT"]
        action = f"reprojected from {crs.authid()}"
    processing.run("native:savefeatures", {
        "INPUT": layer, "OUTPUT": out_gpkg, "LAYER_NAME": name,
        "ACTION_ON_EXISTING_FILE": 0 if first else 1})
    first = False
    log.append((name, action, layer.featureCount()))
    print(f"{name:<30} {action:<32} {layer.featureCount():>7}")

Breakdown: Passing the PROJ pipeline string as OPERATION makes native:reprojectlayer use exactly the chosen transformation, independent of project settings; the transform context on the Processing context provides the same guarantee for any other step. Layers already in the target CRS are copied, not reprojected, so their coordinates are untouched. Writing every layer as a table in one GeoPackage gives a single, self-contained deliverable. Duplicate names across sources would overwrite each other — make names unique, for example by prefixing the source folder, before writing.

Keep styles with the layers

Reprojection creates new layers without the styling the originals had. GeoPackage can store styles alongside tables, so they can be copied across.

for name, uri, crs in inventory:
    src = QgsVectorLayer(uri, name, "ogr")
    qml = Path(uri.split("|")[0]).with_suffix(".qml")
    dst = QgsVectorLayer(f"{out_gpkg}|layername={name}", name, "ogr")
    if qml.exists():
        dst.loadNamedStyle(str(qml))
    else:
        dst.setRenderer(src.renderer().clone())
    dst.saveStyleToDatabase(name, "default style", True, "")

Breakdown: A sidecar QML next to the source file takes priority; otherwise the source layer's current renderer is cloned. saveStyleToDatabase writes the style into the GeoPackage's layer_styles table and marks it as default, so the layer opens styled anywhere. More on database-stored styles in saving styles to a database.

Verify positions with a check point

A wrong transformation shifts data by metres, which is invisible at most zoom levels. A point with known coordinates in both systems — a survey benchmark, a trig point — verifies the result numerically.

A check point in both systemsA benchmark with published coordinates in the source system is transformed with the chosen operation and compared with its published coordinates in the target system. A difference within the expected accuracy, here a few centimetres, confirms the transformation; a difference of metres means a weaker transformation or a missing grid was used.Known point, known answerbenchmarksource coordspublishedtransformchosen operationcomparetarget coordsΔ ≤ accuracy?

from qgis.core import QgsCoordinateTransform, QgsPointXY

xf = QgsCoordinateTransform(source, target, context)
benchmark_src = QgsPointXY(3_478_912.410, 5_421_337.880)       # published DHDN GK3
benchmark_dst = QgsPointXY(478_853.973, 5_419_603.112)         # published ETRS89 UTM32
got = xf.transform(benchmark_src)
delta = got.distance(benchmark_dst)
print(f"check point difference: {delta:.3f} m")
assert delta < 0.10, "transformation less accurate than expected"

Breakdown: Transforming a benchmark with the same context the batch used and comparing with its published target coordinates tests the whole chain. A difference of a few centimetres confirms a grid-based transformation; a metre or more means a Helmert fallback was used. Use a benchmark inside the data's area, since grid-based accuracy varies across a country. The coordinates above are illustrative; take real ones from your mapping agency.

Write the log

A batch that changed coordinates of a hundred layers should leave a record: what was done to each layer, with which transformation, and how many features went through.

import csv
from datetime import datetime

with open(out_gpkg.replace(".gpkg", "_log.csv"), "w", newline="", encoding="utf-8") as fh:
    w = csv.writer(fh)
    w.writerow(["layer", "action", "features", "target", "operation", "run"])
    for name, action, count in log:
        w.writerow([name, action, count, target.authid(), best.name,
                    datetime.now().isoformat(timespec="seconds")])

Breakdown: The log answers the questions that come up months later — was this layer transformed, with what — without anyone having to reconstruct the run. Including the operation name records the datum transformation, which is the detail most often lost and most often responsible for small, puzzling offsets between datasets.

QGIS version compatibility

native:reprojectlayer with the OPERATION parameter, QgsDatumTransform.operations and transform contexts are available on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. querySublayers on provider metadata exists from 3.22. The list of operations depends on the PROJ version and installed grids.

Troubleshooting

  • Layers end up hundreds of metres off. A source CRS was wrong; fix the label before reprojecting.
  • Differences of a metre or two between machines. Different default transformations; pin the operation.
  • Layers overwrite each other in the GeoPackage. Duplicate layer names; make them unique.
  • Styles are lost. Copy QML or renderers and save them to the GeoPackage.

Conclusion

Inventory every layer's CRS, fix missing or wrong labels first, list and pin one datum transformation, reproject with that operation into one GeoPackage while skipping layers already in the target, carry styles across, verify with a benchmark, and log every layer with the transformation used.

Frequently Asked Questions

Do I need to reproject if QGIS displays everything correctly? For display, no. For analysis, exchange and storage, a single CRS avoids repeated on-the-fly transformations and their subtle differences.

Can I reproject in place? Not safely. Write new outputs and keep the originals until verified.

What about layers in a PostGIS database? Use ST_Transform in SQL, or reproject with Processing and import, as in importing a layer into PostGIS.

How do I handle sources in several different CRSs? List operations for each source–target pair, pin one per pair in the transform context, and pass the matching pipeline when reprojecting each layer. The log should then record the operation per layer rather than one for the whole run.

Does reprojection change attribute values? No. Only geometry changes — but stored area or length fields become stale and should be recomputed.