首页 > 编程语言 >工程坐标转换方法C#代码实现

工程坐标转换方法C#代码实现

时间：2022-12-05 13:46:06浏览次数：64  
标签：set get C# double private 转换方法 坐标 public Math

工程坐标转换方法C#代码实现  

 

目录

 

1. 前言

在前面的文章中系统的阐述了工程坐标的转换类别和转换的方法。关于转换代码实现，有很多的类库：

这里针对GPS接收的WGS84椭球的经纬度转换为地方坐标系的问题，利用C#，对工程坐标转换方法和步骤做出详细的解答。不基于任何类库和函数库，也未使用矩阵库，可以便利的将代码移植到任何语言。

2. 计算总体框架

根据上一篇文章中对七参数、四参数、高程拟合在坐标转换的作用和使用条件的阐述，我们可以将上一篇文章第7节的总结图，按照计算的流程重新绘制。
image

根据上图可知，预将WGS84椭球的GPS坐标需要经过5次转换。其中，

  1. 转换1、转换3在charlee44的博客：大地经纬度坐标与地心地固坐标的转换中详细讲解了，并且有C++代码的实现，利用C#重构即可。
  2. 转换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。结果如下图：
image

作者：Aidan 出处：http://www.cnblogs.com/AidanLee/

-------------------------------------------

个性签名：独学而无友，则孤陋而寡闻。做一个灵魂有趣的人！

如果觉得这篇文章对你有小小的帮助的话，记得在右下角点个“推荐”哦，博主在此感谢！

 

万水千山总是情，打赏一分行不行，所以如果你心情还比较高兴，也是可以扫码打赏博主，哈哈哈(っ•̀ω•́)っ✎⁾⁾！

标签：set,get,C#,double,private,转换方法,坐标,public,Math
From： https://www.cnblogs.com/sexintercourse/p/16952061.html

相关文章

  • HACKTHEBOX VM LIST
    LinuxBoxes:WindowsBoxes:MorechallengingthanOSCP,butgoodpractice:LamelegacyJeeves[Windows]brainfuckBlueBart[Windows]shockerDevel......
  • PHP Codeigniter 中的_remap方法简介
    在Codeigniter中，其实有个容易忽略掉的不错的控制器层处理的方法，那就是remap,这里简单介绍下。其实在CI的URL控制方式表示中，比如：example.com/inde......
  • FP6397S5 高效、高频同步DC-DC降压变频器IC
    FP6397是一种高效、高频同步DC-DC降压变频器。100%占空比功能提供了低退出操作，延长了便携式系统的电池寿命。内部同步开关提高了效率，并消除了对外部肖特基二极管的需要。......
  • 《DFZU2EG_4EV MPSoc之FPGA开发指南》第七章 Verilog HDL语法
    VerilogHDL语法​VerilogHDL（HardwareDescriptionLanguage）是在用途最广泛的C语言的基础上发展起来的一种硬件描述语言，具有灵活性高、易学易用等特点。VerilogHDL可以在......
  • make distclean
    makedistcleanmakedistcleanmakedistcleanmakedistcleanopenwrt清理menuconfigkernel_menuconfig的选择 makeclean清理bin目录下生成的烧录镜像文件makedi......
  • Redis 利用 incr 和 expire 来限流, 并发导致过期时间失效问题
    当某一个接口需要限流时,可以采用redis的incr来递增,记录访问次数,以及expire来设置失效时间.大概的代码如下:r=redis.Redis.connect()key="linyk3"......
  • 推荐两个不错的GC分析工具
    推荐两个不错的GC分析工具，一个是在线的：​​https://gceasy.io/​​，另外一个是​​https://github.com/chewiebug/GCViewer​​上的，把这个下载下来maven打包，mvncleanpackag......
  • (收藏）javascript变量提升
    函数和变量的声明总是会被解析器悄悄地被“提升”到方法体的最顶部上面是javascript变量提升的概念。为了更好地理解“变量提升”，我们先来看一段代码：(function(){cons......
  • scrapy3在python2,python3共存下的使用
    因为安装了PYTHON2,PYTHON3，之前的SCRAPY在PYTHON2下是可以的，但在3下运行失败，关联的还是2，原来要在PYTHON3运行的时候，要用全路径：执行scrapy命令(假设......
  • 题解 [AGC047C] Product Modulo
    显然不能暴力算两两的乘积，而积取模而结果不取模提示我们模数肯定有用。所有为\(0\)的\(a_i\)对答案不会产生任何贡献，可以直接删除，下文不再考虑这种情况。同时我们约定......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99