
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "generated/gallery/slicing_ndcube.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_generated_gallery_slicing_ndcube.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_generated_gallery_slicing_ndcube.py:


==============================
Examples of cropping an NDCube
==============================

One of the powerful aspects of having coordinate-aware data stored as an
`~ndcube.NDCube` is the ability to crop and slice the data and coordinates in a
standardised and easy way.

For example, there may be a region of interest you would like to crop out along a certain dimension
of your cube. In this example, this method to slice an `~ndcube.NDCube` are illustrated.

.. GENERATED FROM PYTHON SOURCE LINES 13-21

.. code-block:: default

    import astropy.wcs
    import numpy as np
    from astropy import units as u
    from astropy.coordinates import SkyCoord, SpectralCoord
    from sunpy.coordinates import frames

    from ndcube import NDCube








.. GENERATED FROM PYTHON SOURCE LINES 22-30

Let's begin by creating an example `~ndcube.NDCube` object.
For this case, we'll generate an `~ndcube.NDCube` that consists of 3 dimensions
space-space-wavelength. This is analogous to an observation including images in multiple
wavelengths.
Lets first define a 3-D numpy array, and then define the WCS information that describes the
data. We'll just create an array of random numbers, and a WCS which consists of the coordinate information
which in this case will be in Helioprojective (i.e. an observation of the Sun) in latitude and longitude,
and over several wavelengths in the range of 10.2 - 11 angstrom.

.. GENERATED FROM PYTHON SOURCE LINES 30-44

.. code-block:: default


    # Define the data of random numbers. Here the spatial dimensions are (45, 45) and the wavelength dimension is 5.
    data = np.random.rand(5, 45, 45)
    # Define the WCS
    wcs = astropy.wcs.WCS(naxis=3)
    wcs.wcs.ctype = 'HPLT-TAN', 'HPLN-TAN', "WAVE"
    wcs.wcs.cunit = 'arcsec', 'arcsec', 'Angstrom'
    wcs.wcs.cdelt = 10, 10, 0.2
    wcs.wcs.crpix = 2, 2, 0
    wcs.wcs.crval = 1, 1, 10
    wcs.wcs.cname = 'HPC lat', 'HPC lon', 'wavelength'
    # Instantiate the `~ndcube.NDCube`
    example_cube = NDCube(data, wcs=wcs)








.. GENERATED FROM PYTHON SOURCE LINES 45-49

So we now have created an `~ndcube.NDCube` named ``example_cube``.
You may have noticed that the order of the WCS is reversed to the array order - this
is normal convention, and something to remember throughout.
Now let's first inspect the cube.

.. GENERATED FROM PYTHON SOURCE LINES 51-52

Here we can inspect the cube by plotting it

.. GENERATED FROM PYTHON SOURCE LINES 52-54

.. code-block:: default

    example_cube.plot()




.. image-sg:: /generated/gallery/images/sphx_glr_slicing_ndcube_001.png
   :alt: slicing ndcube
   :srcset: /generated/gallery/images/sphx_glr_slicing_ndcube_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <mpl_animators.wcs.ArrayAnimatorWCS object at 0x7fe425512390>



.. GENERATED FROM PYTHON SOURCE LINES 55-56

We can also inspect the dimensions of the cube:

.. GENERATED FROM PYTHON SOURCE LINES 56-58

.. code-block:: default

    example_cube.dimensions





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <Quantity [ 5., 45., 45.] pix>



.. GENERATED FROM PYTHON SOURCE LINES 59-60

We can also inspect the world coordinates for all array elements:

.. GENERATED FROM PYTHON SOURCE LINES 60-62

.. code-block:: default

    example_cube.axis_world_coords()





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    (<SkyCoord (Helioprojective: obstime=None, rsun=695700.0 km, observer=None): (Tx, Ty) in arcsec
        [[( -8.99999999,  -8.99999998), ( -8.99999999,   1.        ),
          ( -8.99999999,  10.99999998), ..., ( -9.00000009, 410.99945954),
          ( -9.00000009, 420.99941904), ( -9.00000009, 430.99937657)],
         [(  1.        ,  -8.99999999), (  1.        ,   1.        ),
          (  1.        ,  10.99999999), ..., (  1.        , 410.99946002),
          (  1.        , 420.99941954), (  1.        , 430.99937708)],
         [( 10.99999999,  -8.99999998), ( 10.99999999,   1.        ),
          ( 10.99999999,  10.99999998), ..., ( 11.00000009, 410.99945954),
          ( 11.00000009, 420.99941904), ( 11.00000009, 430.99937657)],
         ...,
         [(410.99945993,  -8.99998221), (410.99946002,   0.99999802),
          (410.99946012,  10.99997826), ..., (410.99946397, 410.99864807),
          (410.99946407, 420.99858784), (410.99946417, 430.99852562)],
         [(420.99941944,  -8.99998133), (420.99941954,   0.99999793),
          (420.99941964,  10.99997719), ..., (420.99942359, 410.99860798),
          (420.99942369, 420.99854677), (420.99942379, 430.99848358)],
         [(430.99937698,  -8.99998044), (430.99937708,   0.99999783),
          (430.99937719,  10.99997609), ..., (430.99938123, 410.99856693),
          (430.99938133, 420.99850472), (430.99938143, 430.99844053)]]>, <SpectralCoord [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>)



.. GENERATED FROM PYTHON SOURCE LINES 63-69

Slicing and cropping the cube
-----------------------------
An `~ndcube.NDCube` can be sliced and cropped both by array indexing (similar to the way a numpy
array in indexed) or by real world coordinates. When we use array indices we say we
are "slicing" the cube. When we use world coordinates we say we are "cropping" the cube.
Let's begin by slicing by array index.

.. GENERATED FROM PYTHON SOURCE LINES 71-75

Slicing by array index
----------------------
To slice the ``example_cube`` so that we extract only one wavelength, we can do:
by indexing as such

.. GENERATED FROM PYTHON SOURCE LINES 75-90

.. code-block:: default


    sliced_cube = example_cube[1, :, :]
    # here we can see we are left with a 2-D cube which is an image at one wavelength.
    sliced_cube.dimensions

    # We can also index a region of interest of the cube at a particular wavelength.
    # Again note that we are slicing here based on the ``array`` index rather than cropping by
    # real world value

    sliced_cube = example_cube[1, 10:20, 20:40]
    sliced_cube.dimensions

    # Now we can inspect the sliced cube, and see it's now a smaller region of interest.
    sliced_cube.plot()




.. image-sg:: /generated/gallery/images/sphx_glr_slicing_ndcube_002.png
   :alt: slicing ndcube
   :srcset: /generated/gallery/images/sphx_glr_slicing_ndcube_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <WCSAxesSubplot: >



.. GENERATED FROM PYTHON SOURCE LINES 91-98

Cropping cube using world coordinate values using :meth:`ndcube.NDCube.crop`
----------------------------------------------------------------------------
In many cases it's more useful to crop a cube to a region of interest based
on real world coordinates such as points in space or over some spectral
range. This is achieved by the :meth:`ndcube.NDCube.crop` method which takes high-level astropy coordinate objects,
such as `~astropy.coordinates.SkyCoord`. :meth:`ndcube.NDCube.crop` returns the smallest cube
in array-index space that contains all the passed points.

.. GENERATED FROM PYTHON SOURCE LINES 100-108

Let's first define some points over which to crop the ``example_cube``. The points are
defined as iterables of scalar high-level coordinate objects. We must provide the same
number of objects in each tuple as required by the WCS to describe all the world axes.
In this example, this means we need to provide a `~astropy.coordinates.SkyCoord` and
`~astropy.coordinates.SpectralCoord` in each iterable. However, if we don't want to crop
by one of the coordinate types, e.g. wavelength, we can replace the corresponding
high-level coordinate object in each iterable with ``None``.
Let's first define two points in space (lat and long) we want to crop but keep all wavelengths:

.. GENERATED FROM PYTHON SOURCE LINES 108-114

.. code-block:: default


    point1 = [SkyCoord(0*u.arcsec, 0*u.arcsec, frame=frames.Helioprojective), None]
    point2 = [SkyCoord(200*u.arcsec, 100*u.arcsec, frame=frames.Helioprojective), None]

    cropped_cube = example_cube.crop(point1, point2)








.. GENERATED FROM PYTHON SOURCE LINES 115-116

Similar to before, we can inspect the dimensions of the sliced cube:

.. GENERATED FROM PYTHON SOURCE LINES 116-119

.. code-block:: default


    cropped_cube.dimensions





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <Quantity [ 5., 21., 11.] pix>



.. GENERATED FROM PYTHON SOURCE LINES 120-121

and we can visualize it:

.. GENERATED FROM PYTHON SOURCE LINES 121-124

.. code-block:: default


    cropped_cube.plot()




.. image-sg:: /generated/gallery/images/sphx_glr_slicing_ndcube_003.png
   :alt: slicing ndcube
   :srcset: /generated/gallery/images/sphx_glr_slicing_ndcube_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <mpl_animators.wcs.ArrayAnimatorWCS object at 0x7fe42458b320>



.. GENERATED FROM PYTHON SOURCE LINES 125-127

Now let's say we also want to crop out the image that includes a specific wavelength.
Let's define a new point, and then include it with the first two we passed to :meth:`~ndcube.NDCube.crop`.

.. GENERATED FROM PYTHON SOURCE LINES 127-132

.. code-block:: default


    point3 = [None, SpectralCoord(10.2*u.angstrom)]

    cropped_cube = example_cube.crop(point1, point2, point3)








.. GENERATED FROM PYTHON SOURCE LINES 133-134

we can inspect the dimensions of the cropped cube:

.. GENERATED FROM PYTHON SOURCE LINES 134-137

.. code-block:: default


    cropped_cube.dimensions





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <Quantity [21., 11.] pix>



.. GENERATED FROM PYTHON SOURCE LINES 138-139

and again visualize it:

.. GENERATED FROM PYTHON SOURCE LINES 139-149

.. code-block:: default


    cropped_cube.plot()

    # Now let's say we instead want to crop over a wavelength range.
    # Let's define a new point, and then include it with the first two we passed to
    # :meth:`~ndcube.NDCube.crop`.
    point4 = [None, SpectralCoord(10.6*u.angstrom)]

    cropped_cube = example_cube.crop(point1, point2, point3, point4)




.. image-sg:: /generated/gallery/images/sphx_glr_slicing_ndcube_004.png
   :alt: slicing ndcube
   :srcset: /generated/gallery/images/sphx_glr_slicing_ndcube_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 150-151

Check dimensions:

.. GENERATED FROM PYTHON SOURCE LINES 151-154

.. code-block:: default


    cropped_cube.dimensions





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <Quantity [ 3., 21., 11.] pix>



.. GENERATED FROM PYTHON SOURCE LINES 155-158

Here we can just see how powerful this can be to easily crop over different world coordinates.
In fact we can make this simpler still by combining the SkyCoords and SpectralCoords into
just two points:

.. GENERATED FROM PYTHON SOURCE LINES 158-164

.. code-block:: default


    point5 = [SkyCoord(0*u.arcsec, 0*u.arcsec, frame=frames.Helioprojective), SpectralCoord(10.2*u.angstrom)]
    point6 = [SkyCoord(200*u.arcsec, 100*u.arcsec, frame=frames.Helioprojective), SpectralCoord(10.6*u.angstrom)]

    cropped_cube = example_cube.crop(point5, point6)
    cropped_cube.dimensions




.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <Quantity [ 3., 21., 11.] pix>




.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.634 seconds)


.. _sphx_glr_download_generated_gallery_slicing_ndcube.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: slicing_ndcube.py <slicing_ndcube.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: slicing_ndcube.ipynb <slicing_ndcube.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
