GDAL - Warp pattern for the inverted Y axis coordinate system (Apple MapKit) gives the image upside down

Question

I have an image of a map that I want to bind to geolocation and overlay on the built-in map on iOS. I use GDAL programmatically.

Step 1 (OK) - geo-referencing (works great)

I am currently calculating a GeoTransform (6 coefficients) from three ground control points to georealize the image, and this gives me the correct 6 coefficients.

Step 2 (PROBLEM) - Warp Image + Get GeoTransform for the new converted image

Image is upside down! This is because the target coordinate system (Apple’s own coordinate system in MapKit) has an inverted y axis, which increases as you move south.

Question

How can I get GDAL to warp the image correctly (and at the same time give me the correct GeoTransform to go with it)?

What i tried

I changed the 5/6 value in the original GeoTransform to warp. This gives the correct image deformation, but the new GeoTransform is erroneous.

CURRENT CODE

- (WarpResultC*)warpImageWithGeoTransform:(NSArray<NSNumber*>*)geoTransformArray sourceFile:(NSString*)inFilepath destinationFile:(NSString*)outFilepath
{
    GDALAllRegister();

    GDALDriverH hDriver;
    GDALDataType eDT;
    GDALDatasetH hDstDS;
    GDALDatasetH hSrcDS;

    // Open the source file.
    hSrcDS = GDALOpen( inFilepath.UTF8String, GA_ReadOnly );
    CPLAssert( hSrcDS != NULL );

    // Set the GeoTransform on the source image
    // HERE IS WHERE I NEED NEGATIVE VALUES OF 4 & 5 TO GET A PROPER IMAGE
    double geoTransform[] = { geoTransformArray[0].doubleValue, geoTransformArray[1].doubleValue, geoTransformArray[2].doubleValue, geoTransformArray[3].doubleValue, -geoTransformArray[4].doubleValue, -geoTransformArray[5].doubleValue };
    GDALSetGeoTransform(hSrcDS, geoTransform);

    // Create output with same datatype as first input band.
    eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));

    // Get output driver (GeoTIFF format)
    hDriver = GDALGetDriverByName( "GTiff" );
    CPLAssert( hDriver != NULL );

    // Create a transformer that maps from source pixel/line coordinates
    // to destination georeferenced coordinates (not destination
    // pixel line).  We do that by omitting the destination dataset
    // handle (setting it to NULL).
    void *hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, NULL, NULL, NULL, FALSE, 0, 1 );
    CPLAssert( hTransformArg != NULL );

    // Get approximate output georeferenced bounds and resolution for file.
    double adfDstGeoTransform[6];
    int nPixels=0, nLines=0;
    CPLErr eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );
    CPLAssert( eErr == CE_None );
    GDALDestroyGenImgProjTransformer( hTransformArg );

    // Create the output file.
    hDstDS = GDALCreate( hDriver, outFilepath.UTF8String, nPixels, nLines, 4, eDT, NULL );
    CPLAssert( hDstDS != NULL );

    // Write out the projection definition.
    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

    // Copy the color table, if required.
    GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS, 1) );
    if( hCT != NULL )
        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS, 1), hCT );

    // Setup warp options.
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = hDstDS;

    /* -------------------------------------------------------------------- */
    /*      Do we have a source alpha band?                                 */
    /* -------------------------------------------------------------------- */
    bool enableSrcAlpha = GDALGetRasterColorInterpretation( GDALGetRasterBand(hSrcDS, GDALGetRasterCount(hSrcDS) )) == GCI_AlphaBand;
    if(enableSrcAlpha) { printf( "Using band %d of source image as alpha.\n", GDALGetRasterCount(hSrcDS) ); }

    /* -------------------------------------------------------------------- */
    /*      Setup band mapping.                                             */
    /* -------------------------------------------------------------------- */
    if(enableSrcAlpha)
        psWarpOptions->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
    else
        psWarpOptions->nBandCount = GDALGetRasterCount(hSrcDS);

    psWarpOptions->panSrcBands = (int *) CPLMalloc(psWarpOptions->nBandCount*sizeof(int));
    psWarpOptions->panDstBands = (int *) CPLMalloc(psWarpOptions->nBandCount*sizeof(int));

    for( int i = 0; i < psWarpOptions->nBandCount; i++ )
    {
        psWarpOptions->panSrcBands[i] = i+1;
        psWarpOptions->panDstBands[i] = i+1;
    }

    /* -------------------------------------------------------------------- */
    /*      Setup alpha bands used if any.                                  */
    /* -------------------------------------------------------------------- */
    if( enableSrcAlpha )
        psWarpOptions->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);

    psWarpOptions->nDstAlphaBand = GDALGetRasterCount(hDstDS);

    psWarpOptions->pfnProgress = GDALTermProgress;

    // Establish reprojection transformer.
    psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( hSrcDS, NULL, hDstDS, NULL, FALSE, 0.0, 1 );
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;
    oOperation.Initialize( psWarpOptions );
    CPLErr warpError = oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );
    CPLAssert( warpError == CE_None );

    GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );
    GDALClose( hDstDS );
    GDALClose( hSrcDS );

    WarpResultC* warpResultC = [WarpResultC new];
    warpResultC.geoTransformValues = @[@(adfDstGeoTransform[0]), @(adfDstGeoTransform[1]), @(adfDstGeoTransform[2]), @(adfDstGeoTransform[3]), @(adfDstGeoTransform[4]), @(adfDstGeoTransform[5])];
    warpResultC.newX = nPixels;
    warpResultC.newY = nLines;

    return warpResultC;
}

Answer 1

You do this using the GDALCreateGenImgProjTransformer . From the documentation: Create a transformer that maps the source pixel / line coordinates to the target coordinates with reference (not the end pixel line). We do this by omitting the handle to the target dataset (setting it to NULL).

gdal_alg.h: .

, / / . GCP / ( , ).

.

- / , GCP. GCP, GDALGCPTransform().

- , . GDALReprojectionTransform().

- . , , , GCP. GCP, GDALGCPTransform(). , hdstDS NULL .