Ближайший сосед Python - координаты

Я хотел проверить, правильно ли использую дерево KD scipy, потому что оно работает медленнее, чем простой брутфорс.

У меня было три вопроса по этому поводу:

Q1.

Если я создам следующие тестовые данные:

nplen = 1000000
# WGS84 lat/long
point = [51.349,-0.19]
# This contains WGS84 lat/long
points = np.ndarray.tolist(np.column_stack(
        [np.round(np.random.randn(nplen)+51,5),
         np.round(np.random.randn(nplen),5)]))

И создайте три функции:

def kd_test(points,point):
    """ KD Tree"""
    return points[spatial.KDTree(points).query(point)[1]]

def ckd_test(points,point):
    """ C implementation of KD Tree"""
    return points[spatial.cKDTree(points).query(point)[1]]

def closest_math(points,point):
    """ Simple angle"""
    return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points))[1:3]   

Однако я ожидал, что дерево cKD будет самым быстрым - запустив это:

print("Co-ordinate: ", f(points,point))
print("Index: ", points.index(list(f(points,point))))
%timeit f(points,point)

Время результата - простой метод перебора работает быстрее:

closest_math: 1 loops, best of 3: 3.59 s per loop
ckd_test: 1 loops, best of 3: 13.5 s per loop
kd_test: 1 loops, best of 3: 30.9 s per loop

Это потому, что я как-то неправильно его использую?

Q2.

Я бы предположил, что даже для того, чтобы получить рейтинг (а не расстояние) ближайших точек, все же необходимо спроецировать данные. Однако кажется, что спроецированные и непроецированные точки дают мне одного и того же ближайшего соседа:

def proj_list(points,
              inproj = Proj(init='epsg:4326'),
              outproj = Proj(init='epsg:27700')):
    """ Projected geo coordinates"""
    return [list(transform(inproj,outproj,x,y)) for y,x in points]
proj_points = proj_list(points)
proj_point = proj_list([point])[0]

Это просто потому, что мой разброс точек недостаточно велик, чтобы вносить искажения? Я выполнял повторный запуск несколько раз и все равно получил тот же индекс из возвращаемых спроектированных и непроектированных списков.

Q3.

Является ли обычно более быстрым проецирование точек (как указано выше) и вычисление расстояния гипотенузы по сравнению с вычислением расстояния гаверсинуса или винсента на (непроектируемых) широте / долготе? И какой вариант будет более точным? Я провел небольшой тест:

from math import *
def haversine(origin,
              destination):
    """
    Find distance between a pair of lat/lng coordinates
    """
    lat1, lon1, lat2, lon2 = map(radians, [origin[0],origin[1],destination[0],destination[1]])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    c = 2 * asin(sqrt(a))
    r = 6371000  # Metres
    return (c * r)

def closest_math_unproj(points,point):
    """ Haversine on unprojected """
    return (min((haversine(point,pt),pt[0],pt[1]) for pt in points))

def closest_math_proj(points,point):
    """ Simple angle since projected"""
    return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points)) 

Результаты:

введите описание изображения здесь

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

Тестирование этого с помощью онлайн-расчета винсенти, по-видимому, является планируемым -координаты - это то, что нужно:

введите описание изображения здесь


person mptevsion    schedule 05.02.2016    source источник
comment
Одно в основном несвязанное предложение: использовать %timeit -n 10 f(points,point), вероятно, удобнее, чем использовать %timeit for x in range(10): f(points,point).   -  person Martin Valgur    schedule 05.02.2016
comment
Кстати, стоит взглянуть на github.com/storpipfugl/pykdtree. Это, вероятно, не решит проблемы эффективности по сравнению с методом грубой силы, но, вероятно, будет немного быстрее, чем реализация scipy по умолчанию.   -  person Martin Valgur    schedule 05.02.2016


Ответы (1)


Q1.

Причина очевидной неэффективности k-d-дерева довольно проста: вы одновременно измеряете как построение, так и запросы к k-d-дереву. Это не то, как вы могли бы или должны использовать k-d дерево: вы должны построить его только один раз. Если вы измеряете только запросы, затрачиваемое время сокращается до нескольких десятков миллисекунд (по сравнению с секундами при использовании подхода грубой силы).

Q2.

Это будет зависеть от пространственного распределения фактически используемых данных и используемой проекции. Могут быть небольшие различия в зависимости от того, насколько эффективна реализация k-d дерева для балансировки построенного дерева. Если вы запрашиваете только одну точку, то результат будет детерминированным и в любом случае не зависит от распределения точек.

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

Q3.

Технически ответ на ваш вопрос тривиален: использование формулы Хаверсина для измерения географических расстояний и точнее, и медленнее. Гарантированность компромисса между точностью и скоростью в значительной степени зависит от вашего варианта использования и пространственного распределения ваших данных (в основном, очевидно, от пространственной протяженности).

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

person Martin Valgur    schedule 05.02.2016
comment
Спасибо, Мартин - это все отвечает. Я просто хотел убедиться, что вы сказали, что формула Хаверсина будет более точной (и, следовательно, формула Винсенти). Что означает, что если точность очень важна, тогда вам подойдет векторизованная формула numpy vincenty? - person mptevsion; 05.02.2016
comment
Извините, что я получаю - если у меня (например) 10 миллионов координат в Великобритании, и моя основная цель - минимизировать ошибку в расстояниях (+ - 1 метр - это отлично), тогда я должен использовать scipy.pdist с векторизованной формулой винсенти вместо проецирования координат с последующим вычислением векторизованного евклидова расстояния? - person mptevsion; 05.02.2016
comment
Ах, прости. Я неправильно прочитал последний вопрос и пропустил, что вы спрашивали о формулах гаверсина или винсенти. Можете передать мой последний ответ. Последний вопрос, вероятно, гораздо лучше подходит для gis.stackexchange.com, чем для SO как такового. - person Martin Valgur; 05.02.2016