Hi,
I am trying to write a program which is supposed to read geotif files, extract
the values I need, then export them in a homemade format.
As my geotif files are on several different UTM zones, I have to change
projections so that all the files I open have the same projection as the "main
tile" (central tile).
So i'm using WarpRegionToBuffer from GDALWarpOperation class with a null
hDstDS, which is supposed to work... but my buffer remains filed with the
default values I assigned to him (-999), even though WarpRegionToBuffer returns
CE_None...
here is my code:
float* MyTile::extractdomain(float rxso, float ryso, float rxne, float ryne,
GDALDataset *maintile)
{
int immai = 0;
int jmmai = 0;
GDALDatasetH hSrcDS, hDstDS;
GDALAllRegister();
hSrcDS = GDALOpen( this->filename.toLatin1().data(), GA_ReadOnly );
// hDstDS = GDALOpen( "out.tif", GA_Update );
hDstDS = 0;
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = 0;
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;
psWarpOptions->eWorkingDataType = GDT_CFloat32;
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( hSrcDS,
GDALGetProjectionRef(hSrcDS),
hDstDS,
maintile->GetProjectionRef(),
FALSE, 0.0, 1 );
//maintile is my central tile, i want to apply its projection to all the other
tiles i'm concatenating
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions );
if (rxso < XSO)
rxso = XSO;
if (ryso < YSO)
ryso = YSO;
if (rxne > XNE)
rxne = XNE;
if (ryne > YNE)
ryne = YNE;
immai = ABS((rxne - rxso)) / dx;
jmmai = ABS((ryne - ryso)) / dy;
try
{
this->ch2d = new float[immai * jmmai];
if (ch2d)
{
for (int i = 0 ; i < immai * jmmai ; i++)
ch2d[i] = -999;
}
}
catch(...)
{
//could not allocate memory
return (0);
}
int nXOff=int((rxso-dx/2.-rxso)/dx);
int nYOff=int((YNE-(ryne+(-dy)/2.))/(dy));
int nXSize=((rxne-rxso+dx)/dx)+1;
int nYSize=((ryne-ryso+(-dy))/(dy))+1;
int xs;
int ys;
double adfGeoTransform[6];
GDALSuggestedWarpOutput(hSrcDS, psWarpOptions->pfnTransformer,
psWarpOptions->pTransformerArg, adfGeoTransform, &xs, &ys);
int error = oOperation.WarpRegionToBuffer(0, 0, nXSize, nYSize,
ch2d,
GDT_CFloat32, nXOff, nYOff, nXSize, nYSize);
if (error == CE_None)
printf("okcool\n");
else
printf("okpascool\n");
int size = sizeof(float);
for (int i = 0 ; i < immai * jmmai ; i++)
{
if (ch2d[i] != -999)
{
//small test to see if I have != values from default...
printf("roflflflfl");
}
}
GDALClose(hSrcDS);
GDALClose(hDstDS);
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );
return (this->ch2d);
}
By the way, I noticed a floating exception in
GDALWarpOperation::WarpRegionToBuffer() if you don't set a value to your
GDALWarpOptions eWorkingDataType attribute, since nWordSize will be null and
the integer overflow test will divide by zero : "nSrcXSize * nSrcYSize >
INT_MAX / (nWordSize * psOptions->nBandCount)"
Any hint / help will be really appreciated!
Thanks,
--
Matthieu
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
