Dst帧未严格包含在Src帧中时，GDALWarpRegionToBuffer和切片

如何解决Dst帧未严格包含在Src帧中时，GDALWarpRegionToBuffer和切片

我目前正在使用gdal api C / C ++，并且遇到gdal扭曲区域以缓冲功能（WarpRegionToBuffer）的问题。

当我的目标数据集未严格包含在我的源数据集的框架中时，不应使用随机数据填充不应包含任何数据值的区域（请参见随附的out_code.tif）。但是，也使用WarpRegionToBuffer函数的gdalwarp命令行功能似乎没有此问题。

1 /这是我使用的代码：

#include <iostream>
#include <string>
#include <vector>

#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,ytl,-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,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,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;
}

结果： output image of previous sample of code (as png,as I can't enclose TIF img)

GdalWarp命令行：

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

命令行结果： output image of previous command line (as png,as I can't enclose TIF img)

您能否帮助我发现我使用GDAL C / C ++ API出了什么问题，以使其具有与gdalwarp命令行类似的行为？ gdalwarp中可能有一种算法，可以在调用WarpRegionToBuffer之前计算目标帧中有用像素的掩码，但我没有找到它。 我非常感谢您在此问题上的帮助！

最诚挚的问候