[gdal-dev] transform input image to UTM

[Jorge Martin]

Hello,

       I want to project an input Lat/Long image into UTM. I need to do this
automatic, so I have to do it in C++. I have found a source code to
implement this issue. I have paste below the source code that I have used.
     I When this code is running, appears some times the followin error
messages:

       "ERROR 1: no system list, errno: 42" and
      "ERROR 6: SetColorTable() not supported for multi-sample TIFF files."

     What is it means?
When I open the output image there are some strange black lines crossing the
output image that not appear in the input image. Is there any reason for
this behaviour?


Many thanks,

Jorge


//----------------------------------SOUCE
CODE----------------------------------------


void Reproject()
{
 // ------------ Example of initializing an image into the GDALDataset
object

    GDALDataset  *poDataset;

    // -------------- Reprojection example

 GDALDatasetH  hSrcDS, hDstDS; //, hRefDS;

    // Open input and output files.

    GDALAllRegister();

    hSrcDS = GDALOpen(
"/home/jjmf/sample_images/t1.09155.AV_MI.143.250m.tif", GA_ReadOnly );

    poDataset = (GDALDataset *) hSrcDS;
    if( poDataset == NULL )
    {

        cout<<"The input file does not exist: "<< endl;

        exit(1);
    }

    GDALDriverH hDriver;
    GDALDataType eDT;

    // Open the source file.

    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.

    const char *pszSrcWKT = NULL;
    char *pszDstWKT;

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



    OGRSpatialReference oSRS;

    oSRS.SetUTM(12, TRUE);

    oSRS.SetWellKnownGeogCS("EPSG:4258");
    //oSRS.importFromProj4("+proj=latlo +datum=WGS84");

    oSRS.exportToWkt( &pszDstWKT );

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

    GDALDestroyGenImgProjTransformer( hTransformArg );

    // Create the output file.
    char *apszOptions[] = {"INTERLEAVE=PIXEL", NULL};

    hDstDS = GDALCreate(hDriver, "/home/jjmf/sample_images/out_new.tif",
nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, apszOptions);

    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));
    GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);


    hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,2));
    GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,2), hCT);


    hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,3));
    GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,3), hCT);


 cout << "Setup warp options ..."<<endl;
    // Setup warp options.

    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = hDstDS;
    psWarpOptions->nBandCount = 3;
    psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) *
psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) *
psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;

    psWarpOptions->pfnProgress = GDALTermProgress;

 cout << "Establishing transformer ..."<<endl;
    // 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.

 cout << "Executing ..."<<endl;
    GDALWarpOperation oOperation;

    // Initialize transformation operation
 oOperation.Initialize(psWarpOptions);

 // Perform the transformation
    oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS),
GDALGetRasterYSize(hDstDS));

 // Clean up the intermediate objects
    GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
    GDALDestroyWarpOptions(psWarpOptions);


 cout <<"Done ..."<<endl;

 // Close the input and output files
    GDALClose(hDstDS);
    GDALClose(hSrcDS);



}

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