title: Indexing a large 3D HDF5 dataset for subsetting based on 2D condition
body: <p>I have a large 3D HDF5 dataset that represents location (X,Y) and time for a certain variable. Next, I have a 2D numpy array containing a classification for the same (X,Y) location. What I would like to achieve is that I can extract all time series from the 3D HDF5 dataset that fall under a certain class in the 2D array.</p>
  
  <p>Here's my example:</p>
  
  <pre><code>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
  &lt;HDF5 dataset "NDVI": shape (1319, 2063, 53), type "&lt;f4"&gt;
  
  # Let's make a random 1319 x 2063 classification containing class numbers 0-4
  classification = np.random.randint(5, size=(1319, 2063))
  </code></pre>
  
  <p>Now we have our 3D HDF5 dataset and a 2D classification. Let's look for pixels that fall under class number '3'</p>
  
  <pre><code># Look for the X,Y locations that have class number '3'
  idx = np.where(classification == 3)
  </code></pre>
  
  <p>This returns me a tuple of size 2 that contains the X,Y pairs that match the condition, and in my random example the amount of pairs is 544433. How should I use now this <code>idx</code> variable to create a 2D array of size (544433,53) that contains the 544433 time series of the pixels that have classification class number '3'?</p>
  
  <p>I did some testing with fancy indexing and pure 3D numpy arrays and this example would work just fine:</p>
  
  <pre><code>subset = 3D_numpy_array[idx[0],idx[1],:]
  </code></pre>
  
  <p>However, the HDF5 dataset is too large to convert to a numpy array; when I'm trying to use the same indexing method directly on the HDF5 dataset:</p>
  
  <pre><code># Try to use fancy indexing directly on HDF5 dataset
  NDVI_subset = np.array(NDVI_data[idx[0],idx[1],:])
  </code></pre>
  
  <p>It throws me an error:</p>
  
  <pre><code>Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
  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
  </code></pre>
  
  <p>Another thing I tried is to <code>np.repeat</code> the classification array in the 3rd dimension to create a 3D array that matches the shape of the HDF5 dataset. The <code>idx</code> variable than gets a tuple of size 3:</p>
  
  <pre><code>classification_3D = np.repeat(np.reshape(classification,(1319,2063,1)),53,axis=2)
  idx = np.where(classification == 3)
  </code></pre>
  
  <p>But the following statement than throws the exact same error:</p>
  
  <pre><code>NDVI_subset = np.array(NDVI_data[idx])
  </code></pre>
  
  <p>Is this because the HDF5 dataset works differently compared to a pure numpy array? The documentation does say "Selection coordinates must be given in increasing order"</p>
  
  <p>Does anyone in that case have a suggestion on how I could get this to work without having to read in the full HDF5 dataset into memory (which does not work)?
  Thank you very much!</p>