It's easier to understand what np.vstack, np.hstack and np.dstack* do by looking at the .shape attribute of the output array.

Using your two example arrays:

print(a.shape, b.shape)
# (3, 2) (3, 2)

• np.vstack concatenates along the first dimension...

  print(np.vstack((a, b)).shape)
  # (6, 2)
  

• np.hstack concatenates along the second dimension...

  print(np.hstack((a, b)).shape)
  # (3, 4)
  

• and np.dstack concatenates along the third dimension.

  print(np.dstack((a, b)).shape)
  # (3, 2, 2)
  

Since a and b are both two dimensional, np.dstack expands them by inserting a third dimension of size 1. This is equivalent to indexing them in the third dimension with np.newaxis (or alternatively, None) like this:

print(a[:, :, np.newaxis].shape)
# (3, 2, 1)

If c = np.dstack((a, b)), then c[:, :, 0] == a and c[:, :, 1] == b.

You could do the same operation more explicitly using np.concatenate like this:

print(np.concatenate((a[..., None], b[..., None]), axis=2).shape)
# (3, 2, 2)

* Importing the entire contents of a module into your global namespace using import * is considered bad practice for several reasons. The idiomatic way is to import numpy as np.

After some joint effort with otterb, we concluded that preallocating of the array is the way to go. Apparently the performance killing bottleneck was the array layout with the image number (n) being the fastest changing index. If we make n the first index of the array (which will default to the "C" ordering: first index changest slowest, last index changes fastest) we get the best performance:

from time import time
import numpy as np

# Create 100 images of the same dimention 256x512 (8-bit). 
# In reality, each image comes from a different file
img = np.random.randint(0,255,(100, 256, 512))

# counter-intuitive but preallocation is slightly slower
stacked = np.empty((100, 256, 512))
t0 = time()
for n in range(100):
    stacked[n] = img[n]
print time()-t0  
print stacked.shape

There are two possible reasons for the error:

1. I have a feeling that your script may not be selecting the correct bands. You should sort the input bands using sorted(glob.glob("*.TIF")) before converting them into numpy arrays.
2. You should rescale your bands / numpy arrays to uint8 format before plotting.

Here is an example that worked for me:

import glob
import numpy as np
import gdalnumeric
import matplotlib.pyplot as plt

BandList = glob.glob('*.TIF')
# returns: ['B4.TIF', 'B6.TIF', 'B3.TIF', 'B5.TIF', 'B1.TIF', 'B7.TIF', 'B2.TIF']

BandList = sorted(BandList)
# returns: ['B1.TIF', 'B2.TIF', 'B3.TIF', 'B4.TIF', 'B5.TIF', 'B6.TIF', 'B7.TIF']

RGB_Image = np.dstack((gdalnumeric.LoadFile(BandList[3]), gdalnumeric.LoadFile(BandList[2]), gdalnumeric.LoadFile(BandList[1]))).astype('uint8')

plt.figure(figsize=(3, 3))
plt.imshow(RGB_Image, interpolation='nearest')
plt.savefig('TestImage.png', format='png', dpi=300, bbox_inches='tight')
plt.close()

And here is the output image: