[Gdal-dev] Problems when warping to visualise with ERDAS
Juan Cicuendez
jicicuendez at gmv.es
Tue Jul 27 05:22:25 EDT 2004
Thanks for your help Frank,
Sorry to post this again but the followup option didn't seem to work. Testing
your suggestion to fix the interleave issue in order to visualise the files
with ERDAS I tried the following option for the GeoTiff generation:
char *pszOptions = "INTERLEAVE=PIXEL"
hDstDS = GDALCreate( hDriver, outputfilename, nPixels, nLines, bands, eDT,
&pszOptions );
But it didn't work. How can I pass that options to the creation method? I also
tried to find documentation about the GeoTiff driver but I could not find it?
is it available anywhere?
I also put up all the code:
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
More information about the Gdal-dev
mailing list