Source code for pyklip.instruments.SPHERE

import os, subprocess
import astropy.io.fits as fits
from astropy import wcs
import numpy as np
import scipy.ndimage as ndimage

import pyklip
from pyklip.instruments.Instrument import Data

[docs]class Ifs(Data): """ A sequence of SPHERE IFS Data. Args: data_cube: FITS file with a 4D-cube (Nfiles, Nwvs, Ny, Nx) with all IFS coronagraphic data Also read spectral cubes and 2D images that have been saved using savedata(). psf_cube: FITS file with a 3-D (Nwvs, Ny, Nx) PSF cube info_fits: FITS file with a table in the 1st ext hdr with parallactic angle info wavelenegth_info: FITS file with a 1-D array (Nwvs) of the wavelength sol'n of a cube psf_cube_size: size of the psf cube to save (length along 1 dimension) nan_mask_boxsize: size of box centered around any pixel <= 0 to mask as NaNs IWA: inner working angle of the data in arcsecs Attributes: input: Array of shape (N,y,x) for N images of shape (y,x) centers: Array of shape (N,2) for N centers in the format [x_cent, y_cent] filenums: Array of size N for the numerical index to map data to file that was passed in filenames: Array of size N for the actual filepath of the file that corresponds to the data ifs_rdp: Reduction algorithm used to obtain the input data. PAs: Array of N for the parallactic angle rotation of the target (used for ADI) [in degrees] wvs: Array of N wavelengths of the images (used for SDI) [in microns]. For polarization data, defaults to "None" IWA: a floating point scalar (not array). Specifies to inner working angle in pixels output: Array of shape (b, len(files), len(uniq_wvs), y, x) where b is the number of different KL basis cutoffs psfs: Spectral cube of size (Nwv, psfy, psfx) where psf_cube_size defines the size of psfy, psfx. psf_center: [x, y] location of the center of the PSF for a frame in self.psfs flipx: True by default. Determines whether a relfection about the x axis is necessary to rotate image North-up East left nfiles: number of datacubes nwvs: number of wavelengths """ # class initialization # Astrometric calibration: Maire et al. 2016 # North angle not used, images are only rotated by parallactic angle and pupil offset. True north correction needs to be applied to any astrometry # north_offset = -102.18 # (pupil offset + true north offset) who knows on the sign on this angle platescale = 0.007462 # Coonstructor def __init__(self, data_cube, psf_cube, info_fits, wavelength_info, keepslices=None, psf_cube_size=21, nan_mask_boxsize=9, IWA=0.15, object_name = None, disable_minimum_filter = False, zeros2nans=False,subtract_psf_background=False): super(Ifs, self).__init__() # read in the data with fits.open(data_cube) as hdulist: self._input = hdulist[0].data # If 4D cube, Nfiles, Nwvs, Ny, Nx for vigan and Nwvs, Nfiles, Nx, Ny for sphere-dc # Read headers to be saved when using savedata. Vigan's or SPHERE-DC code don't include headers but pyklip does it # parameters or for the the location of injected planets. self.prihdrs = hdulist[0].header if np.size(self.input.shape) == 4: try: self.prihdrs['HIERARCH ESO PRO REC1 ID'] # The SPHERE DC has headers with details of the reduction. We can use that information to differentiate between the two types of reduction except KeyError: self.prihdrs['HIERARCH ESO PRO REC1 ID'] = False # If the keyword is missing, then we know it is data from Arthur Vigan's reduction if self.prihdrs['HIERARCH ESO PRO REC1 ID'] == 'sph_ifs_science_dr': ifs_rdp = "sphere-dc" # Set the reduction process self.input = np.swapaxes(self.input,0,1) # Swap the axes between the wavelengths and rotations for sphere-dc elif self.prihdrs['HIERARCH ESO PRO REC1 ID'] == False or self.prihdrs['HIERARCH ESO PRO REC1 ID'] == "F": # determine which version of vigan pipeline it is, based on the info_fits file # info_fits is a table for the old IDL pipeline, but a simple vector in the new Python pipeline with fits.open(info_fits) as hdulist: if len(hdulist) > 1: # old IDL pipeline ifs_rdp = "vigan-idl" # Set the reduction process else: self.input = np.swapaxes(self.input,0,1) # Swap the axes between the wavelengths and rotations for vigan python pipeline ifs_rdp = "vigan-python" # Set the reduction process self._ifs_rdp = ifs_rdp # Store the reduction process self._filenums = np.repeat(np.arange(self.input.shape[0]), self.input.shape[1]) self.nfiles = self.input.shape[0] self.nwvs = self.input.shape[1] # collapse files with wavelengths self.input = self.input.reshape(self.nfiles*self.nwvs, self.input.shape[2], self.input.shape[3]) # zeros are nans, and anything adjacent to a pixel less than zero is 0. if not disable_minimum_filter: input_minfilter = ndimage.minimum_filter(self.input, (0, nan_mask_boxsize, nan_mask_boxsize)) self.input[np.where(input_minfilter <= 0)] = np.nan if zeros2nans: from scipy.signal import correlate2d self.input[np.where(self.input == 0)] = np.nan square = np.ones((nan_mask_boxsize,nan_mask_boxsize)) for k in range(self.input.shape[0]): tmp_input = correlate2d(self.input[k,:,:],square,mode="same") self.input[k,:,:][np.where(np.isnan(tmp_input))] = np.nan # centers are at dim/2 self._centers = np.array([[img.shape[1]/2., img.shape[0]/2.] for img in self.input]) elif np.size(self.input.shape) == 3: # If spectral data cube. self._filenums = np.zeros(self.input.shape[0]) self.nfiles = 1 self.nwvs = self.input.shape[0] # centers are at dim/2 self.centers = np.array([[img.shape[1]/2., img.shape[0]/2.] for img in self.input]) elif np.size(self.input.shape) == 2: # If 2D images like SNR maps. self._filenums = 0 self.nfiles = 1 self.nwvs = 1 self.centers = np.array([[self.input.shape[1]/2., self.input.shape[0]/2.]]) # read in the psf cube with fits.open(psf_cube) as hdulist: psf_cube = hdulist[0].data # Nwvs, Ny, Nx if np.size(psf_cube.shape) == 4: # If more than 1 PSF was taken during observation if ifs_rdp == "vigan-idl": self.psfs = np.median(psf_cube, axis=0) # Take the median of the three PSFs elif ifs_rdp == "vigan-python": self.psfs = np.median(psf_cube, axis=1) # Take the median of the three PSFs elif ifs_rdp == "sphere-dc": self.psfs = np.median(psf_cube, axis=1) # Take the median of the three PSFs else: self.psfs = psf_cube self.psfs_center = [self.psfs.shape[2]/2., self.psfs.shape[1]/2.] # (x,y) # background subtraction if subtract_psf_background: tmp_size = 65 pixelsbefore = tmp_size//2 pixelsafter = tmp_size - pixelsbefore background_tmp = np.copy(self.psfs[:, int(self.psfs_center[1]-pixelsbefore):int(self.psfs_center[1]+pixelsafter), int(self.psfs_center[0]-pixelsbefore):int(self.psfs_center[0]+pixelsafter)]) x_grid, y_grid = np.meshgrid(np.arange(tmp_size)-tmp_size//2, np.arange(tmp_size)-tmp_size//2) r_grid = np.sqrt(x_grid**2 +y_grid**2) psf_mask = np.zeros((tmp_size,tmp_size))*np.nan psf_mask[np.where((20<r_grid)*(r_grid<30))] = 1 background_tmp = background_tmp*psf_mask[None,:,:] # trim the cube pixelsbefore = psf_cube_size//2 pixelsafter = psf_cube_size - pixelsbefore self.psfs = np.copy(self.psfs[:, int(self.psfs_center[1]-pixelsbefore):int(self.psfs_center[1]+pixelsafter), int(self.psfs_center[0]-pixelsbefore):int(self.psfs_center[0]+pixelsafter)]) self.psfs_center = [psf_cube_size/2., psf_cube_size/2.] if subtract_psf_background: self.psfs = self.psfs - np.nanmedian(background_tmp,axis=(1,2))[:,None,None] # read in wavelength solution with fits.open(wavelength_info) as hdulist: self._wvs = hdulist[0].data # repeat for all Nfile cubes self._wvs = np.tile(self.wvs, self.nfiles) # read in PA info among other things. The two reduction processes handle PAs differently and therefore must be extracted differently with fits.open(info_fits) as hdulist: if ifs_rdp == "vigan-idl": metadata = hdulist[1].data self._PAs = np.repeat(metadata["PA"] + metadata['PUPOFF'], self.nwvs) self._filenames = np.repeat(metadata["FILE"], self.nwvs) elif ifs_rdp == "vigan-python": self._PAs = hdulist[0].data self._PAs = np.repeat(self.PAs, self.nwvs) self._filenames = np.repeat(data_cube, self.nwvs) elif ifs_rdp == "sphere-dc": self._PAs = -hdulist[0].data # The SPHERE DC inverts the angles and we must correct for it self._PAs = np.repeat(self.PAs, self.nwvs) self._filenames = [data_cube] self._filenames = np.repeat(data_cube, self.nwvs) # we don't need to flip x for North Up East left self.flipx = False # I have no idea self.IWA = IWA / Ifs.platescale # 0.15" IWA # Creating WCS info for SPHERE self.wcs = [] for vert_angle in self.PAs: w = wcs.WCS() vert_angle = np.radians(vert_angle) pc = np.array([[(-1)*np.cos(vert_angle), (-1)*-np.sin(vert_angle)],[np.sin(vert_angle), np.cos(vert_angle)]]) cdmatrix = pc * self.platescale /3600. w.wcs.cd = cdmatrix self.wcs.append(w) self.wcs = np.array(self.wcs) # self.wcs = np.array([None for _ in range(self.nfiles * self.nwvs)]) self._output = None # The definition of psfs_wvs requires that no wavelengths has been skipped in the input files # But it works with keepslices self.psfs_wvs = np.unique(self.wvs) self.star_peaks = np.nanmax(self.psfs,axis=(1,2)) self.dn_per_contrast = np.squeeze(np.array([self.star_peaks[np.where(self.psfs_wvs==wv)[0]] for wv in self.wvs])) if np.size(self.input.shape) == 2: self.wvs = 0 self.dn_per_contrast = 0 if keepslices is not None: self.input = self.input[keepslices,:,:] self.nfiles = self.input.shape[0]//self.nwvs # only works if keepslices select whole cubes self.nwvs = self.input.shape[0]//self.nfiles # only works if keepslices select whole cubes self.filenums = self.filenums[keepslices] self.centers = self.centers[keepslices] self.wvs = self.wvs[keepslices] self.PAs = self.PAs[keepslices] self.filenames = self.filenames[keepslices] self.dn_per_contrast = self.dn_per_contrast[keepslices] self.wcs = self.wcs[keepslices] # Required for automatically querying Simbad for the spectral type of the star. self.object_name = object_name # self.object_name = os.path.basename(data_cube).split("-")[0].split("_")[0] ################################ ### Instance Required Fields ### ################################ @property def input(self): return self._input @input.setter def input(self, newval): self._input = newval @property def ifs_rdp(self): return self._ifs_rdp @ifs_rdp.setter def ifs_rdp(self, newval): self._ifs_rdp = newval @property def centers(self): return self._centers @centers.setter def centers(self, newval): self._centers = newval @property def filenums(self): return self._filenums @filenums.setter def filenums(self, newval): self._filenums = newval @property def filenames(self): return self._filenames @filenames.setter def filenames(self, newval): self._filenames = newval @property def PAs(self): return self._PAs @PAs.setter def PAs(self, newval): self._PAs = newval @property def wvs(self): return self._wvs @wvs.setter def wvs(self, newval): self._wvs = newval @property def wcs(self): return self._wcs @wcs.setter def wcs(self, newval): self._wcs = newval @property def IWA(self): return self._IWA @IWA.setter def IWA(self, newval): self._IWA = newval @property def output(self): return self._output @output.setter def output(self, newval): self._output = newval ############### ### Methods ### ###############
[docs] def readdata(self, filepaths): """ Reads in the data from the files in the filelist and writes them to fields """ pass
[docs] def savedata(self, filepath, data, klipparams=None, filetype="", zaxis=None , more_keywords=None): """ Save SPHERE Data. Note: In principle, the function only works inside klip_dataset(). In order to use it outside of klip_dataset, you need to define the follwing attributes: dataset.output_centers = dataset.centers Args: filepath: path to file to output data: 2D or 3D data to save klipparams: a string of klip parameters filetype: filetype of the object (e.g. "KL Mode Cube", "PSF Subtracted Spectral Cube") zaxis: a list of values for the zaxis of the datacub (for KL mode cubes currently) more_keywords (dictionary) : a dictionary {key: value, key:value} of header keywords and values which will written into the primary header """ hdulist = fits.HDUList() hdulist.append(fits.PrimaryHDU(data=data,header=self.prihdrs)) # save all the files we used in the reduction # we'll assume you used all the input files # remove duplicates from list if self.ifs_rdp == "vigan-idl" or self.ifs_rdp == "vigan-python": filenames = np.unique(self.filenames) nfiles = np.size(filenames) hdulist[0].header["DRPNFILE"] = (nfiles, "Num raw files used in pyKLIP") for i, filename in enumerate(filenames): hdulist[0].header["FILE_{0}".format(i)] = filename + '.fits' # write out psf subtraction parameters # get pyKLIP revision number pykliproot = os.path.dirname(os.path.dirname(os.path.realpath(__file__))) # the universal_newline argument is just so python3 returns a string instead of bytes # this will probably come to bite me later try: pyklipver = pyklip.__version__ except: pyklipver = "unknown" hdulist[0].header['PSFSUB'] = ("pyKLIP", "PSF Subtraction Algo") hdulist[0].header.add_history("Reduced with pyKLIP using commit {0}".format(pyklipver)) hdulist[0].header['CREATOR'] = "pyKLIP-{0}".format(pyklipver) # store commit number for pyklip hdulist[0].header['pyklipv'] = (pyklipver, "pyKLIP version that was used") if klipparams is not None: hdulist[0].header['PSFPARAM'] = (klipparams, "KLIP parameters") hdulist[0].header.add_history("pyKLIP reduction with parameters {0}".format(klipparams)) # store extra keywords in header if more_keywords is not None: for hdr_key in more_keywords: hdulist[0].header[hdr_key] = more_keywords[hdr_key] # write z axis units if necessary if zaxis is not None: # Writing a KL mode Cube if "KL Mode" in filetype: hdulist[0].header['CTYPE3'] = 'KLMODES' # write them individually for i, klmode in enumerate(zaxis): hdulist[0].header['KLMODE{0}'.format(i)] = (klmode, "KL Mode of slice {0}".format(i)) hdulist[0].header['CUNIT3'] = "N/A" hdulist[0].header['CRVAL3'] = 1 hdulist[0].header['CRPIX3'] = 1. hdulist[0].header['CD3_3'] = 1. if "Spectral" in filetype: uniquewvs = np.unique(self.wvs) # do spectral stuff instead # because wavelength solutoin is nonlinear, we're not going to store it here hdulist[0].header['CTYPE3'] = 'WAVE' hdulist[0].header['CUNIT3'] = "N/A" hdulist[0].header['CRPIX3'] = 1. hdulist[0].header['CRVAL3'] = 0 hdulist[0].header['CD3_3'] = 1 # write it out instead for i, wv in enumerate(uniquewvs): hdulist[0].header['WV{0}'.format(i)] = (wv, "Wavelength of slice {0}".format(i)) center = self.output_centers[0] hdulist[0].header.update({'PSFCENTX': center[0], 'PSFCENTY': center[1]}) hdulist[0].header.update({'CRPIX1': center[0], 'CRPIX2': center[1]}) hdulist[0].header.add_history("Image recentered to {0}".format(str(center))) try: hdulist.writeto(filepath, overwrite=True) except TypeError: hdulist.writeto(filepath, clobber=True) hdulist.close()
[docs] def calibrate_output(self, img, spectral=False, units="contrast"): """ Calibrates the flux of an output image. Can either be a broadband image or a spectral cube depending on if the spectral flag is set. Args: img: unclaibrated image. If spectral is not set, this can either be a 2-D or 3-D broadband image where the last two dimensions are [y,x] If specetral is True, this is a 3-D spectral cube with shape [wv,y,x] spectral: if True, this is a spectral datacube. Otherwise, it is a broadband image. units: currently only support "contrast" w.r.t central star Return: img: calibrated image of the same shape (this is the same object as the input!!!) """ if units == "contrast": if spectral: # spectral cube, each slice needs it's own calibration numwvs = img.shape[0] print(self.dn_per_contrast.shape) print(img.shape) img /= self.dn_per_contrast[:numwvs, None, None] else: # broadband image img /= np.nanmean(self.dn_per_contrast) self.flux_units = "contrast" return img
[docs]class Irdis(Data): """ A sequence of SPHERE IRDIS Data. Args: data_cube: FITS file with a 4D-cube (Nfiles, Nwvs, Ny, Nx) with all IFS coronagraphic data Also read spectral cubes and 2D images that have been saved using savedata(). psf_cube: FITS file with a 3-D (Nwvs, Ny, Nx) PSF cube If None, psf_cube = data_cube.replace("cube_coro","cube_psf") info_fits: FITS file with a table in the 1st ext hdr with parallactic angle info If None, info_fits = data_cube.replace("cube_coro","info") wavelength_str: string to specifiy the band (e.g. "H2H3", "K1K2") psf_cube_size: size of the psf cube to save (length along 1 dimension) IWA: inner working angle of the data in arcsecs OWA: outer working angle of the data in arcsecs Attributes: input: Array of shape (N,y,x) for N images of shape (y,x) centers: Array of shape (N,2) for N centers in the format [x_cent, y_cent] filenums: Array of size N for the numerical index to map data to file that was passed in filenames: Array of size N for the actual filepath of the file that corresponds to the data PAs: Array of N for the parallactic angle rotation of the target (used for ADI) [in degrees] wvs: Array of N wavelengths of the images (used for SDI) [in microns]. For polarization data, defaults to "None" IWA: a floating point scalar (not array). Specifies to inner working angle in pixels OWA: a floating point scalar (not array). Specifies to out working angle in pixels output: Array of shape (b, len(files), len(uniq_wvs), y, x) where b is the number of different KL basis cutoffs psfs: Spectral cube of size (2, psfy, psfx) where psf_cube_size defines the size of psfy, psfx. psf_center: [x, y] location of the center of the PSF for a frame in self.psfs flipx: True by default. Determines whether a relfection about the x axis is necessary to rotate image North-up East left nfiles: number of datacubes prihdrs: SPHERE headers if reduced by Sphere data center nwvs: number of wavelengths (i.e. 2 for dual band imaging) """ # class initialization # Astrometric calibration: Maire et al. 2016 # North angle not used, images are only rotated by parallactic angle and pupil offset. True north correction needs to be applied to any astrometry # north_offset = -1.75 # who knows on the sign on this angle platescale = 0.012255 # dual band imaging central wavelengths wavelengths = {"Y2Y3": (1.022, 1.076), "J2J3": (1.190, 1.273), "H2H3": (1.593, 1.667), "H3H4": (1.667, 1.733), "K1K2": (2.110, 2.251), "B_Y": (1.043, 1.043), "B_J": (1.245, 1.245), "B_H": (1.625, 1.625), "B_Ks": (2.182, 2.182)} # Coonstructor def __init__(self, data_cube, psf_cube, info_fits, wavelength_str, psf_cube_size=21, IWA=0.08, OWA=None, keepslices=None): super(Irdis, self).__init__() # read in the data with fits.open(data_cube) as hdulist: self._input = hdulist[0].data # 4D cube, Nfiles, Nwvs, Ny, Nx # Read headers to be saved when using savedata. Vigan's code doesn't include headers but pyklip does it # parameters or for the the location of injected planets. self._filenames = [data_cube] self.prihdrs = hdulist[0].header if np.size(self.input.shape) == 4: # If 4D if ('PIXSCAL' in self.prihdrs): irdis_rdp = "sphere-dc" # Set the reduction process self.input = np.swapaxes(self.input,0,1) # Swap the axes between the wavelengths and rotations for sphere-dc else: with fits.open(info_fits) as hdulist: if len(hdulist) > 1: # old IDL pipeline irdis_rdp = "vigan-idl" # Set the reduction process else: irdis_rdp = "vigan-python" # Set the reduction process self.input = np.swapaxes(self.input,0,1) # Swap the axes between the wavelengths and rotations for vigan python pipeline self.irdis_rdp = irdis_rdp # Store the reduction process self._filenums = np.repeat(np.arange(self.input.shape[0]), self.input.shape[1]) self.nfiles = self.input.shape[0] self.nwvs = self.input.shape[1] # collapse files with wavelengths self.input = self.input.reshape(self.nfiles*self.nwvs, self.input.shape[2], self.input.shape[3]) if OWA is not None: # Trim cube to region of interest to speed up processing if (OWA / Irdis.platescale) < (0.45*self.input.shape[1]): # Only do it if the OWA in pixels is a bit less than half the image size trim_px = int((self.input.shape[1]/2.0) - (OWA / Irdis.platescale)) self.input = self.input[:, trim_px:-trim_px, trim_px:-trim_px] # centers are at dim/2 self._centers = np.array([[img.shape[1]/2., img.shape[0]/2.] for img in self.input]) elif np.size(self.input.shape) == 3: # If spectral data cube. self._filenums = np.zeros(self.input.shape[0]) self.nfiles = 1 self.nwvs = self.input.shape[0] # centers are at dim/2 self.centers = np.array([[img.shape[1]/2., img.shape[0]/2.] for img in self.input]) elif np.size(self.input.shape) == 2: # If 2D images like SNR maps. self._filenums = 0 self.nfiles = 1 self.nwvs = 1 self.centers = np.array([[self.input.shape[1]/2., self.input.shape[0]/2.]]) # read in the psf cube with fits.open(psf_cube) as hdulist: psf_cube = hdulist[0].data # Nwvs, Ny, Nx if np.size(psf_cube.shape) == 4: if irdis_rdp == "vigan-idl": self.psfs = np.median(psf_cube, axis=0) elif irdis_rdp == "vigan-python": self.psfs = np.median(psf_cube, axis=1) # Take the median of the three PSFs elif irdis_rdp == "sphere-dc": self.psfs = np.median(psf_cube, axis=1) # Take the median of the three PSFs else: self.psfs = psf_cube # multiple PSF sequences were taken. Collapse them and take the average self.psfs_center = [self.psfs.shape[2]/2., self.psfs.shape[1]/2.] # (x,y) # trim the cube pixelsbefore = psf_cube_size//2 pixelsafter = psf_cube_size - pixelsbefore self.psfs = np.copy(self.psfs[:, int(self.psfs_center[1]-pixelsbefore):int(self.psfs_center[1]+pixelsafter), int(self.psfs_center[0]-pixelsbefore):int(self.psfs_center[0]+pixelsafter)]) self.psfs_center = [psf_cube_size/2., psf_cube_size/2.] db_wvs = Irdis.wavelengths[wavelength_str] self._wvs = np.tile(db_wvs, self.nfiles) # read in PA info among other things with fits.open(info_fits) as hdulist: if irdis_rdp == "vigan-idl": metadata = hdulist[1].data self._PAs = np.repeat(metadata["PA"] + metadata['PUPOFF'], self.nwvs) self._filenames = np.repeat(metadata["FILE"], self.nwvs) elif irdis_rdp == "vigan-python": self._PAs = hdulist[0].data self._PAs = np.repeat(self.PAs, self.nwvs) self._filenames = np.repeat(data_cube, self.nwvs) elif irdis_rdp == "sphere-dc": self._PAs = -hdulist[0].data # The SPHERE DC inverts the angles and we must correct for it self._PAs = np.repeat(self.PAs, self.nwvs) self._filenames = np.repeat(data_cube, self.nwvs) # we don't need to flip x for North Up East left self.flipx = False # I have no idea self.IWA = IWA / Irdis.platescale # 0.2" IWA self.OWA = OWA / Irdis.platescale if OWA is not None else None # Creating WCS info for SPHERE self.wcs = [] for vert_angle in self.PAs: w = wcs.WCS() vert_angle = np.radians(vert_angle) pc = np.array([[(-1)*np.cos(vert_angle), (-1)*-np.sin(vert_angle)],[np.sin(vert_angle), np.cos(vert_angle)]]) cdmatrix = pc * self.platescale /3600. w.wcs.cd = cdmatrix self.wcs.append(w) self.wcs = np.array(self.wcs) # self.wcs = np.array([None for _ in range(self.nfiles * self.nwvs)]) self._output = None # The definition of psfs_wvs requires that no wavelengths has been skipped in the input files # But it works with keepslices self.psfs_wvs = np.unique(self.wvs) self.star_peaks = np.nanmax(self.psfs,axis=(1,2)) self.dn_per_contrast = np.squeeze(np.array([self.star_peaks[np.where(self.psfs_wvs==wv)[0]] for wv in self.wvs])) if np.size(self.input.shape) == 2: self.wvs = 0 self.dn_per_contrast = 0 if keepslices is not None: self.input = self.input[keepslices,:,:] self.nfiles = self.input.shape[0]//self.nwvs # only works if keepslices select whole cubes self.nwvs = self.input.shape[0]//self.nfiles # only works if keepslices select whole cubes self.filenums = self.filenums[keepslices] self.centers = self.centers[keepslices] self.wvs = self.wvs[keepslices] self.PAs = self.PAs[keepslices] self.filenames = self.filenames[keepslices] self.dn_per_contrast = self.dn_per_contrast[keepslices] self.wcs = self.wcs[keepslices] # Required for automatically querying Simbad for the spectral type of the star. self.object_name = os.path.basename(data_cube).split("_")[0] ################################ ### Instance Required Fields ### ################################ @property def input(self): return self._input @input.setter def input(self, newval): self._input = newval @property def centers(self): return self._centers @centers.setter def centers(self, newval): self._centers = newval @property def filenums(self): return self._filenums @filenums.setter def filenums(self, newval): self._filenums = newval @property def filenames(self): return self._filenames @filenames.setter def filenames(self, newval): self._filenames = newval @property def PAs(self): return self._PAs @PAs.setter def PAs(self, newval): self._PAs = newval @property def wvs(self): return self._wvs @wvs.setter def wvs(self, newval): self._wvs = newval @property def wcs(self): return self._wcs @wcs.setter def wcs(self, newval): self._wcs = newval @property def IWA(self): return self._IWA @IWA.setter def IWA(self, newval): self._IWA = newval @property def OWA(self): return self._OWA @OWA.setter def OWA(self, newval): self._OWA = newval @property def output(self): return self._output @output.setter def output(self, newval): self._output = newval ############### ### Methods ### ###############
[docs] def readdata(self, filepaths): """ Reads in the data from the files in the filelist and writes them to fields """ pass
[docs] def savedata(self, filepath, data, klipparams=None, filetype="", zaxis=None , more_keywords=None): """ Save SPHERE Data. Note: In principle, the function only works inside klip_dataset(). In order to use it outside of klip_dataset, you need to define the follwing attribute: dataset.output_centers = dataset.centers Args: filepath: path to file to output data: 2D or 3D data to save klipparams: a string of klip parameters filetype: filetype of the object (e.g. "KL Mode Cube", "PSF Subtracted Spectral Cube") zaxis: a list of values for the zaxis of the datacub (for KL mode cubes currently) more_keywords (dictionary) : a dictionary {key: value, key:value} of header keywords and values which will written into the primary header """ hdulist = fits.HDUList() hdulist.append(fits.PrimaryHDU(data=data)) # save all the files we used in the reduction # we'll assume you used all the input files # remove duplicates from list if self.irdis_rdp == "vigan-python": filenames = np.unique(self.filenames) nfiles = np.size(filenames) hdulist[0].header["DRPNFILE"] = (nfiles, "Num raw files used in pyKLIP") for i, filename in enumerate(filenames): hdulist[0].header["FILE_{0}".format(i)] = filename + '.fits' # write out psf subtraction parameters # get pyKLIP revision number pykliproot = os.path.dirname(os.path.dirname(os.path.realpath(__file__))) # the universal_newline argument is just so python3 returns a string instead of bytes # this will probably come to bite me later try: pyklipver = pyklip.__version__ except: pyklipver = "unknown" hdulist[0].header['PSFSUB'] = ("pyKLIP", "PSF Subtraction Algo") hdulist[0].header.add_history("Reduced with pyKLIP using commit {0}".format(pyklipver)) hdulist[0].header['CREATOR'] = "pyKLIP-{0}".format(pyklipver) # store commit number for pyklip hdulist[0].header['pyklipv'] = (pyklipver, "pyKLIP version that was used") if klipparams is not None: hdulist[0].header['PSFPARAM'] = (klipparams, "KLIP parameters") hdulist[0].header.add_history("pyKLIP reduction with parameters {0}".format(klipparams)) # store extra keywords in header if more_keywords is not None: for hdr_key in more_keywords: hdulist[0].header[hdr_key] = more_keywords[hdr_key] # write z axis units if necessary if zaxis is not None: # Writing a KL mode Cube if "KL Mode" in filetype: hdulist[0].header['CTYPE3'] = 'KLMODES' # write them individually for i, klmode in enumerate(zaxis): hdulist[0].header['KLMODE{0}'.format(i)] = (klmode, "KL Mode of slice {0}".format(i)) hdulist[0].header['CUNIT3'] = "N/A" hdulist[0].header['CRVAL3'] = 1 hdulist[0].header['CRPIX3'] = 1. hdulist[0].header['CD3_3'] = 1. if "Spectral" in filetype: uniquewvs = np.sort(np.unique(self.wvs)) # do spectral stuff instead # because wavelength solutoin is nonlinear, we're not going to store it here hdulist[0].header['CTYPE3'] = 'WAVE' hdulist[0].header['CUNIT3'] = "microns" hdulist[0].header['CRPIX3'] = 1. # there's only 2 wvs in DBI hdulist[0].header['CRVAL3'] = uniquewvs[0] hdulist[0].header['CD3_3'] = uniquewvs[1] - uniquewvs[0] # write it out instead center = self.output_centers[0] hdulist[0].header.update({'PSFCENTX': center[0], 'PSFCENTY': center[1]}) hdulist[0].header.update({'CRPIX1': center[0], 'CRPIX2': center[1]}) hdulist[0].header.add_history("Image recentered to {0}".format(str(center))) try: hdulist.writeto(filepath, overwrite=True) except TypeError: hdulist.writeto(filepath, clobber=True) hdulist.close()
[docs] def calibrate_output(self, img, spectral=False, units="contrast"): """ Calibrates the flux of an output image. Can either be a broadband image or a spectral cube depending on if the spectral flag is set. Args: img: unclaibrated image. If spectral is not set, this can either be a 2-D or 3-D broadband image where the last two dimensions are [y,x] If specetral is True, this is a 3-D spectral cube with shape [wv,y,x] spectral: if True, this is a spectral datacube. Otherwise, it is a broadband image. units: currently only support "contrast" w.r.t central star Return: img: calibrated image of the same shape (this is the same object as the input!!!) """ # return img if units == "contrast": if spectral: # spectral cube, each slice needs it's own calibration numwvs = img.shape[0] img /= self.dn_per_contrast[:numwvs, None, None] else: # broadband image img /= np.nanmean(self.dn_per_contrast) self.flux_units = "contrast" return img