[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.czWed Jul 25 03:30:53 EDT 2007
- Previous message: [Gdal-dev] Create Polygon in OGR / C#
- Next message: [Gdal-dev] How to correct use GDALGCPTransform
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
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;
}
-
Previous message: [Gdal-dev] Create Polygon in OGR / C#
-
Next message: [Gdal-dev] How to correct use GDALGCPTransform
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the Gdal-dev mailing list
[Gdal-dev] How to correct use GDALGCPTransform
Frank Warmerdam warmerdam at pobox.comWed Jul 25 09:46:04 EDT 2007
- Previous message: [Gdal-dev] How to correct use GDALGCPTransform
- Next message: [Gdal-dev] Reproject with Mapserver vs gdalwarp
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
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] How to correct use GDALGCPTransform
- Next message: [Gdal-dev] Reproject with Mapserver vs gdalwarp
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the Gdal-dev mailing list