.. _nddata_utils:

Image Utilities
***************

Overview
========

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

.. _cutout_images:

2D Cutout Images
================

Getting Started
---------------

The `~astropy.nddata.utils.Cutout2D` class can be used to create a
postage stamp cutout image from a 2D array. If an optional
`~astropy.wcs.WCS` object is input to
`~astropy.nddata.utils.Cutout2D`, then the
`~astropy.nddata.utils.Cutout2D` object will contain an updated
`~astropy.wcs.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 :ref:`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:

.. doctest-skip::

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

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    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)
    plt.imshow(data, origin='lower')

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 `~astropy.units.Quantity` object::

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

Or contain `~astropy.units.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 `~astropy.units.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
`~astropy.nddata.utils.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:

.. doctest-skip::

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

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Gaussian2D
    from astropy.nddata import Cutout2D
    y, x = np.mgrid[0:500, 0:500]
    data = Gaussian2D(1, 50, 100, 10, 5, theta=0.5)(x, y)
    position = (49.7, 100.1)
    cutout = Cutout2D(data, position, (41, 51))
    plt.imshow(cutout.data, origin='lower')

The cutout object can plot its bounding box on the original data using
the :meth:`~astropy.nddata.utils.Cutout2D.plot_on_original` method:

.. doctest-skip::

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

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Gaussian2D
    from astropy.nddata import Cutout2D
    y, x = np.mgrid[0:500, 0:500]
    data = Gaussian2D(1, 50, 100, 10, 5, theta=0.5)(x, y)
    position = (49.7, 100.1)
    size = (41, 51)
    cutout = Cutout2D(data, position, size)
    plt.imshow(data, origin='lower')
    cutout.plot_on_original(color='white')

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))  # doctest: +FLOAT_CMP
    ((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 `~astropy.nddata.utils.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
`~astropy.nddata.utils.PartialOverlapError` is raised. In all modes,
non-overlapping arrays will raise a
`~astropy.nddata.utils.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)  # doctest: +FLOAT_CMP
    [[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)  # doctest: +FLOAT_CMP
    [[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:

.. doctest-skip::

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


2D Cutout from a `~astropy.coordinates.SkyCoord` Position
---------------------------------------------------------

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

First, we define a `~astropy.coordinates.SkyCoord` position and a
`~astropy.wcs.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
`~astropy.coordinates.SkyCoord` position and ``wcs`` object::

    >>> cutout = Cutout2D(data, position, (30, 40), wcs=wcs)
    >>> plt.imshow(cutout.data, origin='lower')   # doctest: +SKIP

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Gaussian2D
    from astropy.nddata import Cutout2D
    from astropy.coordinates import SkyCoord
    from astropy.wcs import WCS
    y, x = np.mgrid[0:500, 0:500]
    data = Gaussian2D(1, 50, 100, 10, 5, theta=0.5)(x, y)
    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]
    cutout = Cutout2D(data, position, (30, 40), wcs=wcs)
    plt.imshow(cutout.data, origin='lower')

The ``wcs`` attribute of the `~astropy.nddata.utils.Cutout2D` object now
contains the propagated `~astropy.wcs.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)    # doctest: +FLOAT_CMP
    <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)    # doctest: +FLOAT_CMP
    <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
`~astropy.units.Quantity` in angular units (e.g., degrees, arcminutes,
arcseconds, etc.). For this case, a `~astropy.wcs.WCS` object must be
input via the ``wcs`` keyword.

For this example, we will use the data, `~astropy.coordinates.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')   # doctest: +SKIP

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Gaussian2D
    from astropy.nddata import Cutout2D
    from astropy.coordinates import SkyCoord
    from astropy.wcs import WCS
    from astropy import units as u
    y, x = np.mgrid[0:500, 0:500]
    data = Gaussian2D(1, 50, 100, 10, 5, theta=0.5)(x, y)
    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]
    size = u.Quantity((1.5, 2.5), u.arcsec)
    cutout = Cutout2D(data, position, size, wcs=wcs)
    plt.imshow(cutout.data, origin='lower')


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

A `~astropy.nddata.utils.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
`~astropy.nddata.utils.Cutout2D` object is then saved to a new FITS file with
the updated WCS for the cutout region.

.. literalinclude:: examples/cutout2d_tofits.py
   :language: python