With readability, np.choose definitely looks great.

If performance is of essence, you can calculate the linear indices and then use np.take or use a flattened version with .ravel() and extract those specific elements from val_arr. The implementation would look something like this -

def linidx_take(val_arr,z_indices):

    # Get number of columns and rows in values array
     _,nC,nR = val_arr.shape

     # Get linear indices and thus extract elements with np.take
    idx = nC*nR*z_indices + nR*np.arange(nR)[:,None] + np.arange(nC)
    return np.take(val_arr,idx) # Or val_arr.ravel()[idx]

Runtime tests and verify results -

Ogrid based solution from here is made into a generic version for these tests, like so :

In [182]: def ogrid_based(val_arr,z_indices):
     ...:   v_shp = val_arr.shape
     ...:   y,x = np.ogrid[0:v_shp[1], 0:v_shp[2]]
     ...:   return val_arr[z_indices, y, x]
     ...: 

Case #1: Smaller datasize

In [183]: val_arr = np.random.rand(30,30,30)
     ...: z_indices = np.random.randint(0,30,(30,30))
     ...: 

In [184]: np.allclose(z_indices.choose(val_arr),ogrid_based(val_arr,z_indices))
Out[184]: True

In [185]: np.allclose(z_indices.choose(val_arr),linidx_take(val_arr,z_indices))
Out[185]: True

In [187]: %timeit z_indices.choose(val_arr)
1000 loops, best of 3: 230 µs per loop

In [188]: %timeit ogrid_based(val_arr,z_indices)
10000 loops, best of 3: 54.1 µs per loop

In [189]: %timeit linidx_take(val_arr,z_indices)
10000 loops, best of 3: 30.3 µs per loop

Case #2: Bigger datasize

In [191]: val_arr = np.random.rand(300,300,300)
     ...: z_indices = np.random.randint(0,300,(300,300))
     ...: 

In [192]: z_indices.choose(val_arr) # Seems like there is some limitation here with bigger arrays.
Traceback (most recent call last):

  File "<ipython-input-192-10c3bb600361>", line 1, in <module>
    z_indices.choose(val_arr)

ValueError: Need between 2 and (32) array objects (inclusive).


In [194]: np.allclose(linidx_take(val_arr,z_indices),ogrid_based(val_arr,z_indices))
Out[194]: True

In [195]: %timeit ogrid_based(val_arr,z_indices)
100 loops, best of 3: 3.67 ms per loop

In [196]: %timeit linidx_take(val_arr,z_indices)
100 loops, best of 3: 2.04 ms per loop

If I understand correctly, you are looking for advances indexing of first dimension of B. You can use np.indices to create the indices required for the other two dimensions of B and use advanced indexing:

idx = np.indices(A.shape)
C = B[A,idx[0],idx[1]]

Example:

B = np.random.rand(10,20,30)
A = np.array([[   1,    0,    0,    1,    0],
       [   0,    2,    0,    2,    4],
       [   1,    7,    0,    5,    5],
       [   2,    1,    0,    1,    7],
       [   0,    7,    0,    1,    4],
       [   0,    0,    3,    6,    0]])

print(C[4,1]==B[7,4,1])
#True

The question is somewhat ambiguous, but if the intent is to increment some unknown array B at indices (0,0,0), (1,0,1), (3,1,0), and (3,1,1), then the following should be fine:

B[(A.ravel(), ) + np.unravel_index(np.arange(np.prod(A.shape)), A.shape)] += increment

For example:

A = np.array([[0,1],[3,3]])
B = np.zeros((4,2,2), dtype=int)
increment = 1

B[(A.ravel(), ) + np.unravel_index(np.arange(np.prod(A.shape)), A.shape)] += increment

>>> B
array([[[1, 0],
        [0, 0]],

       [[0, 1],
        [0, 0]],

       [[0, 0],
        [0, 0]],

       [[0, 0],
        [1, 1]]])

Another way of doing the same thing is:

w, h = A.shape
indices = (A.ravel(),) + tuple(np.mgrid[:w, :h].reshape(2, -1))

# then
B[indices] += increment

In [107]: arr_3d = np.arange(2*3*4).reshape(2,3,4)                                                           
In [108]: arr_2d = np.array(([3,2,0], [2,3,2]))                                                              
In [109]: arr_2d.shape                                                                                       
Out[109]: (2, 3)
In [110]: arr_3d[[[0],[1]],[0,1,2],arr_2d]                                                                   
Out[110]: 
array([[ 3,  6,  8],
       [14, 19, 22]])

[[0],[1]],[0,1,2] broadcast with each other to index a (2,3) block, the same size as `arr_2d.

ix_ can be used to construct those 2 indices:

In [114]: I,J = np.ix_(range(2), range(3))                                                                   
In [115]: I,J                                                                                                
Out[115]: 
(array([[0],
        [1]]), array([[0, 1, 2]]))
In [116]: arr_3d[I, J, arr_2d]                                                                               
Out[116]: 
array([[ 3,  6,  8],
       [14, 19, 22]])

You can use take array method:

import numpy

a = numpy.array([[5,6,7,8,9],[10,11,12,14,15]])

idx = numpy.random.randint(5, size=(3, a.shape[0], a.shape[1]))

print a.take(idx)

If your matrices get much bigger than 3x3, to the point that memory involved in np.ogrid is an issue, and if your indexes remain binary, you could also do:

np.where(a, b[1], b[0])

But other than that corner case (or if you like code golfing one-liners) the other answer is probably better.