how to add the gcp to a tif and make it to a geotiff?
use gdal_translate -gcp we could do this.
(in gdal_translate.cpp）
else if( EQUAL(argv[i],"-gcp") && i < argc - 4 )
{
/* -gcp pixel line easting northing [elev] */
nGCPCount++;
pasGCPs = (GDAL_GCP *)
CPLRealloc( pasGCPs, sizeof(GDAL_GCP) * nGCPCount );
GDALInitGCPs( 1, pasGCPs + nGCPCount - 1 );
pasGCPs[nGCPCount-1].dfGCPPixel = atof(argv[++i]);
pasGCPs[nGCPCount-1].dfGCPLine = atof(argv[++i]);
pasGCPs[nGCPCount-1].dfGCPX = atof(argv[++i]);
pasGCPs[nGCPCount-1].dfGCPY = atof(argv[++i]);
if( argv[i+1] != NULL
&& (atof(argv[i+1]) != 0.0 || argv[i+1][0] == '0') )
pasGCPs[nGCPCount-1].dfGCPZ = atof(argv[++i]);
/* should set id and info? */
}
I wrtie in my code ：
char* szFilter = "TIFF(*.tif)|*.tif|BMP (*.bmp)|*.bmp|JPG (*.jpg)|
*.jpg|Image(*.img)|*.img||";
CFileDialog dlg(FALSE,TEXT(szFilter),NULL,OFN_HIDEREADONLY |
OFN_OVERWRITEPROMPT,TEXT(szFilter),this);
if(dlg.DoModal()==IDOK)
{
CString strSaveName = dlg.GetPathName(); //获得保存的完整路径
GDALAllRegister();
GDALDataset* poDataset;
poDataset = (GDALDataset *) GDALOpen(filename , GA_Update );
if( !poDataset == NULL )
{
VRTDataset *poVDS;//创建虚拟数据集
int nRasterXSize = poDataset->GetRasterXSize();
int nRasterYSize = poDataset->GetRasterYSize();
poVDS = (VRTDataset *) VRTCreate( nRasterXSize, nRasterYSize );
char * pszOutputSRS = "WGS_84";
poVDS->SetProjection( pszOutputSRS );//SRS设为新导入的SRS
if(m_bChangeCoords)
{
double adfGeoTransform[6];
adfGeoTransform[0] = 20; //m_dTLX;
adfGeoTransform[1] = 2; //m_dCellX;
adfGeoTransform[2] = 0;
adfGeoTransform[3] = 20; //m_dTLY;
adfGeoTransform[4] = 0;
adfGeoTransform[5] = -2; //m_dCellY;
poVDS->SetGeoTransform( adfGeoTransform );// }
GDAL_GCP * pasGCPs = NULL;
int mGcpCount = 4;
//gcplist
pasGCPs = (GDAL_GCP *)CPLRealloc( pasGCPs, sizeof(GDAL_GCP) *
mGcpCount );
pasGCPs[0].dfGCPPixel = 0;
pasGCPs[0].dfGCPLine = 0;
pasGCPs[0].dfGCPX = 863909;
pasGCPs[0].dfGCPY = 993283;
pasGCPs[0].dfGCPZ = 0;
pasGCPs[0].pszId = "1";
pasGCPs[0].pszInfo = "info 1";
pasGCPs[1].dfGCPPixel = 0;
pasGCPs[1].dfGCPLine = 201;
pasGCPs[1].dfGCPX = 863909;
pasGCPs[1].dfGCPY = 1024899;
pasGCPs[1].dfGCPZ = 0;
pasGCPs[1].pszId = "2";
pasGCPs[1].pszInfo = "info 2";
pasGCPs[2].dfGCPPixel = 201;
pasGCPs[2].dfGCPLine = 201;
pasGCPs[2].dfGCPX = 809988;
pasGCPs[2].dfGCPY = 1024899;
pasGCPs[2].dfGCPZ = 0;
pasGCPs[2].pszId = "3";
pasGCPs[2].pszInfo = "info 3";
pasGCPs[3].dfGCPPixel = 201;
pasGCPs[3].dfGCPLine = 0;
pasGCPs[3].dfGCPX = 809988;
pasGCPs[3].dfGCPY = 0;
pasGCPs[3].dfGCPZ = 0;
pasGCPs[3].pszId = "4";
pasGCPs[3].pszInfo = "info 4";
char * m_GcpProjection = "WGS_84";
poVDS->SetGCPs(mGcpCount,pasGCPs,m_GcpProjection);
if
(poVDS->SetGCPs(mGcpCount,pasGCPs,m_GcpProjection)==CE_None) //说明
setgcp是成功的,但是这里并未更改tif，
{
AfxMessageBox("ok");
}
double madfGeoTransform[6] = {444720, 2.5, 0, 3751320, 0, -2.5};/
char * src_adfGeoTransform = poDataset->GetGeoTransform();
poVDS->SetGeoTransform( madfGeoTransform);
poVDS->SetProjection(mProjection);
poVDS->SetGCPs(mGcpCount,pasGCPs,m_GcpProjection);
poVDS->SetMetadata( poDataset->GetMetadata( ) );//元数据保持旧的
char * mGcpProjection;
poVDS->GetGCPProjection();
AfxMessageBox(poVDS->GetGCPProjection());
AfxMessageBox(poVDS->GetGCPCount());
AfxMessageBox(m_GcpProjection);
int i = 0;
int nBandCount = poDataset->GetRasterCount();
for( i = 0; i < nBandCount; i++ )
{
VRTSourcedRasterBand *poVRTBand;
GDALRasterBand *poSrcBand;
GDALDataType eBandType;
poSrcBand = poDataset->GetRasterBand(i+1);
eBandType = poSrcBand->GetRasterDataType();
// Create this band.
poVDS->AddBand( eBandType, NULL );
poVRTBand = (VRTSourcedRasterBand *)
poVDS->GetRasterBand( i+1 );
poVRTBand->AddSimpleSource( poSrcBand, 0, 0,nRasterXSize,
nRasterYSize,0,0,nRasterXSize, nRasterYSize);
/*
--------------------------------------------------------------------
*/
/* copy over some other information of
interest. */
/*
--------------------------------------------------------------------
*/
poVRTBand->CopyCommonInfoFrom( poSrcBand );
}//for
// Write to the output file using CopyCreate().
GDALDriver * poDriver = NULL;
poDriver = poDataset->GetDriver();
GDALDataset* pNewDs = poDriver->CreateCopy( (LPCTSTR)strSaveName,
poVDS, FALSE, NULL,NULL,NULL);
//GDALDataset* mNewDs = RasterIO()
if( pNewDs!= NULL)
delete pNewDs;
delete poVDS;
CPLFree( pszOutputSRS );
}//if !poDataset == NULL
if(poDataset != NULL)
delete poDataset;
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
