DOC: warn about the indexing='xy' default in np.meshgrid

The default of np.meshgrid is indexing='xy'. This is not intuitive - intuitive would be: indexing='ij'. I have been bitten by this default behavior many times since I started using numpy in 2014 and as recently as a month ago. In the last couple of weeks, I've already answered two Stackoverflow questions that related to this 1 and 2.

Reproducing code example:

import numpy as np
print(np.meshgrid(np.arange(2), np.arange(3), np.arange(4))[0].shape)

Generates output

(3, 2, 4)

Most people would expect (2, 3, 4) unless they read the documentation very carefully. The documentation of the indexing parameter:

indexing{‘xy’, ‘ij’}, optional -- Cartesian (‘xy’, default) or matrix (‘ij’) indexing of output. See Notes for more details.

is not helpful in my opinion. To this day, I don't really get what it the options have to do with "cartesian" or "matrix"; all I can remember is: "'ij' is the one I need.". I understand that the default needs to stay there for compatibility and that the arbitrary "swap the order of the first two axes" is likely inherited from how Matlab does it. But at least, the documentation could be more clear, for example:

indexing{‘xy’, ‘ij’}, optional -- Cartesian (‘xy’, default) or matrix (‘ij’) indexing of output. With 'xy' indexing, the order of the first two axes will be swapped. See Notes for more details.

It could even be considered to provide an alternative function that behaves like np.meshgrid(indexing='ij'). I know that np.mgrid sort of does, but its syntax is rather esoteric.

Numpy/Python version information:

1.18.1 3.7.6 (default, Jan 8 2020, 19:59:22)
[GCC 7.3.0]

Start with the np.meshgrid documentation, especially the indexing parameter description and its Notes section. Clarify that the default 'xy' indexing swaps the first two axes, while preserving the existing default and explaining the difference from 'ij'; the documentation should make the resulting axis order understandable.