Source code for eispy.cube.eis_cube

'''EIS spectral cube definitions'''
import numpy as np

from astropy.io import fits
from astropy.nddata import StdDevUncertainty as sdu
from astropy.time import Time
from astropy.wcs import WCS

# from eispy.eis_spectral_cube import EISSpectralCube
# from eispy.calibration.constants import missing

from ndcube import NDCube

import re

__all__ = ['read', 'EISCube', 'EISObservation', 'EISObservationL2']


[docs]def read(filename, er_filename=None): """ Reads in a given .fits file. Parameters ---------- filename: string Location of the FITS file. er_filename: string Location of the error FITS file. Returns ------- EISObservation """ hdulist = fits.open(name=filename) errlist = fits.open(er_filename) if er_filename else None # TODO: Make sure each cube has a correct wcs. hdulist[1].verify('fix') wavelengths = [c.name for c in hdulist[1].columns if c.dim is not None] data = [hdulist[1].data[wav] for wav in wavelengths] errs = [errlist[1].data[wav] if errlist is not None else None for wav in wavelengths] is_l2 = hdulist[0].header['DATA_LEV'] == 2 cubes = [] if is_l2: vmaps, wmaps = [], [] for i in range(len(data)): window = i + 1 header = _dictionarize_header(hdulist[1].header, hdulist[0].header, window) uncertainty = sdu(errs[i]) if errlist else None if is_l2: vmaps.append(data[i][:, :, 0].T) wmaps.append(data[i][:, :, 1].T) data[i] = data[i][:, :, 2:] header['NAXIS1'], header['NAXIS2'], header['NAXIS3'] = data[i].shape wcs = WCS(header=header, naxis=3) data[i] = data[i].T cubes += [EISCube(data[i], wcs, uncertainty=uncertainty, meta=header)] primary_header = _clean(hdulist[0].header) if is_l2: return EISObservationL2(wavelengths, cubes, primary_header, vmaps, wmaps) else: return EISObservation(wavelengths, cubes, primary_header)
[docs]class EISObservation: """ A single EIS observation (that comes from a single .fits file). This class stores a series of `EISCubes`, one for each wavelength that was observed. Parameters ---------- wavelengths : list of str List of wavelengths observed. cubes : list of EISCube List of data cubes for each wavelength. primary_header : dict Primary data header for this observation. """ def __init__(self, wavelengths, cubes, primary_header): self._cubes = dict(zip(wavelengths, cubes)) self._header = primary_header def __getitem__(self, wavelength): """ Return the cube corresponding to *wavelength*. """ return self.cubes[wavelength] @property def wavelengths(self): """ List of wavelengths. """ return list(self._cubes.keys()) @property def cubes(self): """ List of data cubes. """ return self._cubes @property def obs_starttime(self): """ Observation start time. """ return Time(self._header['DATE-OBS'])
[docs]class EISObservationL2(EISObservation): """ An EIS observation derived from an L2 data product. In addition to spectral intensities, also contains rough estimates of the velocity and width maps. Parameters ---------- wavelengths : list of str List of wavelengths observed. cubes : list of EISCube List of data cubes for each wavelength. """ def __init__(self, wavelengths, cubes, primary_header, vmaps, wmaps): super().__init__(wavelengths, cubes, primary_header) self._vmaps = dict(zip(wavelengths, vmaps)) self._wmaps = dict(zip(wavelengths, wmaps))
[docs] def velocity_map(wavelength): return self._vmaps[wavelength]
[docs] def width_map(wavelength): return self._wmaps[wavelength]
[docs]class EISCube(NDCube): ''' EIS Cube subclass. References ---------- For an overview of the mission http://solarb.mssl.ucl.ac.uk/SolarB/ ''' @property def total_intensity(self): """ The intensity summed over an entire spectral window. """ data = np.sum(self.data, axis=0) wcs = self.wcs.dropaxis(2) return NDCube(data, wcs)
def _is_in_window(key, window): ''' Checks if a given key forms part of the specified spectral window. Parameters ---------- key: str The key to be validated window: int The desired window ''' end = re.findall(r'\d+$', key) # finds numbers at the end of the key if len(end) == 0: return True else: return window == int(end[0]) def _dictionarize_header(data_header, primary_header, window): ''' Combines the given FITS primary header and the bintable header for a specified window into a dictionary. Parameters ---------- data_header: astropy.io.fits.Header object, dict, or dict-like object. secondary header to be pruned for the specified window primary_header: astropy.io.fits.Header object, dict, or dict-like object. The main FITS file header window: int The window to be chosen out of the data header. ''' ph = dict(primary_header) dh = {} for key in data_header: if _is_in_window(key, window): newkey = re.sub(r'\d+$', '', key) dh[newkey] = data_header[key] dh.update(ph) dh['CRPIX3'] = 1 dh['CRVAL3'] = dh['TWAVE'] dh.pop('COMMENT', '') dh.pop('NAXIS1', '') return _clean(dh) def _clean(header): # TODO: find a way to identify cubes containing time """ Fixes non-standard or deprecated CTYPEn FITS keywords. Parameters ---------- header : astropy.io.fits.Header The header to be cleaned. """ header['CTYPE1'] = 'HPLN-TAN' # Helioprojective longitude, TAN projection header['CTYPE2'] = 'HPLT-TAN' # Helioprojective latitude, TAN projection header['CTYPE3'] = 'WAVE ' # Wavelength axis, default (TAB) projection header['NAXIS'] = 3 header['DATE-OBS'] = header.pop('DATE_OBS') # Drop the non-numbered keys that are already stored in the numbered keys for key in ['CRPIX', 'CRVAL', 'CDELT', 'CUNIT', 'CTYPE', 'CROTA']: if key in header: header.pop(key) return header