Implementation of Gaussian Mixture Model for clustering when dealing with multidimensional hyperspectral data in python
If I understand your type of data correctly, what you have is essentially an image with $800\times800 =640000$ pixels and each pixel has a vector of 4 values.

If this stands, I suppose you could then transform your data to a $640000\times4$ matrix, so as to conform with scikit-learn's data representation schema of inputting matrices of shape ($\#samples\times\#features$) and then you could use the GMM class implemented by the package.

And you could transform the data back if needed.

P.S.: On a second thought, I don't know if the spatial distribution of pixels plays a role on the generation of the GMM components. Maybe try it out and let us know?

Ok, after you have trained your model you can ask for it to predict the mixture model from which each pixel in your dataset has been generated.

That is if you run again:
$clusters=gmm.predict(newdata)$ you will get an array with values $0-2$, denoting the cluster that each pixel belongs to. You can, then map these values in the image at hand, applying a mask of different color for each pixel in order to see if your model correctly identifies the three classes.

Here is an example of what I mean:

from scipy.misc import imread, imshow
from sklearn.datasets import load_sample_image

X = load_sample_image('china.jpg')
old_shape = X.shape
X = X.reshape(-1,3)
gmm = mixture.GMM(covariance_type='full', n_components=5)
gmm.fit(X)
clusters = gmm.predict(X)
clusters = clusters.reshape(old_shape[0], old_shape[1])
imshow(clusters)