GDAL - 将图像扭曲到倒置的 y 轴坐标系（Apple MapKit）使图像倒置

Posted 2023-03-23

【中文标题】GDAL - 将图像扭曲到倒置的 y 轴坐标系（Apple MapKit）使图像倒置

【英文标题】：GDAL - Warp image to inverted y-axis coordinate system (Apple MapKit) gives image upside down

【发布时间】：2016-06-24 02:11:54

【问题描述】：

我有一张地图图像，我想在 ios 的内置地图上进行地理参考和叠加。我正在以编程方式使用 GDAL。

第 1 步（确定）- 地理参考（工作正常）

目前我正在计算来自三个地面控制点的 GeoTransform（6 个系数）来地理参考图像，它给了我正确的 6 个系数。

第 2 步（问题）- 扭曲图像 + 为新的变换图像获取 GeoTransform

图像变得颠倒了！这是因为目标坐标系（Apple 在 MapKit 中自己的坐标系）具有反转的 y 轴，随着您向南移动，它会增加。

问题

我怎样才能让 GDAL 正确地扭曲图像（同时给我一个正确的 GeoTransform 来配合它）？

我的尝试

我在变形之前将原始 GeoTransform 中的值更改为 5/6。这给出了正确的图像扭曲，但新的 GeoTransform 是错误的。

当前代码

- (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;

【问题讨论】：

检查this question 也许你可以翻转图像 zi = zi[::-1,:] 或者你在here这样的地方缺少一个减号

感谢您的提示。我可以翻转图像（请参阅问题中的详细信息），但我做不到，最后仍然有正确的地理变换。

【参考方案1】：

您可以使用GDALCreateGenImgProjTransformer 来执行此操作。从文档： 创建一个转换器，将源像素/线坐标映射到目标地理参考坐标（不是目标像素线）。我们通过省略目标数据集句柄（将其设置为 NULL）来做到这一点。

来自gdal_alg.h的与您的任务相关的其他信息： 创建图像到图像转换器。

此函数创建一个转换对象，该对象将一个图像上的像素/线坐标映射到另一图像上的像素/线坐标。这些图像可能会在不同的坐标系中进行地理参考，并且可以使用 GCP 在其像素/线坐标和地理参考坐标之间进行映射（与应使用其地理变换的默认假设相反）。

此转换器可能执行三个级联转换。

第一阶段是从源图像像素/线坐标到源图像地理参考坐标，可以使用地理变换完成，或者如果未使用从 GCP 派生的多项式模型进行定义。如果使用 GCP，则此阶段使用 GDALGCPTransform() 完成。

第二阶段是将投影从源坐标系更改为目标坐标系，假设它们不同。这是在内部使用 GDALReprojectionTransform() 完成的。

第三阶段是从目标图像地理参考坐标转换为目标图像坐标。这是使用目标图像地理变换完成的，或者如果不可用，则使用从 GCP 派生的多项式模型。如果使用 GCP，则此阶段使用 GDALGCPTransform() 完成。 如果在创建转换时 hDstDS 为 NULL，则跳过此阶段。

【讨论】：

感谢您的回答，但这就是我正在做的。我正在使用您链接到的那个页面作为基础。问题仍然是由于地理参考坐标系颠倒了，图像被颠倒了。

您是否尝试将 NULL 设置为 hDstDS？ hTransformArg = GDALCreateGenImgProjTransformer(hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, FALSE, 0, 1);

感谢您的帮助，但我不明白您所写的内容是否对我有帮助。我在问题中添加了一些代码，以便您了解我今天所做的事情，那么解释您的意思可能会更容易？谢谢！