wgs84转cgcs2000 C++实现

最新推荐文章于 2025-06-04 18:16:15 发布

原创最新推荐文章于 2025-06-04 18:16:15 发布 · 1.6k 阅读

15 ·

CC 4.0 BY-SA版权

文章标签：

#c++ #开发语言

C++ 专栏收录该内容

18 篇文章

订阅专栏

本文介绍了一个名为GaussProDirect的函数，用于根据经纬度计算在WGS84坐标系统下的高斯投影坐标，作者还提供了自己的实现代码，以及使用示例和精度评估。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

网上很多都不能用，自己实现个

#include <stdio.h>
//#include <windows.h>
#include <iostream>
#include <string>
#include <vector>
#include <stack>
#include <map>
#include <iomanip>
using namespace std;

double Merdian(double EthR0, double Ethe12, double Lat)
{
    // REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.
    double S0 = EthR0 * (1 - Ethe12);
    double e2 = Ethe12;
    double e4 = e2 * e2;
    double e6 = e4 * e2;
    double e8 = e6 * e2;
    double e10 = e8 * e2;
    double e12 = e10 * e2;
    double A1 = 1 + 3 * e2 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e12 / 1048576;
    double B1 = 3 * e2 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e12 / 524288;
    double C1 = 15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e12 / 8388608;
    double D1 = 35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e12 / 1048576;
    double E1 = 315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e12 / 8388608;
    double F1 = 693 * e10 / 1310720 + 9009 * e12 / 5242880;
    double G1 = 1001 * e12 / 8388608;
    double X0 = S0 * (A1 * Lat - B1 * sin(2 * Lat) + C1 * sin(4 * Lat) - D1 * sin(6 * Lat) + 
        E1 * sin(8 * Lat) - F1 * sin(10 * Lat) + G1 * sin(12 * Lat));
    return X0;
}


void GaussProDirect(double longitude, double latitude, double* X, double* Y)
{
    // INPUT // Units of longitude and latitude is RAD (°)
    // REF[1]// 程鹏飞,等.2000国家大地坐标系实用宝典[M].北京:测绘出版社,2008.
    // REF[2]// 孔祥元,郭际明.大地测量学基础(第二版)[M].武汉:武汉大学出版社,2010.
    // Longitude of the Earth central meridian
    double MerLon = 111;
    // Extract the longitudeand latitude
    //Pos = Coord;
    double EthD2R = 0.0174532925199433;// pi / 180
    double Lat = latitude * EthD2R;
    double Lon = longitude * EthD2R - MerLon * EthD2R;
    // Earth orientation parameters of WGS 84 Coordinate System
    double EthR0 = 6378137.0;
    double Ethf = 1 / 298.257223563;//298.257222101; //cgcs2000大地坐标系
    double EthRp = 6356752.3142452;// R0* (1 - f);
    // First Eccentricityand its Squared
    double Ethe12 = 0.006694379990141; // (2f - f * f);
    double Ethe11 = 0.081819190842622;// sqrt(2f - f * f);
    // Second Eccentricityand its Squared
    double Ethe22 = 0.00673949674227643; // (2f - f * f) / (1 + f * f - 2 * f);
    double Ethe21 = 0.08209443794969570;// sqrt((2f - f * f) / (1 + f * f - 2 * f));
    // 高斯投影正算公式
    double RN = EthR0 / sqrt(1 - Ethe12 * sin(Lat) * sin(Lat));
    double Lon2 = Lon * Lon; 
    double Lon4 = Lon2 * Lon2;
    double tnLat = tan(Lat); 
    double tn2Lat = tnLat * tnLat; 
    double tn4Lat = tn2Lat * tn2Lat;
    double csLat = cos(Lat); 
    double cs2Lat = csLat * csLat; 
    double cs4Lat = cs2Lat * cs2Lat;
    double Eta2 = Ethe22 * cs2Lat; 
    double Eta4 = Eta2 * Eta2;
    double NTBLP = RN * tnLat * cs2Lat * Lon2;
    double COE1 = (5 - tn2Lat + 9 * Eta2 + 4 * Eta4) * cs2Lat * Lon2 / 24;
    double COE2 = (61 - 58 * tn2Lat + tn4Lat) * cs4Lat * Lon4 / 720;
    double x = Merdian(EthR0, Ethe12, Lat) + NTBLP * (0.5 + COE1 + COE2);
    double NBLP = RN * csLat * Lon;
    double COE3 = (1 - tn2Lat + Eta2) * cs2Lat * Lon2 / 6;
    double COE4 = (5 - 18 * tn2Lat + tn4Lat + 14 * Eta2 - 58 * tn2Lat * Eta2) * cs4Lat * Lon4 / 120;
    double y = NBLP * (1 + COE3 + COE4) + 500000;
    X = &x; 
    Y = &y;
    printf("GaussProDirect   %lf   %lf\r\n", x, y);

}

int main() {
    //31:49:29   117:55:04
    double longitude = 117.917778;
    double latitude = 31.824722;
    double* X = new double,* Y= new double;
    GaussProDirect(longitude, latitude, X, Y);
    return 0;
}