寻求有关本地机器或集群(Python，R，JavaScript，任何语言(算法的帮助。
我有一个带有坐标的位置列表。
# R script
n <- 10
set.seed(1)
index <- paste0("id_",c(1:n))
lat <- runif(n, 32.0, 41)
lon <- runif(n, 84, 112)*(-1)
values <- as.integer(runif(n, 50, 100))
df <- data.frame(index, lat, lon, values, stringsAsFactors = FALSE)
names(df) <- c('loc_id','lat','lon', 'value')
loc_id lat lon value
1 id_1 34.38958 -89.76729 96
2 id_2 35.34912 -88.94359 60
3 id_3 37.15568 -103.23664 82
4 id_4 40.17387 -94.75490 56
5 id_5 33.81514 -105.55556 63
6 id_6 40.08551 -97.93558 69
7 id_7 40.50208 -104.09332 50
8 id_8 37.94718 -111.77337 69
9 id_9 37.66203 -94.64099 93
10 id_10 32.55608 -105.76847 67
我需要为表格中的每个位置找到 3 个壁橱位置。
这是我在 R 中的代码：
# R script
require(dplyr)
require(geosphere)
start.time <- Sys.time()
d1 <- df
sample <- 999999999999
distances <- list("init1" = sample, "init2" = sample, "init3" = sample)
d1$distances <- apply(d1, 1, function(x){distances})
n_rows = nrow(d1)
for (i in 1:(n_rows-1)) {
# current location
dot1 <- c(d1$lon[i], d1$lat[i])
for (k in (i+1):n_rows) {
# next location
dot2 <- c(d1$lon[k], d1$lat[k])
# distance between locations
meters_between <- as.integer(distm(dot1, dot2, fun = distHaversine))
# updating current location distances
distances <- d1$distances[[i]]
distances[d1$loc_id[k]] <- meters_between
d1$distances[[i]] <- distances[order(unlist(distances), decreasing=FALSE)][1:3]
# updating next location distances
distances <- d1$distances[[k]]
distances[d1$loc_id[i]] <- meters_between
d1$distances[[k]] <- distances[order(unlist(distances), decreasing=FALSE)][1:3]
}
}
但这需要太多时间：
# [1] "For 10 rows and 45 iterations takes 0.124729156494141 sec. Average sec 0.00277175903320313 per row."
# [1] "For 100 rows and 4950 iterations takes 2.54944682121277 sec. Average sec 0.000515039761861165 per row."
# [1] "For 200 rows and 19900 iterations takes 10.1178169250488 sec. Average sec 0.000508433011308986 per row."
# [1] "For 500 rows and 124750 iterations takes 73.7151870727539 sec. Average sec 0.000590903303188408 per row."
我在 Python 中做了同样的事情：
# Python script
import pandas as pd
import numpy as np
n = 10
np.random.seed(1)
data_m = np.random.uniform(0, 5, 5)
data = {'loc_id':range(1, n+1),
'lat':np.random.uniform(32, 41, n),
'lon':np.random.uniform(84, 112, n)*(-1),
'values':np.random.randint(50, 100, n)}
df = pd.DataFrame(data)[['loc_id', 'lat', 'lon', 'values']]
df['loc_id'] = df['loc_id'].apply(lambda x: 'id_{0}'.format(x))
df = df.reset_index().drop('index', axis = 1).set_index('loc_id')
from geopy.distance import distance
from datetime import datetime
start_time = datetime.now()
sample = 999999999999
df['distances'] = np.nan
df['distances'] = df['distances'].apply(lambda x: [{'init1': sample}, {'init2': sample}, {'init3': sample}])
n_rows = len(df)
rows_done = 0
for i, row_i in df.head(n_rows-1).iterrows():
dot1 = (row_i['lat'], row_i['lon'])
rows_done = rows_done + 1
for k, row_k in df.tail(n_rows-rows_done).iterrows():
dot2 = (row_k['lat'], row_k['lon'])
meters_between = int(distance(dot1,dot2).meters)
distances = df.at[i, 'distances']
distances.append({k: meters_between})
distances_sorted = sorted(distances, key=lambda x: x[next(iter(x))])[:3]
df.at[i, 'distances'] = distances_sorted
distances = df.at[k, 'distances']
distances.append({i: meters_between})
distances_sorted = sorted(distances, key=lambda x: x[next(iter(x))])[:3]
df.at[k, 'distances'] = distances_sorted
print df
几乎相同的性能。
有人知道是否有更好的方法吗？在我的任务中，必须为90000个位置完成。甚至想到了Hadoop/MpRc/Spark，但不知道如何在分布式模式下做。
我很高兴听到任何想法或建议。
如果欧几里得距离没问题，那么nn2使用 kd 树和 C 代码，所以它应该很快：
library(RANN)
nn2(df[2:3], k = 4)
在我的不是特别快的笔记本电脑上，处理 n = 10,000 行总共需要 0.06 到 0.11 秒，90,000 行总共需要 1.00 到 1.25 秒。
我可以通过scipy提供python解决方案
from scipy.spatial import distance
from geopy.distance import vincenty
v=distance.cdist(df[['lat','lon']].values,df[['lat','lon']].values,lambda u, v: vincenty(u, v).kilometers)
np.sort(v,axis=1)[:,1:4]
Out[1033]:
array([[384.09948155, 468.15944729, 545.41393271],
[270.07677993, 397.21974571, 659.96238603],
[384.09948155, 397.21974571, 619.616239 ],
[203.07302273, 483.54687912, 741.21396029],
[203.07302273, 444.49156394, 659.96238603],
[437.31308598, 468.15944729, 494.91879983],
[494.91879983, 695.91437812, 697.27399161],
[270.07677993, 444.49156394, 483.54687912],
[530.54946479, 626.29467739, 695.91437812],
[437.31308598, 545.41393271, 697.27399161]])
以下是使用C++ 和我的库解决此问题的方法 GeographicLib(版本 1.47 或更高版本(。 这使用真正的椭球测地线 距离和 有利位置树 以优化对最近邻的搜索。
#include <exception>
#include <vector>
#include <fstream>
#include <string>
#include <GeographicLib/NearestNeighbor.hpp>
#include <GeographicLib/Geodesic.hpp>
using namespace std;
using namespace GeographicLib;
// A structure to hold a geographic coordinate.
struct pos {
string id;
double lat, lon;
pos(const string& _id = "", double _lat = 0, double _lon = 0) :
id(_id), lat(_lat), lon(_lon) {}
};
// A class to compute the distance between 2 positions.
class DistanceCalculator {
private:
Geodesic _geod;
public:
explicit DistanceCalculator(const Geodesic& geod) : _geod(geod) {}
double operator() (const pos& a, const pos& b) const {
double d;
_geod.Inverse(a.lat, a.lon, b.lat, b.lon, d);
if ( !(d >= 0) )
// Catch illegal positions which result in d = NaN
throw GeographicErr("distance doesn't satisfy d >= 0");
return d;
}
};
int main() {
try {
// Read in pts
vector<pos> pts;
string id;
double lat, lon;
{
ifstream is("pts.txt"); // lines of "id lat lon"
if (!is.good())
throw GeographicErr("pts.txt not readable");
while (is >> id >> lon >> lat)
pts.push_back(pos(id, lat, lon));
if (pts.size() == 0)
throw GeographicErr("need at least one location");
}
// Define a distance function object
DistanceCalculator distance(Geodesic::WGS84());
// Create NearestNeighbor object
NearestNeighbor<double, pos, DistanceCalculator>
ptsset(pts, distance);
vector<int> ind;
int n = 3; // Find 3 nearest neighbors
for (unsigned i = 0; i < pts.size(); ++i) {
ptsset.Search(pts, distance, pts[i], ind,
n, numeric_limits<double>::max(),
// exclude the point itself
0.0);
if (ind.size() != n)
throw GeographicErr("unexpected number of results");
cout << pts[i].id;
for (unsigned j = 0; j < ind.size(); ++j)
cout << " " << pts[ind[j]].id;
cout << "n";
}
int setupcost, numsearches, searchcost, mincost, maxcost;
double mean, sd;
ptsset.Statistics(setupcost, numsearches, searchcost,
mincost, maxcost, mean, sd);
long long
totcost = setupcost + searchcost,
exhaustivecost = ((pts.size() - 1) * pts.size())/2;
cerr
<< "Number of distance calculations = " << totcost << "n"
<< "With an exhaustive search = " << exhaustivecost << "n"
<< "Ratio = " << double(totcost) / exhaustivecost << "n"
<< "Efficiency improvement = "
<< 100 * (1 - double(totcost) / exhaustivecost) << "%n";
}
catch (const exception& e) {
cerr << "Caught exception: " << e.what() << "n";
return 1;
}
}
这在一组点(以"id lat lon"的形式(中读取pts.txt， 将它们放在 VP 树中。 然后对于每个点，它查找 3 个最近的点 邻居并打印邻居的 ID 和 ID(排名依据 距离(。
编译它，例如，
g++ -O3 -o nearest nearest.cpp -lGeographic
如果 pts.txt 包含 90000 个点，则计算在 在我家电脑上做大约 6 秒(或每点 70 μs(大约 3380000 距离后 计算。 这 效率约为蛮力计算的 1200 倍 (执行所有N(N− 1(/2 距离计算(。
你可以通过使用粗略的来加快速度("几"的倍数( 距离的近似值(例如，球形或欧几里得(;只 适当地修改距离计算器类。 例如，这个 距离计算器的版本返回球面距离 度：
// A class to compute the spherical distance between 2 positions.
class DistanceCalculator {
public:
explicit DistanceCalculator(const Geodesic& /*geod*/) {}
double operator() (const pos& a, const pos& b) const {
double sphia, cphia, sphib, cphib, somgab, comgab;
Math::sincosd(a.lat, sphia, cphia);
Math::sincosd(b.lat, sphib, cphib);
Math::sincosd(Math::AngDiff(a.lon, b.lon), somgab, comgab);
return Math::atan2d(Math::hypot(cphia * sphib - sphia * cphib * comgab,
cphib * somgab),
sphia * sphib + cphia * cphib * comgab);
}
};
但是现在您有额外的负担，即确保近似 已经足够好了。 我建议只使用正确的测地线距离 首先。
给出了 GeographicLib 中 VP 树实现的详细信息 这里。