LBS 球面距离公式

博客首页 » LBS 球面距离公式

发布于 29 Aug 2013 08:19
标签 blog lbs
WGS_84_reference_frame.png
工作中需要考虑WGS 84坐标体系中的球面距离,许多文章都是只有结果没有原理。终于找到一个描述球面距离公式的文章

关于坐标系

可以参考以下文章,就不赘述了。
http://baike.baidu.com/view/261065.htm
http://wenku.baidu.com/view/4a9c432de2bd960590c67775.html
http://en.wikipedia.org/wiki/World_Geodetic_System
http://ja.wikipedia.org/wiki/%E6%B8%AC%E5%9C%B0%E7%B3%BB

球面距离公式

http://tech.idv2.com/2011/06/17/location-search/

附近地点搜索初探

附近地点搜索,顾名思义,就是搜索用户附近有哪些地点。随着GPS和带有GPS功能的移动设备的普及, 附近地点搜索也变得炙手可热。不过在网上却很少有这方面的讨论。本文的方法并不算最好, 但足以应付一般的应用了。

本文中,数据库采用Oracle,语言采用python。理论上别的数据库和语言也没问题, 但我们要在经纬度上设置两个索引,所以如果你的数据库不支持索引,或者不支持在一个查询中使用两个索引, 那就只能想别的办法了。

球面最短距离公式

球面上任意两点之间的距离计算公式可以参考维基百科上的下述文章,这里就不再赘述了。

Great-circle distance
http://en.wikipedia.org/wiki/Great-circle_distance

Haversine formula
http://en.wikipedia.org/wiki/Haversine_formula

值得一提的是,维基百科推荐使用Haversine公式,理由是Great-circle distance公式用到了大量余弦函数, 而两点间距离很短时(比如地球表面上相距几百米的两点),余弦函数会得出0.999…的结果, 会导致较大的舍入误差。而Haversine公式采用了正弦函数,即使距离很小,也能保持足够的有效数字。 以前采用三角函数表计算时的确会有这个问题,但经过实际验证,采用计算机来计算时,两个公式的区别不大。 稳妥起见,这里还是采用Haversine公式。

haversin(d/R) = haversin(φ2 - φ1) + cos(φ1) cos(φ2) haversin(⊿λ)

distance-haversin-distance.png

其中

haversin(θ) = sin(θ/2)^2 = (1- cos(θ))/2

distance-haversin.png

R为地球半径,可取平均值 6371.137km
φ1, φ2 表示两点的纬度;
Δλ 表示两点经度的差值。
距离计算函数
下面就是计算球面间两点(lat0, lng)-(lat1, lng1)之间距离的函数。

from math import sin, asin, cos, radians, fabs, sqrt
 
EARTH_RADIUS=6371           # 地球平均半径,6371km
 
def hav(theta):
    s = sin(theta / 2)
    return s * s
 
def get_distance_hav(lat0, lng0, lat1, lng1):
    "用haversine公式计算球面两点间的距离。"
    # 经纬度转换成弧度
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lng0 = radians(lng0)
    lng1 = radians(lng1)
 
    dlng = fabs(lng0 - lng1)
    dlat = fabs(lat0 - lat1)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))
 
    return distance
范围搜索算法

在庞大的地理数据库中搜索地点,索引是很重要的。但是,我们的需求是搜索附近地点, 例如,坐标(39.91, 116.37)附近500米内有什么地点?搜索条件是地点坐标与当前坐标之间的距离, 显然是无法应用索引的。

那么换个思路:首先算出“给定坐标附近500米”这个范围的坐标范围。 虽然它是个圆,但我们可以先求出该圆的外接正方形,然后拿正方形的经纬度范围去搜索数据库。

distance-map.png

如图,红色部分为要求的搜索范围,绿色部分为实际搜索范围。

先来求东西两侧的的范围边界。在haversin公式中令φ1 = φ2,可得

⊿λ = 2 arcsin (sin (d/2R / cosφ)

distance-lng.png

写成python代码就是

dlng = 2 * asin(sin(distance / (2 * EARTH_RADIUS)) / cos(lat))
dlng = degrees(dlng)        # 弧度转换成角度

然后求南北两侧的范围边界,在haversin公式中令 Δλ = 0,可得

⊿φ = d/R

distance-lat.png

写成python代码就是

dlat = distance / EARTH_RADIUS
dlng = degrees(dlat)     # 弧度转换成角度

这样,根据当前点坐标,我们可以得出搜索范围为

left-top : (lat + dlat, lng - dlng)
right-top : (lat + dlat, lng + dlng)
left-bottom : (lat - dlat, lng - dlng)
right-bottom: (lat - dlat, lng + dlng)
然后利用这个范围构造SQL语句,即可实现范围查询:

SELECT * FROM place WHERE lat > lat1 AND lat < lat2 AND lng > lng1 AND lng < lng2;

在lat和lng列上建立索引,能从一定程度上提高范围查询的效率。

不过,这样查询到的地点是正方形范围内的地点,一些结果与当前点的距离可能会超出给定的距离。 如果要求严格,可以遍历结果并计算与当前点之间的距离,并过滤掉不符合要求的结果。

总结

附近地点搜索条件是距离,而数据库中一般只保存地点的经纬度,因此无法直接查询。 本文将距离转化成经纬度范围,利用经纬度上的索引,提高查询效率。

Javascript实现
rad = function(x) {return x*Math.PI/180;}
 
distHaversine = function(p1, p2) {
  var R = 6371; // earth's mean radius in km
  var dLat  = rad(p2.lat() - p1.lat());
  var dLong = rad(p2.lng() - p1.lng());
 
  var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
          Math.cos(rad(p1.lat())) * Math.cos(rad(p2.lat())) * Math.sin(dLong/2) * Math.sin(dLong/2);
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
  var d = R * c;
 
  return d.toFixed(3);
}

C#的实现

计算公式如下:
  说明:
  1、公式中经纬度均用弧度表示;
  2、Lat1 Lng1 分别表示A点经、纬度,Lat2 Lng2 分别表示B点经纬度;
  3、a=Lng1 -Lng2 为两点纬度之差 b=Lat1 – Lat2 为两点经度之差;
  4、6378.137为地球半径,单位为公里;
  c#实现函数如下:

 public double DistanceOfTwoPoints(double lat1, double lng1, double lat2, double lng2)
 {
    double radLng1 = lng1 * Math.PI / 180.0;
    double radLng2 = lng2 * Math.PI / 180.0;
    double a = radLng1 - radLng2;
    double b = (lat1-lat2) * Math.PI / 180.0;
    double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a/2), 2) +
    Math.Cos(radLng1) * Math.Cos(radLng2) * Math.Pow(Math.Sin(b/2), 2)))* 6378.137; 
    s =Math.Round(s * 10000) / 10000;
    return s;
 }

返回结果单位为公里。
 public class DistanceAlgorithm
{
    const double PIx = 3.141592653589793;
    const double RADIO = 6378.16;

    /// <summary>
    /// This class cannot be instantiated.
    /// </summary>
    private DistanceAlgorithm() { }

    /// <summary>
    /// Convert degrees to Radians
    /// </summary>
    /// <param name="x">Degrees</param>
    /// <returns>The equivalent in radians</returns>
    public static double Radians(double x)
    {
        return x * PIx / 180;
    }

    /// <summary>
    /// Calculate the distance between two places.
    /// </summary>
    /// <param name="lon1"></param>
    /// <param name="lat1"></param>
    /// <param name="lon2"></param>
    /// <param name="lat2"></param>
    /// <returns></returns>
    public static double DistanceBetweenPlaces(
        double lon1,
        double lat1,
        double lon2,
        double lat2)
    {
        double dlon =  Radians(lon2 - lon1);
        double dlat =  Radians(lat2 - lat1);

        double a = (Math.Sin(dlat / 2) * Math.Sin(dlat / 2)) + Math.Cos(Radians(lat1)) * Math.Cos(Radians(lat2)) * (Math.Sin(dlon / 2) * Math.Sin(dlon / 2));
        double angle = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
        return (angle * RADIO) * 0.62137;//distance in miles
    }

}

Reference:
http://stackoverflow.com/questions/1502590/calculate-distance-between-two-points-in-google-maps-v3


本页面的文字允许在知识共享 署名-相同方式共享 3.0协议和GNU自由文档许可证下修改和再使用,仅有一个特殊要求,请用链接方式注明文章引用出处及作者。请协助维护作者合法权益。


系列文章

文章列表

  • LBS 球面距离公式

这篇文章对你有帮助吗,投个票吧?

rating: 0+x

留下你的评论

Add a New Comment