Choosing a Transformation API: PROJ, pyproj, or Manual Grids

Deciding how you invoke a datum transformation is a design decision that outlives any single batch, and it belongs to the wider discipline of Batch Transformation & Automation for Cadastral Coordinate Pipelines: the same corner of ground can be moved from one realization to another by a raw PROJ pipeline string, by the high-level pyproj.Transformer API, or by reading the grid and interpolating the shift yourself — and for cadastral work the three routes differ far more in auditability, control, and reproducibility than in the final centimetre they land on. A parcel boundary that a tribunal may inspect years later is not served by “whatever operation PROJ happened to pick”; it is served by knowing, and being able to reproduce, the exact steps that were applied. This page frames that choice precisely. It routes a concrete use case to one of the three approaches, transforms one point three ways so the numbers can be compared side by side, and sets out the trade-offs on the axes that matter to a survey deliverable rather than to a web map.

The three approaches are not competitors so much as a ladder of abstraction. pyproj.Transformer.from_crs sits at the top: you name a source and target coordinate reference system and let PROJ choose the operation. A +proj=pipeline string sits in the middle: you name every step, every grid, and the order they execute, and PROJ merely runs what you wrote. Manual grid interpolation sits at the bottom: you open the NTv2 grid, find the four surrounding nodes, and compute the bilinear shift yourself, inspecting each node as you go. Higher on the ladder means less code and more convenience; lower means more control and more that you must get right. The whole of this page is about knowing which rung a given job needs.

Figure — routing a use case to manual grids, an explicit PROJ pipeline, or the high-level pyproj Transformer.

The Three Approaches, Precisely Defined #

Before comparing them it is worth pinning down exactly what each one is, because the differences are easy to blur.

Raw PROJ pipeline string. PROJ exposes a small language for chaining coordinate operations: +proj=pipeline +step +proj=... +step +proj=hgridshift +grids=.... Every conversion, every grid file, and the order in which they run is written out explicitly. pyproj can execute such a string through Transformer.from_pipeline, so you keep Python ergonomics while surrendering none of the specificity. The string itself becomes a text artifact you can commit, diff, and cite in a submission.

High-level pyproj.Transformer API. Transformer.from_crs(src, dst) asks PROJ’s operation factory to search for the best available coordinate operation between two CRSs, ranking candidates by accuracy and grid availability. It is concise and fast to write, and for the overwhelming majority of GIS work it is exactly right. Its one liability for cadastral work is that “best available” is resolved at run time against the installed PROJ data, so the chosen operation can change when the environment changes. The relationship between the explicit and the searched form is examined in depth in PROJ pipeline strings vs the pyproj Transformer API.

Manual grid interpolation. Here you bypass PROJ’s transformation engine for the grid step and read the NTv2 .gsb yourself, locate the four nodes around the query point, and evaluate the bilinear shift in your own code. This is the most work and the narrowest tool, but it is the only one that lets you inspect each contributing node, weight or reject nodes by their accuracy columns, and attach a per-point uncertainty. The mechanics of reading that grid are covered in understanding NTv2 grid shift files in Python, and the specific cases where doing it by hand is justified are laid out in when to use manual grid interpolation over pyproj.

The bilinear kernel that both PROJ and a manual reader apply at each grid cell is identical in form. For a query at fractional position $(t_x, t_y)$ inside a cell whose four corner shifts are $s_{00}, s_{10}, s_{01}, s_{11}$, the interpolated shift is

Because the arithmetic is the same, a correctly written manual reader agrees with PROJ to well within a millimetre; the value of writing it yourself is never a different answer but a fully inspectable one.

Comparison Matrix #

The following table scores the three approaches on the axes that decide a cadastral deliverable. “Control” means how precisely you dictate the operation; “audit” means how easily the exact computation can be reproduced and defended.

The pattern is clear. from_crs wins on ease; the pipeline string wins on the audit-and-reproducibility axis while keeping the engine’s speed; manual interpolation wins only where you genuinely need to open the grid, and it pays for that with complexity and throughput. For most certified cadastral pipelines the pinned pipeline string is the sweet spot — full control at native speed — with the high-level API reserved for exploratory work and manual interpolation reserved for the narrow QA and custom-grid cases.

A Decision Procedure #

Work the axes in order of how expensive they are to get wrong. The following procedure resolves the choice in four checks, and each check maps onto a branch of the decision tree above.

Step 1 — Ask whether you must see inside the grid #

If the deliverable requires inspecting the individual node shifts, weighting or rejecting nodes by their published accuracy, or embedding a per-point uncertainty derived from the accuracy columns, no engine-level API will expose that — choose manual interpolation and stop. This is a small fraction of jobs, but a real one, and it is the first question because it is the only one whose answer forces the most laborious route.

from dataclasses import dataclass


@dataclass(frozen=True)
class ApiChoice:
    approach: str          # "manual" | "pipeline" | "transformer"
    reason: str


def needs_node_inspection(
    inspect_nodes: bool,
    embed_per_node_uncertainty: bool,
    custom_unregistered_grid: bool,
) -> ApiChoice | None:
    """Return a manual-interpolation choice when the job requires opening the
    grid, else None so the caller falls through to the next question."""
    if inspect_nodes or embed_per_node_uncertainty or custom_unregistered_grid:
        return ApiChoice(
            approach="manual",
            reason="per-node inspection, embedded uncertainty, or a custom grid",
        )
    return None

Step 2 — Ask whether the operation must be pinned for audit #

If the result will be certified, submitted to an agency, or may need to be reproduced identically years later on a different machine, the operation must not be chosen by a run-time search — write an explicit pipeline string so the exact steps and grids are fixed in text.

def needs_pinned_operation(
    certified_deliverable: bool,
    reproduce_across_versions: bool,
) -> ApiChoice | None:
    """Return a pipeline-string choice when the exact operation must be fixed
    and reproducible, else None to fall through to the throughput question."""
    if certified_deliverable or reproduce_across_versions:
        return ApiChoice(
            approach="pipeline",
            reason="operation must be pinned and reproducible for audit",
        )
    return None

Step 3 — Default to the high-level API #

If neither of the first two checks fired, the job is interactive, exploratory, or a bulk transform where the standard best operation is acceptable. The concise Transformer.from_crs is the right tool; it is vectorized and cached, so throughput is a non-issue on the same C engine the pipeline uses.

def resolve_api(
    inspect_nodes: bool = False,
    embed_per_node_uncertainty: bool = False,
    custom_unregistered_grid: bool = False,
    certified_deliverable: bool = False,
    reproduce_across_versions: bool = False,
) -> ApiChoice:
    """Route a use case to exactly one of the three approaches by asking the
    questions in ascending order of how costly they are to answer wrongly."""
    manual = needs_node_inspection(
        inspect_nodes, embed_per_node_uncertainty, custom_unregistered_grid
    )
    if manual is not None:
        return manual
    pinned = needs_pinned_operation(certified_deliverable, reproduce_across_versions)
    if pinned is not None:
        return pinned
    return ApiChoice(
        approach="transformer",
        reason="standard best operation is acceptable; favour convenience",
    )

Step 4 — Confirm the routing on a representative case #

Feed the router the flags for a real job before you commit to an implementation. A certified NAD27→NAD83 cadastral batch, for instance, routes to the pipeline string:

choice = resolve_api(certified_deliverable=True)
print(choice.approach, "->", choice.reason)
# pipeline -> operation must be pinned and reproducible for audit

Parameter and Return Reference #

Worked Example: One Point, Three Ways #

Take a single control point in eastern Canada and move it from NAD27 (EPSG:4267) to NAD83 (EPSG:4269) three ways, then compare. All three should agree to within the grid’s own accuracy — the point of the exercise is that they do, so the choice between them is about process, not answer. Note always_xy=True throughout: PROJ’s authority axis order for geographic CRSs is latitude-then-longitude, and forgetting this is the single most common source of a silently transposed coordinate.

import pyproj

# A point near Ottawa, in longitude, latitude order (always_xy convention).
lon_nad27: float = -75.6972
lat_nad27: float = 45.4215

# --- Approach 1: high-level Transformer, PROJ picks the operation ---
searched = pyproj.Transformer.from_crs("EPSG:4267", "EPSG:4269", always_xy=True)
lon_a, lat_a = searched.transform(lon_nad27, lat_nad27)

# --- Approach 2: explicit pipeline string, we name the grid step ---
pipeline_str: str = (
    "+proj=pipeline "
    "+step +proj=hgridshift +grids=ca_nrc_ntv2_0.tif "  # named, pinned grid
    "+step +proj=noop"
)
pinned = pyproj.Transformer.from_pipeline(pipeline_str)
lon_b, lat_b = pinned.transform(lon_nad27, lat_nad27)

For the manual route we reuse the NTv2 reader described in understanding NTv2 grid shift files in Python, then apply the bilinear kernel ourselves. The function below is self-contained given a grid object exposing lat_shift/lon_shift surfaces and their extents in arc-seconds:

import numpy as np


def manual_hgridshift(
    lon_deg: float,
    lat_deg: float,
    lat_grid_sec: np.ndarray,      # latitude-shift surface, arc-seconds, float64
    lon_grid_sec: np.ndarray,      # longitude-shift surface, arc-seconds, +W
    lat_min_sec: float,
    lon_min_sec_west: float,       # NTv2 longitude is positive west
    lat_inc_sec: float,
    lon_inc_sec: float,
) -> tuple[float, float]:
    """Bilinearly interpolate an NTv2 shift and apply it, returning the shifted
    (lon, lat) in degrees. Mirrors PROJ's hgridshift arithmetic exactly."""
    lat_sec = lat_deg * 3600.0
    lon_west_sec = -lon_deg * 3600.0        # +east degrees -> +west arc-seconds
    fj = (lat_sec - lat_min_sec) / lat_inc_sec
    fi = (lon_west_sec - lon_min_sec_west) / lon_inc_sec
    j0, i0 = int(np.floor(fj)), int(np.floor(fi))
    ty, tx = fj - j0, fi - i0               # fractional cell position in [0, 1)

    def bilinear(surface: np.ndarray) -> float:
        s00 = surface[j0, i0]
        s10 = surface[j0, i0 + 1]
        s01 = surface[j0 + 1, i0]
        s11 = surface[j0 + 1, i0 + 1]
        return float(
            (1 - tx) * (1 - ty) * s00 + tx * (1 - ty) * s10
            + (1 - tx) * ty * s01 + tx * ty * s11
        )

    d_lat_sec = bilinear(lat_grid_sec)
    d_lon_sec_west = bilinear(lon_grid_sec)
    lat_out = (lat_sec + d_lat_sec) / 3600.0
    lon_out = -(lon_west_sec + d_lon_sec_west) / 3600.0   # back to +east degrees
    return lon_out, lat_out

Running all three against the same grid data yields three coordinates that agree to the sub-millimetre level. The instructive comparison is not the numbers but the artifacts each route leaves behind:

def compare_three(
    a: tuple[float, float],
    b: tuple[float, float],
    c: tuple[float, float],
) -> dict[str, float]:
    """Report the maximum pairwise coordinate spread, in metres, treating one
    degree of latitude as ~111,320 m for an order-of-magnitude check."""
    deg_to_m = 111_320.0
    lons = [a[0], b[0], c[0]]
    lats = [a[1], b[1], c[1]]
    spread_deg = max(max(lons) - min(lons), max(lats) - min(lats))
    return {"max_spread_m": round(spread_deg * deg_to_m, 4)}


# Expected: a spread far below the survey tolerance — the routes agree.
# The decision between them is governed by audit, not accuracy.

The three approaches converge on the answer; they diverge on what you can show. The searched transform leaves only its result. The pipeline string leaves a citable text record of the exact operation. The manual route leaves the four node values and their accuracies. Choose the route whose artifact your deliverable requires.

Verification and Residual Analysis #

Whichever route you pick, gate it the same way: compare the transformed coordinate against an independently surveyed monument and emit a structured record. The residual against control is the planar separation

and the gate does not care which API produced N_out, E_out — a virtue, because it lets you swap approaches without touching the acceptance test.

import json
import logging
import math

logger = logging.getLogger("api_choice")


def gate_against_control(
    out_en: tuple[float, float],
    control_en: tuple[float, float],
    approach: str,
    tolerance_m: float = 0.020,     # 20 mm horizontal cadastral gate
) -> dict[str, object]:
    """Compare a transformed easting/northing against a control monument and
    emit an audit record tagged with the API that produced it."""
    residual = math.hypot(out_en[0] - control_en[0], out_en[1] - control_en[1])
    record: dict[str, object] = {
        "approach": approach,
        "residual_m": round(residual, 4),
        "tolerance_m": tolerance_m,
        "passed": residual <= tolerance_m,
    }
    logger.info("api_choice_verify %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(
            f"{approach}: residual {residual:.4f} m exceeds {tolerance_m} m"
        )
    return record


rec = gate_against_control((5031234.561, 445678.902), (5031234.567, 445678.910),
                           "pipeline")
assert rec["passed"]

Because the record carries the approach field, an audit trail records not only whether the point met tolerance but which API computed it — the reproducibility metadata that separates a defensible cadastral pipeline from a black box. The broader submission workflow that consumes these records is developed across Batch Transformation & Automation for Cadastral Coordinate Pipelines.

Troubleshooting #

Frequently Asked Questions #

Related #

• PROJ pipeline strings vs the pyproj Transformer API — the explicit-versus-searched decision in full detail.
• When to use manual grid interpolation over pyproj — the narrow cases that justify reading the grid yourself.
• Understanding NTv2 grid shift files in Python — the reader the manual route is built on.
• Batch Transformation & Automation for Cadastral Coordinate Pipelines — the parent reference on automating and certifying transformation batches.