WCS ヘッダーを含む FITS ファイルがあると仮定すると、次のことができます。
#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 には、銀河座標の全空のごく一部に対応する 2D numpy 配列 (NY x NX) が含まれます。
Mollweide プロジェクションを使用して次の散布図にそのスカイマップをオーバーレイしたい場合、その 2d numpy を WCS と呼ばれる healpix 形式に変換する最良の方法は何でしょうか?
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()
ここで、経度と緯度は 2 つの 2D 配列 (numpy.meshgrid で生成) であり、関心のあるいくつかの天体物理ソースの銀河座標が含まれています。
healpix ピクセルを座標に変換して、スカイマップで使用できるものと一致させ、そこから値を補間することもできると思いますが、もっと簡単でエレガントで正確なものがあるはずですよね?