使用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 bDstToSrc		TRUE为地心转经纬度，FALSE为经纬度转地心
* @param nPointCount	点个数
* @param x				X坐标序列
* @param y				Y坐标序列
* @param z				Z坐标序列
* @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/