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数组并提取数据呢？

我通过一个C库从一个形状文件中读取数据，该库将数据以常规python列表的形式提供给我。我明白如果这个库直接给我numpy数组中的坐标，就不需要“填充”numpy数组了。但不幸的是，我的数据起点是一个普通的python列表。更重要的是，我想了解如何在python中快速填充numpy数组。

说明

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

在我的实际应用中，我有一个坐标点的形状文件，并且我有一个API可以获取给定对象的点。大约有20万个对象。所以我反复调用一个函数GetShapeCoords( i )来获取对象i的坐标。这个函数返回一个列表的列表，每个子列表是经纬度对的列表，之所以是列表的列表是因为有些对象是多部分的（即多边形）。然后，在我原来的代码中，当我读取每个对象的点时，我通过调用一个常规的python函数对每个点进行变换，然后使用PIL绘制变换后的点。整个过程大约花了20秒来绘制所有20万个多边形。虽然不算太糟，但还有很大的改进空间。我注意到至少有一半的20秒花在了变换逻辑上，所以我想用numpy来处理。我的原始实现是逐个读取对象，并将所有子列表中的点追加到一个大的numpy数组中，然后我可以在numpy中进行数学运算。

所以，我现在明白，简单地将整个python列表传递给numpy是设置大数组的正确方法。但在我的情况下，我一次只读取一个对象。所以我可以做的一件事是将点保存在一个大的python列表的列表中。当我以这种方式编译了一些大量对象的点（比如说，10000个对象）后，我可以简单地将这个庞大的列表分配给numpy。

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

(a) 真的可以将那个大而不规则的列表的列表的列表传递给numpy，并且快速处理吗？

(b) 然后我想对那个庞大树的叶子中的所有点进行变换。有什么表达式可以让numpy“进入每个子列表，然后进入每个子子列表，然后对那些子子列表中的每个坐标对进行操作，比如将第一个（经度坐标）乘以0.5”？我可以这样做吗？

(c) 最后，我需要将那些变换后的坐标提取出来以便绘制。

Winston下面的回答似乎给出了我如何使用itertools来做到这一点的一些提示。我想做的基本上和Winston做的差不多，都是将列表扁平化。但我不能简单地将其扁平化。当我去绘制数据时，我需要知道一个多边形何时结束，另一个何时开始。所以，我想如果能有一种快速标记每个多边形结束的方法（例如，用一个特殊的坐标对如(-1000, -1000)），那我就可以像Winston的回答中那样使用itertools进行扁平化，然后在numpy中进行变换。然后我需要使用PIL从点到点进行绘制，在这里我想我需要将修改后的numpy数组重新分配回python列表，然后在这个列表中使用常规的python循环进行绘制。这似乎是我最好的选择，除非我写一个C模块来一步到位处理所有读取和绘制的工作。

Answer 1

使用numpy数组的主要目的是尽量避免使用for循环。自己写for循环会导致代码运行得很慢，而使用numpy数组可以利用一些预定义的向量化函数，这样不仅速度快，而且更简单！

所以，如果你想把一个列表转换成数组，可以使用：

point_buffer = np.array(point_list)

如果这个列表里的元素是像(lat, lon)这样的坐标对，那么转换后就会变成一个有两列的数组。

有了这个numpy数组，你可以很方便地一次性处理所有元素。比如，如果你想把每对坐标的第一个元素乘以0.5，就可以简单地这样做（假设第一个元素在第一列）：

point_buffer[:,0] * 0.5

Answer 2

这样会更快：

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：关于数据的层次结构

如果我理解得没错，数据总是一个列表的列表的列表（对象 - 多边形 - 坐标），那么我会采取这样的方式：将数据简化到最低维度，形成一个正方形数组（在这个例子中是二维的），并用一个单独的数组来跟踪更高层级的索引。这基本上是实现了Winston的想法，使用numpy.fromiter和itertools链对象。唯一增加的想法是分支索引。

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)

Answer 3

你提到你的数据是“坐标的列表的列表的列表”。从这点来看，我猜你的数据提取可能是这样的：

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基本上是针对形状统一的数据设计的。