代码改变世界

[Gdal-dev]用GCPs纠正影像的一些资源

2007-09-11 08:46  flyingfish  阅读(2292)  评论(0)  编辑  收藏  举报
从GDAL邮件列表中收集的一些资料.

[Gdal-dev] How to correct use GDALGCPTransform

Ivan Mjartan ivan.mjartan at geovap.cz
Wed Jul 25 03:30:53 EDT 2007

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;

}

More information about the Gdal-dev mailing list

[Gdal-dev] How to correct use GDALGCPTransform

Frank Warmerdam warmerdam at pobox.com
Wed Jul 25 09:46:04 EDT 2007 
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
More information about the Gdal-dev mailing list