加速处理 500 万行坐标数据

Posted 2023-03-05

【中文标题】加速处理 500 万行坐标数据

【英文标题】：speeding up processing 5 million rows of coordinate data

【发布时间】：2016-03-14 02:00:39

【问题描述】：

我有一个包含两列（纬度、经度）的 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 文件或其他内容），这样内部循环就不必查看所有其他点，但是我必须弄清楚如何制作确保每个部分的边界都与其相邻部分的边界进行了检查，这似乎过于复杂，恐怕这比它的价值更令人头疼。

那么有什么方法可以加快速度吗？

【问题讨论】：

我不是地理专家，但我最初的直觉反应是通过消除目标周围 5x5 平方英里之外的点来对数据进行分类。这应该可以通过计算 5 英里是多少弧秒并简单地比较数字来实现....

你能给我们一个样本数据集吗？

这里的低调是如何使这不是 O(n) + O(n^2) 操作...我还应该指出，这种事情是旅行推销员问题。

【参考方案1】：

让我们看看你在做什么。

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

现在，你能告诉我列表是否排序了吗？因为如果它已排序，我们肯定想知道这一点。排序是很有价值的信息，尤其是当您处理 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 英里内是否有一个点。让我们考虑一个 10 平方英里的正方形，以您的第 1 点为中心。这是一个以 (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）。然后，我快速检查我们是否已经知道这一点与另一点接近。

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

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

未排序列表

但是如果列表没有排序呢？

@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 链接提供了一种算法）。

建议

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

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

【讨论】：

基本上，我假设计算“真实”距离的成本很高。所以我在中心点周围定义了一个“盒子”并测试另一个点是否在“盒子”内。如果是这样，那么值得计算更昂贵的真实距离，看看两者是否真的在 5 英里范围内。

请记住，即使在外循环中，您也会这样做 500 万次。简单添加的目的是避免函数调用以及幕后发生的任何其他事情。看看你是否可以计算出一个数值来添加 - 这将是你可以合理进行的最快速度。如果愚蠢的比较通过，您可以进行“正式”比较。

一点谷歌搜索将我带到Wikipedia，它有一个很好的尺寸表。可悲的是，它们以公里为单位，但这只是乘法。更糟糕的是，它们的大小各不相同。似乎一个很好的近似值是度数约为 111 公里减去纬度。 （因此，纬度 0 为 111 公里，纬度 80 为 31 公里。不精确，但计算简单且方向正确。）

我还发现了一篇文章，声称人口> 500的最北端的常驻村刚好在78度以上。更多的谷歌搜索导致***上的“最南端城市”文章表明向南 55 度是一个很好的边界。同样，您似乎不必太担心极端情况。

当然，纬度的大小是相当恒定的。我认为大约是每分钟 1 英里，所以 6 英里将是 6 分钟，或 1/10 度

【参考方案2】：

正如Python calculate lots of distances quickly 的回答所指出的，这是 k-D 树的经典用例。

另一种方法是使用扫线算法，如How do I match similar coordinates using Python?的答案所示

这是适用于您的问题的扫描线算法。在我的笔记本电脑上，运行 5M 随机点需要

import itertools as it
import operator as op
import sortedcontainers     # handy library on Pypi
import time

from collections import namedtuple
from math import cos, degrees, pi, radians, sqrt
from random import sample, uniform

Point = namedtuple("Point", "lat long has_close_neighbor")

miles_per_degree = 69

number_of_points = 5000000
data = [Point(uniform( -88.0,  88.0),     # lat
              uniform(-180.0, 180.0),     # long
              True
             )
        for _ in range(number_of_points)
       ]

start = time.time()
# Note: lat is first in Point, so data is sorted by .lat then .long.
data.sort()

print(time.time() - start)

# Parameter that determines the size of a sliding lattitude window
# and therefore how close two points need to be to be to get flagged.
threshold = 5.0  # miles
lat_span = threshold / miles_per_degree
coarse_threshold = (.98 * threshold)**2

# Sliding lattitude window.  Within the window, observations are
# ordered by longitude.
window = sortedcontainers.SortedListWithKey(key=op.attrgetter('long'))

# lag_pt is the 'southernmost' point within the sliding window.
point = iter(data)
lag_pt = next(point)

milepost = len(data)//10

# lead_pt is the 'northernmost' point in the sliding window.
for i, lead_pt in enumerate(data):
    if i == milepost:
        print('.', end=' ')
        milepost += len(data)//10

    # Dec of lead_obs represents the leading edge of window.
    window.add(lead_pt)

    # Remove observations further than the trailing edge of window.
    while lead_pt.lat - lag_pt.lat > lat_span:
        window.discard(lag_pt)
        lag_pt = next(point)

    # Calculate 'east-west' width of window_size at dec of lead_obs
    long_span = lat_span / cos(radians(lead_pt.lat))
    east_long = lead_pt.long + long_span
    west_long = lead_pt.long - long_span

    # Check all observations in the sliding window within
    # long_span of lead_pt.
    for other_pt in window.irange_key(west_long, east_long):

        if other_pt != lead_pt:
            # lead_pt is at the top center of a box 2 * long_span wide by 
            # 1 * long_span tall.  other_pt is is in that box. If desired, 
            # put additional fine-grained 'closeness' tests here. 

            # coarse check if any pts within 80% of threshold distance
            # then don't need to check distance to any more neighbors
            average_lat = (other_pt.lat + lead_pt.lat) / 2
            delta_lat   = other_pt.lat - lead_pt.lat
            delta_long  = (other_pt.long - lead_pt.long)/cos(radians(average_lat))

            if delta_lat**2 + delta_long**2 <= coarse_threshold:
                break

            # put vincenty test here
            #if 0 < vincenty(lead_pt, other_pt).miles <= close_limit:
            #    break

    else:
        data[i] = data[i]._replace(has_close_neighbor=False)

print()      
print(time.time() - start)

【讨论】：

这是一个很好的技术答案，但我恐怕有点过头了！不过，我一定会尽可能地研究它并在我的努力中使用它，谢谢

【参考方案3】：

如果您按纬度 (n log(n)) 对列表进行排序，并且这些点大致均匀分布，则每个点将在 5 英里范围内减少到大约 1000 个点（餐巾纸数学，不准确）。通过仅查看纬度附近的点，运行时间从 n^2 变为 n*log(n)+.0004n^2。希望这可以加快速度。

【讨论】：

再想一想，我把它归结为 nlog(n)+nlog(.0004n)+.00000144n^2 （给予或接受）通过分类到 5 英里纬度 (n*log(n)) 的桶中，然后按纬度对每个桶进行排序。然后对于每个 bin，您比较 9 个相邻 bin 中的点。对于 500000 点，这会导致与我的第一次优化类似的加速，使其达到约 500000 倍。玩得开心。

“直到 2016 年 3 月 15 日，他才意识到他真的应该在数学课上更加努力......”

对它背后的数字的解释真的很好，谢谢你

没问题（我的女朋友在我做橡皮鸭编程时问了一个非常好的问题，我只是做了数学看看它是如何工作的。）

这里是这些不同方程含义的可视化（x 是点数 y 是时间单位）desmos.com/calculator/uz9ll1vrfw

【参考方案4】：

我会试试pandas。 Pandas 是为高效处理大量数据而设计的。无论如何，这可能有助于提高 csv 部分的效率。但从它的声音来看，你已经给自己解决了一个本质上效率低下的问题。您取第 1 点并将其与其他 4,999,999 点进行比较。然后取第 2 点并将其与其他 4,999,998 点进行比较，依此类推。算一算。这就是您正在进行的 12.5 万亿次 比较。如果您每秒可以进行 1,000,000 次比较，那就是 144 天的计算。如果您每秒可以进行 10,000,000 次比较，那就是 14 天。对于直接 python 中的添加，10,000,000 次操作可能需要大约 1.1 秒，但我怀疑您的比较与添加操作一样快。所以至少给它一两个星期。

或者，你可以想出一个替代算法，虽然我没有任何特别的想法。

【讨论】：

我想知道...我想知道是否可以对某些数据子集使用 K-means 聚类算法，然后使用该学习器对其余数据进行分类？一篇有趣的论文可能会沿着这条道路开始：pdfs.semanticscholar.org/2837/…

【参考方案5】：

我会分三步重做算法：

使用大圆距离，假设误差为 1%，因此让 limit 等于 1.01*limit。

将大圆距离编码为内联函数，这个测试应该很快

你会得到一些误报，你可以用 vincenty 进一步测试

【参考方案6】：

由 Oscar Smith 生成的更好的解决方案。你有一个 csv 文件，只是在 excel 中对其进行了排序，它非常有效）。然后在你的程序中使用二分搜索来查找 5 英里范围内的城市（你可以对二分搜索方法做一些小的改动，这样如果它找到一个满足你条件的城市就会中断）。 另一个改进是设置地图以在您发现一个城市在另一个城市内时记住这对城市。例如，当您发现城市 A 在城市 B 5 英里范围内时，使用 Map 存储该对（B 是键，A 是值）。所以下次遇到B，先在Map里搜索一下，如果有对应的值，就不用再查了。但它可能会使用更多的内存，所以要关心它。希望对你有帮助。

【参考方案7】：

这只是第一次通过，但到目前为止，通过使用great_circle() 而不是vincinty()，我已经将它加快了一半，并清理了其他一些东西。区别解释here，准确率损失约0.17%：

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


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


def isCloseToAnother(pointA, points):
    for pointB in points:
        dist = great_circle(pointA, pointB).miles
        if 0 < dist <= CLOSE_LIMIT:  # (0, close_limit]
            return True

    return False


    with open('geo_input.csv', 'r') as geo_csv:
        reader = csv.reader(geo_csv)
        next(reader, None)  # skip the headers

        geo_points = sorted(map(lambda x: CustomGeoPoint(x[0], x[1]), reader))

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

        # for every point, look at every point until one is found within a mile
        for point in geo_points:
            point.close_to_another_point = isCloseToAnother(point, geo_points)
            writer.writerow([point.latitude, point.longitude,
                             point.close_to_another_point])

我会进一步改进。

之前：

$ time python geo.py

real    0m5.765s
user    0m5.675s
sys     0m0.048s

之后：

$ time python geo.py

real    0m2.816s
user    0m2.716s
sys     0m0.041s

【讨论】：

在您返回并修复错误的 0.2% 后，时间是如何影响的？如果它需要完美，我可以想象验证需要一段时间。

抱歉，我的意思是每个距离的准确度都会降低0.17%。

【参考方案8】：

这个问题可以通过VP tree 解决。这些允许查询数据 距离是服从三角不等式的度量。

VP 树相对于 k-D 树的最大优势在于它们可以盲目 适用于世界任何地方的地理数据，无需担心 关于将其投影到合适的二维空间。此外，真正的测地线 可以使用距离（无需担心之间的差异 测地线距离和投影中的距离）。

这是我的测试：在 世界。将它们放入 VP 树中。

遍历所有点，查询VP树以找到任何邻居a 距离 in (0km, 10km] away. (0km 不包括在这个集合中以避免 找到的查询点。）计算没有这样的点的数量 邻居（在我的例子中是 229573）。

建立 VP 树的成本 = 5000000 * 20 距离计算。

查询成本 = 5000000 * 23 次距离计算。

设置和查询时间为 5m 7s。

我使用 C++ 和 GeographicLib 来计算距离，但是 该算法当然可以用任何语言实现，这里是 python version of GeographicLib。

附录：实现此方法的 C++ 代码给出了here。