Python: ускорение географического сравнения

Question

Я написал некоторый код, который включает в себя вложенный цикл, в котором внутренний цикл выполняется около 1,5 миллионов раз. У меня есть функция в этом цикле, которую я пытаюсь оптимизировать. Я проделал некоторую работу и получил некоторые результаты, но мне нужно немного информации, чтобы проверить, насколько разумно то, что я делаю.

Некоторый фон:

У меня есть две коллекции географических точек (широта, долгота), одна относительно небольшая коллекция и одна относительно огромная коллекция. Для каждой точки в маленькой коллекции мне нужно найти ближайшую точку в большой коллекции.

Очевидный способ сделать это - использовать формулу haversine. Преимущество здесь в том, что расстояния точно точны.

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

def haversine(point1, point2):
    """Gives the distance between two points on earth.
    """
    earth_radius_miles = 3956
    lat1, lon1 = (radians(coord) for coord in point1)
    lat2, lon2 = (radians(coord) for coord in point2)
    dlat, dlon = (lat2 - lat1, lon2 - lon1)
    a = sin(dlat/2.0)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.0)**2
    great_circle_distance = 2 * asin(min(1,sqrt(a)))
    d = earth_radius_miles * great_circle_distance
    return d

Тем не менее, выполнение этого 1,5 миллиона раз занимает около 9 секунд на моей машине (в зависимости от времени). Поскольку точное расстояние не имеет значения, мне нужно только найти ближайшую точку, поэтому я решил попробовать другие функции.

Простая реализация теоремы Пифагора дает мне ускорение примерно на 30%. Думая, что я могу сделать лучше, я написал следующее:

def dumb(point1, point2):
    lat1, lon1 = point1
    lat2, lon2 = point2
    d = abs((lat2 - lat1) + (lon2 - lon1))

, что дает мне коэффициент улучшения в 10 раз. Однако теперь я беспокоюсь, что это не сохранит неравенство треугольника.

Итак, мой последний вопрос состоит из двух частей: я хотел бы иметь функцию, которая работает так же быстро, как dumb, но все же будет правильной. Будет ли dumb работать? Если нет, какие-либо предложения о том, как улучшить мою функцию haversine?

Answer 1

Это расчет, который действительно хорош в numpy .Вместо того, чтобы зацикливаться на всем большом наборе координат, вы можете вычислить расстояние между одной точкой и всем набором данных в одном вычислении.С моими тестами, приведенными ниже, вы можете получить увеличение скорости на порядок.

Вот некоторые временные тесты с вашим методом haversine, вашим методом dumb (не совсем уверен, что это делает) и моим методом тупой хаверсин,Он вычисляет расстояние между двумя точками - одной в Вирджинии и одной в Калифорнии, которые находятся на расстоянии 2293 миль.

from math import radians, sin, cos, asin, sqrt, pi, atan2
import numpy as np
import itertools

earth_radius_miles = 3956.0

def haversine(point1, point2):
    """Gives the distance between two points on earth.
    """
    lat1, lon1 = (radians(coord) for coord in point1)
    lat2, lon2 = (radians(coord) for coord in point2)
    dlat, dlon = (lat2 - lat1, lon2 - lon1)
    a = sin(dlat/2.0)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.0)**2
    great_circle_distance = 2 * asin(min(1,sqrt(a)))
    d = earth_radius_miles * great_circle_distance
    return d

def dumb(point1, point2):
    lat1, lon1 = point1
    lat2, lon2 = point2
    d = abs((lat2 - lat1) + (lon2 - lon1))
    return d

def get_shortest_in(needle, haystack):
    """needle is a single (lat,long) tuple.
        haystack is a numpy array to find the point in
        that has the shortest distance to needle
    """
    dlat = np.radians(haystack[:,0]) - radians(needle[0])
    dlon = np.radians(haystack[:,1]) - radians(needle[1])
    a = np.square(np.sin(dlat/2.0)) + cos(radians(needle[0])) * np.cos(np.radians(haystack[:,0])) * np.square(np.sin(dlon/2.0))
    great_circle_distance = 2 * np.arcsin(np.minimum(np.sqrt(a), np.repeat(1, len(a))))
    d = earth_radius_miles * great_circle_distance
    return np.min(d)


x = (37.160316546736745, -78.75)
y = (39.095962936305476, -121.2890625)

def dohaversine():
    for i in xrange(100000):
        haversine(x,y)

def dodumb():
    for i in xrange(100000):
        dumb(x,y)

lots = np.array(list(itertools.repeat(y, 100000)))
def donumpy():
    get_shortest_in(x, lots)

from timeit import Timer
print 'haversine distance =', haversine(x,y), 'time =',
print Timer("dohaversine()", "from __main__ import dohaversine").timeit(100)
print 'dumb distance =', dumb(x,y), 'time =',
print Timer("dodumb()", "from __main__ import dodumb").timeit(100)
print 'numpy distance =', get_shortest_in(x, lots), 'time =',
print Timer("donumpy()", "from __main__ import donumpy").timeit(100)

И вот что он печатает:

haversine distance = 2293.13242188 time = 44.2363960743
dumb distance = 40.6034161104 time = 5.58199882507
numpy distance = 2293.13242188 time = 1.54996609688

Метод numpy занимает 1,55 секунд, чтобы вычислить то же количество вычислений расстояния, сколько требуется 44,24 секунд для вычисления с помощью метода вашей функции.Вероятно, вы могли бы добиться большего ускорения, объединив некоторые из простых функций в один оператор, но он стал бы длинной, трудной для чтения строкой.

Answer 2

Вы можете рассмотреть какое-то графическое хеширование, то есть быстро найти ближайшие точки и затем рассчитать их. Например, вы можете создать равномерную сетку и распределить точки (из большой коллекции) в ячейках, созданных сеткой.

Теперь, имея точку из небольшой коллекции, вам нужно будет обработать гораздо меньшее количество точек (т. Е. Только в соответствующих ячейках)

Answer 3

abs(lat2 - lat1) + abs(lon2 - lon1) является 1-нормой или манхэттен-метрикой, и поэтому неравенство треугольника выполнено.

Answer 4

Формула, которую вы написали (d = abs (lat2-lat1) + (lon2-lon1)), НЕ сохраняет неравенство треугольника: если вы найдете lat, lon, для которого d равно min, вы не найдете ближайшую точку,но точка, ближайшая к двум диагональным прямым линиям, пересекающим точку, которую вы проверяете!

Я думаю, вам следует упорядочить большое количество точек по широте и долготе (это означает: (1,1), (1,2), (1,3) ... (2,1), (2,2) и т. Д. Затем используйте метод наводчика, чтобы найти некоторые из ближайших точек по широте и долготе (это должно быть очень быстро,это займет время процессора, пропорциональное ln2 (n), где n - количество точек.) Вы можете сделать это легко, например: выберите все точки в квадрате 10x10 вокруг точки, которую вы собираетесь проверить, этоозначает: найти все точки от -10 до +10 в латах (метод наводчика) и снова те, которые от -10 до +10 в долготе (метод наводчика). Теперь у вас есть действительно небольшое количество данных, обрабатываемых процессом,и это должно быть очень быстро!

Answer 5

У меня была похожая проблема, и я решил запустить функцию Cython.На моем MBP 2008 он может выполнять около 1,2 млн итераций в секунду.Проверка типа ускоряет еще на 25%.Без сомнения, возможна дальнейшая оптимизация (за счет ясности).

Вы также можете проверить функцию scipy.spatial.distance.cdist.

from libc.math cimport sin, cos, acos

def distance(float lat1, float lng1, float lat2, float lng2):
    if lat1 is None or lat2 is None or lng1 is None or lng2 is None: return None
    cdef float phi1
    cdef float phi2
    cdef float theta1
    cdef float theta2
    cdef float c
    cdef float arc

    phi1 = (90.0 - lat1)*0.0174532925
    phi2 = (90.0 - lat2)*0.0174532925
    theta1 = lng1*0.0174532925
    theta2 = lng2*0.0174532925

    c = (sin(phi1)*sin(phi2)*cos(theta1 - theta2) + cos(phi1)*cos(phi2))
    arc = acos( c )
    return arc*6371

Answer 6

Самый быстрый способ сделать это - избежать вычисления функции для каждой пары точек, при условии, что ваша относительно небольшая коллекция не очень мала.

Есть несколько баз данных, которые имеют гео-индексы, которые вы можете использовать(mysql, oracle, mongodb ..) или реализуйте что-нибудь самостоятельно.

Вы можете использовать python-geohash .Для каждого документа в меньшей коллекции вам нужно быстро найти набор документов в большой коллекции, которые имеют хэш от geohash.neighbors для самого длинного размера хеша, который соответствует.Вам понадобится использовать соответствующую структуру данных для поиска, или это будет медленно.

Для нахождения расстояния между точками ошибка простого подхода увеличивается с увеличением расстояния между точками, а также зависит отширота.См., Например, http://www.movable -type.co.uk / scripts / gis-faq-5.1.html .