# 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.

Note

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.

## Interpolation¶

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:

- The name of a FITS file
- An
`HDUList`

object - An image HDU object such as a
`PrimaryHDU`

,`ImageHDU`

, or`CompImageHDU`

instance - A tuple where the first element is a
`ndarray`

and the second element is either a`WCS`

or a`Header`

object

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 astropy.io import fits
>>> hdu = fits.open('http://data.astropy.org/galactic_center/gc_msx_e.fits')[0]
Downloading http://data.astropy.org/galactic_center/gc_msx_e.fits [Done]
```

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

```
>>> hdu.header['CTYPE1']
'GLON-CAR'
```

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

## Drizzling¶

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).