5

我有一个包含超过 500 万行地理位置数据的两列（纬度、经度）的 csv 文件。我需要识别不在列表中任何其他点 5 英里范围内的点，并将所有内容输出回另一个具有额外列 ( CloseToAnotherPoint) 的 CSV，即True是否有另一个点在 5 英里范围内，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 文件或其他内容），以便内部循环不必查看所有其他点，但随后我必须弄清楚如何确保边界每个部分都对照其相邻部分的边界进行检查，这似乎过于复杂，恐怕这比它的价值更令人头疼。

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

4

8 回答 8

4

正如Python 计算大量距离的答案很快指出的那样，这是 kD 树的经典用例。

另一种方法是使用扫描线算法，如如何使用 Python 匹配相似坐标的答案所示？

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

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)
于 2016-03-14T02:40:23.347 回答
4

让我们看看你在做什么。

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

    现在，你能告诉我列表是否排序了吗？因为如果它被排序，我们肯定想知道。排序是很有价值的信息，尤其是当您处理 500 万件事情时。

  2. 你遍历所有的geo_points. 据你说，那是500万。

  3. 在外部循环中，您再次循环所有 500 万个geo_points。

  4. 您计算两个循环项目之间的距离（以英里为单位）。

  5. 如果距离小于您的阈值，则在第一个点记录该信息，并停止内部循环。

  6. 当内部循环停止时，您将有关外部循环项的信息写入 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)只需让外循环 go和内循环 go ，您就可以将问题空间减少一半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）。然后，我快速检查我们是否已经知道这一点与另一点接近。

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

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

未排序列表

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

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

建议

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

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

于 2016-03-14T04:46:54.440 回答
3

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

于 2016-03-14T02:20:37.070 回答
2

我会试试熊猫。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 秒，但我怀疑您的比较与添加操作一样快。所以至少给它一两个星期。

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

于 2016-03-14T02:32:15.007 回答
1

我将分三步重做算法：

  1. 使用大圆距离，并假设 1% 的误差，所以让 limit 等于 1.01*limit。

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

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

于 2016-03-14T02:20:38.087 回答
1

这只是第一次通过，但到目前为止，我已经通过使用great_circle()而不是加快了一半vincinty()，并清理了其他一些东西。此处解释了差异，准确性损失约为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
于 2016-03-14T03:05:13.013 回答
1

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

于 2016-03-14T03:12:27.257 回答
1

这个问题可以用VP树来解决。这些允许查询距离是遵守三角不等式的度量的数据。

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

这是我的测试：在世界上随机均匀地生成 500 万个点。将这些放入 VP 树中。

循环遍历所有点，查询VP树以找到距离为（0km，10km]的任何邻居。（0km不包括在此集合中，以避免找到查询点。）计算没有这样邻居的点的数量（在我的情况下是 229573）。

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

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

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

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

附录：这里给出了实现这种方法的 C++ 代码。

于 2016-03-27T09:04:40.893 回答