Fix bug with artifacts in reproject_to_healpix

reproject/healpix/core.py

@@ -2,6 +2,8 @@
 from astropy.coordinates import SkyCoord
 from astropy_healpix import HEALPix, npix_to_nside
 
+from reproject.array_utils import map_coordinates
+
 __all__ = ["healpix_to_image", "image_to_healpix"]
 
 ORDER = {}
@@ -121,7 +123,6 @@ def image_to_healpix(data, wcs_in, coord_system_out, nside, order="bilinear", ne
         no coverage or valid values in the input image, while values of 1
         indicate valid values.
     """
-    from scipy.ndimage import map_coordinates
 
     hp = HEALPix(nside=nside, order="nested" if nested else "ring")
 
@@ -134,13 +135,19 @@ def image_to_healpix(data, wcs_in, coord_system_out, nside, order="bilinear", ne
     world_out = SkyCoord(lon_out, lat_out, frame=coord_system_out)
 
     # Look up pixels in input WCS
-    yinds, xinds = wcs_in.world_to_pixel(world_out)
+    xinds, yinds = wcs_in.world_to_pixel(world_out)
 
     # Interpolate
 
     if isinstance(order, str):
         order = ORDER[order]
 
-    healpix_data = map_coordinates(data, [xinds, yinds], order=order, mode="constant", cval=np.nan)
+    healpix_data = map_coordinates(
+        data,
+        np.array([yinds, xinds]),
+        order=order,
+        mode="constant",
+        cval=np.nan,
+    )
 
     return healpix_data, (~np.isnan(healpix_data)).astype(float)

reproject/healpix/tests/test_healpix.py

@@ -155,7 +155,7 @@ def test_reproject_to_healpix_input_types(valid_celestial_input_data):
 
     # Make sure there are some valid values
 
-    assert np.sum(~np.isnan(healpix_data_ref)) == 3
+    assert np.sum(~np.isnan(healpix_data_ref)) == 4
 
     np.testing.assert_allclose(healpix_data_ref, healpix_data_test)
     np.testing.assert_allclose(footprint_ref, footprint_test)
@@ -180,3 +180,31 @@ def test_reproject_from_healpix_output_types(valid_celestial_output_projections)
 
     np.testing.assert_allclose(output_ref, output_test)
     np.testing.assert_allclose(footprint_ref, footprint_test)
+
+
+def test_reproject_to_healpix_exact_allsky():
+
+    # Regression test for a bug that caused artifacts in the final image if the
+    # WCS covered the whole sky - this was due to using scipy's map_coordinates
+    # one instead of our built-in one which deals properly with the pixels
+    # around the rim.
+
+    shape_out = (160, 320)
+    wcs = WCS(naxis=2)
+    wcs.wcs.crpix = [(shape_out[1] + 1) / 2, (shape_out[0] + 1) / 2]
+    wcs.wcs.cdelt = np.array([-360.0 / shape_out[1], 180.0 / shape_out[0]])
+    wcs.wcs.crval = [0, 0]
+    wcs.wcs.ctype = ["RA---CAR", "DEC--CAR"]
+
+    array = np.ones(shape_out)
+
+    healpix_array, footprint = reproject_to_healpix(
+        (array, wcs),
+        coord_system_out="galactic",
+        nside=64,
+        nested=False,
+        order="bilinear",
+    )
+
+    assert np.all(footprint > 0)
+    assert not np.any(np.isnan(healpix_array))