# 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)]
>>> timestamps1 = [timestamps0[-1]+timedelta(minutes=i+1) for i in range(data1.shape)]
>>> timestamps2 = [timestamps1[-1]+timedelta(minutes=i+1) for i in range(data2.shape)]
>>> 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[:, :, :] = False
>>> cube_meta = {"Description": "This is example NDCube metadata."}
>>> my_cube0 = NDCube(data0, input_wcs, uncertainty=np.sqrt(data0),
...                          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),
...                          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),
...                          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],
```

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],
```

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¶

Just like `NDCube`, `NDCubeSequence` provide simple but powerful plotting APIs to help users visualize their data. Two plotting methods, `plot` and `plot_as_cube`, are provided which correspond to the sequence and cube-like representations of the data, respectively. These methods allows the sequence to be animated as though it were one contiguous `NDCube`. Both methods have the same API and same kwargs as `ndcube.NDCube.plot`. See documentation for `ndcube.NDCube.plot` for more details. The main substantive difference between them is how the axis inputs relate to dimensionality of the data, i.e. the same way that the inputs to NDCubeSequence slicing and `index_as_cube` differ.

## 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.data for cube in my_sequence.data])
>>> print(data.shape)
(3, 4, 5)
```