Post-processing methods

This page explores a few possibilities on how to post-process the datacubes from CHARIS to extract your science scene. The methods shown here are not yet shared.

You first need to reduce all of your raw data files as described in the reduction section. For example, this command will do this for a particular data set. Here, the data is in the folder “data/HD1160/”, and we reduced everything in the “data/HD1160/reduced/” folder.

extractcube ../*.fits HD1160_broad.ini

A typical post-processing sequence would be:

  1. Reduce all the cubes with extractcube, e.g.:

    extractcube ../CRSA1500???.fits HD1160_broad.ini
    
  2. Recenter all the cubes:

    python get_centroid.py CRSA*_cube.fits
    
  3. Apply ADI:

    python ADI.py CRSA*_recen.fits
    
  4. Apply SDI on the cube obtained from ADI:

    python SDI.py ADI_cube.fits
    

While the functions that are distributed work as command-line tools, all of them can also be called in a regular Python environment, as shown below.

ADI

In [14]:
# inline plotting of figures
%pylab inline --no-import-all
matplotlib.rcParams['image.origin'] = 'lower'
matplotlib.rcParams['image.interpolation'] = 'nearest'
Populating the interactive namespace from numpy and matplotlib
In [1]:
charisfolder = 'data/HD1160_broad/reduced/'
In [4]:
import glob,re
from astropy.io import fits
filenames = glob.glob(charisfolder+'CRSA*_cube.fits')

Stacked cubes before recentering:

In [20]:
stacked_cube = np.zeros(fits.open(filenames[0])[1].data.shape)
for i in range(len(filenames)):
    stacked_cube+=fits.open(filenames[i])[1].data
plt.figure(figsize=(10,10))
plt.imshow(stacked_cube[12,25:-25,25:-25])
Out[20]:
<matplotlib.image.AxesImage at 0x11e625ad0>
_images/PostProcessing_11_1.png

Recentering

First, one needs to recenter all the cubes with respect to each other. The following function recreates a version of the cubes which is properly recentered, and appends “_recen” to the original filename. This uses the satellite spots.

In [ ]:
from get_centroid import allcentroid
allcentroid(filenames,prefix=charisfolder)

Display slice from stacked recentered cube

In [23]:
recen_cubenames = glob.glob(charisfolder+'CRSA*_recen.fits')
stacked_cube = np.zeros(fits.open(recen_cubenames[0])[0].data.shape)
for i in range(len(recen_cubenames)):
    stacked_cube+=fits.open(recen_cubenames[i])[0].data
plt.figure(figsize=(10,10))
plt.imshow(stacked_cube[12,25:-25,25:-25])
Out[23]:
<matplotlib.image.AxesImage at 0x11b64c4d0>
_images/PostProcessing_16_1.png

Apply ADI

This uses a Principal Components Analysis approach to clean up the cubes.

In [ ]:

# this is our ADI function
from ADI import apply_ADI
outim = apply_ADI(filenames,     # List of filenames containing the aligned cubes and their variane
                    minsep=10,   # Radial separations at which to suppress weights
                    maxsep=70,   # in low-rank approximation to the PSF
                    nvec=10,     # Number of principal components to fit at each wavelength
                    ivar=True)   # whether we use the inverse variance weights or not

# export the new cube which has been cleaned by ADI
hdr = fits.open(new_cubenames[0])[0].header # this is just so we have a header to attach to the data
out = fits.HDUList(fits.PrimaryHDU(outim,hdr))
out.writeto(charisfolder+'ADI_recen.fits',clobber=True)

Display same slice as before to show how well the PSF was cleaned

In [13]:
plt.figure(figsize=(10,10))
plt.imshow(outim[12])
Out[13]:
<matplotlib.image.AxesImage at 0x116923210>
_images/PostProcessing_21_1.png

Export & display the mean of the cube across wavelengths

In [11]:
dat = np.mean(outim,axis=0)
out = fits.HDUList(fits.PrimaryHDU(dat))
out.writeto(charisfolder+'ADI_recen_mean.fits',clobber=True)
plt.figure(figsize=(10,10))
plt.imshow(dat)
Populating the interactive namespace from numpy and matplotlib
Out[11]:
<matplotlib.image.AxesImage at 0x116a87950>
_images/PostProcessing_23_2.png

Also exporting the median and the sum, just to check things

In [ ]:
out = fits.HDUList(fits.PrimaryHDU(np.median(outim[1:-1],axis=0)))
out.writeto(charisfolder+'ADI_recen_med.fits',clobber=True)
out = fits.HDUList(fits.PrimaryHDU(np.sum(outim[1:-1],axis=0)))
out.writeto(charisfolder+'ADI_recen_sum.fits',clobber=True)

SDI

SDI is then performed on the ADI result.

In [ ]:
from SDI import apply_SDI
sdi = apply_SDI(charisfolder+'ADI_recen.fits',
                minsep=10,
                maxsep=70,
                nphi=10)
out = fits.HDUList(fits.PrimaryHDU(sdi))
out.writeto(charisfolder+'SDI_ADI_recen.fits',clobber=True)
out = fits.HDUList(fits.PrimaryHDU(np.mean(sdi,axis=0)))
out.writeto(charisfolder+'SDI_ADI_recen_mean.fits',clobber=True)
out = fits.HDUList(fits.PrimaryHDU(np.median(sdi[1:-1],axis=0)))
out.writeto(charisfolder+'SDI_ADI_recen_med.fits',clobber=True)

Display the same slice after ADI+SDI

In [22]:
sdi = fits.open(charisfolder+'SDI_ADI_recen.fits')[0].data
dat = np.mean(sdi,axis=0)
out = fits.HDUList(fits.PrimaryHDU(dat))
out.writeto(charisfolder+'ADI_recen_mean.fits',clobber=True)
plt.figure(figsize=(10,10))
plt.imshow(sdi[12])
Out[22]:
<matplotlib.image.AxesImage at 0x11e76b690>
_images/PostProcessing_30_1.png
In [ ]: