https://blog.csdn.net/zhouxuguang236/article/details/17468171
看到一篇講的很詳細的文章，和大家分享一下，感謝原創部落客。
一、影像投影轉換的概念
影像投影轉換就是将一個地理坐标系統轉換到另一個坐标系統，如果在同一個橢球基準面下的轉換就是嚴密的轉換，如果在同一個橢球體不同基準面的轉換是不嚴密的，不同橢球體之間的轉換是不嚴密的，這就需要用到七參數、三參數等方法。需要兩個不同坐标系統下公共點坐标求得系數。例如北京54和WG4-84坐标下的同一點的經緯度或者是經過投影後的平面坐标也是不同的。那麼影像投影主要分為哪些步驟呢？說白了，就三個步驟，第一，坐标轉換；第二，影像的重采樣，最後就是寫入到新檔案中。
首先來說第一步，坐标轉換需要轉換四個坐标，也就是四個角點。也許有人說兩個點就夠了，左下點和右上點。在此，我告訴你，這是錯誤的。投影轉換後這個四個角點組成的矩形那麼很有可能就不是矩形了，如果你取兩個點做轉換那麼後面的影像投影轉換後的範圍就不正确了。
或者再有人問，我怎麼知道影像的四個角點的坐标啊？這個很簡單，通過仿射變換系數，它就是影像的像素坐标（行列号）和地理坐标之間進行關聯的系數。一般是六個參數。在GDAL中可以通過以下這個函數來獲得，如果影像有仿射變換系數的話。如果沒有仿射變換系數但是有控制點的話也能解算出系數；如果都沒有，說明這幅影像是沒有地理參考的，那麼狠遺憾的告訴你，這個影像就不能做投影轉換。還有就是，如果你這幅影像沒有投影的話也就不能做投影轉換了，因為你根本不知道這幅影像的投影是啥。
| CPLErr GDALDataset::GetGeoTransform | ( | double * | padfTransform | ) |
其中padfTransform就存儲了這六個參數。這是個六個double型的數的數組。
在向北向的圖像中，padfTransform[1]代表像素寬度，padfTransform[5]代表像素高度。圖像左上角的坐标是（padfTransform[0],padfTransform[3]），adfGeoTransform[1] X方向也就是橫向的分辨率大小，padfTransform [2] 旋轉系數，如果為0，就是标準的正北向圖像，padfTransform [4] 旋轉系數，如果為0，就是标準的正北向圖像，知道了這兩個參數的意義，那麼我們就可以得到四個角點的地理坐标了。
在正北向的圖像中，四個角點的坐标計算如下：
//計算源圖像的MBR
double dbX[4];
double dbY[4];
double dbZ[4] = {0,0,0,0};
dbX[0] = adfDstGeoTransform[0]; //左上角點坐标
dbY[0] = adfDstGeoTransform[3];
//右上角坐标
dbX[1] = adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1];
dbY[1] = adfDstGeoTransform[3];
//右下角點坐标
dbX[2]= adfDstGeoTransform[0]+ nXsize*adfDstGeoTransform[1]+
nYsize*adfDstGeoTransform[2];
dbY[2] = adfDstGeoTransform[3] + nXsize*adfDstGeoTransform[4]
+ nYsize*adfDstGeoTransform[5];
//左下角坐标
dbX[3] = adfDstGeoTransform[0];
dbY[3] = adfDstGeoTransform[3] + nYsize*adfDstGeoTransform[5];
這樣我們就找到了需要參與投影轉換的坐标點了，下一步就是坐标轉換，坐标轉換的過程通過GDAL的接口實作，其底層依賴了PROJ4地圖投影開源類庫。對于不同的橢球體之間變換需要用到三參數布爾莎或者七參數布爾莎模型，具體過程就是首先将經緯度大地坐标轉換為地心坐标系下的空間直角坐标，然後用布爾莎模型計算，最後将計算後的結果重新轉換到目标地理坐标系統下的經緯度大地坐标。有了需要轉換的坐标後，我們将對上述四個點的坐标進行變換，其函數如下：
bool TranformCoordsOGR(char* pszSrcWkt,char* pszDestWkt, int nCount,double* x,double* y,double* z ,double *dfParaSrc,double* dfParaDst,int nParaCount) { //建立OGR的空間參考系 OGRSpatialReference oSrcSrs; //源坐标系統 OGRSpatialReference oDestSrs; //目的坐标系統 oSrcSrs.importFromWkt(&pszSrcWkt); oDestSrs.importFromWkt(&pszDestWkt); int nSameGeoCS = oSrcSrs.IsSameGeogCS(&oDestSrs); //相同的橢球基準面，則進行轉換 if (nSameGeoCS) { OGRCoordinateTransformation *poCT = NULL; poCT = OGRCreateCoordinateTransformation( &oSrcSrs,&oDestSrs ); if (NULL == poCT) { return false; } int nFlag = poCT->Transform(nCount,x,y,z); if (nFlag) { OGRCoordinateTransformation::DestroyCT(poCT); return true; } return false; } else //不同的橢球體基準面,要設定七參數或者三參數 { int nFlag = 0; //如果是地理坐标系，直接轉換為空間直角坐标 OGRErr err = 0; double dbAsrc = 0; double dbBsrc = 0; double dbEsrc = 0; dbAsrc = oSrcSrs.GetSemiMajor(&err); dbBsrc = oSrcSrs.GetSemiMinor(&err); dbEsrc = 1-pow((dbBsrc/dbAsrc),2.0); if (oSrcSrs.IsProjected()) { OGRSpatialReference* poTmpSRS = oSrcSrs.CloneGeogCS(); OGRCoordinateTransformation *poCTTmp = NULL; poCTTmp = OGRCreateCoordinateTransformation( &oSrcSrs,poTmpSRS ); nFlag = poCTTmp->Transform(nCount,x,y,z); if (!nFlag) { OGRCoordinateTransformation::DestroyCT(poCTTmp); return false; } OGRCoordinateTransformation::DestroyCT(poCTTmp); //pj_geodetic_to_geocentric(dbAsrc,dbEsrc,nCount,0,x,y,z); } //将經緯度坐标轉換為空間直角坐标 CGeoEllipse geoEllipse(dbAsrc,dbBsrc); double dbX = 0; double dbY = 0; double dbZ = 0; for (int i = 0; i < nCount; i ++) { geoEllipse.BLH_XYZ(x[i],y[i],z[i],dbX,dbY,dbZ); x[i] = dbX; y[i] = dbY; z[i] = dbZ; } //七參數模型 vector<double> vecX; vector<double> vecY; vector<double> vecZ; for (int i = 0; i < nCount; i ++) { double dbX = dfParaSrc[0] + (1+dfParaSrc[6])*(x[i]+dfParaSrc[5]*y[i]-dfParaSrc[4]*z[i]); vecX.push_back(dbX); double dbY = dfParaSrc[1] + (1+dfParaSrc[6])*(-dfParaSrc[5]*x[i]+y[i]+dfParaSrc[3]*z[i]); vecY.push_back(dbY); double dbZ = dfParaSrc[2] + (1+dfParaSrc[6])*(dfParaSrc[4]*x[i]-dfParaSrc[3]*y[i]+z[i]); vecZ.push_back(dbZ); } memcpy(x,&vecX[0],sizeof(double)*nCount); memcpy(y,&vecY[0],sizeof(double)*nCount); memcpy(z,&vecZ[0],sizeof(double)*nCount); double dbAdst = 0; double dbBdst = 0; double dbEdst = 0; dbAdst = oDestSrs.GetSemiMajor(&err); dbBdst = oDestSrs.GetSemiMinor(&err); //再将空間直角坐标轉換為地理坐标，即經緯度坐标 CGeoEllipse geoEllipse1(dbAdst,dbBdst); for (int i = 0; i < nCount; i ++) { geoEllipse1.XYZ_BLH(x[i],y[i],z[i],dbX,dbY,dbZ); x[i] = dbX; y[i] = dbY; z[i] = dbZ; } if (oDestSrs.IsProjected()) { const char* pszProjName = oDestSrs.GetAttrValue("PROJECTION"); OGRSpatialReference* poTmpSRS = oDestSrs.CloneGeogCS(); int nZone = oDestSrs.GetUTMZone(); char* pszTmp; poTmpSRS->exportToWkt(&pszTmp); OGRCoordinateTransformation *poCTTmp = NULL; poCTTmp = OGRCreateCoordinateTransformation( poTmpSRS,&oDestSrs ); if (NULL == poCTTmp) { //MessageBox(NULL,_T("失敗"),_T("提示"),MB_OK); return false; } nFlag = poCTTmp->Transform(nCount,x,y,z); if (!nFlag) { OGRCoordinateTransformation::DestroyCT(poCTTmp); return false; } OGRCoordinateTransformation::DestroyCT(poCTTmp); return true; } return true; } return false; }
上述代碼中有空間直角坐标和大地坐标之間的變換，這個是我自己寫的，讀者也可以使用PROJ中的接口進行變換。
第二步就是影像重采樣了，重采樣就是通過原始影像的像素值内插得到新到采樣點上的像素值。這個可以直接用GDAL中重采樣接口來完成。
第三步就不用詳細說了，一般投影轉換後需要将投影後的影像寫入的新檔案中，直接用GDAL的讀寫接口來完成。
二、GDAL影像投影轉換的過程中分辨率和仿射變換系數的估算
上一節已經完成了點的投影轉換，那麼我們現在就要估算投影後的像素分辨率大小和仿射變換系數了。
如果投影變換前是投影坐标系統，投影轉換後也是投影坐标系統，或者說另外一種情況：投影變換前是地理坐标系統，投影變換後也是地理坐标系統，并且坐标的機關都一緻的，那麼分辨率大小基本上沒變換，可以用變換前的分辨率大小。如果變換前是地理坐标系統，投影變換後是投影坐标系統，假設地理坐标系統以度為機關，投影坐标系統以米為機關，那麼投影後的像素大小可以這樣估計，因為經線上一個緯度的距離大約是111km，那麼變換後的分辨率可以由原始分辨率乘以111000；相反的話，如果變換前是投影坐标系統，投影變換後是地理坐标系統，假設地理坐标系統以度為機關，投影坐标系統以米為機關，同理投影後的像素大小可以這樣估計，變換後的分辨率可以由原始分辨率除以111000。
仿射變換系數這樣也就可以确定了，左上角的坐标就是最小x值，最大y值，在變換後的四個角點坐标中比較獲得，分辨率上一段也講了如何獲得。對于正北向的圖像這就夠了。 然後行列數就用變換後的四個角點組成的區域的MBR的寬度除以橫向分辨率得到列數，高度除以縱向分辨率得到行數。
具體的代碼片段如下：
//轉換為PROJ4結構
projPJ pj_SourceProjection = NULL, pj_DestinationProjection = NULL;
pj_SourceProjection = pj_init_plus(pszSrcProj);
pj_DestinationProjection = pj_init_plus(pszDestProj);
if (pj_is_latlong( pj_SourceProjection ) && !pj_is_latlong(pj_DestinationProjection))
{
dbRes = dbRes * 111000;
}
else if (!pj_is_latlong( pj_SourceProjection ) &&
pj_is_latlong(pj_DestinationProjection))
{
dbRes = dbRes / 111000;
}
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]);
//估算行列号
adfDstGeoTransform[0] = dbMinx; //左上角點坐标
adfDstGeoTransform[3] = dbMaxy;
//padfTransform[1] 像素寬度, padfTransform[5]像素高度
adfDstGeoTransform[1] = dbRes;
adfDstGeoTransform[5] = -dbRes;
//估算行列數
nPixels = ceil(fabs(dbMaxx-dbMinx)/dbRes);
nLines = ceil(fabs(dbMaxy-dbMiny)/dbRes);
三、投影後處理
投影後處理主要就是重采樣和寫入資料到新檔案中了。
重采樣不多說，主要有最鄰近、雙線性插值、立方卷積法等。可以直接調用GDAL接口，寫入也就不多說了。
這個在網上有很多這樣的代碼，包括李民錄大哥的部落格等。不過還是将代碼共享出來
// 建立輸出檔案
hDstDS = GDALCreate(hDriver,info.m_strOutputFile.c_str(), nPixels, nLines,GDALGetRasterCount(hSrcDS),eDT,NULL);
CPLAssert( hDstDS != NULL );
// 寫入投影
GDALSetProjection( hDstDS,info.m_strOutWkt.c_str());
GDALSetGeoTransform( hDstDS,adfDstGeoTransform );
// 複制顔色表，如果有的話
GDALColorTableH hCT;
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1));
if( hCT != NULL )
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1),hCT );
ProgressCreateEx( &hPolygonizeProgress,_T(""),TRUE,FALSE );
ProgressBeginEx( hPolygonizeProgress,PROGRESS_MODE_PERCENT, 100);
ProgressSetStepTitleEx( hPolygonizeProgress, _T("圖像重投影"));
// 建立變換選項
GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS =hSrcDS;
psWarpOptions->hDstDS =hDstDS;
int nBandCount = GDALGetRasterCount(hSrcDS);
psWarpOptions->nBandCount = nBandCount;
psWarpOptions->panSrcBands =
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
for (int i = 0; i < nBandCount; i ++)
{
psWarpOptions->panSrcBands[i] = i+1;
}
psWarpOptions->panDstBands =
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
for (int i = 0; i < nBandCount; i ++)
{
psWarpOptions->panDstBands[i] = i+1;
}
psWarpOptions->pfnProgress = ImageReProjectProgress;
psWarpOptions->eResampleAlg = info.m_enResampleAlg;
// 建立重投影變換函數
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( hSrcDS,
GDALGetProjectionRef(hSrcDS),
hDstDS,
GDALGetProjectionRef(hDstDS),
FALSE,0.0, 1 );
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// 初始化并且執行變換操作
GDALWarpOperation oOperation;
oOperation.Initialize(psWarpOptions );
oOperation.ChunkAndWarpImage(0,0,GDALGetRasterXSize(hDstDS),GDALGetRasterYSize(hDstDS));
ProgressEndEx( hPolygonizeProgress );
ProgressCloseEx( &hPolygonizeProgress );
hPolygonizeProgress=NULL;
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
![如何将高程資料(等高線)轉換成xi'an80或者beijing54坐标系[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)