出处:GDAL邮件列表。下边是按照C格式写的结构比较清晰足以说明问题,就不再贴自己写的了。
[Gdal-dev] How to correct use GDALGCPTransform
Ivan Mjartan
ivan.mjartan at geovap.cz
Thu Aug 9 05:42:59 EDT 2007
- Previous message: [Gdal-dev] OGR Sqlite driver
- Next message: [Gdal-dev] How to correct use GDALGCPTransform
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
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)
{
GDALDatasetH hSrcDS,hVDS; double adfDstGeoTransform[6]; hDriver = GDALGetDriverByName( "VRT" ); hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
GDALDriverH hDriver; 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/
- Previous message: [Gdal-dev] OGR Sqlite driver
- Next message: [Gdal-dev] How to correct use GDALGCPTransform
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the Gdal-dev mailing list
注:鉴于关注此问题的朋友较多,上传一些段供参考。