大地坐标表示顺序为：纬度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 mathdef 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;}
能极大的简化计算。
