Combining a celestial WCS with a wavelength axis

The goal of this example is to construct a spectral-image cube of AIA images at different wavelengths.

This will showcase how to add an arbitrarily spaced wavelength dimension to a celestial WCS.

import matplotlib.pyplot as plt

import astropy.units as u

import sunpy.data.sample
import sunpy.map

from ndcube import NDCube
from ndcube.extra_coords import QuantityTableCoordinate
from ndcube.wcs.wrappers import CompoundLowLevelWCS

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/latest/lib/python3.12/site-packages/sunpy/util/sysinfo.py:111: SunpyUserWarning: Importing sunpy.map without its extra dependencies may result in errors.
The following packages are not installed:
['reproject>=0.10.0']
To install sunpy with these dependencies use `pip install sunpy[map]` or `pip install sunpy[all]` for all extras.
If you installed sunpy via conda, please report this to the community channel: https://matrix.to/#/#sunpy:openastronomy.org
  warn_user(f"Importing sunpy.{extras} without its extra dependencies may result in errors.\n"

We will use the sample data that sunpy provides to construct a sequence of AIA image files for different wavelengths using sunpy.map.Map.

aia_files = [sunpy.data.sample.AIA_094_IMAGE,
             sunpy.data.sample.AIA_131_IMAGE,
             sunpy.data.sample.AIA_171_IMAGE,
             sunpy.data.sample.AIA_193_IMAGE,
             sunpy.data.sample.AIA_211_IMAGE,
             sunpy.data.sample.AIA_304_IMAGE,
             sunpy.data.sample.AIA_335_IMAGE,
             sunpy.data.sample.AIA_1600_IMAGE]
# `sequence=True` causes a sequence of maps to be returned, one for each image file.
sequence_of_maps = sunpy.map.Map(aia_files, sequence=True)
# Sort the maps in the sequence in order of wavelength.
sequence_of_maps.maps = sorted(sequence_of_maps.maps, key=lambda m: m.wavelength)

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063305_0094_lowres.fits:   0%|          | 0.00/971k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.82file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.82file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063301_0131_lowres.fits:   0%|          | 0.00/962k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.58file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.58file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063302_0171_lowres.fits:   0%|          | 0.00/973k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  5.11file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  5.10file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063307_0193_lowres.fits:   0%|          | 0.00/1.00M [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  5.89file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  5.87file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063302_0211_lowres.fits:   0%|          | 0.00/988k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.96file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.96file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063334_0304_lowres.fits:   0%|          | 0.00/965k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.58file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.58file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063303_0335_lowres.fits:   0%|          | 0.00/979k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.46file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.46file/s]

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

AIA20110607_063305_1600_lowres.fits:   0%|          | 0.00/930k [00:00<?, ?B/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.26file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.26file/s]
/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/latest/lib/python3.12/site-packages/astropy/io/fits/hdu/image.py:634: VerifyWarning: Invalid 'BLANK' keyword in header.  The 'BLANK' keyword is only applicable to integer data, and will be ignored in this HDU.
  warnings.warn(msg, VerifyWarning)

Using an astropy.units.Quantity of the wavelengths of the images, we can construct a 1D lookup-table WCS via QuantityTableCoordinate. This is then combined with the celestial WCS into a single 3D WCS via CompoundLowLevelWCS.

wavelengths = u.Quantity([m.wavelength for m in sequence_of_maps])
wavelengths_wcs = QuantityTableCoordinate(wavelengths, physical_types="em.wl", names="wavelength").wcs
cube_wcs = CompoundLowLevelWCS(wavelengths_wcs, sequence_of_maps[0].wcs)

Combine the new 3D WCS with the stack of AIA images using ndcube.NDCube. Note that because we set the wavelength to the first axis in the WCS (cube_wcs), the final data cube is stacked such that wavelength corresponds to the array axis which is last. This is due to the convention that WCS axis ordering is reversed compared to data array axis ordering.

my_cube = NDCube(sequence_of_maps.as_array(), wcs=cube_wcs)
# Produce an interactive plot of the spectral-image stack.
my_cube.plot(plot_axes=['y', 'x', None])

plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/ndcube/checkouts/latest/examples/creating_even_spaced_wavelength_visualisation.py:56: SunpyDeprecationWarning: The as_array function is deprecated and will be removed in sunpy 7.1. Use the data and mask properties instead.
  my_cube = NDCube(sequence_of_maps.as_array(), wcs=cube_wcs)
/home/docs/checkouts/readthedocs.org/user_builds/ndcube/conda/latest/lib/python3.12/site-packages/astropy/units/format/generic.py:603: UnitsWarning: The unit 'Angstrom' has been deprecated in the VOUnit standard. Suggested: 0.1nm.
  warnings.warn(message, UnitsWarning)