astropy / reproject

Python-based Astronomical image reprojection :milky_way: - maintainer @astrofrog
https://reproject.readthedocs.io
BSD 3-Clause "New" or "Revised" License
105 stars 64 forks source link

Nans appearing in reproject_to_healpix #284

Open runburg opened 2 years ago

runburg commented 2 years ago

I might be misunderstanding something, but I find that when using reproject_to_healpix, the returned healpy map has exactly nside number of NaN pixels. I'm trying to reproject the Fermi galactic background map and I persistently get NaN pixels in the reprojection where there are none before the projection.

I can add an MWE and a sample file if needed.

astrofrog commented 2 years ago

A MWE would be helpful, thanks!

runburg commented 2 years ago

The Fermi background file is here: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html

from astropy.io import fits
from astropy.wcs import WCS
import reproject

# this is the Fermi galactic background map
gal_bg = fits.open('gll_iem_v07.fits')

wcs = WCS(gal_bg[0].header).dropaxis(2)

hpx_map, footprint = reproject.reproject_to_healpix((gal_bg[0].data[0], wcs), 'galactic', nside=nside)

Here hpx_map has nside number of NaNs and footprint has nside number of 0s. Let me know if this is unclear or doesn't work!

lfarinaa commented 5 months ago

A kind reminder that this is not fixed. These pixels correspond to the nside pixels that run along the meridian that bisects the rhomboid at the "back" of the healpix tesellation. See fig: image I would say that this could be a symptom that whatever backend is used by reproject_to_healpix (I think this would be map_coordinates) doesn't understand that the image wraps at -180/180 deg, and so the centers of these pixels appear to be beyond the domain of the image you are reprojecting.

As a stopgap, I'm re-interpolating from the neighboring reprojected pixels:

badPix=np.argwhere(np.isnan(hMap))
badPixNeigh=hp.pixelfunc.get_all_neighbours(nsideGalFlux,badPix)
for i in np.arange(0,nside):
    hMap[badPix[i]]=np.nanmean(hMap[badPixNeigh[:,i]])

(Renormalize afterwards) Obviously, this in not a great solution. It would be great if someone could take a quick look at this.