python - 加速处理500万行坐标数据

Question

我有一个包含两列(纬度，经度)的csv文件，其中包含超过500万行地理位置数据。
我需要识别不在列表中任何其他点的5英里之内的点，并将所有内容输出回另一个具有额外列(CloseToAnotherPoint)的CSV中，如果还有另一点位于5英里之内，则为True和False如果没有。

这是我当前使用geopy的解决方案(不进行任何网络调用，仅使用该函数来计算距离):

from geopy.point import Point
from geopy.distance import vincenty
import csv


class CustomGeoPoint(object):
    def __init__(self, latitude, longitude):
        self.location = Point(latitude, longitude)
        self.close_to_another_point = False


try:
    output = open('output.csv','w')
    writer = csv.writer(output, delimiter = ',', quoting=csv.QUOTE_ALL)
    writer.writerow(['Latitude', 'Longitude', 'CloseToAnotherPoint'])

    # 5 miles
    close_limit = 5
    geo_points = []

    with open('geo_input.csv', newline='') as geo_csv:
        reader = csv.reader(geo_csv)
        next(reader, None) # skip the headers
        for row in reader:
            geo_points.append(CustomGeoPoint(row[0], row[1]))

    # for every point, look at every point until one is found within 5 miles
    for geo_point in geo_points:
        for geo_point2 in geo_points:
            dist = vincenty(geo_point.location, geo_point2.location).miles
            if 0 < dist <= close_limit: # (0,close_limit]
                geo_point.close_to_another_point = True
                break
        writer.writerow([geo_point.location.latitude, geo_point.location.longitude,
                         geo_point.close_to_another_point])

finally:
    output.close()

从您可能看到的结果可以看出，该解决方案非常慢。实际上，它是如此之慢，以至于我让它运行了3天，但仍然没有完成!

我曾考虑过尝试将数据分成多个块(多个CSV文件或其他内容)，以使内部循环不必查看所有其他点，但是随后我必须弄清楚如何确保边框每个部分中的每个部分都与相邻部分的边界相对，这似乎太复杂了，我担心这将是一个令人头疼的问题，而不是值得的。

那么关于如何使其更快的任何指示呢？

Answer 1

让我们看看你在做什么。

您将所有要点读入名为geo_points的列表中。

现在，您能告诉我列表是否已排序吗？因为如果排序，我们肯定是想要知道这一点。排序是有值(value)的信息，尤其是当您要处理500万种东西时。

您遍历所有geo_points。根据您的说法，这就是500万。

在外部循环中，您将再次遍历所有500万个geo_points。

您可以计算两个循环项之间的距离(以英里为单位)。

如果距离小于阈值，则在第一个点上记录该信息，然后停止内循环。

当内部循环停止时，将有关外部循环项目的信息写入CSV文件。

注意几件事。首先，您在外循环中循环了500万次。然后，您在内部循环中循环了500万次。

这就是O(n²)的意思。

下次您看到有人在谈论“哦，这是O(log n)，而另一件事是O(n log n)”时，请记住这种体验-您正在运行n²算法，其中n在这种情况下为5,000,000。糟透了，邓尼特？

无论如何，您有一些问题。

问题1:您最终将比较每个点与自身。它们的距离应为零，这意味着它们都将被标记为在任何距离阈值之内。如果您的程序完成，则所有单元格都将标记为True。

问题2:当您将点1与点12345进行比较，并且它们之间的距离都在阈值范围内时，您正在记录有关点1的信息。但是，您不会记录关于另一点的相同信息。您知道点＃12345(geo_point2)反身在点＃1的阈值内，但您不必写下来。因此，您错过了跳过500万次比较的机会。

问题3:如果比较点1和点2，但它们不在阈值距离之内，那么将点2和点1比较时会发生什么？您的内部循环每次都是从列表的开头开始的，但是您知道已经将列表的开始与列表的末尾进行了比较。您只需将外循环设为i in range(0, 5million)而将内循环设为j in range(i+1, 5million)，就可以将问题空间减少一半。

答案？

考虑您在平面上的经度和纬度。您想知道5英里内是否有一个点。让我们考虑一个以1号点为中心的10平方英里的正方形。那是一个以(X1，Y1)为中心的正方形，左上角为(X1-5miles，Y1 + 5miles)，右下角为(X1 + 5miles，Y1-5miles)。现在，如果某个点在该正方形内，则可能不在您的＃1点的5英里之内。但是您可以打赌，如果它在那个广场外，那么它就在5英里之外。

正如@SeverinPappadeaux指出的那样，像地球这样的椭球体上的距离与平面上的距离并不完全相同。但是那又怎样呢？将您的正方形设置得更大一些，以允许差异，然后继续!

排序列表

这就是为什么排序很重要的原因。如果所有点都按X，Y(或Y，然后X-等等)排序，并且您知道了，那么您真的可以加快速度。因为您可以在X(或Y)坐标太大时停止扫描，而不必经过500万个点。

那将如何工作？与之前相同，只是您的内部循环将进行如下检查:

five_miles = ... # Whatever math, plus an error allowance!
list_len = len(geo_points) # Don't call this 5 million times

for i, pi in enumerate(geo_points):

    if pi.close_to_another_point:
        continue   # Remember if close to an earlier point

    pi0max = pi[0] + five_miles
    pi1min = pi[1] - five_miles
    pi1max = pi[1] + five_miles

    for j in range(i+1, list_len):
        pj = geo_points[j]
        # Assumes geo_points is sorted on [0] then [1]
        if pj[0] > pi0max:
            # Can't possibly be close enough, nor any later points
            break
        if pj[1] < pi1min or pj[1] > pi1max:
            # Can't be close enough, but a later point might be
            continue

        # Now do "real" comparison using accurate functions.
        if ...:
            pi.close_to_another_point = True
            pj.close_to_another_point = True
            break

我在那做什么首先，我将一些数字输入局部变量。然后，我使用enumerate给我一个i值和对外部点的引用。 (您称为geo_point)。然后，我迅速检查我们是否已经知道，这一点与另一点很接近。

如果没有，我们将不得不扫描。因此，我只扫描列表中的“较晚”点，因为我知道外部循环会扫描较早的点，并且我绝对不希望将某个点与其自身进行比较。我正在使用一些临时变量来缓存涉及外循环的计算结果。在内部循环中，我对临时对象进行了一些愚蠢的比较。他们无法告诉我这两个点是否彼此接近，但是我可以检查它们是否确实没有接近并且向前跳过。

最后，如果简单检查通过，则继续进行昂贵的检查。如果确实通过了检查，请确保在两个点上都记录结果，因此以后可以跳过第二点。

未排序列表

但是，如果列表未排序怎么办？

@RootTwo将您指向一棵kD树(其中D代表“维度”，而k代表“2”)。如果您已经知道二叉搜索树，这个想法真的很简单:您遍历维度，比较树中偶数级的X和比较奇数级的Y(反之亦然)。这个想法是这样的:

def insert_node(node, treenode, depth=0):
    dimension = depth % 2  # even/odd -> lat/long
    dn = node.coord[dimension]
    dt = treenode.coord[dimension]

    if dn < dt:
        # go left
        if treenode.left is None:
            treenode.left = node
        else:
            insert_node(node, treenode.left, depth+1)
    else:
        # go right
        if treenode.right is None:
            treenode.right = node
        else:
            insert_node(node, treenode.right, depth+1)

这会怎么办？这将为您提供可搜索的树，可以在O(log n)时间中插入点。这意味着整个列表的O(n log n)，比n平方好! (500万的对数基数2基本上是23。因此n log n是500万乘以23，而500万乘以500万!)

这也意味着您可以进行有针对性的搜索。由于树是有序的，因此查找“关闭”点相当简单(@RootTwo的Wikipedia链接提供了一种算法)。

建议

我的建议是，如果需要的话，只写代码对列表进行排序。它更容易编写，也更易于手工检查，并且是一次单独的传递，您只需要进行一次即可。

对列表进行排序后，请尝试上面显示的方法。它接近您所做的事情，应该易于理解和编码。