Add option to instantiate an NDCube with array indices in pixel order

Start at NDCube instantiation and review the linked proof-of-concept branch for how an index_order="pixel" option changes array shape, axis metadata, and slicing. Define the comprehensive behavior needed beyond the examples, then verify that pixel-ordered cubes support the demonstrated construction, indexing, and plotting workflows.

Describe the feature

I find it maddening that NDCube bakes in that the array indices are in the reverse order of the WCS pixel axes. This is a natural consequence of the natural (row-major) ordering of NumPy arrays, where the array indices start with the largest strides. However, that means a user will always have to think about whether the API is asking for input in array-index ordering or in WCS-pixel ordering.

Take this cube from the documentation:

import astropy.wcs
import matplotlib.pyplot as plt
import numpy as np

from ndcube import NDCube

# Define data array.
data = np.random.rand(3, 4, 5)
# Define WCS transformations in an astropy WCS object.
wcs = astropy.wcs.WCS(naxis=3)
wcs.wcs.ctype = 'WAVE', 'HPLT-TAN', 'HPLN-TAN'
wcs.wcs.cunit = 'Angstrom', 'deg', 'deg'
wcs.wcs.cdelt = 0.2, 0.5, 0.4
wcs.wcs.crpix = 0, 2, 2
wcs.wcs.crval = 10, 0.5, 1
wcs.wcs.cname = 'wavelength', 'HPC lat', 'HPC lon'

# Now instantiate the NDCube
my_cube = NDCube(data, wcs=wcs)

Wavelength is the first WCS pixel axis, and is indexed in the first position in the attached WCS:

>>> my_cube.wcs.axis_type_names[0]
'wavelength'

but for natural NumPy ordering, wavelength comes last, so it is indexed in the third position in NDCube's own API:

>>> my_cube.array_axis_physical_types[2]
('em.wl',)
>>> my_cube[:, :, 0]  # extract the 2D HPC map for the first wavelength
<ndcube.ndcube.NDCube object at 0x00000191F8D6D150>
NDCube
------
Shape: (3, 4)
Physical Types of Axes: [('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon'), ('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon')]
Unit: None
Data Type: float64
>>> my_cube[:, :, 0].plot()
>>> plt.show()

Proposed solution

I propose that a keyword argument be added to NDCube instantiation that enables having the array indices be in the same order as the WCS pixel axes. For such a NDCube, the user doesn't have to think about which part of the API uses which ordering.

Here is an example using this proof-of-concept branch to create a new version of the above cube:

>>> new_cube = NDCube(data.T, wcs=wcs, index_order="pixel")  # new keyword argument
>>> new_cube.shape
(5, 4, 3)
>>> new_cube.array_axis_physical_types[0]  # wavelength is now indexed in the first position
('em.wl',)
>>> new_cube[0, :, :]
<ndcube.ndcube.NDCube object at 0x0000021C2FFECE90>
NDCube
------
Shape: (4, 3)
Physical Types of Axes: [('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon'), ('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon')]
Unit: None
Data Type: float64
>>> new_cube[0, :, :].plot()
>>> plt.show()

(This proof of concept is far from comprehensive, but it is able to do the above.)