加速处理500万行坐标数据

zer*_*pha 5 python csv geopy

我有一个包含两列(纬度,经度)的csv文件,其中包含超过500万行的地理位置数据.我需要识别与列表中任何其他点不在5英里范围内的点,并将所有内容输出回另一个具有额外列的CSV(CloseToAnotherPoint),True如果有另一个点在5英里内,并且False如果有"T.

这是我当前使用的解决方案geopy(不进行任何Web调用,只使用函数来计算距离):

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()
Run Code Online (Sandbox Code Playgroud)

正如您可能从中看到的那样,这种解决方案非常缓慢.实际上这么慢,我让它运行了3天,但仍然没有完成!

我曾经考虑过将数据拆分成块(多个CSV文件或其他东西),这样内部循环就不必查看其他所有点,但是我必须弄清楚如何确保边界每个部分都按照相邻部分的边界进行检查,这看起来过于复杂,我担心它会比它的价值更令人头痛.

那么关于如何加快速度的任何指针呢?

agh*_*ast 4

让我们看看你在做什么。

\n\n
    \n
  1. 您将所有点读入名为 的列表中geo_points。

    \n\n

    现在,你能告诉我列表是否已排序吗?因为如果它被排序,我们肯定想知道。排序是很有价值的信息,尤其是当您处理 500 万个东西时。

  2. \n
  3. 您循环遍历所有geo_points. 据你说,那是 500 万。

  4. \n
  5. 在外循环中,您再次循环所有 500 万个geo_points。

  6. \n
  7. 您计算两个循环项目之间的距离(以英里为单位)。

  8. \n
  9. 如果距离小于阈值,则在第一个点上记录该信息,并停止内循环。

  10. \n
  11. 当内部循环停止时,您将有关外部循环项的信息写入 CSV 文件。

  12. \n
\n\n

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

\n\n

这就是 O(n\xc2\xb2) 的意思。

\n\n

下次当你看到有人谈论“哦,这是 O(log n) 但其他事情是 O(n log n)”时,请记住这个经历 - 你正在运行一个 n\xc2\xb2 算法,其中 n 在这个案例是5,000,000。糟透了,知道吗?

\n\n

无论如何,你有一些问题。

\n\n

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

\n\n

问题 2:当您将点 #1 与点 #12345 进行比较,并且它们彼此之间的距离在阈值范围内时,您正在记录有关点 #1 的信息。但您不会记录另一点的相同信息。你知道点 #12345 (geo_point2) 反射性地位于点 #1 的阈值内,但您没有写下来。所以您错过了跳过超过 500 万次比较的机会。

\n\n

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

\n\n

答案?

\n\n

考虑平面上的纬度和经度。您想知道 5 英里内是否有一个点。让我们考虑一下以第 1 点为中心的 10 英里正方形。这是一个以 (X1, Y1) 为中心的正方形,左上角位于 (X1 - 5 英里,Y1 + 5 英里),右下角位于 (X1 + 5 英里,Y1 - 5 英里)。现在,如果一个点在那个正方形内,它可能不是您的 1 点的 5 英里范围内。但你可以打赌,如果它在那个广场之外,那么距离就超过 5 英里。

\n\n

正如@SeverinPappadeaux 指出的那样,像地球这样的球体上的距离与平面上的距离并不完全相同。但那又怎样呢?将你的方块设置得大一点以允许差异,然后继续!

\n\n

排序列表

\n\n

这就是为什么排序很重要。如果所有点都按 X 排序,然后按 Y(或 Y,然后 X - 无论如何)排序,并且您知道这一点,那么您确实可以加快速度。因为当 X(或 Y)坐标变得太大时,您可以简单地停止扫描,并且不必经过 500 万个点。

\n\n

那会如何运作呢?和以前一样,除了你的内部循环会有一些像这样的检查:

\n\n
five_miles = ... # Whatever math, plus an error allowance!\nlist_len = len(geo_points) # Don\'t call this 5 million times\n\nfor i, pi in enumerate(geo_points):\n\n    if pi.close_to_another_point:\n        continue   # Remember if close to an earlier point\n\n    pi0max = pi[0] + five_miles\n    pi1min = pi[1] - five_miles\n    pi1max = pi[1] + five_miles\n\n    for j in range(i+1, list_len):\n        pj = geo_points[j]\n        # Assumes geo_points is sorted on [0] then [1]\n        if pj[0] > pi0max:\n            # Can\'t possibly be close enough, nor any later points\n            break\n        if pj[1] < pi1min or pj[1] > pi1max:\n            # Can\'t be close enough, but a later point might be\n            continue\n\n        # Now do "real" comparison using accurate functions.\n        if ...:\n            pi.close_to_another_point = True\n            pj.close_to_another_point = True\n            break\n
Run Code Online (Sandbox Code Playgroud)\n\n

我在那里做什么?首先,我将一些数字放入局部变量中。然后我用来enumerate给我一个i值和对外部点的引用。(你所谓的geo_point)。然后,我快速检查我们是否已经知道这一点与另一点接近。

\n\n

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

\n\n

最后,如果简单的检查通过了,那么就继续进行昂贵的检查。如果检查确实通过了,请务必记录这两点的结果,以便我们稍后可以跳过第二点。

\n\n

未排序列表

\n\n

但如果列表未排序怎么办?

\n\n

@RootTwo 将您指向 kD 树(其中 D 代表“维度”,k 在本例中为“2”)。如果您已经了解二叉搜索树,那么这个想法非常简单:循环遍历维度,比较树中偶数级别的 X 并比较奇数级别的 Y(反之亦然)。这个想法是这样的:

\n\n
def insert_node(node, treenode, depth=0):\n    dimension = depth % 2  # even/odd -> lat/long\n    dn = node.coord[dimension]\n    dt = treenode.coord[dimension]\n\n    if dn < dt:\n        # go left\n        if treenode.left is None:\n            treenode.left = node\n        else:\n            insert_node(node, treenode.left, depth+1)\n    else:\n        # go right\n        if treenode.right is None:\n            treenode.right = node\n        else:\n            insert_node(node, treenode.right, depth+1)\n
Run Code Online (Sandbox Code Playgroud)\n\n

这会做什么?这将为您提供一个可搜索的树,其中可以在 O(log n) 时间内插入点。这意味着整个列表的复杂度为 O(n log n),这比 n 的平方要好得多!(500 万的对数底数 2 基本上是 23。所以 n log n 是 500 万乘以 23,而 500 万乘以 500 万!)

\n\n

这也意味着您可以进行有针对性的搜索。由于树是有序的,因此寻找“接近”点相当简单(@RootTwo 的维基百科链接提供了一种算法)。

\n\n

建议

\n\n

我的建议是,如果需要的话,只需编写代码来对列表进行排序。它更容易书写,也更容易手动检查,而且它是一张单独的通行证,您只需进行一次。

\n\n

对列表进行排序后,请尝试我上面展示的方法。它与您正在做的事情很接近,并且应该很容易让您理解和编码。

\n