astro py:docs

- Index
- Modules

### Navigation

- next »
- « previous |
- Astropy v1.1.1 »
- N-dimensional datasets (`astropy.nddata`) »
- astropy.nddata.utils

# astropy.nddata.utils

## 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')

(Source code, png, hires.png, pdf)

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

(Source code, png, hires.png, 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')

(Source code, png, hires.png, pdf)

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.value, position.dec.value]
    >>> 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')

(Source code, png, hires.png, pdf)

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

(Source code, png, hires.png, pdf)

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