工程坐标转换方法C#代码实现
目录
- 1. 前言
- 2. 计算总体框架
- 3. C#代码实现
- 3.1 整体类的构建
- 3.2 椭球参数赋值
- 3.3 转换1、3(大地经纬度坐标与地心地固坐标的转换)
- 3.4 投影转换
- 3.5 转换2的实现(三参数、七参数)
- 3.6 转换5的实现(四参数+高程拟合)
- 3.7 调用过程
- 3.7.1 一步法
- 3.7.2 两步法
- 4. 总结
1. 前言
在前面的文章中系统的阐述了工程坐标的转换类别和转换的方法。关于转换代码实现,有很多的类库:
- GDAL
- SharpProj - Providing OSGEO PROJ for .Net (Core)
- ProjNet (for GeoAPI)
这里针对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/1000000private 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;}// 西安 80else if (type == EllipsoidType.西安80){this.a = 6378140;this.f = 1 / 298.257;}// 北京 54else 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。结果如下图:
工程坐标转换方法C#代码实现相关推荐
- 经纬度坐标转换为工程坐标
1. 绪论 在施工和工程测量时,经常需要将GPS坐标转换为工程中所使用的坐标.在Bing上检索到的类似问题,基本描述为两个坐标系的转换,但实际上并非如此. 本文将详细解释转换过程和转换方法. 1.1 ...
- 使用Eclipse可以方便的统计工程或文件的代码行数,
使用Eclipse可以方便的统计工程或文件的代码行数,方法如下: 1.点击要统计的项目或许文件夹,在菜单栏点击Search,然后点击File... 2.选中正则表达式(Regular expressi ...
- BestMPRBaseVtk-003-修改工程,搬运官方代码并尝试理解-2
修改工程,搬运官方代码并尝试理解-2 接着上篇,今天接着搞官方的源码,太难了,真想回家卖红薯去,你们说靠卖红薯可以养家糊口吗? 文章目录 修改工程,搬运官方代码并尝试理解-2 1 setRende ...
- 机器学习和特征工程理论与python代码实现 晓物智联
文章来源于:http://www.52phm.cn/blog/detail/23 最初来源于本人的kesci专栏 课题:特征工程理论及代码实现 日期:2019.9.21 作者:小知同学 描述:本篇比较 ...
- 相机计算坐标公式_相机采样点的坐标转换方法与流程
本发明涉及数字图像处理与计算机视觉领域,尤其涉及一种相机采样点的坐标转换方法. 背景技术: 广义上说振动是指描述系统状态的参量(如位移.电压)在其基准值上下交替变化的过程:狭义的指机械振动,即力学系统 ...
- iOS小技能:-fobjc-arc和 -fno-objc-arc 的使用(在非ARC工程中集成ARC代码、在ARC工程中集成非ARC的第三方代码)
文章目录 前言 I ARC 简介 1.1 ARC的规则 1.2 OC中有强参照strong和弱参照weak. 1.3 ARC只能工作于OC. 前言 在非ARC工程中集成ARC代码: 使用-fobjc- ...
- 【JS】1070- 8个工程必备的JavaScript代码片段(建议添加到项目中)
8个工程必备的JavaScript代码片段,听过这样起博客标题可以提高阅读量.???? 最近写博客好累,让8月征文活动搞的,今天水一篇好了,麻烦不要给我点赞,不想看到消息通知的小红点. 1. 获取文件 ...
- 16个工程必备的JavaScript代码片段
作者:_红领巾 https://juejin.cn/post/7000919400249294862 1. 下载一个excel文档 同时适用于word,ppt等浏览器不会默认执行预览的文档,也可以用于 ...
- 8个工程必备的JavaScript代码片段(建议添加到项目中)
点击上方 前端瓶子君,关注公众号 回复算法,加入前端编程面试算法每日一题群 8个工程必备的JavaScript代码片段,听过这样起博客标题可以提高阅读量.???? 最近写博客好累,让8月征文活动搞的, ...
最新文章
- 查看临界区等待线程数量
- 干货 | 深入仓储管理系统你需要了解的15件事
- 2.Hadoop的学习(Ubuntu的目录及权限)
- python爬新闻并保存csv_用python爬取内容怎么存入 csv 文件中
- aws lambda_AWS Lambda事件源映射:使您的触发器混乱无序
- 列出一个目录中所有文件及大小
- java中stack集合框架
- linux7基础——给用户添加sudo权限
- 第十周学习总结--助教
- 接入Google Play SDK
- 中国农用喷雾机市场趋势报告、技术动态创新及市场预测
- [译] 用 Swift 创建自定义的键盘
- Java 生成二维码实战
- 1 dhcp服务器的配置文件,Linux1 DHCP服务器配置 主配置文件(dhcpd.conf)
- AJAX框架大全 (AJAX Frameworks)
- 20190223深信服测试一面回顾
- 计算机职业规划500字中专,计算机中专生职业规划范文500字中专生职业生涯规划书范文.doc...
- win10下git报fatal: open /dev/null or dup failed解决办法(附null.sys文件下载)
- hive 正则表达式 过滤字符串里的中文
- 欧美是怎么做创新的?
热门文章
- 最好的护眼灯是什么牌子?央视315护眼灯合格名单
- 从青腾爬向宇宙:科技巨头与少年科学家的故事
- 网络教学资源平台设计与实现--公告发布系统数据表
- python爬虫爬图片教程_Python爬虫爬图片需要什么
- 纽曼u盘无法打开的解决方法 量产工具
- python word操作添加超链接_使用pythondocx在MSWord中添加超链接
- WINPE U盘版制作-老毛桃版本
- android stm32 蓝牙模块,STM32+USART+蓝牙模块(BT04)
- Realme GT Neo2T ROOT 解锁BL教程
- 24V转3.3V电路设计