高斯正反算—投影坐标转大地坐标、大地坐标转投影坐标(附有完整代码及测试结果)

        本文意在介绍高斯正反算的基本原理和代码实现!针对高斯正反算网上给出的方法很多,但是我试了之后发现多少都有些问题:或公式原理问题、或精度问题!

        通过查找资料对其进行了总结与测试!原理!代码!测试结果!本文都一一给出,此外本文还对常见坐标系:北京54、西安80、WGS84、CGCS2000等坐标系的高斯转换都做出了实现并使用 QT 进行封装可视化!

        对比公式和代码实现可快速理解高斯正反算!

        精度可达 0.0001!堪称最全公式原理总结!

        坐标转换整套流程包括:像素坐标转投影坐标、投影坐标转大地坐标、大地坐标转空间直角坐标、七参数转换、空间直角坐标转大地坐标、大地坐标转投影坐标、投影坐标转像素坐标; 本人均已实现,且每一个环节都已经过测试、如有需要欢迎在下方留言评论!!!

一、常用椭球参数

北京54坐标系

西安80坐标系

WGS坐标系

CGC2000坐标系

a

6378245.0000000000

6378140.0000000000

6378137.0000000000

6378137.0000000000

b

6356863.0187730473

6356755.2881575287

6356752.3142

6356752.314

f

1/298.3

1/298.257

1/298.257223563

1/298.257222101

c

6399698.9017827110

6399596.6519880105

6399593.6258

6399593.6259

e1

0.006693421622966

0.006694384999588

0.00669437999013

0.00669438002290

e2

0.006738525414683

0.006739501819473

0.00673949674227

0.00673949677548

二、高斯正反算原理 

1、高斯正算(大地坐标转投影坐标) 

 2、高斯反算(投影坐标转大地坐标)

三、高斯正反算代码实现 

#include<iostream>
#include<cmath>
#include "stdio.h"

#define pi 3.141592653589793238463
#define p0 206264.8062470963551564

//wgs84参考椭球
const double e = 0.00669438002290;
const double e1 = 0.00673949677548;
const double b = 6356752.3141;
const double a = 6378137.0;

using namespace std;


//大地坐标转投影坐标
void DadiPoint2ProjectPoint(double B, double L)
{
    //把度转化为弧度
    B = B * pi / 180;
    L = L * pi / 180;

    double N, t, n, c, V, Xz, m1, m2, m3, m4, m5, m6, a0, a2, a4, a6, a8, M0, M2, M4, M6, M8, x0, y0, l;

    int L_num;
    double L_center;

    //中央子午线经度,6°带
    L_num = (int)(L * 180 / pi / 6.0) + 1;
    L_center = 6 * L_num - 3;

    //中央子午线经度,3°带
    //L_num = (int)(L * 180 / pi / 3.0 + 0.5);
    //L_center = 3 * L_num;				  

    l = (L / pi * 180 - L_center) * 3600; //求带号、中央经线、经差

    M0 = a * (1 - e);
    M2 = 3.0 / 2.0 * e * M0;
    M4 = 5.0 / 4.0 * e * M2;
    M6 = 7.0 / 6.0 * e * M4;
    M8 = 9.0 / 8.0 * e * M6;

    a0 = M0 + M2 / 2.0 + 3.0 / 8.0 * M4 + 5.0 / 16.0 * M6 + 35.0 / 128.0 * M8;
    a2 = M2 / 2.0 + M4 / 2 + 15.0 / 32.0 * M6 + 7.0 / 16.0 * M8;
    a4 = M4 / 8.0 + 3.0 / 16.0 * M6 + 7.0 / 32.0 * M8;
    a6 = M6 / 32.0 + M8 / 16.0;
    a8 = M8 / 128.0;

    Xz = a0 * B - a2 / 2.0 * sin(2 * B) + a4 / 4.0 * sin(4 * B) - a6 / 6.0 * sin(6 * B) + a8 / 8.0 * sin(8 * B);  //计算子午线弧长
    c = a * a / b;
    V = sqrt(1 + e1 * cos(B) * cos(B));
    N = c / V;
    t = tan(B);
    n = e1 * cos(B) * cos(B);

    m1 = N * cos(B);
    m2 = N / 2.0 * sin(B) * cos(B);
    m3 = N / 6.0 * pow(cos(B), 3) * (1 - t * t + n);
    m4 = N / 24.0 * sin(B) * pow(cos(B), 3) * (5 - t * t + 9 * n);
    m5 = N / 120.0 * pow(cos(B), 5) * (5 - 18 * t * t + pow(t, 4) + 14 * n - 58 * n * t * t);
    m6 = N / 720.0 * sin(B) * pow(cos(B), 5) * (61 - 58 * t * t + pow(t, 4));
    x0 = Xz + m2 * l * l / pow(p0, 2) + m4 * pow(l, 4) / pow(p0, 4) + m6 * pow(l, 6) / pow(p0, 6);
    y0 = m1 * l / p0 + m3 * pow(l, 3) / pow(p0, 3) + m5 * pow(l, 5) / pow(p0, 5);   //计算x y坐标

    double x = x0;
    //double y = y0 + 500000 + 1000000 * L_num;    //化为国家统一坐标
    double y = y0 + 500000;     //化为国家统一坐标

    cout << "方法一 x=" << x << endl;
    cout << "方法一 y=" << y << endl;
}


//投影坐标转大地坐标
void ProjectPoint2DadiPoint(double x, double y, double l0)
{
    //l0为中央经度
    double Bf, B0, FBf, M, N, V, t, n, c, y1, n1, n2, n3, n4, n5, n6, a0, a2, a4, a6, M0, M2, M4, M6, M8, l;

    int L_num, L_center;

    L_num = (int)(x / 1000000.0);
    y1 = y - 500000;
    //y1 = y - 500000 - L_num * 1000000;

    //L_center = ((L_num + 1) * 6 - 3)*pi*180;		//中央子午线经度,6°带
    //cout<<"L_center="<<L_center<<endl;
    //L_center = L_num * 3;			//中央子午线经度,3°带

    M0 = a * (1 - e);
    M2 = 3.0 / 2.0 * e * M0;
    M4 = 5.0 / 4.0 * e * M2;
    M6 = 7.0 / 6.0 * e * M4;
    M8 = 9.0 / 8.0 * e * M6;

    a0 = M0 + M2 / 2.0 + 3.0 / 8.0 * M4 + 5.0 / 16.0 * M6 + 35.0 / 128.0 * M8;
    a2 = M2 / 2.0 + M4 / 2 + 15.0 / 32.0 * M6 + 7.0 / 16.0 * M8;
    a4 = M4 / 8.0 + 3.0 / 16.0 * M6 + 7.0 / 32.0 * M8;
    a6 = M6 / 32.0 + M8 / 16.0;

    cout << "a0=" << a0 << endl;
    cout << "a2=" << a2 << endl;
    cout << "a4=" << a4 << endl;
    cout << "a6=" << a6 << endl;

    Bf = x / a0;
    B0 = Bf;
    cout<<"B0="<<B0<<endl;
    
     cout<<"sin(2 * B0)="<<sin(2 * B0)/2<<endl;

    while ((fabs(Bf - B0) > 0.0000001) || (B0 == Bf))
    {
        B0 = Bf;
        FBf = -a2 / 2.0 * sin(2 * B0) + a4 / 4.0 * sin(4 * B0) - a6 / 6.0 * sin(6 * B0);
        Bf = (x - FBf) / a0;
    }    //迭代求数值为x坐标的子午线弧长对应的底点纬度

    cout<<"Bf="<<Bf<<endl;

    t = tan(Bf);                            //一样
    c = a * a / b;
    V = sqrt(1 + e1 * cos(Bf) * cos(Bf));   //一样
    N = c / V;                              //一样
    M = c / pow(V, 3);                      //一样
    n = e1 * cos(Bf) * cos(Bf);             //一样(为n的平方)

    n1 = 1 / (N * cos(Bf));
    n2 = -t / (2.0 * M * N);
    n3 = -(1 + 2 * t * t + n) / (6.0 * pow(N, 3) * cos(Bf));
    n4 = t * (5 + 3 * t * t + n - 9 * n * t * t) / (24.0 * M * pow(N, 3));
    n5 = (5 + 28 * t * t + 24 * pow(t, 4) + 6 * n + 8 * n * t * t) / (120.0 * pow(N, 5) * cos(Bf));
    n6 = -t * (61 + 90 * t * t + 45 * pow(t, 4)) / (720.0 * M * pow(N, 5));

    //秒
    double B = (Bf + n2 * y1 * y1 + n4 * pow(y1, 4) + n6 * pow(y1, 6)) / pi * 180;
    
    double L0=l0;

    l = n1 * y1 + n3 * pow(y1, 3) + n5 * pow(y1, 5);
    //double L = L_center + l / pi * 180;    //反算得大地经纬度
    double L = L0 + l / pi * 180;    //反算得大地经纬度

    cout << "方法一 B=" << B << endl;
    cout << "方法一 L=" << L << endl;
}

说明:高斯正算中的输入为度; 

 四、高斯正反算结果

1、已知旧坐标系投影坐标数据

 2、高斯反算(投影坐标转大地坐标)

 3、高斯正算(大地坐标转投影坐标)

 说明:高斯正反算经度由上述转换结果对比可知;

常见坐标系:北京54、西安80、WGS84、CGCS2000等坐标系的高斯转换都做出了实现并使用 QT 进行封装可视化:高斯正反算设计实现!!!-C++文档类资源-CSDN下载

坐标转换整套流程包括:像素坐标转投影坐标、投影坐标转大地坐标、大地坐标转空间直角坐标、七参数转换、空间直角坐标转大地坐标、大地坐标转投影坐标、投影坐标转像素坐标; 本人均已实现,且每一个环节都已经过测试、如有需要欢迎在下方留言评论!!!

如果需要代码包,请在评论区留言!!! 

如果需要代码包,请在评论区留言!!! 

如果需要代码包,请在评论区留言!!! ​​​​​​​

  • 37
    点赞
  • 193
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 48
    评论
高斯/UTM投影是地图制图中常用的投影方式,它通过将地球表面分成若干个椭球体投影带来将地球表面坐标换成直角坐标坐标或将直角坐标坐标换成地球表面坐标的功能。正向换指将地球表面坐标换为直角坐标坐标,反向换则指将直角坐标坐标换为地球表面坐标正反代码包含以下主要步骤: 1. 读取投影带号,确定投影系统和椭球体参数; 2. 对经纬度进行基准面换,将大地坐标换为空间直角坐标; 3. 计投影带宽度,确定中央经线; 4. 计投影坐标系原点; 5. 根据直角坐标坐标和椭球体参数计高斯倍带投影坐标; 6. 根据投影坐标和椭球体参数计经纬度坐标。 正向换的代码示例如下: ```python import math def LLtoUTM(lat, lon): # 读取UTM带号以及椭球体参数 zone = math.floor((lon + 180) / 6) + 1 a = 6378137 f = 1/298.257223563 k0 = 0.9996 # 基准面换 e2 = f * (2 - f) n = f / (2 - f) rho = a * (1 - e2) / pow(1 - e2 * pow(math.sin(lat), 2), 1.5) nu = a / pow(1 - e2 * pow(math.sin(lat), 2), 0.5) psi = nu / rho sin_t = math.sin(lat) cos_t = math.cos(lat) eta2 = e2 * pow(cos_t, 2) # 计投影带宽度 lon0 = (zone - 1) * 6 - 180 + 3 x0 = 500000 y0 = 0 # 计投影坐标系原点 m = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * lat - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * math.sin(2 * lat) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * math.sin(4 * lat) - (35 * e2 * e2 * e2 / 3072) * math.sin(6 * lat)) y = k0 * (m - y0) + k0 * nu * sin_t * ((lon - lon0) * math.pi / 180 + (1 - eta2) * math.sin(2 * lat) / 2 + (5 - 18 * eta2 + eta2 * eta2 + 72 * eta2 * eta2 * eta2) * math.sin(4 * lat) / 24 + (61 - 58 * eta2 + eta2 * eta2 + 600 * eta2 * eta2 * eta2) * math.sin(6 * lat) / 720) x = k0 * nu * cos_t * ((lon - lon0) * math.pi / 180 + (math.sin(lat) / cos_t) * (eta2 / 2 + (5 - eta2 + 9 * eta2 * eta2 + 4 * eta2 * eta2 * eta2) * pow(math.sin(lat), 2) / 24 + (61 - 58 * eta2 + eta2 * eta2 + 600 * eta2 * eta2 * eta2) * pow(math.sin(lat), 4) / 720)) return x, y ``` 反向换的代码示例如下: ```python def UTMtoLL(x, y): # 读取UTM带号以及椭球体参数 a = 6378137 f = 1/298.257223563 k0 = 0.9996 e2 = f * (2 - f) n = f / (2 - f) rho = a * (1 - e2) / pow(1 - e2 * pow(math.sin(lat), 2), 1.5) nu = a / pow(1 - e2 * pow(math.sin(lat), 2), 0.5) psi = nu / rho # 计投影带宽度 w = y / k0 z = (w / a / psi + 5 * (1 - n + 9 * n * n / 4 - 61 * n * n * n / 64) * math.sin(2 * lat) / 16 + (61 * n * n / 32 - 90 * n * n * n / 64) * math.sin(4 * lat) / 32 + (495 * n * n * n / 512) * math.sin(6 * lat)) / math.cos(lat) lat = (w / a - psi * z * math.sin(lat) * (1 + z * z / 6 + (5 - 18 * z * z + z * z * z * z + 72 * z * z * z * z * z *z) / 120)) * 180 / math.pi lon = (zone - 1) * 6 - 180 + 3 + z * math.fmod(1, 400000) / (k0 * nu * math.cos(lat) * math.pi / 180) return lat, lon ``` 以上是高斯/UTM投影坐标正反代码的简单示例。需要注意的是,以上代码并不是完整的实现代码,仅供参考,实际应用仍需要根据具体情况进行适当的修改和调整。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 48
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

数据库内核

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值