GDAL座標轉換

1、引言

最近研究了一下GIS、測繪學的座標轉換的問題，感受大部分資料專業性太強，上來就是一通專業性論述；但感受對於相關從業者來講，其實沒必要了解那麼多背景知識的；就經過GDAL這個工具，來簡單總結下座標轉換相關的問題。
GDAL座標轉換其實也是調用proj4來實現，可是proj4有個特別麻煩的地方，就是座標系描述的部分特別繁複，須要對專業知識有必定的瞭解。使用GDAL則相對簡單不少。

2、地理座標系

地理座標系就是常說的經緯度座標系，好比用GPS直接獲取的座標就是在地理座標系下獲取的。一個真實座標不管怎麼變換，必定會有地理座標系做爲基準，也必定有能夠轉換出來的經緯度座標。以國內的狀況來講，經常使用到的地理座標系有WGS84，beijing54，xian80，CGCS2000這四種。
GDAL能夠像proj4那樣自定義座標系，也能夠僅經過字符串定義一些經常使用的座標系，但本文認爲最方便的仍是經過EPSG數據庫定義的代碼來定義一個地理座標系統；畢竟不少時候須要兼容的地理座標系不少，所有一個個自定義座標系太麻煩。

OGRSpatialReference spatialReference;
//spatialReference.importFromEPSG(4326);//WGS84
//spatialReference.importFromEPSG(4214);//BeiJing54
spatialReference.importFromEPSG(4610);//XIAN80
//spatialReference.importFromEPSG(4490);//CGCS2000

char *pszWKT = nullptr;
spatialReference.exportToPrettyWkt(&pszWKT);
cout << pszWKT << endl;
CPLFree(pszWKT);
pszWKT = nullptr;

如上演示了GDAL定義地理座標系並輸出座標系的信息。四種不一樣的座標系只須要改變相應的代碼編號，用起來十分方便。最終打印輸出了的xian80座標系信息：

在這裏必定要注意，使用這種方式定義地理座標系必定要經過配置GDAL_DATA路徑，不然控制檯會輸出錯誤信息：

CPLSetConfigOption("GDAL_DATA", "gdaldata");

「gdaldata」表示一個路徑（這裏用的是相對路徑，固然也能夠設置成絕對路徑），是GDAL編譯完成後會生成的一個目錄，裏面記錄了各類座標系的參數文件。例如進入這個文件夾，用Excel打開pcs.csv這個文件，就能夠看到各類座標系的參數了。

除了這種方法，也能夠在環境變量中設置GDAL_DATA變量來實現。

3、投影平面座標系

經緯度座標是曲面上的座標，曲面上的座標投影到平面，不一樣的投影方式就會產生不一樣的平面座標；即便是同一種投影方式，不一樣的參數獲得的平面座標也會不一樣。也就是說，地理座標系下的經緯度座標與投影座標系下的平面座標，是一對多的關係而不是一對一的關係。以國內的狀況來講，常常用到的投影有橫軸墨卡託投影，高斯-克呂格投影和UTM投影。
在Global Mapper中三種投影設置方式以下：



能夠看出，高斯-克呂格投影和UTM投影其實都是橫軸墨卡託投影，二者都是經過設置帶號（Zone）來實現設置橫軸墨卡託投影的具體參數（Parameters）。在GDAL裏面，高斯-克呂格投影就是經過設置橫軸墨卡託投影來實現的。以下演示了一個xian80座標系，3度帶帶號38的橫軸墨卡託投影。

OGRSpatialReference spatialReference;
spatialReference.importFromEPSG(4610);          //XIAN80
spatialReference.SetTM(0, 114, 1.0, 38500000, 0);

設置橫軸墨卡託投影的函數SetTM()有六個參數：

參數	設置
dfCenterLat	通常爲0
dfCenterLong	中央經線，決定了是哪一投影帶
dfScale	通常爲1.0，UTM設置爲0.9996
dfFalseEasting	通常爲500000，高斯-克呂格前面加上帶號
dfFalseNorthing	通常爲0

用以前方法獲得座標系信息並輸出，信息以下：

4、座標轉換

定義好座標系以後，就能夠進行座標轉換了。以下爲xian80地理座標系下某點（113.6，38.8）用高斯-克呂格投影到帶平面座標系：

OGRSpatialReference* pLonLat = spatialReference.CloneGeogCS();
OGRCoordinateTransformation* LonLat2XY = OGRCreateCoordinateTransformation(pLonLat, &spatialReference);
if (!LonLat2XY)
{       
    return 1;
}
    
double x = 113.6;
double y = 38.8;
printf("經緯度座標：%.9lf\t%.9lf\n", x, y);
if (!LonLat2XY->Transform(1, &x, &y))
{
    return 1;
}
printf("平面座標：%.9lf\t%.9lf\n", x, y);

OGRCoordinateTransformation::DestroyCT(LonLat2XY);
LonLat2XY = nullptr;

這裏經過複製以前定義的高斯-克呂格投影平面座標系獲得相同的地理座標系（固然也能夠自定義新座標系），而後使用OGRCoordinateTransformation::Transform()方法來進行座標轉換。最後的輸出結果爲：

經過Global Mapper的座標轉換工具來驗證結果是否正確：


比較能夠發現二者轉換的結果基本一致。除此以外，將平面座標逆投影到地理座標也是能夠的，只須要在OGRCreateCoordinateTransformation（）的時候顛倒下順序便可。

5、其餘

1.GDAL默認不編譯proj.4，使用的時候須要添加proj.4的支持。
2.同一地理座標系的投影轉換是嚴密的，但不一樣地理座標系之間須要先轉換到地心立體座標系，而後經過七參數轉換。
3.能夠根據座標值選擇正確的分帶，使用這個分帶的上下幾個分帶進行投影問題也不是很大。可是分帶差距太大可能有問題，由於投影公式只能在必定範圍內有效；即不一樣的軟件對的時候結果比較一致，錯的時候結果可能不一樣。
以上三個問題有時間再作進一步的研究和總結。

6、參考文獻

1.GDAL源碼剖析（十一）之OGR投影說明
2.墨卡託投影、高斯-克呂格投影、UTM投影及我國分帶方法
3.GDAL庫學習筆記(五):座標系之間的轉化
4.GIS座標轉換庫Proj.4的使用
5.GDAL影像投影轉換