问题描述
我目前正在使用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之前计算目标帧中有用像素的掩码,但我没有找到它。 我非常感谢您在此问题上的帮助!
最诚挚的问候
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)