Spectral cube modifications (apodization, interpolation etc.)

You can use the following functions to modify your cube data and export it to a new hdf5 cube (ORCS compatible) or a fits cube.

Even if possible, you may want to use a level 3 cube since level 3 cubes contains complex data which allows for perfect interpolation, apodization etc. Older cubes contains only the real spectrum.

Note that these functions can take a very long time to run :)

from orcs.process import SpectralCube
import orb.utils.io
import logging
import numpy as np

Apodization (exported to a new hdf5 cube that can be used with ORCS)

With a level 2, 2.5 or 3 cube

cube = SpectralCube('M1_2022_SN3.merged.cm1.hdf5')

cube_path = 'NGC6888_North_SN3.merged.cm1.1.0.hdf5'
# You must first make a **hardcopy** of the original cube and rename it (here as outcube.hdf5)
export_path = 'outcube.hdf5'

def apodize_cube(cube_path, export_path, coeff):
    cube = SpectralCube(cube_path)
    if not np.iscomplex(cube[1000,1000,0]):
        print('WARNING: cube is not complex, opt for a level 3 cube')

    logging.getLogger().setLevel(logging.CRITICAL)
    fakespec = cube.get_spectrum(0,0, r=0)
    with cube.open_hdf5() as f:
        with orb.utils.io.open_hdf5(export_path, 'a') as fout:
            for ii in range(cube.dimx):
                print('\r', ii, '/', cube.dimx, end='')
                for ij in range(cube.dimy):
                    ispec = cube[ii,ij,:]
                    fakespec.data = ispec
                    # here is the apodization function. You can change this function and export it
                    # to an ORCS-compatible HDF5 cube as long as you don't change the axis of the data
                    # (no interpolation)
                    fakespecA = fakespec.apodize(coeff)
                    # e.g. if you want to take the absolute value instead of the real part :
                    # fakespecA.data = np.abs(fakespecA.data).astype(complex)

                    # you may now write the new spectrum in place of the old data
                    fout['data'][ii,ij,:] = fakespecA.data
    logging.getLogger().setLevel(logging.INFO)

apodize_cube(cube_path, export_path, 2)

logging.getLogger().setLevel(logging.INFO)

With a level 1 cube

cube_path = 'M57_SN3.merged.cm1.1.0.hdf5'
# You must first make a **hardcopy** of the original cube and rename it (here as outcube.hdf5)
export_path = 'outcube_old.hdf5'

cube = SpectralCube(cube_path)

def apodize_old_cube(cube_path, export_path, coeff):
    cube = SpectralCube(cube_path)
    if not np.iscomplex(cube[1000,1000,0]):
        print('WARNING: cube is not complex, opt for a level 3 cube')
    fakespec = cube.get_spectrum(0,0, r=0)
    logging.getLogger().setLevel(logging.CRITICAL)
    with cube.open_hdf5() as f:
        with orb.utils.io.open_hdf5(export_path, 'a') as fout:
            for iquad in range(9):
                quadpath = cube.oldcube._get_hdf5_quad_data_path(iquad)
                quad = f[quadpath]
                quadout = fout[quadpath]
                print('quadrant', iquad, '/ 9 =======')
                for ii in range(quad.shape[0]):
                    print('\r', ii, '/', cube.shape[0], end='')
                    for ij in range(quad.shape[1]):
                        ispec = quad[ii,ij,:]
                        fakespec.data = ispec
                        fakespec = fakespec.apodize(coeff)
                        quadout[ii,ij,:] = fakespec.data.real
    logging.getLogger().setLevel(logging.INFO)

apodize_old_cube(cube_path, export_path, 2)

Interpolation on a new axis (here a loglinear axis)

import orb.utils.io
import orb.utils.spectrum

cube = SpectralCube('NGC6888_North_SN3.merged.cm1.1.0.hdf5')

cube_axis = cube.get_axis(1000, 1000)

# create the projection axis
loglin_axis = 10**np.linspace(np.log10(cube_axis.data[0]), np.log10(cube_axis.data[-1]), 1000)

export_path = 'outcube.fits'

# get flambda depending on the level of the cube, to export the flux calibration with the data
if not np.iscomplex(cube[1000,1000,0]):
    print('WARNING: cube is not complex, opt for a level 3 cube')

flambda = np.ones(cube.dimz, dtype=float)
if cube.get_level() >= 3:
    if 'flambda' in cube.params:
        flambda = cube.params.flambda

    if np.size(flambda) == 1:
        flambda *= np.ones(cube.dimz, dtype=float)
    elif np.size(flambda) != cube.dimz:
        logging.warning('bad flux calibration, output will not be flux calibrated')
        flambda = np.ones(cube.dimz, dtype=float)


# write data to a HDF5 file (impossible to do with a FITS file without loading the whole data in memory)
with orb.utils.io.open_hdf5('test.hdf5', 'w') as fout:
    fout.create_dataset('/data', shape=(cube.dimx, cube.dimy, len(loglin_axis)), dtype=float)
    for ii in range(cube.dimx):
        print('\r', ii, '/', cube.dimx, end='')
        for ij in range(cube.dimy):
            ispec = cube[ii,ij,:]
            newspec = orb.utils.spectrum.project(ispec, cube_axis.data, loglin_axis, 10).astype(float)
            fout['data'][ii,ij,:] = newspec


# export the HDF5 to a FITS cube
with orb.utils.io.open_hdf5('test.hdf5', 'r') as f:
    shdu = orb.utils.io.open_large_fits_file(cube, export_path)

    for iz in range(f['data'].shape[2]):
        print('\r', iz, 'Exporting frame {}'.format(iz), end='')
        shdu.write(f['data'][:,:,iz].real.astype(np.float32).T * flambda[iz])

    shdu.close()