参考:
https://www.cnblogs.com/rainbow70626/p/7989907.html
https://www.cnblogs.com/DHUtoBUAA/p/6706642.html
http://cdn.hujiulong.com/geohey/blog/mercator/play.html
方法一(不含类):
mercator.h文件:
#pragma once // Add this line before including math.h: #define _USE_MATH_DEFINES // Additions for MS Windows compilation: #ifndef M_PI #define M_PI acos(-1.0) #endif #ifndef M_PI_2 #define M_PI_2 1.57079632679489661922 #endif #include <math.h> /* * Mercator transformation * accounts for the fact that the earth is not a sphere, but a spheroid */ #define D_R (M_PI / 180.0) #define R_D (180.0 / M_PI) #define R_MAJOR 6378137.0 #define R_MINOR 6356752.3142 #define RATIO (R_MINOR/R_MAJOR) #define ECCENT (sqrt(1.0 - (RATIO * RATIO))) #define COM (0.5 * ECCENT) inline double fmin(double x, double y) { return(x < y ? x : y); } inline double fmax(double x, double y) { return(x > y ? x : y); } static double deg_rad(double ang) { return ang * D_R; } double merc_x(double lon) { return R_MAJOR * deg_rad(lon);//墨卡托投影到平面的横坐标 } double merc_y(double lat) { lat = fmin(89.5, fmax(lat, -89.5)); double phi = deg_rad(lat); double sinphi = sin(phi); double con = ECCENT * sinphi; con = pow((1.0 - con) / (1.0 + con), COM); double ts = tan(0.5 * (M_PI * 0.5 - phi)) / con; return 0 - R_MAJOR * log(ts); } static double rad_deg(double ang) { return ang * R_D; } double merc_lon(double x) { return rad_deg(x) / R_MAJOR;//返回经度 } double merc_lat(double y) { double ts = exp(-y / R_MAJOR); double phi = M_PI_2 - 2 * atan(ts); double dphi = 1.0; int i; for (i = 0; fabs(dphi) > 0.000000001 && i < 15; i++) { double con = ECCENT * sin(phi); dphi = M_PI_2 - 2 * atan(ts * pow((1.0 - con) / (1.0 + con), COM)) - phi; phi += dphi; } return rad_deg(phi);//返回纬度 }
mercator.cpp文件:
/* 本文件来自于:https://wiki.openstreetmap.org/wiki/Mercator */ #include"mercator.h" #include<iostream> using namespace std; int main() { cout <<"输出横坐标:"<< merc_x(120) << endl; cout << "输出纵坐标:" << merc_y(60) << endl; cout << "输出经度:" << merc_lon(654321) << endl; cout << "输出纬度:" << merc_lat(123456) << endl; system("pause"); return 0; }
方法二(用类写):
mercator.h文件:
#pragma once #include<cmath> //圆周率 const double PI = 3.1415926535897932; //角度到弧度的转换 double DegreeToRad(double degree) { return PI*((double)degree / (double)180); } //弧度到角度的转换 double RadToDegree(double rad) { return (180 * rad) / PI; } class MercatorProj { public: MercatorProj(); void SetAB(double& a, double& b); void SetB0(double b0); void SetL0(double l0); int FromProj(double X, double Y, double& B, double& L); int ToProj(double B, double L, double &X, double &Y); ~MercatorProj(); private: int __IterativeTimes; //反向转换程序中的迭代次数 double __IterativeValue; //反向转换程序中的迭代初始值 double __A; //椭球体长半轴,米 double __B; //椭球体短半轴,米 double __B0; //标准纬度,弧度 double __L0; //原点经度,弧度 }; //构造函数中赋予参数默认值 MercatorProj::MercatorProj() :__IterativeTimes(15), __IterativeValue(0), __A(0), __B(0),__B0(0), __L0(0) { } //设定__A与__B void MercatorProj::SetAB(double& a, double& b)//原程序的参数写作double a, double b { if (a <= 0 || b <= 0) { return; } __A = a; __B = b; } //设定__B0 void MercatorProj::SetB0(double b0) { if (b0<-PI / 2 || b0>PI / 2) { return; } __B0 = b0;//设定标准纬度,标准纬度线和原点经度线的交点组成了投影后的平面坐标的原点 } //设定__L0 void MercatorProj::SetL0(double l0) { if (l0<-PI || l0>PI) { return; } __L0 = l0;//设定原点经度 } /******************************************* 投影正向转换程序 double B: 维度,弧度 double L: 经度,弧度 double& X: 纵向直角坐标 double& Y: 横向直角坐标 *******************************************/ int MercatorProj::ToProj(double B, double L, double &X, double &Y) { double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp; double E = exp(1); if (L<-PI || L>PI || B<-PI / 2 || B>PI / 2) { return 1; } if (__A <= 0 || __B <= 0) { return 1; } f = (__A - __B) / __A; dtemp = 1 - (__B / __A)*(__B / __A); if (dtemp<0) { return 1; } e = sqrt(dtemp); dtemp = (__A / __B)*(__A / __B) - 1; if (dtemp<0) { return 1; } e_ = sqrt(dtemp); NB0 = ((__A*__A) / __B) / sqrt(1 + e_*e_*cos(__B0)*cos(__B0)); //NB0=N, 法线长(《大地测量学基础》第二版第103页,或109、110页) K = NB0*cos(__B0);//平行圈半径 Y = K*(L - __L0); X = K*log(tan(PI / 4 + B / 2)*pow((1 - e*sin(B)) / (1 + e*sin(B)), e / 2)); return 0; } /******************************************* 投影反向转换程序 double X: 纵向直角坐标 double Y: 横向直角坐标 double& B: 维度,弧度 double& L: 经度,弧度 *******************************************/ int MercatorProj::FromProj(double X, double Y, double& B, double& L) { double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp; double E = exp(1); if (__A <= 0 || __B <= 0) { return 1; } f = (__A - __B) / __A; dtemp = 1 - (__B / __A)*(__B / __A); if (dtemp<0) { return 1; } e = sqrt(dtemp); dtemp = (__A / __B)*(__A / __B) - 1; if (dtemp<0) { return 1; } e_ = sqrt(dtemp); NB0 = ((__A*__A) / __B) / sqrt(1 + e_*e_*cos(__B0)*cos(__B0)); K = NB0*cos(__B0); L = Y / K + __L0; B = __IterativeValue; for (int i = 0; i<__IterativeTimes; i++) { B = PI / 2 - 2 * atan(pow(E, (-X / K))*pow(E, (e / 2)*log((1 - e*sin(B)) / (1 + e*sin(B))))); } return 0; } MercatorProj::~MercatorProj() { }
mercator.cpp文件:
#include"mercator_projection.h" #include<iostream> using namespace std; //测试主函数: int main() { MercatorProj m_mp; double b0, l0; double latS, lgtS, latD, lgtD; b0 = 30; //b0 = 0; l0 = 0; latS = 60;//测试数据,纬度60度 lgtS = 120;//测试数据,经度120度 double a = 6378137;//长半轴 double b = 6356752.3142;//短半轴 m_mp.SetAB(a, b); // WGS 84 m_mp.SetB0(DegreeToRad(b0)); m_mp.SetL0(DegreeToRad(l0)); m_mp.ToProj(DegreeToRad(latS), DegreeToRad(lgtS), latD, lgtD); cout << "X="<<latD << " " <<"Y="<< lgtD << endl; // 7248377.351067:11578353.630128 latS = 123456;//测试数据 lgtS = 654321;//测试数据 m_mp.FromProj(latS, lgtS, latD, lgtD); latD = RadToDegree(latD); lgtD = RadToDegree(lgtD); cout << "B="<<latD << " " <<"L="<< lgtD << endl; // 1.288032: 6.781493 system("pause"); return 0; }