Skip to main content
`Ctrl`+`K`

# Image Utilities

## Overview

The `astropy.nddata.utils` module includes general utility functions for array operations.

## 2D Cutout Images

### Getting Started

The `Cutout2D` class can be used to create a postage stamp cutout image from a 2D array. If an optional `WCS` object is input to `Cutout2D`, then the `Cutout2D` object will contain an updated `WCS` corresponding to the cutout array.

First, we simulate a single source on a 2D data array. If you would like to simulate many sources, see Efficient Model Rendering with Bounding Boxes.

Note: The pair convention is different for size and position! The position is specified as (x,y), but the size is specified as (y,x).

    >>> import numpy as np
    >>> from astropy.modeling.models import Gaussian2D
    >>> y, x = np.mgrid[0:500, 0:500]
    >>> data = Gaussian2D(1, 50, 100, 10, 5, theta=0.5)(x, y)

Now, we can display the image:

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots()
    >>> ax.imshow(data, origin='lower')

(`png`, `svg`, `pdf`)

Next we can create a cutout for the single object in this image. We create a cutout centered at position `(x, y) = (49.7, 100.1)` with a size of `(ny, nx) = (41, 51)` pixels:

    >>> from astropy.nddata import Cutout2D
    >>> from astropy import units as u
    >>> position = (49.7, 100.1)
    >>> size = (41, 51)     # pixels
    >>> cutout = Cutout2D(data, position, size)

The `size` keyword can also be a `Quantity` object:

    >>> size = u.Quantity((41, 51), u.pixel)
    >>> cutout = Cutout2D(data, position, size)

Or contain `Quantity` objects:

    >>> size = (41*u.pixel, 51*u.pixel)
    >>> cutout = Cutout2D(data, position, size)

A square cutout image can be generated by passing an integer or a scalar `Quantity`:

    >>> size = 41
    >>> cutout2 = Cutout2D(data, position, size)

    >>> size = 41 * u.pixel
    >>> cutout2 = Cutout2D(data, position, size)

The cutout array is stored in the `data` attribute of the `Cutout2D` instance. If the `copy` keyword is `False` (default), then `cutout.data` will be a view into the original `data` array. If `copy=True`, then `cutout.data` will hold a copy of the original `data`. Now we display the cutout image:

    >>> cutout = Cutout2D(data, position, (41, 51))
    >>> fig, ax = plt.subplots()
    >>> ax.imshow(cutout.data, origin='lower')

(`png`, `svg`, `pdf`)

The cutout object can plot its bounding box on the original data using the `plot_on_original()` method:

    >>> plt.imshow(data, origin='lower')
    >>> cutout.plot_on_original(color='white')

(`png`, `svg`, `pdf`)

Many properties of the cutout array are also stored as attributes, including:

    >>> # shape of the cutout array
    >>> print(cutout.shape)
    (41, 51)

    >>> # rounded pixel index of the input position
    >>> print(cutout.position_original)
    (50, 100)

    >>> # corresponding position in the cutout array
    >>> print(cutout.position_cutout)
    (25, 20)

    >>> # (non-rounded) input position in both the original and cutout arrays
    >>> print((cutout.input_position_original, cutout.input_position_cutout))
    ((49.7, 100.1), (24.700000000000003, 20.099999999999994))

    >>> # the origin pixel in both arrays
    >>> print((cutout.origin_original, cutout.origin_cutout))
    ((25, 80), (0, 0))

    >>> # tuple of slice objects for the original array
    >>> print(cutout.slices_original)
    (slice(80, 121, None), slice(25, 76, None))

    >>> # tuple of slice objects for the cutout array
    >>> print(cutout.slices_cutout)
    (slice(0, 41, None), slice(0, 51, None))

There are also two `Cutout2D` methods to convert pixel positions between the original and cutout arrays:

    >>> print(cutout.to_original_position((2, 1)))
    (27, 81)

    >>> print(cutout.to_cutout_position((27, 81)))
    (2, 1)

### 2D Cutout Modes

There are three modes for creating cutout arrays: `'trim'`, `'partial'`, and `'strict'`. For the `'partial'` and `'trim'` modes, a partial overlap of the cutout array and the input `data` array is sufficient. For the `'strict'` mode, the cutout array has to be fully contained within the `data` array, otherwise an `PartialOverlapError` is raised. In all modes, non-overlapping arrays will raise a `NoOverlapError`. In `'partial'` mode, positions in the cutout array that do not overlap with the `data` array will be filled with `fill_value`. In `'trim'` mode only the overlapping elements are returned, thus the resulting cutout array may be smaller than the requested `size`.

The default uses `mode='trim'`, which can result in cutout arrays that are smaller than the requested `size`:

    >>> data2 = np.arange(20.).reshape(5, 4)
    >>> cutout1 = Cutout2D(data2, (0, 0), (3, 3), mode='trim')
    >>> print(cutout1.data)
    [[0. 1.]
     [4. 5.]]
    >>> print(cutout1.shape)
    (2, 2)
    >>> print((cutout1.position_original, cutout1.position_cutout))
    ((0, 0), (0, 0))

With `mode='partial'`, the cutout will never be trimmed. Instead it will be filled with `fill_value` (the default is `numpy.nan`) if the cutout is not fully contained in the data array:

    >>> cutout2 = Cutout2D(data2, (0, 0), (3, 3), mode='partial')
    >>> print(cutout2.data)
    [[nan nan nan]
     [nan  0.  1.]
     [nan  4.  5.]]

Note that for the `'partial'` mode, the positions (and several other attributes) are calculated for on the valid (non-filled) cutout values:

    >>> print((cutout2.position_original, cutout2.position_cutout))
    ((0, 0), (1, 1))
    >>> print((cutout2.origin_original, cutout2.origin_cutout))
    ((0, 0), (1, 1))
    >>> print(cutout2.slices_original)
    (slice(0, 2, None), slice(0, 2, None))
    >>> print(cutout2.slices_cutout)
    (slice(1, 3, None), slice(1, 3, None))

Using `mode='strict'` will raise an exception if the cutout is not fully contained in the data array:

    >>> cutout3 = Cutout2D(data2, (0, 0), (3, 3), mode='strict')
    PartialOverlapError: Arrays overlap only partially.

### 2D Cutout from a `SkyCoord` Position

The input `position` can also be specified as a `SkyCoord`, in which case a `WCS` object must be input via the `wcs` keyword.

First, we define a `SkyCoord` position and a `WCS` object for our data (usually this would come from your FITS header):

    >>> from astropy.coordinates import SkyCoord
    >>> from astropy.wcs import WCS
    >>> position = SkyCoord('13h11m29.96s -01d19m18.7s', frame='icrs')
    >>> wcs = WCS(naxis=2)
    >>> rho = np.pi / 3.
    >>> scale = 0.05 / 3600.
    >>> wcs.wcs.cd = [[scale*np.cos(rho), -scale*np.sin(rho)],
    ...               [scale*np.sin(rho), scale*np.cos(rho)]]
    >>> wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
    >>> wcs.wcs.crval = [position.ra.to_value(u.deg),
    ...                  position.dec.to_value(u.deg)]
    >>> wcs.wcs.crpix = [50, 100]

Now we can create the cutout array using the `SkyCoord` position and `wcs` object:

    >>> cutout = Cutout2D(data, position, (30, 40), wcs=wcs)
    >>> fig, ax = plt.subplots()
    >>> ax.imshow(cutout.data, origin='lower')

(`png`, `svg`, `pdf`)

The `wcs` attribute of the `Cutout2D` object now contains the propagated `WCS` for the cutout array. Now we can find the sky coordinates for a given pixel in the cutout array. Note that we need to use the `cutout.wcs` object for the cutout positions:

    >>> from astropy.wcs.utils import pixel_to_skycoord
    >>> x_cutout, y_cutout = (5, 10)
    >>> pixel_to_skycoord(x_cutout, y_cutout, cutout.wcs)
    <SkyCoord (ICRS): (ra, dec) in deg
        ( 197.8747893, -1.32207626)>

We now find the corresponding pixel in the original `data` array and its sky coordinates:

    >>> x_data, y_data = cutout.to_original_position((x_cutout, y_cutout))
    >>> pixel_to_skycoord(x_data, y_data, wcs)
    <SkyCoord (ICRS): (ra, dec) in deg
        ( 197.8747893, -1.32207626)>

As expected, the sky coordinates in the original `data` and the cutout array agree.

### 2D Cutout Using an Angular `size`

The input `size` can also be specified as a `Quantity` in angular units (e.g., degrees, arcminutes, arcseconds, etc.). For this case, a `WCS` object must be input via the `wcs` keyword.

For this example, we will use the data, `SkyCoord` position, and `wcs` object from above to create a cutout with size 1.5 x 2.5 arcseconds:

    >>> size = u.Quantity((1.5, 2.5), u.arcsec)
    >>> cutout = Cutout2D(data, position, size, wcs=wcs)
    >>> fig, ax = plt.subplots()
    >>> ax.imshow(cutout.data, origin='lower')

(`png`, `svg`, `pdf`)

## Saving a 2D Cutout to a FITS File with an Updated WCS

A `Cutout2D` object can be saved to a FITS file, including the updated WCS object for the cutout region. In this example, we download an example FITS image and create a cutout image. The resulting `Cutout2D` object is then saved to a new FITS file with the updated WCS for the cutout region.

    # Download an example FITS file, create a 2D cutout, and save it to a
    # new FITS file, including the updated cutout WCS.
    from astropy.io import fits
    from astropy.nddata import Cutout2D
    from astropy.utils.data import download_file
    from astropy.wcs import WCS

    def download_image_save_cutout(url, position, size):
        # Download the image
        filename = download_file(url)

        # Load the image and the WCS
        hdu = fits.open(filename)[0]
        wcs = WCS(hdu.header)

        # Make the cutout, including the WCS
        cutout = Cutout2D(hdu.data, position=position, size=size, wcs=wcs)

        # Put the cutout image in the FITS HDU
        hdu.data = cutout.data

        # Update the FITS header with the cutout WCS
        hdu.header.update(cutout.wcs.to_header())

        # Write the cutout to a new FITS file
        cutout_filename = "example_cutout.fits"
        hdu.writeto(cutout_filename, overwrite=True)

    if __name__ == "__main__":
        url = "https://astropy.stsci.edu/data/photometry/spitzer_example_image.fits"

        position = (500, 300)
        size = (400, 400)
        download_image_save_cutout(url, position, size)