Why is my custom transformation pipeline not applied with GDALReprojectImage?

Question

I am reprojecting a GeoTIFF image using GDAL's GDALReprojectImage while applying a custom transformation pipeline defined using PROJ's +proj=pipeline. The code produces an output image, but the transformation defined in the pipeline is not applied—the output aligns with the default reprojection instead.

The pipeline is correctly parsed, as seen in the debug logs (PROJ_TRACE), but the transformation steps are ignored. When testing the pipeline with gdalwarp, the results are correct.

Here’s the custom pipeline I am trying to apply (no need to analyze it, it is correct):

+proj=pipeline
+step +inv +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 +alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
+step +inv +proj=hgridshift +grids=Slovakia_JTSK03_to_JTSK.gsb
+step +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 +alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
+step +inv +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 +alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
+step +proj=push +v_3
+step +proj=cart +ellps=bessel
+step +proj=helmert +x=485.021 +y=169.465 +z=483.839 +rx=-7.786342 +ry=-4.397554 +rz=-4.102655 +s=0 +convention=coordinate_frame
+step +inv +proj=cart +ellps=WGS84
+step +proj=pop +v_3
+step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84

I am using gdal-sys in Rust, but the implementation in C would be very similar. Here’s the relevant part of my code:

unsafe {
    let source_wkt = CString::new("EPSG:8353").unwrap();
    let target_wkt = CString::new("EPSG:3857").unwrap();

    let warp_options = GDALCreateWarpOptions();

    let option1 = CString::new(format!("COORDINATE_OPERATION={pipeline}")).unwrap();
    let mut options: Vec<*mut i8> = vec![option1.into_raw(), ptr::null_mut()];

    let gen_img_proj_transformer = GDALCreateGenImgProjTransformer2(
        source_ds.c_dataset(),
        target_dataset.c_dataset(),
        options.as_mut_ptr(),
    );

    if gen_img_proj_transformer.is_null() {
        panic!("Failed to create image projection transformer");
    }

    (*warp_options).pTransformerArg = gen_img_proj_transformer;
    (*warp_options).pfnTransformer = Some(GDALGenImgProjTransform);
    (*warp_options).eResampleAlg = GDALResampleAlg::GRA_Lanczos;
    (*warp_options).hSrcDS = source_ds.c_dataset();
    (*warp_options).hDstDS = target_dataset.c_dataset();

    let warp_result = GDALReprojectImage(
        source_ds.c_dataset(),
        source_wkt.as_ptr(),
        target_dataset.c_dataset(),
        target_wkt.as_ptr(),
        GDALResampleAlg::GRA_Lanczos,
        0.0,
        0.0,
        None,
        ptr::null_mut(),
        warp_options,
    );

    drop(CString::from_raw(options[0]));
    GDALDestroyWarpOptions(warp_options);
    GDALDestroyGenImgProjTransformer(gen_img_proj_transformer);

    assert!(warp_result == CPLErr::CE_None, "Reprojection failed");
}

Debugging Details

The pipeline is correctly parsed, as seen in the PROJ debug logs:

PROJ_TRACE: pipeline: Pipeline: init - inv, 9
PROJ_TRACE: pipeline:     proj=krovak
PROJ_TRACE: pipeline:     lat_0=49.5
PROJ_TRACE: pipeline:     lon_0=24.8333333333333
...
PROJ_TRACE: pipeline: Pipeline: Building arg list for step no. 1
PROJ_TRACE: pipeline: Pipeline: init - inv, 3
PROJ_TRACE: pipeline:     proj=hgridshift
PROJ_TRACE: pipeline:     grids=Slovakia_JTSK03_to_JTSK.gsb
PROJ_TRACE: hgridshift: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
PROJ_TRACE: hgridshift: pj_ellipsoid - final:    ellps=GRS80

I’ve confirmed that gdalwarp with the same parameters (e.g., -s_srs, -t_srs, -ct) produces the correct output:

gdalwarp -s_srs EPSG:8353 -t_srs EPSG:3857 -ct "+proj=pipeline ..." source.tif warped.tif

No errors are produced during execution, either in logs or as return values.

---

Question:

Why is the custom transformation pipeline being ignored in GDALReprojectImage? Is there an additional step or configuration required to ensure the pipeline is applied correctly? Should I use another GDAL function, such as GDALWarp, for this purpose?

I’m looking for insights from GDAL developers or users who have experience with custom transformations using PROJ pipelines. Any guidance would be greatly appreciated.

Answer 1

After some investigation and experimenting with the GDAL API, I found that using GDALReprojectImage directly wasn't sufficient for applying a custom transformation pipeline (+proj=pipeline). The function internally overrides the pfnTransformer, which prevents the custom pipeline from being applied.

To solve the issue, I used GDALCreateWarpOperation and GDALChunkAndWarpImage to perform the reprojection manually. Here's the working code in Rust:

unsafe {
    // Create warp options
    let warp_options = GDALCreateWarpOptions();

    // Define the custom pipeline
    let option1 = CString::new(format!("COORDINATE_OPERATION={pipeline}")).unwrap();
    let mut options: Vec<*mut i8> = vec![option1.into_raw(), ptr::null_mut()];

    // Create the transformer
    let gen_img_proj_transformer = GDALCreateGenImgProjTransformer2(
        source_ds.c_dataset(),
        target_ds.c_dataset(),
        options.as_mut_ptr(),
    );
    if gen_img_proj_transformer.is_null() {
        panic!("Failed to create image projection transformer");
    }

    // Configure warp options
    (*warp_options).pTransformerArg = gen_img_proj_transformer;
    (*warp_options).pfnTransformer = Some(GDALGenImgProjTransform);
    (*warp_options).eResampleAlg = GDALResampleAlg::GRA_Lanczos;
    (*warp_options).hSrcDS = source_ds.c_dataset();
    (*warp_options).hDstDS = target_ds.c_dataset();
    (*warp_options).nDstAlphaBand = 0;
    (*warp_options).nSrcAlphaBand = 0;

    // Initialize default band mapping
    GDALWarpInitDefaultBandMapping(warp_options, 3);

    // Create the warp operation
    let warp_operation = GDALCreateWarpOperation(warp_options);
    assert!(!warp_operation.is_null(), "Failed to create warp operation");

    // Perform the warp on the full image
    let result = GDALChunkAndWarpImage(
        warp_operation,
        0,                                // Start X
        0,                                // Start Y
        GDALGetRasterXSize(target_ds),    // Width
        GDALGetRasterYSize(target_ds),    // Height
    );
    assert!(
        result == CPLErr::CE_None,
        "ChunkAndWarpImage failed with error code: {:?}",
        result
    );

    // Clean up
    GDALDestroyWarpOperation(warp_operation);
    GDALDestroyWarpOptions(warp_options);
    GDALDestroyGenImgProjTransformer(gen_img_proj_transformer);
    drop(CString::from_raw(options[0]));
}