astropy / astropy-healpix

BSD-licensed HEALPix for Astropy - maintained by @astrofrog and @lpsinger
https://astropy-healpix.readthedocs.io
BSD 3-Clause "New" or "Revised" License
50 stars 22 forks source link

Document HEALPix projection with Astropy #99

Open cdeil opened 6 years ago

cdeil commented 6 years ago

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:

  1. How can I use astropy.wcs.WCS to achieve the HEALPix projection?
  2. Why does the HEALPix projection yield values in the range X = {0, 360} and Y={-90, 90}, i.e. values in the same range as LON = {0, 360} and LAT={-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).

astrofrog commented 6 years ago

I've always been confused by the HPX projection too, so I agree some clarification docs would be good.