GDAL坐標轉換

一、引言

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

二、地理坐標系

地理坐標系就是常說的經緯度坐標系，比如用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變量來實現。

三、投影平面坐標系

經緯度坐標是曲面上的坐標，曲面上的坐標投影到平面，不同的投影方式就會產生不同的平面坐標；即使是同一種投影方式，不同的參數得到的平面坐標也會不同。也就是說，地理坐標系下的經緯度坐標與投影坐標系下的平面坐標，是一對多的關系而不是一對一的關系。以國內的情況來說，經常用到的投影有橫軸墨卡托投影，高斯-克呂格投影和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

用之前方法得到坐標系信息並輸出，信息如下：

四、坐標轉換

定義好坐標系之后，就可以進行坐標轉換了。如下為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（）的時候顛倒下順序即可。

五、其他

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

六、參考文獻

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