Hi, I'm getting a bit lost in numpy's fancy indexing, using argsort on a slice of a 4d arrays, and then using the resulting indices to sort other slices of that same array, preferably in one go. I have a data[d, b, y, x] array (for the curious: date, band, y, x: satellite image) Then I create a sorting index from band 0: idx = data[:, 0, :, :].argsort(axis=0) Now I want to sort all 7 bands accoring to this index, without looping, or at least avoid looping as far as possible. Something like sortedData = data[idx], but then for each 'b'. However, I cant get multidim indexing using the argsort result to work. Any hints? Hard to explain. I hope someone understands it, otherwise please let me know and I'll try to refrase :) Cheers, Vincent Schut.