Embedding a 1D Flare Spectra in a 2D Image#

We have a 1D spectra as a function of wavelength. To convert this to a 3D data cube (2 spatial dimensions, 1 wavelength dimension), we need to:

  1. Choose a spatial location for the flare

  2. Create an array to embed spectra in with dimensions \((N_\lambda,2N+1,2N+1)\), where \(N\) is the flare kernel size in pixel space

  3. Create a celestial WCS

  4. Interpolate spectra to constant wavelength bins

  5. Create the combined WCS

  6. Create spectral cube using blurred data and combined WCS

[1]:

import pathlib

import asdf
import ndcube
import numpy as np
import astropy.io.fits
import astropy.units as u
import astropy.wcs
from fiasco.io import Parser
from scipy.interpolate import PchipInterpolator
import sunpy.map
import sunpy.io._fits as sunpy_fits

from sunpy.coordinates import get_earth, Helioprojective
from astropy.coordinates import SkyCoord
from astropy.wcs.utils import wcs_to_celestial_frame

from synthesizAR.instruments.util import add_wave_keys_to_header

from mocksipipeline.util import read_data_cube
from mocksipipeline.instrument.configuration import moxsi_slot

[4]:

goes_class = ['a1', 'b1', 'b7', 'm1', 'm5', 'x5']
flare_spectra = []
with asdf.open('../data/reference_spectra/caspi_flare_spectra.asdf', copy_arrays=True) as af:
    for gc in goes_class:
        flare_spectra.append((gc, ndcube.NDCube(af.tree[gc]['data'], wcs=af.tree[gc]['wcs'])))
flare_spectra = ndcube.NDCollection(flare_spectra, aligned_axes=(0,))

WARNING: UnitsWarning: The unit 'Angstrom' has been deprecated in the VOUnit standard. Suggested: 0.1nm. [astropy.units.format.utils]

Interpolate each of the spectra to a common wavelength axis because we need them to be sampled regularly

[5]:

wave_min = 0.1 * u.AA
wave_max = 80 * u.AA
flare_spectra_interp = {}
for k in flare_spectra:
    energy_axis = flare_spectra[k].axis_world_coords(0)[0]
    wavelength_axis = energy_axis.to('AA', equivalencies=u.spectral())
    delta_wave = np.min(np.fabs(np.diff(wavelength_axis))) * 50  # Don't need the full resolution
    wavelength_axis_norm = u.Quantity(np.arange(wave_min.to_value('AA'), wave_max.to_value('AA'), delta_wave.to_value('AA')), 'AA')
    f_interp = PchipInterpolator(energy_axis.to_value('keV'), flare_spectra[k].data, extrapolate=True)
    data_interp = f_interp(wavelength_axis_norm.to_value('keV', equivalencies=u.spectral())) * flare_spectra[k].unit
    flare_spectra_interp[k] = (wavelength_axis_norm, data_interp)

[6]:

earth_observer = get_earth('2020-01-01')
ref_coord = SkyCoord(0, 0, unit='arcsec', frame=Helioprojective(observer=earth_observer))
array_size = (30, 130)  # this should always be even
flare_spectra_interp_cubes = {}
for k in flare_spectra_interp:
    wavelength, data = flare_spectra_interp[k]
    data_array = np.zeros(wavelength.shape+array_size) * data.unit
    data_array[:, array_size[0]//2-1:array_size[0]//2+1, array_size[1]//2-1:array_size[1]//2+1] = data[:, None, None]
    header = sunpy.map.make_fitswcs_header(
        array_size,
        ref_coord,
        scale=moxsi_slot.optical_design.spatial_plate_scale,
        unit=data.unit,
    )
    header = add_wave_keys_to_header(wavelength, header)
    wcs = astropy.wcs.WCS(header=header)
    flare_spectra_interp_cubes[k] = ndcube.NDCube(data_array, wcs=wcs, meta=header)

WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]

[7]:

results_root = pathlib.Path('/Users/wtbarnes/Documents/codes/mocksipipeline/pipeline/results/')

[9]:

for k in flare_spectra_interp_cubes:
    results_dir = results_root / f'caspi_spectra_{k}'
    results_dir.mkdir(exist_ok=True)
    sunpy_fits.write(results_dir / 'spectral_cube.fits',
                     flare_spectra_interp_cubes[k].data,
                     flare_spectra_interp_cubes[k].meta,
                     overwrite=True)

WARNING: VerifyWarning: Verification reported errors: [astropy.io.fits.verify]
WARNING: VerifyWarning: HDU 0: [astropy.io.fits.verify]
WARNING: VerifyWarning:     'NAXIS3' card at the wrong place (card 35).  Fixed by moving it to the right place (card 5). [astropy.io.fits.verify]
WARNING: VerifyWarning: Note: astropy.io.fits uses zero-based indexing.
 [astropy.io.fits.verify]

[ ]: