GDALWarpRegionToBuffer и мозаика, когда Dst Frame не содержится строго в Src Frame ⇐ C++

[Anonymous]

В настоящее время я работаю с gdal api C/C++ и столкнулся с проблемой, связанной с функцией преобразования области gdal в буфер (WarpRegionToBuffer).
Когда мой целевой набор данных не является строго содержащийся в кадре моего исходного набора данных, область, где не должно быть значений данных, заполняется случайными данными (см. прилагаемый out_code.tif). Однако функциональность командной строки gdalwarp, которая также использует функцию WarpRegionToBuffer, похоже, не имеет этой проблемы.
1/ Вот код, который я использую:

Код: Выделить всё

#include 
#include 
#include 

#include "gdal.h"
#include "gdalwarper.h"
#include "cpl_conv.h"

int main(void)
{
std::string pathSrc = "in.dt1";

//these datas will be provided by command line
std::string pathDst = "out_code.tif";

double resolutionx = 0.000833333;
double resolutiony = 0.000833333;

//destination corner coordinates: top left (tl) bottom right (br)

float_t xtl = -1;
float_t ytl = 45;
float_t xbr = 2;
float_t ybr = 41;

//tile size defined by user
int tilesizex = 256;
int tilesizey = 256;

float width = std::ceil((xbr - xtl)/resolutionx);
float height = std::ceil((ytl - ybr)/resolutiony);

double adfDstGeoTransform[6] = {xtl, resolutionx, 0, ytl, 0, -resolutiony};

GDALDatasetH  hSrcDS, hDstDS;

// Open input file
GDALAllRegister();
hSrcDS = GDALOpen(pathSrc.c_str(), GA_ReadOnly);
GDALDataType  eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

// Create output file, using same spatial reference as input image, but new geotransform
GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
hDstDS = GDALCreate( hDriver, pathDst.c_str(), width, height, GDALGetRasterCount(hSrcDS), eDT, NULL );

OGRSpatialReference oSRS;
char *pszWKT = NULL;

//force geo projection
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszWKT );
GDALSetProjection( hDstDS, pszWKT );
//Fetches the coefficients for transforming between pixel/line (P,L) raster space,
//and projection coordinates (Xp,Yp) space.
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

// Setup warp options
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
psWarpOptions->nBandCount = 1;
psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
psWarpOptions->panSrcBands[0] = 1;
psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->pfnProgress = GDALTermProgress;

//these datas will be calculated in order to warp tile by tile
//current tile size
int cursizex = 0;
int cursizey = 0;

double nbtilex = std::ceil(width/tilesizex);
double nbtiley = std::ceil(height/tilesizey);

int starttilex = 0;
int starttiley = 0;

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

// Initialize and execute the warp operation on region
GDALWarpOperation oOperation;
oOperation.Initialize(psWarpOptions);

for (int ty = 0; ty < nbtiley; ty++) {
//handle last tile size
//if it last tile change size otherwise keep tilesize
for (int tx = 0; tx < nbtilex;  tx++) {
//if it last tile change size otherwise keep tilesize
starttiley = ty * tilesizey;
starttilex = tx * tilesizex;

cursizex = std::min(starttilex + tilesizex, (int)width) - starttilex;
cursizey = std::min(starttiley + tilesizey, (int)height) - starttiley;

float * buffer = new float[cursizex*cursizey];
memset(buffer, 0, cursizex*cursizey);

//warp source
CPLErr ret = oOperation.WarpRegionToBuffer(
starttilex, starttiley, cursizex, cursizey,
buffer,
eDT);
if (ret != 0) {
CEA_SIMONE_ERROR(CPLGetLastErrorMsg());
throw std::runtime_error("warp error");
}

//write the fuzed tile in dest
ret = GDALRasterIO(GDALGetRasterBand(hDstDS,1),
GF_Write,
starttilex, starttiley, cursizex, cursizey,
buffer, cursizex, cursizey,
eDT,
0, 0);
if (ret != 0) {
CEA_SIMONE_ERROR("raster io write error");
throw std::runtime_error("raster io write error");
}
delete(buffer);
}
}

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

return 0;
}

Результат:
выходное изображение предыдущего примера кода (в формате png, так как я не могу вложить изображение TIF)
Командная строка GdalWarp:

Код: Выделить всё

gdalwarp -te -1 41 2 45 -tr 0.000833333 0.000833333 in.dt1 out_cmd_line.tif

Результат командной строки:
выходное изображение предыдущей командной строки (в формате png, так как я не могу вложить изображение TIF)
Можете ли вы помочь мне найти, что не так с использованием API GDAL C/C++, чтобы поведение было аналогично командной строке gdalwarp? Вероятно, в gdalwarp есть алгоритм, который вычисляет маску полезных пикселей в целевом кадре перед вызовом WarpRegionToBuffer, но я его не нашел.
Буду очень признателен за помощь в решении этой проблемы!
С уважением

Подробнее здесь: https://stackoverflow.com/questions/641 ... in-src-fra