遙感影像中投影坐标系轉換

概述

投影轉換一般是講A投影轉換為B投影，包括地理坐标系和投影坐标系之間的互相轉換等，這裡要提到一點的是不同地理坐标系之間的轉換，不同地理坐标系采用的橢球體和基準面有差異，是以，這之間的轉換還涉及到橢球體之間的轉換，這時候需要考慮到橢球體轉換的問題，一般采用7參數轉換模型或者3參數模型，這種轉換是一種嚴格的投影轉換，轉換之後資料由于橢球體和基準面的不同，導緻坐标會有位置偏差，一般偏差幾十米。有時候為了資料能夠完好疊加，采用不嚴格轉換方法，将待轉換資料的投影資訊中的坐标系參數進行修改，比如UTM28N轉換為WGS84坐标系，轉換之前，将UTM28N坐标系中的基準面通過定義投影修改為WGS84所采用的的基準面，然後進行轉換，這樣的轉換之後前後資料可以完好疊加。

投影轉換條件

• 存在投影坐标系的遙感影像資料；
• GDAL和PROJ三方庫。

資料

輸入資料

UTM等投影坐标系的遙感影像。

輸出資料

WGS84投影坐标系的遙感影像。

投影轉換流程圖

關鍵點分析

在整個流程圖轉換過程中，需要注意以下幾點：

判斷是否存在投影資訊和投影坐标

從GIS資料中擷取投影資訊，通過GDAL中的GetProjectionRef()函數即可擷取；從GIS資料中擷取坐标資訊時就需要注意，有些GIS資料的坐标資訊存儲在影像資料内，其坐标資訊格式為：

1. 左上角X坐标；
2. 地圖單元中的一個象素在X方向上的X分辨率尺度；
3. 旋轉量, 0表示沒有旋轉；
4. 左上角Y坐标；
5. 旋轉量, 0表示沒有旋轉；
6. 地圖單元中的一個象素在Y方向上的Y分辨率尺度的負值。

有些坐标資訊是存儲在世界坐标檔案中，其坐标資訊格式為：

1. 地圖單元中的一個象素在X方向上的X分辨率尺度；
2. 旋轉量, 0表示沒有旋轉；
3. 旋轉量, 0表示沒有旋轉；
4. 地圖單元中的一個象素在Y方向上的Y分辨率尺度的負值；
5. 左上角X坐标；
6. 左上角Y坐标。

// 存儲地理位置的參數資訊
	//  dWorldFileParameters[0] /* 左上角 x */
	//  dWorldFileParameters[1] /* 東西方向一個像素對應的距離 */
	//  dWorldFileParameters[2] /* 旋轉, 0表示上面為北方 */
	//  dWorldFileParameters[3] /* 左上角 y */
	//  dWorldFileParameters[4] /* 旋轉, 0表示上面為北方 */
	//  dWorldFileParameters[5] /* 南北方向一個像素對應的距離 */
        //存儲位置調整，以便和GDAL的參數要求一樣
	dWorldFileParameters[0] = aTempWorldFileParameters[4];
	dWorldFileParameters[1] = aTempWorldFileParameters[0];
	dWorldFileParameters[2] = aTempWorldFileParameters[2];
	dWorldFileParameters[3] = aTempWorldFileParameters[5];
	dWorldFileParameters[4] = aTempWorldFileParameters[1];
	dWorldFileParameters[5] = aTempWorldFileParameters[3];
           

兩者的存儲坐标資訊内容一樣，但是兩者存儲的坐标資訊順序不一樣，為避免應用過程中出現錯誤則需要人工對讀取後的坐标資訊順序進行調整。

擷取投影名稱

擷取GIS投影資訊後，發現無法進行各個投影資訊的比較，比如：

原始影像投影坐标系：PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS  84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32628"]]。

轉換後投影坐标系：GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]。

OGRSpatialReference fRef;
	char *tmp = NULL;
	tmp = (char *)malloc(strlen(pProj) + 1);
	strcpy_s(tmp, strlen(pProj)+1, pProj);
	fRef.importFromWkt(&tmp);
	OGR_SRSNode* pOgrNode = fRef.GetRoot();
	OGR_SRSNode* pOgrNodeChild = pOgrNode->GetChild(0);
	if (pOgrNodeChild == NULL)
	{
		return false;
	}
	const char* pProjName = pOgrNodeChild->GetValue();
           

這時就得建立OGRSpatialReference類中OGR_SRSNode節點類，再通過GetValue函數擷取投影坐标系名稱，比如擷取示例原始影像投影坐标系名稱為WGS 84 / UTM zone 28N，擷取轉換後的投影坐标系名稱為WGS 84，然後對坐标系進行比較就可以滿足要求。

進行投影坐标轉換

OGRSpatialReference fRef, tRef;
	char *tmp = NULL;
	/** 獲得projRef的一份拷貝 **/
	/** 由于projRef是const char*,下面的一個函數不接受，是以需要轉換成非const **/
	tmp = (char *)malloc(strlen(projRef) + 1);
	strcpy_s(tmp, strlen(projRef)+1, projRef);

	/** 設定原始的坐标參數，和test.tif一緻 **/
	fRef.importFromWkt(&tmp);
	/** 設定轉換後的坐标 **/
	tRef.SetWellKnownGeogCS(m_strProjName.c_str());
	
	/** 進行坐标轉換，下面的内容如果不安裝proj将會無法編譯 **/
	OGRCoordinateTransformation *coordTrans;
	coordTrans = OGRCreateCoordinateTransformation(&fRef, &tRef);
	
	coordTrans->Transform(1, &DPt.m_dX, &DPt.m_dY);
           

由OGRSpatialReference類中的importFromWkt函數和SetWellKnownGeogCS函數建立對象，再根據OGRCoordinateTransformation類中Transform函數進行轉換。由OGRCreateCoordinateTransformation建立OGRCoordinateTransformation類對象時需判斷OGRSpatialReference類對象是否建立成功，然後再導入PROJ4.0庫，如果不導入proj4.0将無法建立OGRCoordinateTransformation類對象，進而無法進行坐标轉換。

結論

在C++程式設計中，通過轉換投影可以将各種投影坐标系資料顯示在同一坐标系中，極大友善GIS資料顯示。