Python геокод фильтрации по расстоянию

Question

Мне нужно отфильтровать геокоды по близости к местоположению. Например, я хочу отфильтровать список геокодов ресторанов, чтобы идентифицировать эти рестораны в пределах 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

Примечание: я использую верхнюю норму для расстояния.

Answer 1

Кажется, не было хорошей реализации Python. К счастью, боковая панель «Связанные статьи» - наш друг. Эта статья SO указывает на превосходную статью , в которой приведены математические и Java-реализации. Фактическая функция, которая вам требуется, довольно коротка и встроена в мой код Python ниже. Протестировано в показанной степени. Прочитайте предупреждения в комментариях.

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

Answer 2

Вот как вы рассчитываете расстояние между широтой / длинной парой, используя формулу haversine:

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» (также в км) относительно вашего порога.Если вам нужно что-то еще, кроме км, отрегулируйте радиус.

Извините, я не могу дать вам выпадающее решение, но я не понимаю ваш скелет кода (см. Комментарий).

Также обратите внимание, что в наши дни вы, вероятно, хотите использовать сферический закон косинусов, а не хаверсин.Преимущества в числовой стабильности больше не стоят того, и это чертовски просто для понимания, кодирования и использования.

Answer 3

Ответ Джона Мачина мне очень помог.Есть только небольшая ошибка: широта и долгота поменялись местами на boundigbox:

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

, это решает проблему.Причина в том, что долготы не меняют свое расстояние на градус, а широты.Их расстояние зависит от долготы.

Answer 4

Если вы храните данные в MongoDB, он выполняет хорошо индексируемые поиски геолокации и превосходит приведенные выше решения на чистом Python, поскольку он будет выполнять оптимизацию для вас.

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