Make reproject_to_healpix Compatible with Multiple images with the same WCS

[parsotat]

reproject_to_healpix can be made faster if it allows for the same input format as retroject_interp, where states:

If the data array contains more dimensions than are described by the input header or WCS, the extra dimensions (assumed to be the first dimensions) are taken to represent multiple images with the same coordinate information. The coordinate transformation will be computed once and then each image will be reprojected, offering a speedup over reprojecting each image individually.

I have multiple images, formatted in a numpy array where each index of the first array corresponds to different energy ranges, in which each image has the same WCS. I currently have to iterate over the first axis to get the HEALPIX map for each energy range and then combine all the maps appropriately. Since I know that these images all correspond to the same WCS, it would be ideal to save the computation of the coordinate transformation.

To be explicit my code looks something like this:

new_array = np.zeros((energy_nbins, healpix_npixels))
for e in range(energy_nbins):
    array, footprint = reproject_to_healpix(
                        (image[e, :, :], wcs), coordsys,
                        nside=healpix_nside)
   new_array[e, :] = array

And it would be great for the reprojection call to be:

new_array= reproject_to_healpix((image, wcs), coordsys, nside=healpix_nside, footprint=False)

Currently, my version of reproject is 0.14.0.