python - 有没有比 np.where 更好的方法使用条件语句来索引大型数组？

Question

抱歉文字墙。我试图缩短它，但我认为一切都会对愿意阅读全部内容的人有所帮助。

我有 xyz 点云，我正在尝试将它们网格化为 DEM(数字高程模型，如果您不熟悉 - 只是高程的二维数组)。我的 DEM 需要具有比点云低得多的分辨率(相反，它们不需要几乎一样高的分辨率)，因此每个 DEM 单元中的点云中有大约 10 个点。我的代码迭代最终 DEM 的行和列，并根据提供的 XYZ 坐标计算每个像元的高程值。为此，我指定网格大小和间距，然后代码计算每个网格单元的 x 和 y 最大值和最小值。然后，它找到 x 和 y 位于该像元的最大值和最小值内的所有 z 值，拒绝异常值并取剩余值的平均值来确定该像元的最终值。 xyz 值存储在如下所示的数组中:

xyz = np.array([[x1, y1, z1],
                [x2, y2, z2],
                [x3, y3, z3]])       #     with tens of thousands of xyz cominations

我当前的方法涉及创建 X 和 Y 中的像元边界值列表，然后将边界内具有 x 和 y 值的所有 z 值添加到高程列表中:

dx = 0.5           # cell size in meters

xmin = np.min(xyz[:,0])
xmax = np.max(xyz[:,0])
ymin = np.min(xyz[:,1])
ymax = np.max(xyz[:,1])
xoffset = ((xmax-xmin) % dx)/2.0
yoffset = ((ymax-ymin) % dx)/2.0


xlims = np.arange(xmin+xoffset, xmax, dx)   # list of grid cell limits in x
ylims = np.arange(ymin+yoffset, ymax, dx)   # list of grid cell limits in y

DEM = np.empty((len(ylims)-1, len(xlims)-1), 'float')  # declares output array

for i in range(DEM.shape[0]):           # iterate over rows of final DEM
    for j in range(DEM.shape[1]):       # iterate over columns of final DEM
        bottom = ylims[i]               
        top = ylims[i+1]                
        left = xlims[j]                 # these rows just pick minimum and
        right = xlims[j+1]              # maximum of cell [i,j]

        elevations = xyz[np.where(((xyz[:,0] > left) &
                            (xyz[:,0] < right)) & 
                            ((xyz[:,1] > bottom) & 
                            (xyz[:,1] < top)))[0]][:,2]
        elevations = reject_outliers(elevations)

        if len(elevations) == 0:
            elevation = np.nan
        else:
            elevation = np.mean(elevations)

        DEM[i,j] = elevation

这可行，但我必须做数百个 DEM，每个 DEM 都有数十万个点，所以如果我这样做，我需要等待一周的时间来等待我的计算机通过此操作。对我来说它也显得很笨重。有没有办法简化这个过程？

Answer 1

一种可能性是找到这些点适合一次的段，然后循环遍历这些组，这样您就不必不断地屏蔽相同的元素，即使它们只属于一个段。

实现此目的的一种方法是使用 numpy 的内置 searchsorted:

xmin = np.min(xyz[:,0])
xmax = np.max(xyz[:,0])
ymin = np.min(xyz[:,1])
ymax = np.max(xyz[:,1])
xoffset = ((xmax-xmin) % dx)/2.0
yoffset = ((ymax-ymin) % dx)/2.0

xlims = np.arange(xmin+xoffset, xmax, dx)   # list of grid cell limits in x
ylims = np.arange(ymin+yoffset, ymax, dx)   # list of grid cell limits in y

DEM = np.empty((len(ylims) - 1, len(xlims) - 1), dtype=float)  # declares output array

# Find the bins that each point fit into
x_bins = np.searchsorted(xlims, xyz[:, 0]) - 1
y_bins = np.searchsorted(ylims, xyz[:, 1]) - 1

for i in range(DEM.shape[0]):           # iterate over rows of final DEM
    y_mask = y_bins == i
    for j in range(DEM.shape[1]):
        elevations = xyz[y_mask & (x_bins == j), 2]
        elevations = reject_outliers(elevations)
        if len(elevations) == 0:
            elevations = np.nan
        else:
            elevations = np.mean(elevations)

        DEM[i, j] = elevations

当我使用 timeit 和 xyz = randn(100000, 3) 和 dx = 0.1 分析已经列出的替代方案时(在将reject_outliers定义为 lambda x: x 之后)，我得到了以下时间:

1. (您的)where 方法:10 个循环，3 个循环中最好的:每个循环 6.69 秒
2. (dermon 的)逐步方法:10 个循环，3 次最佳:每个循环 16.7 秒
3. 搜索排序方法:10 个循环，3 个循环中最好的:每个循环 2.11 秒

但是，如果您愿意使用 Pandas，您可以修改代码，以便将双 for 循环替换为 Pandas 出色的 groupby 功能:

elevation_df = pd.DataFrame({'x_bins': x_bins, 'y_bins': y_bins, 'z': xyz[:, 2]})
for x_y_bins, data in elevation_df.groupby(['x_bins', 'y_bins']):
    elevations = reject_outliers(data['z'])
    elevations = data['z']
    if len(elevations) == 0:
        elevations = np.nan
    else:
        elevations = np.mean(elevations)
    if 0 <= x_y_bins[1] < DEM.shape[0] and 0 <= x_y_bins[0] < DEM.shape[1]: 
        DEM[x_y_bins[1], x_y_bins[0]] = elevations

这几乎将时间缩短了 2 倍(10 个循环，最好是 3 个:每个循环 1.11 秒)。

我还应该注意到，由于您的 np.arange 命令，您似乎已经排除了范围中的某些点。在上面，我假设您打算排除这些点，但如果您想包含您可以使用的所有数据:

xlims = np.arange(xmin-xoffset, xmax+dx, dx)
ylims = np.arange(ymin-yoffset, ymax+dx, dx)

如果您使用这些范围，您可以将我之前的 for 循环修改为:

for i in np.unique(y_bins):
    y_mask = y_bins == i
    for j in np.unique(x_bins[y_mask]):

这将我之前的搜索排序示例的 timeit 结果降低到 10 个循环，最好是 3 次:每个循环 1.56 秒，这至少更接近 Pandas groupby。