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