Эффективно находить индексы ближайших точек на непрямоугольной 2D-сетке

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

Я написал функцию, которая делает это, но она ДЕЙСТВИТЕЛЬНО медленная.

def find_indices(lon,lat,x,y):
    lonlat = np.dstack([lon,lat])
    delta = np.abs(lonlat-[x,y])
    ij_1d = np.linalg.norm(delta,axis=2).argmin()
    i,j = np.unravel_index(ij_1d,lon.shape)
    return i,j

ind = [find_indices(lon,lat,p*) for p in points]

Я почти уверен, что в numpy/scipy есть лучшее (и более быстрое) решение. Я уже довольно много гуглил, но ответ до сих пор ускользал от меня.

Есть предложения по эффективному поиску индексов соответствующих (ближайших) точек?

PS: этот вопрос возник из другого (нажмите).

Изменить: Решение

Основываясь на ответе @Cong Ma, я нашел следующее решение:

def find_indices(points,lon,lat,tree=None):
    if tree is None:
        lon,lat = lon.T,lat.T
        lonlat = np.column_stack((lon.ravel(),lat.ravel()))
        tree = sp.spatial.cKDTree(lonlat)
    dist,idx = tree.query(points,k=1)
    ind = np.column_stack(np.unravel_index(idx,lon.shape))
    return [(i,j) for i,j in ind]

Чтобы представить это решение, а также решение из ответа Дивакара, вот некоторые тайминги функции, в которой я использую find_indices (и где это узкое место с точки зрения скорости) (см. ссылку выше):

spatial_contour_frequency/pil0                :   331.9553
spatial_contour_frequency/pil1                :   104.5771
spatial_contour_frequency/pil2                :     2.3629
spatial_contour_frequency/pil3                :     0.3287

pil0 — мой первоначальный подход, pil1 — Divakar, а pil2/pil3 — окончательное решение выше, где дерево создается «на лету» в pil2 (т. е. для каждой итерации цикла, в котором вызывается find_indices) и только один раз в pil3 ( подробности в другой теме). Несмотря на то, что усовершенствование Divakar моего первоначального подхода дает мне 3-кратное ускорение, cKDTree выводит это на совершенно новый уровень с еще 50-кратным ускорением! А перенос создания дерева из функции делает работу еще быстрее.


python numpy scipy performance interpolation
person flotzilla    schedule 02.10.2015    source источник
comment
В вашем скрипте вы создаете новое дерево при каждом вызове find_indices. Если ваша сетка фиксирована между вызовами, вы тратите время на построение одного и того же дерева снова и снова. На самом деле построение этого дерева является самым медленным вызовом этой функции.   -  person Cong Ma    schedule 05.10.2015
comment
Да, я заметил, это то, над чем я сейчас работаю. ;) Я обновлю ответ соответственно. Спасибо за замечание.   -  person flotzilla    schedule 05.10.2015

Ответы (2)


arrow_upward
4
arrow_downward

Если точки достаточно локализованы, вы можете попробовать непосредственно реализацию scipy.spatial cKDTree, как обсуждалось мной в другом сообщении. Этот пост был об интерполяции, но вы можете игнорировать это и просто использовать часть запроса.

тл;др версия:

Прочтите документацию по scipy.sptial.cKDTree< /а>. Создайте дерево, передав (n, m)-образный объект numpy ndarray в инициализатор, и дерево будет создано из n m-мерных координат.

tree = scipy.spatial.cKDTree(array_of_coordinates)

После этого используйте tree.query() для получения k-го ближайшего соседа (возможно, с аппроксимацией и распараллеливанием, см. документы) или используйте tree.query_ball_point() для поиска всех соседей в пределах заданного допуска расстояния.

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

person Cong Ma    schedule 02.10.2015

arrow_upward
1
arrow_downward

Вот общий векторизованный подход с использованием scipy.spatial.distance.cdist. -

import scipy

# Stack lon and lat arrays as columns to form a Nx2 array, where is N is grid**2
lonlat = np.column_stack((lon.ravel(),lat.ravel()))

# Get the distances and get the argmin across the entire N length
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)

# Get the indices corresponding to grid's shape as the final output
ind = np.column_stack((np.unravel_index(idx,lon.shape))).tolist()

Пробный запуск -

In [161]: lon
Out[161]: 
array([[-11.   ,  -7.82 ,  -4.52 ,  -1.18 ,   2.19 ],
       [-12.   ,  -8.65 ,  -5.21 ,  -1.71 ,   1.81 ],
       [-13.   ,  -9.53 ,  -5.94 ,  -2.29 ,   1.41 ],
       [-14.1  ,  -0.04 ,  -6.74 ,  -2.91 ,   0.976]])

In [162]: lat
Out[162]: 
array([[-11.2  ,  -7.82 ,  -4.51 ,  -1.18 ,   2.19 ],
       [-12.   ,  -8.63 ,  -5.27 ,  -1.71 ,   1.81 ],
       [-13.2  ,  -9.52 ,  -5.96 ,  -2.29 ,   1.41 ],
       [-14.3  ,  -0.06 ,  -6.75 ,  -2.91 ,   0.973]])

In [163]: lonlat = np.column_stack((lon.ravel(),lat.ravel()))

In [164]: idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)

In [165]: np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
Out[165]: [[0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [3, 3]]

Тесты времени выполнения -

Определить функции:

def find_indices(lon,lat,x,y):
    lonlat = np.dstack([lon,lat])
    delta = np.abs(lonlat-[x,y])
    ij_1d = np.linalg.norm(delta,axis=2).argmin()
    i,j = np.unravel_index(ij_1d,lon.shape)
    return i,j

def loopy_app(lon,lat,pts):
    return [find_indices(lon,lat,pts[i,0],pts[i,1]) for i in range(pts.shape[0])]

def vectorized_app(lon,lat,points):
    lonlat = np.column_stack((lon.ravel(),lat.ravel()))
    idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
    return np.column_stack((np.unravel_index(idx,lon.shape))).tolist()

Тайминги:

In [179]: lon = np.random.rand(100,100)

In [180]: lat = np.random.rand(100,100)

In [181]: points = np.random.rand(50,2)

In [182]: %timeit loopy_app(lon,lat,points)
10 loops, best of 3: 47 ms per loop

In [183]: %timeit vectorized_app(lon,lat,points)
10 loops, best of 3: 16.6 ms per loop

Для повышения производительности можно использовать np.concatenate место np.column_stack.

person Divakar    schedule 03.10.2015