工程坐标转换方法C#代码实现
1. 前言
在前面的文章中系统的阐述了工程坐标的转换类别和转换的方法。关于转换代码实现,有很多的类库:
这里针对GPS接收的WGS84椭球的经纬度转换为地方坐标系的问题,利用C#,对工程坐标转换方法和步骤做出详细的解答。不基于任何类库和函数库,也未使用矩阵库,可以便利的将代码移植到任何语言。
2. 计算总体框架
根据上一篇文章中对七参数、四参数、高程拟合在坐标转换的作用和使用条件的阐述,我们可以将上一篇文章第7节的总结图,按照计算的流程重新绘制。
根据上图可知,预将WGS84椭球的GPS坐标需要经过5次转换。其中,
- 转换1、转换3在charlee44的博客:大地经纬度坐标与地心地固坐标的转换中详细讲解了,并且有C++代码的实现,利用C#重构即可。
- 转换2、转换5,以及他们的组合,在我的上一篇文章(工程)坐标转换类别和方法也详细的讲解了。
因此,根据计算原理,直接可以利用C#代码实现。
3. C#代码实现
3.1 整体类的构建
5个转换是对点的操作,不妨构建自定义点类MyPoint
,在这个类中定义转换方法。在实现转换方法之前,需要定义数据属性,以承载转换参数和转换数据。代码框架如下:
internal class MyPoint | |
{ | |
// 定义椭球类型。这里仅列举了4中国内常见的椭球类型 | |
// 国际椭球可以增加自行定义 | |
public enum EllipsoidType | |
{ | |
WGS84, | |
CGCS2000, | |
西安80, | |
北京54 | |
} | |
//大地坐标经度、维度、高度 | |
public double L { get; set; } | |
public double B { get; set; } | |
public double H { get; set; } | |
//空间坐标系 | |
public double X { get; set; } | |
public double Y { get; set; } | |
public double Z { get; set; } | |
//七参数转换后的空间坐标 | |
public double X2 { get; set; } | |
public double Y2 { get; set; } | |
public double Z2 { get; set; } | |
private double a = 0, f = 0, b = 0, e = 0, e2 = 0; //椭球参数 | |
private readonly double rho = 180 / Math.PI; | |
private readonly double d2r = Math.PI / 180; | |
public double Xs { get; set; } | |
public double Ys { get; set; } | |
public double Hs { get; set; } | |
//七参数 三个线性平移量-单位米 三个旋转平移量-十进制秒为单位(运算时注意转换为度) 比例因子-单位百万分率 (ppm) | |
//测量队给出的七参数单位与计算的单位不同,要进行单位转化 1 秒=0.0000048481373323 弧度 | |
//尺度因子有两种单位的表示形式,一种结果约为1,如1.0000045,用k表示; | |
//另一种就是ppm的表示形式,稍微比1大一点,如4.5,用m表示。k=m/1000000 | |
private double dx = 0, dy = 0, dz = 0, rx = 0, ry = 0, rz = 0, m = 0, k = 0; | |
} |
3.2 椭球参数赋值
常见的椭球参数值在我的文章经纬度坐标转换为工程坐标可以找到,这里选取与上述代码对应的4类椭球,并在上述MyPoint
类中增加函数EllipsoidParam(EllipsoidType type)
。
/// <summary> | |
/// 椭球参数设置 | |
/// </summary> | |
/// <param name="type">椭球类型</param> | |
private void EllipsoidParam(EllipsoidType type) | |
{ | |
// CGCS2000 椭球参数 | |
if (type == EllipsoidType.CGCS2000) | |
{ | |
this.a = 6378137; | |
this.f = 1 / 298.257222101; | |
} | |
// 西安 80 | |
else if (type == EllipsoidType.西安80) | |
{ | |
this.a = 6378140; | |
this.f = 1 / 298.257; | |
} | |
// 北京 54 | |
else if (type == EllipsoidType.北京54) | |
{ | |
this.a = 6378245; | |
this.f = 1 / 298.3; | |
} | |
// WGS-84 | |
else | |
{ | |
this.a = 6378137; | |
this.f = 1 / 298.257223563; | |
} | |
this.b = this.a * (1 - this.f); | |
this.e = Math.Sqrt(this.a * this.a - this.b * this.b) / this.a; //第一偏心率 | |
this.e2 = Math.Sqrt(this.a * this.a - this.b * this.b) / this.b; //第二偏心率 | |
} |
3.3 转换1、3(大地经纬度坐标与地心地固坐标的转换)
charlee44的博客有C++代码的实现,现在利用C#重构即可。上述MyPoint
类中增加BLH2XYZ(EllipsoidType type)
和XYZ2BLH(EllipsoidType type)
两个函数。
/// <summary> | |
/// 经纬度坐标转空间直角坐标 | |
/// </summary> | |
/// <param name="type">椭球类型</param> | |
public void BLH2XYZ(EllipsoidType type = EllipsoidType.WGS84) | |
{ | |
EllipsoidParam(type); | |
double sB = Math.Sin(this.B * d2r); | |
double cB = Math.Cos(this.B * d2r); | |
double sL = Math.Sin(this.L * d2r); | |
double cL = Math.Cos(this.L * d2r); | |
double N = this.a / (Math.Sqrt(1 - this.e * this.e * sB * sB)); | |
this.X = (N + this.H) * cB * cL; | |
this.Y = (N + this.H) * cB * sL; | |
this.Z = (N * (1 - this.e * this.e) + this.H) * sB; | |
this.X2 = this.X; | |
this.Y2 = this.Y; | |
this.Z2 = this.Z; | |
} | |
/// <summary> | |
/// 空间直角坐标转经纬度坐标 | |
/// </summary> | |
/// <param name="type">椭球类型</param> | |
public void XYZ2BLH(EllipsoidType type) | |
{ | |
EllipsoidParam(type); | |
// 这里转出来的B L是弧度 | |
this.L = Math.Atan(this.Y2 / this.X2) + Math.PI; | |
this.L = this.L * 180 / Math.PI; | |
// B需要迭代计算 | |
double B2 = Math.Atan(Z2 / Math.Sqrt(X2 * X2 + Y2 * Y2)); | |
double B1; | |
double N; | |
while (true) | |
{ | |
N = a / Math.Sqrt(1 - f * (2 - f) * Math.Sin(B2) * Math.Sin(B2)); | |
B1 = Math.Atan((Z2 + N * f * (2 - f) * Math.Sin(B2)) / Math.Sqrt(X2 * X2 + Y2 * Y2)); | |
if (Math.Abs(B1 - B2) < 1e-12) | |
break; | |
B2 = B1; | |
} | |
this.B = B2 * 180 / Math.PI; | |
double sB = Math.Sin(this.B * d2r); | |
double cB = Math.Cos(this.B * d2r); | |
this.H = this.Z2 / sB - N * (1 - this.e * this.e); | |
} |
3.4 投影转换
此处仅实现了常见的高斯-克里格投影。上述MyPoint
类中增加GaussProjection(EllipsoidType type, ProjectionSetting prjSetting)
函数。
/// <summary> | |
/// 利用高斯投影将指定椭球类型的经纬度坐标转为投影坐标 | |
/// </summary> | |
/// <param name="type">椭球类型</param> | |
/// <param name="prjSetting">投影设置实例</param> | |
public void GaussProjection(EllipsoidType type, ProjectionSetting prjSetting) | |
{ | |
this.EllipsoidParam(type); | |
double l = (this.L - prjSetting.CenterL) / this.rho; | |
double cB = Math.Cos(this.B * this.d2r); | |
double sB = Math.Sin(this.B * this.d2r); | |
double s2b = Math.Sin(this.B * this.d2r * 2); | |
double s4b = Math.Sin(this.B * this.d2r * 4); | |
double s6b = Math.Sin(this.B * this.d2r * 6); | |
double s8b = Math.Sin(this.B * this.d2r * 8); | |
double N = this.a / Math.Sqrt(1 - this.e * this.e * sB * sB); // 卯酉圈曲率半径 | |
double t = Math.Tan(this.B * this.d2r); | |
double eta = this.e2 * cB; | |
double m0 = this.a * (1 - this.e * this.e); | |
double m2 = 3.0 / 2.0 * this.e * this.e * m0; | |
double m4 = 5.0 / 4.0 * this.e * this.e * m2; | |
double m6 = 7.0 / 6.0 * this.e * this.e * m4; | |
double m8 = 9.0 / 8.0 * this.e * this.e * m6; | |
double a0 = m0 + 1.0 / 2.0 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8; | |
double a2 = 1.0 / 2.0 * m2 + 1.0 / 2.0 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8; | |
double a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8; | |
double a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8; | |
double a8 = 1.0 / 128.0 * m8; | |
// X1为自赤道量起的子午线弧长 | |
double X1 = a0 * (this.B * this.d2r) - 1.0 / 2.0 * a2 * s2b + 1.0 / 4.0 * a4 * s4b - 1.0 / 6.0 * a6 * s6b + 1.0 / 8.0 * a8 * s8b; | |
this.Xs = X1 + N / 2 * t * cB * cB * l * l + N / 24 * t * (5 - t * t + 9 * Math.Pow(eta, 2) + 4 * Math.Pow(eta, 4)) * Math.Pow(cB, 4) * Math.Pow(l, 4) | |
+ N / 720 * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(cB, 6) * Math.Pow(l, 6); | |
this.Ys = N * cB * l + N / 6 * (1 - t * t + eta * eta) * Math.Pow(cB, 3) * Math.Pow(l, 3) | |
+ N / 120 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * eta * eta * t * t) * Math.Pow(cB, 5) * Math.Pow(l, 5); | |
this.Hs = this.H; | |
// 假东 假北偏移 | |
this.Xs += prjSetting.PseudoNorth; | |
this.Ys += prjSetting.PseudoEast; | |
} | |
其中,ProjectionSetting
是一个投影参数设置类,独立于MyPoint
类,用于设定中央经线、东偏等投影参数。
internal class ProjectionSetting | |
{ | |
private double _centerL; | |
public double CenterL | |
{ | |
get { return _centerL; } | |
set { _centerL = value; } | |
} | |
private double _centerB; | |
public double CenterB | |
{ | |
get { return _centerB; } | |
set { _centerB = value; } | |
} | |
private double _pseudoEast; | |
public double PseudoEast | |
{ | |
get { return _pseudoEast; } | |
set { _pseudoEast = value; } | |
} | |
private double _pseudoNorth; | |
public double PseudoNorth | |
{ | |
get { return _pseudoNorth; } | |
set { _pseudoNorth = value; } | |
} | |
private double _prjScale; | |
public double PrjScale | |
{ | |
get { return _prjScale; } | |
set { _prjScale = value; } | |
} | |
/// <summary> | |
/// 设置全部的投影参数 | |
/// </summary> | |
/// <param name="centerL"></param> | |
/// <param name="centerB"></param> | |
/// <param name="pseudoEast"></param> | |
/// <param name="pseudoNorth"></param> | |
/// <param name="prjScale"></param> | |
public ProjectionSetting(double centerL, double centerB, | |
double pseudoEast, double pseudoNorth, | |
double prjScale) | |
{ | |
CenterL = centerL; | |
CenterB = centerB; | |
PseudoEast = pseudoEast; | |
PseudoNorth = pseudoNorth; | |
PrjScale = prjScale; | |
} | |
/// <summary> | |
/// 仅设置中央经线和东偏 | |
/// </summary> | |
/// <param name="centerL"></param> | |
/// <param name="pseudoEast"></param> | |
public ProjectionSetting(double centerL, double pseudoEast) | |
{ | |
CenterL = centerL; | |
CenterB = 0.0; | |
PseudoEast = pseudoEast; | |
PseudoNorth = 0.0; | |
PrjScale = 1.0; | |
} | |
/// <summary> | |
/// 默认常用投影参数,中央经线120°,东偏500000 | |
/// </summary> | |
public ProjectionSetting() | |
{ | |
CenterL = 120.0; | |
CenterB = 0.0; | |
PseudoEast = 500000; | |
PseudoNorth = 0.0; | |
PrjScale = 1.0; | |
} | |
} |
3.5 转换2的实现(三参数、七参数)
上述MyPoint
类中增加SevenParamTrans(Datum7Paras datum7Paras)
和TreeParamTrans(Datum3Paras datum3Paras)
函数。
/// <summary> | |
/// 利用7参数进行坐标系之间转换 | |
/// </summary> | |
/// <param name="datum7Paras">7参数实例</param> | |
public void SevenParamTrans(Datum7Paras datum7Paras) | |
{ | |
this.dx = datum7Paras.Dx; | |
this.dy = datum7Paras.Dy; | |
this.dz = datum7Paras.Dz; | |
this.rx = datum7Paras.Rx * 0.0000048481373323; //1 秒=0.0000048481373323 弧度 | |
this.ry = datum7Paras.Ry * 0.0000048481373323; | |
this.rz = datum7Paras.Rz * 0.0000048481373323; | |
this.m = datum7Paras.PPM; | |
this.k = this.m / 1000000; | |
this.X2 = (1 + k) * (this.X + this.rz * this.Y - this.ry * this.Z) + this.dx; | |
this.Y2 = (1 + k) * (-this.rz * this.X + this.Y + this.rx * this.Z) + this.dy; | |
this.Z2 = (1 + k) * (this.ry * this.X - this.rx * this.Y + this.Z) + this.dz; | |
} | |
/// <summary> | |
/// 利用3参数进行坐标系之间转换 | |
/// </summary> | |
/// <param name="datum3Paras">3参数实例</param> | |
public void TreeParamTrans(Datum3Paras datum3Paras) | |
{ | |
this.dx = datum3Paras.Dx; | |
this.dy = datum3Paras.Dy; | |
this.dz = datum3Paras.Dz; | |
this.X2 = this.X + this.dx; | |
this.Y2 = this.Y + this.dy; | |
this.Z2 = this.Z + this.dz; | |
} | |
Datum3Paras
和Datum7Paras
是独立于MyPoint
类,用于设定坐标转换参数。
/// <summary> | |
/// 7参数 | |
/// </summary> | |
internal class Datum7Paras | |
{ | |
private double _dx; | |
public double Dx | |
{ | |
get { return _dx; } | |
set { _dx = value; } | |
} | |
private double _dy; | |
public double Dy | |
{ | |
get { return _dy; } | |
set { _dy = value; } | |
} | |
private double _dz; | |
public double Dz | |
{ | |
get { return _dz; } | |
set { _dz = value; } | |
} | |
private double _rx; | |
public double Rx | |
{ | |
get { return _rx; } | |
set { _rx = value; } | |
} | |
private double _ry; | |
public double Ry | |
{ | |
get { return _ry; } | |
set { _ry = value; } | |
} | |
private double _rz; | |
public double Rz | |
{ | |
get { return _rz; } | |
set { _rz = value; } | |
} | |
private double _ppm; | |
public double PPM | |
{ | |
get { return _ppm; } | |
set { _ppm = value; } | |
} | |
public Datum7Paras(double dx, double dy, double dz, | |
double rx, double ry, double rz, | |
double ppm) | |
{ | |
_dx = dx; | |
_dy = dy; | |
_dz = dz; | |
_rx = rx; | |
_ry = ry; | |
_rz = rz; | |
_ppm = ppm; | |
} | |
} |
internal class Datum3Paras | |
{ | |
private double _dx; | |
public double Dx | |
{ | |
get { return _dx; } | |
set { _dx = value; } | |
} | |
private double _dy; | |
public double Dy | |
{ | |
get { return _dy; } | |
set { _dy = value; } | |
} | |
private double _dz; | |
public double Dz | |
{ | |
get { return _dz; } | |
set { _dz = value; } | |
} | |
public Datum3Paras(double dx, double dy, double dz) | |
{ | |
Dx = dx; | |
Dy = dy; | |
Dz = dz; | |
} | |
} |
3.6 转换5的实现(四参数+高程拟合)
上述MyPoint
类中增加Transform4Para(Trans4Paras transPara)
函数。此处,高程拟合仅实现了已知一个测点的固定改正差。
/// <summary> | |
/// 投影坐标获取后,进一步利用4参数转换坐标 | |
/// </summary> | |
/// <param name="transPara"></param> | |
public void Transform4Para(Trans4Paras transPara) | |
{ | |
var X1 = transPara.Dx; | |
var Y1 = transPara.Dy; | |
var cosAngle = Math.Cos(transPara.A); | |
var sinAngle = Math.Sin(transPara.A); | |
X1 += transPara.K * (cosAngle * this.Xs - sinAngle * this.Ys); | |
Y1 += transPara.K * (sinAngle * this.Xs + cosAngle * this.Ys); | |
this.Xs = X1; | |
this.Ys = Y1; | |
// 固定改正差 | |
this.Hs += transPara.Dh; | |
} |
Trans4Paras
是独立于MyPoint
类,用于设定坐标转换参数。
internal class Trans4Paras | |
{ | |
private double _dx; | |
public double Dx | |
{ | |
get { return _dx; } | |
set { _dx = value; } | |
} | |
private double _dy; | |
public double Dy | |
{ | |
get { return _dy; } | |
set { _dy = value; } | |
} | |
private double _a; | |
public double A | |
{ | |
get { return _a; } | |
set { _a = value; } | |
} | |
private double _k; | |
public double K | |
{ | |
get { return _k; } | |
set { _k = value; } | |
} | |
private double _dh; | |
public double Dh | |
{ | |
get { return _dh; } | |
set { _dh = value; } | |
} | |
public Trans4Paras(double dx, double dy, double a, double k, double dh) | |
{ | |
Dx = dx; | |
Dy = dy; | |
A = a; | |
K = k; | |
Dh = dh; | |
} | |
public Trans4Paras() | |
{ | |
} | |
} |
3.7 调用过程
里面的参数,因为保密原因,做出了随机更改,实际使用时可根据自己情况赋值。
3.7.1 一步法
// 实例化计算参数 | |
MyPoint p = new MyPoint();. | |
p.L=113.256; | |
p.B=31.565; | |
p.H=5.216; | |
// 经纬度转空间坐标 | |
p.BLH2XYZ(); | |
// 实例化七参数 | |
Datum7Paras datum7Paras = new Datum7Paras( | |
489.2994563566, 141.1525159753, 15.74421120568, | |
-0.164423, 4.141573, -4.808299, | |
-6.56482989958); | |
p.SevenParamTrans(datum7Paras); | |
// 空间坐标转回经纬度 | |
p.XYZ2BLH(EllipsoidType.WGS84); | |
// 高斯投影 经纬度转平面坐标 | |
// 实例化投影参数类 | |
ProjectionSetting projectionSetting = new ProjectionSetting(120,500000); | |
p.GaussProjection(EllipsoidType.WGS84, projectionSetting); | |
3.7.2 两步法
// 实例化计算参数 | |
MyPoint p = new MyPoint();. | |
p.SetLBH(113.256,31.565,5.216); | |
// 经纬度转空间坐标 | |
p.BLH2XYZ(); | |
// 实例化七参数 | |
Datum7Paras datum7Paras = new Datum7Paras( | |
489.2994563566, 141.1525159753, 15.74421120568, | |
-0.164423, 4.141573, -4.808299, | |
-6.56482989958); | |
p.SevenParamTrans(datum7Paras); | |
// 空间坐标转回经纬度 | |
p.XYZ2BLH(EllipsoidType.WGS84); | |
// 高斯投影 经纬度转平面坐标 | |
// 实例化投影参数类 | |
ProjectionSetting projectionSetting = new ProjectionSetting(120,500000); | |
p.GaussProjection(EllipsoidType.WGS84, projectionSetting); | |
Trans4Paras transformPara = new(6456.15957352521, -134618.390707439, 0.011104964500129, 1.00002537583871, 5.788); | |
p.Transform4Para(transformPara); | |
4. 总结
至此,关于工程坐标系转化,即GPS接收的WGS84椭球的经纬度转换为地方坐标系的问题,基本全部实现。代码正确性和准确性的验证是与 南方GPS工具箱做对比。例如,采用上述的一步法,在设定好坐标、7参数、投影参数后,计算发现,与南方GPS工具箱在y方向偏差1mm。结果如下图:
-------------------------------------------
个性签名:独学而无友,则孤陋而寡闻。做一个灵魂有趣的人!
如果觉得这篇文章对你有小小的帮助的话,记得在右下角点个“推荐”哦,博主在此感谢!
万水千山总是情,打赏一分行不行,所以如果你心情还比较高兴,也是可以扫码打赏博主,哈哈哈(っ•̀ω•́)っ✎⁾⁾!