地心空间(直角)坐标系－－定义为原点O与地球质心重合，Z轴指向地球北极，X轴指向格林尼治子午面与地球赤道的交点，Y轴垂直于XOZ平面构成右手坐标系。地心空间直角坐标系是坐标系的一种，测量学上用于描述任一点的位置。
地心大地坐标系－－定义为地球椭球的中心与地球质心(质量中心)重合，椭球的短轴与地球自转轴重合。地心大地经度L，是过地面点的椭球子午面与格林尼治天文台子午面的夹角；地心大地纬度B，是过点的椭球法线(与参考椭球面正交的直线)和椭球赤道面的夹角；大地高H，是地面点沿椭球法线到地球椭球面的距离。
如下图所示，P点的坐标如果使用XYZ表示，就是地心直角坐标，如果使用BLH表示就是地心大地坐标。
地心直角坐标系一般用来描述卫星位置较多，比如SPOT5卫星的位置。对于SPOT5的遥感影像，里面的dim文件中含有描述卫星位置和速度的项。里面卫星的位置都是使用地心直角坐标系来进行描述，比如下面的DIM文件片段：
-2.0394400196e+06
4.2728461045e+06
5.4215671770e+06
-5.0095518940e+02
5.8130406670e+03
-4.7582155460e+03
-10-14T03:16:27.000000
-2.0531141383e+06
4.445277e+06
5.2762355325e+06
-4.1067833320e+02
5.6762657790e+03
-4.9297861590e+03
-10-14T03:16:57.000000
从上面的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
{
//弧度转度
x[i]*=RAD_TO_DEG;
y[i]*=RAD_TO_DEG;
}
}
else
{
for(int i=0; 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<
}
cout<
GeocentLonLatTransform(NULL, FALSE, 12, pGeoccsX, pGeoccsY, pGeoccsZ, NULL);
for(int i=0; i<12;i++)
{
cout.setf(ios_base::fixed);//设置cout为定点输出格式(设置当前流为小数形式输出)
cout<
}
return 0;
}上面测试代码执行输出的结果如下图。需要说明的是，上面的代码编译和执行时需要PROJ4库的支持。
参考资料：
