I'm thinking this should be

dout = self.view(type=ndarray)[indx]

because it is possible that the _baseclass overrides the indexing.

I see you logic (especially after the troubles with masked matrices for .flat...; indeed, cc @seberg) but one of the main points of this PR is actually that, to allow the baseclass to have its say in getting and setting. For astropy's Quantity, getting a proper slice would work fine, since that gets turned back into a MaskedArray and is updated (which is important as it copies the Quantity's unit). But if one gets a single item, ndarray's __getitem__ leaves one with an np.float64 object, which, if not masked, will be returned as is, i.e., one has lost the unit (or, more generally, removed the baseclass aspect of the item).

Indeed, the above is why we have a Quantity.__getitem__ : we return a single item as a a̶r̶r̶a̶y̶ ̶s̶c̶a̶l̶a̶r̶ 0-dimensional array instead (and have an internal property that remembers that this is a scalar).

Not sure what to think of this, except that it seems like this "masked matrix" stuff was always broken with indexing... It likely makes sense to actually change this, the finalization stage should attempt to fix the mask shape (needs test with matrix). But shape fixing in the data array and not having the same thing in the mask array, still seems like asking for trouble to me in any case (i.e. the whole reason why the .flat was needed)

It does seem to make sense to allow the _baseclass to override the indexing. I did check scipy with your changes and didn't see anything the looked related to this. Hmm... could you add a note here and we'll give this a try.

Looks like the code in MaskedArray could use a cleanup, if only for clarity.

Curious that _data and data are the same thing, likewise _mask and mask. I wonder why that was done? Backward compatibility? Intent to deprecate one?

Masked array could probably use a rewrite using __numpy_ufunc__ anyway when it is time...

Isn't data a property and _data the real array? Though I think you are probably right and it is identical. If it doesn't, it would be better if .data would return a view on _data and not the object itself.

About the matrix stuff/subclasses that change the shape... This makes enough trouble in numpy, getting it to work right in ma is probably hell and buggy at best at this time (just my guess on things, I never even used masked arrays ^^). For what I care, I would be fine with deprecating all that .flat usage. Do we really need to support subclasses that mess with the shape inside masked arrays?

@charris - It is indeed odd that _data and data are the same thing. For _mask and mask it looks more like the usual private property shielded by a getter and setter.
@seberg - definitely would be great to move to __numpy_ufunc__... but it does need someone to sit down and be careful...

Y'all, sorry for not being able to be more responsive on all these problems. The _data is indeed the private property accessing the underlying ndarray (just like _mask stores the boolean ndarray corresponding to, you guessed it). Neither _data or _mask are supposed to be accessible to the lay person, hence the data and mask properties. Well, hence the mask: the data was introduced to give a quick access to the ndarray in cases where you'd want to drop the MaskedArray part for a while, perform faster computations and pull the mask back at the end.
In the problem at hand, the objective is to get the underlying data in its 'least-adulterated' form: @charris is right, it should be `dout=self.view(...)[indx].

@pierregm - the problem with doing dout=self.view(np.ndarray)[indx] is that, as mentioned above, this presumes that a single item from a subclass can be represented as a np.float64 scalar. This, however, is false at least for astropy's Quantity, where one does need to keep the unit.
Note that the problem is worse for __setitem__: the baseclass may not allow one to set parts to just anything, it may need to check that, e.g., the data put in has the same class (or, for Quantity, that the units are consistent, and convert if necessary).

OK, this line has pretty much everything the bothers me ;) Note the here _data is viewed as _baseclass, which may, but I'm not sure, have a different shape than self._data in the matrix case or some other odd case. Same thing above, where ndarray.view(self, ndarray) may be different from self.data, the latter may actually be a matrix type. I'd also like to know what the difference between self.data and self._data is, likewise self.mask and self._mask.

Dealing with masked arrays is like trying to sneak up on a tiger while creeping across eggshells. Disaster threatens at every move.

@charris - yes, it is weird, but it may be that this predates the following lines

    def _get_data(self):
        """Return the current data, as a view of the original
        underlying data.

        """
        return ndarray.view(self, self._baseclass)
    _data = property(fget=_get_data)
    data = property(fget=_get_data)

Similarly,

    def _get_mask(self):
        """Return the current mask.

        """
        # We could try to force a reshape, but that wouldn't work in some cases.
#        return self._mask.reshape(self.shape)
        return self._mask
    mask = property(fget=_get_mask, fset=__setmask__, doc="Mask")

So, this would seem to be safe... Or is there a MaskedMatrix class or so where parts get overwritten? My larger worry is that this was all written at a time when people were not as used to writing lots and lots of test cases...

Turns out the code actually postdates those functions you quote by about two years. Two different authors though, so we may not be the first struggling with the code ;) The mask part does look safe, although I'd like to know what the comment on line 2006 means.

I don't feel very confidence reviewing masked array code, there are a lot of subtleties that I don't understand yet...