基于2D条件对大型3D HDF5数据集进行子集化索引

Question

我有一个大的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变量创建一个包含分类类号为“3”的像素的544433时间序列的二维大小数组（544433,53）？

我做了一些测试用花哨的索引和纯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在第三分类数组用于创建与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数据集到内存中（这是行不通的）？ 非常感谢！

Answer 1

h5py中的高级/花式索引几乎没有np.ndarray一般。

设置一个小的测试情况：

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],...]，但重复指数

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) 

我可以选择所有奇数值值或更改顺序会产生错误，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数组执行fancier索引。

在我odd values情况下，我只需要做：

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