Examples of cropping an NDCube

One of the powerful aspects of having coordinate-aware data stored as an 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 are illustrated.

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

Let’s begin by creating an example NDCube object. For this case, we’ll generate an 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.

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

So we now have created an 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.

Here we can inspect the cube by plotting it

example_cube.plot()

<mpl_animators.wcs.ArrayAnimatorWCS object at 0x7fcdcdfaf010>

We can also inspect the dimensions of the cube:

example_cube.dimensions

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

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

example_cube.axis_world_coords()

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord
  warnings.warn(

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

Slicing and cropping the cube

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

Slicing by array index

To slice the example_cube so that we extract only one wavelength, we can do: by indexing as such

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

<WCSAxes: >

Cropping cube using world coordinate values using 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 ndcube.NDCube.crop() method which takes high-level astropy coordinate objects, such as SkyCoord. ndcube.NDCube.crop() returns the smallest cube in array-index space that contains all the passed points.

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 SkyCoord and 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:

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)

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord
  warnings.warn(

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

cropped_cube.dimensions

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

and we can visualize it:

cropped_cube.plot()

<mpl_animators.wcs.ArrayAnimatorWCS object at 0x7fcdbffd9f90>

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 crop().

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

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

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:671: AstropyUserWarning: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change
  warnings.warn(

we can inspect the dimensions of the cropped cube:

cropped_cube.dimensions

<Quantity [21., 11.] pix>

and again visualize it:

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)

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:671: AstropyUserWarning: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change
  warnings.warn(

Check dimensions:

cropped_cube.dimensions

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

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:

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

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/stable/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:671: AstropyUserWarning: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change
  warnings.warn(

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