使用Ceres实现WGS84到GCJ-02坐标相互转换

• 前言

最近偶然才知道有GCJ-02这个坐标系，作为测绘从业人员，直到现在才知道我国还存在这么一个坐标系，实在感觉有些对不起专业。不过我主要做航摄，平时也主要接触的是真实大地坐标，航摄成图的后续处理也从未关注过，之前确实没有接触过这一坐标系统。

目前WGS84到GCJ-02坐标的转换算法似乎已经是公开的，公开搜索到了好几个WGS84到GCJ-02的转换代码，核心算法都是一致的，大家似乎使用了同一加密公式。假设这一公式确实是李成民实现WGS84到GCJ-02坐标转换使用的精确公式，那么从加密过程看，使用的公式也只是普通的多项式，并没有使用其他特殊处理。相应的反解时，基于加密公式列非线性方程，使用牛顿迭代法完全可以获取非线性方程的数值解。理论上给定迭代收敛条件，完全可以获得较高精度的迭代数值解。

即使网上流传的加密算法使用的加密公式并非原版，只是近似加密公式。只要李成民在对WGS84坐标加密时并没有做特殊处理，两个坐标系的坐标变换还是通过多项式变换，高精度的纠偏算法还是可以实现的。主要是现在有ceres这一神器，使用ceres的数值求导完全可以解决这一问题。将GCJ-02坐标当做观测值列误差方程，调用国家保密插件配合ceres数值求导计算residuals就行。

在知乎跟踪这个问题时，我表达了类似的看法，有知友不相信我的结论，于是这里我借助ceres实现WGS84到GCJ-02坐标大相互转换。因为我手上也没有国家保密插件，只能当网上公开的加密算法使用的加密公式是精确公式。

• 公开的加、解密算法

这里我收集了2个版本的加解密算法。

第一个是来自博客园的文章《GCJ-02火星坐标系和WGS-84坐标系转换关系》，代码如下：

# -*- coding: utf-8 -*-
import json
import math

x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 扁率

def wgs84togcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return lng, lat
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [mglng, mglat]


def gcj02towgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return lng, lat
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [lng * 2 - mglng, lat * 2 - mglat]


def transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
        0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
        0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    if lng < 72.004 or lng > 137.8347:
        return True
    if lat < 0.8293 or lat > 55.8271:
        return True
    return False


if __name__ == '__main__':
    [lng,lat]=[114.061202,22.529388]
    [dstlng, dstlat] = gcj02towgs84(lng, lat)
    print(dstlng, dstlat)

第二个版本是知乎上一篇文章《【对错代码分辨】wgs84坐标系和火星坐标系的转换中的对与错》引用的QGIS相关插件作者geohey提供的转换代码：

# -*- coding: utf-8 -*-
##########################################################################################
"""
/***************************************************************************
 OffsetWGS84Core
                                 A QGIS plugin
 Class with methods for geometry and attributes processing
                              -------------------
        begin                : 2016-10-11
        git sha              : $Format:%H$
        copyright            : (C) 2017 by sshuair
        email                : [email protected]
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from __future__ import print_function
##########################################################################################
from builtins import zip
import math
from math import sin, cos, sqrt, fabs, atan2
from math import pi as PI
# from numba import jit


# =================================================sshuair==================================