Python地理编码按距离过滤答案

【问题标题】：Python geocode filtering by distancePython地理编码按距离过滤
【发布时间】：2011-03-12 01:38:54
【问题描述】：

我需要过滤地理编码以接近某个位置。例如，我想过滤餐厅地理编码列表，以识别我当前位置 10 英里范围内的餐厅。

有人能指点我一个将距离转换为纬度和经度增量的函数吗？例如：

class GeoCode(object):
   """Simple class to store geocode as lat, lng attributes."""
   def __init__(self, lat=0, lng=0, tag=None):
      self.lat = lat
      self.lng = lng
      self.tag = None

def distance_to_deltas(geocode, max_distance):
   """Given a geocode and a distance, provides dlat, dlng
      such that

         |geocode.lat - dlat| <= max_distance
         |geocode.lng - dlng| <= max_distance
   """
   # implementation
   # uses inverse Haversine, or other function?
   return dlat, dlng

注意：我使用的是最高标准来表示距离。

【问题讨论】：

对不起，我不明白。您是否希望 inverse_havesine 返回一个采用“其他”参数并返回 True 或 False 的可调用对象？还是您打算以其他方式传递“其他”？
(1) “有人能指点我吗”：someone == google (2) “提供 dlat, dlng 这样”的东西没有提到 dlat, dlng -- 请编辑你的问题。 (3) 什么是“距离度量的最高范数”？
@约翰·马钦。当然，谷歌也是你了解最高标准的朋友。
@Ranieri 抱歉，编辑了函数以澄清 dlat、dlng。逆半正弦应该给出相对于参考地理编码的最大 lat/lng 增量，以定义参考地理编码周围的 max_distance 正方形区域。
@AndrewB：我知道什么是“最高标准”。我问“什么是“距离度量的最高范数”。此外，“逆haversine”是一个与haversine函数相反的函数（@ 987654321@）；您滥用该术语来描述您的问题，这在某些感觉更容易的 (point1, point2) -> 距离问题的反面（可以使用 hasrsines 或其他方法计算）。

标签： python geocoding

【解决方案1】：

似乎没有一个好的 Python 实现。幸运的是，SO“相关文章”侧边栏是我们的朋友。 This SO article 指向一个excellent article，它给出了数学和Java 实现。您需要的实际功能相当短，并且嵌入在下面的 Python 代码中。测试到显示的程度。阅读 cmets 中的警告。

from math import sin, cos, asin, sqrt, degrees, radians

Earth_radius_km = 6371.0
RADIUS = Earth_radius_km

def haversine(angle_radians):
    return sin(angle_radians / 2.0) ** 2

def inverse_haversine(h):
    return 2 * asin(sqrt(h)) # radians

def distance_between_points(lat1, lon1, lat2, lon2):
    # all args are in degrees
    # WARNING: loss of absolute precision when points are near-antipodal
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    dlat = lat2 - lat1
    dlon = radians(lon2 - lon1)
    h = haversine(dlat) + cos(lat1) * cos(lat2) * haversine(dlon)
    return RADIUS * inverse_haversine(h)

def bounding_box(lat, lon, distance):
    # Input and output lats/longs are in degrees.
    # Distance arg must be in same units as RADIUS.
    # Returns (dlat, dlon) such that
    # no points outside lat +/- dlat or outside lon +/- dlon
    # are <= "distance" from the (lat, lon) point.
    # Derived from: http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
    # WARNING: problems if North/South Pole is in circle of interest
    # WARNING: problems if longitude meridian +/-180 degrees intersects circle of interest
    # See quoted article for how to detect and overcome the above problems.
    # Note: the result is independent of the longitude of the central point, so the
    # "lon" arg is not used.
    dlat = distance / RADIUS
    dlon = asin(sin(dlat) / cos(radians(lat)))
    return degrees(dlat), degrees(dlon)

if __name__ == "__main__":

    # Examples from Jan Matuschek's article

    def test(lat, lon, dist):
        print "test bounding box", lat, lon, dist
        dlat, dlon = bounding_box(lat, lon, dist)
        print "dlat, dlon degrees", dlat, dlon
        print "lat min/max rads", map(radians, (lat - dlat, lat + dlat))
        print "lon min/max rads", map(radians, (lon - dlon, lon + dlon))

    print "liberty to eiffel"
    print distance_between_points(40.6892, -74.0444, 48.8583, 2.2945) # about 5837 km
    print
    print "calc min/max lat/lon"
    degs = map(degrees, (1.3963, -0.6981))
    test(*degs, dist=1000)
    print
    degs = map(degrees, (1.3963, -0.6981, 1.4618, -1.6021))
    print degs, "distance", distance_between_points(*degs) # 872 km

【讨论】：

【解决方案2】：

这是您使用半正弦公式计算纬度/经度对之间距离的方法：

import math 

R = 6371 # km
dLat = (lat2-lat1) # Make sure it's in radians, not degrees
dLon = (lon2-lon1) # Idem 
a = math.sin(dLat/2) * math.sin(dLat/2) +
    math.cos(lat1) * math.cos(lat2) * 
    math.sin(dLon/2) * math.sin(dLon/2) 
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) 
d = R * c;

现在根据您的阈值测试“d”（也以公里为单位）是微不足道的。如果你想要公里以外的东西，请调整半径。

很抱歉，我无法为您提供即插即用的解决方案，但我不了解您的代码框架（请参阅评论）。

另外请注意，现在您可能希望使用余弦球面定律而不是Haversine。数值稳定性的优势已经不值得了，而且理解、编码和使用都非常简单。

【讨论】：

是的，我有半正弦公式。我正在寻找反转的实现。
修复了代码 sn-p。基本上我想计算地理编码的 max_distance 范围内的最大/最小纬度/经度值。所以 inverse_havesine 应该返回一个 dlat, dlng 元组，它们是仍然在范围内的 +/- 增量。这有帮助吗？

【解决方案3】：

如果您将数据存储在 MongoDB 中，它会为您提供很好的索引地理定位搜索，并且优于上述纯 Python 解决方案，因为它会为您处理优化。

http://www.mongodb.org/display/DOCS/Geospatial+Indexing

【讨论】：

克里斯·哥伦布不会用这个装备“发现”美国：“当前的实现假设一个平坦地球的理想模型，这意味着纬度 (y) 和经度 (x) 的弧度表示到处都是一样的距离。”
是的，那将成为两极附近的限制。中心点的位置 N 或 S 将被视为比应有的更近，而位置 E 或 W 将被视为更远。

【解决方案4】：

John Machin 的回答对我帮助很大。只是有个小错误：boundigbox中的经纬度互换了：

dlon = distance / RADIUS
dlat = asin(sin(dlon) / cos(radians(lon)))
return degrees(dlat), degrees(dlon)

这解决了问题。原因是经度不会改变每度的距离 - 但纬度会改变。它们的距离取决于经度。

【讨论】：