坐标系转换算法与投影算法

火星坐标

火星坐标是国家测绘局为了国家安全在原始坐标的基础上进行偏移得到的坐标，基本国内的电子地图、导航设备都是采用的这一坐标系或在这一坐标的基础上进行二次加密得到的。 火星坐标的真实名称应该是 GCJ-02 坐标。基本上所有的国内的电子地图采用的都是火星坐标系甚至 Google 地图中国部分都特意为中国政府做了偏移。

百度坐标

火星坐标是在国际标准坐标 WGS-84 上进行的一次加密，由于国内的电子地图都要至少使用火星坐标进行一次加密，百度直接就任性一些，直接自己又研究了一套加密算法，来了个 二次加密，这就是我们所熟知的百度坐标(BD-09)。

WGS-84 坐标

WGS-84 坐标是一个国际的标准，一般卫星导航，原始的 GPS 设备中的数据都是采用这一坐标系。 国外的 Google 地图、OSM 等采用的都是这一坐标。

Node.js 库推荐

一个提供了百度坐标（BD09）、国测局坐标（火星坐标，GCJ02）、和 WGS84 坐标系之间的转换的工具模块。

coordtransform

墨卡托投影

墨卡托投影（Mercator Projection）— 等角正轴切圆柱投影

由荷兰地图学家墨卡托（G. Mercator）于 1569 年创拟，为地图投影方法中影响最大的，又称正轴等角圆柱投影。

假设地球被围在一中空的圆柱里，其基准纬线与圆柱相切（赤道）接触，然后再假想地球中心有一盏灯，把球面上的图形投影到圆柱体上，再把圆柱体展开，这就是一幅选定基准纬线上的“墨卡托投影”绘制出的地图。

墨卡托投影特点：

1、没有角度变形，由每一点向各方向的长度比相等；

2、经纬线都是平行直线，且相交成直角；经线间隔相等，纬线间隔从基准纬线处向两极逐渐增大。

3、长度和面积变形明显，但基准纬线处无变形，变形从基准纬线处向两极变形逐渐增大，但因为它具有各个方向均等扩大的特性，保持了方向和相互位置关系的正确。

在地图上保持方向和角度的正确是墨卡托投影的优点，墨卡托投影地图常用作航海图和航空图，如果循着墨卡托投影图上两点间的直线航行，方向不变可以一直到达目的地，因此它对船舰在航行中定位、确定航向都具有有利条件，给航海者带来很大方便。

百度地图和 Google Maps 使用的投影方法都是墨卡托投影。

经纬度转墨卡托投影坐标

/*
WGS84经纬度转墨卡托投影坐标，貌似1px相当于1m
平面坐标x = 经度*20037508.34/180
平面坐标y = log（tan（（90+纬度）*PI/360））/（PI/180）*20037508.34/180
*/
export function lonlatToMercator(lon, lat) {
  let x = (lon * 20037508.34) / 180;
  let y = Math.log(Math.tan(((90 + lat) * Math.PI) / 360)) / (Math.PI / 180);
  y = (y * 20037508.34) / 180;
  return { x, y };
}

/*
墨卡托投影坐标转WGS84经纬度坐标
经度 = 平面坐标x/20037508.34*180
纬度 = 180/(PI*(2*atan(exp(平面坐标y/20037508.34*180*PI/180))-PI/2)
*/
export function mercatorToLonlat(x, y) {
  let lon = (x / 20037508.34) * 180;
  let lat = (y / 20037508.34) * 180;
  lat = (180 / Math.PI) * (2 * Math.atan(Math.exp((lat * Math.PI) / 180)) - Math.PI / 2);
  return { lon, lat };
}

UTM 投影

UTM 投影全称为“通用横轴墨卡托投影”(UNIVERSAL TRANSVERSE MERCATOR PROJECTION)，是一种“等角横轴割圆柱投影”，椭圆柱割地球于南纬 80 度、北纬 84 度两条等高圈，投影后两条相割的经线上没有变形，而中央经线上长度比 0.9996。 UTM 投影是为了全球战争需要创建的，美国于 1948 年完成这种通用投影系统的计算。 与高斯-克吕格投影相似，该投影角度没有变形，中央经线为直线，且为投影的对称轴，中央经线的比例因子取 0.9996 是为了保证离中央经线左右约 180km 处有两条不失真的标准经线。 UTM 投影分带方法与高斯-克吕格投影相似，将北纬 84 度至南纬 80 度之间按经度分为 60 个带,每带 6 度.从西经 180 度起算,两条标准经线距中央经线为 180KM 左右,中央经线比例系数为 0.9996.

高斯-克吕格投影

高斯-克吕格投影（Gauss - Kruger projection）——等角横轴切椭圆柱投影

由德国数学家、物理学家、天文学家高斯于 19 世纪 20 年代拟定，后经德国大地测量学家克吕格于 1912 年对投影公式加以补充，故称为高斯-克吕格投影。

高斯-克吕格投影即等角横切椭圆柱投影。假想用一个椭圆柱横切于地球椭球体的某一经线上，这条与圆柱面相切的经线，称中央经线。以中央经线为投影的对称轴，将东西各 3° 或 1°30′的两条子午线所夹经差 6° 或 3° 的带状地区按数学法则、投影法则投影到圆柱面上，再展开成平面，即高斯-克吕格投影，简称高斯投影。这个狭长的带状的经纬线网叫做高斯-克吕格投影带。

高斯-克吕格投影特点：

1、中央子午线是直线，其长度不变形；其他子午线是凹向中央子午线的弧线，并以中央子午线为对称轴；

2、赤道线是直线，但有长度变形；其他纬线为凸向赤道的弧线，并以赤道为对称轴；

3、经线和纬线投影后仍然保持正交；

4、离开中央子午线越远，变形越大。

若采用分带投影的方法，可使投影边缘的变形不致过大。我国各种大、中比例尺地形图采用了不同的高斯-克吕格投影带。其中大于 1:1 万的地形图采用 3° 带；1:2.5 万至 1:50 万的地形图采用 6° 带。

正算公式

public double[] gaussBLtoXY(double longitude,double latitude) {
    int ProjNo = 0;

    // 带宽
    int ZoneWide = 6;

    double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
    double a, f, e2, ee, NN, T, C, A, M, iPI;

    // 3.1415926535898/180.0;
    iPI = 0.0174532925199433;

    // 54年北京坐标系参数
    a = 6378245.0;
    f = 1.0 / 298.3;

    // 80年西安坐标系参数
    // a=6378140.0;
    // f=1/298.257;

    ProjNo = (int) (longitude / ZoneWide);
    longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
    longitude0 = longitude0 * iPI;

    // 经度转换为弧度
    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;
    return new double[] { xval, yval };
}

反算公式

public double[] gaussXYtoBL(double X,double Y) {
    int ProjNo;
    int ZoneWide; // //带宽
    double[] output = new double[2];
    double longitude1, latitude1, longitude0, X0, Y0, xval, yval;// latitude0,
    double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
    iPI = 0.0174532925199433; // //3.1415926535898/180.0;
    a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数
    //a = 6378140.0;
    //f = 1 / 298.257; // 80年西安坐标系参数
    // 54年北京坐标系参数

    ZoneWide = 6; // //6度带宽
    ProjNo = (int) (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;
}

兰伯特等角圆锥投影

兰伯特等角圆锥投影也称兰勃脱正形圆锥投影，该投影的微分圆投影后仍为圆形。经线为辐射直线，纬线为同心圆圆弧。指定两条标准纬度线 Q1，Q2，在这两条纬度线上没有长度变形,即 M=N=1。此种投影也叫等角割圆锥投影，可用来编制中，小比例尺地图。等角圆锥投影有广泛的应用，特别适宜于作为中纬度处沿纬度线伸展的制图区域之投影，投影后经线为辐射直线,纬度线为同心圆圆弧。我国的分省图，即为两条标准纬度线为 Q1=25 度，Q2=45 度的兰伯特等角圆锥投影。1962 年以后，百万分一地图采用了等角圆锥投影(南纬度 80 度，北纬度 84 度)，极区附近，采用等角方位投影(极球面投影)。

地图分幅为:

纬度 60 以下，纬度差 4 经差 6 度分幅

纬度 60-76，纬度差 4 经差 12 度分幅

纬度 76-84，纬度差 4 经差 24 度分幅

纬度 84-88，纬度差 4 经差 36 度分幅

88-90 仍为一幅图

每幅图内两条标准纬线的纬度:

Q1=QS+40 分 (南纬度) Q2=QN-40 分(北纬度)

投影后经线是辐射直线，东西图幅可完全拼接，南北图幅有裂隙。

我国采用等角割圆锥，Q1=PHIS+35 分 Q2=PHIN-35 分

好文要顶 关注我 收藏该文

高斯-克吕格投影与 UTM 投影的异同

高斯-克吕格投影与 UTM 投影都是横轴墨卡托投影的变种，目前一些国外的软件或国外进口仪器的配套软件往往不支持高斯-克吕格投影，但支持 UTM 投影，因此常有把 UTM 投影当作高斯-克吕格投影的现象。

从投影几何方式看，高斯-克吕格投影是“等角横切椭圆柱投影”，投影后中央经线保持长度不变，即比例系数为 1；UTM 投影是“等角横轴割圆柱投影”，圆柱割地球于南纬 80 度、北纬 84 度两条等高圈，投影后两条割线上没有变形，中央经线上长度比 0.9996。

从计算结果看，两者主要差别在比例因子上，高斯-克吕格投影中央经线上的比例系数为 1，UTM 投影为 0.9996，高斯-克吕格投影与 UTM 投影可近似采用 X[UTM]=0.9996*X[高斯]，Y[UTM]=0.9996*Y[高斯]，进行坐标转换（注意：如坐标纵轴西移了 500000 米，转换时必须将 Y 值减去 500000 乘上比例因子后再加 500000）。

	输入坐标（度）	高斯投影（米）	UTM 投影（米）	Xutm=0.9996*X 高斯 / Yutm=0.9996*Y 高斯
纬度值（X）	32	3543600.9	3542183.5	3543600.9*0.9996≈3542183.5
经度值（Y）	121	21310996.8	311072.4	（310996.8-500000）*0.9996+500000≈311072.4

从分带方式看，两者的分带起点不同，高斯-克吕格投影自 0 度子午线起每隔经差 6 度自西向东分带，第 1 带的中央经度为 3°；UTM 投影自西经 180° 起每隔经差 6 度自西向东分带，第 1 带的中央经度为-177°，因此高斯-克吕格投影的第 1 带是 UTM 的第 31 带。

此外，两投影的东伪偏移（False Easting）都是 500 公里，高斯-克吕格投影北伪偏移（False_Northing）为零，UTM 北半球投影北伪偏移（False Northing）为零，南半球则为 10000 公里。