• 墨卡托投影代码


    参考:

    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;
    }
  • 相关阅读:
    阿里云短信服务工具类
    vue.config.js
    elementui Tree 树形控件增删改查
    vue 实时显示年月日时分秒星期上下午
    1553:【例 2】暗的连锁
    CF825G Tree Queries
    最短母串
    寻找好串
    无限链计数
    异或运算
  • 原文地址:https://www.cnblogs.com/yibeimingyue/p/10290485.html
Copyright © 2020-2023  润新知