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, let’s simulate a single source on a 2D data array. If you would like to simulate many sources, see Efficient Model Rendering with Bounding Boxes.

>>> 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, let’s display the image:

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

(png, svg, pdf)

../_images/utils-1.png

Now let’s 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) = (40, 50) pixels:

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

The size keyword can also be a Quantity object:

>>> size = u.Quantity((40, 50), u.pixel)
>>> cutout = Cutout2D(data, position, size)

or contain Quantity objects:

>>> size = (40*u.pixel, 50*u.pixel)
>>> cutout = Cutout2D(data, position, size)

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

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

>>> size = 40 * 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. Let’s display the cutout image:

>>> cutout = Cutout2D(data, position, (40, 50))
>>> plt.imshow(cutout.data, origin='lower')

(png, svg, pdf)

../_images/utils-2.png

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)

../_images/utils-3.png

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

>>> # shape of the cutout array
>>> print(cutout.shape)
(40, 50)

>>> # 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, 120, None), slice(25, 75, None))

>>> # tuple of slice objects for the cutout array
>>> print(cutout.slices_cutout)
(slice(0, 40, None), slice(0, 50, 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, let’s 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 let’s create the cutout array using the SkyCoord position and wcs object:

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

(png, svg, pdf)

../_images/utils-4.png

The wcs attribute of the Cutout2D object now contains the propagated WCS for the cutout array. Let’s 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’ll 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)
>>> plt.imshow(cutout.data, origin='lower')   

(png, svg, pdf)

../_images/utils-5.png

Reference/API

astropy.nddata.utils Module

This module includes helper functions for array operations.

Functions

extract_array(array_large, shape, position) Extract a smaller array of the given shape and position from a larger array.
add_array(array_large, array_small, position) Add a smaller array at a given position in a larger array.
subpixel_indices(position, subsampling) Convert decimal points to indices, given a subsampling factor.
overlap_slices(large_array_shape, …[, mode]) Get slices for the overlapping part of a small and a large array.
block_reduce(data, block_size[, func]) Downsample a data array by applying a function to local blocks.
block_replicate(data, block_size[, conserve_sum]) Upsample a data array by block replication.

Classes

NoOverlapError Raised when determining the overlap of non-overlapping arrays.
PartialOverlapError Raised when arrays only partially overlap.
Cutout2D(data, position, size[, wcs, mode, …]) Create a cutout object from a 2D array.