Falko's suggestion to use contourf works with a bit of finagling. It's a bit limited since at least my version of contourf has a few bugs where it sometimes renders one of the planes in front of other planes it should be behind, but for now only plotting either the three front or three back sides of the cube will do:

import numpy as np
import math
import matplotlib.pyplot as plot
import mpl_toolkits.mplot3d.axes3d as axes3d

def cube_marginals(cube, normalize=False):
    c_fcn = np.mean if normalize else np.sum
    xy = c_fcn(cube, axis=0)
    xz = c_fcn(cube, axis=1)
    yz = c_fcn(cube, axis=2)
    return(xy,xz,yz)

def plotcube(cube,x=None,y=None,z=None,normalize=False,plot_front=False):
    """Use contourf to plot cube marginals"""
    (Z,Y,X) = cube.shape
    (xy,xz,yz) = cube_marginals(cube,normalize=normalize)
    if x == None: x = np.arange(X)
    if y == None: y = np.arange(Y)
    if z == None: z = np.arange(Z)

    fig = plot.figure()
    ax = fig.gca(projection='3d')

    # draw edge marginal surfaces
    offsets = (Z-1,0,X-1) if plot_front else (0, Y-1, 0)
    cset = ax.contourf(x[None,:].repeat(Y,axis=0), y[:,None].repeat(X,axis=1), xy, zdir='z', offset=offsets[0], cmap=plot.cm.coolwarm, alpha=0.75)
    cset = ax.contourf(x[None,:].repeat(Z,axis=0), xz, z[:,None].repeat(X,axis=1), zdir='y', offset=offsets[1], cmap=plot.cm.coolwarm, alpha=0.75)
    cset = ax.contourf(yz, y[None,:].repeat(Z,axis=0), z[:,None].repeat(Y,axis=1), zdir='x', offset=offsets[2], cmap=plot.cm.coolwarm, alpha=0.75)

    # draw wire cube to aid visualization
    ax.plot([0,X-1,X-1,0,0],[0,0,Y-1,Y-1,0],[0,0,0,0,0],'k-')
    ax.plot([0,X-1,X-1,0,0],[0,0,Y-1,Y-1,0],[Z-1,Z-1,Z-1,Z-1,Z-1],'k-')
    ax.plot([0,0],[0,0],[0,Z-1],'k-')
    ax.plot([X-1,X-1],[0,0],[0,Z-1],'k-')
    ax.plot([X-1,X-1],[Y-1,Y-1],[0,Z-1],'k-')
    ax.plot([0,0],[Y-1,Y-1],[0,Z-1],'k-')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    plot.show()

plot_front=True plot_front=False Other data (not shown)

Nested matrix produced this:

\documentclass{article}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{xcolor}
\begin{document}
\begin{align*}
\colorlet{tmp}{.}\color{red}\left[\color{tmp}
\begin{matrix}
\left[
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}\right]\,& \color{red}\raisebox{-1.5em}{\scalebox{2}{,}}\\
\\
%
\left[
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}
\begin{bmatrix}
1&0&0\\
0&1&0\\
0&0&1
\end{bmatrix}\right]&
\end{matrix}
\vphantom{\begin{array}{c}~\\[8.5em]\end{array}}
\color{red}\right]
\end{align*}
\end{document}

The anatomy of an ndarray in NumPy looks like this red cube below: (source: Physics Dept, Cornell Uni)

Once you leave the 2D space and enter 3D or higher dimensional spaces, the concept of rows and columns doesn't make much sense anymore. But still you can intuitively understand 3D arrays. For instance, considering your example:

In [41]: b
Out[41]: 
array([[[ 1,  2,  3],
        [ 4,  5,  6]],

       [[ 7,  8,  9],
        [10, 11, 12]]])

In [42]: b.shape
Out[42]: (2, 2, 3)

Here the shape of b is (2, 2, 3). You can think about it like, we've two (2x3) matrices stacked to form a 3D array. To access the first matrix you index into the array b like b[0] and to access the second matrix, you index into the array b like b[1].

# gives you the 2D array (i.e. matrix) at position `0`
In [43]: b[0]
Out[43]: 
array([[1, 2, 3],
       [4, 5, 6]])


# gives you the 2D array (i.e. matrix) at position 1
In [44]: b[1]
Out[44]: 
array([[ 7,  8,  9],
       [10, 11, 12]])

However, if you enter 4D space or higher, it will be very hard to make any sense out of the arrays itself since we humans have hard time visualizing 4D and more dimensions. So, one would rather just consider the ndarray.shape attribute and work with it.

More information about how we build higher dimensional arrays using (nested) lists:

For 1D arrays, the array constructor needs a sequence (tuple, list, etc) but conventionally list is used.

In [51]: oneD = np.array([1, 2, 3,])    
In [52]: oneD.shape
Out[52]: (3,)

For 2D arrays, it's list of lists but can also be tuple of lists or tuple of tuples etc:

In [53]: twoD = np.array([[1, 2, 3], [4, 5, 6]])
In [54]: twoD.shape
Out[54]: (2, 3)

For 3D arrays, it's list of lists of lists:

In [55]: threeD = np.array([[[1, 2, 3], [2, 3, 4]], [[5, 6, 7], [6, 7, 8]]])

In [56]: threeD.shape
Out[56]: (2, 2, 3)

P.S. Internally, the ndarray is stored in a memory block as shown in the below picture. (source: Enthought)

For better performance, avoid calling ax.scatter multiple times, if possible. Instead, pack all the x,y,z coordinates and colors into 1D arrays (or lists), then call ax.scatter once:

ax.scatter(x, y, z, c=volume.ravel())

The problem (in terms of both CPU time and memory) grows as size**3, where size is the side length of the cube.

Moreover, ax.scatter will try to render all size**3 points without regard to the fact that most of those points are obscured by those on the outer shell.

It would help to reduce the number of points in volume -- perhaps by summarizing or resampling/interpolating it in some way -- before rendering it.

We can also reduce the CPU and memory required from O(size**3) to O(size**2) by only plotting the outer shell:

import functools
import itertools as IT
import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def cartesian_product_broadcasted(*arrays):
    """
    http://stackoverflow.com/a/11146645/190597 (senderle)
    """
    broadcastable = np.ix_(*arrays)
    broadcasted = np.broadcast_arrays(*broadcastable)
    dtype = np.result_type(*arrays)
    rows, cols = functools.reduce(np.multiply, broadcasted[0].shape), len(broadcasted)
    out = np.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

# @profile  # used with `python -m memory_profiler script.py` to measure memory usage
def main():
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection='3d')

    size = 512
    volume = np.random.rand(size, size, size)
    x, y, z = cartesian_product_broadcasted(*[np.arange(size, dtype='int16')]*3).T
    mask = ((x == 0) | (x == size-1) 
            | (y == 0) | (y == size-1) 
            | (z == 0) | (z == size-1))
    x = x[mask]
    y = y[mask]
    z = z[mask]
    volume = volume.ravel()[mask]

    ax.scatter(x, y, z, c=volume, cmap=plt.get_cmap('Greys'))
    plt.show()

if __name__ == '__main__':
    main()

But note that even when plotting only the outer shell, to achieve a plot with size=512 we still need around 1.3 GiB of memory. Also beware that even if you have enough total memory but, due to a lack of RAM, the program uses swap space, then the overall speed of the program will slow down dramatically. If you find yourself in this situation, then the only solution is to find a smarter way to render an acceptable image using fewer points, or to buy more RAM.

This answer creates a 3D barplot as an alternative method for representing your data. This will allow you to visualize the surface with each column representing a single z-value. It is also possible to view a flattened heat map of the array by rotating the figure to display the 'bottom'.

NOTE: I used a different, larger sample set in order to create the examples. This is just to show that the method works for more complex inputs.

Example Output

This is a barplot of a larger sample set. You can clearly see the data being represented even though the 'surface' of the data is not flat.

This is a bottom-up view of the figure. It represents a heatmap of your input data.

Code

from matplotlib import cbook
from matplotlib import cm
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np

# Load and format data
dem = cbook.get_sample_data('jacksboro_fault_dem.npz', np_load=True)
z = dem['elevation']
nrows, ncols = z.shape
x = np.linspace(dem['xmin'], dem['xmax'], ncols)
y = np.linspace(dem['ymin'], dem['ymax'], nrows)
x, y = np.meshgrid(x, y)

region = np.s_[5:50, 5:50]
x, y, z = x[region].ravel(), y[region].ravel(), z[region].ravel()

# parameters for bar3d(...) func
bottom = np.full_like(z, np.min(z))
width = (np.max(x)-np.min(x))/np.sqrt(np.shape(x)[0])
depth = (np.max(y)-np.min(y))/np.sqrt(np.shape(y)[0])

# creating color_values for the figure
offset = z + np.abs(z.min())
fracs = offset.astype(float)/offset.max()
norm = colors.Normalize(fracs.min(), fracs.max())
color_values = cm.jet(norm(fracs.tolist()))

# Set up and display the plot
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.bar3d(x,y,bottom,width,depth,z, color=color_values, shade=True)

## Used to view the figure from the bottom
# ax.view_init(270, 0)

plt.show()

Does this help?

Documentation

3D Bar Charts from Matplotlib

Another example using the sample data from Matplotlib

If you have a dset like that, and you want to just get the 1 values, you could use nonzero, which "returns a tuple of arrays, one for each dimension of a, containing the indices of the non-zero elements in that dimension.".

For example, we can make a simple 3d array:

>>> import numpy
>>> numpy.random.seed(29)
>>> d = numpy.random.randint(0, 2, size=(3,3,3))
>>> d
array([[[1, 1, 0],
        [1, 0, 0],
        [0, 1, 1]],

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

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

and find where the nonzero elements are located:

>>> d.nonzero()
(array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2]), array([0, 0, 1, 2, 2, 0, 0, 1, 2, 2, 0, 0, 1, 2]), array([0, 1, 0, 1, 2, 1, 2, 0, 1, 2, 0, 1, 1, 2]))
>>> z,x,y = d.nonzero()

If we wanted a more complicated cut, we could have done something like (d > 3.4).nonzero() or something, as True has an integer value of 1 and counts as nonzero.

Finally, we plot:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, -z, zdir='z', c= 'red')
plt.savefig("demo.png")

giving