CODEWITHSUNDEEP
pythonnumpymatplotlibastropyhealpy
Mireia
Mireia

Reputation: 11

FITS with WCS to healpy

Assume I have a FITS file that contains a WCS header, so that I can do:

#import healpy as hp
#import astropy.io.fits as pyfits
#from astropy.wcs import WCS

listofhdus = pyfits.open(FITS)
wcs = WCS(listofhdus[0].header)

listofhdus[0].data would contain a 2D numpy array (NY times NX) which correspond to a small portion of the full sky in Galactic coordinates.

What would be the best way to convert that 2d numpy into healpix format, known the WCS, if I want to include overlay that skymap into the following scatter plot with Mollweide projection?

NPIX = hp.nside2npix(512)
m = 0*np.arange(NPIX)
hp.mollview(m, title="test")
hp.projscatter(longitude, latitude, lonlat=True, 
               coord='G',marker='+',color='black')
hp.graticule()

where longitude and latitude are two 2D arrays (generated with numpy.meshgrid) with Galactic coordinates of some astrophysical sources I am interested on?.

I guess I could try to translate my healpix pixels into coordinates, somehow match them with the ones available in my skymap and interpolate the values from there, but there must be something easier, more elegant and precise, right?

Upvotes: 1

Views: 705

Answers (1)

Daniel Lenz
Daniel Lenz

Reputation: 3877

You should check out the reproject, which takes care of all your reprojection needs to/from HEALPix.

In particular, check out reproject_to_healpix.

If you need a better control of your interpolation routine (using a convolution-based gridder, ensuring flux conservation, performance), I'd suggest you check out cygrid, which has some examples for HEALPix as well.

Upvotes: 1

Related Questions