title: something like partial/selective np.ndenumerate?
body: <p>I have read this <a href="https://stackoverflow.com/questions/9570050/numpy-ndarray-enumeration-over-a-proper-subset-of-the-dimensions">post</a> and suspect I need to understand indexing better. I would like to do something like <code>np.ndenumerate</code>, but (in this particular case) only over the first three dimensions, returning a series of vectors:</p>
  
  <pre><code>x = np.linspace(0., 1., 4)
  y = np.linspace(0., 1., 5)
  z = np.linspace(0., 1., 2)
  X, Y, Z = np.meshgrid(x, y, z)   # define a grid
  F = np.zeros(X.shape + (3,))     # a 3D vector field 
  F = np.random.rand(5*4*2*3).reshape(5,4,2,3)   # added this later just to give non-zero for testing
  
  thetas = np.linspace(0, 2.*np.pi, 21)[:-1]    # drop the duplicate
  v = np.array([(np.cos(theta), np.sin(theta), 0.0) for theta in thetas])
  
  for tup, vec in magical_enumerate(F, axis=(0,1,2)):  # it's magic! (need correct expression)
      F(tup) = np.cross(v, vec).sum(axis=0)     # sum over 20 vectors in v
  </code></pre>
  
  <p>Is there a way to do this without lots of loops or list interpretations? The grids will be large, so numpythony and speed are appreciated. What about non-sequential dimensions (e.g. axis=(0,2))?  Thanks! </p>