墨卡托投影代码

转载

mb5ff9812a47cf3 2019-01-19 09:24:00

文章标签 #define #include 测试数据迭代 ios 文章分类 代码人生

方法一（不含类）：

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 << "\t" <<"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 << "\t" <<"L="<< lgtD << endl;
    // 1.288032: 6.781493
    system("pause");
    return 0;
}