一、椭球
1、常用椭球参数
克拉索夫斯基椭球 | 1975年国际椭球 | WGS-84椭球 | CGCS2000(2000中国大地坐标系) | |
---|---|---|---|---|
a | 6378245(m) | 6378140(m) | 6378137(m) | 6378137(m) |
b | 6356863.0187730473(m) | 6356755.2881575287(m) | 6356752.3142 (m) | 6356752.3141(m) |
c | 6399698.9017827110(m) | 6399596.6519880105(m) | 6399593.6258(m) | 6399593.6259(m) |
α | 1/298.3 | 1/298.257 | 1/298.257223563 | 1/298.257222101 |
e2 | 0.006693421622966 | 0.006694384999588 | 0.00669437999013 | 0.00669438002290 |
e’2 | 0.006738525414683 | 0.006739501819473 | 0.00673949674227 | 0.00673949677548 |
决定旋转椭球的大小和形状只要知道以上参数中的两个就可(至少包含一个长度元素)
我国的北京54坐标系是基于克拉索夫斯基椭球建立的,
西安80坐标系是基于1975年国际椭球建立,
GPS定位系统是基于WGS-84椭球,
CGCS2000坐标系是全球地心坐标系在我国的体现,坐标系原点是地球质心(包括大气和海洋),
Z轴指向国际时间局1984年定义的CTP,X轴为起始子午面与通过原点且正交于Z轴的赤道面,构成右手地心地固直角坐标系
2、
C#构建椭球类
class Ellipsoid
{
//长半轴
public double a;
//短半轴
public double b;
//扁率
public double f;
//第一偏心率平方
public double ec;
//第二偏心率平方
public double ecc;
/// <summary>
///
/// </summary>
/// <param name="a为长半轴"></param>
/// <param name="Invf为扁率"></param>
public Ellipsoid(double a,double Invf)
{
this.a = a;
f = 1 / Invf;
evaluate();
}
public Ellipsoid()
{
a = 6378137.000;
f = 1 / 298.3;
evaluate();
}
void evaluate()
{
double b = a - a * f;
ec = 1 - Math.Pow(b, 2) / Math.Pow(a, 2);
ecc= ec / (1-ec);
}
}
二、大地坐标与空间直角坐标转换
1、转换类
class Transform
{
public double a;
public double ec;
public double ecc;
public Transform(Ellipsoid ellipsoid)
{
a = ellipsoid.a;
ec = ellipsoid.ec;
ecc = ellipsoid.ecc;
}
public PointXYZ BLHToXYZ(PointBLH bLH)
{
double B = bLH.B;
double L = bLH.L;
double H = bLH.H;
//辅助计算公式
double W = Math.Sqrt(1 - ec * Math.Pow(Math.Sin(B), 2));
double N = a / W;
PointXYZ xYZ = new PointXYZ();
xYZ.X = (N + H) * Math.Cos(B) * Math.Cos(L);
xYZ.Y = (N + H) * Math.Cos(B) * Math.Sin(L);
xYZ.Z = (N * (1 - ec) + H) * Math.Sin(B);
return xYZ;
}
public PointBLH XYZToBLH(PointXYZ xYZ)
{
double X = xYZ.X;
double Y = xYZ.Y;
double Z = xYZ.Z;
PointBLH bLH = new PointBLH();
bLH.L = Math.Atan(Y / X);
if (X < 0)
bLH.L += Math.PI;
var r = Math.Sqrt(X * X + Y * Y);
var B1 = Math.Atan(Z / r);
double B2;
while (true)
{
var W1 = Math.Sqrt(1 - ec * (Math.Sin(B1) * Math.Sin(B1)));
var N1 = a / W1;
B2 = Math.Atan((Z + N1 * ec * Math.Sin(B1)) / r);
if (Math.Abs(B2 - B1) <= 0.0000000001)
break;
B1 = B2;
}
bLH.B = B2;
var W = Math.Sqrt(1 - ec * (Math.Sin(B2) * Math.Sin(B2)));
var N = a / W;
bLH.H = r / Math.Cos(B2) - N;
return bLH;
}
}
把输入的点构造成Mypoint类
class PointBLH
{
public double B;
public double L;
public double H;
}
class PointXYZ
{
public double X;
public double Y;
public double Z;
}
2、大地坐标系与空间之间坐标系转换公式
大地坐标换算到空间直角坐标
P点投影到椭球面上,大地高为H
空间直角坐标换算到大地坐标
解算大地纬度的时候需要进行迭代计算,
保证前后两次迭代<e-10时,可保证计算精度达到0.0001″