百度坐标系下经纬度距离的计算

最近的一个项目中有用到经纬度距离的计算，数据库中存储的是百度的经纬度。由于先前学习过一些 地图经纬度坐标系的知识，所以心中产生了困惑：使用随机偏移过的经纬度地址计算出来的距离是否是正确的？重新梳理了一些基本知识：

地球坐标 (WGS84)

• 国际标准，从 GPS 设备中取出的数据的坐标系
• 国际地图提供商使用的坐标系

火星坐标 (GCJ-02)也叫国测局坐标系

• 中国标准，从国行移动设备中定位获取的坐标数据使用这个坐标系
• 国家规定： 国内出版的各种地图系统（包括电子形式），必须至少采用GCJ-02对地理位置进行首次加密。

百度坐标 (BD-09)

• 百度标准，百度 SDK，百度地图，Geocoding 使用(本来就乱了，百度又在火星坐标上来个二次加密)

从设备获取经纬度（GPS）坐标

• 如果使用的是百度sdk那么可以获得百度坐标（bd09）或者火星坐标（GCJ02),默认是bd09
• 如果使用的是ios的原生定位库，那么获得的坐标是WGS84
• 如果使用的是高德sdk,那么获取的坐标是GCJ02

更多参考信息：

• https://en.wikipedia.org/wiki/Restrictions_on_geographic_data_in_China
• https://zh.wikipedia.org/w/index.php?oldid=44719670

最准确的距离计算应当采用WGS84的坐标，采用火星坐标系和百度坐标从原理上就可能存在偏差。但是是否可以对经纬度直接进行计算呢，考虑到百度地图有测量距离的功能，于是拔了他的代码用来研究：

代码来源： http://api.map.baidu.com/getscript?v=2.0

var j = void 0,
        p = null,
H.prototype.nb = function(a) {
        return a && this.lat == a.lat && this.lng == a.lng
};

getDistance: function(a, b) {
            if (a && b) {
                if (a.nb(b)) return 0;
                var c = 0,
                    c = R.To(a, b);
                if (c === p || c === j) c = 0;
                return c
            }
        },

nb: function(a) {
            return !(a instanceof eb) || this.xj() ? q : this.se().nb(a.se()) && this.of().nb(a.of())
        },

function R() {}
    R.prototype = new gc;
    x.extend(R, {
        To: function(a, b) {
            if (!a || !b) return 0;
            a.lng = this.ND(a.lng, -180, 180);
            a.lat = this.RD(a.lat, -74, 74);
            b.lng = this.ND(b.lng, -180, 180);
            b.lat = this.RD(b.lat, -74, 74);
            return this.Re(this.Sk(a.lng), this.Sk(b.lng), this.Sk(a.lat), this.Sk(b.lat))
        },
Re: function(a, b, c, d) {
            return this.jP * Math.acos(Math.sin(c) * Math.sin(d) + Math.cos(c) * Math.cos(d) * Math.cos(b - a))
        },
		jP: 6370996.81,

		RD: function(a, b, c) {
            b != p && (a = Math.max(a, b));
            c != p && (a = Math.min(a, c));
            return a
        },
        ND: function(a, b, c) {
            for (; a > c;) a -= c - b;
            for (; a < b;) a += c - b;
            return a
        }
		Sk: function(a) {
            return Math.PI * a / 180
        },        
});

使用Python对上述代码进行了重写：

def get_distance_bd09(point_a, point_b):
    """
    算法来源：http://developer.baidu.com/map/jsdemo.htm#a6_1
    :param pointA: {lat:29.490295, lng:106.486654}
    :param pointB: {lat:29.615467, lng:106.581515}
    :return:米
    """
    Radius = 6370996.81 #球半径

    if (point_a and point_b):
        if point_a["lat"] == point_b["lat"] and point_a["lng"] == point_b["lng"]:
            distance = 0
        else:
            a_lat = point_a["lat"] * math.pi / 180
            a_lng = point_a["lng"] * math.pi / 180
            b_lat = point_b["lat"] * math.pi / 180
            b_lng = point_b["lng"] * math.pi / 180
            # print(a_lng,b_lng,a_lat,b_lat)
            distance = Radius * math.acos(math.sin(a_lat) * math.sin(b_lat) + math.cos(a_lat) * math.cos(b_lat) * math.cos(b_lng - a_lng))
        print(distance)

发现此计算方法与wgs84的算法完全一致。

def get_distance_wgs84(point_a, point_b):
    """
    算法来源：https://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula
    :param pointA: {lat:29.490295, lng:106.486654}
    :param pointB: {lat:29.615467, lng:106.581515}
    :return:米
    """
    Radius = 6370996.81 #球半径

    d_lat = math.radians(point_b["lat"] - point_a["lat"])
    d_lng = math.radians(point_b["lng"] - point_a["lng"])

    a = math.sin(d_lat/2) * math.sin(d_lat/2) + \
        math.cos(math.radians(point_a["lat"])) * math.cos(math.radians(point_b["lat"])) * \
        math.sin(d_lng/2) * math.sin(d_lng/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = Radius * c
    return d

也就是说百度忽略了经纬度的偏移，可能原因是经纬度距离越是距离的近，算出的距离偏差越小，越是距离远，算出的偏差才会远。

补充知识

地理空间距离计算方法较多，可以分为两类：

• 球面模型，这种模型将地球看成一个标准球体，球面上两点之间的最短距离即大圆弧长，这种方法使用较广
• 椭球模型，该模型最贴近真实地球，精度也最高，但计算较为复杂

基于球面模型的地理空间距离计算公式，将地球看成圆球，假设地球上有A(ja,wa)，B(jb,wb)两点（注：ja和jb分别是A和B的经度，wa和wb分别是A和B的纬度），A和B两点的球面距离就是AB的弧长，AB弧长=R*角AOB（注：角AOB是A跟B的夹角，O是地球的球心，R是地球半径，约为6367000米）。如何求出角AOB呢？可以先求AOB的最大边AB的长度，再根据余弦定律可以求夹角。

如果业务场景仅仅是在一个城市范围内进行距离计算，也就是说两个点之间的距离一般不会超过200多千米。由于范围小，可以认为经线和纬线是垂直的，如图所示，要求A（116.8，39,78）和B（116.9，39.68）两点的距离，我们可以先求出南北方向距离AM，然后求出东西方向距离BM，最后求矩形对角线距离，即sqrt(AMAM + BMBM)。

public static double distanceSimplify(double lat1, double lng1, double lat2, double lng2, double[] a) {
     double dx = lng1 - lng2; // 经度差值
     double dy = lat1 - lat2; // 纬度差值
     double b = (lat1 + lat2) / 2.0; // 平均纬度
     double Lx = toRadians(dx) * 6367000.0* Math.cos(toRadians(b)); // 东西距离
     double Ly = 6367000.0 * toRadians(dy); // 南北距离
     return Math.sqrt(Lx * Lx + Ly * Ly);  // 用平面的矩形对角距离公式计算总距离
    }
}

参考内容：

• https://www.biaodianfu.com/lbs-precision.html
• https://tech.meituan.com/lucene-distance.html
• http://www.movable-type.co.uk/scripts/gis-faq-5.1.html