NIRSPEC Data Analysis

Contents:

class nsdata.aperture

IRAF aperture class.

nsdata.bfixpix(data, badmask, n=4, retdat=False)

Replace pixels flagged as nonzero in a bad-pixel mask with the average of their nearest four good neighboring pixels.

INPUTS:

data : numpy array (two-dimensional)

badmask : numpy array (same shape as data)

OPTIONAL_INPUTS:  

n : int

number of nearby, good pixels to average over

retdat : bool

If True, return an array instead of replacing-in-place and do _not_ modify input array data. This is always True if a 1D array is input!

RETURNS:

another numpy array (if retdat is True)

TO_DO:

Implement new approach of Popowicz+2013 (http://arxiv.org/abs/1309.4224)

nsdata.cleanec(input, output, **kw)

Clean raw echelleogram files of bad pixels.

Similar to IRAF’s ‘cosmicrays’ or ‘crmed’ tasks, but this only looks along the spectrum’s dispersion direction. The less parallel your spectra are to the pixel grid, the less well this may work; in this case try increasing the ‘window’ parameter somewhat.

INPUTS:

file : (str or numpy.array)

input echellogram file to be cleaned.

DEFAULT_OPTIONS:  

dict(dispaxis=0, npasses=1, mapout=None, verbose=False,

threshold=300, nsigma=15, window=25, clobber=False, hdr=None)

nsdata.collapse_objlist(object_list, keyword, suffix='', doarray=False)

Given a LIST of identical objects, collapse the objects to form a list of only a single keyword.

aplist = nsdata.getirafap('ap1')
cenlist = nsdata.collapse_objlist(aplist, 'keyword')

If you want to get fancy, you can add a suffix to the keyword (e.g., to subscript):

newlist = nsdata.collapse_objlist(object_list, 'keyword', suffix='[0]')

nsdata.cutoffmask(filename, cutoff=[0, inf], writeto=None, clobber=True)

Create a simple mask from a FITS frame or array, based on whether

its values are above or below specified cutoff values.

INPUT:

filename: str or numpy array – frame filename or data

str – filename of FITS frame to load numpy array – data frame

OPTIONS:

cutoff : list – of form [lower, higher]; values below ‘lower’ or

above ‘higher’ will be flagged as bad.

writeto : str – filename to write output to. If None, returns the array.

RETURNS:

if writeto is None:

a 2D boolean numpy array: True for bad pixels, False for other pixels.

else:

returns None

nsdata.dark_correct(rawpre, procpre, sdark, postfix, startnum, endnum, numplaces=4)

Use similar parameters as nsdata.filelist() to correct frames for dark current. Frames will have a ‘d’ appended to their filenames.

EXAMPLE:

nsdata.dark_correct('/raw/targ_', '/proc/targ_', '/proc/sflat', '', 10, 20)

2008-06-11 16:53 IJC: Created

nsdata.darkbpmap(filelist, clipsigma=5, sigma=25, writeto=None, clobber=False, verbose=False)

Use dark frames to construct a bad pixel map based on unusually

variable pixels.

INPUT:

filelist: str, list, or 3D numpy array – dark frame filenames or data

str – IRAF-style file list (beginning with ‘@’ symbol) list – Python-style list of strs numpy array – 3D (L*M*N) stack of L dark frames

OPTIONS:

clipsigma : scalar – significance threshold for removing transient

(cosmic-ray-like) events

sigma : scalar – significance threshold for greater-than-average

pixel variability.

writeto : str – filename to write output to. If None, returns the array.

RETURNS:

if writeto is None:

a 2D boolean numpy array: True for bad pixels, False for other pixels.

else:

returns None

nsdata.dispeval(c, olim, xlim, shift=0, function='chebyshev', fullxlim=None)

Evaluate an IRAF-generated dispersion function (for now, Chebyshev only). Needs the input “C” matrix from ECIDENTIFY, as well as the limits of orders and pixels:

w = nsdata.dispeval(C, [32,37], [1,1024])

It can also be used in conjunction with nsdata.getdisp():

fn = 'ec/ecmar21s0165s'
d  = nsdata.getdisp(fn)
w  = nsdata.dispeval(d[0], d[1], d[2], shift=d[3])

SEE ALSO:nsdata.getdisp(), nsdata.interp_spec() NOTE:May not be correct for multi-order (echelle) dispersion solutions!

2008-06-24 11:48 IJC

nsdata.divdata(data, op='median', axis=None, badval=nan, returndata=False)

Take the mean/median along a specified direction and divide

the rest of the data by it.

EXAMPLE:

p   = [[1,2,3], [4,5,7]]
q1 =  nsdata.divdata(p, op='mean',   axis=0)
q2 =  nsdata.divdata(p, op='median', axis=1)

Gives:

q1: [[-1.5, -1.5, -2], [1.5, 1.5, 2]]

q2: [[ -1, 0, 1], [ -1, 0, 2]]

KEYWORDS:

*op: operation to perform; either ‘median’ (DEFAULT),

‘mean’, or ‘none’ (i.e., divide by 1)

*axis: axis along which to perform ‘op’; if None (DEFAULT),

‘op’ is performed on the entire data set as a whole.

*badval: value to replace any residual nan/inf values with.

DEFAULT is nan. Makes two passes, pre- and post-division.

*returndata: Whether to also return the data series by which the

division was performed. DEFAULT is False.

REQUIREMENTS:

Numerical analysis routines

SEE ALSO:

subdata()

nsdata.envMet(filename, tz=-10, planet=None, date=None, ignore='***')

Read the envMet.arT Keck telemetry file.

INPUT:

filename – str.

OPTIONAL_INPUTS:  

tz: offset (in hours) from GMT (HST = -10)

planet: an.planet object to compute orbital phase & HJD (if desired)

date: observation run date; grab time tags from this to use for

averaging over (not yet implemented!!!)

nsdata.file2list(filelist, prefix='', suffix='')

Convert an IRAF-like data file list of images to a Python list.

newlist = nsdata.file2list('filelist', prefix='', suffix='')

A prefix or suffix can also be appended to each filename. This is useful if, e.g., you need an explicit file extension identifier.

SEE ALSO:wfilelist(), filelist(), sfilelist()

nsdata.filelist(prefix, postfix, numlist, numplaces=4)

2008-05-27 17:12 IJC: Create a list of filenames based on a specified prefix, and a list of numeric values. You can also specify the number of digits in the filenames; default is 4.

EXAMPLE:

files = filelist(prefix, postfix, numlist, numplaces=4)

SEE ALSO:

wfilelist(), sfilelist(), file2list()

nsdata.find_features(spec, width=4, absorption=True)

Scan through a spectrum and attempt to identify features based on the sign of derivatives. Return the indices of the centers of the features.

EXAMPLE:

ind = find_features(spec, width=4, absorption=True)

2008-06-20 10:06 IJC: Created

nsdata.fix_quadnoise(*args, **kw)

Fix the 8-row coherent patterns in each quadrant of NIRSPEC using

linear least-squares after removing outliers.

INPUTS:

file – a filename or list of FITS files. The file suffix

‘.fits’ is appended if the file cannot be found. If this parameter begins with the ‘@’ (‘at’) symbol, it is interpreted as an IRAF file list.

OPTIONAL_INPUTS:  

prefix – prefix to add to the fixed files.

clobber – overwrite existing files

verbose – boolean flag for more output printed to the screen

OUTPUTS:

none

EXAMPLE:

fix_quadnoise(file, verbose=False, clobber=True) fix_quadnoise(‘mar21s0165.fits’)

nsdata.gd2jd(datestr)

Convert a string Gregorian date into a Julian date using Pylab.

If no time is given (i.e., only a date), then noon is assumed. Timezones can be given, but UTC is assumed otherwise.

EXAMPLES:

print gd2jd('Aug 11 2007')   #---------------> 2454324.5
print gd2jd('Aug 11 2007, 12:00 PST')  #-----> 2454324.29167
print gd2jd('12:00 PM, January 1, 2000')  #--> 2451545.0

REQUIREMENTS:

matplotlib

SEE ALSO:

jd2gd()

nsdata.getdisp(filename, mode='echelle')

Read the most recent dispersion function values from an IRAF dispersion database file:

D = nsdata.getdisp('ec/ecmar21s0165s')

where:

D[0] = coefficient matrix, C_mn

D[1] = [min_order, max_order]

D[2] = [min_pix, max_pix]

D[3] = pixel shift applied to the fit (units of PIXELS!)

D[4] = type: 1=chebyshev, 2=legendre

It is designed for use in conjunction with NSDATA.DISPEVAL:

fn = 'ec/ecmar21s0165s'
d  = nsdata.getdisp(fn)
w  = nsdata.dispeval(d[0], d[1], d[2], shift=d[3])

SEE ALSO:nsdata.dispeval()

2008-06-24 14:06 IJC

nsdata.getirafap(filename, filelist=False, verbose=False)

Get various parameters from an IRAF aperture file, or a list thereof

ap = getirafap(filename, filelist=False)
ap = getirafap(listname, filelist=True )

Input:filename / listname – a string filelistlist – set to True if input is a list of files Output:ap – IRAF ‘aperture’ object (or a list of these objects)

nsdata.getval(filename, key, *ext, **kw)

Get a keyword’s value from a header in a FITS file, or a list of files.

Syntax is the same as pyfits.getval:

@type filename: string or list

@param filename: input FITS file name, or list of filenames

@type key: string

@param key: keyword name

@param ext: The rest of the arguments are for extension specification.

See L{getdata} for explanations/examples.

@return: keyword value

@rtype: string, integer, or float

An extra keyword is path=’’ – a path to prepend to the filenames.

nsdata.imshow(data, x=, []y=, []aspect='auto', interpolation='nearest', cmap=None, vmin=, []vmax=[])

Version of pylab’s IMSHOW with my own defaults:

imshow(data, aspect='auto', interpolation='nearest', cmap=cm.gray, vmin=[], vmax=[])

Other IMSHOW options are default, but a new one exists:

x= and y= let you set the axes values by passing in the x and y coordinates.

nsdata.initobs(date, **kw)

Initialize variables for Nirspec data analysis.

INPUT:

date – a string of type YYYYMMMDD (e.g., 2008mar21 or 2008jun15a)

OPTIONAL_INPUT:

remote=False – changes processed-data directory as specified

within.

interp=True – whether to load wavelength-interpolated

(upsampled) spectra (if True) or the raw 1024-pixel spectra (if False)

OUTPUT:

a tuple containing the following values, in order:

planet – name of planet for use in analysis.planet()

datalist – list of data file numbers to analyse

_proc – processed data directory

wavefilename – filename of the wavelength solution

starmodelfilename – path and filename of stellar model

planetmodelfilename – path and filename of planet model

aplist – a list of the IRAF aperture filenames (for use with

nsdata.getirafap)

telluric – the FITS file spectrum of the telluric/A0V spectrum

n_aperture – number of echelle apertures for this setup

filter – NIRSPEC filter used

prefix – filename prefix

calnod – whether calibrator stars were nodded.

rowfix – list of four lists; which rows to fix in each of

four quadrants (see FIX_QUADNOISE).

nsdata.interp_spec(filename, w, w_interp, suffix='int', k=1, badval=0, clobber=True, verbose=False)

Reinterpolate a spectrum from a FITS file with a given wavelength calibration into a given wavelength grid.

result = nsdata.interp_spec(filename, wl, wl_new, k=1, badval=0, clobber=True)

Inputs:

filename – a FITS file, a list of FITS files, or a file list

(‘.fits’ optional)

wavelengths – a wavelength map of the FITS files to be reinterpolated

wl_new – the new wavelength map (preferably from nsdata.wl_grid)

Output:

result – 0 if something went wrong; 1 otherwise.

Example:

d_cor = nsdata.getdisp(‘ec/ecscifile’)

w_old = nsdata.dispeval(d[0], d[1], d[2], shift=d[3])

w_new = nsdata.wl_grid(w_old, 0.075)

result = nsdata.interp_spec(‘mar21s0161s’, w_old, w_new)

SEE ALSO:

getdisp(), dispeval(), wl_grid()

nsdata.intstr(num, numplaces=4)

A simple function to map an input number into a string padded with zeros (default 4). Syntax is: out = intstr(6, numplaces=4) –> 0006

2008-05-27 17:12 IJC: Created

nsdata.isfits(filename)

Test (CRUDELY!) whether a file is FITS format or not.

result = nsdata.isfits(filename)

Inputs:

filename – a string representing a filename

Outputs:

result – 0 if not a FITS file;

1 if the explicitly passed filename is a FITS file; 2 if conditions for 1 are not met, but

‘filename.fits’ is a FITS file.

If ‘filename’ does not exist, also checks filename+’.fits’.

nsdata.jd2gd(juldat)

Convert a numerial Julian date into a Gregorian date using Pylab.

Timezone returned will be UTC.

EXAMPLES:

print jd2gd(2454324.5)  #--> 2007-08-12 00:00:00
print jd2gd(2451545)    #--> 2000-01-01 12:00:00

SEE ALSO:

gd2jd()

nsdata.linespec(loc, ew, win, **kw)

Create a delta-function line spectrum based on a wavelength grid and a list of line locations and equivalent widths.

INPUTS:

loc – location of lines in the emission frame of reference

ew – equivalent widths of lines, in units of wavelength grid.

Positive values are emission lines.

w_in – wavelength grid in the emission frame, with values

monotonically increasing (best if it is linearly spaced)

All inputs should be lists or one-dimensional arrays of scalars

OPTIONAL_INPUTS:  

cont=None – set continuum values in the emission frame;

nearest=False – if True, use full pixels instead of partial

verbose=False – if True, print out various messages

OUTPUTS:

s – delta-function line spectrum, with a continuum level of zero

EXAMPLE:

(NEEDS TO BE UPDATED!):

w   = [2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7]
loc = [2.1, 2.35, 2.62]
ew  = [0.1, .02, .01]
s = linespec(loc, ew, w)
print s  #  --->  [0, 1, 0, 0.1, 0.1, 0, 0.08, 0.02]

NOTE:

This may give incorrect results for saturated lines.

nsdata.mederr(data, ntrials=1000, mode='median')

Return the median or mode, and the 68.3% error on that quantity, using bootstrapping.

nsdata.nAir(vaclam, T=293.15, P=100000.0, fco2=0.0004, pph2o=0.0)

Return the index of refraction of air at a given wavelength.

INPUTS:

vaclam: scalar or Numpy array

Vacuum wavelength (in microns) at which to calculate n

T : scalar

temperature in Kelvin

P : scalar

pressure in Pascals

fc02 : scalar

carbon dioxide content, as a fraction of the total atmosphere

pph2o : scalar

water vapor partial pressure, in Pascals

REFERENCE:

specfi Boensch and Potulski, 1998 Metrologia 35 133

nsdata.nAir_old(vaclam, T=288.15, P=101325.0)

Return the index of refraction of air at a given wavelength.

INPUTS:

vaclam – scalar or Numpy array – Vacuum wavelength (in

microns) at which to calculate n

T – scalar – temperature in Kelvin

P – scalar – pressure in Pascals.

NOTES:

This assumes a dry atmosphere with 0.03% CO2 by volume (“standard air”)

REFERENCE:

81st CRC Handbook (c.f. Edlen 1966)

nsdata.posofend(str1, str2)

returns the position immediately _after_ the end of the occurence of str2 in str1. If str2 is not in str1, return -1.

nsdata.preprocess(*args, **kw)

Basic processing of NIRSPEC data: set JD and HJD fields, fix bad-row ‘quadnoise’, clean cosmic rays, flatten and remove bad pixels.

INPUTS:

input =

output =

Input and output must both be specified, and must both be different files.

OPTIONAL_INPUTS:  

qfix = True

qpref = ‘’ – string to preface all quad-fixed frames

flat = None – flat field for iraf.ccdproc

mask = None – bad pixel mask for iraf.ccdproc

clobber = False – overwrite file if input and output files are same

verbose = False

cleanec = False – run nsdata.cleanec (a “poor-man’s iraf.cosmicrays”)

cthreshold = 300 – threshold for nsdata.cleanec

csigma = 20 – sigma threshold for nsdata.cleanec

cwindow = 25 – window size for nsdata.cleanec

cleancr = False – run iraf.cosmicrays

rthreshold = 300 – threshold for iraf.cosmicrays

rratio = 5 – fluxratio threshold for iraf.cosmicrays

Any inputs set to ‘None’ disable that part of the processing.

EXAMPLE:

ns.preprocess('@'+rawcal, '@'+proccal, qfix=, qpref=                        flat=, mask=, clobber=, verbose=)

Note that although the syntax is similar to Pyraf tasks, only Python Booleans or their equivalent should be used for flags (i.e., use True or False instead of ‘yes’ or ‘no’

nsdata.quadd(arg1, arg2)

Add two arguments in quadrature.

print quadd(0.1, 0.2)       -->  0.2236
print sqrt(0.1**2 + 0.2**2) -->  0.2236

nsdata.readMagiqlog(filename, dat=['guidestats'])

Read specified types of data from the Keck Magiq telemetry log.

filename : str – logfile name

dat : list of str – types of data to return

guidestats – centroid x/y, fwhm, star & skyflux

nsdata.readarm(filename)

Read a Keck .arM telemetry file.

So far, tested only with envAut.arM

nsdata.readart(filename, ignore='***')

Read a Keck .arT telemetry file.

nsdata.readfile(fn, cols=None)

Read data from an space- or tab-delimited ASCII file.

nsdata.readnstemps(filename)

Read the Keck telemetry file “nirspecTemps”

nsdata.repval(data, badval, newval)

Replace all occurrences of one value with another value. This handles nan and inf as well.

EXAMPLE:

To set all ‘nan’ values to zero, just type:

C = repval(data, nan, 0)

nsdata.setjd(filename, **kw)

Set JD and HJD fields in a FITS file, based on the UTC date and time header keywords.

INPUTS:

filename : str

OPTIONS:

ra=None : Right Ascension in decimal degrees. If None, use “RA” header key

dec=None : Declination in decimal degrees. If None, use “DEC” header key

EXAMPLE:

dict(date='date-obs', time='UTC', epoch='equinox', jd='JD', hjd='hjd',             verbose=False, ra=None, dec=None)

This does not take your position on Earth into account, so it must be less accurate than ~0.1 sec.

nsdata.sfilelist(prefix, postfix, numlist, numplaces=4, delim=', ')

2008-05-27 17:12 IJC: Create a delimited string of filenames based on a specified prefix, and a list of numeric values. You can also specify the number of digits in the numeric part filenames; default is 4. Default delimiter is a comma.

EXAMPLE:

files = sfilelist(prefix, postfix, numlist, numplaces=4, delim=',')

nsdata.showfits(filename)

2008-05-28 13:15 IJC: Routine to load and show a FITS file.

showfits(‘blah.fits’)

nsdata.strl2f(filename, strl, clobber=True, EOL='\n')

Write a list of strings to the specified filename.

INPUTS:filename: string, name of file to write to strl: list, to be written to specified file.

Returns the filename

Note:this is only designed for single-depth lists (i.e., no doubly-deep string lists).

nsdata.subab(afiles, bfiles, outfiles, clobber=False)

Take raw A & B nod files from a set of observations and subtract them; output the difference images.

INPUTS:

afiles – list of A-position filenames

bfiles – list of B-position filenames

outfiles – list of output filename

NOTE:

All input lists will be truncated to the length of the shortest input list.

For now, the output file has the header of th A-file of the pair; this is not optimal and should be fixed to record both A & B headers!

EXAMPLE:

import nsdata as ns inita = ns.initobs(‘2008jul12A’) initb = ns.initobs(‘2008jul12B’) afiles = inita[12][0][1:-4] bfiles = initb[12][0][:-4] outfiles = [(os.path.split(bfn)[0] + ‘/%s-%s’ % (os.path.split(afn)[1], os.path.split(bfn)[1])).replace(‘.fits’, ‘’) + ‘.fits’ for afn, bfn in zip(afiles, bfiles)] ns.subab(afiles, bfiles, outfiles)

nsdata.subdata(data, op='median', axis=None, returndata=False)

Take the mean/median along a specified direction and subtract it from the rest of the data.

EXAMPLE:

p   = [[1,2,3], [4,5,7]]
q1 =  nsdata.subdata(p, op='mean',   axis=0)
q2 =  nsdata.subdata(p, op='median', axis=1)

Gives:

q1: [[-1.5, -1.5, -2], [1.5, 1.5, 2]]

q2: [[ -1, 0, 1], [ -1, 0, 2]]

KEYWORDS:

*op: operation to perform; either ‘median’ (DEFAULT) or ‘mean’

*axis: axis along which to perform ‘op’; if None (DEFAULT),

‘op’ is performed on the entire data set as a whole.

*returndata: Whether to also return the data series by which the

division was performed. DEFAULT is False.

REQUIREMENTS:Numerical analysis routines SEE ALSO:divdata()

nsdata.wfilelist(prefix, postfix, numlist, numplaces=4, tempname='wfilelist_py.tmp')

2008-05-27 17:12 IJC: Create an ASCII file of a list of filenames based on a specified prefix, starting number, and ending number. You can also specify the number of digits in the filenames; default is 4.

If the file already exists, it is overwritten.

EXAMPLE:

filelist(prefix, postfix, numlist, numplaces=4, tempname='wfilelist_py.tmp')

SEE ALSO:

filelist(), sfilelist(), file2list()

nsdata.wl_grid(w, dispersion, method='log', mkhdr=False, verbose=False)

Generate a linear or logarithmic wavelength map with constant dispersion.

w_new = wl_grid(w_old, dispersion, method='log')

INPUTS:

wl_old – array; the original wavelength map. If composed

of N echelle orders of length M, should be shape N x M

dispersion – maximum desired dispersion (wavelength units

per pixel).

OPTIONAL_INPUTS:  

method – str; ‘log’ (DEFAULT) or ‘linear’ mkhdr – bool; False (DEFAULT) or True. If True, output

FITS header keys for each of the N echelle orders

verbose – bool; False (DEFAULT) or True.

OUTPUTS:

w_new – the new wavelength map, with constant dispersion hdr – list of dict, each containing FITS wavelength headers

EXAMPLE1:

w = array([1, 1.1, 1.2])
w2, h = wl_grid(w, 0.05, method='log', mkhdr=True)

Which gives: w2 = CRVAL1 * CDELT1**X, where X={0,...,M-1}

EXAMPLE2:

w = array([1.00, 1.10, 1.21])
w2, h = wl_grid(w, 0.05, method='linear', mkhdr=True)

Which gives: w2 = CRVAL1 + CDELT1*X, where X={0,...,M-1}

SEE ALSO:

interp_spec()

nsdata.write_exptime(filename, coadds='coadds')

Read ‘itime’ and ‘coadds’ from a specified file, and write into it an “exptime’ header keyword (IRAF likes it this way). If the filename does not have a ‘.fits’ extension, write_exptime will attempt to add one in order to find the file.

EXAMPLE:

nsdata.write_exptime('blah.fits')

nsdata.wspectext(inputspec)

Write a FITS-file spectrum to a column-format ASCII file.

nsdata.wspectext('filelist')   ## Not yet working
nsdata.wspectext('target.fits')
nsdata.wspectext(specdata)     ## Not yet working

Inputs:

If 3D: The first band is the data; the next-to-last is the

error; the last is the wavelength

Outputs: TBW

2008-07-07 10:44 IJC