高斯平面坐标与大地坐标转换

最新推荐文章于 2025-04-30 15:24:54 发布

小白菜鸡_

最新推荐文章于 2025-04-30 15:24:54 发布

阅读量1.3k

点赞数 16

CC 4.0 BY-SA版权

文章标签： c++

本文链接：https://blog.youkuaiyun.com/iioo_ooo_ioo/article/details/142151068

文章目录

包含两部分：
- 高斯平面坐标xy转为大地坐标BL。
- 大地坐标BL转为高斯平面坐标xy。
用C++语言编写，能够读取txt坐标文件，并将结果导出到新的txt文件。
椭球模型为WGS84，可以根据实际需要更改参数。
不足之处请指正！

张华海.应用大地测量学(第4版).中国矿业大学出版社.2016

1. 高斯平面坐标xy转为大地坐标BL

convert.h

#pragma once
#include<math.h>
#include<iomanip>
#include <cmath>
using namespace std;

#define PI  3.14159265358979323846;

typedef struct
{
    double pos1;
    double pos2;
}position;

//WGS84椭球参数
double a = 6378137.0;
double f = 1 / 298.257223563;

double e2 = 2 * f - f * f;   //第一偏心率平方
double e_2 = e2 * (1 + e2);  //第二偏心率平方
double c = a / sqrt(1 - e2);
double beita0 = 1 - 3 * e_2 / 4 + 45 * e_2 * e_2 / 64 - 175 * e_2 * e_2 * e_2 / 256 + 11025 * e_2 * e_2 * e_2 * e_2 / 16384;
double beita2 = beita0 - 1;
double beita4 = 15 * e_2 * e_2 / 32 - 175 * e_2 * e_2 * e_2 / 384 + 3675 * e_2 * e_2 * e_2 * e_2 / 8192;
double beita6 = -35 * e_2 * e_2 * e_2 / 96 + 735 * e_2 * e_2 * e_2 * e_2 / 2048;
double beita8 = 315 * e_2 * e_2 * e_2 * e_2 / 1024;

double FxB(double B)
{
    double temp = (c * beita4 + (c * beita6 + c * beita8 * cos(B) * cos(B)) * cos(B) * cos(B)) * cos(B) * cos(B);
    return c * beita2 + temp * sin(B) * cos(B);
}

double FxBl(double B, double l)
{
    double W = sqrt(1 - e2 * sin(B) * sin(B));
    double N = a / W;  //卯酉圈曲率半径
    double t2 = tan(B) * tan(B);
    double yita2 = e_2 * cos(B) * cos(B);

    double a2 = N * sin(B) * cos(B) / 2;
    double a4 = N * sin(B) * cos(B) * cos(B) * cos(B) * (5 - t2 + 9 * yita2 + 4 * yita2 * yita2) / 24;
    double a6 = N * sin(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (61 - 58 * t2 + t2 * t2) / 720;

    return a2 * l * l + a4 * l * l * l * l + a6 * l * l * l * l * l * l;
}

double FyBl(double B, double l)
{
    double W = sqrt(1 - e2 * sin(B) * sin(B));
    double N = a / W;  //卯酉圈曲率半径
    double t2 = tan(B) * tan(B);
    double yita2 = e_2 * cos(B) * cos(B);

    double a3 = N * cos(B) * cos(B) * cos(B) * (1 - t2 + yita2) / 6;
    double a5 = N * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (5 - 18 * t2 + t2 * t2 + 14 * yita2 - 58 * t2 * yita2) / 120;

    return a3 * l * l * l + a5 * l * l * l * l * l;
}



//坐标反算。高斯平面坐标 x、y -> 经纬度 B(latitude)、L(longitude)
//y = N * 1000000 + 500000 + y
position xy2BL(double x, double y)
{
    //带宽，可选6或者3
    int zoneWide = 6;

    //带号
    int zoneNum = 0;
    int temp = (int)floor(log(abs((double)y)) / log((double)10)) + 1;
    zoneNum = (int)(y / pow((double)10, temp - 2));

    //中央子午线经度
    int L0 = 0;
    if (zoneWide == 6)
        L0 = 6 * zoneNum - 3;
    else
        L0 = 3 * zoneNum;


    //迭代初值
    y = y - (double)zoneNum * 1000000 - 500000;
    double B = x / (c * beita0);
    double a1 = a * cos(B) / sqrt(1 - e2 * sin(B) * sin(B));
    double l = y / a1;

    //迭代
    double R2D = 180 / PI;
    double tempB = 0;
    double templ = 0;
    while (abs(B - tempB) * R2D > 0.00001 || abs(l - templ) * R2D > 0.00001)
    {
        tempB = B;
        templ = l;

        B = (x - FxB(B) - FxBl(B, l)) / (c * beita0);
        a1 = a * cos(B) / sqrt(1 - e2 * sin(B) * sin(B));
        l = (y - FyBl(B, l)) / a1;
    }

    position BL;
    BL.pos1 = B*R2D;
    BL.pos2 = l*R2D + L0;
    return BL;
}

main.cpp

#include<iostream>
#include<fstream>
#include<string>
#include <vector> 
#include "convert.h"

int rolnum = 32420;
double xy_[32420][2];
double BL_[32420][2];

int main()
{
    double xx = 0.0, yy = 0.0;

    ifstream infile;    //输入流
    ofstream outfile;   //输出流
    infile.open("E:/xy.txt", ios::in);  //打开文件

    //读取文件直到末尾，数据存在 BLH 中
    int i = 0;
    while (!infile.eof())
    {
        infile >> xx >> yy;
        xy_[i][0] = xx;
        xy_[i][1] = yy;
        i++;
    }
    infile.close();   //关闭文件

    int k = 0;
    position temp;
    for (k = 0; k < rolnum; k++)
    {
        temp = xy2BL(xy_[k][0], xy_[k][1]);
        BL_[k][0] = temp.pos1;
        BL_[k][1] = temp.pos2;
    }

    //每次写都定位的文件结尾，不会丢失原来的内容用app，用out则会丢失原来的内容
    int j = 0;
    outfile.open("E:/BL.txt", ios::out);
    for (j = 0; j < rolnum; ++j)
        outfile << setprecision(11) << BL_[j][0] << "\t" << setprecision(12) << BL_[j][1] << endl;
    outfile.close();

    cout << "**************** finished！****************" << endl;

	return 0;
}

2. 大地坐标BL转为高斯平面坐标xy

convert.h

#include<math.h>
#include<iomanip>
using namespace std;

typedef struct
{
    double pos1;
    double pos2;
}position;

//WGS84椭球参数
double a = 6378137.0;
double f = 1 / 298.257223563;

//坐标正算。经纬度 B(latitude)、L(longitude) -> 高斯平面坐标 x、y
position BL2xy(double latitude, double longitude)
{
    //带宽，可选6或者3
    int zoneWide = 6;

    //带号
    int zoneNum;
    if (zoneWide == 6)
        zoneNum = (int)longitude / 6 + 1;
    else
        zoneNum = (int)longitude / 3;

    //中央子午线经度
    int L0;
    if (zoneWide == 6)
        L0 = 6 * zoneNum - 3;
    else
        L0 = 3 * zoneNum;

    //经差
    double l = longitude - L0;
    
    //角度->弧度
    double D2R = 3.1415926535898 / 180.0;
    latitude = latitude * D2R;
    longitude = longitude * D2R;
    l = l * D2R;

    double e2 = 2 * f - f * f;   //第一偏心率平方
    double e_2 = e2 * (1 + e2);  //第二偏心率平方
    double W = sqrt(1 - e2 * sin(latitude) * sin(latitude));
    double N = a / W;  //卯酉圈曲率半径
    double t2 = tan(latitude) * tan(latitude);
    double yita2 = e_2 * cos(latitude) * cos(latitude);

    double c = a / sqrt(1 - e2);
    double beita0 = 1 - 3 * e_2 / 4 + 45 * e_2 * e_2 / 64 - 175 * e_2 * e_2 * e_2 / 256 + 11025 * e_2 * e_2 * e_2 * e_2 / 16384;
    double beita2 = beita0 - 1;
    double beita4 = 15 * e_2 * e_2 / 32 - 175 * e_2 * e_2 * e_2 / 384 + 3675 * e_2 * e_2 * e_2 * e_2 / 8192;
    double beita6 = -35 * e_2 * e_2 * e_2 / 96 + 735 * e_2 * e_2 * e_2 * e_2 / 2048;
    double beita8 = 315 * e_2 * e_2 * e_2 * e_2 / 1024;
    double X = c * (beita0 * latitude + sin(latitude) * (beita2 * cos(latitude)
        + beita4 * cos(latitude) * cos(latitude) * cos(latitude)
        + beita6 * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude)
        + beita8 * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude)));

    double x, y;
    x = X 
        + l * l * N * sin(latitude) * cos(latitude) / 2 
        + l * l * l * l * N * sin(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * (5 - t2 + 9 * yita2 + 4 * yita2 * yita2) / 24
        + l * l * l * l * l * l * N * sin(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * (61 - 58 * t2 + t2 * t2) / 720;

    y = N * cos(latitude) * l
        + (N * cos(latitude) * cos(latitude) * cos(latitude) * (1 - t2 + yita2) * l * l * l) / 6
        + (N * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * (5 - 18 * t2 + t2 * t2 + 14 * yita2 - 58 * t2 * yita2) * l * l * l * l * l) / 120
        + (N * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * cos(latitude) * (61 - 479 * t2 + 179 * t2 * t2 - t2 * t2 * t2) * l * l * l * l * l * l * l) / 5040;

    position xy;
    xy.pos1 = x;
    xy.pos2 = zoneNum * 1000000 + 500000 + y;
    return xy;
}

main.cpp

#include<iostream>
#include<fstream>
#include<string>
#include <vector> 
#include "convert.h"

int rolnum = 32420;
double BLH[32420][3];
double xy_[32420][2];

int main() 
{
    double bb = 0.0, ll = 0.0, hh = 0.0;
    
    ifstream infile;    //输入流
    ofstream outfile;   //输出流
    infile.open("E:/BLH.txt", ios::in);  //打开文件

    //读取文件直到末尾，数据存在 BLH 中
    int i = 0;
    while (!infile.eof())
    {
        infile >> bb >> ll >> hh;
        BLH[i][0] = bb;
        BLH[i][1] = ll;
        BLH[i][2] = hh;
        i++;
    }
    infile.close();   //关闭文件

    //计算 x、y，存入 xy_z 中
    int k = 0;
    position temp;
    for (k = 0; k < rolnum; k++)
    {
        temp = BL2xy(BLH[k][0], BLH[k][1]);
        xy_[k][0] = temp.pos1;
        xy_[k][1] = temp.pos2;
    }

    //每次写都定位的文件结尾，不会丢失原来的内容用app，用out则会丢失原来的内容
    int j = 0;
    outfile.open("E:/xy.txt", ios::out);
    for (j = 0; j < rolnum; ++j)
        outfile << setprecision(13) << xy_[j][0] << "\t" << setprecision(14) << xy_[j][1] << endl;
    outfile.close();

    cout <<"**************** finished！****************" << endl;
    return 0;
}