使用GDAL對遙感影像進行投影轉換
1.在進行投影轉換前需要做的工作：
首先編譯GDAL和proj4庫,GDAL中進行投影轉換依賴的開源庫是proj4，可以到網上下載已經編譯好的包含proj4的GDAL，CSDN上就有很多資源。接下來將編譯好的GDAL和PROJ4的lib和include路徑分別加在VS2010的“附加庫目錄”，和“附加包含目錄”中；在聯結器->附加依賴項中新增編譯好的proj.lib檔案和gdal_i.lib。然後將proj.dll拷到相應的debug檔案下。
環境設定好之後，接下來進行投影轉換的第一步。
2.投影轉換：
我做的是同一橢球體基準面下的轉換，將投影座標轉化為地理座標。只要把原理和環境弄清楚了，相信其他的轉換也很簡單。
OGRSpatialReference Raster_spf;
OGRSpatialReference *Dest_Raster_spf;
Raster_spf=m_pDataset->GetProjectionRef();
if (Raster_spf.IsProjected())
{
char*pszSrcWKT=NULL;
Raster_spf.exportToWkt(&pszSrcWKT);
Dest_Raster_spf=Raster_spf.CloneGeogCS();
char *pszDestWKT=NULL;
Dest_Raster_spf->exportToPrettyWkt(&pszDestWKT);
OGRCoordinateTransformation *poCT;
poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);//反過來就是地理座標轉投影座標
if (NULL==poCT)
{
QMessageBox::warning(this, tr("Failed"), tr("poCT is null!"), tr("OK"));
}
int ndbFlag=poCT->Transform(4,dbX,dbY,NULL);
if (ndbFlag)
{
OGRCoordinateTransformation::DestroyCT(poCT);
}
其中dbX和dbY儲存的是影像四個頂點座標。
3,轉換後處理
首先，重新計算地理六引數和行列號
dbRes = dbRes / 111000; //投影座標轉地理座標解析度需除以111000m,反過來就得乘以111000m.原理參考上述部落格
double dbMinx = 0;
double dbMaxx = 0;
double dbMiny = 0;
double dbMaxy = 0;
dbMinx = min(min(min(dbX[0],dbX[1]),dbX[2]),dbX[3]);
dbMaxx = max(max(max(dbX[0],dbX[1]),dbX[2]),dbX[3]);
dbMiny = min(min(min(dbY[0],dbY[1]),dbY[2]),dbY[3]);
dbMaxy = max(max(max(dbY[0],dbY[1]),dbY[2]),dbY[3]);
//估算行列號
padfTransform[0] = dbMinx; //左上角點座標
padfTransform[3] = dbMaxy;
//padfTransform[1] 畫素寬度, padfTransform[5]畫素高度
padfTransform[1] = dbRes;
padfTransform[5] = -dbRes;
//估算行列數
int nPixels = ceil(fabs(dbMaxx-dbMinx)/dbRes);
int nLines = ceil(fabs(dbMaxy-dbMiny)/dbRes);
接下來，利用gdalwarp來重取樣並將其寫入新的影像中。
GDALDriver* pGDalDriver = NULL;//用於creat的驅動
GDALDataset * poDataset; //結果影像
GDALResampleAlg eResampleMethod=GRA_Bilinear; //預設設定取樣方式為雙線性插值，時間會稍短一些
int iDstWidth=nPixels;
int iDstHeight=nLines;
const char *pszFormat;
pszFormat= m_pDataset->GetDriver()->GetDescription();
pGDalDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
QString FilenameD;
FilenameD=QFileDialog::getSaveFileName(this,tr("Save Image"),"",tr("Images (*.tif)"));
if(FilenameD.isEmpty())
{ return; }
poDataset=pGDalDriver->Create(FilenameD.toStdString().c_str(),iDstWidth,iDstHeight,dataBands,dataType,NULL);
poDataset->SetGeoTransform(padfTransform);
poDataset->SetProjection(pszDestWKT);
//構造座標轉換關係
void *hTransformArg=NULL;
hTransformArg=GDALCreateGenImgProjTransformer((GDALDatasetH)m_pDataset,pszSrcWKT,(GDALDatasetH)poDataset,
pszDestWKT,FALSE,0.0,1);
GDALTransformerFunc pfnTransformer=GDALGenImgProjTransform;
//構造gdalwarp的變化選項
GDALWarpOptions *psWO=GDALCreateWarpOptions();
psWO->papszWarpOptions=CSLDuplicate(NULL);
psWO->eWorkingDataType=dataType;
psWO->eResampleAlg=eResampleMethod;
psWO->hSrcDS=(GDALDatasetH)m_pDataset;
psWO->hDstDS=(GDALDatasetH)poDataset;
psWO->pfnTransformer=pfnTransformer;
psWO->pTransformerArg=hTransformArg;
psWO->nBandCount=dataBands;
psWO->panSrcBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));
psWO->panDstBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));
psWO->panSrcBands[0]=1;
psWO->panDstBands[0]=1;
//建立GDALWarp執行物件並使用GDALWarpOptions來進行初始化
GDALWarpOperation oWo;
oWo.Initialize(psWO);
//執行處理,轉換會耗點時間，請耐心等待
oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);
//釋放和關閉檔案
GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
GDALDestroyWarpOptions(psWO);
GDALClose((GDALDatasetH)m_pDataset);
GDALClose((GDALDatasetH)poDataset);
到此，投影轉換的整個過程就結束了。之前遇到的問題是在轉換過程中執行poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);得到的poCT為NULL，這肯定就是跟proj4庫有關了，後來才發現是環境設定的問題。還有在除錯時，執行到oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);就不動了，原以為是轉換除了問題，但跟蹤發現結果影像中相關的資訊都得到了，應該不會有錯。後來發現是重取樣過程比較慢，需要時間，但最鄰近和雙線性插值可能需要的時間稍短些。