WGS84GPS坐标与大地坐标系互相转化

大地坐标表示顺序为：纬度B、经度L、海拔H
空间直角坐标表示方法：X、Y、Z

进行地图投影等价于将BLH转换为XYZ，涉及到两个椭球体以及坐标系之间的转换，本文主要讨论WGS84坐标系下BLH与XYZ相互转换的办法。

WGS84坐标系的参数：
已知长半轴a=6378137 ，
WGS84椭球扁率：f=1/298.257223563
椭球扁率f=（a-b）/a ,可求出b，
椭球第一偏心率：ee=（aa-bb）/aa，带入可求得e
W为第一辅助系数：W=√(1-e²sin²B) B为纬度
N为椭圆曲率半径：N=a/W

用python代码写出：

import math
def get_coordinate(latitude,longitude,altitude):
	B=math.radians(latitude);
	L=math.radians(longitude);
	H=altitude;
	f=1/298.257223563;
	r=6378137;
	b=r*(1-f);
	e=math.sqrt(2*f-f*f);
	N=r/math.sqrt(1-e*e*math.sin(B)*math.sin(B));
	data =[(N+H)*math.cos(B)*math.cos(L),(N+H)*math.cos(B)*math.sin(L),(N*(1-e*e)+H)*math.sin(B)];
	return data;

java语言写成了一个方法：

	public static double[] get_coordinate(double latitude,double longitude,double altitude){
		double B=Math.toRadians(latitude) ,L=Math.toRadians(longitude) ,H=altitude,x,y,z;
		double f=1/298.257223563,r=6378137;
		double b=r*(1-f),e=Math.sqrt(2*f-f*f);
		double N=r/Math.sqrt(1-e*e*Math.sin(B)*Math.sin(B));
		x=(N+H)*Math.cos(B)*Math.cos(L);
		y=(N+H)*Math.cos(B)*Math.sin(L);
		z=(N*(1-e*e)+H)*Math.sin(B);
		double[] data ={x,y,z};
		return data;
	}

这里完成了WGS84坐标系转大地坐标系，但是大地坐标系反解WGS84坐标系较为复杂

我们可以直接通过XY求出L经度，但是对于纬度和海拔求解较为复杂。
这里我们可以通过一种简单的方法将计算最简化。我们将N视为不变，来计算前后GPS坐标。程序如下：
python:

def return_gps(m,n,q):
	f=1/298.257223563;a=6378137;
	b=a*(1-f);e=math.sqrt(2*f-f*f);
	L=math.atan(n/m)+math.pi;
	B=math.atan(math.cos(L)*q/m/(1-N*e*e/(N+40.0038)));
	H=(m/math.cos(L)/math.cos(B)-N+n/math.sin(L)/math.cos(B)-N)/2;
	L1=L/math.pi*180;B1=B/math.pi*180;
	data = [B1,L1,H];
	return data;

Java：

	public static double[] returnlocation(double m,double n,double q){
		double f=1/298.257223563,a=6378137;
		double b=a*(1-f),e=Math.sqrt(2*f-f*f);
		double L=Math.atan(n/m)+Math.PI;
		double B=Math.atan(Math.cos(L)*q/m/(1-N*e*e/(N+40.0038)));
		double H=(m/Math.cos(L)/Math.cos(B)-N+n/Math.sin(L)/Math.cos(B)N)/2;
		double L1=L/Math.PI*180,B1=B/Math.PI*180;
		double[] data = {B1,L1,H};
		return data;
	}

能极大的简化计算。