NDCube

NDCube is the fundamental class of the ndcube package and is designed to handle data contained in a single N-D array described by a single set of WCS transformations. NDCube is subclassed from astropy.nddata.NDData and so inherits the same attributes for data, wcs, uncertainty, mask, meta, and unit. The WCS object contained in the .wcs attribute is subclassed from astropy.wcs.WCS and contains a few additional attributes to enable to keep track of its relationship to the data.

Initialization

To initialize the most basic NDCube object, all you need is a numpy.ndarray containing the data, and an astropy.wcs.WCS object describing the transformation from array-element space to real world coordinates. Let’s create a 3-D array of data with shape (3, 4, 5) where every value is 1:

>>> import numpy as np
>>> data = np.ones((3, 4, 5))

Now let’s create an astropy.wcs.WCS object describing the translation from the array element coordinates to real world coordinates. Let the first data axis be helioprojective longitude, the second be helioprojective latitude, and the third be wavelength. Note that due to (confusing) convention, the order of the axes in the WCS object is reversed relative to the data array.

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

Now that we have a data array and a corresponding WCS object, we can create an NDCube instance by doing:

>>> from ndcube import NDCube
>>> my_cube = NDCube(data, input_wcs)

The data array is stored in the mycube.data attribute while the WCS object is stored in the my_cube.wcs attribute. However, when manipulating/slicing the data is it better to slice the object as a whole. (See section on Slicing.) So the .data attribute should only be used to access a specific value(s) in the data. Another thing to note is that as part of the initialization, the WCS object is converted from an astropy.wcs.WCS to an ndcube.utils.wcs.WCS object which has some additional features for tracking “missing axes”, etc. (See section on Missing Axes.)

Thanks to the fact that NDCube is subclassed from astropy.nddata.NDData, you can also supply additional data to the NDCube instance. These include: metadata (dict or dict-like) located at NDCube.meta; a data mask (boolean numpy.ndarray) located at NDCube.mask marking, for example, reliable and unreliable pixels; an uncertainty array (numpy.ndarray) located at NDCube.uncertainty describing the uncertainty of each data array value; and a unit (astropy.units.Unit or unit str). For example:

>>> mask = np.zeros_like(my_cube.data, dtype=bool)
>>> meta = {"Description": "This is example NDCube metadata."}
>>> my_cube = NDCube(data, input_wcs, uncertainty=np.sqrt(data),
...                         mask=mask, meta=meta, unit=None)
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]

N.B. The above warning is 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.

Dimensions

NDCube has useful properties for inspecting its data shape and axis types, dimensions and world_axis_physical_types:

>>> my_cube.dimensions
<Quantity [3., 4., 5.] pix>
>>> my_cube.world_axis_physical_types
('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')

dimensions returns an Quantity of pixel units giving the length of each dimension in the NDCube while world_axis_physical_types returns an iterable of strings denoting the type of physical property represented by each axis. The axis names are in accordance with the International Virtual Observatory Alliance (IVOA) UCD1+ controlled vocabulary. Here the shape and axis types are given in data order, not WCS order.

Slicing

Arguably NDCube’s most powerful capability is its slicing. Slicing an NDCube instance using the standard slicing notation allows users to access sub-regions of their data while simultaneously slicing not only the other array attributes (e.g. uncertainty, mask, etc.) but also the WCS object. This ensures that even though the data array has changed size and shape, each array element will still correspond to the same real world coordinates as they did before. An example of how to slice a 3-D NDCube object is:

>>> my_cube_roi = my_cube[3:5, 10:100, 30:37]

Slicing can also reduce the dimension of an NDCube, e.g.:

>>> my_2d_cube = my_cube[0, 10:100, 30:37]

In addition to slicing by index, NDCube supports a basic version of slicing/indexing by real world coordinates via the crop_by_coords method. This takes a list of astropy.units.Quantity instances representing the minimum real world coordinates of the region of interest in each dimension. The order of the coordinates must be the same as the order of the data axes. A second iterable of Quantity must also be provided which gives the widths of the region of interest in each data axis:

>>> from astropy.units import Quantity
>>> my_cube_roi = my_cube.crop_by_coords(
... [Quantity(0.7, unit="deg"), Quantity(1.3e-5, unit="deg"), Quantity(1.04e-9, unit="m")],
... [Quantity(0.6, unit="deg"), Quantity(1., unit="deg"), Quantity(0.08e-9, unit="m")])

This method does not rebin or interpolate the data if the region of interest does not perfectly map onto the array’s “pixel” grid. Instead it translates from real world to pixel coordinates and rounds to the nearest integer pixel before indexing/slicing the NDCube instance. Therefore it should be noted that slightly different inputs to this method can result in the same output.

Missing Axes

Some WCS axis types are coupled. For example, the helioprojective latitude and longitude of the Sun as viewed by a camera on a satellite orbiting Earth do not map independently to the pixel grid. Instead, the longitude changes as we move vertically along the same x-position if that single x-position is aligned anywhere other than perfectly north-south along the Sun’s central meridian. The analagous is true of the latitude for any y-pixel position not perfectly aligned with the Sun’s equator. Therefore, knowledge of both the latitude and longitude must be known to derive the pixel position along a single spatial axis and vice versa.

However, there are occasions when a data array may only contain one spatial axis, e.g. data from a slit-spectrograph. In this case, simply extracting the corresponding latitude or longitude axis from the WCS object would cause the translations to break.

To deal with this scenario, NDCube supports “missing” WCS axes. An additional attribute is added to the WCS object (NDCube.wcs.missing_axis) which is a list of bool type indicating which WCS axes do not have a corresponding data axis. This allows translation information on coupled axes to persist even if the data axes do not. This feature also makes it possible for NDCube to seamlessly reduce the data dimensionality via slicing. In the majority of cases a user will not need to worry about this feature. But it is useful to be aware of as many of the coordinate transformation functionalities of NDCube are only made possible by the missing axis feature.

Extra Coordinates

In the case of some datasets, there may be additional translations between the array elements and real world coordinates that are not included in the WCS. Consider a 3-D data cube from a rastering slit-spectrograph instrument. The first axis corresponds to the x-position of the slit as it steps across a region of interest in a given pattern. The second corresponds to latitude along the slit. And the third axis corresponds to wavelength. However, the first axis also corresponds to time, as it takes time for the slit to move and then take another exposure. It would be very useful to have the measurement times also associated with the x-axis. However, the WCS may only handle one translation per axis.

Fortunately, NDCube has a solution to this. Values at integer (pixel) steps along an axis can be stored within the object and accessed via the extra_coords property. To attach extra coordinates to an NDCube instance, provide an iterable of tuples of the form (str, int, Quantity or array-like) during instantiation. The 0th entry gives the name of the coordinate, the 1st entry gives the data axis to which the extra coordinate corresponds, and the 2nd entry gives the value of that coordinate at each pixel along the axis. So to add timestamps along the 0th axis of my_cube we do:

>>> from datetime import datetime, timedelta
>>> # Define our timestamps.  Must be same length as data axis.
>>> axis_length = int(my_cube.dimensions[0].value)
>>> timestamps = [datetime(2000, 1, 1)+timedelta(minutes=i)
...               for i in range(axis_length)]
>>> extra_coords_input = [("time", 0, timestamps)]
>>> # Generate NDCube as above, except now set extra_coords kwarg.
>>> my_cube = NDCube(data, input_wcs, uncertainty=np.sqrt(data),
...                  mask=mask, meta=meta, unit=None,
...                  extra_coords=extra_coords_input)
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]

The extra_coords property returns a dictionary where each key is a coordinate name entered by the user. The value of each key is itself another dictionary with keys 'axis' and 'value' giving the corresponding data axis number and coordinate value at each pixel as supplied by the user:

>>> my_cube.extra_coords 
{'time': {'axis': 0, 'value': [datetime.datetime(2000, 1, 1, 0, 0), datetime.datetime(2000, 1, 1, 0, 1), datetime.datetime(2000, 1, 1, 0, 2)]}}

Just like the data array and the WCS object, the extra coordinates are sliced automatically when the NDCube instance is sliced. So if we take the first slice of my_cube in the 0th axis, the extra time coordinate will only contain the value from that slice.:

>>> my_cube[0].extra_coords 
{'time': {'axis': None, 'value': datetime.datetime(2000, 1, 1, 0, 0)}}

Note that the axis value is now None because the dimensionality of the NDCube has been reduced via the slicing:

>>> my_cube[0].dimensions
<Quantity [4., 5.] pix>

and so the time extra coordinate no longer corresponds to a data axis. This would not have been the case if we had done the slicing so the length of the 0th axis was >1:

>>> my_cube[0:2].dimensions
<Quantity [2., 4., 5.] pix>
>>> my_cube[0:2].extra_coords 
{'time': {'value': [datetime.datetime(2000, 1, 1, 0, 0), datetime.datetime(2000, 1, 1, 0, 1)], 'axis': 0}}

Plotting

To quickly and easily visualize N-D data, NDCube provides a simple-to-use, yet powerful plotting method, plot, which produces a sensible visualization based on the dimensionality of the data. It is intended to be a useful quicklook tool and not a replacement for high quality plots or animations, e.g. for publications. The plot method can be called very simply, like so:

>>> my_cube.plot() 

The type of visualization returned depends on the dimensionality of the data within the NDCube object. For 1-D data a line plot is produced, similar to matplotlib.pyplot.plot. For 2-D data, an image is produced similar to that of matplotlib.pyplot.imshow. While for a >2-D data, a sunpy.visualization.imageanimator.ImageAnimatorWCS object is returned. This displays a 2-D image with sliders for each additional dimension which allow the user to animate through the different values of each dimension and see the effect in the 2-D image.

No args are required. The necessary information to generate the plot is derived from the data and metadata in the NDCube itself. Setting the x and y ranges of the plot can be done simply by indexing the NDCube object itself to the desired region of interest and then calling the plot method, e.g.:

>>> my_cube[0, 10:100, :].plot() 

In addition, some optional kwargs can be used to customize the plot. The axis_ranges kwarg can be used to set the axes ticklabels. See the ImageAnimatorWCS documentation for more detail. However, if this is not set, the axis ticklabels are automatically derived in real world coordinates from the WCS object within the NDCube.

By default the final two data dimensions are used for the plot axes in 2-D or greater visualizations, but this can be set by the user using the images_axes kwarg:

>>> my_cube.plot(image_axes=[0,1]) 

where the first entry in the list gives the index of the data index to go on the x-axis, and the second entry gives the index of the data axis to go on the y-axis.

In addition, the units of the axes or the data can be set by the unit_x_axis, unit_y_axis, unit kwargs. However, if not set, these are derived from the NDCube wcs and unit attributes.

Coordinate Transformations

The fundamental point the WCS system is the ability to easily translate between pixel and real world coordinates. For this purpose, NDCube provides convenience wrappers for the better known astropy functions, astropy.wcs.WCS.all_pix2world and astropy.wcs.WCS.all_world2pix. These are pixel_to_world, world_to_pixel, and axis_world_coords. It is highly recommended that when using NDCube these convenience wrappers are used rather than the original astropy functions for a few reasons. For example, they can track house-keeping data, are aware of “missing” WCS axis, are unit-aware, etc.

To use pixel_to_world, simply input Quantity objects with pixel units. Each Quantity corresponds to an axis so the number of Quantity objects should equal the number of data axes. Also, the order of the quantities should correspond to the data axes’ order, not the WCS order. The nth element of each Quantity describes the pixel coordinate in that axis. For example, if we wanted to transform the pixel coordinates of the pixel (2, 3, 4) in my_cube we would do:

>>> import astropy.units as u
>>> real_world_coords = my_cube.pixel_to_world(
... Quantity([2], unit=u.pix), Quantity([3], unit=u.pix), Quantity([4], unit=u.pix))

To convert two pixels with pixel coordinates (2, 3, 4) and (5, 6, 7), we would call pixel_to_world like so:

>>> real_world_coords = my_cube.pixel_to_world(
... Quantity([2, 5], unit=u.pix), Quantity([3, 6], unit=u.pix), Quantity([4, 7], unit=u.pix))

As can be seen, since each Quantity describes a different pixel coordinate of the same number of pixels, the lengths of each Quantity must be the same.

pixel_to_world returns a similar list of Quantities to those that were input, except that they are now in real world coordinates:

>>> real_world_coords
[<Quantity [1.40006967, 2.6002542 ] deg>, <Quantity [1.49986193, 2.99724799] deg>, <Quantity [1.10e-09, 1.16e-09] m>]

The exact units used are defined within the NDCube instance’s WCS object. Once again, the coordinates of the nth pixel is given by the nth element of each of the Quantity objects returned.

Using world_to_pixel to convert real world coordinates to pixel coordinates is exactly the same, but in reverse. This time the input Quantity objects must be in real world coordinates compatible with those defined in the NDCube instance’s WCS object. The output is a list of Quantity objects in pixel unit.:

>>> pixel_coords = my_cube.world_to_pixel(
... Quantity(1.40006967, unit="deg"), Quantity(1.49986193, unit="deg"),
...  Quantity(1.10000000e-09,  unit="m"))
>>> pixel_coords
[<Quantity 2.00000001 pix>, <Quantity 3. pix>, <Quantity 4. pix>]

Note that both pixel_to_pixel and world_to_pixel can handle non-integer pixels. Moreover, they can also handle pixel beyond the bounds of the NDCube and even negative pixels. This is because the WCS translations should be valid anywhere in space, and not just within the field of view of the NDCube. This capability has many useful applications, for example, in comparing observations from different instruments with overlapping fields of view.

There are times however, when you only want to know the real world coordinates of the NDCube field of view. To make this easy, NDCube has a another coordinate transformation method axis_world_coords. This method returns the real world coordinates for each pixel along a given data axis. So in the case of my_cube, if we wanted the wavelength axis we could call:

>>> my_cube.axis_world_coords(2)
<Quantity [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>

Note we set axes to 2 since axes is defined in data axis order. We can also define the axis using any unique substring from the axis names defined in ndcube.NDCube.world_axis_physical_types:

>>> my_cube.world_axis_physical_types
('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl')
>>> # Since 'wl' is unique to the wavelength axis name, let's use that.
>>> my_cube.axis_world_coords('wl')
<Quantity [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>

Notice how this returns the same result as when we set axes to the corresponding data axis number.

As discussed above, some WCS axes are not independent. For those axes, axis_world_coords returns a Quantity with the same number of dimensions as dependent axes. For example, helioprojective longitude and latitude are dependent. Therefore if we ask for longitude, we will get back a 2D Quantity with the same shape as the longitude x latitude axes lengths. For example:

>>> longitude = my_cube.axis_world_coords('lon')
>>> my_cube.dimensions
<Quantity [3., 4., 5.] pix>
>>> longitude.shape
(3, 4)
>>> longitude
<Quantity [[0.60002173, 0.59999127, 0.5999608 , 0.59993033],
           [1.        , 1.        , 1.        , 1.        ],
           [1.39997827, 1.40000873, 1.4000392 , 1.40006967]] deg>

It is also possible to request more than one axis’s world coordinates by setting axes to an iterable of data axis number and/or axis type strings.:

>>> my_cube.axis_world_coords(2, 'lon')
(<Quantity [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>,
 <Quantity [[0.60002173, 0.59999127, 0.5999608 , 0.59993033],
            [1.        , 1.        , 1.        , 1.        ],
            [1.39997827, 1.40000873, 1.4000392 , 1.40006967]] deg>)

Notice that the axes’ coordinates have been returned in the same order in which they were requested.

Finally, if the user wants the world coordinates for all the axes, axes can be set to None, which is in fact the default.:

>>> my_cube.axis_world_coords()
(<Quantity [[0.60002173, 0.59999127, 0.5999608 , 0.59993033],
          [1.        , 1.        , 1.        , 1.        ],
          [1.39997827, 1.40000873, 1.4000392 , 1.40006967]] deg>,
 <Quantity [[1.26915033e-05, 4.99987815e-01, 9.99962939e-01,
             1.49986193e+00],
          [1.26918126e-05, 5.00000000e-01, 9.99987308e-01,
           1.49989848e+00],
          [1.26915033e-05, 4.99987815e-01, 9.99962939e-01,
           1.49986193e+00]] deg>,
 <Quantity [1.02e-09, 1.04e-09, 1.06e-09, 1.08e-09, 1.10e-09] m>)

As stated previously, NDCube is only written to handle single arrays described by single WCS instances. For cases where data is made up of multiple arrays, each described by different WCS translations, ndcube has another class, NDCubeSequence, which will discuss in the next section.