Skip to content

Creating Datacubes with IGRINS data

kfkaplan edited this page Dec 6, 2019 · 8 revisions

Written by Kyle Kaplan April 24, 2018

Previous versions of this code had major issues which are now fixed in the latest version. If you have any questions or encounter issues please contact Kyle Kaplan (kkaplan@usra.edu).

Table of Contents

Main

Introduction

Using this datacube building code, IGRINS slit-scans of extended emission line objects (i.e. nebulae) are combined into a flux calibrated 3D datacube reminiscent of channel maps used in radio and sub-mm astronomy. Given proper astrometry, these datacubes allow the output and exploration of spatially resolved emission line fluxes and flux ratios which can be combined with data from other instruments.

Install

The datacube code relies on the plotspec.py python library which is designed to process 2D IGRINS emission line data. Download the latest version of plotspec along with the datacube building demo from: https://github.com/kfkaplan/plotspec

Plotspec and the datacube code are designed to run in Python 2.7. Anaconda (https://anaconda.org) is the recommended python environment. Other python setups might work, but people have had issues in the past with other setups.

The code requires following python libraries:

  • astropy (usually comes with anaconda)
  • bottleneck (you probably need to install)
  • pylab (usually comes with python)
  • scipy (usually comes with python)
  • numpy (usually comes with python)
  • matplotlib (usually comes with python)
Make sure the proper paths and settings are set at the top of plotspec.py:
#Global variables user should set
pipeline_path = '/Volumes/home/plp/'
save_path = '/Volumes/home/results/'
#pipeline_path = '/Volumes/IGRINS_Data/plp/' #Paths for running on linux laptop
#save_path = '/Volumes/IGRINS_Data/results/'
#save_path = '/home/kfkaplan/Desktop/results/'
#pipeline_path = '/Volumes/IGRINS_data_Backup/plp/'
#save_path = '/Volumes/IGRINS_data_Backup/results/' #Define path for saving temporary files'
scratch_path = save_path + 'scratch/' #Define a scratch path for saving some temporary files
if not os.path.exists(scratch_path): #Check if directory exists
	print 'Directory '+ scratch_path + ' does not exist.  Making new directory.'
	os.mkdir(scratch_path) #If path does not exist, make directory
#default_wave_pivot = 0.625 #Scale where overlapping orders (in wavelength space) get stitched (0.0 is blue side, 1.0 is red side, 0.5 is in the middle)
default_wave_pivot = 0.85 #Scale where overlapping orders (in wavelength space) get stitched (0.0 is blue side, 1.0 is red side, 0.5 is in the middle)
set_velocity_range =100.0 # +/- km/s for interpolated velocity grid
set_velocity_res = 1.0 #Resolution of velocity grid
#slit_length = 62 #Number of pixels along slit in both H and K bands
slit_length = 100 #Number of pixels along slit in both H and K bands
block = 750 #Block of pixels used for median smoothing, using iteratively bigger multiples of block
cosmic_horizontal_mask = 3 #Number of pixels to median smooth horizontally (in wavelength space) when searching for cosmics
cosmic_horizontal_limit  = 10.0 #Number of times the data must be above it's own median smoothed self to find cosmic rays
cosmic_s2n_min = 2.5 #Minimum S/N needed to flag a pixel as a cosmic ray

A note about code performance:

There is theoretically no limit to how many IGRINS slits can be combined into one large datacube, beyond computing power and memory. Larger cubes take longer to make and eat up large amounts of RAM. The code is designed for functionality, and is not optimized for performance. If you plan to construct large datacubes, a fast computer with plenty of RAM and harddrive space is highly recommended.

Working Example

If you would like to see how the datacube building code works and how to install and run it, you can download a working example of the datacube building code created by Heeyoung Oh (hyoh@utexas.edu) here: https://drive.google.com/drive/folders/1W5Up2Rzn05WhuDOdt4c1gYulJyQPEvxZ?usp=sharing

Please refer to the enclosed readme file for details.

Parts of the datacube building code

The datacube building code consists of the following files:

datacube_demo.py
Script for creating and processing the datacube. It loads the datacube making library into memory, builds the datacube, and saves the output fits files specified by the user.
make_datacube.py
Datacube making python library. Imported by _datacube_demo.py_. User sets important variables at top of the library.
demo_datacube_input.dat
Text input file specifying all the details for the individual IGRINS slits that make up the slit-scan(s) to combine into a datacube.
plotspec.py
Python library imported by make_datacube.py that has the various classes and definitions used to analyze 2D IGRINS data.

Datacube geometry

The datacube geometry is the transformation between the spatial xy spaxel coordinates to angular coordinates on the sky. The geometry determines how to set up the positions for slit-scan observations, calculates the astrometry, and stores the astrometry in the FITS headers.

The datacube plate scale is calculated by dividing the IGRINS slit length in arcseconds by the number of vertical (y direction) spaxels found in the 2D spectra outputted by the IGRINS pipeline. The slit length depends on the telescope IGRINS was on (14.8, 9.3, and 4.9 arcseconds on McDonald, DCT, or Gemini South respectively) and is set near the top of make_datacube.py. The code automatically scales the slit width to match the slit length (a ratio of 1/15):

length_of_slit_on_sky = 15.0 #Length of the slit on the sky, in arcseconds, depends on the telescope IGRINS was on, here it is set for McDonald

The reference slit is the geometric "zero point" for the datacube. The center of the reference slit sets the 0,0 position for the spaxel x,y coordinates around which the positions of all other slits are specified relative to. The position angle (PA) of the reference slit similarly sets the zeropoint angle of the x,y grid. The reference slit's information is set near the top of make_datacube.py:

ra = 83.83205 #RA of center of reference slit pointing in decimal degrees
dec = -5.42417 #Dec of center of reference slit pointing in decimal degrees
reference_pointing_date = 20141125  #Night of pointing where the center of the slit is used to determine the RA and Dec. for the astrometry
reference_pointing_frameno = 120  #Frame number of pointing where center of the slit is used to determine the RA and Dec for the astrometry
The variables ``reference_pointing_date`` and ``reference_pointing_frameno`` refer to the reference slit's night and frame number. The reference slit must be listed in the input file. The RA and Dec. are the right ascension and declination in J2000 decimal degrees for the reference slit's center. The rest of the vital information for the reference slit is stored in the input file ``datacube_demo_input.dat``:
20141125	120	135.0	0.0	0.0	3.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
The 1st and 2nd columns store the night and frame number of the reference slit, and should match the variables reference_pointing_date and reference_pointing_frameno set at the top of make_datacube.py. The 3rd, 4th, and 5th columns store information about the geometry. The 3rd column is the position-angle (PA) of the IGRINS slit on the sky measured in degrees from the north rotated towards the east. The top of the slit would be pointing north if the PA=0 deg. The default PA for the IGRINS slit is 90 deg, or east-west. In this demo, the PA=135 deg for the reference slit. The 4th and 5th columns represent the shift in x and y from the center of the reference slit in units of arcseconds. The code will round fractions of an arcsecond to the nearest whole spaxel. This ddemo has the x and y shifts set to 0 to start from the center of the reference slit.

The positive y-direction always points towards the top of the reference slit's long axis, and the positive x-direction always points to the right of the reference slit's short axis from the perspective where the slit top is oriented upwards. The orientation of the x,y axes depends on the PA of the reference slit. This can be a bit confusing, so it is illustrated in the figures below.

From the perspective that the cardinal directions on the sky are fixed, the orientation of the xy pixel coordinates depend on the PA of the reference slit as follows:

From the perspective where the xy pixel coordinates are fixed and the reference slit center is at 0,0 with the slit top oriented upwards, the orientaiton of the cardinal directions on the sky depend on PA as follows:

After the orientation of the reference slit is set up, the rest of the slit positions in the datacube are specified relative to the reference slit in the input file datacube_demo_input.dat:

20141125	120	135.0	0.0	0.0	3.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20161230	85	45.0	0.0	0.0	7.0	99	99	5.515	5.500	1.0	0.0	H2_demo.dat	10.0	-	1
20141125	119	135.0	-1.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	123	135.0	-2.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	121	135.0	1.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	125	135.0	2.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	133	135.0	3.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1

The first five columns in the input file specify the geometric parameters for each slit in the datacube:

Column 1
Night slit was observed as YYYYMMDD. Slits can be all from the same night, or from different nights, but they must all be taken on the same telescope since the plate scale for the datacube is fixed.
Column 2
Frame number of slit. This is the first frame number specified for an observation when running the IGRINS pipeline.
Column 3
The PA of slit. The PA of the reference slit can be any angle but all other slits must be right angles to the PA of the reference slit. This is due to the fact the datacube building code paints the slits into the cube in x,y spaxel space and does not do any interpolation.
Column 4
x spaxel offset in arcseconds from the center of the reference slit.
Column 5
y spaxel offset in arcseconds from the center of the reference slit.
In this demo, the IGRINS slit is 1 arcsec wide on the 2.7 m telescope at McDonald Observatory, and each slit in the datacube is shifted 1 arcsecond from each adjacent slit. The reference slit PA=135 deg and is parallel to all the other slits, except one. Most of the the other slits have a PA=135 deg, and are offset from the reference slit is in the x direction by up to +/- 3 arcseconds. A single slit is rotated 90 degrees from the reference slit to a PA=90 deg and centered at x,y = 0,0. This perpendicular slit overlaps all the other slits in the datacube and is required to flux calibrate all the slits to each other (the relative flux calibration process is discussed later). The datacube building code will dynamically set the size of the x,y pixel grid to fit all slits in the datacube. The figures below show the geometry of the datacube specified in this demo:

Telluric Correction and Relative Flux Calibration With A0V Standard Stars

Getting scientific results from datacubes relies upon accurate emission line flux ratios, which require a good correction for the telluric absorption lines and good relative flux calibration across the entire wavelength range observed by IGRINS (H & K bands 1.45-2.45 µm). We have chosen to use A0V stars as standards because they have a well known continuum shape, broad H I absorption lines, and weak metal lines in the near-IR. The standard star for a given slit in the datacube should have been observed at a similar airmass and under similar atmospheric conditions as the science data.

For the telluric correction and relative flux calibration of individual IGRINS slits, we follow techniques similar to those used by Spextool (Vacca et al. 2003, http://irtfweb.ifa.hawaii.edu/~spex). These techniques are the same as those followed for ordinary plotspec (non-datacube building) scripts. We assume every A0V star has a continuum shape similar to Vega and construct a synthetic A0V spectrum by modifying the model Vega spectrum vegallpr25.50000resam5 by R. Kurucz (http://kurucz.havard.edu/stars.html). The model Vega spectrum is artificially reddened to match the V and B magnitudes of the A0V standard star, and then depth and broadening of the stellar H I absorption lines in the model are modified to match the lines observed in the A0V standard star spectrum. We divide the IGRINS spectrum by this synthetic modified Vega spectrum to derive the ratio of counts-to-flux at each wavelength to simultaneously get the relative flux calibration and correct for telluric absorption.

To ensure the modified model Vega spectrum matches the observed A0V standard star for each slit in the datacube, the following columns must be modified in the input file datacube_demo_input.dat:

Column 7
Specifies the frame number of the A0V standard star observed in the same night as the science observation for this slit. Multiple slit positions from the same night can share the same standard star provided the airmass and atmospheric conditions are the similar enough.
Column 9
Specifies the V magnitude of the A0V standard star to account for any reddening in its spectrum. The V magnitude can easily be looked up in SIMBAD (http://simbad.u-strasbg.fr/simbad/).
Column 10
Same as column 9 but specifies the B magnitude.
Column 11
Scale of the H I absorption line depths in the A0V standard star vs. the model Vega spectrum. Larger numbers mean deeper H I lines. Smaller numbers mean shallower H I lines. The default value = 1.0 which is equal to the depth of the Vega H I lines.
Column 12
Gaussian smoothing of the H I absorption lines in the model Vega spectrum to fit the broadening of the H I lines in the A0V standard star spectrum. Larger numbers mean more smoothing and broader H I lines. Smaller numbers mean less smoothing and narrower H I lines. The default value = 0.0 which matches the (narrow) H I lines found in the model Vega spectrum.

To see how well the H I lines in the modified Vega spectrum fit the A0V standard star spectrum, set the variable save_checks at the top of make_datacube.py = True:

save_checks = True #Save pdfs of the output of each pointing (used for checking flux calibration, ect.)
For every slit listed in the input file, setting save_checks = True creates a file called check_flux_calib_NIGHT_FRAME.pdf inside a directory named after the datacube (set at the top of make_datacube.py which itself is situated inside the save_path directoryset at the top of plotspec.py:

Each check_flux_calib_NIGHT_FRAME.pdf displays the following on the first page.

This is a zoom in of four Brackett H I absorption lines in the A0V standard star spectrum. The grey colored plot-line is the non-blaze corrected IGIRNS spectrum of the A0V star. Here the blaze curves for the individual IGRINS orders are quite evident. The narrow absorption features (narrower than the H I lines) are telluric absorption lines. The black plot-line is the same spectrum as the grey plot-line but with the H I lines from the modified Vega spectrum divided out. The blue plot-lines are the "estimated" continuum flux for the A0V star determined by fitting and averaging the what the pipeline thinks is the continuum and blaze function found in adjacent orders. The goal here is to properly adjust the depth and broadening of the H I lines (columns 11 and 12 in the input file) to best fit the H I lines in the Vega model to the actual A0V standard star spectrum. The better the fit, the closer the black plot-line match the blue plot-line.

The second page of check_flux_calib_NIGHT_FRAME.pdf looks like this:

This page shows the unmodified model Vega spectrum as the dashed blue plot-line, and the model continuum artificially reddened to match the A0V standard star as the solid black plot-line. There is nothing that really needs adjusting in the input file for this, since the A0V standard's B and V magnitudes (columns 9 and 10 in the input file) are easily obtained from SIMBAD.

You will probably need to iterate a few times, modifying columns 11 and 12 in the input file and checking each slit's check_flux_calib_NIGHT_FRAME.pdf file to best fit the Vega model's H I absorption lines to the actual A0V standard star's spectrum to obtain the best telluric correction and relative flux calibration. For multiple slits observed in the same night, you usually only need to do this for one A0V star which can be shared among multiple slits, but you will need separate standards for slits observed on different nights, or when the atmospheric conditions or airmass change significantly over slits from a single night.

Wavelength Solution, Velocities, and Data Structure

The structure of the datacube consists of three axes: RA, Dec., and velocity with the RA and Dec. stored in J2000 FK5 coordinates and the velocity stored in km/s. The "datacube" actually consists of multiple cubes, one for each spectral emission line, which could be thought of as a fourth axis. The datacube code will output 3D cubes in this structure or images collapsed along the velocity axis for any chosen spectral line (discussed later in this manual).

To build the datacube, the spectral lines must be specified in a line list. Typically, the line list is limited to the needed spectral lines. More lines will increase the computational resources necessary to build the cube. Spectral line lists are stored in the line_lists directory as .dat files which consist of two tab separated columns like the following example:

1.5884880	H I 4-14
2.12183	H2 1-0 S(1)
2.166120	H I 4-7
2.1986	[Kr III]
2.2865	[Se IV]
The first column is the wavelength of the line in µm and the second column is a short text label which is used in the datacube building code to identify that individual line. Information about the spectral lines for each slit is specified in the following columns in the datacube input file:
Column 13
Name of the spectral line list to use.
Column 14
Shift in velocity (in km/s) for all spectral lines for each slit. This is the barycentric velocity correction.
Column 15
Specifies a file that shifts the velocity of each individual spectral line. Normally this is not required so leave this column as a dash, but it can be used to correct for imperfect wavelength solutions or other abnormalities.
In order to calculate velocities, the wavelength solution must be known. In the input file:
Column 8
Frame number for the wavelength solution to use for each slit in the datacube. The wavelength solution can be calculated from the telluric lines in an A0V standard star spectrum, the OH emission lines in a sky spectrum, or even the ThAr arclamp for older data when IGRINS had its calibration unit installed. In practice, the best wavelength solutions seem to be obtained when first "registering the sky" using OH sky lines from a sky spectrum, which is used as an initial guess to calculate the wavelength solution from telluric absorption lines in the A0V standard star spectrum. Blocks of slits observed in the same night typically use the same wavelength solution.
The datacube building code takes the wavelength solution and linearly interpolates the 2D spectrum for each spectral line in each slit from position-wavelength space to position-velocity space. The default range in velocities is set to +/- 100 km/s with a resolution of 1 km/s per spaxel. This can be adjusted at the top of make_datacube.py with the following variables:
velocity_range = 100.0 #Set range of velocity axis in datacube
velocity_res = 1.0 #Set size of pixel in velocity space
Increasing the size of the velocity range or making the velocity resolution finer will proportionally increase the computational time and memory requirements to run the code and the disk space required to store the resulting fits files.

For an individual slit, interpolation for a spectral line from wavelength to velocity space looks like the following example of the 1-0 S(1) molecular hydrogen (H2) line seen in the Orion Bar where the position is the spatial axis along the IGRINS slit and the velocity is the specified velocity range that the spectral line is linearly interpolated to:

The above example can be thought of as a slice through the datacube consisting of this single slit.

Datacube flux calibration

Good science with a datacube requires the spatially resolved emission line flux ratios to vary smoothly across the entire cube, properly reflecting the actual spatially resolved flux ratios from the observed object in the sky. The throughput of the atmosphere, telescope, and instrument is never exactly the same from observation to observation, as the airmass and atmospheric conditions are constantly changing. This is especially true with passing clouds. Slits from different nights could have widely varying throughput. Even exposure times might vary. The solution to these problems is to carefully construct the datacube in such a way that the relative flux between every slit in the cube is scaled to match the other slits as they are placed into the cube one after another. We are not talking about the relative flux calibration between different wavelengths which is explained above using A0V standard stars, but the relative flux calibration between the different slits themselves in the datacube.

Even if you assume the throughput for each slit in the cube is identical, you must account differences in exposure times. This is set in the following column in the input file:

Column 6
Set the relative exposure time between each slit in the datacube. If all exposure times are identical, set to 1.0. If a slit has double the exposure time, set to 2.0. If a slit has half the exposure time, set to 0.5. The code will scale the flux in each slit by the inverse of the number placed in this column. It is possible to use this scaling to manually flux calibrate the slits, if you don't want to do it automatically.
To use the automatic flux calibration, set _flux_calibrate_ at the top of _make_datacube.py_ to True:
flux_calibrate = True #Turn on or off relative flux calibration
If no flux calibration is done, overlapping slits are still averaged together to maximize the signal-to-noise.

The flux calibration between all the slits in the datacube is done by summing up the flux in a specific emission line for a slit being placed into the cube and scaling the flux in that slit to match the flux of previously placed spaxels in the cube. The emission line should be bright enough that is observed across in all the slits that make up the datacube. The emission line should ideally be be free of telluric absorption and other sources of contamination. A range of velocities is specified over which the chosen line is summed over. A signal-to-noise threshold is specified, below which spaxels are not be used for the scaling of the fluxes. All this is set at the top of _make_datacube.py:

flux_calibration_line = '1-0 S(1)' #Name of emission line to use for flux calibration, should be something bright!
flux_calibration_velocity_range = 10.0 #Range of velocities to collapse for flux calibration in +/- km/s
flux_calibration_s2n_thresh = 3.0 #Threshold for pixel S/N to be used for flux calibration

We will start by discussing placing a single slit-scan into the datacube. First we flux calibrate Slit 2 against Slit 1. More complicated configurations and the use of blocks will be discussed later. Slit 1 is the first slit listed in the input file, and is the first one placed into the datacube. Slit 1 serves as the basis that Slit 2 will be flux calibrated to match. Slit 1 should have a high S/N and be positioned so that Slit 2 and all other subsequent slits spatially overlap with Slit 1. A good choice for the Slit 1 position and PA to have it centered and perpendicular to the other slits in the slit-scan. Slit 2 is the next slit, after Slit 1, listed in the input file. Before Slit 2 is placed into the datacube, it must be flux calibrated to match Slit 1. The code finds all spatially overlapping spaxels above the set S/N threshold for the specified emission line. The flux in the overlapping spaxels for both Slit 1 and Slit 2 is summed for the chosen emission line. The ratio of the summed fluxes is then used to scale the flux of all spaxels for all emission lines in Slit 2 to match the flux of Slit 1. Slit 3, and all subsequent slits, are then placed into the cube and flux calibrated against previously placed slits following the order that the slits are listed in the input file. Each slit needs to overlap at least one previously placed slit for the flux calibration to properly work. In this way, a relative flux calibrated datacube is built up for a single slit scan. You must be mindful of the order you list the slits in the input file.

Sets of slits can be organized into groups or blocks. Each block is treated as it's own mini-datacube and the slits within a block are flux calibrated against each other, independent of the slits in other blocks. Once all blocks are individually flux calibrated, the blocks themselves are placed into the datacube in order of their block numbers, and flux calibrated against each other, similar to how the individual slits are flux calibrated. Each block must spatially overlap, at least a little bit, previously placed blocks for the flux calibration to work properly. Blocks are a useful way to combine multiple overlapping slit-scans (even those taken on different nights), as illustrated in the diagram below. In this way. multiple slit-scans can be combined to create large datacubes.

To turn on using blocks for the flux calibration, set the use_blocks variable at the top of make_datacube.py to True. Set it to False to not use blocks:

use_blocks = False #Use blocks for flux calibration (set in input file)?  If not using blocks, turn it off here to speed up datacube building
Note that using blocks is more processor and memory intensive so if you are not using them, set _use_blocks_ to False to speed up running the code.

The block number each slit belongs to is specified in the following column of the input file:

Column 16
Flux calibration block number that each slit belongs to. All the slits in a single block are flux calibrated against each other in the order listed in the input file, and then the blocks themselves are flux calibrated against one another in the order of their block numbers. If there are no separate blocks, you usually just leave the block number equal to 1 in this column.
It is possible to get creative with how you place slits into the datacube to do your flux calibration. Be mindful of the slit order and block numbers and make sure there is a enough S/N in overlapping spaxels to achieve a good relative flux calibration across the entire datacube.

Run the datacube building code and save your FITS files

Once everything is set up in plotspec.py, make_datacube.py, and demo_datacube_input.dat, it is time to run the code and create output fits files. Running the datacube building code is controlled with the script datacube_demo.py. The example demo script is not very long. We will examine each line of code and describe what it does.

Import the needed python libraries. Here make_datacube.py is imported as cubelib.

#Test demo script for make_datacube.py library
from scipy import *
import make_datacube as cubelib #Import library to make datacubes

Set the working directory to save all output fits files and the velocity range (in km/s) over which to sum and collapse the 3D emission line data stored in the datacube into 2D images.

workdir =  '/Volumes/IGRINS_Data/datacube_demo/' #Set to where you want to save resulting fits files
vrange = [-10.0,10.0] #Velocity range

Create the datacube object which stores everything about the datacube.

cube = cubelib.data() #Create datacube object

Interpolate and fill in small spatial gaps between slits in the datacube. Comment out to turn off the gap filling.

cube.fill_gaps() #Fill in nans, this is optional, comments out if you don't want to do this

Save a fits file storing a "cube" (i.e. channel map) of the 3D information (RA, Dec., velocity) for a specified emission line. The data structure is described in more detail earlier in this manual under "Wavelength Solution, Velocities, and Data Structure". The output fits file can be opened and viewed in DS9 and other datacube visualization software. The example below saves the 3D data for the H2 1-0 S(1) line.

cube.savecube('1-0 S(1)', workdir+'1-0_S(1)_cube.fits') #Save datacube of an emission line

Save 2D image (RA, Dec.) for the specified emission line. The velocity axis is collapsed by summed over the specified velocity range to create the image. Otherwise, the format of the output fits file is identical to the 3D one described above. The example below saves an image of the 2D flux for the H2 1-0 S(1) line. cube.saveimage('1-0 S(1)', workdir+'1-0_S(1)_img.fits', vrange=vrange) #Save image of an emission line in the velocity range "vrange", here set to +/- 10 km/s

Same as 2D image above (RA, Dec.) but saves the line flux ratio map between the two specified emission lines instead of the flux of a single line. The example below save the H2 2-1 S(1)/1-0 S(1) line ratio.

cube.saveratio('2-1 S(1)', '1-0 S(1)', vrange=vrange,  fname=workdir+'ratio_21s1_10s1.fits') #Test save ratio maps

An arbitrary number of cubes or 2D images can be outputted in this way, especially if many emission lines are used. Once the script datacube_demo.py has been properly set up, it is run like any other python program.

python demo_datacube.py

Extra Settings

Below are descriptions of extra settings that can be found near the top of make_datacube.py.

Apply continuum subtraction to the 2D spectrum for each slit. Useful to remove continuum sources such as stars or background scattered light contaminating the emission lines. The continuum is fit and subtracted with a 2D running median filter along the wavelength axis of each slit.

subtract_continuum = False #Turn on or off continuum subtraction

Fill nan pixels along the 2D spectrum of each slit. Useful for filling holes created by masked bad pixels or cosmic rays along the velocity axis of the datacubes. The nan valued pixels are replaced by the median of nearby pixels along the wavelength axis of a slit's 2D spectrum.

fill_nans = True #Turn on or off filling nans

Quick Reference

Set various parameters in make_datacube.py

#~~~~~~~~~~~~~~~MODIFY THESE PARAMETERS FOR IGRINS DATA~~~~~~~~~~~~~~~~
save.name('Datacube Name') #Name to save results as
#~~~~~~~~~~~~~~~~~~~~~READ INPUT FILE~~~~~~~~~~~~~~~~~~~~~~~~~~~
input_file = 'demo_datacube_input.dat'
input = ascii.read(input_file)
n_pointings = input['Date'].size
#~~~~~~~~~~~~~~~MODIFY THESE PARAMETERS FOR ASTROMETRY~~~~~~~~~~~~~~~~
velocity_range = 100.0 #Set range of velocity axis in datacube
velocity_res = 1.0 #Set size of pixel in velocity space
ra = 83.83205 #RA of center of reference slit pointing in decimal degrees
dec = -5.42417 #Dec of center of reference slit pointing in decimal degrees
reference_pointing_date = 20141125  #Night of pointing where the center of the slit is used to determine the RA and Dec. for the astrometry
reference_pointing_frameno = 120  #Frame number of pointing where center of the slit is used to determine the RA and Dec for the astrometry
flux_calibrate = True #Turn on or off relative flux calibration
subtract_continuum = False #Turn on or off continuum subtraction
use_blocks = False #Use blocks for flux calibration (set in input file)?  If not using blocks, turn it off here to speed up datacube building
fill_nans = True #Turn on or off filling nans
save_checks = True #Save pdfs of the output of each pointing (used for checking flux calibration, ect.)
flux_calibration_line = '1-0 S(1)' #Name of emission line to use for flux calibration, should be something bright!
flux_calibration_velocity_range = 10.0 #Range of velocities to collapse for flux calibration in +/- km/s
flux_calibration_s2n_thresh = 3.0 #Threshold for pixel S/N to be used for flux calibration
length_of_slit_on_sky = 15.0 #Length of the slit on the sky, in arcseconds, depends on the telescope IGRINS was on, here it is set for McDonald

Setup input file demo_datacube_input.dat

## This file contains all the necessary information to build a datacube from many different IGRINS pointings
# It can easily handle observations taken on different nights, but the user must take care to ensure the parameter
# for each pointing is correct.
#
# 1: Date = civil date of observation in YYYMMDD (ex. 20150722), name of directory for raw data files
# 2: Frame # = number of first frame of a single pointing
# 3: PA = Position angle, the PA of the first frame is the PA of the entire datacube, other PAs can only be different by 90 degrees
# 4: X shift = Shift in x-direction perpendicular to the first slit, in arcseconds, typically this is in 1 arcsecond steps
# 5: Y shift = Shift in y-direction parallel to the first slit, in arcseconds, typically 0 but can change depending on pointing error or multiple maps being combined
# 6: Exp = Relative exposure, use to correct for differences in exposure time, or for relative flux calibration, Example:: 2.0 means exposure time was double the rest of the
# 7: Std # = number of frame for A0V standard star to use
# 8: Wave # = number of frame for wavelength solution (typically you want to use the sky frame used wavelength solution correction by OH lines, but ThAr frame can be used)
# 9: A0V V = V magnitude of A0V standard star, used for reddening correction
# 10: A0V B = B magnitude of A0V standard star, used for reddening correction
# 11: A0V HI Scale = Scale of H I absorption lines in synthetic A0V stanadard star spectrum relative to Vega, check output "check_flux_calib.pdf" in results directory
# 12: A0V smoothing = Gaussian smoothing of H I absorption lines in synthetic A0V stanadard star spectrum in wavelength space, check output "check_flux_calib.pdf" in results directory
# 13: Spectral line file = File name of line list to use
# 14: Delta v = Shift in velocity (relative to observatory) of spectral lines for target
# 15: v shift file = File storing shifts of individual lines, can be used to correct imperfect wavelength calibration, set to '-' to ignore this
# 16: Block definining how independant flux calibrations are done.  Blocks are first flux calibrated against the pointings inside each block then the blocks are flux calibrated against each other, in the input order
#
#1			2	3	4	5	6	7	8		9		10		11	12					13	14 15 16
Date 	FrameNo	PA	X_Shift	Y_Shift	Exp	StdNo WaveNo	A0V_V	A0V_B	A0V_HI_Scale	A0V_Smooth	Spec_Line_File	Delta_V	V_Shift_File	Block
#Actual datacube
20141125	120	135.0	0.0	0.0	3.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20161230	85	45.0	0.0	0.0	7.0	99	99	5.515	5.500	1.0	0.0	H2_demo.dat	10.0	-	1
20141125	119	135.0	-1.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	123	135.0	-2.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	121	135.0	1.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	125	135.0	2.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1
20141125	133	135.0	3.0	0.0	1.0	126	118	6.410	6.412	0.85	1.0	H2_demo.dat	8.0	-	1

Setup script demo_datacube.py to run datacube building code and save results

#Test demo script for make_datacube.py library
from scipy import *
import make_datacube as cubelib #Import library to make datacubes

workdir =  '/Volumes/IGRINS_Data/datacube_demo/' #Set to where you want to save resulting fits files
vrange = [-10.0,10.0] #Velocity range

#Demo of saving files from datacube
cube = cubelib.data() #Create datacube object
cube.fill_gaps() #Fill in nans, this is optional, comments out if you don't want to do this
cube.savecube('1-0 S(1)', workdir+'1-0_S(1)_cube.fits') #Save datacube of an emission line
cube.saveimage('1-0 S(1)', workdir+'1-0_S(1)_img.fits', vrange=vrange) #Save image of an emission line in the velocity range "vrange", here set to +/- 10 km/s
cube.saveratio('2-1 S(1)', '1-0 S(1)', vrange=vrange,  fname=workdir+'ratio_21s1_10s1.fits') #Test save ratio maps

Run the script:

python demo_datacube.py