NDCubeSequence

NDCubeSequence is a class for handling multiple NDCube objects as though they were one contiguous data set. Another way of thinking about it is that NDCubeSequence provides the ability to manipulate a data set described by multiple separate WCS transformations.

Regarding implementation, an NDCubeSequence instance is effectively a list of NDCube instances with some helper methods attached.

Initialization

To initialize the most basic NDCubeSequence object, all you need is a list of NDCube instances. So let us first define three 3-D NDCubes for slit-spectrograph data as we did in the NDCube section of this tutorial. First we define the data arrays and WCS objects:

>>> # Define data for cubes
>>> import numpy as np
>>> data0 = np.ones((3, 4, 5))
>>> data1 = data0 * 2
>>> data2 = data1 * 2

>>> # Define WCS object for all cubes.
>>> import astropy.wcs
>>> wcs_input_dict = {
... 'CTYPE1': 'WAVE    ', 'CUNIT1': 'Angstrom', 'CDELT1': 0.2, 'CRPIX1': 0, 'CRVAL1': 10, 'NAXIS1': 5,
... 'CTYPE2': 'HPLT-TAN', 'CUNIT2': 'deg', 'CDELT2': 0.5, 'CRPIX2': 2, 'CRVAL2': 0.5, 'NAXIS2': 4,
... 'CTYPE3': 'HPLN-TAN', 'CUNIT3': 'deg', 'CDELT3': 0.4, 'CRPIX3': 2, 'CRVAL3': 1, 'NAXIS3': 3}
>>> input_wcs = astropy.wcs.WCS(wcs_input_dict)

Let’s also define an extra coordinate of time assigned to the 0th cube data axis and another label coordinate assigned to the cubes as wholes. (See NDCube section of this guide of more detail.) Let the slices along the 0th axis be separated by one minute and the slices in preceding cube are followed directly in time by the slices in the next:

>>> from datetime import datetime, timedelta
>>> timestamps0 = [datetime(2000, 1, 1)+timedelta(minutes=i) for i in range(data0.shape[0])]
>>> timestamps1 = [timestamps0[-1]+timedelta(minutes=i+1) for i in range(data1.shape[0])]
>>> timestamps2 = [timestamps1[-1]+timedelta(minutes=i+1) for i in range(data2.shape[0])]
>>> extra_coords_input0 = [("time", 0, timestamps0), ("label", None, "hello")]
>>> extra_coords_input1 = [("time", 0, timestamps1), ("label", None, "world")]
>>> extra_coords_input2 = [("time", 0, timestamps2), ("label", None, "!")]

Now we can define our cubes.

>>> from ndcube import NDCube
>>> from ndcube import NDCubeSequence
>>> # Define a mask such that all array elements are unmasked.
>>> mask = np.empty(data0.shape, dtype=object)
>>> mask[:, :, :] = False
>>> cube_meta = {"Description": "This is example NDCube metadata."}
>>> my_cube0 = NDCube(data0, input_wcs, uncertainty=np.sqrt(data0),
...                          mask=mask, meta=cube_meta, unit=None,
...                          extra_coords=extra_coords_input0)
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
>>> my_cube1 = NDCube(data1, input_wcs, uncertainty=np.sqrt(data1),
...                          mask=mask, meta=cube_meta, unit=None,
...                          extra_coords=extra_coords_input1)
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
>>> my_cube2 = NDCube(data2, input_wcs, uncertainty=np.sqrt(data2),
...                          mask=mask, meta=cube_meta, unit=None,
...                          extra_coords=extra_coords_input2)
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]

N.B. The above warnings are due to the fact that astropy.nddata.uncertainty is recommended to have an uncertainty_type attribute giving a string describing the type of uncertainty. However, this is not required. Also note that due to laziness, we have used the same WCS translations in each NDCube instance above. However, it would be more common for each NDCube instance to have a different WCS, and in that case the usefulness of NDCubeSequence is more pronounced. Nonetheless, this case can still be used to adequately demonstrate the capabilities of NDCubeSequence.

Finally, creating an NDCubeSequence becomes is simple:

>>> my_sequence = NDCubeSequence([my_cube0, my_cube1, my_cube2])

While, each NDCube in the NDCubeSequence can have its own meta, it is also possible to supply additional metadata upon initialization of the NDCubeSequence. This metadata may be common to all sub-cubes or is specific to the sequence rather than the sub-cubes. This metadata is input as a dictionary:

>>> my_sequence_metadata = {"Description": "This is some sample NDCubeSequence metadata."}
>>> my_sequence = NDCubeSequence([my_cube0, my_cube1, my_cube2],
...                              meta=my_sequence_metadata)

and stored in the my_sequence.meta attribute. Meanwhile, the NDCube instances are stored in my_sequence.data. However, analgously to NDCube, it is strongly advised that the data is manipulated by slicing the NDCubeSequence rather than more manually delving into the .data attribute. For more explanation, see the section on Slicing.

Common Axis

It is possible (although not required) to set a common axis of the NDCubeSequence. A common axis is defined as the axis of the sub-cubes parallel to the axis of the sequence.

For example, assume the 0th axis of the sub-cubes, my_cube0, my_cube1 and my_cube2 in the NDCubeSequence, my_sequence, represent time as we have indicated by setting the time extra coordinate. In this case, my_cube0 represents observations taken from a period directly before my_cube1 and my_cube2 and the sub-cubes are ordered chronologically in the sequence. Then moving along the 0th axis of one sub-cube and moving along the sequence axis from one cube to the next both represent movement in time. The difference is simply the size of the steps. Therefore it can be said that the 0th axis of the sub-cubes is common to the sequence.

To define a common axis, set the kwarg during intialization of the NDCubeSequence to the desired data axis number:

>>> my_sequence = NDCubeSequence([my_cube0, my_cube1, my_cube2],
...                              meta=my_sequence_metadata, common_axis=0)

Defining a common axis enables the full range of the NDCubeSequence features to be utilized including ndcube.NDCubeSequence.plot, ndcube.NDCubeSequence.common_axis_extra_coords, and ndcube.NDCubeSequence.index_as_cube. See following sections for more details on these features.

Dimensions

Analagous to ndcube.NDCube.dimensions, there is also a ndcube.NDCubeSequence.dimensions property for easily inspecting the shape of an NDCubeSequence instance:

>>> my_sequence.dimensions
(<Quantity 3. pix>, <Quantity 3. pix>, <Quantity 4. pix>, <Quantity 5. pix>)

Slightly differently to ndcube.NDCube.dimensions, ndcube.NDCubeSequence.dimensions returns a tuple of astropy.units.Quantity instances with pixel units, giving the length of each axis. This is in constrast to the single Quantity returned by NDCube. This is because NDCubeSequence supports sub-cubes of different lengths along the common axis if it is set. In that case, the corresponding quantity in the dimensions tuple will have a length greater than 1 and list the length of each sub-cube along the common axis.

Equivalent to ndcube.NDCube.world_axis_physical_types, ndcube.NDCubeSequence.world_axis_physical_types returns a tuple of the physical axis types. The same IVOA UCD1+ controlled words are used for the cube axes as is used in ndcube.NDCube.world_axis_physical_types. The sequence axis is given the label 'meta.obs.sequence' as it is the IVOA UCD1+ controlled word that best describes it. To call, simply do:

>>> my_sequence.world_axis_physical_types
('meta.obs.sequence', 'custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

Slicing

As with NDCube, slicing an NDCubeSequence using the standard slicing API simulataneously slices the data arrays, WCS objects, masks, uncertainty arrays, etc. in each relevant sub-cube. For example, say we have three NDCubes in an NDCubeSequence, each of shape (3, 4, 5). Say we want to obtain a region of interest between the 1st and 2nd pixels (inclusive) in the 2nd dimension and 1st and 3rd pixels (inclusive) in the 3rd dimension of the 0th slice along the 0th axis in only the 1st (not 0th) and 2nd sub-cubes in the sequence. This would be a cumbersome slicing operation if treating the sub-cubes independently. (This would be made even worse without the power of NDCube where the data arrays, WCS objects, masks, uncertainty arrays, etc. would all have to be sliced independently!) However, with NDCubeSequence this becomes as simple as indexing a single array:

>>> regions_of_interest_in_sequence = my_sequence[1:3, 0, 1:3, 1:4]
>>> regions_of_interest_in_sequence.dimensions
(<Quantity 2. pix>, <Quantity 2. pix>, <Quantity 3. pix>)
>>> regions_of_interest_in_sequence.world_axis_physical_types
('meta.obs.sequence', 'custom:pos.helioprojective.lat', 'em.wl')

This will return a new NDCubeSequence with 2 2-D NDCubes, one for each region of interest from the 3rd slice along the 0th axis in each original sub-cube. If our regions of interest only came from a single sub-cube - say the 0th and 1st slices along the 0th axis in the 1st sub-cube - an NDCube is returned:

>>> roi_from_single_subcube = my_sequence[1, 0:2, 1:3, 1:4]
>>> roi_from_single_subcube.dimensions
<Quantity [2., 2., 3.] pix>
>>> roi_from_single_subcube.world_axis_physical_types
('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

If a common axis has been defined for the NDCubeSequence one can think of it as a contiguous data set with different sections along the common axis described by different WCS translations. Therefore it would be useful to be able to index the sequence as though it were one single cube. This can be achieved with the ndcube.NDCubeSequence.index_as_cube property. In our above example, my_sequence has a shape of (<Quantity 3. pix>, <Quantity 3.0 pix>, <Quantity 4.0 pix>, <Quantity 5.0 pix>) and a common axis of 0. Therefore we can think of my_sequence as a having an effective cube-like shape of (<Quantity 9.0 pix>, <Quantity 4.0 pix>, <Quantity 5.0 pix>) where the first sub-cube extends along the 0th cube-like axis from 0 to 3, the second from 3 to 6 and the third from 6 to 9. Say we want to extract the same region of interest as above, i.e. my_sequence[1, 0:2, 1:3, 1:4]. Then this can be acheived by entering:

>>> roi_from_single_subcube = my_sequence.index_as_cube[3:5, 1:3, 1:4]
>>> roi_from_single_subcube.dimensions
<Quantity [2., 2., 3.] pix>
>>> roi_from_single_subcube.world_axis_physical_types
('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

In this case the entire region came from a single sub-cube. However, index_as_cube also works when the region of interest spans multiple sub-cubes in the sequence. Say we want the same region of interest in the 2nd and 3rd cube dimensions from the final slice along the 0th cube axis of the 0th sub-cube, the whole 1st sub-cube and the 0th slice of the 2nd sub-cube. In cube-like indexing this corresponds to slices 2 to 7 along to the 0th cube axis:

>>> roi_across_subcubes = my_sequence.index_as_cube[2:7, 1:3, 1:4]
>>> roi_across_subcubes.dimensions
(<Quantity 3. pix>, <Quantity [1., 3., 1.] pix>, <Quantity 2. pix>, <Quantity 3. pix>)
>>> roi_across_subcubes.world_axis_physical_types
('meta.obs.sequence', 'custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

Notice that since the sub-cubes are now of different lengths along the common axis, the corresponding Quantity gives the lengths of each cube individually. See section on Dimensions for more detail.

Cube-like Dimensions

To help with handling an NDCubeSequence with a common axis as if it were a single cube, there exist cube-like equivalents of the dimensions and world_axis_physical_types methods. They are intuitively named cube_like_dimensions and cube_like_world_axis_physical_types. These give the lengths and physical types of the axes as if the data were stored in a single NDCube. So in the case of my_sequence, with three sub-cubes, each with a length of 3 along the common axis, we get:

>>> my_sequence.cube_like_dimensions
<Quantity [9., 4., 5.] pix>
>>> my_sequence.cube_like_world_axis_physical_types
('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

Note that cube_like_dimensions returns a single Quantity in pixel units, as if it were ndcube.NDCube.dimensions. This is in contrast to ndcube.NDCubeSequence.dimensions that returns a tuple of Quantity.

Common Axis Extra Coordinates

If a common axis is defined, it may be useful to view the extra coordinates along that common axis defined by each of the sub-cube extra_coords as if the NDCubeSequence were one contiguous Cube. This can be done using the common_axis_extra_coords property:

>>> my_sequence.common_axis_extra_coords
{'time': array([datetime.datetime(2000, 1, 1, 0, 0),
      datetime.datetime(2000, 1, 1, 0, 1),
      datetime.datetime(2000, 1, 1, 0, 2),
      datetime.datetime(2000, 1, 1, 0, 3),
      datetime.datetime(2000, 1, 1, 0, 4),
      datetime.datetime(2000, 1, 1, 0, 5),
      datetime.datetime(2000, 1, 1, 0, 6),
      datetime.datetime(2000, 1, 1, 0, 7),
      datetime.datetime(2000, 1, 1, 0, 8)], dtype=object)}

This returns a dictionary where each key gives the name of a coordinate. The value of each key is the values of that coordinate at each pixel along the common axis. Since all these coordinates must be along the common axis, it is not necessary to supply axis information as it is with ndcube.NDCube.extra_coords making ndcube.NDCubeSequence.common_axis_extra_coords simpler. Because this property has a functional form and calculates the dictionary each time from the constituent sub-cubes’ ndcube.NDCube.extra_coords attributes, ndcube.NDCubeSequence.common_axis_extra_coords is effectively sliced when the NDCubeSequence is sliced, e.g.:

>>> my_sequence[1:3].common_axis_extra_coords
{'time': array([datetime.datetime(2000, 1, 1, 0, 3),
      datetime.datetime(2000, 1, 1, 0, 4),
      datetime.datetime(2000, 1, 1, 0, 5),
      datetime.datetime(2000, 1, 1, 0, 6),
      datetime.datetime(2000, 1, 1, 0, 7),
      datetime.datetime(2000, 1, 1, 0, 8)], dtype=object)}

Sequence Axis Extra Coordinates

Analgous to common_axis_extra_coords, it is also possible to access the extra coordinates that are not assigned to any NDCube data axis via the ndcube.NDCubeSequence.sequence_axis_extra_coords property. Whereas common_axis_extra_coords returns all the extra coords with an 'axis' value equal to the common axis, sequence_axis_extra_coords returns all extra coords with an 'axis' value of None. Another way of thinking about this when there is no common axis set, is that they are assigned to the sequence axis. Hence the property’s name.:

>>> my_sequence.sequence_axis_extra_coords
{'label': array(['hello', 'world', '!'], dtype=object)}

Plotting

The ndcube.NDCubeSequence.plot method allows the sequence to be animated as though it were one contiguous NDCube. It has the same API and same kwargs as ndcube.NDCube.plot. See documentation for ndcube.NDCube.plot for more details.

Explode Along Axis

During analysis of some data - say of a stack of images - it may be necessary to make some different fine-pointing adjustments to each image that isn’t accounted for the in the original WCS translations, e.g. due to satellite wobble. If these changes are not describable with a single WCS object, it may be desirable to break up the N-D sub-cubes of an NDCubeSequence into an sequence of sub-cubes with dimension N-1. This would enable a separate WCS object to be associated with each image and hence allow individual pointing adjustments.

Rather than manually dividing the datacubes up and deriving the corresponding WCS object for each exposure, NDCubeSequence provides a useful method, explode_along_axis. To call it, simply provide the number of the data cube axis along which you wish to break up the sub-cubes:

>>> exploded_sequence = my_sequence.explode_along_axis(0)

Assuming we are using the same my_sequence as above, with dimensions.shape (<Quantity 3.0 pix>, <Quantity 3.0 pix>, <Quantity 4.0 pix>, <Quantity 5.0 pix>), the exploded_sequence will be an NDCubeSequence of nine 2-D NDCubes each with shape (<Quantity 4.0 pix>, <Quantity 5.0 pix>).:

>>> # Check old and new shapes of the squence
>>> my_sequence.dimensions
(<Quantity 3. pix>, <Quantity 3. pix>, <Quantity 4. pix>, <Quantity 5. pix>)
>>> exploded_sequence.dimensions
(<Quantity 9. pix>, <Quantity 4. pix>, <Quantity 5. pix>)

Note that any cube axis can be input. A common axis need not be defined.

Extracting Data Arrays

It is possible that you may have some procedures that are designed to operate on arrays instead of NDCubeSequence objects. “Therefore it may be useful to extract the data (or other array-like information such as uncertainty or mask) in the NDCubeSequence into a single ndarray. A succinct way of doing this operation is using python’s list comprehension features.

In the above examples we defined the my_sequence NDCubeSequence object.:

>>> # Print dimensions of my_sequence as a reminder
>>> print(my_sequence.dimensions)
(<Quantity 3. pix>, <Quantity 3. pix>, <Quantity 4. pix>, <Quantity 5. pix>)

In this section we will use this object to demonstrate extracting data arrays from NDCubeSequence objects. For example, say we wanted to make a 4D array out of the data arrays within the NDCubes of my_sequence.:

>>> # Make a single 4D array of data in sequence.
>>> data = np.stack([cube.data for cube in my_sequence.data])
>>> print(data.shape)
(3, 3, 4, 5)

If instead, we want to define a 3D array where every NDCube in the NDCubeSequence is appended together, we can use numpy’s vstack function:

>>> # Make a 3D array
>>> data = np.vstack([cube.data for cube in my_sequence.data])
>>> print(data.shape)
(9, 4, 5)

Finally, we can also create 3D arrays by slicing NDCubeSequence objects. Here we slice the NDCubeSequence along the fastest-changing dimension:

>>> # Slice sequence to make 3D array
>>> data = np.stack([cube[2].data for cube in my_sequence.data])
>>> print(data.shape)
(3, 4, 5)