Regular celestial images and cubes

One of the most common types of data to reproject are celestial images or n-dimensional data (such as spectral cubes) where two of the axes are celestial. There are several existing algorithms that can be used to reproject such data:

  • Interpolation (such as nearest-neighbor, bilinear, biquadratic interpolation and so on). This is the fastest algorithm and is suited to common use cases, but it is important to note that it is not flux conserving, and will not return optimal results if the input and output pixel sizes are very different.
  • Drizzling, which consists of determining the exact overlap fraction of pixels, and optionally allows pixels to be rescaled before reprojection. A description of the algorithm can be found in Fruchter and Hook (2002). This method is more accurate than interpolation but is only suitable for images where the field of view is small so that pixels are well approximated by rectangles in world coordinates. This is slower but more accurate than interpolation for small fields of view.
  • Computing the exact overlap of pixels on the sky by treating them as four-sided spherical polygons on the sky and computing spherical polygon intersection. This is essentially an exact form of drizzling, and should be appropriate for any field of view. However, this comes at a significant performance cost. This is the algorithm used by the Montage package, and we have implemented it here using the same core algorithm.

Currently, this package implements interpolation and spherical polygon intersection.


The reprojection/resampling is always done assuming that the image is in surface brightness units. For example, if you have an image with a constant value of 1, reprojecting the image to an image with twice as high resolution will result in an image where all pixels are all 1. This limitation is due to the interpolation algorithms (the fact that interpolation can be used implicitly assumes that the pixel values can be interpolated which is only the case with surface brightness units). If you have an image in flux units, first convert it to surface brightness, then use the functions described below. In future, we will provide a convenience function to return the area of all the pixels to make it easier.


Reprojection using interpolation can be done using the high-level reproject_interp() function:

>>> from reproject import reproject_interp

This function takes two main arguments. The first argument is the image to reproject, together with WCS information about the image. This can be either:

In the case of a FITS file or an HDUList object, if there is more than one Header-Data Unit (HDU), the hdu_in argument is also required and should be set to the ID or the name of the HDU to use.

The second argument is the WCS information for the output image, which should be specified either as a WCS or a Header instance. If this is specified as a WCS instance, the shape_out argument to reproject_interp() should also be specified, and be given the shape of the output image using the Numpy (ny, nx) convention (this is because WCS, unlike Header, does not contain information about image size).

As an example, we start off by opening a FITS file using Astropy:

>>> from import fits
>>> hdu ='')[0]    
Downloading [Done]

The image is currently using a Plate Carée projection:

>>> hdu.header['CTYPE1']   

We can create a new header using a Gnomonic projection:

>>> new_header = hdu.header.copy()   
>>> new_header['CTYPE1'] = 'GLON-TAN'   
>>> new_header['CTYPE2'] = 'GLAT-TAN'   

And finally we can call the reproject_interp() function to reproject the image:

>>> from reproject import reproject_interp
>>> new_image, footprint = reproject_interp(hdu, new_header)   

The reproject_interp() function returns two arrays - the first is the reprojected input image, and the second is a ‘footprint’ array which shows the fraction of overlap of the input image on the output image grid. This footprint is 0 for output pixels that fall outside the input image, 1 for output pixels that fall inside the input image. For more information about footprint arrays, see the Footprint arrays section.

We can then easily write out the reprojected image to a new FITS file:

>>> fits.writeto('reprojected_image.fits', new_image, new_header)   

The order of the interpolation can be controlled by setting the order= argument to either an integer or a string giving the order of the interpolation. Supported strings include:

  • 'nearest-neighbor': zeroth order interpolation
  • 'bilinear': fisst order interpolation
  • 'biquadratic': second order interpolation
  • 'bicubic': third order interpolation


Support for the drizzle algorithm will be implemented in future versions.

Spherical Polygon Intersection

Exact reprojection using the spherical polygon intersection can be done using the high-level reproject_exact() function:

>>> from reproject import reproject_exact

The two first arguments, the input data and the output projection, should be specified as for the reproject_interp() function described in Interpolation. In addition, an optional parallel= option can be used to control whether to parallelize the reprojection, and if so how many cores to use (see reproject_exact() for more details). For this algorithm, the footprint array returned gives the exact fractional overlap of new pixels with the original image (see Footprint arrays for more details).