Рассчитать двумерное попарное расстояние на большом многомерном трехмерном массиве

Question

У меня есть массив numpy 3 миллиона точек в виде [pt_id, x, y, z]. Цель состоит в том, чтобы вернуть все пары точек, у которых евклидово расстояние два числа min_d и max_d.

Евклидово расстояние между x и y, а не на z. Тем не менее, я хотел бы сохранить массив с атрибутами pt_id_from, pt_id_to, distance.

Я использую расст SciPy для вычисления расстояния:

import scipy.spatial.distance 
coords_arr = np.array([['pt1', 2452130.000, 7278106.000, 25.000], 
         ['pt2', 2479539.000, 7287455.000, 4.900], 
         ['pt3', 2479626.000, 7287458.000, 10.000], 
         ['pt4', 2484097.000, 7292784.000, 8.800], 
         ['pt5', 2484106.000, 7293079.000, 7.300], 
         ['pt6', 2484095.000, 7292891.000, 11.100]]) 

dists = scipy.spatial.distance.pdist(coords_arr[:,1:3], 'euclidean') 
np.savetxt('test.out', scipy.spatial.distance.squareform(dists), delimiter=',') 

Что я должен сделать, чтобы вернуть массив формы: [pt_id_from, pt_id_to, distance]?

Answer 1

Ну, ['pt1', 'pt2', distance_as_number] не совсем возможно. Наиболее близким к смешанным типам данных является структурированный массив, но тогда вы не можете делать такие вещи, как result[:2,0]. Вам нужно будет индексировать имена полей и индексы массивов отдельно, например: result[['a','b']][0].

Вот мое решение:

import numpy as np 
import scipy.spatial.distance 

coords_arr = np.array([['pt1', 2452130.000, 7278106.000, 25.000], 
         ['pt2', 2479539.000, 7287455.000, 4.900], 
         ['pt3', 2479626.000, 7287458.000, 10.000], 
         ['pt4', 2484097.000, 7292784.000, 8.800], 
         ['pt5', 2484106.000, 7293079.000, 7.300], 
         ['pt6', 2484095.000, 7292891.000, 11.100]]) 

dists = scipy.spatial.distance.pdist(coords_arr[:,1:3], 'euclidean') 

# Create a shortcut for `coords_arr.shape[0]` which is basically 
# the total amount of points, hence `n` 
n = coords_arr.shape[0] 

# `a` and `b` contain the indices of the points which were used to compute the 
# distances in dists. In this example: 
# a = [0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4] 
# b = [1, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 5, 4, 5, 5] 
a = np.arange(n).repeat(np.arange(n-1, -1, -1)) 
b = np.hstack([range(x, n) for x in xrange(1, n)]) 

min_d = 1000 
max_d = 10000 

# Find out which distances are in range. 
in_range = np.less_equal(min_d, dists) & np.less_equal(dists, max_d) 

# Define the datatype of the structured array which will be the result. 
dtype = [('a', '<f8', (3,)), ('b', '<f8', (3,)), ('dist', '<f8')] 

# Create an empty array. We fill it later because it makes the code cleaner. 
# Its size is given by the sum over `in_range` which is possible 
# since True and False are equivalent to 1 and 0. 
result = np.empty(np.sum(in_range), dtype=dtype) 

# Fill the resulting array. 
result['a'] = coords_arr[a[in_range], 1:4] 
result['b'] = coords_arr[b[in_range], 1:4] 
result['dist'] = dists[in_range] 

print(result) 

# In caste you don't want a structured array at all, this is what you can do: 
result = np.hstack([coords_arr[a[in_range],1:], 
        coords_arr[b[in_range],1:], 
        dists[in_range, None]]).astype('<f8') 
print(result) 

Структурированный массив:

[([2479539.0, 7287455.0, 4.9], [2484097.0, 7292784.0, 8.8], 7012.389393067102) 
([2479539.0, 7287455.0, 4.9], [2484106.0, 7293079.0, 7.3], 7244.7819152821985) 
([2479539.0, 7287455.0, 4.9], [2484095.0, 7292891.0, 11.1], 7092.75912462844) 
([2479626.0, 7287458.0, 10.0], [2484097.0, 7292784.0, 8.8], 6953.856268287403) 
([2479626.0, 7287458.0, 10.0], [2484106.0, 7293079.0, 7.3], 7187.909362255481) 
([2479626.0, 7287458.0, 10.0], [2484095.0, 7292891.0, 11.1], 7034.873843929257)] 

ndarray:

[[2479539.0, 7287455.0, 4.9, 2484097.0, 7292784.0, 8.8, 7012.3893], 
[2479539.0, 7287455.0, 4.9, 2484106.0, 7293079.0, 7.3, 7244.7819], 
[2479539.0, 7287455.0, 4.9, 2484095.0, 7292891.0, 11.1, 7092.7591], 
[2479626.0, 7287458.0, 10.0, 2484097.0, 7292784.0, 8.8, 6953.8562], 
[2479626.0, 7287458.0, 10.0, 2484106.0, 7293079.0, 7.3, 7187.9093], 
[2479626.0, 7287458.0, 10.0, 2484095.0, 7292891.0, 11.1, 7034.8738]] 

Answer 2

Вы можете использовать np.where, чтобы получить координаты расстояний в пределах диапазона, а затем сгенерировать новый список в вашем формате, фильтруя одинаковые пары. Например:

>>> import scipy.spatial.distance 
>>> import numpy as np 
>>> coords_arr = np.array([['pt1', 2452130.000, 7278106.000, 25.000], 
...      ['pt2', 2479539.000, 7287455.000, 4.900], 
...      ['pt3', 2479626.000, 7287458.000, 10.000], 
...      ['pt4', 2484097.000, 7292784.000, 8.800], 
...      ['pt5', 2484106.000, 7293079.000, 7.300], 
...      ['pt6', 2484095.000, 7292891.000, 11.100]]) 
>>> 
>>> dists = scipy.spatial.distance.pdist(coords_arr[:,1:3], 'euclidean') 
>>> dists = scipy.spatial.distance.squareform(dists) 
>>> x, y = np.where((dists >= 8000) & (dists <= 30000)) 
>>> [(coords_arr[x[i]][0], coords_arr[y[i]][0], dists[y[i]][x[i]]) for i in xrange(len(x)) if x[i] < y[i]] 
[('pt1', 'pt2', 28959.576688895162), ('pt1', 'pt3', 29042.897927032005)] 

Answer 3

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

n = coords_arr.shape[0] # number of points 
D = scipy.spatial.distance.squareform(dists) # distance matrix 

data = [] 
for i, j in itertools.combinations(range(n), 2): 
    pt_a = coords_arr[i, 0] 
    pt_b = coords_arr[j, 0] 
    d_ab = D[i,j] 
    data.append([pt_a, pt_b, d_ab]) 

result_arr = np.array(data) 

Если память является проблемой, вы можете захотеть изменить расстояние поиск с помощью огромной матрицы D в поиске значения непосредственно в dists с использованием индекса i и j.