为基于 2D 条件的子集索引大型 3D HDF5 数据集

Posted 2023-03-11

【中文标题】为基于 2D 条件的子集索引大型 3D HDF5 数据集

【英文标题】：Indexing a large 3D HDF5 dataset for subsetting based on 2D condition

【发布时间】：2016-12-10 06:06:26

【问题描述】：

我有一个大型 3D HDF5 数据集，表示某个变量的位置 (X,Y) 和时间。接下来，我有一个 2D numpy 数组，其中包含相同 (X,Y) 位置的分类。我想要实现的是，我可以从 3D HDF5 数据集中提取属于 2D 数组中某个类的所有时间序列。

这是我的例子：

import numpy as np
import h5py

# Open the HDF5 dataset
NDVI_file = 'NDVI_values.hdf5'
f_NDVI = h5py.File(NDVI_file,'r')
NDVI_data = f_NDVI["NDVI"]

# See what's in the dataset
NDVI_data
<HDF5 dataset "NDVI": shape (1319, 2063, 53), type "<f4">

# Let's make a random 1319 x 2063 classification containing class numbers 0-4
classification = np.random.randint(5, size=(1319, 2063))

现在我们有了 3D HDF5 数据集和 2D 分类。让我们寻找属于“3”类的像素

# Look for the X,Y locations that have class number '3'
idx = np.where(classification == 3)

这将返回一个大小为 2 的元组，其中包含与条件匹配的 X、Y 对，在我的随机示例中，对的数量是 544433。我现在应该如何使用这个 idx 变量来创建二维数组大小 (544433,53) 包含分类类别号为“3”的像素的 544433 个时间序列？

我用花哨的索引和纯 3D numpy 数组做了一些测试，这个例子可以正常工作：

subset = 3D_numpy_array[idx[0],idx[1],:]

但是，HDF5 数据集太大，无法转换为 numpy 数组；当我尝试直接在 HDF5 数据集上使用相同的索引方法时：

# Try to use fancy indexing directly on HDF5 dataset
NDVI_subset = np.array(NDVI_data[idx[0],idx[1],:])

它给我一个错误：

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "h5py\_objects.pyx", line 54, in h5py._objects.with_phil.wrapper     (C:\aroot\work\h5py\_objects.c:2584)
File "h5py\_objects.pyx", line 55, in h5py._objects.with_phil.wrapper (C:\aroot\work\h5py\_objects.c:2543)
File "C:\Users\vtrichtk\AppData\Local\Continuum\Anaconda2\lib\site-packages\h5py\_hl\dataset.py", line 431, in __getitem__
selection = sel.select(self.shape, args, dsid=self.id)
File "C:\Users\vtrichtk\AppData\Local\Continuum\Anaconda2\lib\site-packages\h5py\_hl\selections.py", line 95, in select
sel[args]
File "C:\Users\vtrichtk\AppData\Local\Continuum\Anaconda2\lib\site-packages\h5py\_hl\selections.py", line 429, in __getitem__
raise TypeError("Indexing elements must be in increasing order")
TypeError: Indexing elements must be in increasing order

我尝试的另一件事是np.repeat 3 维中的分类数组，以创建与 HDF5 数据集形状匹配的 3D 数组。 idx 变量得到一个大小为 3 的元组：

classification_3D = np.repeat(np.reshape(classification,(1319,2063,1)),53,axis=2)
idx = np.where(classification == 3)

但以下语句会引发完全相同的错误：

NDVI_subset = np.array(NDVI_data[idx])

这是因为 HDF5 数据集与纯 numpy 数组的工作方式不同吗？文档确实说“选择坐标必须按递增顺序给出”

在这种情况下，是否有人建议我如何使其工作而无需将完整的 HDF5 数据集读入内存（这不起作用）？ 非常感谢！

【问题讨论】：

h5py doc 对高级索引或花哨索引有何评论？我会研究它，然后设置一个小得多的测试用例，我可以在移动两个 3d 之前在 2d 数组上测试这种索引。我可以在哪里打印所有值。当然，H5 在索引方面可能会受到更多限制。

docs.h5py.org/en/latest/high/dataset.html#fancy-indexing

【参考方案1】：

h5py 中的高级/精美索引不如np.ndarray 那样通用。

设置一个小测试用例：

import h5py
f=h5py.File('test.h5','w')
dset=f.create_dataset('data',(5,3,2),dtype='i')
dset[...]=np.arange(5*3*2).reshape(5,3,2)
x=np.arange(5*3*2).reshape(5,3,2)

ind=np.where(x%2)

我可以选择所有奇数：

In [202]: ind
Out[202]: 
(array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4], dtype=int32),
 array([0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2], dtype=int32),
 array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int32))

In [203]: x[ind]
Out[203]: array([ 1,  3,  5,  7,  9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29])
In [204]: dset[ind]
...
TypeError: Indexing elements must be in increasing order

我可以使用如下列表对单个维度进行索引：dset[[1,2,3],...]，但重复索引值或更改顺序会产生错误，dset[[1,1,2,2],...] 或 dset[[2,1,0],...]。 dset[:,[0,1],:] 没问题。

几个切片都可以，dset[0:3,1:3,:]，或者切片和列表，dset[0:3,[1,2],:]。

但是 2 个列表 dset[[0,1,2],[1,2],:] 产生了一个

TypeError: Only one indexing vector or array is currently allowed for advanced selection

所以np.where 的索引元组在几个方面是错误的。

我不知道这有多少是h5存储的约束，有多少只是h5py模块中的不完整开发。也许两者兼而有之。

因此，您需要从文件中加载更简单的块，并对生成的 numpy 数组执行更高级的索引。

在我的odd values 案例中，我只需要这样做：

In [225]: dset[:,:,1]
Out[225]: 
array([[ 1,  3,  5],
       [ 7,  9, 11],
       [13, 15, 17],
       [19, 21, 23],
       [25, 27, 29]])

【讨论】：

感谢分享这个例子，它很好地阐明了 h5py 索引的局限性。我看到了您针对odd values 案例的解决方案，但我担心在我的案例中没有这样简单的解决方案。我必须分别为每个时间步切出 2D 数组，然后对 Numpy 数组进行花哨的索引，但这似乎不是一个快速而优雅的解决方案。