经纬度与高斯-克吕格平面坐标转换

原作者链接：https://blog.csdn.net/jianyi7659/article/details/7583339

前言
支持将地理坐标（经纬度坐标）转换到高斯-克吕格投影下的平面坐标，如北京54平面坐标、西安80平面坐标、大地2000平面坐标。

LatLong2GaussKruger.js

/**
 * WGS-84椭球经纬度坐标转高斯-克吕格平面坐标
 */

    
        //基本变量定义
var a = 0;//'椭球体长半轴
var b = 0;// '椭球体短半轴
var f = 0;//'扁率
var e = 0;// '第一偏心率
var e1 = 0; //'第二偏心率

var FE= 0;//'东偏移
var FN=0;//'北偏移
var L0=0;//'中央经度
var W0=0;//'原点纬线
var k0=0;//'比例因子
        // PI: 3.14159265358979,

        /**
         * 用于初始化椭球参数
         * @param {*} TuoqiuCanshu 枚举类型，提供北京54、西安80、WGS84、CGCS2000椭球参数
         * @param {*} CentralMeridian 中央经线
         * @param {*} OriginLatitude 原点纬度，如果是标准的分幅，则该参数是0
         * @param {*} EastOffset 东偏移
         * @param {*} NorthOffset 北偏移
         */
        function constructor(TuoqiuCanshu, CentralMeridian, OriginLatitude, EastOffset, NorthOffset) {
            /*
             *  Canshu
             *  Beijing54 = 0;Krassovsky （北京54采用） 6378245 6356863.0188
                Xian80 = 1;IAG 75（西安80采用） 6378140 6356755.2882
                WGS84 = 2;WGS 84 6378137 6356752.3142
                CGCS2000 = 3
            */
            console.log(TuoqiuCanshu)
           if (TuoqiuCanshu == 0)//北京五四
            {
                a = 6378245;
                b = 6356863.0188;
            }
            if (TuoqiuCanshu == 1)// '西安八零
            {
                a = 6378140;
                b = 6356755.2882;
            }
            if (TuoqiuCanshu == 2)//'WGS84
            {
                a = 6378137;
                b = 6356752.3142;
            }
            if (TuoqiuCanshu == 3)//'CGCS2000
            {
                a = 6378137;
                b = 6356752.314140356;
            }
            f = (a - b) / a;//扁率
            console.log(f)
            //e = Math.sqrt(1 - Math.pow((b / a) ,2));//'第一偏心率
            e = Math.sqrt(2 * f - Math.pow(f, 2));//'第一偏心率

            //eq = Math.sqrt(Math.pow((a / b) , 2) - 1);//'第二偏心率
            e1 = e / Math.sqrt(1 - Math.pow(e, 2));//'第二偏心率
            L0 = CentralMeridian;//中央经
            W0 = OriginLatitude;//原点纬线
            k0 = 1;//'比例因子
            FE = EastOffset;//东偏移
            FN = NorthOffset;//北偏移
        }
        /**
         * 经纬度转高斯-克吕格平面坐标
         * @param {*} J 经度
         * @param {*} W 维度
         */
function JWgetGK(J, W) {
            
            //'给出经纬度坐标，转换为高克投影坐标
            var resultP = {};
            var BR = (W - W0) * Math.PI / 180;//纬度弧长
            var lo = (J - L0) * Math.PI / 180; //经差弧度
             
            var N = a / Math.sqrt(1 - Math.pow((e * Math.sin(BR)), 2)) //卯酉圈曲率半径
          
            //求解参数s
            var B0;
            var B2;
            var B4;
            var B6;
            var B8;
    var C = Math.pow(a, 2) / b;
            B0 = 1 - 3 * Math.pow(e1, 2) / 4 + 45 * Math.pow(e1, 4) / 64 - 175 * Math.pow(e1, 6) / 256 + 11025 * Math.pow(e1, 8) / 16384;
    B2 = B0 - 1
    
            B4 = 15 / 32 * Math.pow(e1, 4) - 175 / 384 * Math.pow(e1, 6) + 3675 / 8192 * Math.pow(e1, 8);
            B6 = 0 - 35 / 96 * Math.pow(e1, 6) + 735 / 2048 * Math.pow(e1, 8);
            B8 = 315 / 1024 * Math.pow(e1, 8);
            var s = C * (B0 * BR + Math.sin(BR) * (B2 * Math.cos(BR) + B4 * Math.pow((Math.cos(BR)), 3) + B6 * Math.pow((Math.cos(BR)), 5) + B8 * Math.pow((Math.cos(BR)), 7)));
          
            var t = Math.tan(BR);
            var g = e1 * Math.cos(BR);

            var XR = s + Math.pow(lo, 2) / 2 * N * Math.sin(BR) * Math.cos(BR) + Math.pow(lo, 4) * N * Math.sin(BR) * Math.pow((Math.cos(BR)), 3) / 24 * (5 - Math.pow(t, 2) + 9 * Math.pow(g, 2) + 4 * Math.pow(g, 4)) + Math.pow(lo, 6) * N * Math.sin(BR) * Math.pow((Math.cos(BR)), 5) * (61 - 58 * Math.pow(t, 2) + Math.pow(t, 4)) / 720;
            var YR = lo * N * Math.cos(BR) + Math.pow(lo, 3) * N / 6 * Math.pow((Math.cos(BR)), 3) * (1 - Math.pow(t, 2) + Math.pow(g, 2)) + Math.pow(lo, 5) * N / 120 * Math.pow((Math.cos(BR)), 5) * (5 - 18 * Math.pow(t, 2) + Math.pow(t, 4) + 14 * Math.pow(g, 2) - 58 * Math.pow(g, 2) * Math.pow(t, 2));

            resultP.x = YR + FE;
    resultP.y = XR + FN;
  
            return resultP;
        }

 /**
         带带号计算方法
         */
function JWgetGK1(longitude, latitude) {
    var resultP = {};
    var ProjNo = 0;
    var ZoneWide; //带宽
    var longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
    var a, f, e2, ee, NN, T, C, A, M, iPI;
    iPI = 0.0174532925199433;  //3.1415926535898 / 180.0;
    ZoneWide = 6; //*度带宽 
    // a = 6378245.0; f = 1.0 / 298.3; //54年北京坐标系参数 
    a = 6378140.0; f = 1 / 298.257; //80年西安坐标系参数 
    ProjNo = parseInt(longitude / ZoneWide);
    longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
    longitude0 = longitude0 * iPI;
    latitude0 = 0;
    longitude1 = longitude * iPI; //经度转换为弧度
    latitude1 = latitude * iPI; //纬度转换为弧度
    e2 = 2 * f - f * f;
    ee = e2 / (1.0 - e2);
    NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1) * Math.sin(latitude1));
    T = Math.tan(latitude1) * Math.tan(latitude1);
    C = ee * Math.cos(latitude1) * Math.cos(latitude1);
    A = (longitude1 - longitude0) * Math.cos(latitude1);
    M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude1));
    xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14 * C - 58 * ee) * A * A * A * A * A / 120);
    yval = M + NN * Math.tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 270 * C - 330 * ee) * A * A * A * A * A * A / 720);
    X0 = 1000000 * (ProjNo + 1) + 500000;
    Y0 = 0;
    xval = xval + X0; yval = yval + Y0;
    resultP.x = xval;
    resultP.y = yval;

    return resultP;
}


        /**
         * 高斯-克吕格平面转经纬度坐标
         * @param {*} X  
         * @param {*} Y 
         */
        function GKgetJW(X, Y) {
            //'给出高克投影坐标，转换为经纬度坐标
            var resultP = {};
            var El1 = (1 - Math.sqrt(1 - Math.pow(e, 2))) / (1 + Math.sqrt(1 - Math.pow(e, 2)));

            var Mf = (Y - FN) / k0;//真实坐标值

            var Q = Mf / (a * (1 - Math.pow(e, 2) / 4 - 3 * Math.pow(e, 4) / 64 - 5 * Math.pow(e, 6) / 256));//角度

            Bf = Q + (3 * El1 / 2 - 27 * Math.pow(El1, 3) / 32) * Math.sin(2 * Q) + (21 * Math.pow(El1, 2) / 16 - 55 * Math.pow(El1, 4) / 32) * Math.sin(4 * Q) + (151 * Math.pow(El1, 3) / 96) * Math.sin(6 * Q) + 1097 / 512 * Math.pow(El1, 4) * Math.sin(8 * Q);
            Rf = a * (1 - Math.pow(e, 2)) / Math.sqrt(Math.pow((1 - Math.pow((e * Math.sin(Bf)), 2)), 3));
            Nf = a / Math.sqrt(1 - Math.pow((e * Math.sin(Bf)), 2));//卯酉圈曲率半径
            Tf = Math.pow((Math.tan(Bf)), 2);
            D = (X - FE) / (k0 * Nf);

            Cf = Math.pow(e1, 2) * Math.pow((Math.cos(Bf)), 2);

            var B = Bf - Nf * Math.tan(Bf) / Rf * (Math.pow(D, 2) / 2 - (5 + 3 * Tf + 10 * Cf - 9 * Tf * Cf - 4 * Math.pow(Cf, 2) - 9 * Math.pow(e1, 2)) * Math.pow(D, 4) / 24 + (61 + 90 * Tf + 45 * Math.pow(Tf, 2) - 256 * Math.pow(e1, 2) - 3 * Math.pow(Cf, 2)) * Math.pow(D, 6) / 720);
            var L = CentralMeridian * Math.PI / 180 + 1 / Math.cos(Bf) * (D - (1 + 2 * Tf + Cf) * Math.pow(D, 3) / 6 + (5 - 2 * Cf + 28 * Tf - 3 * Math.pow(Cf, 2) + 8 * Math.pow(e1, 2) + 24 * Math.pow(Tf, 2)) * Math.pow(D, 5) / 120);

            var Bangle = B * 180 / Math.PI;
            var Langle = L * 180 / Math.PI;

            resultP.x = Langle;//RW * 180 / Math.PI;
            resultP.y = Bangle + W0;//RJ * 180 / Math.PI;
            return resultP;
        }
   


export default {constructor,JWgetGK};

使用示例

import LatLong2GaussKruger from "../Main/earthContainer/LatLong2GaussKruger"


       LatLong2GaussKruger.constructor(2,120,0,500000,0)
        console.log("经纬度坐标：120.335194,29.313418")
        let coord = LatLong2GaussKruger.JWgetGK(120.335194,29.313418);
        console.log("转换后：" + coord.x + "," + coord.y);
![效果](https://img-blog.csdnimg.cn/20201214223038100.png#pic_center)