坐标转换-大地转高斯平面&平面坐标转换

最新推荐文章于 2024-05-09 16:58:25 发布

alan1914

最新推荐文章于 2024-05-09 16:58:25 发布

阅读量1.6w

点赞数 12

分类专栏： Java 文章标签： java 高斯正反算平面坐标转换

本文链接：https://blog.csdn.net/db0n15/article/details/84567846

版权

Java 专栏收录该内容

13 篇文章 0 订阅

订阅专栏

楔子

以前呢，总感觉大学老师教的都没用，基本上都用不上。得了，这两天碰上了。
我们公司给甲方做了个小程序，因为业务原因必须使用深圳独立坐标系。
一个业务需求-导航。因为我们是获取手机的gps坐标，我们起先是使用的甲方提供的坐标转换服务，发现84大地转到深圳独立坐标（深圳高斯平面坐标）时精度误差极大。于是甲方给我们提供了部分控制点信息，让我们自己去完成坐标的转换，不再使用甲方提供的服务。
一把眼泪，大一学的地图学和测量学早把这些东西忘光了。
一步步来吧，随让自己当初不好好学呢。

数据

49990.201  63839.358           2521448.83   496112.77             113.86475295492812            22.79054504954949
49839.360  97644.067           2511136.78   476638.23             113.86988381063978            22.779607662363563

注：这是这里修改了数据反正大概长这样大概位数是这么多至于具体的值。我已经改得面目全非了
用于四参数计算已经足够了

解决思路

这些都是什么鬼数据？
这些数据都怎么转换？
转换公式是什么?
转换代码如何编写？

ps:折腾了快一天后，事实证明，一开始设计的解决问题思路是对的
转换逻辑 84大地坐标系 <–>84平面直角坐标系 <–>深圳平面直角坐标系

数据分析

// 深圳高斯平面数据(至于其他的参数就不给了)
49990.201  63839.358
// 84高斯平面数据
 2521448.83   496112.77
// 84大地坐标
113.86475295492812  22.79054504954949

软件转换+手工验算

首先呢，理论知识都忘的差不多了，还是先计算清楚把。
这里使用的coord软件
在这里插入图片描述

第一步平面坐标到平面坐标的转换

设置->计算四参数
输入上面两组数据源单击 -> 计算
单击 ->导出得到如下四参数
在这里插入图片描述

注:这里是深圳平面转84平面的四参数
如何要获取 84平面转深圳平面的四参数重新计算即可

公式如下

在这里插入图片描述

代码如下

   /**
     * 深圳平面坐标转84平面
     *
     * @param x 纵坐标
     * @param y 横坐标
     * @return
     */
    public double[] convert54to84(double x, double y) {
        double[] coord = {0.0, 0.0};
        double X;
        double Y;

        // Coord gm 控制点计算结果
        // 核心参数
        final double DX = 2470000.759939;
        final double DY = 390000.400211;
        final double T = -0.0174591152;
        final double K = 0.999887111016;

        double a = DX;
        double b = DY;
        double t = T;
        double k = K;

		// 核心代码
        X = a + x * k * Math.cos(t) - y * k * Math.sin(t);
        Y = b + y * k * Math.cos(t) + x * k * Math.sin(t);

        coord[0] = X;
        coord[1] = Y;
        return coord;
    }

注：这里84平面转深圳平面是一样的只是核心参数不同罢了
核心参数已经改的面目全非反正大概数据长这样罢了

第二步大地坐标于高斯平面坐标的转换

其实这个用专业的术语是高斯正反算

高斯正算
大地坐标转平面坐标
高斯反算
平面坐标转大地坐标

注：别看了，我正儿八经也没全看懂，主要是在网上找到了现成的代码。具体要看公式的话这里看肯定是不行的，建议下本大地测量学的书，那里会手把手教你如何计算出来。记得大一还是大二的时候老师让我们算过，我可能算过也可能没算过，反正快4年了我忘了。

代码


/**
 * 由经纬度反算成高斯投影坐标(高斯正算)
 *
 * @param longitude
 * @param latitude
 */
public double[] BLToGauss(double longitude, double latitude) {

    int ProjNo = 0;

    // 带宽
    int ZoneWide = 3;

    double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
    double a, f, e2, ee, NN, T, C, A, M, iPI;

    // 3.1415926535898/180.0;
    iPI = 0.0174532925199433;

    // 84年北京坐标系参数
    a = 6378137;
    f = 1.0 / 298.2572236;

    ProjNo = 0;
    longitude0 = ProjNo * ZoneWide ;
    longitude0 = 114 * iPI;

    // 经度转换为弧度
    longitude1 = longitude * iPI;

    // 纬度转换为弧度
    latitude1 = latitude * iPI;

    e2 = 2 * f - f * f;
    ee = e2 * (1.0 - e2);
    NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1) * Math.sin(latitude1));
    T = Math.tan(latitude1) * Math.tan(latitude1);
    C = ee * Math.cos(latitude1) * Math.cos(latitude1);
    A = (longitude1 - longitude0) * Math.cos(latitude1);
    M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude1));
    xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
    yval = M + NN * Math.tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
    X0 = 1000000 * (ProjNo ) + 500000;
    Y0 = 0;
    xval = xval + X0;
    yval = yval + Y0;
    return new double[] { xval, yval };
}

/**
 * 由高斯投影坐标反算成经纬度 84坐标系
 *
 * @param X
 * @param Y
 * @return
 */
public double[] GaussToBL(double X, double Y) {

    int ProjNo;
    // 带宽
    int ZoneWide;
    double[] output = new double[2];
    // latitude0,
    double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
    double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
    //3.1415926535898/180.0;
    iPI = 0.0174532925199433;
    a = 6378137;
    f = 1 / 298.2572236;
    // 3度带宽
    ZoneWide = 3;
    // 查找带号
    ProjNo = 0;
    // 中央经线
    longitude0 = 114 * Math.PI / 180;

    X0 = ProjNo * 1000000L + 500000L;
    Y0 = 0;
    xval = X - X0;
    // 带内大地坐标
    yval = Y - Y0;
    e2 = 2 * f - f * f;
    e1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2));
    ee = e2 / (1 - e2);
    M = yval;
    u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
    fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.sin(4 * u) + (151 * e1 * e1 * e1 / 96) * Math.sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.sin(8 * u);
    C = ee * Math.cos(fai) * Math.cos(fai);
    T = Math.tan(fai) * Math.tan(fai);
    NN = a / Math.sqrt(1.0 - e2 * Math.sin(fai) * Math.sin(fai));
    R = a * (1 - e2) / Math.sqrt((1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)));
    D = xval / NN;
    // 计算经度(Longitude) 纬度(Latitude)
    longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / Math.cos(fai);
    latitude1 = fai - (NN * Math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
    // 转换为度 DD
    output[0] = longitude1 / iPI;
    output[1] = latitude1 / iPI;
    return output;
}

注：我在这里埋了个坑，网上找的代码是6度带的。我这里需要的是三度带的，而且我的业务区域很小，只需要114度的中央经线即可，我就偷懒写死了。有会计算三度带的小伙伴可以@我。我从书上看的三度带的计算公式很简单，问题是套在代码里不好使。
回答：今天研究了一下，发现我的y坐标少了两位。所以我套用网上的代码不好使。
我这里的y坐标（横坐标）496112.77
真实的y坐标(横坐标)应该是 38496112.77

Es6版本高斯正反算坐标转换和四参数平面坐标转换

注意：四参数平面坐标转换已经将核心参数改的面目全非了根据自己的控制点填写 const参数

/**
 * 深圳54平面坐标转84平面
 *
 * @param x 纵坐标
 * @param y 横坐标
 * @return
 */
function convert54GaussTo84Gauss (x, y) {
  // Coord gm 控制点计算结果
  const DX = 2470000.000
  const DY = 390000.000
  const T = -0.01900
  const K = 0.90000000

  let coord = [0.0, 0.0]
  let X
  let Y

  let a = DX
  let b = DY
  let t = T
  let k = K

  X = a + x * k * Math.cos(t) - y * k * Math.sin(t)
  Y = b + y * k * Math.cos(t) + x * k * Math.sin(t)

  coord[0] = X
  coord[1] = Y
  console.log('转换前：' + x + ',' + y + '\t转换后：' + coord)
  return coord
}

/**
 * 84高斯平面坐标转深圳54高斯平面坐标
 *
 * @param x 纵坐标
 * @param y 横坐标
 * @return
 */
function convert84GaussTo54Gauss (x, y) {
  // Coord gm 控制点计算结果
  const DX = -2470000.000
  const DY = -410000.000
  const T = 0.01900
  const K = 0.90000000


  let coord = [0.0, 0.0]
  let X
  let Y

  let a = DX
  let b = DY
  let t = T
  let k = K

  X = a + x * k * Math.cos(t) - y * k * Math.sin(t)
  Y = b + y * k * Math.cos(t) + x * k * Math.sin(t)

  coord[0] = X
  coord[1] = Y
  console.log('转换前：' + x + ',' + y + '\t转换后：' + coord)
  return coord
}

/**
 * 由经纬度反算成高斯投影坐标
 *
 * @param longitude 经度
 * @param latitude 纬度
 */
function convert84BLToGauss (longitude, latitude) {
  // 带宽
  const ZONE_WITH = 3
  // 带号
  // 三度带计算公式 带号 = （经度 + 1.5度） /  带宽
  // 注意 3度带和6度带的 带号和中央经线的计算方式不同
  const PROJ_NO = Math.floor(((longitude + 1.5) / ZONE_WITH))
  // 中央经线（弧度制）
  // 三度带计算公式： 中央经线(弧度制) = 带号 * 带宽 * (π/180)
  const longitude0 = (PROJ_NO * ZONE_WITH) * (Math.PI / 180)
  // 长半径（84）
  const a = 6378137
  // 偏率（84）
  const f = 1 / 298.2572236

  let X0, Y0, xval, yval
  let e2, ee, NN, T, C, A, M

  // 经度转换为弧度
  longitude = longitude * (Math.PI / 180)
  // 纬度转换为弧度
  latitude = latitude * (Math.PI / 180)

  e2 = 2 * f - f * f
  ee = e2 * (1.0 - e2)
  NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude) * Math.sin(latitude))
  T = Math.tan(latitude) * Math.tan(latitude)
  C = ee * Math.cos(latitude) * Math.cos(latitude)
  A = (longitude - longitude0) * Math.cos(latitude)
  M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.sin(2 * latitude) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude))
  xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120)
  yval = M + NN * Math.tan(latitude) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720)
  X0 = 1000000 * PROJ_NO + 500000
  Y0 = 0
  xval = xval + X0
  yval = yval + Y0
  let coord = [xval, yval]
  console.log('转换前：' + longitude * 180 / Math.PI + ',' + latitude * 180 / Math.PI + '\t转换后：' + coord)
  return coord
}

/**
 * 由高斯投影坐标反算成经纬度 84坐标系
 *
 * @param x 纵坐标
 * @param y 横坐标
 * @return
 */
function convert84GaussToBL (x, y) {
  // 带宽
  const ZONE_WITH = 3
  // 带号
  // 三度带计算公式 带号 = 横坐标前两位
  // 注意 3度带和6度带的 带号和中央经线的计算方式不同
  const PROJ_NO = Math.floor(y / 1000000)
  // 中央经线（弧度制）
  // 三度带计算公式： 中央经线(弧度制) = 带号 * 带宽 * (π/180)
  const longitude0 = (PROJ_NO * ZONE_WITH) * (Math.PI / 180)
  // 长半径（84）
  const a = 6378137
  // 偏率（84）
  const f = 1 / 298.2572236

  let longitude, latitude, X0, Y0, xval, yval
  let e1, e2, ee, NN, T, C, M, D, R, u, fai

  X0 = 0
  Y0 = PROJ_NO * 1000000 + 500000
  xval = x - X0
  yval = y - Y0
  e2 = 2 * f - f * f
  e1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2))
  ee = e2 / (1 - e2)
  M = xval
  u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256))
  fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.sin(4 * u) + (151 * e1 * e1 * e1 / 96) * Math.sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.sin(8 * u)
  C = ee * Math.cos(fai) * Math.cos(fai)
  T = Math.tan(fai) * Math.tan(fai)
  NN = a / Math.sqrt(1.0 - e2 * Math.sin(fai) * Math.sin(fai))
  R = a * (1 - e2) / Math.sqrt((1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)) * (1 - e2 * Math.sin(fai) * Math.sin(fai)))
  D = yval / NN
  // 计算经度(Longitude) 纬度(Latitude)
  longitude = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / Math.cos(fai)
  latitude = fai - (NN * Math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720)
  // 转换为度
  let coord = [0.0, 0.0]
  coord[0] = longitude * 180 / Math.PI
  coord[1] = latitude * 180 / Math.PI
  console.log('转换前：' + x + ',' + y + '\t转换后：' + coord)
  return coord
}

console
在这里插入图片描述

这里参考了一些其他的博文：

https://blog.csdn.net/u010534192/article/details/77196064 坐标转换流程与公式七参数四参数
https://blog.csdn.net/m0_37970224/article/details/79677868 坐标高斯正反算

alan1914

关注

12
点赞
踩
94

收藏

觉得还不错? 一键收藏
6
评论
坐标转换-大地转高斯平面&平面坐标转换

楔子以前呢，总感觉大学老师教的都没用，基本上都用不上。得了，这两天碰上了。我们公司给甲方做了个小程序，因为业务原因必须使用深圳独立坐标系。一个业务需求-导航。因为我们是获取手机的gps坐标，我们起先是使用的甲方提供的坐标转换服务，发现84大地转到深圳独立坐标（深圳高斯平面坐标）时精度误差极大。于是甲方给我们提供了部分控制点信息，让我们自己去完成坐标的转换，不再使用甲方提供的服务。一把眼泪，...
复制链接

扫一扫