Artefacts appearing during reprojecting back to healpix

[Wasim04]

Hi,

I was trying to use the reproject_to_healpix function. It did reproject the distribution I anticipated. However, a 'chain' like '--' line artefact appeared. Please let me know how to remove this. following is the small script, to reproduce the issue.

nside = 64
npix = hp.nside2npix(nside)
healpix_map = np.arange(npix)

# Define the output WCS (World Coordinate System) for the flat 2D projection
shape_out = (160, 320)
wcs = WCS(naxis=2)
wcs.wcs.crpix = [shape_out[1] / 2, shape_out[0] / 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"]

# testMap is a probability distribution map. shaped as (1, 49152).

# Project HEALPix map to a flat 2D image
flat_image, footprint = reproject_from_healpix(
    (testMap[0], 'galactic'), 
    wcs, 
    shape_out=shape_out,
    nested = True,
    order='bilinear'
)

plt.imshow(flat_image, origin='lower')
plt.title('Flat 2D Projection')
plt.colorbar()
plt.show()

# Reproject the flat image back to a HEALPix map
healpix_map_reprojected, _ = reproject_to_healpix(
    (flat_image, wcs), 
    coord_system_out='galactic', 
    nside=nside,
    nested = True,
    order='bilinear' 
)

hp.mollview(healpix_map_reprojected, title="Reprojected HEALPix Map")
plt.show()

[Wasim04]

Hi @astrofrog, Kindly accept my apology for tagging you. I am wondering if you could kindly check the above problem and advise me a solution.

Many Thanks Wasim

[astrofrog]

@Wasim04 - sorry for the delay in getting back to you, would you be able to share the input file? (for example via a dropbox or google drive or other link)

[Wasim04]

Hi Thomas, Thanks for your response. Here is an example input file. It's just a nside 64 probability map. https://github.com/Wasim04/test_repo

Many Thanks Wasim

[astrofrog]

@Wasim04 - the artifacts are due to the fact some of the output HEALPIX cells don't have a match in the input image (and are therefore set to NaN) - you can make the artifacts disappear by changing the reproject_to_healpix call with:

healpix_map_reprojected, footprint = reproject_to_healpix(

to capture the footprint, then e.g.

healpix_map_reprojected[footprint == 0] = 0

Note that I had to change nested to False to get reasonable results in mollview.

Why are some pixels not overlapping in the first place? I think the definition of the WCS is slightly wrong, as for example if I try and get the coordinates of the corner of the image the first value should be 180,-90 and the second should be defined, not NaN:

In [13]: print(wcs.pixel_to_world(-0.5, -0.5))
    ...: print(wcs.pixel_to_world(319.5, 159.5))
<SkyCoord (ICRS): (ra, dec) in deg
    (179.4375, -89.4375)>
<SkyCoord (ICRS): (ra, dec) in deg
    (nan, nan)>

I think the crpix line should be:

wcs.wcs.crpix = [(shape_out[1]+1) / 2, (shape_out[0]+1) / 2]

which yields sensible results:

In [19]: print(wcs.pixel_to_world(-0.5, -0.5))
    ...: print(wcs.pixel_to_world(319.5, 159.5))
<SkyCoord (ICRS): (ra, dec) in deg
    (180., -90.)>
<SkyCoord (ICRS): (ra, dec) in deg
    (180., 90.)>

I still get artifacts in this case, though not if I change coord_system_out='icrs' interestingly. I will need to investigate this further.

[astrofrog]

This should be fixed with https://github.com/astropy/reproject/pull/459