地心空间（直角）坐标系－－定义为原点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>-10-14T03:16:27.000000</TIME></Point><Point><Location><X>-2.0531141383e+06</X><Y>4.445277e+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>-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、/view/1114841.htm
2、/view/1115634.htm
3、/view/284430.htm
4、/zh-cn/help/main/10.1/index.html#/na/003r0000002s000000/
