[gdal-dev] Warp results in wrong colours

[Nik Sands]

My application is doing a warp using GDAL which works fine for some source images, but with others results in an image with the colours all wrong (eg, water areas are red instead of blue).  I'm guessing that the RGB(A) (bands?) are getting mixed up somehow, but I really don't know where to start looking after my web searches produced nothing useful.

The problem only occurs for some images.  So far the only thing that the problem source images have in common is that they were all converted from ECW to GeoTIFF by 'gdal_translate' (with no options).  The source GeoTIFF files have all the correct colours before the warp.

The code I'm using is based on the GDAL warp API tutorial and is included below.

I'd be grateful if somebody could point me in the right direction.

Thanks in advance,
Nik.



		void *hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 );
		
		if ( hTransformArg == NULL )
		{
			NSLog(@"Failed to create transformation.");
			return NO;
		}
		
		// Get approximate output georeferenced bounds and resolution for file.
		
		double adfDstGeoTransform[6];
		int nPixels=0, nLines=0;
		
		if ( GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines ) != CE_None )
		{
			NSLog(@"Failed to get suggested warp output.");
			return NO;
		}
		
		GDALDestroyGenImgProjTransformer( hTransformArg );
		
		// Create the output file.
		
		NSLog(@"source band count:  %d", GDALGetRasterCount(hSrcDS));
		
		hDstDS = GDALCreate( hDriver, [dstPath cStringUsingEncoding:NSASCIIStringEncoding], nPixels, nLines, GDALGetRasterCount(hSrcDS) + 1, eDT, NULL );
		
		if ( hDstDS == NULL )
		{
			NSLog(@"Failed to open destination file '%@' with GDALCreate().", [dstPath lastPathComponent]);
			return NO;
		}
		
		// Write out the projection definition.
		
		GDALSetProjection( hDstDS, pszDstWKT );
		GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
		
		// Copy the color table, if required.
		
		GDALColorTableH hCT;
		
		hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
		if( hCT != NULL )
			GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
		
		// 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 = reportWarpProgress;
		psWarpOptions->nDstAlphaBand = GDALGetRasterCount(hSrcDS) + 1;
		
		// 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.
		
		GDALWarpOperationH oOperation = GDALCreateWarpOperation( psWarpOptions );
		
		GDALChunkAndWarpImage( oOperation, 0, 0,
							  GDALGetRasterXSize( hDstDS ),
							  GDALGetRasterYSize( hDstDS ) );