[gdal-dev] GDALWarpOperation::WarpRegionToBuffer() doesn't write in my buffer

[bosow at club-internet.fr]

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







-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20101119/cc2b3413/attachment-0001.html