使用PROJ4庫將地心直角座標（XYZ）轉為地心大地座標（BLH）

地心空間（直角）座標系－－定義為原點O與地球質心重合，Z軸指向地球北極，X軸指向格林尼治子午面與地球赤道的交點，Y軸垂直於XOZ平面構成右手座標系。地心空間直角座標系是座標系的一種，測量學上用於描述任一點的位置。

地心大地座標系－－定義為地球橢球的中心與地球質心（品質中心）重合，橢球的短軸與地球自轉軸重合。地心大地經度L，是過地面點的橢球子午面與格林尼治天文檯子午面的夾角；地心大地緯度B，是過點的橢球法線（與參考橢球面正交的直線）和橢球赤道面的夾角；大地高H，是地面點沿橢球法線到地球橢球面的距離。

如所示，P點的座標如果使用XYZ表示，就是地心直角座標，如果使用BLH表示就是地心大地座標。

地心直角座標系一般用來描述衛星位置較多，比如SPOT5衛星的位置。對於SPOT5的遙感影像，裡面的dim檔案中含有描述衛星位置和速度的項。裡面衛星的位置都是使用地心直角座標系來進行描述，比如下面的DIM檔案片段：

<Point><Location><X>-2.0394400196e+06</X><Y>4.2728461045e+06</Y><Z>5.4215671770e+06</Z></Location><Velocity><X>-5.0095518940e+02</X><Y>5.8130406670e+03</Y><Z>-4.7582155460e+03</Z></Velocity><TIME>2008-10-14T03:16:27.000000</TIME></Point><Point><Location><X>-2.0531141383e+06</X><Y>4.4452004277e+06</Y><Z>5.2762355325e+06</Z></Location><Velocity><X>-4.1067833320e+02</X><Y>5.6762657790e+03</Y><Z>-4.9297861590e+03</Z></Velocity><TIME>2008-10-14T03:16:57.000000</TIME></Point>

從上面的dim檔案片段中可以看出，在某一時刻的衛星位置是使用地心直角座標系表示，大多數時候是需要將上面的座標轉為地心大地座標，也就是經緯度和大地高表示的座標。下面就如何使用PROJ4庫來進行轉換進行說明。座標轉換核心函數如下：

/*** 批量將WGS84地心座標系轉為WGS84經緯度座標* @param pTransformArg轉換參數，設定為NULL，設定這個參數是方便用GDAL的函數指標* @param bDstToSrcTRUE為地心轉經緯度，FALSE為經緯度轉地心* @param nPointCount點個數* @param xX座標序列* @param yY座標序列* @param zZ座標序列* @param panSuccess轉換就誒過標記序列* @return 成功執行傳回值為true，否則傳回值為false*/int GeocentLonLatTransform(void *pTransformArg, int bDstToSrc, int nPointCount,    double *x, double *y, double *z, int *panSuccess){if (panSuccess != NULL)memset(panSuccess, FALSE, nPointCount);// 地心座標系const char* geoccs="+proj=geocent +datum=WGS84";// 經緯度，WGS84基準const char* latlon="+proj=latlong +datum=WGS84";projPJ pjGeoccs, pjLatlon; //初始化當前投影對象if(!(pjGeoccs= pj_init_plus(geoccs)))return FALSE;if(!(pjLatlon= pj_init_plus(latlon)))return FALSE;if (bDstToSrc){int iRev = pj_transform(pjGeoccs, pjLatlon, nPointCount, 1, x, y, z);if (iRev != 0)return FALSE;for(int i=0; i<nPointCount; i++){//弧度轉度x[i]*=RAD_TO_DEG;y[i]*=RAD_TO_DEG;}}else{for(int i=0; i<nPointCount; i++){//度轉弧度x[i]*=DEG_TO_RAD;y[i]*=DEG_TO_RAD;}int iRev = pj_transform(pjLatlon, pjGeoccs, nPointCount, 1, x, y, z);if (iRev != 0)return FALSE;}pj_free(pjGeoccs);pj_free(pjLatlon);return TRUE;}

下面我們再編寫一個函數來調用上面的函數進行測試。測試代碼如下，測試中一共使用了12組點，分別進行正變換和逆變換，從逆變換的結果與原始點對比發現，座標與輸入的一致。

int GeoCent2LLH() {double pGeoccsX[12]={-2.3825143026e+06,-953076.900000,-968629.800000,-984133.100000,-999587.000000,-1014989.400000,-1030337.600000,-1045628.000000,-1060860.500000,-1076032.900000,-1091144.700000,-1106195.200000};double pGeoccsY[12]={4.0316337093e+06,-6542517.500000,-6560998.500000,-6578987.500000,-6596481.500000,-6613479.000000,-6629982.500000,-6645987.000000,-6661486.000000,-6676487.000000,-6690984.500000,-6704978.500000}; double pGeoccsZ[12]={5.4665429711e+06,2453130.200000,2397415.200000,2341526.000000,2285467.000000,2229241.500000,2172853.500000,2116305.200000,2059601.200000,2002746.600000,1945745.600000,1888602.700000}; GeocentLonLatTransform(NULL, TRUE, 12, pGeoccsX, pGeoccsY, pGeoccsZ, NULL);for(int i=0; i<12;i++){cout.setf(ios_base::fixed);//設定cout為定點輸出格式(設定當前流為小數形式輸出)cout<<"經緯度:    "<<pGeoccsX[i]<<"    "<<pGeoccsY[i]<<"    "<<pGeoccsZ[i]<<endl;}cout<<"\n\n";GeocentLonLatTransform(NULL, FALSE, 12, pGeoccsX, pGeoccsY, pGeoccsZ, NULL);for(int i=0; i<12;i++){cout.setf(ios_base::fixed);//設定cout為定點輸出格式(設定當前流為小數形式輸出)cout<<"地心座標:    "<<pGeoccsX[i]<<"    "<<pGeoccsY[i]<<"    "<<pGeoccsZ[i]<<endl;}return 0;}

上面測試代碼執行輸出的結果如。需要說明的是，上面的代碼編譯和執行時需要PROJ4庫的支援。

參考資料：

1、http://baike.baidu.cn/view/1114841.htm

2、http://baike.baidu.cn/view/1115634.htm

3、http://baike.baidu.com.cn/view/284430.htm

4、http://resources.arcgis.com/zh-cn/help/main/10.1/index.html#/na/003r0000002s000000/