It's now possible to index dask arrays using advanced Numpy indexing, which means one can do e.g.:

dask_array.vindex[ipix, jpix]

so this means that for the nearest neighbor method in reproject_interp we could just do this directly on input dask arrays. Similarly, for bilinear interpolation one would need to fetch four indices for each output pixel, but it seems doable.

We should investigate to see if the performance implications of this are prohibitive, but as long as it's not too slow it might be a good way and the only way to avoid dumping gigantic arrays to disk.

If we could do this it would make things easier for reprojecting from HiPS to standard projections.