[Gdal-dev] How to correct use GDALGCPTransform

[Ivan Mjartan]

Hi, Frank and all .

Thanks for previous reply. It helps me very much.

So I was little bit testing and I found solution for me, I am using 
GDALGenImgProjTransform and virtual dataset like you said.

It works fine but I have question Is it true when I use virtual datatset 
that whole raster is loading in to memory? Maybe I have to save input raster 
into intermediate tiff file (to set GCP into input raster, it can be jpg or 
other), because my input can by very large, just when raster is loading 
whole.

Thanks very much.



If someone will need to solve similar problem my code looks like follow.

BTW: Maybe sorry I really don't know how to correct reply into my thread.


GDALDatasetH CreateVirtualDatasetFromSourceFile(GDAL_GCP *pasGCPList,int 
numPoints,char *pszSrcFilename,const char *pszFormat){
 GDALDriverH hDriver;
 GDALDatasetH hSrcDS,hVDS;
 double adfDstGeoTransform[6];
 hDriver = GDALGetDriverByName( "VRT" );

 hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );

 if( hSrcDS == NULL ){
  return NULL;
 }

 hVDS = GDALCreateCopy( hDriver,"", hSrcDS,FALSE, NULL, NULL, NULL );

 GDALClose( hSrcDS );


 //nastavim GCPs body
 GDALSetGCPs(hVDS,numPoints,pasGCPList,"");
 GDALSetProjection( hVDS, "");

 GDALGCPsToGeoTransform(numPoints,pasGCPList,adfDstGeoTransform,TRUE);
 GDALSetGeoTransform( hVDS, adfDstGeoTransform );
 return hVDS;
}


GDALDatasetH GDALWarpCreateOutput(GDALDatasetH hSrcDS, const char 
*papszDesFile, const char *pszFormat,int nOrder)
{
    GDALDriverH hDriver;
    GDALDatasetH hDstDS;
    void *hTransformArg;
    GDALColorTableH hCT = NULL;
    int nDstBandCount = 0;
 int nPixels, nLines;
 double adfDstGeoTransform[6];
 GDALDataType eDT;
 char **papszCreateOptions = NULL;
    hDriver = GDALGetDriverByName( pszFormat );


 if( hSrcDS == NULL )
        return NULL;

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


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

    if( hCT != NULL )
    {
        hCT = GDALCloneColorTable( hCT );
    }

    hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, "",NULL, 
"",TRUE, 1000.0, nOrder );

    if( hTransformArg == NULL )
        return NULL;

 GDALSuggestedWarpOutput( hSrcDS,GDALGenImgProjTransform, 
hTransformArg,adfDstGeoTransform, &nPixels, &nLines );

 GDALDestroyGenImgProjTransformer( hTransformArg );
 papszCreateOptions = CSLSetNameValue(papszCreateOptions,"TFW", "YES");

    hDstDS = GDALCreate( hDriver, papszDesFile, nPixels, 
nLines,nDstBandCount, eDT,papszCreateOptions);

    if( hDstDS == NULL )
        return NULL;

    GDALSetProjection( hDstDS, "" );
    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

    if( hCT != NULL )
    {
        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
        GDALDestroyColorTable( hCT );
    }

    return hDstDS;
}

C_GDAL4_API int warpImage(double* rasterXArray,double* rasterYArray,double* 
realXArray,double* realYArray, int numPoints,char *pszSrcFilename,char 
*pszDstFilename)
{
    GDALDatasetH hDstDS;
 GDALDatasetH hSrcDS;
    const char         *pszFormat = "GTiff";

    int                 i,j;
    void               *hTransformArg;
    GDALTransformerFunc pfnTransformer = NULL;
    GDALDataType        eOutputType = GDT_Unknown, eWorkingType = 
GDT_Unknown;
    GDALResampleAlg     eResampleAlg = GRA_NearestNeighbour;
 GDAL_GCP *pasGCPList = NULL;


 //create gcp list
 pasGCPList = (GDAL_GCP *) malloc((numPoints)*sizeof(GDAL_GCP) );

 for(j=0 ; j<numPoints;j++ )
 {
    pasGCPList[j].dfGCPPixel = rasterXArray[j];
    pasGCPList[j].dfGCPLine  = rasterYArray[j];
    pasGCPList[j].dfGCPX = realXArray[j];
    pasGCPList[j].dfGCPY = realYArray[j];
    pasGCPList[j].pszId = (char *) malloc( 10*sizeof(char) );
   pasGCPList[j].pszInfo = (char *) malloc( 10*sizeof(char) );
    sprintf(pasGCPList[j].pszId,"");
   sprintf(pasGCPList[j].pszInfo,"");
 }



    GDALAllRegister();

 hSrcDS = 
CreateVirtualDatasetFromSourceFile(pasGCPList,numPoints,pszSrcFilename,pszFormat);
 hDstDS = GDALWarpCreateOutput(hSrcDS, pszDstFilename,pszFormat,0);

    if( hSrcDS == NULL )
        return 1;

 if( hDstDS == NULL )
        return 1;

/* -------------------------------------------------------------------- */
/*      Create a transformation object from the source to               */
/*      destination coordinate system.                                  */
/* -------------------------------------------------------------------- */
    hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, "",hDstDS, 
"",TRUE, 1000.0,0);

    if( hTransformArg == NULL )
         return 1;

    pfnTransformer = GDALGenImgProjTransform;
    GDALWarpOptions *psWO = GDALCreateWarpOptions();

    psWO->eWorkingDataType = eWorkingType;
    psWO->eResampleAlg = eResampleAlg;

    psWO->hSrcDS = hSrcDS;
    psWO->hDstDS = hDstDS;

    psWO->pfnTransformer = pfnTransformer;
    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;
    }

    GDALWarpOperation oWO;

    if( oWO.Initialize( psWO ) == CE_None )
    {
  oWO.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( 
hDstDS ),GDALGetRasterYSize( hDstDS ) );
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */

    if( hTransformArg != NULL )
        GDALDestroyGenImgProjTransformer( hTransformArg );

    GDALDestroyWarpOptions( psWO );


 //for ( j=0; i<numPoints;++j ) {
 // free( pasGCPList[j].pszId );
 // free( pasGCPList[j].pszInfo );
 //}
 //free( pasGCPList );

  GDALClose( hSrcDS );
    GDALClose( hDstDS );

    GDALDestroyDriverManager();

    return 0;
}



----- Original Message ----- 
From: "Frank Warmerdam" <warmerdam at pobox.com>
To: "Ivan Mjartan" <ivan.mjartan at geovap.cz>
Cc: <gdal-dev at lists.maptools.org>
Sent: Wednesday, July 25, 2007 3:46 PM
Subject: Re: [Gdal-dev] How to correct use GDALGCPTransform


> Ivan Mjartan wrote:
>> 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 .)
>
> Ivan,
>
> I believe your problem is that the warper needs a transformer that goes
> from input pixel/line coordinates to destination pixel/line coordinates 
> (and
> the reverse).  But the GDALGCPTransformer only does source pixel/line to
> source georeferenced coordinates.  If you look at GenImgProjTransform() in
> gdal/alg/gdaltransformer.cpp you may note it actually performs several 
> steps.
> One from source pixel/line to source georef.  Then potentially reproject
> to destination pixel/line.  Then dest georef to dest pixel/line.
>
> You could implement your own similar transformer somewhat similar to
> GDALGenImgProjTransform, but I think the easier way is just to associate
> the GCPs with the source image.  If it is impractical to write the image
> to a new TIFF file with gdal_translate, and with the GCPs associated, then
> another approach without intermediate files would be to create an 
> in-memory
> VRT dataset using GDALCreateCopy() from the source image, and then call
> SetGCPs() on that VRT dataset.
>
> Then use the VRT as an input with the GenImgProj transformer.
>
> BTW, to avoid having a VRT written to disk, just use the filename "" for
> it when calling CreateCopy().
>
> BTW, the way you included your code (tab expansion, extra blank lines)
> made it nearly unreadable.
>
> Good work on all you have already learned about the warper!
>
> Best regards,
> -- 
> ---------------------------------------+--------------------------------------
> I set the clouds in motion - turn up   | Frank Warmerdam, 
> warmerdam at pobox.com
> light and sound - activate the windows | http://pobox.com/~warmerdam
> and watch the world go round - Rush    | President OSGeo, http://osgeo.org
>
>
>