Векторизация для вычисления многих расстояний

Я новичок в numpy / pandas и векторизованном вычислении. Я выполняю задачу данных, где у меня есть два набора данных. Набор данных 1 содержит список мест с их долготой и широтой и переменную A. В наборе данных 2 также содержится список мест с их долготой и широтой. Для каждого места в наборе данных 1 я хотел бы рассчитать его расстояния ко всем местам в наборе данных 2, но я хотел бы получить только количество мест в наборе данных 2, которые меньше, чем значение переменной A. Обратите внимание, что оба наборы данных очень большие, поэтому мне нужно использовать векторизованные операции для ускорения вычисления.

Например, мой набор данных1 может выглядеть следующим образом:

id lon lat varA 1 20.11 19.88 100 2 20.87 18.65 90 3 18.99 20.75 120 

и мой dataset2 может выглядеть следующим образом:

 placeid lon lat a 18.75 20.77 b 19.77 22.56 c 20.86 23.76 d 17.55 20.74 

Тогда для id == 1 в наборе данных1, я хотел бы рассчитать его расстояния ко всем четырем точкам (a, c, c, d) в наборе данных2, и я хотел бы подсчитать, сколько расстояний меньше соответствующих значение varA. Например, рассчитанные четыре расстояния составляют 90, 70, 120, 110, а varA равно 100. Затем значение должно быть 2.

У меня уже есть векторизованная функция для вычисления расстояния между двумя парами координат. Предположим, что функция (haversine (x, y)) реализована правильно, у меня есть следующий код.

 dataset2['count'] = dataset1.apply(lambda x: haversine(x['lon'],x['lat'],dataset2['lon'], dataset2['lat']).shape[0], axis = 1) 

Однако это дает общее количество строк, но не те, которые удовлетворяют моим требованиям.

Кто-нибудь сможет указать мне, как заставить код работать?

3 Solutions collect form web for “Векторизация для вычисления многих расстояний”

Если вы можете проецировать координаты на локальную проекцию (например, UTM ), которая довольно прямолинейна с pyproj и, как правило, более благоприятна, чем lon / lat для измерения, тогда существует намного более быстрый способ использования scipy.spatial . Ни один из df['something'] = df.apply(...) и np.vectorize() не являются действительно векторизованными, под капотом они используют цикл.

 ds1 id lon lat varA 0 1 20.11 19.88 100 1 2 20.87 18.65 90 2 3 18.99 20.75 120 ds2 placeid lon lat 0 a 18.75 20.77 1 b 19.77 22.56 2 c 20.86 23.76 3 d 17.55 20.74 from scipy.spatial import distance # gey coordinates of each set of points as numpy array coords_a = ds1.values[:,(1,2)] coords_b = ds2.values[:, (1,2)] coords_a #out: array([[ 20.11, 19.88], # [ 20.87, 18.65], # [ 18.99, 20.75]]) distances = distance.cdist(coords_a, coords_b) #out: array([[ 1.62533074, 2.70148108, 3.95182236, 2.70059253], # [ 2.99813275, 4.06178532, 5.11000978, 3.92307278], # [ 0.24083189, 1.97091349, 3.54358575, 1.44003472]]) 

distances фактически являются расстоянием между каждой парой точек. coords_a.shape(3, 2) а coords_b.shape(4, 2) , поэтому результат равен (3,4) . Метрика по умолчанию для np.distance является eculidean , но есть и другие показатели. Для этого примера предположим, что vara :

 vara = np.array([2,4.5,2]) 

(вместо 100 90 120 ). Нам нужно определить, какое значение на distances в строке 1 меньше 2 , в строке два меньше, чем 4.5 , …, одним из способов решения этой проблемы является вычитание каждого значения из vara из соответствующей строки (обратите внимание, что мы должны изменить размер vara ) :

 vara.resize(3,1) res = res - vara #out: array([[-0.37466926, 0.70148108, 1.95182236, 0.70059253], # [-1.50186725, -0.43821468, 0.61000978, -0.57692722], # [-1.75916811, -0.02908651, 1.54358575, -0.55996528]]) 

то установление положительных значений на нуль и положительное положительное значение даст нам окончательный массив:

 res[res>0] = 0 res = np.absolute(res) #out: array([[ 0.37466926, 0. , 0. , 0. ], # [ 1.50186725, 0.43821468, 0. , 0.57692722], # [ 1.75916811, 0.02908651, 0. , 0.55996528]]) 

Теперь, чтобы суммировать по каждой строке:

 sum_ = res.sum(axis=1) #out: array([ 0.37466926, 2.51700915, 2.34821989]) 

и подсчитывать элементы в каждой строке:

 count = np.count_nonzero(res, axis=1) #out: array([1, 3, 3]) 

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

 x, y = np.mgrid[0:4, 0:4] points = zip(x.ravel(), y.ravel()) tree = spatial.cKDTree(points) tree.query_ball_point([2, 0], 1) [4, 8, 9, 12] 

query_ball_point() находит все точки на расстоянии r точки (ов) x, и это удивительно быстро.

последнее замечание: не используйте эти алгоритмы с входом lon / lat, особенно если ваша область интересов далека от экватора, потому что ошибка может стать огромной.

ОБНОВИТЬ:

Чтобы спроецировать координаты, вам необходимо преобразовать из WGS84 (lon/lat) в соответствующий UTM . Чтобы узнать, какую utm-зону вам следует проектировать, используйте epsg.io.

 lon = -122.67598 lat = 45.52168 WGS84 = "+init=EPSG:4326" EPSG3740 = "+init=EPSG:3740" Proj_to_EPSG3740 = pyproj.Proj(EPSG3740) Proj_to_EPSG3740(lon,lat) # out: (525304.9265963673, 5040956.147893889) 

Вы можете сделать df.apply() и использовать Proj_to_... для проекта df.

IIUC:

Источник DF:

 In [160]: d1 Out[160]: id lon lat varA 0 1 20.11 19.88 100 1 2 20.87 18.65 90 2 3 18.99 20.75 120 In [161]: d2 Out[161]: placeid lon lat 0 a 18.75 20.77 1 b 19.77 22.56 2 c 20.86 23.76 3 d 17.55 20.74 

haversine функция haversine :

 def haversine(lat1, lon1, lat2, lon2, to_radians=True, earth_radius=6371): if to_radians: lat1, lon1, lat2, lon2 = pd.np.radians([lat1, lon1, lat2, lon2]) a = pd.np.sin((lat2-lat1)/2.0)**2 + \ pd.np.cos(lat1) * pd.np.cos(lat2) * pd.np.sin((lon2-lon1)/2.0)**2 return earth_radius * 2 * pd.np.arcsin(np.sqrt(a)) 

Решение:

 x = d2.assign(x=1) \ .merge(d1.loc[d1['id']==1, ['lat','lon']].assign(x=1), on='x', suffixes=['','2']) \ .drop(['x'], 1) x['dist'] = haversine(x.lat, x.lon, x.lat2, x.lon2) 

выходы:

 In [163]: x Out[163]: placeid lon lat lat2 lon2 dist 0 a 18.75 20.77 19.88 20.11 172.924852 1 b 19.77 22.56 19.88 20.11 300.078600 2 c 20.86 23.76 19.88 20.11 438.324033 3 d 17.55 20.74 19.88 20.11 283.565975 

фильтрация:

 In [164]: x.loc[x.dist < d1.loc[d1['id']==1, 'varA'].iat[0]] Out[164]: Empty DataFrame Columns: [placeid, lon, lat, lat2, lon2, dist] Index: [] 

давайте изменим d1 , поэтому несколько строк удовлетворяют критериям:

 In [171]: d1.loc[0, 'varA'] = 350 In [172]: d1 Out[172]: id lon lat varA 0 1 20.11 19.88 350 # changed: 100 --> 350 1 2 20.87 18.65 90 2 3 18.99 20.75 120 In [173]: x.loc[x.dist < d1.loc[d1['id']==1, 'varA'].iat[0]] Out[173]: placeid lon lat lat2 lon2 dist 0 a 18.75 20.77 19.88 20.11 172.924852 1 b 19.77 22.56 19.88 20.11 300.078600 3 d 17.55 20.74 19.88 20.11 283.565975 

Используйте scipy.spatial.distance.cdist с вашим пользовательским алгоритмом расстояния, как metric

 h = lambda u, v: haversine(u['lon'], u['lat'], v['lon'], v['lat']) dist_mtx = scipy.spatial.distance.cdist(dataset1, dataset2, metric = h) 

Затем, чтобы проверить номер в области, просто передайте его

 dataset2['count'] = np.sum(dataset1['A'][:, None] > dist_mtx, axis = 0) 
 
Interesting Posts for Van-Lav

AttributeError: объект 'module' не имеет атрибута 'request'

Отладка кода Conda и Visual Studio

Почему вы должны избегать использования «is» и «is not»?

Как установить версию Python по умолчанию на 3.3 в OS X?

Как вы даже можете дать вход (FST) FST? Где идет выход?

TextCtrl обеспечивает исключение из исключения в wxPython

установка cPickle с помощью python 3.5

Программа, которая открывает текстовый файл, подсчитывает количество слов и сообщает верхние N слов, упорядоченных по количеству раз, когда они появляются в файле?

как динамически создавать переменную класса в python

Python: тестирование для None, тестирование для логического значения

У Python есть неизменный список?

Разделить строки с несколькими разделителями?

Печать строки unicode в python независимо от среды

pandas: Как получить метод .to_string () для выравнивания заголовков столбцов со значениями столбцов?

Начальная практика Python?

Python - лучший язык программирования в мире.