WGS84与大地2000坐标转换(Java,C#,Dart)

一、坐标转换的必要性

平面坐标在道路测绘,隧道测量,农业建筑业等室外勘测等方面有着广泛的应用,各行业基本都会涉及到移动端测量之后不能满足屏幕坐标,所以需要经纬度的转换,移动端勘测结果都是WGS84坐标或者GCJ-02格式坐标,而实际工程项目中却需要的是平面坐标[CGCS2000]或[北京坐标]或[西安坐标系],转换过程都是通过高斯投影来计算结果。高斯投影毕竟是将一个不规则的椭圆投影必定会变形,因此在实际应用中,要考虑尽量减小误差。而移动端进行高斯投影的结果误差太大而无法用到实际应用中。所以通过分析高斯投影过程中产生误差的决定因素来减少误杀对于所有行业都显的很必要。

二、我国三大常用坐标系

我国三大常用坐标系区别

• 北京54
• 西安80
• WGS-84

1、北京54坐标系(BJZ54)

北京54坐标系为参心大地坐标系
北京54坐标系，属三心坐标系，长轴6378245m，短轴6356863，扁率1/298. 3;

2、西安80坐标系

1978年4月在西安召开全国天文大地网平差会议，确定重新定位，建立我国新的坐标系。(即1985国家高程基准)。西安80坐标系，属三心坐标系，长轴6378140m。

3、WGS-84坐标系(World Geodetic System) 是-种国际上采用的地心坐标系。

坐标原点为地球质心，其地心空间直角坐标系的Z轴指向国际时间局(BIH)WGS84坐标系，长轴6378137. 000m，短轴6356752.314，扁率1/298. 257223563。由于采用的椭球基准不一样，并且由于投影的局限性,使的全国各地并不存在一至的转换参数。对于这种转换由于量较大，有条件的话，一般都采用GPS联测己知点，应用GPS软件自动完成坐标的转换。当然若条件不许可，且有足够的重
合点，也可以进行人工解算。

4、2000 国家大地坐标系

英文缩写为CGCS2000。2000国家大地坐标系是全球地心坐标系在我国的具
体体现，其原点为包括海洋和大气的整个地球的质量中心。2000国家大地坐标
系采用的地球椭球参数如下:
长半轴a=6378137m， 扁率f=1/298.257222101,
地心引力常数GM=3. 986004418X 1014m3s-2
自转角速度w=7.292115X 10-5rads- 1

参考资料:

2000国家大地坐标系
中央子午线经度的计算
WGS84 经纬度坐标到北京 54 高斯投影坐标的转换

5，注意点:

1.明确是WGS8和那个坐标系转换，确定参数值影响很大

大地2000的参数:
double a = 6378137; //椭球长半轴
double b = 6356752.3142451795; //椭球短半轴
double e = 0.081819190842621; //第一偏心率
double eC = 0.0820944379496957; //第二偏心率
对应的其他坐标系自己可以百度。

2.明确自己项目点经纬度所在的带号计算出带所在的中央带经度

这里自己可以百度看看，我国有两种分带，3度带，6度带。
360度=3度X. => X=120份
360度=6度Y => Y=60份
360度=2π弧度
每度=2π/360=π/180

地球是圆的，如上图东西分界线为0度，左右分为东经西经，如果是三分带，那么可以知道，-1.5度(即西经)和 1.5度(即东经)之间的夹角就是3度。1.5到4.5之间和-4.5度(西经)到-1.5(西经)之间是关于子午线对称的，所以经度纬度相等的有两个。我们大地坐标一般避免经度为负都会将平面地图坐标系向左平移500公里。如下图。

所谓的项目点所在的中央经度。例如3分带划分，东经1.5度和4.5度之间就是3，同样3.5和6.5之间就是5, 我们项目点东经3.53226,3.74,3.2334348,4.02323,5.023232…那么它的中央经度为带的中点 3。转换为计算机代码,Math.roun(经度/3)*3.

Dart代码

import 'dart:math';
import 'package:luhenchang_plugin/map/bean/latlng.dart';
import 'bean/tuple.dart';
/**
 * 版权：渤海新能 版权所有
 *
 * @author feiWang
 * 版本：1.5
 * 创建日期：2020/10/16
 * 描述：AndroidDaggerStudy
 * E-mail : 1276998208@qq.com
 * CSDN:https://blog.csdn.net/m0_37667770/article
 * GitHub:https://github.com/luhenchang
 */
class LonLatAndXY {
   
  //@将WGS84经纬度转为大地2000坐标。我们是国家电网项目数据很精确的了。
  //@param B 纬度
  //@param L 经度
  //@param degree //
  //@param withBand 默认=false
  //@return
  Tuple gps84ToXY(double B, double L, double degree, {
   withBand = false}) {
   
    List<double> xy = new List();
    xy.add(0);
    xy.add(0);
    double a = 6378137; //椭球长半轴
    double b = 6356752.3142451795; //椭球短半轴
    double e = 0.081819190842621; //第一偏心率
    double eC = 0.0820944379496957; //第二偏心率

    double L0 = 0; //中央子午线经度
    int n = 0; //带号
    if (degree == 6) {
   
      //6度
      n = ((L + degree / 2) / degree).round().toInt();
      L0 = degree * n - degree / 2;
    } else {
   
      //3度
      n = (L / degree).round().toInt();
      L0 = degree * n;
    }
    //开始计算
    double radB = B * pi / 180; //纬度(弧度)
    double radL = L * pi / 180; //经度(弧度)
    double deltaL = (L - L0) * pi / 180; //经度差(弧度)
    double N = a * a / b / sqrt(1 + eC * eC * cos(radB) * cos(radB));
    double C1 = 1.0 +
        3.0 / 4 * e * e +
        45.0 / 64 * pow(e, 4) +
        175.0 / 256 * pow(e, 6) +
        11025.0 / 16384 * pow(e, 8);
    double C2 = 3.0 / 4 * e * e +
        15.0 / 16 * pow(e, 4) +
        525.0 / 512 * pow(e, 6) +
        2205.0 / 2048 * pow(e, 8);
    double C3 = 15.0 / 64 * pow(e, 4) +
        105.0 / 256 * pow(e, 6) +
        2205.0 / 4096 * pow(e, 8);
    double C4 = 35.0 / 512 * pow(e, 6) + 315.0 / 2048 * pow(e, 8);
    double C5 = 315.0 / 131072 * pow(e, 8);
    double t = tan(radB);
    double eta = eC * cos(radB);
    double X = a *
        (1 - e * e) *
        (C1 * radB -
            C2 * sin(2 * radB) / 2 +
            C3 * sin(4 * radB) / 4 -
            C4 * sin(6 * radB) / 6 +
            C5 * sin(8 * radB));

    xy[0] = X +
        N *
            sin(radB) *
            cos(radB) *
            pow(deltaL, 2) *
            (1 +
                pow(deltaL * cos(radB), 2) *
                    (5 - t * t + 9 * eta * eta + 4 * pow(eta, 4)) /
                    12 +
                pow(deltaL * cos(radB), 4) *
                    (61 - 58 * t * t + pow(t, 4)) /
                    360) /
            2;
    xy[1] = N *
            deltaL *
            cos(radB) *
            (1 +
                pow(deltaL * cos(radB), 2) * (1 - t * t + eta * eta) / 6 +
                pow(deltaL * cos(radB), 4) *
                    (5 -
                        18 * t * t +
                        pow(t, 4) -
                        14 * eta * eta -
                        58 * eta * eta * t * t) /
                    120) +
        500000; // +n * 1000000;
    return new Tuple(xy[0], xy[1]);
  }

  //@将大地2000转为WGS84
  //高斯投影反算为大地平面。
  // x，y ，高斯平面坐标点
  //L0 通过经纬度来获取中央带所在带的角度
  //return B纬度 , L经度
  LatLng xyTowgs84(double x, double y, double L0) {
   
    //中央子午线经度
    //WGS-84   椭球体参数
    double a = 6378137.0; //major semi axis
    double efang = 0.0066943799901413; //square of e
    double e2fang = 0.0067394967422764; //suqre of e2

    y = y - 500000;

    //主曲率计算
    double m0, m2, m4, m6, m8;
    m0 = a * (1 - efang);
    m2 = 3.0 / 2.0 * efang * m0;
    m4 = efang * m2 * 5.0 / 4.0;
    m6 = efang * m4 * 7.0 / 6.0;
    m8 = efang * m6 * 9.0 / 8.0;

    //子午线曲率计算
    double a0, a2, a4, a6, a8;
    a0 = m0 + m2 / 2.0 + m4 * 3.0 / 8.0 + m6 * 5.0 / 16.0 + m8 * 35.0 / 128.0;
    a2 = m2 / 2.0 + m4 / 2.0 + m6 * 15.0 / 32.0 + m8 * 7.0 / 16.0;
    a4 = m4 / 8.0 + m6 * 3.0 / 16.0 + m8 * 7.0 / 32.0;
    a6 = m6 / 32.0 + m8 / 16.0;
    a8 = m8 / 128.0;

    double X = x;
    double FBf = 0;
    double Bf0 = X / a0, Bf1 = 0;

    //计算Bf的值，直到满足条件
    while ((Bf0 - Bf1) >= 0.0001) {
   
      Bf1 = Bf0;
      FBf = -a2*sin(2*Bf0)/2 +a4*