Combining images into mosaics

Warning

The mosaicking functionality in the reproject package is currently experimental, so use with care and please report issues at https://github.com/astropy/reproject

The reproject.mosaicking sub-package includes helper functions for constructing mosaics from multiple images. These are find_optimal_celestial_wcs(), which can be used to construct a single optimal WCS/shape that overlaps with multiple images, and reproject_and_coadd(), which given images and a target WCS/shape will reproject all the images then combine them into a mosaic. We describe these in the sections below.

For the examples on this page we will use the PyVO module to retrieve tiles from the 2MASS survey around the M17 region:

.. doctest-requires:: pyvo
>>> from astropy.io import fits
>>> from astropy.coordinates import SkyCoord
>>> from pyvo.dal import imagesearch
>>> pos = SkyCoord.from_name('M17')
>>> table = imagesearch('https://irsa.ipac.caltech.edu/cgi-bin/2MASS/IM/nph-im_sia?type=at&ds=asky&',
...                    pos, size=0.25).to_table()
>>> table = table[(table['band'].astype('S') == 'K') & (table['format'].astype('S') == 'image/fits')]
>>> m17_hdus =  [fits.open(url)[0] for url in table['download'].astype('S')]

Computing an optimal WCS

Basic usage

Given a series of images, the find_optimal_celestial_wcs() function can be used to find an output WCS and shape (i.e. an output header) that overlaps with all the inpute images. Note that you don’t necessarily need to use this if you already know the final header or WCS you want to use for the images - in this case you can skip straight to Reprojecting and co-adding images into a mosaic.

You can optionally provide options to try and constrain the solution, as we will see below. To start off, let’s consider the simplest example, which is to call find_optimal_celestial_wcs() with the files downloaded above, but no additional information:

.. doctest-requires:: pyvo
>>> from reproject.mosaicking import find_optimal_celestial_wcs
>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus)

The first argument to find_optimal_celestial_wcs() should be a list where each element is either a filename, an HDU object (e.g. PrimaryHDU or ImageHDU), an HDUList object, or a tuple of (array, wcs). In the example above, we have passed a list of HDUs. We can now look at the output WCS and shape:

.. doctest-requires:: pyvo
>>> wcs_out.to_header()  
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =      900.07981909504 / Pixel coordinate of reference point
CRPIX2  =      1099.9484609446 / Pixel coordinate of reference point
CDELT1  =     -0.0002777777845 / [deg] Coordinate increment at reference point
CDELT2  =      0.0002777777845 / [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'           / Right ascension, gnomonic projection
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection
CRVAL1  =      275.18472258448 / [deg] Coordinate value at reference point
CRVAL2  =     -16.141349044263 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =     -16.141349044263 / [deg] Native latitude of celestial pole
...
>>> shape_out
(2201, 1800)

Coordinate system

By default, the coordinate system of the first file is used, and the final WCS is set up so that North (in that coordinate system) is up. In the case above, the images are in equatorial coordinates, so the final WCS is also in equatorial coordinates. We can force the output WCS to instead be in Galactic coordinates by setting the frame= argument to a coordinate frame object such as Galactic or one of the string shortcuts defined in astropy (e.g. 'fk5', 'galactic', etc.):

.. doctest-requires:: pyvo
>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus,
...                                                 frame='galactic')

the resulting WCS is then in Galactic coordinates:

.. doctest-requires:: pyvo
>>> wcs_out.to_header()  
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =      1214.1034417971 / Pixel coordinate of reference point
CRPIX2  =      1310.1351675461 / Pixel coordinate of reference point
CDELT1  =     -0.0002777777845 / [deg] Coordinate increment at reference point
CDELT2  =      0.0002777777845 / [deg] Coordinate increment at reference point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'GLON-TAN'           / galactic longitude, gnomonic projection
CTYPE2  = 'GLAT-TAN'           / galactic latitude, gnomonic projection
CRVAL1  =      15.116960053834 / [deg] Coordinate value at reference point
CRVAL2  =    -0.72166403860488 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =    -0.72166403860488 / [deg] Native latitude of celestial pole
...
>>> shape_out
(2623, 2424)

Orientation

As mentioned above, by default the image will be lined up so that North is up, but this may not always be optimal because if the mosaic is rotated compared to North, there may be a lot of empty space in the final mosaic. The auto_rotate option can therefore be used to find the optimal rotation for the WCS that minimizes the area of the final mosaic file:

>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus,
...                                                 frame='galactic',
...                                                 auto_rotate=True)

Note that this requires Shapely 1.6 or later to be installed. We can now look at the final WCS and shape:

>>> wcs_out.to_header()  
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =      1102.3949574309 / Pixel coordinate of reference point
CRPIX2  =      900.46211361965 / Pixel coordinate of reference point
PC1_1   =     0.88188439264557 / Coordinate transformation matrix element
PC1_2   =     0.47146571244169 / Coordinate transformation matrix element
PC2_1   =    -0.47146571244169 / Coordinate transformation matrix element
PC2_2   =     0.88188439264557 / Coordinate transformation matrix element
CDELT1  =     -0.0002777777845 / [deg] Coordinate increment at reference point
CDELT2  =      0.0002777777845 / [deg] Coordinate increment at reference point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'GLON-TAN'           / galactic longitude, gnomonic projection
CTYPE2  = 'GLAT-TAN'           / galactic latitude, gnomonic projection
CRVAL1  =      15.116960053834 / [deg] Coordinate value at reference point
CRVAL2  =    -0.72166403860488 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =    -0.72166403860488 / [deg] Native latitude of celestial pole
...
>>> shape_out
(1800, 2201)

As expected, the optimal shape is smaller than was returned previously.

Pixel resolution

By default, the final mosaic will have the pixel resolution (i.e. the pixel scale along the pixel axes) of the highest resolution input image, but this can be overriden using the resolution= keyword argument:

.. doctest-requires:: pyvo
>>> from astropy import units as u
>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus,
...                                                 resolution=1.5 * u.arcsec)

Projection and reference coordinate

Finally, you can customize the projection to use as well as the reference coordinate. To change the projection from the default (which is the gnomonic projection, or TAN), you can use the projection= keyword argument, which should be set to a valid three-letter FITS-WCS projection code:

.. doctest-requires:: pyvo
>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus,
...                                                 projection='CAR')

To customize the reference coordinate (where the projection is centered) you can set the reference= keyword argument to an astropy SkyCoord object:

.. doctest-requires:: pyvo
>>> from astropy.coordinates import SkyCoord
>>> coord = SkyCoord.from_name('M17')
>>> wcs_out, shape_out = find_optimal_celestial_wcs(m17_hdus,
...                                                 reference=coord)

Reprojecting and co-adding images into a mosaic

Assuming that you have a set of images that you want to combine into a mosaic, as well as a target header or WCS and shape (which you either determined independently, or with Computing an optimal WCS), you can make use of the reproject_and_coadd() function to produce the mosaic:

.. doctest-requires:: pyvo
>>> from reproject import reproject_interp
>>> from reproject.mosaicking import reproject_and_coadd
>>> array, footprint = reproject_and_coadd(m17_hdus,
...                                        wcs_out, shape_out=shape_out,
...                                        reproject_function=reproject_interp)

The first argument to reproject_and_coadd() should be a list where each element is either a filename, an HDU object (e.g. PrimaryHDU or ImageHDU), an HDUList object, or a tuple of (array, wcs).

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

Finally, the reproject_function should be used to specify which function to use to reproject individual tiles - this should be either reproject_interp() or reproject_exact() - with the latter being slower but more accurate (and necessary for flux conservation). Keyword arguments for these functions (e.g. order for reproject_interp()) can be passed as keyword arguments to reproject_and_coadd().

The example above will return an array which is the mosaic itself, and a footprint, which shows how many input images contributed to each output pixel. We can take a look at the output:

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 8))
ax1 = plt.subplot(1, 2, 1)
im1 = ax1.imshow(array, origin='lower', vmin=600, vmax=800)
ax1.set_title('Mosaic')
ax2 = plt.subplot(1, 2, 2)
im2 = ax2.imshow(footprint, origin='lower')
ax2.set_title('Footprint')

(png, svg, pdf)

_images/mosaicking-2.png

In some cases, including the above example, each tile that was used to compute the mosaic has an arbitrary offset due e.g. to different observing conditions. The reproject_and_coadd() includes an option to match the backgrounds (assuming a constant additive offset in each image):

.. doctest-requires:: pyvo
>>> array_bgmatch, _ = reproject_and_coadd(m17_hdus,
...                                        wcs_out, shape_out=shape_out,
...                                        reproject_function=reproject_interp,
...                                        match_background=True)

By adjusting the stretch, we can see the difference more clearly between the mosaic made with background matching and that made without - the one without shows vertical striping, especially on the left.

import numpy as np
import matplotlib.pyplot as plt

ax1 = plt.subplot(1, 2, 1)
im1 = ax1.imshow(array, origin='lower', vmin=635, vmax=660)
ax1.set_title('No background matching')
ax2 = plt.subplot(1, 2, 2)
im2 = ax2.imshow(array_bgmatch, origin='lower', vmin=635, vmax=660)
ax2.set_title('Background matching')

(png, svg, pdf)

_images/mosaicking-4.png

The background matching works by finding all overlapping pairs of images and determining the median difference for each pair, then using a stochastic gradient descent method to find the optimal additive corrections (a positive or negative constant for each image) to minimize differences. We additionally place the constraint that the average correction should be zero, but since there’s no reason that the average correction should be exactly zero, you should be aware that the final mosaic may be offset from the absolute surface brightness/flux by a constant additive factor. The algorithm should be robust for many situations and does not currently have any exposed options for fine tuning.