First hack at resampling an NDCube

ndcube_resampling.py

import numpy as np
import matplotlib.pyplot as plt

from ndcube import NDCube
from astropy.wcs import WCS
from reproject import reproject_interp


data = np.random.rand(3, 3, 30)

# Replace some data to make it look good on a graph
data[0, 0, :] = [0.01, 0.0, 0.01, 0.03, 0.04, 0.07, 0.12, 0.15, 0.28,
                 0.3, 0.52, 0.55, 0.78, 0.73, 0.87, 0.93, 0.98, 0.67, 0.69,
                 0.69, 0.42, 0.37, 0.3, 0.15, 0.15, 0.07, 0.05, 0.02, 0.0, 0.01]

wcs = WCS(naxis=3)
wcs.wcs.ctype = 'WAVE', 'HPLT-TAN', 'HPLN-TAN'
wcs.wcs.cunit = 'Angstrom', 'deg', 'deg'
wcs.wcs.cdelt = 20, 0.5, 0.4
wcs.wcs.crpix = 1, 2, 2
wcs.wcs.crval = 5000, 0.5, 1
wcs.wcs.cname = 'wavelength', 'HPC lat', 'HPC lon'

cube = NDCube(data, wcs=wcs)
print(cube.data.shape)    # Prints (3, 3, 30)

SCALING_FACTOR = 10

# Prepare a WCS to be used for resampling
# CDELT1 is modified
resampled_wcs = WCS(naxis=3)
resampled_wcs.wcs.ctype = 'WAVE', 'HPLT-TAN', 'HPLN-TAN'
resampled_wcs.wcs.cunit = 'Angstrom', 'deg', 'deg'
resampled_wcs.wcs.cdelt = 20 / SCALING_FACTOR, 0.5, 0.4
resampled_wcs.wcs.crpix = 1, 2, 2
resampled_wcs.wcs.crval = 5000, 0.5, 1
resampled_wcs.wcs.cname = 'wavelength', 'HPC lat', 'HPC lon'

# Resample using reproject, using the new WCS as a base
resampled_data, _ = reproject_interp(cube, output_projection=resampled_wcs,
                                     shape_out=(3, 3, 30 * SCALING_FACTOR), order='bicubic')

resampled_cube = NDCube(resampled_data, wcs=resampled_wcs)
print(resampled_cube.data.shape)    # Prints (3, 3, 300)

# Extract X and Y from the cube to plot
wavelengths = cube.axis_world_coords('wl')[0].Angstrom
values = cube[0, 0, :].data

resampled_wavelengths = resampled_cube.axis_world_coords('wl')[0].Angstrom
resampled_values = resampled_cube[0, 0, :].data

# Plot and compare
plt.subplots(2, 2, figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Value')
plt.title('Normal')
plt.plot(wavelengths, values)

plt.subplot(1, 2, 2)
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Value')
plt.title('Resampled')
plt.plot(resampled_wavelengths, resampled_values)

plt.savefig('Resampling Test.jpg', dpi=180)
plt.show()