GDAL - Образец деформации для перевернутой системы координат оси Y (Apple MapKit) дает изображение вверх дном

Question

У меня есть изображение карты, на которой я хочу привязать геоданные и наложить на встроенную карту на iOS. Я использую GDAL программно.

Шаг 1 (OK) - Geo ссылки (работает отлично)

В настоящее время я расчет геокоординат (на 6 коэффициентов) из трех опорных точек гео ссылки на изображении, и это дает мне правильное 6 коэффициентов.

Шаг 2 (ПРОБЛЕМА) - Деформация изображение + получить геокоординат для нового преобразованного изображения

Изображение становится с ног на голову! Это связано с тем, что целевая система координат (собственная система координат Apple в MapKit) имеет инвертированную ось y, которая увеличивается по мере продвижения на юг.

Вопрос

Как я могу получить GDAL деформировать изображение правильно (и в то же время дать мне правильный геокоординат пойти с ним)?

То, что я пытался

Я изменил значение 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; 
} 

Answer 1

Вы можете сделать это с помощью GDALCreateGenImgProjTransformer. Из документации: Создайте трансформатор, который отображает исходные координаты пикселя/линии в целевые координаты с привязкой (не конечная пиксельная линия). Мы делаем это, опуская дескриптор целевого набора данных (установив его в NULL).

Другая родственная Информация для вашей задачи из gdal_alg.h: Создать образ трансформера изображения.

Эта функция создает объект преобразования, который отображает из координаты пикселя/линии на одном изображении на координаты пикселя/линии на другом изображении. Изображения могут быть привязаны к географическим привязкам в разных системах координат и могут использовать GCP для сопоставления между их координатами пикселя/линии и координатами с привязкой (в отличие от предположения по умолчанию, что их геотрансформация должна использоваться).

Этот трансформатор потенциально выполняет три конкатенированных преобразования.

Первый этап - это исходные координаты пикселя/линии исходного изображения для координат с привязкой к исходному изображению и может быть выполнен с использованием геотрансформации или если не определен с использованием полиномиальной модели, полученной из опорных точек. Если используются GCP, этот этап выполняется с использованием GDALGCPTransform().

Второй этап - изменение проекций из исходной системы координат в целевую систему координат, если они отличаются. Это выполняется внутренне с помощью GDALReprojectionTransform().

Третий этап - преобразование координат привязки координат назначения к координатам места назначения. Это делается с использованием геотрансформации целевого изображения или, если оно недоступно, с использованием полиномиальной модели, полученной из GCP. Если используются GCP, этот этап выполняется с использованием GDALGCPTransform(). Этот этап пропускается, если hDstDS имеет значение NULL, когда создается преобразование.