[Gdal-dev] Re: Problems when warping

[Juan Cicuendez]

Frank Warmerdam <warmerdam <at> pobox.com> writes:

> 
> Juan Cicuendez wrote:
> > Hi everybody,
> > 
> > I am trying to change the projection of an RGB scanned map. This is my 
first 
> > experience with GDAL and I followed the instructions of the Warping 
tutorial. 
> > In order to process all the three bands I set psWarpOptions->nBandCount = 0;
> > I finally obtained the geotiff file but when I try to visualised it from 
ERDAS 
> > it only recognises the first band, however if I try from ENVI it goes fine 
an 
> > sees all the bands. Looking at the image with windows Imaging application, 
it 
> > looks very strange and this could be a problem with the interleaving of the 
> > data. Has anybody come across this problem before? Does anybody know how to 
fix 
> > this problems in order that the files can be visualised from ERDAS?
> > 
> > I would appreciate if anyone could give me hand, I tried many things but 
> > nothing seemed to work.
> 
> Juan,
> 
> By default GDAL creates TIFF files with PLANARCONFIG_SEPERATE (band 
interleaved)
> but I have observed some versions of Erdas tools and other applications having
> problems with this.  You can pass the INTERLEAVE=PIXEL creation option to the
> GTiff driver when creating a new TIFF file to create a pixel interleaved file.
> 
> Of course, if you could provide a smallish sample file demonstrating the
> problem it would be easier to determine if this is the issue.
> 
> Best regards,
> 
Thanks for your help Frank,

Testing your suggestion. I tried something like:

char *pszOptions = "INTERLEAVE=PIXEL"
hDstDS = GDALCreate( hDriver, outputfilename, nPixels, nLines, bands, eDT, 
&pszOptions ); 

, but of course it didn't work. How can I pass the opticon INTERLEAVE=PIXEL 
when creating the geoTiff file? I haven't found the documentation of the 
GeoTiff driver. Where can I find it?

I put up all the code of the function:

CPLErr change_projection(char *filename, char *outputfilename, char *pszDstWKT)
{

    GDALDriverH hDriver;
    GDALDataType eDT;
    GDALDatasetH hDstDS;
    GDALDatasetH hSrcDS;

    // Open the source file. 
    hSrcDS = GDALOpen( filename, GA_ReadOnly );
    CPLAssert( hSrcDS != NULL );
    
    // Create output with same datatype as first input band. 

    eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

    // Get output driver (GeoTIFF format)

    hDriver = GDALGetDriverByName( "GTiff" );
    CPLAssert( hDriver != NULL );

    // Get Source coordinate system. 

    char *pszSrcWKT = NULL;

    pszSrcWKT = (char *) GDALGetProjectionRef( hSrcDS );
    CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );


    // Create a transformer that maps from source pixel/line coordinates
    // to destination georeferenced coordinates (not destination 
    // pixel line).  We do that by omitting the destination dataset
    // handle (setting it to NULL). 

    void *hTransformArg;

    hTransformArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, 
FALSE, 0.0, 1 );
    CPLAssert( hTransformArg != NULL );

    // Get approximate output georeferenced bounds and resolution for file. 

    double adfDstGeoTransform[6];
    int nPixels=0, nLines=0;
    CPLErr eErr;

    eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, 
hTransformArg, 
                                                adfDstGeoTransform, &nPixels, 
&nLines );
    CPLAssert( eErr == CE_None );

    // Create the output file.
    int bands = GDALGetRasterCount(hSrcDS);
    char *pszOptions = "INTERLEAVE=PIXEL"

    hDstDS = GDALCreate( hDriver, outputfilename, nPixels, nLines, bands, eDT, 
&pszOptions ); 
    
    CPLAssert( hDstDS != NULL );

    // 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 );

    GDALClose( hDstDS );

    hDstDS = (GDALDataset *) GDALOpen( outputfilename, GA_Update );

    
    // Setup warp options. 
    
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = hDstDS;

    psWarpOptions->nBandCount = 0;

    psWarpOptions->pfnProgress = GDALTermProgress;   

    // Establish reprojection transformer. 

    psWarpOptions->pTransformerArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, 
                                         pszSrcWKT, 
                                         hDstDS,
                                         pszDstWKT, 
                                         FALSE, 0.0, 1 );

    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation. 

    GDALWarpOperation oOperation;

    oOperation.Initialize( psWarpOptions );
    oOperation.ChunkAndWarpImage( 0, 0, 
                                  GDALGetRasterXSize( hDstDS ), 
                                  GDALGetRasterYSize( hDstDS ) );

    
    GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );  


    GDALClose( hDstDS );
    GDALClose( hSrcDS );


    return eErr;
}

Cheers,
Juan