Support fancy indexing

Start with NDCube.axis_world_coords and NDCube.wcs.world_to_array_indices, then review sunpy/map/maputils.py and its pixelate_coord_path workflow. Implement the decided option 1 using tabular coordinate transformations for the tuple of slices and index arrays, returning a new NDCube that supports the motivating path-extraction workflow.

Describe the feature

Provide a function that takes an NDCube and a tuple of slice items and/or arrays of array indices and returns a new NDCube. As fancy indexing is not suitable for all underlying data array representations (e.g. dask), this should be a function rather than a method on NDCube.

Proposed solution

The bulk of the work in developing this functionality is building the new WCS. Fancy indexing enables discontinuities in the coordinate transformations. We foresee the following paths forward:

1. Use NDCube.axis_world_coords to get tables of all the coordinates, and then produce a tabular gWCS object using the fancy index inputs. Do this for all axes, even if the original transformation could be maintained in some cases. This makes the implementation easier, but inefficient for large resultant cubes.
2. Only produce tabular transformations for the axes that are fancy indexed. Use CompositeWCS to combine the fancy indexed tabulated transformations with the those unaffected by the fancy indexing. CompositeWCS is slow, but may be more memory efficient in some cases.
3. Develop a new WCS implementation that wraps a FITS-WCS within a gWCS model. @Cadair thinks this is a tractable amount of work, but is probably a longer term goal.

Decision

For the first version of this functionality, pursue option 1.

Motivating Use Case

Ability to produce N-D time/frequency/polarisation-distance plots. This functionality will enable this to be acheived by the following work-flow:

• Use sunpy.map.maputils.pixelate_coord_path to get the world coordinates of all pixels along a path in the spatial plane of an NDCube.
• Use NDCube.wcs.world_to_array_indices to get the array indices of the path
• Build a slice item tuple valid of the cube, e.g. (path_indices_x, path_indices_y, slice(None),...) (assuming the spatial dimensions are the first two array axes)
• Extract pixels along path, retaining other dimensions in new (N-1)-D cube: new_cube = fancy_index_ndcube(old_cube, (path_indices, slice(None),...))