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 ) );
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
