# 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


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.

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>


## 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>


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

In [ ]:


# this is our ADI function
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))


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>


Export & display the mean of the cube across wavelengths

In [11]:

dat = np.mean(outim,axis=0)
out = fits.HDUList(fits.PrimaryHDU(dat))
plt.figure(figsize=(10,10))
plt.imshow(dat)

Populating the interactive namespace from numpy and matplotlib

Out[11]:

<matplotlib.image.AxesImage at 0x116a87950>


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 = fits.HDUList(fits.PrimaryHDU(np.sum(outim[1:-1],axis=0)))



# SDI¶

SDI is then performed on the ADI result.

In [ ]:

from SDI import apply_SDI
minsep=10,
maxsep=70,
nphi=10)
out = fits.HDUList(fits.PrimaryHDU(sdi))
out = fits.HDUList(fits.PrimaryHDU(np.mean(sdi,axis=0)))
out = fits.HDUList(fits.PrimaryHDU(np.median(sdi[1:-1],axis=0)))


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[22]:

<matplotlib.image.AxesImage at 0x11e76b690>

In [ ]: