大地坐标系简介

大地坐标系是大地测量中以参考椭球面为基准面建立起来的坐标系，也就是地理坐标系，地面点的位置用大地经度L、大地纬度B和大地高度H表示，单位为度分秒。其中有三类比较常用的大地坐标系统，即地心坐标系统、参心坐标系统和地方独立坐标系统。

介绍不同坐标系统之前，首先简要介绍下经纬度的概念。

大地经纬度示意图：（B为纬度，L为经度，H为高度）

某地的大地经度是本初子午线与本地子午线的夹角，如图中的经度。本初子午线是指通过英国伦敦格林尼治天文台原址的那条经线，也称为0°经线。大地经度分为东经和西经，东经为正，西经为负，东西各为180度，每15度时间相差一个小时。

某地的大地纬度是将地球表面点投射到地球椭球体上，投影点法线方向与赤道平面的夹角为该点的大地纬度，即垂直于该地椭球面的线与赤道面的夹角。

1、地心坐标系 地心坐标系是以地球质心为原点建立的空间直角坐标系，或以球心与地球质心重合的地球椭球面为基准面所建立的大地坐标系，通常分为地心空间直角坐标系(以x，y，z为其坐标元素)和地心大地坐标系(以B，L，H为其坐标元素)。

常用的有WGS84、CGCS2000等。（主要是选的椭球体不一样，所以有不同的坐标系统）

其中WGS84是美国国家影像制图局NIMA在WGS72基础上改进的、全球统一的地心坐标系，是现有应用于导航、精确大地测量和地图制图的最好的全球大地参考系（GPS星历的坐标系，于2002年修正）。

注：地心坐标系又名协议地球坐标系，与GPS的瞬时地球坐标系要对应起来。

2、参心坐标系 WGS是针对全球建立的坐标系，针对某一地区可能误差较大，为了让地图能尽可能描述某个地区的地形，便建立所谓的“参心坐标系”。为了贴合某地区，参心坐标系的中心并不在地球质心，是以所选椭球体圆心为中心，X、Y、Z三轴的定义与地心坐标系一致。

我国常用的参心坐标系有BJZ54（原）、西安80和BJZ54（北京54）（新）等。

3、常见坐标系简介 1、WGS坐标系 WGS84坐标系全称World Geodetic System - 1984，是为了解决GPS定位而产生的全球统一的一个地心坐标系。

● 椭球体：WGS84椭球 ● 长半径=6378137m ● 短半径b=6356752.3142m ● 第一偏心率平方e2=0.00669437999013 ● 第二偏心率平方e2=0.006739496742227 ● 扁率=1/298.257223563

2、BJZ54坐标系 北京54参心坐标系，它是以克拉索夫斯基椭球为基础，经局部平差后产生的坐标系。1954年北京坐标系可以认为是前苏联1942年坐标系的延伸。它的原点不在北京而是在前苏联的普尔科沃。

● 椭球体：克拉索夫斯基椭球 ● 长半径=6378245m ● 短半径b=6356863.0188m ● 第一偏心率平方e2=0.006693421622 ● 扁率=1/298.3

3、西安80坐标系 西安80坐标系是指1980年西安参心坐标系，该坐标系的大地原点设在我国中部的陕西省泾阳县永乐镇，位于西安市西北方向约60公里，又简称西安大地原点。基准面为黄海水面。 ● 椭球体：IAG椭球 ● 长半径 a=6378140m ● 短半径 b=6356755m ● 第一偏心率平方e2=0.00669438499959 ● 扁率=1/298.25722101

4、CGCS2000坐标系 2000国家大地坐标系是全球地心坐标系在我国的具体体现，是我国当前最新的国家大地坐标系，其全称为China Geodetic Coordinate System 2000，其原点为包括海洋和大气的整个地球的质量中心。

● 椭球体：CGCS2000坐标系 ● 长半轴 a=6378137m ● 短半轴 b=6356752.314m ● 第一偏心率平方e2=0.00669438002290 ● 扁率=1/298.257222101

js大地坐标与经纬度坐标互转

//高斯投影坐标反算成经纬度
function GaussToBL(X,Y){
    let ProjNo;
    let ZoneWide; 带宽
    let output = new Array(2);
    let longitude1,latitude1, longitude0, X0,Y0, xval,yval;//latitude0,
    let e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
    iPI = 3.14159265358979324/180.0; 3.1415926535898/180.0;
    // a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数
    a=6378140.0; f=1.0/298.257; //80年西安坐标系参数
    ZoneWide = 6; 6度带宽
    ProjNo = parseInt(X/1000000) ; //查找带号
    longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;
    longitude0 = longitude0 * iPI ; //中央经线


    X0 = ProjNo*1000000+500000;
    Y0 = 0;
    xval = X-X0; yval = Y-Y0; //带内大地坐标
    e2 = 2*f-f*f;
    e1 = (1.0-Math.sqrt(1-e2))/(1.0+Math.sqrt(1-e2));
    ee = e2/(1-e2);
    M = yval;
    u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
    fai = u+(3*e1/2-27*e1*e1*e1/32)*Math.sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(4*u) +(151*e1*e1*e1/96)*
        Math.sin(6*u)+(1097*e1*e1*e1*e1/512)*Math.sin(8*u);
    C = ee*Math.cos(fai)*Math.cos(fai);
    T = Math.tan(fai)*Math.tan(fai);
    NN = a/Math.sqrt(1.0-e2*Math.sin(fai)*Math.sin(fai));
    R = a*(1-e2)/Math.sqrt((1-e2*Math.sin(fai)*Math.sin(fai))*(1-e2*Math.sin(fai)*Math.sin(fai))*(1-e2*Math.sin
    (fai)*Math.sin(fai)));
    D = xval/NN;
    //计算经度(Longitude) 纬度(Latitude)
    longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D
        *D*D*D*D/120)/Math.cos(fai);
    latitude1 = fai -(NN*Math.tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24
        +(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);
    //转换为度 DD
    output[0] = longitude1 / iPI;
    output[1] = latitude1 / iPI;
    return output;
}

//经纬度=>高斯投影
function BLToGauss(longitude, latitude){
    let ProjNo=0;
    let ZoneWide; 带宽
    let ret=Array(2);
    let longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
    let a,f, e2,ee, NN, T,C,A, M, iPI;
    iPI = 0.0174532925199433; 3.1415926535898/180.0;
    ZoneWide = 6; 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+72*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+600*C-330*ee)*A*A*A*A*A*A/720);
    X0 = 1000000*(ProjNo+1)+500000;
    Y0 = 0;
    xval = xval+X0; yval = yval+Y0;
    ret[0]=xval;
    ret[1]=yval;

    return ret;
}

export default {
    GaussToBL,
    BLToGauss,
}