The reprojected healpix images are displayed by DS9 with an incorrect right ascension value.
The problem seems to stem from the parameter CRVAL1 = 180.0 used to set the _targer_header_.
By setting CRVAL1 = 0, DS9 correctly displays the image.

This problem does not occur using Aladin Sky Atlas.

One common use case is for people to reproject an image so it's aligned with the celestial north. We should provide a function that given a header or WCS returns a header or WCS that is aligned with the north for that coordinate system. It would also be great to be able to specify which coordinate system use (e.g. ICRS north, Galactic north, etc.).

I was getting a crash executing the example from the docs here:

from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
hdu1 = fits.open(get_pkg_data_filename('galactic_center/gc_2mass_k.fits'))[0]
hdu2 = fits.open(get_pkg_data_filename('galactic_center/gc_msx_e.fits'))[0]
from reproject import reproject_interp, reproject_exact
array_interp, footprint_interp = reproject_interp(hdu2, hdu1.header)
array_exact, footprint_exact = reproject_exact(hdu2, hdu1.header)

The crash report is here: https://gist.github.com/cdeil/cda42d1c5357e053b1c2

This is with Python 3.4 from Macports and from an IPython notebook.
I don't know what the problem is.
The IPython console output mentions a CFITSIO version mismatch:

$ ipython notebook
[I 22:10:11.655 NotebookApp] Using MathJax from CDN: https://cdn.mathjax.org/mathjax/latest/MathJax.js
[W 22:10:11.699 NotebookApp] Terminals not available (error was No module named 'terminado')
[I 22:10:11.699 NotebookApp] Serving notebooks from local directory: /private/tmp
[I 22:10:11.699 NotebookApp] 0 active kernels 
[I 22:10:11.699 NotebookApp] The IPython Notebook is running at: http://localhost:8888/
[I 22:10:11.699 NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation).
[I 22:10:29.713 NotebookApp] Creating new notebook in 
[I 22:10:30.557 NotebookApp] Kernel started: 1e1521ba-8ed9-47a1-ad7c-dcc47599478e
[I 22:11:42.558 NotebookApp] KernelRestarter: restarting kernel (1/5)
WARNING:root:kernel 1e1521ba-8ed9-47a1-ad7c-dcc47599478e restarted
[I 22:12:31.523 NotebookApp] Saving file at /Untitled.ipynb
[I 22:14:31.521 NotebookApp] Saving file at /Untitled.ipynb
[I 22:16:31.518 NotebookApp] Saving file at /Untitled.ipynb
[I 22:18:32.513 NotebookApp] Saving file at /Untitled.ipynb

WARNING: version mismatch between CFITSIO header (v3.34) and linked library (v3.35).

[I 22:20:32.511 NotebookApp] Saving file at /Untitled.ipynb
[I 22:24:33.508 NotebookApp] Saving file at /Untitled.ipynb
[I 22:26:34.508 NotebookApp] Saving file at /Untitled.ipynb
[I 22:28:35.501 NotebookApp] Saving file at /Untitled.ipynb
Assertion failed: ok (mailbox.cpp:84)
Assertion failed: ok (mailbox.cpp:84)
[I 22:30:15.567 NotebookApp] KernelRestarter: restarting kernel (1/5)
WARNING:root:kernel 1e1521ba-8ed9-47a1-ad7c-dcc47599478e restarted

I've tried a few times, but I'm not able to reproduce the issue.
Just thought I'd mention it here in case someone else sees a similar crash ...

Continuing the discussion from https://github.com/astrofrog/reproject/pull/67#issuecomment-98984750 and #3, I think it would be nice if the high-level reproject API contained a flag conserve_flux (or something like that) that's False by default and when set does this ("sb" = "surface brightness):

solid_angle_image = compute_solid_angle_image(flux_image)
sb_image = flux_image / solid_angle_image
sb_image2 = reproject(sb_image)
solid_angle_image2 = compute_solid_angle_image(sb_image2)
flux_image2 = solid_angle_image2 * sb_image2

@astrofrog – Do you think this is a common enough usecase to expose in the high-level API functions or is there some better way?

By writing wrappers that take 2-d arrays (Nx4) we can make sure it is possible to compute areas/overlaps for many pixels/polygons.

Aside from #61 and #69, I'm not sure if anything else needs to be done before an initial release. Many of the other issues are not critical for 0.1 and I'd like to have this package see some real-world testing. I'm planning on having frequent releases.

The docs are now working and plots are showing: http://reproject.readthedocs.org/en/latest/

@cdeil - is there anything that you think should be done before a first release?

For the 0.1 release, we need to make sure that the API is at least close to what we'd like in the long term. This is also a reminder that I need to add support for HEALPIX via the high-level reproject interface, and we also need to make sure that we expose some of the lower level functions for more fine-tuned control.

I get this warning in the build:

/usr/bin/clang -fno-strict-aliasing -fno-common -dynamic -pipe -Os -fwrapv -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/include -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c reproject/overlapArea.c -o build/temp.macosx-10.9-x86_64-2.7/reproject/overlapArea.o
reproject/overlapArea.c:1262:7: warning: 'finite' is deprecated: first deprecated in OS X 10.9 [-Wdeprecated-declarations]
   if(mNaN(area) || area < 0.)
      ^
reproject/mNaN.h:7:30: note: expanded from macro 'mNaN'
#define mNaN(x) isnan(x) || !finite(x)
                             ^
/usr/include/math.h:718:12: note: 'finite' declared here
extern int finite(double) __OSX_AVAILABLE_BUT_DEPRECATED(__MAC_10_0, __MAC_10_9, __IPHONE_NA, __IPHONE_NA);
           ^
1 warning generated.

I'll have to look around what the proper and equivalent way of handling NaNs in this way is unless someone else knows.
Tests with input NaNs should be in place first.

Can you make this independent from hyperion?

$ python setup.py build
Traceback (most recent call last):
  File "setup.py", line 15, in <module>
    from hyperion.testing.helper import HyperionTest
ImportError: No module named hyperion.testing.helper

Or maybe put it in photutils or astropy?

With Python 2.6 I get this error:

$ python2.6 setup.py test
running test
running build
running build_py
running build_ext
Traceback (most recent call last):
  File "/Users/deil/code/python-reprojection/runtests.py", line 2841, in <module>
    do_exec(entry, locals())
  File "<string>", line 1, in do_exec
  File "<string>", line 1, in <module>
  File "py._apipkg", line 159, in __getattribute__
  File "py._apipkg", line 144, in getmod
  File "py._apipkg", line 37, in importobj
  File "/Users/deil/code/python-reprojection/runtests.py", line 2817, in load_module
    do_exec(co, module.__dict__)
  File "<string>", line 1, in do_exec
  File "pytest", line 14, in <module>
  File "/Users/deil/code/python-reprojection/runtests.py", line 2817, in load_module
    do_exec(co, module.__dict__)
  File "<string>", line 1, in do_exec
  File "_pytest.config", line 376, in <module>
  File "py._std", line 15, in __getattr__
AttributeError: py.std: could not import argparse

Should we bundle argparse or import the astropy version?

It would be of great utility a command-line function to reproject an HEALPix image
just specifying the required resolution of the pixels in degrees.

I'm proposing the function below that would fit in a next release, if it is deemed appropriate.

Bests,
Giuseppe

from __future__ import division
import math
from astropy.utils.data import get_pkg_data_filename
from astropy.io import fits
from reproject import reproject_from_healpix
from astropy.wcs import WCS
from astropy.coordinates import FK5

class header_reproject:
     def __init__ (self, CTYPE, CDELT, input_image, output_image):
          """
          Quick reprojection of a LVC Healpix sky map using the plate scale in degrees per pixel (CDELT).
Parameters
          ----------

          CTYPE : string
                projection --- e.g. "CAR": Cartesian projection; "MOL": Mollweide projection; "AIT": AITOFF projection

          CDELT : float
                the plate scale in degrees per pixel; for square pixels, CDELT1 = CDELT2

          input_image : string
                input LVC HEALPix sky map

          output_image : string
                clobber=True

          Default parameters :
               header :
                    NAXIS    =  2 
                    NAXIS1   = (360/CDELT)
                    NAXIS2   = (180/CDELT)
                    CRPIX1   = [(360/CDELT)+1]/2
                    CRVAL1   = 180.0
                    CDELT1   = -CDELT
                    CUNIT1   = 'deg'
                    CRPIX2   = [(180/CDELT)+1]/2
                    CRVAL2   = 0.0
                    CDELT2   = CDELT
                    CUNIT2   = 'deg '
                    COORDSYS = 'icrs'

               nested : True
                    NESTED pixel ordering is assume

          Returns : 
                 the reprojected image in fits format.

          Example :
            ***save the full code as LVC_reproject.py***
               from LVC_reproject import reproject
               reproject.grid_size(CTYPE='CAR',CDELT=1,input_image='allsky/ligo_simulated.fits.gz',
                                             output_image="reproject_1.fits")
          """
          self.CTYPE=CTYPE
          self.CDELT=CDELT
          self.input_image=input_image
          self.output_image=output_image

     def grid_size(self,CTYPE,CDELT,input_image,output_image):
          self.CTYPE=CTYPE
          reference_header['CTYPE1']="RA---"+CTYPE
          reference_header['CTYPE2']="DEC--"+CTYPE

          self.CDELT=CDELT
          reference_header['CDELT1']=-CDELT
          reference_header['CDELT2']=CDELT
          x_axis=math.ceil(360/CDELT)
          y_axis=math.ceil(180/CDELT)

          reference_header['NAXIS1']=x_axis
          reference_header['CRPIX1']=(int(x_axis)+1)/2

          reference_header['NAXIS2']=y_axis
          reference_header['CRPIX2']=(int(y_axis)+1)/2

          self.input_image=input_image
          self.output_image=output_image


          def write_image(self):
               filename_LVC = get_pkg_data_filename(input_image)
               hdu_LVC = fits.open(filename_LVC)[1]

               array, footprint = reproject_from_healpix(filename_LVC, reference_header)

               reference_wcs = WCS(reference_header)
               array, footprint = reproject_from_healpix(filename_LVC, reference_wcs,
                                          shape_out=(240,480))
               data = hdu_LVC.data['PROB']
               array, footprint = reproject_from_healpix((data, 'icrs'),
                                           reference_header, nested=True)
               data = hdu_LVC.data['PROB']
               array, footprint = reproject_from_healpix((data, 'icrs'),
                                           reference_header, nested=True)

               array, footprint = reproject_from_healpix((data, FK5(equinox='J2010')),
                                          reference_header, nested=True)

               fits.writeto(output_image,array,reference_header,clobber=True)

          write_image(self)

#loading header parameters to grid_size
file_header="""
NAXIS   = 2 
NAXIS1  = 360
NAXIS2  = 180
CTYPE1  = ''
CRPIX1  = 180
CRVAL1  = 180.0
CDELT1  = -1.0 
CUNIT1  = 'deg'
CTYPE2  = ''
CRPIX2  = 90
CRVAL2  = 0.0
CDELT2  =  1.0
CUNIT2  = 'deg '
COORDSYS= 'icrs'
"""
reference_header=fits.Header.fromstring(file_header, sep='\n')

reproject=header_reproject(CTYPE='',CDELT="",input_image='',output_image="")

It looks like currently the C code is passing global variables to functions, which I think will not be thread-safe.

For example, it should detect whether the WCS is >2-dimensional, whether there is a celestial component, and whether the celestial WCS is independent.

During the build on i386, the test_reproject_celestial_3d test failed:

_________________________ test_reproject_celestial_3d __________________________

    def test_reproject_celestial_3d():
        [...]
        out_full, foot_full = _reproject_full(array_in, wcs_in, wcs_out, (3, 160, 170))
        out_celestial, foot_celestial = _reproject_celestial(array_in, wcs_in, wcs_out, (3, 160, 170))
>       np.testing.assert_allclose(out_full, out_celestial)

reproject/interpolation/tests/test_core.py:391: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x_id = array([[[ True,  True,  True, ...,  True,  True,  True],
        [ True,  True...  True],
        [ True,  True,  True, ...,  True,  True,  True]]], dtype=bool)
y_id = array([[[False, False, False, ...,  True,  True,  True],
        [False, False...  True],
        [ True,  True,  True, ...,  True,  True,  True]]], dtype=bool)
hasval = 'nan'

    def chk_same_position(x_id, y_id, hasval='nan'):
        try:
            assert_array_equal(x_id, y_id)
        except AssertionError:
            msg = build_err_msg([x, y],
                                err_msg + '\nx and y %s location mismatch:'
                                % (hasval), verbose=verbose, header=header,
                                names=('x', 'y'), precision=precision)
>           raise AssertionError(msg)
E           AssertionError: 
E           Not equal to tolerance rtol=1e-07, atol=0
E           
E           x and y nan location mismatch:
E            x: array([[[ nan,  nan,  nan, ...,  nan,  nan,  nan],
E                   [ nan,  nan,  nan, ...,  nan,  nan,  nan],
E                   [ nan,  nan,  nan, ...,  nan,  nan,  nan],...
E            y: array([[[  1.,   1.,   1., ...,  nan,  nan,  nan],
E                   [  1.,   1.,   1., ...,  nan,  nan,  nan],
E                   [  1.,   1.,   1., ...,  nan,  nan,  nan],...

/usr/lib/python2.7/dist-packages/numpy/testing/utils.py:680: AssertionError

Full build log here.

It is unclear whether the package is licensed under "BSD-2" or "BSD-3": there are many files that claim to be under BSD-3, but then refer to LICENSE.rst which is actually BSD-2. This makes it difficult to upload the package to Debian; it was just rejected because of the incomplete license I did.

Could I ask to clean this up and to put all files either under BSD-2 or BSD-3?

Do we want to support 32-bit floats as input / output?

Could be useful to avoid making 64-bit copies where very large 32-bit input images are used.

How would this be implemented?
Worth the effort?

reprojection should accept multiple hdus and then figure out the least common multiple image and reproject all the hdus into that, either individually or to a single final mosaic.

Basically do what montage.mosaic does.

@cdeil - with the current master, the following returns zero:

import numpy as np
from reproject import compute_overlap

ilon = [[  3e-5, 0.00, 0.00, 3e-5 ]]
ilat = [[  8e-8, 8e-5, 1e-4, 1e-4 ]]
olon = [[  1e-4, 7e-6, 7e-6, 1e-4 ]]
olat = [[  0.00, 0.00, 1e-4, 1e-4 ]]
olon = np.roll(olon, 2)[::-1]
olat = np.roll(olat, 2)[::-1]

print compute_overlap(ilon, ilat, olon, olat)

if I replace the .c file with the original one before modifications (and keep the rest of the files on master), I get

(array([  2.83328916e-12]), array([ 1.]))

Any idea what is happening here and what change might have caused this?

Annoyingly, the following doesn't work, either with the old file, or the new file:

import numpy as np
from reproject import compute_overlap

ilon = [[  3e-5, 0.00, 0.00, 3e-5 ]]
ilat = [[  8e-8, 8e-5, 1e-4, 1e-4 ]]
olon = [[  1e-4, 7e-6, 7e-6, 1e-4 ]]
olat = [[  0.00, 0.00, 1e-4, 1e-4 ]]

print compute_overlap(ilon, ilat, olon, olat)

This returns zero in any case. Any ideas why this might be? I've been staring at it for two hours and I can't figure it out. Maybe a fresh pair of eyes will see the issue.

Two improvements need to be made to parse_output_projection:

• GIve an error if the shape is an empty tuple (can happen if NAXIS=0)
• Have an else clause to catch other cases where e.g. filename is passed.

WCSLIB includes a utility to go from a HEALPIX representation to a 2d HPX projection in an exact way with no interpolation (from which we can then use the existing celestial reprojection tools in reproject). Switching to using WCSLIB would therefore have the following advantages:

• It is LGPL instead of GPL so https://github.com/astrofrog/reproject/issues/82 is no longer an issue since we could drop the HEALPIX dependency.
• It is bundled in Astropy so we could maybe call the WCSLIB functions from there instead of having it as an extra dependency (maybe @mdboom could say if this is possible?)
• Since transforming to HPX is lossless, we can then use reproject_exact on the result and therefore implement #80 (cc @cdeil)

The downsides are:

• Not clear if there is a conversion tool for HPX -> HEALPIX but in principle there's no reason we couldn't try and contribute that to WCSLIB
• At the moment, this is a command-line utility in WCSLIB so we'd need that functionality to be moved to the library. I've asked Mark if this would be possible.

This is not something we can do immediately because of the two downsides, but I think we could make it work in the mid-term future!

At the moment the compute_overlap and solid_angle functions use degrees for the input angles and return the solid angle in steradians (see #3).

Wouldn't it be better to use radians as input for these functions.
This is what other trigonometric functions like sin do and it's usually faster to have all angles in radians in library code.

I can make the change if you agree.

Could you please add a license?

I didn't look closely yet, but I think I'd like to use this functionality for gammalib, a GPL3 package if possible:
https://cta-redmine.irap.omp.eu/issues/1017#change-3505

We should enable continuous integration on AppVeyor. If anyone wants to open a PR for this, I can flick the switch on the AppVeyor side! (we should just be able to re-use the yaml file from other astropy projects)

Hello,

I downloaded the .zip file of reproject,
and made installation using "python ./setup.py install".
It seems to be successful, but when I import reproject,
I ran into the error message as follows:

import reproject
Traceback (most recent call last):
File "", line 1, in
File "reproject/init.py", line 14, in
from .spherical_intersect import reproject_exact
File "reproject/spherical_intersect/init.py", line 5, in
from .high_level import *
File "reproject/spherical_intersect/high_level.py", line 6, in
from .core import _reproject_celestial
File "reproject/spherical_intersect/core.py", line 11, in
from ._overlap import _compute_overlap
ImportError: No module named _overlap

Is there any clue of what is missing?
Thanks very much.

Best wishes,
Baobab

This is a feature request to add a flux-conserving (and equivalently count-conserving) reprojection method to / from HEALPIX.

For HEALPIX <-> HEALPIX it would be trivial, just multiply / divide by 4 to the power of difference in HEALPIX order. E.g. if one gets a huge HEALPIX file but doesn't need the high resolution, one could downsample. Not sure if this should be in reproject or in healpy. @lpsinger Maybe this exists in healpy and I missed it?

For HEALPIX <-> normal WCS projection image I don't know how to do it. I thought HEALPIX pixel boundaries are great circles and the 4-corner polygon "exact" method from Montage could be used as well, but this page http://healpix.jpl.nasa.gov/html/intronode4.htm explains that pixel boundaries are not great circles on the sphere.
The use case I have for this count- and flux-preserving reprojection would be to be able to go back and forth between HEALPIX and normal projections and be able to do correct photometry of sources e.g. on FERMI data. @lpsinger Easy or hard to implement? Worth the effort?

Could also include Montage itself for comparison.

As proposed by @mdboom in astropy/astropy#2731 (comment) we should return the coverage mask showing the footprint of the old image in the pixels of the new projection.

Do you want to also wrap the functionality in overlapAreaPP.c?
I think we should either do this or remove it.

EDIT: changed title from "Wrap or remove overlapAreaPP.c?" to "Wrap overlapAreaPP.c".

This is low priority, but in the long term it would be nice to be able to enable debugging of the C routines from the Python code.

This would need to take a HDU object since WCS doesn't technically contain the dimensions of the image, only the projection information. I think astropy.wcs.WCS does contain the naxis info, but I've been told that this could be removed at any point since WCS doesn't say anything about the dimensions of the image.

Fun times ahead:

• Matplotlib is <1.4 on RTD but need >= 1.4 for plot directive to work properly (and pip won't update matplotlib on RTD)
• Healpy complains about setuptools not being recent enough: The required version of setuptools (>=14.3.1) is not available

Can't update setuptools ourselves because pip needs to run egg_info on all dependencies before installing anything.

I think I'm just going to switch to static plots in the docs. I don't have time to waste trying to find workarounds for the above.

This isn't really an issue so much as a place to discuss comparison between tools.

I have two reprojection routines in FITS_tools now, both under the label 'hcongrid'. The AST version is capable of transforming between projections (gal<->cel, etc.), while the scipy map_coordinates version requires that the input and output be in the same system. At the moment, the initial routine Tom has written up in #25 is loosely equivalent to hcongrid.

https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/hcongrid.py

Right now it appears that reprojecting by interpolation is only available as polynomial interpolation. Is that an intended restriction? For example, I would love to see Lanczos resampling (due to it's excellent smoothness and theoretically near-ideal S/N properties). But I presume other people have preferred algorithms they might also want to use?

Not much more to add, see: http://reproject.readthedocs.org/en/stable/

As discussed in https://github.com/astrofrog/reproject/issues/58 and at #pyastro15, the reproject functions themselves should stay as-is but we'll want a function that can take multiple HDUs and then figure out the output WCS to use, and later we can then create an easy 'mosaic' function that just wraps the function to find the header, and then calls reproject on all tiles before combining them.

@bsipocz - I'll assume that you are working on this. No rush of course, but just to say that this is your project if you still want to do it :)

Sorry to bother you again...

The license file which is in the git repository does not get packaged into the pypi tarball.

When reprojecting large images, we may run into memory issues for storing the corners or centers of all pixels, so we should have a way to project chunk by chunk (for all reprojection methods)

It would be nice to have a high-level docs page describing which reprojection functions preserve flux and which preserve surface brightness.

This is the #1 confusion for first-time Montage users and @ellisowen and I weren't sure what the current functions in reproject do at first. If I understand correctly the one in interpolation roughly preserves surface brightness (depending on the input image and interpolation parameters) and the one in spherical_intersect preserves flux up to rounding errors. I.e. for different-resolution input / output images the result from the two functions will be very different in normalisation?

This was originally meant to just be a python version of the montage algorithm, but since we want to use it as a general repository for reprojection, we need a better name. What about reproject-utils or reprojectutils?

The method-specific reprojection routines should be imported into the top level, but for that they need to be given different names, such as reproject_interpolation_2d, reproject_drizzle_2d, and reproject_whatever_2d where whatever is a better name than spherical_intersection

A few remaining things to do for the HEALPIX projection:

• @lpsinger - do we need to use sky_to_sky and the healpy rotator class, or can we switch to using a slightly modified version of the function here: https://github.com/astrofrog/reproject/blob/master/reproject/wcs_utils.py#L17 which uses astropy.coordinates?
• Why does the reference header define COORDSYS in the examples? As far as I know this is not part of the FITS standard, and instead we should use RADESYS if the header is in RA/Dec and otherwise the coordinate system is clear from the CTYPE1/CTYPE2. For normal (non-healpix) WCS we can use the existing Astropy function wcs_to_celestial_frame to return a coordinate frame from the WCS (and it will figure out all the equinox etc. stuff correctly)
• I think we want to avoid dealing with filenames in healpix.py and instead have any I/O in the high-level reproject function to keep things simple.
• We need to improve the test coverage, which took a dive down with these changes
• We need to document the HEALPIX reprojection in the main docs.
• The healpix code should be moved to a sub-directory called healpix to keep things tidy

I feel dirty writing this, but I have to re-release under GPL. Healpy is GPL v2 because of HEALPIX being GPL v2, and importing is enough to make reproject GPL. This is really annoying because I was hoping to use reproject in APLpy but now I can't. There are two ways forward:

• Ask the HEALPIX developers to re-release under LGPL, which I think solves this issue
• Write out own implementation of HEALPIX which may not be so hard because we actually only use a very small fraction of the functionality in HEALPIX

UPDATE: as discussed below, not relicensing under GPL since it's just for an optional dependency, so have just added a note that if functions using that dependency are used, it should be treated as GPL.

Given an HDU object and a Header object, project the original image onto the new Header and return a new HDU object.

Once #18 is merged, we can switch to using Cython.

@cdeil - do you want to take a shot at this?

As discussed in #3, we can write a routine directly at the C level to compute polygon areas (rather than relying on compute_overlap as in #3).

This code:

https://github.com/astrofrog/reproject/blob/master/reproject/spherical_intersect/reproject.py#L80

should probably be in Cython or C for speed - there are currently four nested Python loops!

Hi folk,
when I try to project an Healpix image in plate carrée projection,
the coordinate grid in DS9 are shifted - as in attachment image.
The same in Aladin Sky Atlas but to a limited extent.
This for different input values in the header;
I am attaching an example code:

from astropy.io import fits
hdu_ligo = fits.open(filename_ligo)[1]

target_header = fits.Header.fromstring("""
NAXIS = 2
NAXIS1 = 1000
NAXIS2 = 800
CTYPE1 = 'RA---CAR'
CRPIX1 = 500
CRVAL1 = 0
CDELT1 = -0.4
CUNIT1 = 'deg '
CTYPE2 = 'DEC--CAR'
CRPIX2 = 400
CRVAL2 = 0.0
CDELT2 = 0.4
CUNIT2 = 'deg '
COORDSYS= 'icrs '
""", sep='\n')

array, footprint = reproject_from_healpix(filename_ligo, target_header)

from astropy.wcs import WCS
target_wcs = WCS(

target_header)
array, footprint = reproject_from_healpix(filename_ligo, target_wcs,
shape_out=(240,480))

data = hdu_ligo.data['PROB']
array, footprint = reproject_from_healpix((data, 'icrs'),
target_header, nested=True)

from astropy.coordinates import FK5
array, footprint = reproject_from_healpix((data, FK5(equinox='J2010')),
target_header, nested=True)

fits.writeto('ds9_out.fits', array,target_header)

many thanks,
Giuseppe

We should add a small speed benchmark for a few common use cases.

This should probably be done before changing the code to work without global variables as proposed in #8.

cc @perrygreenfield