numpy 数组：快速填充和提取数据

Question

请参阅此问题底部的重要说明。

我正在使用 numpy 来加快一些经度/纬度坐标的处理。不幸的是，我的 numpy“优化”使我的代码运行速度比不使用 numpy 时运行速度快 5 倍。

瓶颈似乎在于用我的数据填充 numpy 数组，然后在我完成数学转换后提取该数据。为了填充数组，我基本上有一个循环，例如：

point_list = GetMyPoints() # returns a long list of ( lon, lat ) coordinate pairs
n = len( point_list )
point_buffer = numpy.empty( ( n, 2 ), numpy.float32 )

for point_index in xrange( 0, n ):
    point_buffer[ point_index ] = point_list[ point_index ]

那个循环，只是在操作之前填充 numpy 数组，非常慢，比没有 numpy 的整个计算要慢得多。 （也就是说，这不仅仅是 python 循环本身的缓慢，而且显然在将每个小数据块从 python 传输到 numpy 时会产生巨大的开销。）在另一端也有类似的缓慢；在处理完 numpy 数组后，我在循环中访问每个修改后的坐标对，再次作为

some_python_tuple = point_buffer[ index ]

再一次，提取数据的循环比没有 numpy 的整个原始计算慢得多。那么，我如何实际填充 numpy 数组并从 numpy 数组中提取数据，而不会破坏首先使用 numpy 的目的？

我正在使用 C 库从形状文件中读取数据，该库将数据作为常规 python 列表传递给我。我知道，如果图书馆将坐标交给我已经在 numpy 数组中，则不需要“填充”numpy 数组。但不幸的是，我使用数据的起点是常规的 python 列表。更重要的是，总的来说，我想了解如何使用 python 中的数据快速填充一个 numpy 数组。

澄清

上面显示的循环实际上过于简单化了。我在这个问题中这样写是因为我想专注于我看到的尝试在循环中缓慢填充 numpy 数组的问题。我现在明白这样做很慢。

在我的实际应用程序中，我拥有的是坐标点的形状文件，并且我有一个 API 来检索给定对象的点。大约有 200,000 个对象。所以我反复调用函数GetShapeCoords( i ) 来获取对象i 的坐标。这将返回一个列表列表，其中每个子列表是 lon/lat 对的列表，它是列表列表的原因是某些对象是多部分的（即多多边形）。然后，在我的原始代码中，当我读取每个对象的点时，我通过调用常规 python 函数对每个点进行转换，然后使用 PIL 绘制转换后的点。整个过程花了大约 20 秒来绘制所有 200,000 个多边形。不可怕，但还有很大的改进空间。我注意到这 20 秒中至少有一半用于转换逻辑，所以我想我会在 numpy.而我最初的实现只是一次读取一个对象，并不断将子列表中的所有点附加到一个大的 numpy 数组中，然后我可以在 numpy 中进行数学运算。

所以，我现在明白，简单地将整个 python 列表传递给 numpy 是设置大数组的正确方法。但就我而言，我一次只读取一个对象。所以我可以做的一件事是在一个列表列表的大 python 列表中继续附加点。然后当我以这种方式编译了大量对象的点（例如，10000 个对象）时，我可以简单地将怪物列表分配给 numpy。

所以我现在的问题是三个部分：

(a) numpy 是不是真的可以把这么大的、不规则形状的列表列表的列表，然后又快又好地吞下去？

(b) 然后我希望能够转换那棵怪物树的叶子中的所有点。例如，让 numpy 的表达式是“进入每个子列表，然后进入每个子子列表，然后对于在这些子子列表中找到的每个坐标对，将第一个（lon 坐标）乘以 0.5”？我可以这样做吗？

(c) 最后，我需要将这些转换后的坐标取回以便绘制它们。

Winston 在下面的回答似乎暗示了我如何使用 itertools 来完成这一切。我想做的很像温斯顿所做的，把列表弄平。但我不能完全把它弄平。当我去绘制数据时，我需要能够知道一个多边形何时停止，下一个多边形何时开始。因此，我认为如果有一种方法可以使用特殊的坐标对（如 (-1000, -1000) 或类似的东西）快速标记每个多边形（即每个子子列表）的结尾，我可以让它工作。然后我可以像 Winston 的回答那样用 itertools 展平，然后在 numpy 中进行转换。然后我需要使用 PIL 从点到点实际绘制，在这里我想我需要将修改后的 numpy 数组重新分配回 python 列表，然后在常规 python 循环中遍历该列表以进行绘图。这似乎是我最好的选择，而不是仅仅编写一个 C 模块来为我一步处理所有的阅读和绘图？

Answer 1

您将您的数据描述为“坐标列表列表”。从这里我猜你的提取看起来像这样：

for x in points:
   for y in x:
       for Z in y:
           # z is a tuple with GPS coordinates

这样做：

# initially, points is a list of lists of lists
points = itertools.chain.from_iterable(points)
# now points is an iterable producing lists
points = itertools.chain.from_iterable(points)
# now points is an iterable producing coordinates
points = itertools.chain.from_iterable(points)
# now points is an iterable producing individual floating points values
data = numpy.fromiter(points, float)
# data is a numpy array containing all the coordinates
data = data.reshape( data.size/2,2)
# data has now been reshaped to be an nx2 array

itertools 和 numpy.fromiter 都是用 c 实现的，而且非常高效。因此，这应该会很快完成转换。

您问题的第二部分并没有真正表明您想对数据做什么。索引 numpy 数组比索引 python 列表慢。您可以通过对数据执行大量操作来获得速度。如果不知道更多关于您正在使用这些数据做什么，很难建议如何修复它。

更新：

我已经使用 itertools 和 numpy 完成了所有工作。我对试图理解此代码造成的任何脑损伤概不负责。

# firstly, we use imap to call GetMyPoints a bunch of times
objects = itertools.imap(GetMyPoints, xrange(100))
# next, we use itertools.chain to flatten it into all of the polygons
polygons = itertools.chain.from_iterable(objects)
# tee gives us two iterators over the polygons
polygons_a, polygons_b = itertools.tee(polygons)
# the lengths will be the length of each polygon
polygon_lengths = itertools.imap(len, polygons_a)
# for the actual points, we'll flatten the polygons into points
points = itertools.chain.from_iterable(polygons_b)
# then we'll flatten the points into values
values = itertools.chain.from_iterable(points)

# package all of that into a numpy array
all_points = numpy.fromiter(values, float)
# reshape the numpy array so we have two values for each coordinate
all_points = all_points.reshape(all_points.size // 2, 2)

# produce an iterator of lengths, but put a zero in front
polygon_positions = itertools.chain([0], polygon_lengths)
# produce another numpy array from this
# however, we take the cumulative sum
# so that each index will be the starting index of a polygon
polygon_positions = numpy.cumsum( numpy.fromiter(polygon_positions, int) )

# now for the transformation
# multiply the first coordinate of every point by *.5
all_points[:,0] *= .5

# now to get it out

# polygon_positions is all of the starting positions
# polygon_postions[1:] is the same, but shifted on forward,
# thus it gives us the end of each slice
# slice makes these all slice objects
slices = itertools.starmap(slice, itertools.izip(polygon_positions, polygon_positions[1:]))
# polygons produces an iterator which uses the slices to fetch
# each polygon
polygons = itertools.imap(all_points.__getitem__, slices)

# just iterate over the polygon normally
# each one will be a slice of the numpy array
for polygon in polygons:
    draw_polygon(polygon)

您可能会发现最好一次处理一个多边形。将每个多边形转换为一个 numpy 数组并对其进行矢量运算。这样做你可能会获得显着的速度优势。将所有数据放入 numpy 可能有点困难。

由于您的数据形状奇特，这比大多数 numpy 的东西更难。 Numpy 几乎假设了一个形状统一的数据世界。

Answer 2

使用 numpy 数组的要点是尽可能避免 for 循环。自己编写 for 循环会导致代码变慢，但使用 numpy 数组，您可以使用预定义的矢量化函数，这些函数更快（也更容易！）。

因此，对于将列表转换为数组，您可以使用：

point_buffer = np.array(point_list)

如果列表中包含(lat, lon) 之类的元素，则会将其转换为具有两列的数组。

使用该 numpy 数组，您可以轻松地一次操作所有元素。例如，要在您的问题中将每个坐标对的第一个元素乘以 0.5，您可以简单地做（假设第一个元素例如在第一列中）：

point_buffer[:,0] * 0.5

Answer 3

这样会更快：

numpy.array(point_buffer, dtype=numpy.float32)

修改数组，而不是列表。如果可能的话，显然最好避免首先创建列表。

编辑 1：分析

这里有一些测试代码，演示了 numpy 将列表转换为数组的效率（很好）。而且我的列表到缓冲区的想法只能与 numpy 相比，而不是更好。

import timeit

setup = '''
import numpy
import itertools
import struct
big_list = numpy.random.random((10000,2)).tolist()'''

old_way = '''
a = numpy.empty(( len(big_list), 2), numpy.float32)
for i,e in enumerate(big_list):
    a[i] = e
'''

normal_way = '''
a = numpy.array(big_list, dtype=numpy.float32)
'''

iter_way = '''
chain = itertools.chain.from_iterable(big_list)
a = numpy.fromiter(chain, dtype=numpy.float32)
'''

my_way = '''
chain = itertools.chain.from_iterable(big_list)
buffer = struct.pack('f'*len(big_list)*2,*chain)
a = numpy.frombuffer(buffer, numpy.float32)
'''

for way in [old_way, normal_way, iter_way, my_way]:
    print timeit.Timer(way, setup).timeit(1)

结果：

0.22445492374
0.00450378469941
0.00523579114088
0.00451488946237

编辑 2：关于数据的分层性质

如果我了解数据始终是列表列表（对象 - 多边形 - 坐标），那么这就是我要采取的方法：将数据减少到创建方形数组的最低维度（2D in这种情况下）并使用单独的数组跟踪更高级别分支的索引。这本质上是 Winston 使用 itertools 链对象的 numpy.fromiter 的想法的实现。唯一增加的想法是分支索引。

import numpy, itertools

# heirarchical list of lists of coord pairs
polys = [numpy.random.random((n,2)).tolist() for n in [5,7,12,6]]

# get the indices of the polygons:
lengs = numpy.array([0]+[len(l) for l in polys])
p_idxs = numpy.add.accumulate(lengs)

# convert the flattend list to an array:
chain = itertools.chain.from_iterable
a = numpy.fromiter(chain(chain(polys)), dtype=numpy.float32).reshape(lengs.sum(), 2)

# transform the coords
a *= .5

# get a transformed polygon (using the indices)
def get_poly(n):
    i0 = p_idxs[n]
    i1 = p_idxs[n+1]
    return a[i0:i1]

print 'poly2', get_poly(2)
print 'poly0', get_poly(0)