Provide entrypoint for specifying alternate versions of `map_coordinates`

`reproject_interp` relies on `scipy.ndimage.map_coordinates` for interpolating the image to the new pixel array. This is potentially quite slow for large images/complicated WCSs. However, there are several alternate implementations of `map_coordinates` that could provide potential speedups, e.g.

- `cupy` implementation of `map_coordinates` when your array is on a GPU: https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.ndimage.map_coordinates.html
- In-progress `dask-image` implementation of `map_coordinates`: https://github.com/dask/dask-image/pull/237

I see that `map_coordinates` is already abstracted away from `scipy.ndimage.map_coordinates` in some sense by placing it in its own function in `array_utils`. Is there a sensible way for providing an entrypoint for swapping out the underlying `map_coordinates` function? One possibility of course would be just passing in a function as an argument from the high-level `reproject_interp` interface all the way down to `map_coordinates`. Is there a cleaner way?

Start by reading the map_coordinates wrapper in array_utils and tracing how reproject_interp reaches it. Compare the existing abstraction with the CuPy and dask-image implementations linked in the issue. Done means a documented, maintainable entrypoint for selecting an alternate implementation without breaking the existing SciPy path.

Written by the indexing model from the issue text.