[Gdal-dev] How to correct use GDALGCPTransform

[Ivan Mjartan]

Hi list,

So I have scan of raster without georeference and I need make geo
position of this raster into SJTSK coordinate (*EPSG:102067*). I need do it 
on the fly. So I am trying use GDALCreateGCPTransformer.

My input is raster it can be tiff, bitmap, jpg and so on . output can be 
Geotiff

In my test I use polynomial order 1 and 3 GCPs points.

I followed tutorial on GDAL page. But I am not success I am doing something 
wrong and I not able find right way.

In the beginning I created blank dataset to get GeoTransform and output size 
of image like in gedalwarp.exe it is correct and I get black Geotiff and 
coordinate and size is correct. (The same when I use gedal_translate.exe to 
set GCPs and than gdalwarp to warp image). But when try use 
ChunkAndWarpImage with same Transformer the output image is still blank.  So 
I also tried debug this method and I find that inside this method is used 
WarpRegion Method but there was wrong parameters of destination image. So I 
try use WarpRegion directly without success L

I thing that I am wrong using GCPTransformer I thing that I have to set 
something in GDALWarpOptions but what? This is for me big mystery. I was 
looking on internet but there is nothing L

Can I see some sample of using GCPTransformer ?

May be I can use GenImgProjTransformer but I don't know how to set GCPs with 
out prior transform to Gtiff set GCPs and saving to Disk. (Input rasters can 
be jpg bitmap .)


So like input I used this image http://80.188.225.114/temp/input.tif

I am using GDAL 1.4.2





Thanks a lot for your help



Ivan Mjartan


the code looks like follow:

Method to get destination dataset



GDALDatasetH GDALWarpCreateOutput( char *papszSrcFiles,char *pszFilename,

                      const char *pszFormat,GDAL_GCP pasGCPList[3],int 
nReqOrder,int bReversed,int nGCPCount)

{

    GDALDriverH hDriver;

    GDALDatasetH hDstDS;

    void *hTransformArg;

    GDALColorTableH hCT = NULL;

    int nDstBandCount = 0;

                        GDALDataType eDT;

                        int nOrder = 1;





                        hDriver = GDALGetDriverByName( pszFormat );



        GDALDatasetH hSrcDS;



        hSrcDS = GDALOpen( papszSrcFiles, GA_ReadOnly );



                        if( hSrcDS == NULL )

            exit( 1 );



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



        nDstBandCount = GDALGetRasterCount(hSrcDS);



        hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );



                        if( hCT != NULL )

        {

            hCT = GDALCloneColorTable( hCT );

        }



                        hTransformArg = GDALCreateGCPTransformer (

                                   nGCPCount,

                                   pasGCPList,

                                   nReqOrder,

                                   bReversed

                        );



        if( hTransformArg == NULL )

            return NULL;



        double adfThisGeoTransform[6];

        int    nThisPixels, nThisLines;



                        CPLErr eErr;



        eErr = GDALSuggestedWarpOutput( hSrcDS,GDALGCPTransform, 
hTransformArg,adfThisGeoTransform,&nThisPixels, &nThisLines);



                        /*if (eErr == CE_None)

                        {

            return NULL;

                        }*/



                        GDALDestroyGCPTransformer(hTransformArg);



                        GDALClose( hSrcDS );



                        hDstDS = GDALCreate( hDriver, pszFilename, 
nThisPixels, nThisLines,

                         nDstBandCount, eDT, NULL);



                        if( hDstDS == NULL )

                                   return NULL;



                        GDALSetProjection( hDstDS, "");

                        GDALSetGeoTransform( hDstDS, adfThisGeoTransform );



                        if( hCT != NULL )

                        {

                                   GDALSetRasterColorTable( 
GDALGetRasterBand(hDstDS,1), hCT );

                                   GDALDestroyColorTable( hCT );

                        }



                        return hDstDS;

}


and my main method looks:


int main( int argc, char ** argv )



{

    GDALDatasetH       hDstDS;

    const char         *pszFormat = "GTiff";

    void               *hTransformArg;

            int i;



    GDALDataType        eOutputType = GDT_Unknown, eWorkingType = 
GDT_Unknown;

    GDALResampleAlg     eResampleAlg = GRA_NearestNeighbour;



            char               **papszWarpOptions = NULL;

            //moje definice



            char *pszSrcFilename = NULL;

            char *pszDstFilename = NULL;

            GDAL_GCP pasGCPList[3];



            int nReqOrder = 1;

            int bReversed = 0;

            int nGCPCount = 3;



/* -------------------------------------------------------------------- */

/*      input and output rasters                         */

/* -------------------------------------------------------------------- */

            //pszSrcFilename = 
CPLStrdup("c:\\Projects\\GDAL\\gdal142\\gedalwarp\\gedalwarp\\Debug\\out_beroGCP.tif");

            //pszSrcFilename = 
CPLStrdup("c:\\Projects\\GDAL\\gdal142\\gedalwarp\\gedalwarp\\Debug\\21.tif");

            pszSrcFilename = 
CPLStrdup("c:\\Projects\\GDAL\\gdal142\\gedalwarp\\gedalwarp\\Debug\\21.tif");

            pszDstFilename = 
CPLStrdup("c:\\Projects\\GDAL\\gdal142\\gedalwarp\\gedalwarp\\Debug\\out_r.tif");



    GDALAllRegister();





            //gcplist

            pasGCPList[0].dfGCPPixel = 0;

            pasGCPList[0].dfGCPLine = 0;

            pasGCPList[0].dfGCPX = -863909;

            pasGCPList[0].dfGCPY = -993283;

            pasGCPList[0].dfGCPZ = 0;

            pasGCPList[0].pszId = "1";

            pasGCPList[0].pszInfo = "info 1";



            pasGCPList[1].dfGCPPixel = 0;

            pasGCPList[1].dfGCPLine = 3403;

            pasGCPList[1].dfGCPX = -863909;

            pasGCPList[1].dfGCPY = -1024899;

            pasGCPList[1].dfGCPZ = 0;

            pasGCPList[1].pszId = "2";

            pasGCPList[1].pszInfo = "info 2";



            pasGCPList[2].dfGCPPixel = 3800;

            pasGCPList[2].dfGCPLine = 3403;

            pasGCPList[2].dfGCPX = -809988;

            pasGCPList[2].dfGCPY = -1024899;

            pasGCPList[2].dfGCPZ = 0;

            pasGCPList[2].pszId = "3";

            pasGCPList[2].pszInfo = "info 3";



            hDstDS = GDALWarpCreateOutput( pszSrcFilename, 
pszDstFilename,pszFormat,pasGCPList,nReqOrder,bReversed,nGCPCount);



    if( hDstDS == NULL )

        exit( 1 );



    GDALDatasetH hSrcDS;



    hSrcDS = GDALOpen(pszSrcFilename, GA_ReadOnly );



    if( hSrcDS == NULL )

         exit( 2 );



/* -------------------------------------------------------------------- */

/*      Create a transformation object                                */

/* -------------------------------------------------------------------- */

        hTransformArg = GDALCreateGCPTransformer (

                                   nGCPCount,

                                   pasGCPList,

                                   nReqOrder,

                                   bReversed

                        );



        if( hTransformArg == NULL )

            return NULL;



/* -------------------------------------------------------------------- */

/*      Setup warp options.                                             */

/* -------------------------------------------------------------------- */

        GDALWarpOptions *psWO = GDALCreateWarpOptions();

                        psWO->pfnProgress = GDALTermProgress;

                        psWO->eWorkingDataType = 
GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

        //psWO->eResampleAlg = eResampleAlg;



        psWO->hSrcDS = hSrcDS;

        psWO->hDstDS = hDstDS;



        psWO->pfnTransformer = GDALGCPTransform;

        psWO->pTransformerArg = hTransformArg;



/* -------------------------------------------------------------------- */

/*      Setup band mapping.                                             */

/* -------------------------------------------------------------------- */



                        psWO->nBandCount = GDALGetRasterCount(hSrcDS);



        psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));

        psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));



        for( i = 0; i < psWO->nBandCount; i++ )

        {

            psWO->panSrcBands[i] = i+1;

            psWO->panDstBands[i] = i+1;

        }



/* -------------------------------------------------------------------- */

/*      Initialize and execute the warp.                                */

/* -------------------------------------------------------------------- */

        GDALWarpOperation oWO;



        if( oWO.Initialize( psWO ) == CE_None )

        {

                                   //if (oWO.ChunkAndWarpImage( 0, 0, 
GDALGetRasterXSize( hDstDS ),GDALGetRasterYSize( hDstDS ) )==CE_Failure){

                                               //printf("Error");



                                   //}

                                   oWO.WarpRegion( 0, 0, 
GDALGetRasterXSize( hDstDS ),GDALGetRasterYSize( hDstDS ),0, 0, 
GDALGetRasterXSize( hSrcDS ),GDALGetRasterYSize( hSrcDS ));

        }



/* -------------------------------------------------------------------- */

/*      Cleanup                                                         */

/* -------------------------------------------------------------------- */



                        if( hTransformArg != NULL ){

                                   GDALDestroyGCPTransformer( 
hTransformArg );

                        }



                        GDALClose( hDstDS );

                        GDALClose( hSrcDS );

        GDALDestroyWarpOptions( psWO );

                        GDALDestroyDriverManager();



    return 0;

}