astropy-healpix
astropy-healpix copied to clipboard
Document HEALPix projection with Astropy
I remember that when first encountering HEALPix I was very confused, because it's a pixelisation and a WCS projection, and I couldn't find docs how to work with the HEALPix projection in Astropy.
Specifically, this doesn't apply the projection, and doesn't give an error, it just passes the values through without change:
>>> from astropy.wcs import WCS
>>> wcs = WCS(key='HPX')
>>> print(wcs)
WCS Keywords
Number of WCS axes: 2
CTYPE : '' ''
CRVAL : 0.0 0.0
CRPIX : 0.0 0.0
PC1_1 PC1_2 : 1.0 0.0
PC2_1 PC2_2 : 0.0 1.0
CDELT : 1.0 1.0
NAXIS : 0 0
>>> print(wcs.wcs_world2pix(42, 43, 1))
[array(42.), array(43.)]
>>> print(wcs.wcs_pix2world(42, 43, 1))
[array(42.), array(43.)]
I did manage to access the projection via the astropy.modeling.projections.Sky2Pix_HEALPix wrapper for that WCS transform:
>>> from astropy.coordinates import Angle
>>> from astropy.modeling.models import Sky2Pix_HEALPix, Pix2Sky_HEALPix
>>> x, y = Sky2Pix_HEALPix()(-148, 65)
>>> print(x, y)
-141.89216598555333 66.14250235769997
>>> print(Angle(x, 'deg').rad, Angle(y, 'deg').rad)
-2.4764854792342104 1.1544044416499766
>>> lon, lat = Pix2Sky_HEALPix()(x, y)
>>> print(lon, lat)
-147.99999999999997 65.0
This is the example from here and the results are identical.
I have two questions:
- How can I use
astropy.wcs.WCS
to achieve the HEALPix projection? - Why does the HEALPix projection yield values in the range
X = {0, 360}
andY={-90, 90}
, i.e. values in the same range asLON = {0, 360}
andLAT={-90, 90}
, making it easy to think(X, Y)
would be angles again (which they aren't)?
For the second question, I think it's a bit of a convention to have X=LON
, another scale factor, e.g. X=LON/360
could have been chosen? And for Y
, is there a fundamental reason why it has the same range as LAT
, or is is a bit by chance that it comes out to Y={-90, 90}
?
Once these two points are clarified, I'd like to make a small docs addition in astropy-healpix
to clarify this, and give an example how one can do the HEALPix projection with Astropy (since it's not exposed in the astropy-healpix API).