Add HEALPIX reprojection #27

keflavich · 2014-01-07T16:13:56Z

Tool to convert HEALPIX images to Galactic or ICRS coordinates, based tightly on @cdeil's function in gammapy. See healpy/healpy#129.

cdeil · 2014-01-07T18:46:05Z

@keflavich Thanks for writing this!

The most common use case is going to be that the user has a FITS file with a HEALPIX HDU.

How about adding a second healpix_hdu_to_image function so that the user doesn't have to read COORDSYS and ORDERING from the header and pass appropriate hpx_coord_system and nest parameters?

Or if it's easy to make a HEALPIX HDU from given data, hpx_coord_system and nest we could only have the healpix_hdu_to_image function?

keflavich · 2014-01-10T08:21:45Z

@cdeil - added, in case you didn't see it. I kept the other functions because I like having different levels of access - no sense creating duplicate HDUs in memory if you've already loaded one (especially with nsides=2048 healpix, where each copy is >1 GB).

@cdeil - any thoughts on generating tests for this? We shouldn't include data files, but I don't have the familiarity with healpix needed to know what would make a good test.

cdeil · 2014-01-10T08:32:15Z

Thanks for adding the extra versions!

About tests: I don't think we need any real input files to test reprojection, we can always generate test images that are constant or with a checkerboard or sin pattern or something.
I barely know anything about HEALPIX and reprojection, but maybe testing this would be useful:

Sample input / output image at some positions (values should be identical up to interpolation / re-gridding errors)
Integrate over some regions (e.g. circles or boxes or full image) (values should be identical up to integration / interpolation / re-gridding errors).

But I think the main thing to test is that the main use cases work, not to test how good the reprojection methods are numerically.

(sorry I don't have time to work on @python-reprojection@ at the moment)

cdeil · 2014-01-21T16:48:48Z

@keflavich Do you have time to add some basic tests?
If not I'll give it a go because I need this and would like to get it into master soon.

keflavich · 2014-01-21T16:52:09Z

Not at the moment, sorry - go ahead and add them if you can.

astrofrog · 2014-07-16T21:55:10Z

So just to be clear, WCSLIB does include the HPX projection now - but does this do something different?

keflavich · 2014-07-16T21:59:37Z

I have no idea. I have no recollection of writing this! And it was only 6 months ago...

astrofrog · 2014-07-16T22:02:17Z

I guess that this still implements reprojection, but I'm wondering whether it duplicates stuff we would want in the main reprojection routine. We should separate out the healpix-specific stuff from the reprojection itself.

astrofrog · 2014-07-22T21:23:46Z

@keflavich @cdeil - do you have an example of a healpix tile we could use to try and see if it currently works with vanilla astropy.wcs?

cdeil · 2014-07-22T21:29:49Z

The healpy repo has example files here and the healpix repo has more test files here.

I won't have time to try this out before next week, so if someone else has time, please go ahead.

DrPhilEvans · 2014-11-18T08:36:57Z

Hi,
This is a really useful package, but I'm having some problems: the RA/Dec that I get out of the reprojection is wrong! I'm basically taking a HEALPIX map and trying to reproject it to an Aitoff projection, and save it. The code works, but if I plot it (in ds9, or using kapteyn) the RA/Decs are clearly wrong. Here are the original (mollview) and reprojected images (plotted in ds9)

And the code (minus imports) is below - can anyone tell me what I'm doing wrong?

Thanks!

naxis1=2000
naxis2=1000
fra=180.0
fdec=0.0
cdelt1 = 360. / naxis1
cdelt2 = 180 / naxis2
crpix1 = (naxis1 + 1) * 0.5
crpix2 = (naxis2 + 1) * 0.5
header = {'NAXIS': 2 ,
'NAXIS1': naxis1 ,
'NAXIS2': naxis2 ,
'CDELT1': cdelt1,
'CRPIX1': crpix1,
'CRVAL1': fra,
'CTYPE1': 'RA---AIT',
'CUNIT1': 'deg ',
'CDELT2': cdelt2,
'CRPIX2': crpix2,
'CRVAL2': fdec,
'CTYPE2': 'DEC--AIT',
'CUNIT2': 'deg',
'COORDSYS': 'ICRS'
}

healpix_data=hp.read_map('hpMap.fits.gz',nest=True)
rep=healpix.healpix_to_image(healpix_data,header,'ICRS', True)
hdu=fits.PrimaryHDU(rep)
for key, value in header.iteritems():
hdu.header[key]=value
hdulist=fits.HDUList([hdu])
hdulist.writeto('test.img', clobber=True)

astrofrog · 2014-11-18T09:14:36Z

@DrPhilEvans - I can look into it today. Could you share the input data with me? If so, can you send it by email to thomas.robitaille@gmail.com?

astrofrog · 2014-11-18T09:16:02Z

@DrPhilEvans - and just to be clear, are you relying on the code from this pull request? (presumably that is where healpix_to_image comes from?)

DrPhilEvans · 2014-11-18T09:18:49Z

@astrofrog I think so. I'm not familiar with Github (I use svn) so I wasn't sure how to pull this - in the end I downloaded the main "reproject" project, and then added the "healpix.py" file via the "Files changed" link at the top of this page.

astrofrog · 2014-11-18T10:02:08Z

@keflavich @cdeil - it would be great if we could brush up this code and add some tests. What I personally think would be nice is if for the high level interface we could just do this via the high level reproject function. That is, if I do:

reproject('somefile.fits', new_header)

and if somefile.fits is a healpix file (or if a healpix HDU is passed) it would seamlessly deal with it. This would prevent us from adding too many high level entry points. Now for converting from images to healpix we need to think about how that would work, but that's not implemented at the moment anyway, so I think we can worry about that later.

I looked into @DrPhilEvans's issue, but it relates entirely to code in this PR (and doesn't use any of the rest of the reproject package) so @keflavich or @cdeil it would be great if you could look into it if you have some time?

keflavich · 2014-11-18T10:04:37Z

Well, I agree with you that reproject should not care about the input format. I can't say I have any time, but I do have a project in mind that should use Planck all-sky data, so I'll see if I can reproject that..

...edit: wait, I wrote this PR!?! It's too easy to lose track of these things.

astrofrog · 2014-11-18T10:08:01Z

Thanks!

keflavich · 2014-11-18T10:42:15Z

I tried reprojecting this file:
http://irsa.ipac.caltech.edu/data/Planck/release_1/all-sky-maps/previews/HFI_SkyMap_217_2048_R1.10_nominal/index.html
(there's a .fits at that link)
to this header:

WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =              468.716 / Pixel coordinate of reference point
CRPIX2  =              1123.74 / Pixel coordinate of reference point
PC1_1   =    -0.00199999986216 / Coordinate transformation matrix element
PC2_2   =     0.00199999986216 / Coordinate transformation matrix element
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'GLON-CAR'           / galactic longitude, plate caree projection
CTYPE2  = 'GLAT-CAR'           / galactic latitude, plate caree projection
CRVAL1  =        3.65245551155 / [deg] Coordinate value at reference point
CRVAL2  =       0.605737090951 / [deg] Coordinate value at reference point
PV2_1   =                  0.0 / CAR projection parameter
LONPOLE =                  0.0 / [deg] Native longitude of celestial pole
LATPOLE =         89.394262909 / [deg] Native latitude of celestial pole
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
NAXIS   =                    2
NAXIS1  =                 3089
NAXIS2  =                 1628
COORDSYS= 'galactic'

and it looks fine (agrees with my other data).

I think the problem is likely the sign flip in latitude I introduced, which was necessary to get the reprojected image to agree with Galactic coordinates and is therefore correct for one coordinate system... but apparently not for the other? I'll investigate a little further.

DrPhilEvans · 2014-11-18T10:44:56Z

@keflavich Thanks.
FYI, the file I'm working on (at the moment) is one I downloaded from here:
http://www.ligo.org/scientists/first2years/2015/compare/12157/bayestar.fits.gz
If it helps.

keflavich · 2014-11-18T10:51:02Z

Well, I've uncovered some bugs, and when I tried to transform to a tangent projection, it just crashed on me:

WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =                  750 / Pixel coordinate of reference point
CRPIX2  =                  750 / Pixel coordinate of reference point
PC1_1   =    -0.00199999986216 / Coordinate transformation matrix element
PC2_2   =     0.00199999986216 / Coordinate transformation matrix element
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'RA---TAN'           / galactic longitude, plate caree projection
CTYPE2  = 'DEC--TAN'           / galactic latitude, plate caree projection
CRVAL1  =            266.55149 / [deg] Coordinate value at reference point
CRVAL2  =           -28.861962 / [deg] Coordinate value at reference point
PV2_1   =                  0.0 / CAR projection parameter
LONPOLE =                  0.0 / [deg] Native longitude of celestial pole
LATPOLE =         89.394262909 / [deg] Native latitude of celestial pole
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
NAXIS   =                    2
NAXIS1  =                 1500
NAXIS2  =                 1500
COORDSYS= 'icrs    '

Error encountered at /Users/adam/repos/Healpix_3.11/src/cxx/Healpix_cxx/healpix_base.cc, line 1178
(function void T_Healpix_Base<I>::get_interpol(const pointing&, fix_arr<I, 4ul>&, fix_arr<double, 4ul>&) const [with I = long int])

invalid theta value

libc++abi.dylib: terminating with uncaught exception of type PlanckError

so this is already getting me pretty far out of my depth.

cdeil · 2014-11-18T10:56:45Z

@keflavich I think healpy should be written such that it's impossible to make it crash with any input. Do you have a self-contained example I can try? If I see the same crash we should file a healpy issue.

cdeil · 2015-04-13T14:14:16Z

I keep getting emails asking for this functionality (because there used to be this old function that never quite worked in Gammapy). I once again don't have time to work on this and am not very good with HEALPIX in the first place.

Some more references that might be useful:

astrofrog · 2015-04-13T14:32:15Z

Sounds like something we could sprint on next week?

cdeil · 2015-04-13T14:40:54Z

Sounds like something we could sprint on next week?

Definitely!
(I'm not sure how useful I'll be concerning the core implementation of the interpolation routine, but we can probably look what the others have done and I can help setting up tests using files from http://www.atnf.csiro.au/people/mcalabre/WCS/ or elsewhere.)

lpsinger · 2015-04-24T15:28:34Z

What is the exact call that is crashing? I'm going to try to debug it...

lpsinger · 2015-04-24T15:37:29Z

reproject/healpix.py

+    >>> healpix_isnested = fits.getheader(healpix_filename,ext=1)['ORDERING'] == 'NESTED'
+    >>> reference_image = fits.open('reference_image.fits')[0]
+    >>> reprojected_data = healpix_to_image(healpix_data, reference_image, healpix_system, nest=healpix_isnested)
+    >>> fits.writeto('new_image.fits', reprojected_data, reference_image.header)


Problem here: by default, hp.read_map reorders the pixels to RING order. I'll make a pull request to fix this...

keflavich · 2015-04-24T15:38:47Z

@lpsinger - if I use the incorrect sign for lat, I get a crash. Otherwise, it works. However, I haven't tested for the RA/Dec one reported above.

This is my incomplete / WIP test code (note that it is not a functioning test and it is missing data):

from astropy.io import fits
from reproject import healpix

header = fits.Header.fromtextfile('hfi_output.hdr')
url = 'http://irsa.ipac.caltech.edu/data/Planck/release_1/all-sky-maps/maps/HFI_SkyMap_217_2048_R1.10_nominal.fits'
url = 'http://www.ligo.org/scientists/first2years/2015/compare/12157/bayestar.fits.gz'
out = healpix.healpix_reproject_file('/Users/adam/Downloads/HFI_SkyMap_217_2048_R1.10_nominal_ZodiCorrected.fits', header)

lpsinger · 2015-04-24T15:40:44Z

Is reproject set up to have its doctests run as part of setup.py test? Can I just put the corrected example code in the docstring?

keflavich · 2015-04-24T15:44:54Z

Yes, I think reproject runs doctests... but I can't prove it...

lpsinger · 2015-04-24T15:45:31Z

This is my incomplete / WIP test code (note that it is not a functioning test and it is missing data):

Can you provide me with a reference header?

lpsinger · 2015-04-24T17:03:28Z

Yes, I think reproject runs doctests... but I can't prove it...

No, it doesn't... I asked on the Astropy list how to turn this on. Anyway, I've gotten healpix_to_image working as of 38ff431. Still need to fix doctests and turn them on.

I am using this HEALPix test file:
http://www.ligo.org/scientists/first2years/2015/compare/12157/bayestar.fits.gz

Here's a plot of the original HEALPix image:

And here is my test script:

import astropy.io.fits
import astropy.wcs
from reproject.healpix import healpix_to_image
import healpy as hp
reference_header = astropy.io.fits.Header()
reference_header.update({
 'COORDSYS': 'icrs',
 'CDELT1': -0.4,
 'CDELT2': 0.4,
 'CRPIX1': 500,
 'CRPIX2': 400,
 'CRVAL1': 180.0,
 'CRVAL2': 0.0,
 'CTYPE1': 'RA---MOL',
 'CTYPE2': 'DEC--MOL',
 'CUNIT1': 'deg',
 'CUNIT2': 'deg',
 'NAXIS': 2,
 'NAXIS1': 1000,
 'NAXIS2': 800})
healpix_filename = 'bayestar.fits.gz'
healpix_data = hp.read_map(healpix_filename, verbose=False)
healpix_coordsys = 'C'
reprojected_data = healpix_to_image(healpix_data, reference_header, healpix_coordsys, nest=False, interp=False)
astropy.io.fits.writeto('new_image.fits', reprojected_data, reference_header)

And here a screen shot of the output image in DS9:

astrofrog · 2015-04-24T22:49:25Z

For the doctests, I think there's a doctest_plus option in setup.cfg? But before this all works properly, I think I'll have to merge in @cdeil's updates to the standard affiliated package files.

bsipocz · 2015-04-24T23:23:47Z

@astrofrog - it doesn't run atm, I can open a PR for it tomorrow.

astrofrog · 2015-04-25T09:36:09Z

@bsipocz - just to check, do you mean that you will open a PR to edit setup.cfg? Or fix https://github.com/astrofrog/reproject/pull/57?

astrofrog · 2015-04-25T09:40:58Z

@lpsinger - thanks for working on this! Could you open a PR to @keflavich's branch so that we can include the commits here?

DrPhilEvans · 2015-04-27T10:06:27Z

This is a minor point, but when I last used the reproject option to make
a full sky plot, my screen filled with error messages, because many of
the pixels in the FITS image have no valid coordinates (since the FITS
image is rectangular, but the aitoff/mollwiede/something projection is
oval). I have a feeling this can't be suppressed without actually
editing the healpy source code itself (I gave up and wrote my own C++
"HealPixToIMG" programme), but since the thread is currently active,
it's worth asking the question: can we get rid of these "error" messages?

On 25/04/15 10:40, Thomas Robitaille wrote:

@lpsinger https://github.com/lpsinger - thanks for working on this!
Could you open a PR to @keflavich https://github.com/keflavich's
branch so that we can include the commits here?

—
Reply to this email directly or view it on GitHub
https://github.com/astrofrog/reproject/pull/27#issuecomment-96169226.

Phil Evans,
Swift Development Scientist
X-ray and Observational Astronomy Group,
University of Leicester

Tel: +44 (0)116 252 5059
Mobile: +44 (0)7780 980240 (work)
Mobile: +44 (0)7974 977723 (personal)

pae9@leicester.ac.uk

http://www.star.le.ac.uk/~pae9
http://www.swift.ac.uk

Follow me as a Swift scientist on Twitter: @swift_phil
http://www.star.le.ac.uk/~pae9/twitter

lpsinger · 2015-04-27T14:59:59Z

@DrPhilEvans, did you try my branch? I'm not getting warning messages.

DrPhilEvans · 2015-04-27T15:01:31Z

Ah, no probably not. I got a bit lost navigating the different trees :S
Can you point me to the right place?

On 27/04/15 15:59, Leo Singer wrote:

@DrPhilEvans https://github.com/DrPhilEvans, did you try my branch?
I'm not getting warning messages.

—
Reply to this email directly or view it on GitHub
https://github.com/astrofrog/reproject/pull/27#issuecomment-96689199.

Phil Evans,
Swift Development Scientist
X-ray and Observational Astronomy Group,
University of Leicester

Tel: +44 (0)116 252 5059
Mobile: +44 (0)7780 980240 (work)
Mobile: +44 (0)7974 977723 (personal)

pae9@leicester.ac.uk

http://www.star.le.ac.uk/~pae9
http://www.swift.ac.uk

Follow me as a Swift scientist on Twitter: @swift_phil
http://www.star.le.ac.uk/~pae9/twitter

astrofrog · 2015-04-27T15:16:52Z

@lpsinger - just to let you know, when you open a pull request to @keflavich's repo, this is the repo to open a PR to: https://github.com/keflavich/python-reprojection (the repo here changed name after @keflavich forked it)

lpsinger · 2015-04-27T15:32:33Z

Ah, no probably not. I got a bit lost navigating the different trees :S
Can you point me to the right place?

$ git remote add lpsinger https://github.com/lpsinger/healpy.git
$ git fetch lpsinger
$ git checkout --track lpsinger/healpix

…ix, healpix_reproject_file Note that there will be some issues with round-trip projection (HEALPix->WCS->HEALPix) if some of the pixels lie precisely on the WCS map boundary.

Working and tested implementation of healpix_to_image, image_to_healpix, and healpix_reproject_file

astrofrog · 2015-04-29T19:10:53Z

@keflavich - can you change line 96 in .travis.yml from:

# - if [[ $SETUP_CMD != egg_info ]]; then $PIP_INSTALL ...; fi

to

- if [[ $SETUP_CMD != egg_info ]]; then $PIP_INSTALL healpy; fi

? If you have no time on this you can also add me (temporarily) as a committer on your python-reprojection repo and I can try and wrap up this PR.

lpsinger · 2015-04-29T19:34:10Z

FYI: possibly related to using healpy.rotator for coordinate conversion... healpy/healpy#247

…to healpix

coveralls · 2015-04-30T06:42:06Z

Coverage decreased (-5.55%) to 80.66% when pulling c4b8f14 on keflavich:healpix into 95a3d43 on astrofrog:master.

astrofrog · 2015-04-30T08:47:24Z

I'm going to go ahead and merge this now to take the burden off @keflavich but @lpsinger, I have some questions I'll raise in a separate issue.

Add HEALPIX reprojection

This was referenced Sep 27, 2014

healpix_to_image healpy/healpy#202

Closed

Add HEALPIX to / from normal image reprojection functions healpy/healpy#129

Closed

Remove healpix_to_image function (moved to reproject repo) gammapy/gammapy#205

Merged

keflavich added 4 commits November 18, 2014 11:06

added HEALPIX code (copied in part from gammapy)

73e362f

sign flip in lat - confusing!

f05423e

add non-interpolated option and new example

b52cd25

added convenience functions for reprojecting files & HDUs

05dec3e

keflavich force-pushed the healpix branch from ad39807 to 05dec3e Compare November 18, 2014 10:06

update to use newest astropy API and remove invalid import

3e7bf94

lpsinger reviewed Apr 24, 2015
View reviewed changes

lpsinger and others added 2 commits April 27, 2015 13:13

Working and tested implementation of healpix_to_image, image_to_healp…

8eff729

…ix, healpix_reproject_file Note that there will be some issues with round-trip projection (HEALPix->WCS->HEALPix) if some of the pixels lie precisely on the WCS map boundary.

Merge pull request #1 from lpsinger/healpix

f731259

Working and tested implementation of healpix_to_image, image_to_healpix, and healpix_reproject_file

keflavich added 3 commits April 29, 2015 22:32

add healpy install to .travis

81dc63e

Merge branch 'healpix' of github.com:keflavich/python-reprojection in…

cdb4170

…to healpix

python 3: zips are iterables

c4b8f14

astrofrog added a commit that referenced this pull request Apr 30, 2015

Merge pull request #27 from keflavich/healpix

799f2af

Add HEALPIX reprojection

astrofrog merged commit 799f2af into astropy:master Apr 30, 2015

Add HEALPIX reprojection #27

Add HEALPIX reprojection #27

Conversation

keflavich commented Jan 7, 2014

cdeil commented Jan 7, 2014

keflavich commented Jan 10, 2014

cdeil commented Jan 10, 2014

cdeil commented Jan 21, 2014

keflavich commented Jan 21, 2014

astrofrog commented Jul 16, 2014

keflavich commented Jul 16, 2014

astrofrog commented Jul 16, 2014

astrofrog commented Jul 22, 2014

cdeil commented Jul 22, 2014

DrPhilEvans commented Nov 18, 2014

astrofrog commented Nov 18, 2014

astrofrog commented Nov 18, 2014

DrPhilEvans commented Nov 18, 2014

astrofrog commented Nov 18, 2014

keflavich commented Nov 18, 2014

astrofrog commented Nov 18, 2014

keflavich commented Nov 18, 2014

DrPhilEvans commented Nov 18, 2014

keflavich commented Nov 18, 2014

cdeil commented Nov 18, 2014

cdeil commented Apr 13, 2015

astrofrog commented Apr 13, 2015

cdeil commented Apr 13, 2015

lpsinger commented Apr 24, 2015

lpsinger Apr 24, 2015

Choose a reason for hiding this comment

keflavich commented Apr 24, 2015

lpsinger commented Apr 24, 2015

keflavich commented Apr 24, 2015

lpsinger commented Apr 24, 2015

lpsinger commented Apr 24, 2015

astrofrog commented Apr 24, 2015

bsipocz commented Apr 24, 2015

astrofrog commented Apr 25, 2015

astrofrog commented Apr 25, 2015

DrPhilEvans commented Apr 27, 2015

lpsinger commented Apr 27, 2015

DrPhilEvans commented Apr 27, 2015

astrofrog commented Apr 27, 2015

lpsinger commented Apr 27, 2015

astrofrog commented Apr 29, 2015

lpsinger commented Apr 29, 2015

coveralls commented Apr 30, 2015

astrofrog commented Apr 30, 2015