title: numpy 3d to 2d transformation based on 2d mask array
body: <p>If I have an ndarray like this:</p>
  
  <pre><code>&gt;&gt;&gt; a = np.arange(27).reshape(3,3,3)
  &gt;&gt;&gt; a
  array([[[ 0,  1,  2],
          [ 3,  4,  5],
          [ 6,  7,  8]],
  
         [[ 9, 10, 11],
          [12, 13, 14],
          [15, 16, 17]],
  
         [[18, 19, 20],
          [21, 22, 23],
          [24, 25, 26]]])
  </code></pre>
  
  <p>I know I can get the maximum along a certain axis using <code>np.max(axis=...)</code>:</p>
  
  <pre><code>&gt;&gt;&gt; a.max(axis=2)
  array([[ 2,  5,  8],
         [11, 14, 17],
         [20, 23, 26]])
  </code></pre>
  
  <p>Alternatively, I could get the <strong>indices</strong> along that axis which correspond to the maximum values from:</p>
  
  <pre><code>&gt;&gt;&gt; indices = a.argmax(axis=2)
  &gt;&gt;&gt; indices
  array([[2, 2, 2],
         [2, 2, 2],
         [2, 2, 2]])
  </code></pre>
  
  <p>My question -- Given the array <code>indices</code> and the array <code>a</code>, is there an elegant way to reproduce the array the array returned by <code>a.max(axis=2)</code>?</p>
  
  <p>This would probably work:</p>
  
  <pre><code>import itertools as it
  import numpy as np
  def apply_mask(field,indices):
      data = np.empty(indices.shape)
  
      #It seems highly likely that there is a more numpy-approved way to do this.
      idx = [range(i) for i in indices.shape]
      for idx_tup,zidx in zip(it.product(*idx),indices.flat):
          data[idx_tup] = field[idx_tup+(zidx,)]
      return data
  </code></pre>
  
  <p>But, it seems pretty hacky/inefficient.  It also doesn't allow for me to use this with any axis other than the "last" axis.  Is there a numpy function (or some use of magical numpy indexing) to make this work?  The naive <code>a[:,:,a.argmax(axis=2)]</code> doesn't work.</p>
  
  <p><strong>UPDATE</strong>:</p>
  
  <p>It seems the following also works (and is a little nicer):</p>
  
  <pre><code>import numpy as np
  def apply_mask(field,indices):
      data = np.empty(indices.shape)
  
      for idx_tup,zidx in np.ndenumerate(indices):
          data[idx_tup] = field[idx_tup+(zidx,)]
  
      return data
  </code></pre>
  
  <hr>
  
  <p>I would like to do this because I would like to extract the indices based on the data in 1 array (typically using <code>argmax(axis=...)</code>) and use those indices to pull data out of a bunch of other (equivalently shaped) arrays.  I'm open to alternative ways to accomplish this (e.g. using boolean masked arrays).  However, I like the "safety" that I get using these "index" arrays.  With this I am guaranteed to have the right number of elements to create a new array which looks like a 2d "slice" through the 3d field.</p>