"""
Time-series photometric reduction for ground-based IR photometry and
spectroscopy.
Assumes a single filter and a regular dither pattern (or 'staring'),
and sky (& dark) removal via subtraction of temporally adjacent frames.
"""
# 2009-12-16 08:36 IJC: Begun.
#from pyraf import iraf
import astrolib
import phot, os
import analysis as an
import numpy as np
import pdb
try:
from astropy.io import fits as pyfits
except:
import pyfits
#_iraf = "/Users/ianc/iraf/"
#iraf.load('fitsutil')
#iraf.load('astutil')
#iraf.unlearn('ccdproc')
#iraf.unlearn('imcombine')
homedir = os.path.expanduser('~')
[docs]class baseObject:
"""Empty object container.
"""
def __init__(self):
return
[docs]def initobs(obsname='20091203k', verbose=False, cal='dome', home=homedir):
"""Initialize variables for ground-based IR data analysis.
OPTIONAL INPUT:
obsname: (str) -- e.g., 20091203k, 20091203h, '20120427_chip2'
cal : str
Valid options are 'dome' and 'sky', in case two sets of
calibration files are available. If not, this option will
have no effect.
OUTPUT:
'obs', an object containing many useful values in fields.
"""
# 2009-12-14 22:48 IJC: Created.
# 2010-01-20 16:16 IJC: Added IRTF/SpeX Dec2010 observations
# 2011-11-09 08:50 IJMC: Added some of the 2011 IRTF/SpeX Gj1214b observations
# 2011-12-16 15:14 IJMC: Added cal keyword, and Subaru/MOIRCS photometry
# 2012-03-15 07:24 IJMC: Added new SpeX case, and some error-catching.
# 2012-04-28 00:37 IJMC: Added 20120427 HAT-P-12 MOIRCS chip 2 data.
# 2013-04-08 09:07 IJMC: Adding new MOSFIRE run, GJ 3470b.
# 2013-06-04 12:04 IJMC: Added FORS2 mode.
# 2013-10-02 19:35 IJMC: Updated KOI-961/Kepler-42 files.
from numpy import tile, sort
from nsdata import gd2jd
if obsname=='20120903':
planet = None
#_raw = homedir+'/proj/transit/data/20111214/'
_proc = homedir+'/proj/browndwarfs/data/120903_O2K/'
_prefix = 'proc_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(7737, 7770, 2)
skydarks = range(7833, 7853, 2) # Corresponds to flat-field frames
domeflats = range(7789, 7803, 2)
domedarks = range(7773, 7787, 2)
sci = range(6943, 7728, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20111214_IC_ch1.fits'
for badframe in [7271, 7361, 7363]:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20140505':
planet = 'HAT-P-12 b'
_raw = homedir+'/proj/transit/data/20140505/'
_proc = homedir+'/proj/transit/data/20140505_proc/'
_prefix = 'proc_m140506_'
_suffix = '.fits'
_rawprefix = 'm140506_'
_rawsuffix = '.fits'
skies = range(9998, 9999)
flats = range(2, 15)
darks = range(15, 28)
sci = range(37, 102) + range(103, 180) + range(181,199)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + '20140505_badpixelmask_v3.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20140321':
planet = 'HAT-P-12 b'
_raw = homedir+'/proj/transit/data/20140321/'
_proc = homedir+'/proj/transit/data/20140321_proc/'
_prefix = 'proc_m140322_'
_suffix = '.fits'
_rawprefix = 'm140322_'
_rawsuffix = '.fits'
skies = range(9998, 9999)
flats = range(207, 227)
darks = range(227, 247)
sci = range(8, 34) + range(36,88) + range(89, 125)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + '20140321_badpixelmask.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20140205':
planet = 'HAT-P-12 b'
_raw = homedir+'/proj/transit/data/20140205/'
_proc = homedir+'/proj/transit/data/20140205_proc/'
_prefix = 'proc_m140205_'
_suffix = '.fits'
_rawprefix = 'm140205_'
_rawsuffix = '.fits'
skies = range(9998, 9999)
flats = range(252, 302)
darks = range(306, 352)
sci = range(32, 77) + range(78, 161)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + '20140205_badpixelmask.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20130530':
planet = 'GJ 436 b'
_raw = homedir+'/proj/transit/data/20130530/'
_proc = homedir+'/proj/transit/data/20130530_proc/'
_prefix = 'proc_fors2_big_'
_suffix = '.fits'
_rawprefix = 'fors2_big_'
_rawsuffix = '.fits'
skies = range(12411, 13012, 60) #(2831, 2891, 2951, 3011)
flats = list(set(range(12351, 13002, 10)) - set(skies))
darks = range(1761, 1802, 10)
sci = range(471, 4592, 10)
print "Some serious problems with my FORS file-naming conventions. Fix it!"
system = 'vlt/fors2'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 5
badpixelmask = _raw + 'XXXX_TBD.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20130408b':
planet = 'KOI-961 c'
_raw = homedir+'/proj/transit/data/20130408/'
_proc = homedir+'/proj/transit/data/20130408_proc/'
_prefix = 'proc_m130407_'
_suffix = '.fits'
_rawprefix = 'm130408_'
_rawsuffix = '.fits'
skies = range(76, 83)
flats = range(56, 76)
darks = range(30,50)
sci = range(615, 667) + range(669,700)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'kep42_badpixelmask_20130408.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20130408a':
planet = 'GJ 3470 b'
_raw = homedir+'/proj/transit/data/20130408/'
_proc = homedir+'/proj/transit/data/20130408_proc/'
_prefix = 'proc_m130407_'
_suffix = '.fits'
_rawprefix = 'm130408_'
_rawsuffix = '.fits'
skies = range(50,56)
flats = range(10,30)
darks = range(30,50)
sci = range(99, 549) #range(99, 500) #549 #range(200, 350) #range(99, 200) #range(99, 200)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpix_10sep2012.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20130329':
planet = 'GJ 3470 b'
_raw = homedir+'/proj/transit/data/20130329/'
_proc = homedir+'/proj/transit/data/20130329_proc/'
_prefix = 'proc_N20130329S_'
_suffix = '.fits'
_rawprefix = 'N20130329S'
_rawsuffix = '.fits'
skies = range(323, 332)
flats = range(292,323)
darks = range(118, 133)
sci = range(157,292)
system = 'gemini/gmosn/r600'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpix_XXXX.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20121010':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20121010/'
_proc = homedir+'/proj/transit/data/20121010_proc/'
_prefix = 'proc_m121010_'
_suffix = '.fits'
_rawprefix = 'm121010_'
_rawsuffix = '.fits'
skies = range(191,221) + range(500,519)
flats = range(165, 170) + range(185,190) + range(577, 626)
darks = range(160, 165) + range(180, 185) + range(528, 577)
sci = range(225, 488)
system = 'keck/mosfire'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpix_10sep2012.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20121107_chip2':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20121107/'
_proc = homedir+'/proj/transit/data/20121107_proc/'
_prefix = 'proc_MCSA0019'
_suffix = '.fits'
_rawprefix = 'MCSA0019'
_rawsuffix = '.fits'
skyflats = range(3686, 3703, 2)
skydarks = range(4132, 4151, 2)
domeflats = range(4008,4027,2) # Corresponds to DOMEFLAT_ON range(3568, 3587, 2)
domedarks = range(3986,4005,2) # Corresponds to DOMEFLAT_OFF range(3542, 3561, 2)
sci = range(3710, 3985, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20120429_IC_ch2_combined.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20120429_wasp3_chip2':
planet = 'WASP-3 b'
_raw = homedir+'/proj/transit/data/20120429/'
_proc = homedir+'/proj/transit/data/20120429_proc/'
_prefix = 'proc_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(9256, 9271, 2)
skydarks = range(9342, 9401, 2)
domeflats = range(8662, 8681, 2) # Corresponds to DOMEFLAT_ON
domedarks = range(8682, 8701, 2) # Corresponds to DOMEFLAT_OFF
sci = range(8960, 9181, 2) + range(9190, 9255, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20120429_IC_ch2_combined.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
if obsname=='20120427_chip2':
planet = 'HAT-P-12 b'
_raw = homedir+'/proj/transit/data/20120427/'
_proc = homedir+'/proj/transit/data/20120427_proc/'
_prefix = 'proc_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(8540, 8575, 2)
skydarks = range(8608, 8653, 2)
domeflats = range(8192, 8211, 2) # Corresponds to DOMEFLAT_ON
domedarks = range(8164, 8183, 2) # Corresponds to DOMEFLAT_OFF
sci = range(8332, 8349, 2) + range(8360, 8539, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20120427_IC_ch2_updated.fits'
for badframe in []:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20091203k':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/raw/20091203g/'
_proc = homedir+'/proj/transit/data/proc/20091203g/'
_prefix = '03DEB'
_suffix = '.FTS'
flats = range(2,43) + range(45,52)
darks = range(52,92)
sci = range(156,620)
system = 'lick/geminib'
timezone = 'PST'
ndither = 8
procsuffix = '_df'
numdigits = 3
elif obsname=='20091203h':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/raw/20091203g/'
_proc = homedir+'/proj/transit/data/proc/20091203g/'
_prefix = '03DEA'
_suffix = '.FTS'
flats = range(2,43)
darks = range(53,113)
sci = range(177,641)
system = 'lick/geminia'
timezone = 'PST'
ndither = 8
procsuffix = '_df'
numdigits = 3
elif obsname=='20091228a':
planet = 'WASP-12 b'
_raw = homedir+'/proj/pcsa/data/raw/2009dec28/'
_proc = homedir+'/proj/pcsa/data/proc/wasp12_2009dec28/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
#sci = range(176,678)
sci = list(np.sort(np.hstack(([176],np.hstack((178+4*np.arange(125),181+4*np.arange(125)))))))
system = 'irtf/spexprism'
timezone = 'HST'
ndither = 2
procsuffix = ''
numdigits = 4
elif obsname=='20091228b':
planet = 'WASP-12 b'
_raw = homedir+'/proj/pcsa/data/raw/2009dec28/'
_proc = homedir+'/proj/pcsa/data/proc/wasp12_2009dec28/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.b.fits'
flats = [-1]
darks = [-1]
#sci = range(176,678)
sci = list(np.sort(np.hstack(([177],np.hstack((179+4*np.arange(125),180+4*np.arange(125)))))))
system = 'irtf/spexprism'
timezone = 'HST'
ndither = 2
procsuffix = ''
numdigits = 4
elif obsname=='20091230a':
planet = 'WASP-12 b'
_raw = homedir+'/proj/pcsa/data/raw/2009dec30/'
_proc = homedir+'/proj/pcsa/data/proc/wasp12_2009dec30/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
#sci = range(1,389)
sci = list(np.sort(np.hstack((3+np.arange(30)*4, 6+np.arange(30)*4, 123, 125+np.arange(66)*4,128+np.arange(66)*4))))
for badframe in [158,159,160,161,162,163,164]:
try:
sci.remove(badframe)
except:
trash = False
system = 'irtf/spexprism'
timezone = 'HST'
ndither = 2
procsuffix = ''
numdigits = 4
elif obsname=='20091230b':
planet = 'WASP-12 b'
_raw = homedir+'/proj/pcsa/data/raw/2009dec30/'
_proc = homedir+'/proj/pcsa/data/proc/wasp12_2009dec30/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.b.fits'
flats = [-1]
darks = [-1]
#sci = range(1,389)
sci = list(np.sort(np.hstack((4+np.arange(30)*4, 5+np.arange(30)*4, 124, 126+np.arange(66)*4,127+np.arange(66)*4))))
for badframe in [158,159,160,161,162,163,164]:
try:
sci.remove(badframe)
except:
trash = False
system = 'irtf/spexprism'
timezone = 'HST'
ndither = 2
procsuffix = ''
numdigits = 4
elif obsname=='20091230guide':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/raw/20091230/guidedog/'
_proc = '<noprocdata>'
_prefix = 'spc'
_suffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(1,1257)
for badframe in []:
try:
sci.remove(badframe)
except:
trash = False
system = 'irtf/guide'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2011aug05':
planet = 'GJ 1214 b'
_raw = homedir+'/proj/pcsa/data/raw/2011aug05/spec/'
_proc = homedir+'/proj/pcsa/data/proc/gj1214_2011aug05/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'filename'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(251,462)
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2011aug05cal':
planet = 'HD 161289'
_raw = homedir+'/proj/pcsa/data/raw/2011aug05/spec/'
_proc = homedir+'/proj/pcsa/data/proc/gj1214_2011aug05/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'filename'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(462,502)
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2011aug16':
planet = 'GJ 1214 b'
_raw = homedir+'/proj/pcsa/data/raw/2011aug16/spec/'
_proc = homedir+'/proj/pcsa/data/proc/gj1214_2011aug16/'
_prefix = 'spectra01_'
_suffix = '.fits'
_rawprefix = 'filename'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(1, 199)
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2011sep12':
planet = 'GJ 1214 b'
_raw = homedir+'/proj/pcsa/data/raw/2011sep12/spec/'
_proc = homedir+'/proj/pcsa/data/proc/gj1214_2011sep12/'
_prefix = 'spectra01_'
_suffix = '.fits'
_rawprefix = 'filename'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(161, 185) + range(190, 353)
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2011oct01':
planet = 'GJ 1214 b'
_raw = homedir+'/proj/pcsa/data/raw/2011oct01/spec/'
_proc = homedir+'/proj/pcsa/data/proc/gj1214_2011oct01/'
_prefix = 'spectra01_'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = range(151, 283)
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2012mar15a':
planet = 'HD 97658 b'
_raw = homedir+'/proj/pcsa/data/raw/2012mar15/spec/'
_proc = homedir+'/proj/pcsa/data/proc/2012mar15/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.a.fits'
flats = [-1]
darks = [-1]
sci = [181, 183, 186, 187, 190, 191, 194, 195, 198, 199, 201, 204, 205, 208, 209, 212, 213, 216, 217, 220]
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='2012mar15b':
planet = 'HD 97658 b'
_raw = homedir+'/proj/pcsa/data/raw/2012mar15/spec/'
_proc = homedir+'/proj/pcsa/data/proc/2012mar15/'
_prefix = 'spectra'
_suffix = '.fits'
_rawprefix = 'spc'
_rawsuffix = '.b.fits'
flats = [-1]
darks = [-1]
sci = [182, 184, 185, 188, 189, 192, 193, 196, 197, 200, 202, 203,
206, 207, 210, 211, 214, 215, 218, 219]
system = 'irtf/spexsxd'
timezone = 'HST'
ndither = 1
procsuffix = ''
numdigits = 4
elif obsname=='20111214_chip1':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20111214/'
_proc = homedir+'/proj/transit/data/20111214_proc/'
_prefix = 'proc_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(7737, 7770, 2)
skydarks = range(7833, 7853, 2) # Corresponds to flat-field frames
domeflats = range(7789, 7803, 2)
domedarks = range(7773, 7787, 2)
sci = range(6943, 7728, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20111214_IC_ch1.fits'
for badframe in [7271, 7361, 7363]:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20111214_chip2':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20111214/'
_proc = homedir+'/proj/transit/data/20111214_proc/'
_prefix = 'proc_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(7738, 7770, 2)
skydarks = range(7834, 7853, 2) # Corresponds to flat-field frames
domeflats = range(7790, 7803, 2)
domedarks = range(7774, 7787, 2)
sci = range(6944, 7729, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20111214_IC_ch2.fits'
for badframe in [7272, 7362, 7364]:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20111214_chip1lin':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20111214/'
_proc = homedir+'/proj/transit/data/20111214_proc/'
_prefix = 'proc_linear_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(7737, 7770, 2)
skydarks = range(7833, 7853, 2) # Corresponds to flat-field frames
domeflats = range(7789, 7803, 2)
domedarks = range(7773, 7787, 2)
sci = range(6943, 7728, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20111214_IC_ch1.fits'
for badframe in [7271, 7361, 7363]:
try:
sci.remove(badframe)
except:
pass
elif obsname=='20111214_chip2lin':
planet = 'WASP-12 b'
_raw = homedir+'/proj/transit/data/20111214/'
_proc = homedir+'/proj/transit/data/20111214_proc/'
_prefix = 'proc_linear_MCSA0018'
_suffix = '.fits'
_rawprefix = 'MCSA0018'
_rawsuffix = '.fits'
skyflats = range(7738, 7770, 2)
skydarks = range(7834, 7853, 2) # Corresponds to flat-field frames
domeflats = range(7790, 7803, 2)
domedarks = range(7774, 7787, 2)
sci = range(6944, 7729, 2)
system = 'subaru/moircs'
timezone = 'UTC'
ndither = 1
procsuffix = ''
numdigits = 4
badpixelmask = _raw + 'badpixmap_20111214_IC_ch2.fits'
for badframe in [7272, 7362, 7364]:
try:
sci.remove(badframe)
except:
pass
masterflat = ''
if cal=='dome':
if 'domeflats' in locals():
flats = domeflats
masterflat = _proc + 'masterdomeflat' + _suffix
if 'domedarks' in locals():
darks = domedarks
elif cal=='sky':
if 'skyflats' in locals():
flats = skyflats
masterflat = _proc + 'masterskyflat' + _suffix
if 'skydarks' in locals():
darks = skydarks
else:
print "'cal' keyword must be either 'dome' or 'sky'!"
if 'skies' not in locals():
skies = []
if masterflat=='':
masterflat = _proc + 'masterflat' + _suffix
if _rawprefix and _rawsuffix:
pass
else:
_rawprefix = _prefix
_rawsuffix = _suffix
if timezone == 'PST':
toffset = 8
elif timezone=='PDT':
toffset = 7
elif timezone=='HST':
toffset = 10
elif timezone=='UTC':
toffset = 0
fmtstr = '%0'+str(numdigits)+'g'
allindices = darks + flats + sci + skies
allindex = np.argsort(allindices)
rawall = [(_raw+_rawprefix+ (fmtstr % num) +_rawsuffix) for num in np.sort(allindices)]
types = np.array(['dark']*len(darks) + ['flat']*len(flats) + \
['sci']*len(sci + ['sky']*len(skies)))[allindex]
rawdark = [(_raw + _rawprefix + (fmtstr % num) + _rawsuffix) for num in darks]
rawflat = [(_raw + _rawprefix + (fmtstr % num) + _rawsuffix) for num in flats]
rawsci = [(_raw + _rawprefix + (fmtstr % num) + _rawsuffix) for num in sci]
rawsky = [(_raw + _rawprefix + (fmtstr % num) + _rawsuffix) for num in skies]
dsubflat = [(_proc + _prefix +(fmtstr % num) +'_d'+ _suffix) for num in flats]
dsubsci = [(_proc + _prefix + (fmtstr % num) +'_d'+ _suffix) for num in sci]
dsuball = [(_proc + _prefix + (fmtstr % num) +'_d'+ _suffix) for num in np.sort(allindices)]
procsci = [(_proc + _prefix + (fmtstr % num) +procsuffix + _suffix) for num in sci]
procsky = [(_proc + _prefix + (fmtstr % num) +procsuffix + _suffix) for num in skies]
if system=='lick/geminia':
gain = 0 # photons (i.e., electrons) per data unit
readnoise = 0 # photons (i.e., electrons) per individual read
naxis1, naxis2 = 256, 256
elif system=='lick/geminib':
gain = 0 # photons (i.e., electrons) per data unit
readnoise = 0 # photons (i.e., electrons) per individual read
naxis1, naxis2 = 256, 256
elif system.find('irtf/spex')>=0:
gain = 13 # photons (i.e., electrons) per data unit
readnoise = 50 # photons (i.e., electrons) per individual read
naxis1,naxis2 = 1024,1024
elif system.find('irtf/guide')>=0:
gain = 14.7 # photons (i.e., electrons) per data unit
readnoise = 65 # photons (i.e., electrons) per individual read
naxis1,naxis2 = 512,512
if system.find('lick/gemini')>=0:
headers = [pyfits.getheader(file) for file in rawall]
timestrs = ['%s %s' % (h['date_obs'], h['time_obs']) \
for h in headers]
jd = np.array([gd2jd(t) for t in timestrs])+toffset/24.
coadds = np.array([h['coadds'] for h in headers])
elif system.find('irtf/spex')>=0:
try:
try:
headers = [pyfits.getheader(file) for file in procsci]
except:
headers = [pyfits.getheader(file, output_verify='ignore') for file in rawsci]
timestrs = ['%s %s' % (h['date_obs'], h['time_obs']) \
for h in headers]
jd = np.array([gd2jd(t) for t in timestrs])
ra = map(phot.hms, [h['RA'] for h in headers])
dec = map(phot.dms, [h['DEC'] for h in headers])
hjd = np.array([astrolib.helio_jd(jdi - 2400000., rai, deci) + 2400000. for \
(jdi, rai, deci) in zip(jd, ra, dec)])
coadds = np.array([h['co_adds'] for h in headers])
except:
jd = None
ra = None
dec = None
hjd = None
coadds = None
headers = [dict()]
elif system.find('irtf/guide')>=0:
headers = [pyfits.getheader(file) for file in rawsci]
timestrs = ['%s %s' % (h['date_obs'], h['time_obs']) \
for h in headers]
jd = np.array([gd2jd(t) for t in timestrs])
coadds = np.array([h['co_adds'] for h in headers])
elif system.find('subaru/moircs')>=0:
try:
headers = [pyfits.getheader(file) for file in procsci]
except:
headers = [pyfits.getheader(file) for file in rawsci]
mjd = [0.5 * (h['mjd-str'] + h['mjd-end']) for h in headers]
jd = np.array(mjd) + 2400000.5
planet_obj = an.getobj(planet)
ra = planet_obj.ra * 15
dec = planet_obj.dec
hjd = np.array([astrolib.helio_jd(mjdi + 0.5, ra, dec) + \
2400000. for mjdi in mjd])
coadds = np.array([h['coadd'] for h in headers])
# Note that this gain is the average of the Chip 1 & 2 values:
gain = 3.4 # photons (i.e., electrons) per data unit
readnoise = 30 # photons (i.e., electrons) per individual read
naxis1,naxis2 = 2048, 2048
elif system.find('keck/mosfire')>=0:
try:
headers = [pyfits.getheader(file) for file in procsci]
except:
headers = [pyfits.getheader(file) for file in rawsci]
timestrs_start = ['%s %s' % (h['date-obs'], h['time-obs']) \
for h in headers]
timestrs_end = ['%s %s' % (h['date-obs'], h['time-end']) \
for h in headers]
frametime = np.array([h['truitime'] for h in headers])
numreads = np.array([h['readdone'] for h in headers])
jd_start = np.array([gd2jd(t) for t in timestrs_start])
jd_end = np.array([gd2jd(t) for t in timestrs_end])
#???midpoint = ((TIME-END - frameTime * (numreads + 1.5) - TIME_OBS) / 2)
#???jd = (jd_end - (frametime * (numreads + 1.5))/86400. - jd_start) / 2.
jd = 0.5*(jd_start + jd_end)
planet_obj = an.getobj(planet)
ra = planet_obj.ra * 15
dec = planet_obj.dec
hjd = np.array([astrolib.helio_jd(jdi - 2400000.0, ra, dec) + \
2400000. for jdi in jd])
coadds = np.array([h['coaddone'] for h in headers])
# Note that this gain is the average of the Chip 1 & 2 values:
gain = headers[0]['sysgain'] # photons (i.e., electrons) per data unit
readnoise = np.sqrt(21.**2 / (0.5*headers[0]['readdone']) + 2.2**2) # photons (i.e., electrons) per individual read
naxis1,naxis2 = 2048, 2048
elif system.find('gemini/gmos')>=0:
try:
headers = [pyfits.getheader(file) for file in procsci]
except:
headers = [pyfits.getheader(file) for file in rawsci]
timestrs_start = ['%s %s' % (h['date-obs'], h['time-obs']) \
for h in headers]
timestrs_end = ['%s %s' % (h['date-obs'], h['utend']) \
for h in headers]
itimes = np.array([h['exptime'] for h in headers])
coadds = np.ones(itimes.size)
jd_start = np.array([gd2jd(t) for t in timestrs_start])
jd_end = np.array([gd2jd(t) for t in timestrs_end])
jd = 0.5*(jd_start + jd_end)
planet_obj = an.getobj(planet)
ra = planet_obj.ra * 15
dec = planet_obj.dec
hjd = np.array([astrolib.helio_jd(jdi - 2400000.0, ra, dec) + \
2400000. for jdi in jd])
# Note that this gain is the average of the Chip 1 & 2 values:
gain = 2.4 # photons (i.e., electrons) per data unit
readnoise = 4.5 # photons (i.e., electrons) per individual read
coadds
elif system.find('vlt/fors2')>=0:
try:
headers = [pyfits.getheader(file) for file in procsci]
except:
headers = [pyfits.getheader(file) for file in rawsci]
try:
for h in headers:
h['AIRMASS'] = 0.5*(h['HIERARCH ESO TEL AIRM END'] + h['HIERARCH ESO TEL AIRM START'])
except:
pass
timestrs_start = [h['DATE'] for h in headers]
jd_start = np.array([gd2jd(t) for t in timestrs_start])
itimes = np.array([h['exptime'] for h in headers])
jd = jd_start + itimes/86400./2.
coadds = np.ones(itimes.size)
planet_obj = an.getobj(planet)
ra = planet_obj.ra * 15
dec = planet_obj.dec
hjd = np.array([astrolib.helio_jd(jdi - 2400000.0, ra, dec) + \
2400000. for jdi in jd])
# Note that this gain is the average of the Chip 1 & 2 values:
gain = headers[0]['HIERARCH ESO DET OUT1 GAIN'] # photons (i.e., electrons) per data unit
readnoise = headers[0]['HIERARCH ESO DET OUT1 RON'] # photons (i.e., electrons) per individual read
if np.all(['time' in h for h in headers]):
itimes = np.array([h['itime'] for h in headers])
elif 'exp1time' in headers[0]:
itimes = np.array([h['exp1time'] for h in headers])
elif np.all(['truitime' in h for h in headers]):
itimes = np.array([h['truitime'] for h in headers])
elif np.all(['exptime' in h for h in headers]):
itimes = np.array([h['exptime'] for h in headers])
else:
itimes = None
obs = baseObject()
obs._raw = _raw
obs._proc = _proc
obs.rawall = np.array(rawall)
obs.rawdark = np.array(rawdark )
obs.rawflat = np.array(rawflat )
obs.rawsci = np.array(rawsci )
obs.dsubflat = np.array(dsubflat)
obs.dsubsci = np.array(dsubsci )
obs.dsuball = np.array(dsuball )
obs.procsci = np.array(procsci )
obs.target = planet
obs.gain = gain
obs.readnoise = readnoise
obs._suffix = _suffix
obs._prefix = _prefix
obs.types = types
obs.jd= jd
obs.itimes= itimes
obs.coadds = coadds
try:
obs.naxis1 = headers[0]['naxis1']
obs.naxis2 = headers[0]['naxis2']
except:
pass
obs.ndither = ndither
obs.system = system
obs.headers = headers
obs.masterflat = masterflat
try:
obs.hjd = hjd
except:
obs.hjd = None
if 'skies' in locals():
obs.rawsky = rawsky
obs.procsky = procsky
try:
obs.badpixelmask = badpixelmask
except:
obs.badpixelmask = None
return obs
[docs]def red_dark(obs, ind, num=np.nan, verbose=False, sigma=3, niter=1, clobber=False):
"""Combine dark frames into a single dark frame:"""
_sdark = obs._proc + ('superdark%03g' % num) + obs._suffix
darklist = obs.rawall[ind]
framecenter = obs.naxis1/2, obs.naxis2/2
framesize = obs.naxis1, obs.naxis2
if verbose: print "dark file list>>" , darklist
darkstack = phot.subreg2(darklist, framecenter, framesize, verbose=verbose)
if verbose: print "loaded dark stack"
superdark = np.zeros(framesize,float)
for ii in range(obs.naxis1):
for jj in range(obs.naxis2):
superdark[ii,jj] = an.meanr(darkstack[:,ii,jj],nsigma=sigma, \
niter=niter)
hdr0 = pyfits.getheader(darklist[0])
hdr0.update('sup-dark', 'super-dark created by IR.py')
pyfits.writeto(_sdark, superdark, header=hdr0, clobber=clobber)
if verbose: print "Done making superdark!"
return
[docs]def red_sub(ifiles, ofiles, subfile, num=np.nan, clobber=False):
"""Subtract one file from a list of others, and save in a new location.
First two inputs are lists of strings, 'subfile' is a string.
"""
subtractor = pyfits.getdata(subfile)
for infile, outfile in zip(ifiles,ofiles):
thisin = pyfits.getdata(infile)
thishdr = pyfits.getheader(infile)
thishdr.update('subtract', '%s minus %s' % (os.path.split(infile)[-1], os.path.split(subfile)[-1]))
pyfits.writeto(outfile, thisin-subtractor, header=thishdr, clobber=clobber)
return
[docs]def red_flat(obs, ind, num=np.nan, verbose=False, sigma=3, niter=1, clobber=False):
"""Dark-correct flats, make normalized super-flat."""
_sdark = obs._proc + ('superdark%03g' % num) + obs._suffix
_sflat = obs._proc + ('superflat%03g' % num) + obs._suffix
rawflat = obs.rawall[ind]
dsubflat = obs.dsuball[ind]
nflat = len(rawflat)
framecenter = obs.naxis1/2, obs.naxis2/2
framesize = obs.naxis1, obs.naxis2
red_sub(rawflat, dsubflat, _sdark, clobber=True)
if verbose: print "dark-subtracted flat frames"
flatstack = phot.subreg(dsubflat, framecenter, framesize, verbose=verbose)
if verbose: print "read in flat stack"
flatmedians = an.amedian(flatstack.reshape((nflat,np.prod(framesize))),1).reshape((nflat,1,1))
superflat = np.zeros(framesize,float)
for ii in range(obs.naxis1):
for jj in range(obs.naxis2):
superflat[ii,jj] = an.meanr(flatstack[:,ii,jj]/flatmedians,nsigma=sigma, \
niter=niter)
if verbose: print "computed median stack"
hdr0 = pyfits.getheader(dsubflat[0])
hdr0.update('sup-flat', 'super-flat created by IR.py')
pyfits.writeto(_sflat, superflat, header=hdr0, clobber=clobber)
if verbose: print "wrote super-flat"
return
[docs]def red_masterflat(obs, nums=[], verbose=False, sigma=3, niter=1, clobber=False):
"""Combine numbered super-flats into a single flat."""
sflats = [obs._proc + ('superflat%03g' % num) + obs._suffix for num in nums]
nflat = len(sflats)
framecenter = obs.naxis1/2, obs.naxis2/2
framesize = obs.naxis1, obs.naxis2
sflatstack = phot.subreg(sflats, framecenter, framesize, verbose=verbose)
if verbose: print "read in super-flat stack"
sflatmedians = an.amedian(sflatstack.reshape((nflat,np.prod(framesize))),1).reshape((nflat,1,1))
masterflat = np.zeros(framesize,float)
for ii in range(obs.naxis1):
for jj in range(obs.naxis2):
masterflat[ii,jj] = an.meanr(sflatstack[:,ii,jj]/sflatmedians,nsigma=sigma, \
niter=niter, verbose=verbose)
if verbose: print "computed master-flat"
hdr0 = pyfits.getheader(sflats[0])
hdr0.update('mas-flat', 'master-flat created by IR.py')
pyfits.writeto(obs.masterflat, masterflat, header=hdr0, clobber=clobber)
if verbose: print "wrote master-flat"
return
[docs]def red_sci(obs, ind, num=np.nan, verbose=False, sigma=3, niter=1):
"""Dark-subtract science frames.
I don't currently use this (Jan 2010)"""
red_sub(rawflat, dsubflat, _sdark, clobber=True)
rawsci = obs.rawall[ind]
dsubsci = obs.dsuball[ind]
nsci = len(rawsci
)
framecenter = obs.naxis1/2, obs.naxis2/2
framesize = obs.naxis1, obs.naxis2
red_sub(rawsci, dsubsci, _sdark, clobber=True)
return
[docs]def ditherskysub(obs, sci_index, clobber=False, verbose=False, flatten=True, skyramp=False):
"""
Subtract sky (and dark) levels from a set of dithered exposures
with equal exposure settings.
For now, compute a robust mean sky frame from each ndither set.
Then subtract them, flat-field them, and output.
NOTE: unlike other reduction tasks, the input index here refers
ONLY to science frames.
"""
from pyraf import iraf
rawsci = obs.rawsci[sci_index]
procsci = obs.procsci[sci_index]
nsci = len(rawsci)
ndither = obs.ndither
masterflat = pyfits.getdata(obs.masterflat)
framecenter = obs.naxis1/2, obs.naxis2/2
framesize = obs.naxis1, obs.naxis2
rawstack = phot.subreg(rawsci, framecenter, framesize)
procstack = np.zeros(rawstack.shape,float)
if verbose: print "Read in raw science stack!"
if (nsci % ndither)<>0:
print "WARNING: Input file list length (%i) not evenly"+ \
"divisible by size of dither set (%i)." \
% (nsci, ndither)
nsets = np.ceil(nsci/ndither)
if skyramp:
print "not tested yet -- and probably quite slow."
j0 = int(obs.jd[0])
for ii in range(nsets):
i0 = ii*ndither
i1 = min(ii*ndither+ndither, nsci)
thesejd = obs.jd[i0:i1]-j0
for kk in range(obs.naxis1):
for ll in range(obs.naxis2):
thisfit = an.polyfitr(thesejd, rawstack[i0:i1,kk,ll], \
1, sigma)
procstack[i0,i1,kk,ll] = rawstack[i0:i1,kk,ll] - \
polyval(thisfit, thesejd)
thissky = an.amedian(rawstack[i0:i1,:,:],0)
procstack[i0:i1,:,:] = rawstack[i0:i1,:,:]-thissky
if verbose: print "Sky-subtracted raw science data"
else:
subtractedSky = np.zeros(nsci,float)
for ii in range(nsets):
i0 = ii*ndither
i1 = min(ii*ndither+ndither, nsci)
thissky = an.amedian(rawstack[i0:i1,:,:],0)
procstack[i0:i1,:,:] = rawstack[i0:i1,:,:]-thissky
subtractedSky[i0:i1] = np.median(thissky.ravel())
if verbose: print "Sky-subtracted raw science data"
if flatten:
procstack = procstack/masterflat
if verbose: print "Flat-fielded science data"
iraf.observatory('set','lick')
for jj in range(nsci):
hdr = pyfits.getheader(rawsci[jj])
hdr.update('sky-sub', subtractedSky[jj])
hdr.update('sky-note','Sky-subtracted using %i-frame dither set' \
% ndither)
hdr.update('flatten','Flattened using %s' \
% os.path.split(obs.masterflat)[1] )
hdr.update('exptime', hdr['coadds']*hdr['itime'])
pyfits.writeto(procsci[jj], procstack[jj,:,:], hdr, clobber=clobber)
setairmassG(proscsci[jj])
if verbose: print "Wrote corrected science data to disk."
return
[docs]def calibrate(obs, verbose=False, clobber=False):
"""
Take raw NIR frames and output sky-subtracted, flat-fielded frames
ready for photometry extraction.
"""
# 2009-12-16 09:02 IJC: Begun.
# initialize and create arrays of filenames, itimes, coadds, all
# raw filenames (passed in via 'obs')
# Create indices with which to extract various observation types:
darkindex = obs.types=='dark'
flatindex = obs.types=='flat'
sciindex = obs.types=='sci'
# From darks and flats, determine a list of overlapping ITIMEs
darktimes = obs.itimes[darkindex]
flattimes = obs.itimes[flatindex]
caltimes = list(set(darktimes).intersection(flattimes))
ncal = len(caltimes)
for ii in range(ncal):
thiscaltime = caltimes[ii]
# combine darks into super-darks ('dark001', 'dark002', etc.) --
# for each dark-frame ITIME
darkindex = (obs.types=='dark') * (obs.itimes==thiscaltime)
red_dark(obs, darkindex, num=ii, verbose=verbose, clobber=clobber)
# dark-subtract flats (_proc+procflat) and combine dark-subtracted
# flats into super-flats ('flat001', 'flat002', etc.) -- for
# each flat-frame ITIME
flatindex = (obs.types=='flat') * (obs.itimes==thiscaltime)
red_flat(obs, flatindex, num=ii, verbose=verbose, clobber=clobber)
# Combine all super-flats into one master flat.
red_masterflat(obs, range(ncal), verbose=verbose, clobber=clobber)
# dark-subtract science frames (_proc+procsci) and flatten
# dark-subtracted science frames
# red_sci()
#iraf.setjd(output, date='date_obs', time='TIME_OBS',
# epoch='equinox', jd='JD', hjd='hjd')
return
[docs]def setairmassG(filename):
"""Use PyRaf tasks to set airmass for Lick/GEMINI camera observations.
Assumes PST times."""
from numpy import pi,sin,cos,arctan2,abs
from pyraf import iraf
import matplotlib.dates as dates
lat = (37. +20.6/60.)*pi/180.
long = (121. +38.2/60.)*pi/180.
epoch = 2000
#alt = 1290 # meters
hdr = pyfits.open(filename,mode='update')
# Get data and time and put them in the correct format
datestr = '%s %s PST' % (hdr[0].header['date_obs'],hdr[0].header['time_obs'])
exptime = hdr[0].header['exptime']
date = dates.num2date(dates.datestr2num(datestr)) # to UTC
UTCdate = ('%i-%02g-%02gT%02g:%02g:%02g' % (date.year, date.month,date.day,date.hour,date.minute,date.second))
# Get RA & dec in decimal degrees
ra = (phot.hms(hdr[0].header['ra']))*pi/180.
ha = (phot.hms(hdr[0].header['ha']))*pi/180. + 2*pi*exptime/86400.
dec= (phot.dms(hdr[0].header['dec']))*pi/180.
# Compute terms for coordinate conversion
term1 = sin(lat)*sin(dec)+cos(lat)*cos(dec)*cos(ha)
term2 = cos(lat)*sin(dec)-sin(lat)*cos(dec)*cos(ha)
term3 = -cos(dec)*sin(ha)
rad = abs(term3 +1j*term2)
az = arctan2(term3,term2)
alt = arctan2(term1, rad)
# Compute airmass
z = pi/2. - alt
airmass = 1./(cos(z) + 0.50572*(96.07995-z*180./pi)**(-1.6364))
hdr[0].header.update('airmass',airmass,comment='airmass at middle of observation')
hdr[0].header.update('date-obs',UTCdate,comment='UTC date/time at start of observation')
hdr[0].header.update('epoch',epoch)
hdr.flush()
hdr.close()
return
##################################################
# Helper functions for moircs calibration routines:
##################################################
[docs]def makemasterframe(list_or_array, nsigma=None):
"""If an array is passed, return it. Otherwise, return a
median stack of the input filename list."""
if hasattr(list_or_array, 'shape') and len(list_or_array.shape)>1:
masterframe = np.array(list_or_array, copy=False)
else:
if nsigma is None:
masterframe = np.median(map(pyfits.getdata, list_or_array), axis=0)
else:
masterframe = an.meanr(map(pyfits.getdata, list_or_array), nsigma=nsigma, axis=0)
return masterframe
[docs]def makeflat(flats, dark, badpix=None, samplemask=None, finalnorm=True):
"Helper function; input is list of FITS filenames and dark (to subtract)."
from pylab import find
nflat = len(flats)
shape = (pyfits.getval(flats[0], 'naxis1'), \
pyfits.getval(flats[0], 'naxis2'))
masterflat = np.zeros((nflat,) + shape, dtype=np.float32)
medvalues = np.zeros((nflat, 1, 1), dtype=np.float32)
if badpix is None:
badpix = np.zeros(shape, dtype=bool)
# Define region to sample (to appropriately normalize frames)
if samplemask is None:
samplemask = True - badpix
else:
samplemask = (True - badpix) * samplemask
these_indices = find(samplemask.ravel())
# Load individual frames and compute normalization factors:
for ii in range(nflat):
masterflat[ii] = (pyfits.getdata(flats[ii]) - dark)
medvalues[ii] = np.median(masterflat[ii].ravel()[these_indices].ravel())
masterflat[ii] /= medvalues[ii]
#pdb.set_trace()
## Redefine masterflat, from a stack to a single frame:
masterflat = np.median(masterflat, 0)
if not finalnorm:
masterflat *= medvalues.mean()
return masterflat
##################################################
# End helper functions
##################################################
[docs]def moircs_photcal(skydark, skyflat, domedark, domeflat, objdark, rawframes, procframes, badpixelmask=None, clobber=False, retmaster=False, skysub=True, twoatatime=False, verbose=False, chip=None, linearize=True):
"""Tanaka-san's prescription for inital calibration of MOIRCS images.
:INPUTS:
skydark : list of str, or 2D Numpy array
List of raw, long sky darks, or a 2D Numpy array Master Sky Dark
skyflat : list of str, or 2D Numpy array
List of raw, long sky flat frames, or a 2D Numpy array Master
Sky Flat. Exposure time should be the same as for skydark.
domedark : list of str, or 2D Numpy array
List of raw, short dome darks, or a 2D Numpy array Master Dome Dark
domeflat : list of str, or 2D Numpy array
List of raw, short dome flats, or a 2D Numpy array Master Dome
Flat. Exposure time per frame should be same as for skydark
objdark : list of str, or 2D Numpy array
List of raw, short darks, or a 2D Numpy array Master Dark.
Per-frame exposure time should be the same as for `rawframes`
rawframes : list of str
List of raw, unprocessed science frames
procframes : list of str
List of processed science frames, to be written to disk
clobber : bool
Whether to overwrite existing FITS files
retmaster : bool
Whether to return master dark and (normalized) flat frames
skysub : bool
Whether to subtract the sky (not sure why you wouldn't want to...)
twoatatime: bool
Whether to correct both chips simultaneously, in 2-frame
pairs. If true, all the file inputs above must be len-2 lists
of filename lists (or 2xN str arrays).
chip : 1 or 2
Which chip to use (to set gain and/or linearity correction)
linearize : bool
Whether or not to correct for detector nonlinearity.
:EXAMPLE:
::
import ir
obs2dome = ir.initobs('20111214_chip2', cal='dome')
obs2sky = ir.initobs('20111214_chip2', cal='sky')
obj2dark = ['%sMCSA00187%i.fits' % (obs2sky._raw, ii) for ii in range(804,833,2)]
outfns2 = [fn2.replace('proc_M', 'proc_linear_M') for fn2 in obs2sky.procsci]
masterframes2 = ir.moircs_photcal(obs2sky.rawdark, obs2sky.rawflat, \
obs2dome.rawdark, obs2dome.rawflat, obj2dark, \
obs2sky.rawsci, outfns2, chip=2, \
clobber=False, retmaster=True, verbose=True)
obs1dome = ir.initobs('20111214_chip1', cal='dome')
obs1sky = ir.initobs('20111214_chip1', cal='sky')
obj1dark = ['%sMCSA00187%i.fits' % (obs1sky._raw, ii) for ii in range(803,833,2)]
outfns1 = [fn1.replace('proc_M', 'proc_linear_M') for fn1 in obs1sky.procsci]
masterframes1 = ir.moircs_photcal(obs1sky.rawdark, obs1sky.rawflat, \
obs1dome.rawdark, obs1dome.rawflat, obj1dark, \
obs1sky.rawsci, outfns1, chip=1, \
clobber=False, retmaster=True, verbose=True)
::
import ir
obs1dome = ir.initobs('20111214_chip1', cal='dome')
obs1sky = ir.initobs('20111214_chip1', cal='sky')
obs2dome = ir.initobs('20111214_chip2', cal='dome')
obs2sky = ir.initobs('20111214_chip2', cal='sky')
obj1dark = ['%sMCSA00187%i.fits' % (obs1sky._raw, ii) for ii in range(803,833,2)]
obj2dark = ['%sMCSA00187%i.fits' % (obs2sky._raw, ii) for ii in range(804,833,2)]
outfns = [[fn1.replace('proc_M', 'proc_simul_M') for fn1 in obs1sky.procsci],
[fn2.replace('proc_M', 'proc_simul_M') for fn2 in obs2sky.procsci]]
ir.moircs_photcal([obs1sky.rawdark, obs2sky.rawdark], \
[obs1sky.rawflat, obs2sky.rawflat], \
[obs1dome.rawdark, obs2dome.rawdark], \
[obs1dome.rawflat, obs2dome.rawflat], \
[obj1dark, obj2dark], \
[obs1sky.rawsci, obs2sky.rawsci], \
outfns, \
clobber=False, retmaster=True, verbose=True, twoatatime=True)
:NOTES:
- Make master dark frames for flat and science observations.
- Subtract master flat-dark frame from each raw flat frame
- Make object masks for each raw data. [Not implemented yet!]
- median-normalize each raw flat frame.
- Median-combine these normalized raw sky data with object
masks. This is master median flat (sky).
- Make dome flat data.
- Subtract 21-sec master dark from each WASP12 raw frame.
- Flatfield them by Dome Flats.
- Apply correction for inherent nonlinearity of IR detectors.
- Measure each sky counts. Multiply the master normalized median
sky frame by the sky counts (i.e. scale the sky for each data)
then subtract that sky frame from the dome-flatfielded
data. Apply the procedure to all frames.
- Subtract the residual global pattern by the plane-fitting
(imsurfit in IRAF, for example) if necessary.
As you may already know, making object mask is a critical part of the
good median-sky data.
"""
# 2011-12-22 13:16 IJMC: Created
# 2012-01-01 21:38 IJMC: Added Skysub option
# 2012-01-02 09:59 IJMC: Added twoatatime option so that sky
# subtraction will be common between both
# chips.
# 2012-06-08 17:03 IJMC: Now apply linearity correction.
# Define helper functions:
if twoatatime:
if badpixelmask is None:
badpixelmask = [None, None]
# Load in a bad pixel mask array, if one is not provided:
if twoatatime:
for ii in range(2):
if badpixelmask[ii] is not None and \
(not hasattr(badpixelmask[ii], 'shape')):
badpixelmask[ii] = pyfits.getdata(badpixelmask[ii])
else:
if badpixelmask is not None and (not hasattr(badpixelmask, 'shape')):
badpixelmask = pyfits.getdata(badpixelmask)
# - Make 300-sec and 21-sec master dark frames.
if verbose: print "Creating master dark frames"
if twoatatime:
masterskydark = [makemasterframe(skydark[ii]) for ii in range(2)]
masterdomedark = [makemasterframe(domedark[ii]) for ii in range(2)]
masterobjdark = [makemasterframe(objdark[ii]) for ii in range(2)]
else:
masterskydark = makemasterframe(skydark)
masterdomedark = makemasterframe(domedark)
masterobjdark = makemasterframe(objdark)
# - Subtract 300sec dark frame from each 300-sec raw data for sky
# flats
# - Make object masks for each raw data.
print "For now, skipping object masking for flat frames..."
# - normalize each raw data.
# - Median-combine these normalized raw sky data with object
# masks. This is master median sky.
if verbose: print "Creating and dark-subtracting master flat frames"
if twoatatime:
masterskyflat = []
for ii in range(2):
if hasattr(skyflat[ii], 'shape') and len(skyflat[ii].shape)>1:
masterskyflat.append( np.array(skyflat[ii], copy=False) )
else:
masterskyflat.append( makeflat(skyflat[ii], masterskydark[ii], \
badpix=badpixelmask[ii]) )
else:
if hasattr(skyflat, 'shape') and len(skyflat.shape)>1:
masterskyflat = np.array(skyflat, copy=False)
else:
masterskyflat = makeflat(skyflat, masterskydark, \
badpix=badpixelmask)
# - Make dome flat data.
if twoatatime:
masterdomeflat = []
for ii in range(2):
if hasattr(domeflat, 'shape') and len(domeflat.shape)>1:
masterdomeflat.append( np.array(domeflat[ii], copy=False) )
else:
masterdomeflat.append( makeflat(domeflat[ii], masterdomedark[ii], \
badpix=badpixelmask[ii]) )
else:
if hasattr(domeflat, 'shape') and len(domeflat.shape)>1:
masterdomeflat = np.array(domeflat, copy=False)
else:
masterdomeflat = makeflat(domeflat, masterdomedark, \
badpix=badpixelmask)
if verbose: print "Calibrating and writing raw science frames to disk"
fmap = makefmap_moircs() # (for linearity correction)
# - Subtract 21-sec master dark from each WASP12 raw data.
# - Flatfield them by Dome Flats.
# - Correct for linearity!
# - Measure each sky counts. Multiply the master normalized median
# sky frame by the sky counts (i.e. scale the sky for each data)
# then subtract that sky frame from the dome-flatfielded
# data. Apply the procedure to all frames.
# - Subtract the residual global pattern by the plane-fitting
# (imsurfit in IRAF, for example) if necessary.
## 2011-12-22 14:26 IJMC: SURFFIT: TBD
if twoatatime:
for rawfile1, procfile1, rawfile2, procfile2 in \
zip(rawframes[0], procframes[0], rawframes[1], procframes[1]):
frame0_1 = pyfits.getdata(rawfile1)
frame0_2 = pyfits.getdata(rawfile2)
hdr_1 = pyfits.getheader(rawfile1)
hdr_2 = pyfits.getheader(rawfile2)
frame1_1 = (frame0_1 - masterobjdark[0]) / masterdomeflat[0]
frame1_2 = (frame0_2 - masterobjdark[1]) / masterdomeflat[1]
hdradd(hdr_1, 'COMMENT', 'Flat- and dark-corrected by MOIRCS_photcal')
hdradd(hdr_1, 'COMMENT', ' using TWOATATIME mode.')
hdradd(hdr_2, 'COMMENT', 'Flat- and dark-corrected by MOIRCS_photcal')
hdradd(hdr_2, 'COMMENT', ' using TWOATATIME mode.')
exptime = hdr_1['EXPTIME']
ncoadd = hdr_1['COADD']
nread = hdr_1['DET-NSMP']
if linearize:
frame2_1, liniter1 = linearity_correct(frame1_1, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=1), retall=True)
frame2_2, liniter2 = linearity_correct(frame1_2, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=2), retall=True)
hdradd(hdr_1, 'COMMENT', 'Linearity-corrected by MOIRCS_photcal.')
hdradd(hdr_2, 'COMMENT', 'Linearity-corrected by MOIRCS_photcal.')
hdradd(hdr_1, 'LIN_ITER', liniter1)
hdradd(hdr_2, 'LIN_ITER', liniter2)
if skysub:
skyval = np.median(np.concatenate((frame2_1, frame2_2)).ravel())
frame2_1 -= masterskyflat[0] * skyval
frame2_2 -= masterskyflat[1] * skyval
hdr_1.update('N_SKYSUB', skyval)
hdr_2.update('N_SKYSUB', skyval)
else:
hdr_1.update('N_SKYSUB', 0)
hdr_2.update('N_SKYSUB', 0)
pyfits.writeto(procfile1, frame2_1.astype(np.float32), header=hdr_1, clobber=clobber)
pyfits.writeto(procfile2, frame2_2.astype(np.float32), header=hdr_2, clobber=clobber)
else:
for rawfile, procfile in zip(rawframes, procframes):
frame0 = pyfits.getdata(rawfile)
hdr0 = pyfits.getheader(rawfile)
frame1 = (frame0 - masterobjdark) / masterdomeflat
hdradd(hdr0, 'COMMENT', 'Dark-subtracted by MOIRCS_photcal.')
hdradd(hdr0, 'COMMENT', 'Flat-fielded by MOIRCS_photcal.')
exptime = hdr0['EXPTIME']
ncoadd = hdr0['COADD']
nread = hdr0['DET-NSMP']
if linearize:
frame2, liniter = linearity_correct(frame1, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=chip), retall=True)
hdradd(hdr0, 'COMMENT', 'Linearity-corrected by MOIRCS_photcal.')
hdradd(hdr0, 'LIN_ITER', liniter)
if skysub:
skyval = np.median(frame2.ravel())
frame2 -= masterskyflat * skyval
hdradd(hdr0, 'COMMENT', 'Sky-subtracted by MOIRCS_photcal.')
hdr0.update('N_SKYSUB', skyval)
else:
hdradd(hdr0, 'COMMENT', 'No sky subtraction performed.')
hdr0.update('N_SKYSUB', 0)
pyfits.writeto(procfile, frame2.astype(np.float32), header=hdr0, clobber=clobber)
if retmaster:
ret = masterskydark, masterskyflat, masterdomedark, masterdomeflat, masterobjdark
else:
ret = None
return ret
[docs]def moircs_speccal(skydark, sky, domedark, domeflat, rawframes, procframes, badpixelmask=None, samplemask=None, clobber=False, retmaster=False, skysub=True, verbose=False, writefiles=True, gain=3.3, readnoise=40, chip=None):
"""Perform inital calibration of MOIRCS spectroscopic data.
:INPUTS:
skydark : list of str, or 2D Numpy array
List of raw, long sky darks, or a 2D Numpy array Master Sky Dark.
The "sky dark" frame is the dark frame taken with the same
integration time as the main science frames ("rawframes") and
the sky frames ("skyflat").
sky : list of str, or 2D Numpy array
List of raw sky frames, or a 2D Numpy array Master Sky
frame. Exposure time should be the same as for "skydark" and
"rawframes".
domedark : list of str, or 2D Numpy array
List of raw, short dome darks, or a 2D Numpy array Master Dome
Dark. In this (spectroscopic) case, the "domedark" files are
the lamp-off spectroscopic files used to measure the thermal
background underlying the main dome files ("domeflat").
domeflat : list of str, or 2D Numpy array
List of raw flat frames (taken with the dome flat lamp ON), or
a 2D Numpy array Master Dome Flat. Exposure time per frame
should be same as for "domedark".
rawframes : list of str
List of raw, unprocessed science frames. Integration times
should be the same as "skydark" and "sky."
procframes : list of str
List of processed science frames, to be written to disk.
badpixelmask : str, or 2D Numpy array
Bad pixel mask: equal to zero at good pixels, and to 1 at bad pixels.
samplemask : 2D Numpy array
Sampling mask tells the sky-scaling algorithms which parts of
the image to use. Equal to 1 at points to be sampled, and
zero elsewhere. If set to None (default), a crude algorithm
will attempt to find these regions based on the "domedark"
frames.
clobber : bool
Whether to overwrite existing FITS files
retmaster : bool
Whether to return master dark and (normalized) flat frames
skysub : bool
Whether to subtract the sky (e.g., you might not want to if
sky conditions changed by a lot, and the subtraction is no
longer valid).
writefiles : bool
Whether to write files (i.e., just to test things out)
chip : 1 or 2 ONLY
Which detector is being reduced (critical for linearity correction)
:EXAMPLE:
::
import ir
odome = ir.initobs('20120427_chip2', cal='dome')
osky = ir.initobs('20120427_chip2', cal='sky')
masterframes = ir.moircs_speccal(osky.rawdark, osky.rawflat, \
odome.rawdark, odome.rawflat, \
odome.rawsci, odome.procsci, \
badpixelmask=odome.badpixelmask, \
clobber=False, retmaster=True, \
verbose=True, skysub=False, writefiles=False, chip=2)
:DESCRIPTION:
- Make master 'science dark' and bad pixel mask from individual
"skydark" frames.
- Make super-dome flat frames (two: lamp on, and lamp off)
- Subtract LAMPOFF from LAMPON frame, creating master,
unnormalized dome (MUD) flat.
- Normalize each spectrum in the MUD flat along the dispersion
direction (e.g., using a low-order polynomial), creating the
master dome flat.
- Make master sky frame by median-combining the raw sky frames
and subtracting the master 'science dark' frame.
- For each raw science frame: subtract 'science dark' frame;
scale and subtract the master sky frame; divide by the master
dome flat.
- Apply linearity correction.
:NOTES:
Unfortunately, MOIRCS reduction is best done as a two-step
process. First, set up the observations in :func:`initobs`
using some standard, generic bad pixel mask (perhaps produced by
:func:`makebadpixelmask_darks`); reduce the data using this
routine. From the reduced, flattened data, generate a second
(higher-fidelity) bad pixel map (using
:func:`makebadpixelmask_science`). Combine these two bad pixel
masks (and perhaps, as well, the appropriate mask from the
MOIRCS instrument website). Insert the proper reference to this
in :func:`initobs`, and then re-reduce the data.
"""
# 2012-04-28 01:14 IJMC: Created from moircs_photcal
# 2012-06-08 21:40 IJMC: Now does linearity corrections.
# 2012-07-06 17:06 IJMC: Added check for which detector user is reducing.
# 2012-08-10 00:04 IJMC: Now save subregion keywords into FITS headers.
from scipy import signal
import spec
import tools
#from pyraf import iraf
twoatatime = False
nsigma = 5
# Define helper functions:
if twoatatime:
if badpixelmask is None:
badpixelmask = [None, None]
def findbrightregions(data, medfilt=9, nbin=200):
mdmax = np.sort(data.ravel())[[int(0.99*data.size)]]
valbins = np.linspace(-mdmax, mdmax, nbin+1)
y, x0 = np.histogram(data.ravel(), valbins)
x = .5 * (x0[1::] + x0[0:-1])
gfit, err_gfit = spec.fitGaussian(y, guess=[y.max()*5, 3, nbin/2, np.median(y)], verbose=verbose)
cutoff_val = np.interp(gfit[2]+gfit[1]*5, range(nbin), x)
return signal.medfilt2d((data > cutoff_val).astype(np.float32), medfilt)
# Load in a bad pixel mask array, if one is not provided:
if twoatatime:
for ii in range(2):
if badpixelmask[ii] is not None and \
(not hasattr(badpixelmask[ii], 'shape')):
badpixelmask[ii] = pyfits.getdata(badpixelmask[ii])
else:
if badpixelmask is not None and (not hasattr(badpixelmask, 'shape')):
badpixelmask = pyfits.getdata(badpixelmask)
fmap = makefmap_moircs() # (for linearity correction)
# - Make master 'science dark' from individual "skydark"
# frames.
#
# - Make super-dome flat frames (two: lamp on, and lamp off)
if verbose: print "Creating master dark frames"
if twoatatime:
masterskydark = [makemasterframe(skydark[ii]) for ii in range(2)]
masterdomedark = [makemasterframe(domedark[ii]) for ii in range(2)]
else:
masterskydark = makemasterframe(skydark)
masterdomedark = makemasterframe(domedark)
# If one has not been passed in, define a sampling region (to
# normalize spectral frames); determine the boundaries of these
# regions (to save in the header).
if samplemask is None:
if twoatatime:
samplemask = [findbrightregions(masterdomedark[ii], medfilt=9) for ii in range(2)]
samplemask_corners = [tools.findRectangles(samplemask[ii], minsepy=50, minsepx=200) for ii in range(2)]
nspec = [s.shape[1] for s in samplemask_corners]
else:
samplemask = findbrightregions(masterdomedark, medfilt=9)
samplemask_corners = tools.findRectangles(samplemask, minsepy=50, minsepx=200)
nspec = samplemask_corners.shape[0]
# Never sample bad pixels:
if twoatatime:
samplemask = [samplemask[ii] * (True - badpixelmask[ii])]
else:
samplemask = samplemask * (True - badpixelmask)
# - Make master sky frame by median-combining the raw sky frames
# and subtracting the master 'science dark' frame.
if twoatatime:
mastersky = []
for ii in range(2):
if hasattr(sky, 'shape') and len(sky.shape)>1:
mastersky.append( np.array(sky[ii], copy=False) )
else:
mastersky.append( makeflat(sky[ii], masterskydark[ii], \
badpix=badpixelmask[ii], \
samplemask=samplemask[ii]) )
else:
if hasattr(sky, 'shape') and len(sky.shape)>1:
mastersky = np.array(sky, copy=False)
else:
mastersky = makeflat(sky, masterskydark, \
badpix=badpixelmask, samplemask=samplemask)
# - Subtract LAMPOFF from LAMPON frame, creating master,
# unnormalized dome (MUD) flat.
if verbose: print "Creating and dark-subtracting master dome flat"
if twoatatime:
mudflat = []
for ii in range(2):
if hasattr(domeflat[ii], 'shape') and len(domeflat[ii].shape)>1:
mudflat.append( np.array(domeflat[ii], copy=False) )
else:
mudflat.append( makeflat(domeflat[ii], masterdomedark[ii], \
badpix=badpixelmask[ii], \
samplemask=samplemask[ii], finalnorm=False) )
else:
if hasattr(domeflat, 'shape') and len(domeflat.shape)>1:
mudflat = np.array(domeflat, copy=False)
else:
mudflat = makeflat(domeflat, masterdomedark, \
badpix=badpixelmask, samplemask=samplemask, finalnorm=False)
# - Normalize each spectrum in the MUD flat along the
# dispersion direction (e.g., using a low-order
# polynomial), creating the master dome flat. PyRAF could
# do this, but it crashes -- so, do it on our own!
masterdomeflat = spec.normalizeSpecFlat(mudflat*gain, nspec=nspec, minsep=100, badpixelmask=badpixelmask)
# - For each raw science frame: subtract 'science dark' frame;
# scale and subtract the master sky frame; divide by the master
# dome flat.
if verbose: print "Calibrating and writing raw science frames to disk"
if twoatatime:
for rawfile1, procfile1, rawfile2, procfile2 in \
zip(rawframes[0], procframes[0], rawframes[1], procframes[1]):
frame0_1 = pyfits.getdata(rawfile1)
frame0_2 = pyfits.getdata(rawfile2)
hdr_1 = pyfits.getheader(rawfile1)
hdr_2 = pyfits.getheader(rawfile2)
frame1_1 = (frame0_1 - masterskydark[0])
frame1_2 = (frame0_2 - masterskydark[1])
frame1_1 /= masterdomeflat[0]
frame1_2 /= masterdomeflat[1]
hdradd(hdr_1, 'COMMENT', 'Flat- and dark-corrected by MOIRCS_speccal')
hdradd(hdr_1, 'COMMENT', ' using TWOATATIME mode.')
hdradd(hdr_2, 'COMMENT', 'Flat- and dark-corrected by MOIRCS_speccal')
hdradd(hdr_2, 'COMMENT', ' using TWOATATIME mode.')
for jj in range(samplemask_corners[0].shape[0]):
hdradd(hdr_1, 'SUBREG%i' % jj, str(samplemask_corners[0][jj]))
for jj in range(samplemask_corners[1].shape[0]):
hdradd(hdr_2, 'SUBREG%i' % jj, str(samplemask_corners[1][jj]))
hdradd(hdr_1, 'NSUBREG', samplemask_corners[0].shape[0])
hdradd(hdr_2, 'NSUBREG', samplemask_corners[1].shape[0])
exptime = hdr_1['EXPTIME']
ncoadd = hdr_1['COADD']
nread = hdr_1['DET-NSMP']
frame2_1, liniter1 = linearity_correct(frame1_1, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=1), retall=True)
frame2_2, liniter2 = linearity_correct(frame1_2, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=2), retall=True)
hdradd(hdr_1, 'COMMENT', 'Linearity-corrected by MOIRCS_speccal.')
hdradd(hdr_1, 'LIN_ITER', liniter1)
hdradd(hdr_2, 'COMMENT', 'Linearity-corrected by MOIRCS_speccal.')
hdradd(hdr_2, 'LIN_ITER', liniter2)
if skysub:
sample_indices = [samplemask[ii].nonzero() for ii in range(2)]
skyval = np.median(np.concatenate((frame2_1[sample_indices[0]], frame2_2[sample_indices[1]])))
frame2_1 -= mastersky[0] * skyval
frame2_2 -= mastersky[1] * skyval
hdradd(hdr_1, 'COMMENT', 'Sky-subtracted by MOIRCS_speccal.')
hdr_1.update('N_SKYSUB', skyval)
hdradd(hdr_2, 'COMMENT', 'Sky-subtracted by MOIRCS_speccal.')
hdr_2.update('N_SKYSUB', skyval)
else:
hdr_1.update('N_SKYSUB', 0)
hdradd(hdr_1, 'COMMENT', 'No scaled-sky subtraction.')
hdr_2.update('N_SKYSUB', 0)
hdradd(hdr_2, 'COMMENT', 'No scaled-sky subtraction.')
if writefiles:
pyfits.writeto(procfile1, frame2_1.astype(np.float32), header=hdr_1, clobber=clobber)
pyfits.writeto(procfile2, frame2_2.astype(np.float32), header=hdr_2, clobber=clobber)
else:
for rawfile, procfile in zip(rawframes, procframes):
frame0 = pyfits.getdata(rawfile)
hdr0 = pyfits.getheader(rawfile)
frame1 = (frame0 - masterskydark)
hdradd(hdr0, 'COMMENT', 'Dark-subtracted by MOIRCS_speccal')
frame1 /= masterdomeflat
hdradd(hdr0, 'COMMENT', 'Flat-corrected by MOIRCS_speccal')
for jj in range(samplemask_corners.shape[0]):
hdradd(hdr0, 'SUBREG%i' % jj, str(samplemask_corners[jj]))
hdradd(hdr0, 'NSUBREG', samplemask_corners.shape[0])
exptime = hdr0['EXPTIME']
ncoadd = hdr0['COADD']
nread = hdr0['DET-NSMP']
if chip is None:
try:
chip = hdr0['DET-ID']
except:
chip = None
if verbose:
print "File %s lacks header keyword DET-ID. Setting chip to None." % rawfile
frame2, liniter = linearity_correct(frame1, nread, ncoadd, 11.8, exptime, fmap, linearity_moircs, verbose=verbose, linfuncargs=dict(chip=chip), retall=True)
hdradd(hdr0, 'COMMENT', 'Linearity-corrected by MOIRCS_speccal.')
hdradd(hdr0, 'LIN_ITER', liniter)
if skysub:
sample_indices = samplemask.nonzero()
skyval = np.median(frame2[sample_indices])
frame2 -= mastersky * skyval
hdradd(hdr0, 'COMMENT', 'Sky-subtracted by MOIRCS_speccal')
hdr0.update('N_SKYSUB', skyval)
else:
hdr0.update('N_SKYSUB', 0)
if writefiles:
pyfits.writeto(procfile, frame2.astype(np.float32), header=hdr0, clobber=clobber)
if retmaster:
ret = masterskydark, mastersky, masterdomedark, masterdomeflat, samplemask
else:
ret = None
if not writefiles:
print "'writefiles' was False, so no files were actually written to disk."
return ret
[docs]def spectral_subregion_mask(flatframedata, minsep=42, maxboundarywidth=15):
"""Find boundaries of each slit's spectrum, using a flat-field frame:
Helper function for :func:`mosfire_speccal` and :func:`moircs_speccal`
"""
# 2013-01-23 08:16 IJMC: Transferred to standalone function
import tools
vertical_gradient = np.diff(flatframedata, axis=0)
avg_gradient = vertical_gradient.mean(1)
gradient_dispersion = an.dumbconf(avg_gradient, .683)[0]
pos_thresh = np.median(avg_gradient) + 5*gradient_dispersion
neg_thresh = np.median(avg_gradient) - 5*gradient_dispersion
# Estimate the spectral boundaries:
peak_up0 = tools.find_peaks(avg_gradient, sep=minsep, thresh=pos_thresh)
peak_down0 = tools.find_peaks(-avg_gradient, sep=minsep, thresh=-neg_thresh)
# Boundaries without a partner are spurious:
peak_up, peak_down = [], []
for pup in peak_up0[1]:
if (np.abs(peak_down0[1] - pup) < maxboundarywidth).any():
peak_up.append(pup)
peak_down.append(peak_down0[1][(np.abs(peak_down0[1] - pup) < maxboundarywidth).nonzero()][0])
boundaries = np.vstack((peak_up, peak_down)).mean(0)
boundaries = np.concatenate(([0], np.sort(boundaries), [peak_down0[1].max()]))
# Define a mask:
mask = np.zeros(flatframedata.shape, dtype=bool)
for ii in range(boundaries.size-1):
ind1 = max(0, int(boundaries[ii] + np.floor(maxboundarywidth/2.)))
ind2 = min(2048, int(boundaries[ii+1] - np.floor(maxboundarywidth/2.)))
mask[ind1:ind2] = True
return mask, boundaries
[docs]def mosfire_speccal(sky, domedark, domeflat, rawframes, procframes, badpixelmask=None, samplemask=None, clobber=False, retmaster=False, skysub=False, verbose=False, writefiles=True, gain=2.15, readnoise=21, linearize=False, flatten=True, polycoef=None, slitjawflatcorrect=True, dodark=True):
"""Perform inital calibration of MOSFIRE spectroscopic data.
:INPUTS:
rawsky : list of str, or 2D Numpy array
List of raw sky frames, or a 2D Numpy array Master Sky
frame. Exposure time should, ideally, be the same as for
"rawframes". [Deprecated -- no longer necessary!]
domedark : list of str, or 2D Numpy array
List of raw, short dome darks, or a 2D Numpy array Master Dome
Dark. In this (spectroscopic) case, the "domedark" files are
the lamp-off spectroscopic files used to measure the thermal
background underlying the main dome files ("domeflat").
domeflat : list of str, or 2D Numpy array
List of raw flat frames (taken with the dome flat lamp ON), or
a 2D Numpy array Master Dome Flat. Exposure time per frame
should be same as for "domedark".
rawframes : list of str
List of raw, unprocessed science frames. Integration times
should be the same as "skydark" and "sky."
procframes : list of str
List of processed science frames, to be written to disk.
badpixelmask : str, or 2D Numpy array
Bad pixel mask: equal to zero at good pixels, and to 1 at bad
pixels. Note that bad pixels are not interpolated-over in this
stage of the reduction; that doesn't happen until spectra
extraction.
samplemask : 2D Numpy array
Sampling mask tells the sky-scaling algorithms which parts of
the image to use. Equal to 1 at points to be sampled, and
zero elsewhere. If set to None (default), a crude algorithm
will attempt to find these regions based on the "domedark"
frames.
clobber : bool
Whether to overwrite existing FITS files
retmaster : bool
Whether to return master dark and (normalized) flat frames
skysub : bool
Whether to subtract the sky (e.g., you might not want to if
sky conditions changed by a lot, and the subtraction is no
longer valid).
writefiles : bool
Whether to write files (i.e., just to test things out)
linearize : bool
Whether to correct for detector nonlinearity.
flatten : bool
Whether to use flat-frames for correction
polycoef : str, None, or sequence
Polynomial coefficients for detector linearization: see
:func:`ir.linearity_correct` and :func:`ir.linearity_mosfire`.
slitjawflatcorrect : bool
Whether to try to correct the dome flats for the uneven slit
width (from the slit jaw mechanisms).
dodark: bool
If True, apply correction for dark frames in 'domedark.'
If False, don't! (i.e., Y or J band observations)
:EXAMPLE:
::
import ir
obs = ir.initobs('20121010')
masterframes = ir.mosfire_speccal(obs.rawsky[0:10], \
obs.rawdark[::2], obs.rawflat[::2], \
obs.rawsci, obs.procsci, \
badpixelmask=obs.badpixelmask, \
clobber=False, retmaster=True, \
verbose=True, skysub=False, \
linearize=True, polycoef='/Users/ianc/proj/transit/data/mosfire_linearity/linearity/mosfire_linearity_cnl_coefficients_bic-optimized.fits')
mastersky, masterdomedark, masterdomeflat, mudflat, samplemask = \
masterframes
:DESCRIPTION:
- Make super-dome flat frames (two: lamp on, and lamp off)
- Subtract LAMPOFF from LAMPON frame, creating master,
unnormalized dome (MUD) flat.
- Normalize each spectrum in the MUD flat along the dispersion
direction (e.g., using a low-order polynomial), creating the
master dome flat.
- Make master sky frame by median-combining the raw sky frames
- For each raw science frame: scale and subtract the master sky
frame (if desired); divide by the master dome flat.
- Apply linearity correction. (if desired)
:NOTES:
Unfortunately, MOSFIRE reduction is best done as a two-step
process. First, set up the observations in :func:`initobs`
using some standard, generic bad pixel mask (perhaps produced by
:func:`makebadpixelmask_darks`); reduce the data using this
routine. From the reduced, flattened data, generate a second
(higher-fidelity) bad pixel map (using
:func:`makebadpixelmask_science`). Combine these two bad pixel
masks (and perhaps, as well, the appropriate mask from the
MOIRCS instrument website). Insert the proper reference to this
in :func:`initobs`, and then re-reduce the data.
"""
# 2012-10-21 20:24 IJMC: Created from moircs_speccal
# 2012-11-05 13:36 IJMC: Added reference pixel computation.
# 2013-04-08 10:55 IJMC: Added 'slitjawflatcorrect' option.
# 2013-10-02 19:49 IJMC: Added 'dodark' option.
from scipy import signal
import spec
import tools
#from pyraf import iraf
nsigma = 5
# Load in a bad pixel mask array, if one is not provided:
if badpixelmask is not None and (not hasattr(badpixelmask, 'shape')):
badpixelmask = pyfits.getdata(badpixelmask)
if linearize:
fmap = makefmap_mosfire() # (for linearity correction)
# - Make super-dome flat frames (two: lamp on, and lamp off)
if dodark:
if verbose: print "Creating master dark frames"
masterdomedark = makemasterframe(domedark)
else:
if verbose: print "No dark frames subtracted."
masterdomedark = 0.
# If one has not been passed in, define a sampling region (to
# normalize spectral frames); determine the boundaries of these
# regions (to save in the header).
if samplemask is None:
if verbose: print "Creating sampling mask"
masterdomeflat0 = makemasterframe(domeflat)
samplemask, boundaries = spectral_subregion_mask(masterdomeflat0)
samplemask *= (masterdomeflat0 > 500)
samplemask_corners = tools.findRectangles(samplemask, minsepy=42, minsepx=200)
del masterdomeflat0
nspec = samplemask_corners.shape[0]
# Never sample bad pixels:
if badpixelmask is None:
badpixelmask = np.zeros(samplemask.shape, dtype=bool)
samplemask = samplemask * (True - badpixelmask)
# - Make master sky frame by median-combining the raw sky frames
# and subtracting the master 'science dark' frame.
if skysub:
if verbose: print "Creating master sky frame."
if hasattr(sky, 'shape') and sky.ndim==2:
mastersky = np.array(sky, copy=False)
else:
mastersky = makeflat(sky, 0., \
badpix=badpixelmask, samplemask=samplemask)
else:
if verbose: print "skysub=False, so not creating master sky frame."
mastersky = None
# - Subtract LAMPOFF from LAMPON frame, creating master,
# unnormalized dome (MUD) flat.
if verbose: print "Creating and dark-subtracting master dome flat"
if hasattr(domeflat, 'shape') and domeflat.ndim==2:
mudflat = np.array(domeflat, copy=False)
else:
mudflat = makeflat(domeflat, masterdomedark, \
badpix=badpixelmask, samplemask=samplemask, finalnorm=False)
# - Normalize each spectrum in the MUD flat along the
# dispersion direction (e.g., using a low-order
# polynomial), creating the master dome flat. PyRAF could
# do this, but it crashes -- so, do it on our own!
#pdb.set_trace()
masterdomeflat = spec.normalizeSpecFlat(samplemask*mudflat*gain, nspec=nspec, minsep=42, badpixelmask=badpixelmask + (mudflat < 500))
if slitjawflatcorrect:
masterdomeflat /= np.median(masterdomeflat, 1).reshape(2048,1)
# - For each raw science frame: scale and subtract the master
# sky frame (if desired); divide by the master dome flat.
if verbose: print "Calibrating and writing raw science frames to disk"
for rawfile, procfile in zip(rawframes, procframes):
hdr0 = pyfits.getheader(rawfile)
# Compute reference pixel values, if possible:
if ('sampmode' in hdr0 and hdr0['sampmode']==2) or \
('smpmode0' in hdr0 and hdr0['smpmode0']==2):
refpix_values, frame1 = mosfire_refpixvalues(rawfile, retframe=True, mode='cds')
for refpix, refpixval in enumerate(refpix_values):
hdradd(hdr0, 'REFPIX00', refpixval)
else:
print "SAMPMODE is not CDS: cannot compute reference pixel levels."
frame1 = pyfits.getdata(rawfile)
if flatten:
frame1 /= masterdomeflat
hdradd(hdr0, 'COMMENT', 'Flat-corrected by MOSFIRE_speccal')
if slitjawflatcorrect:
hdradd(hdr0, 'COMMENT', 'Option -slitjawflatcorrect- used for master dome flat.')
else:
hdradd(hdr0, 'COMMENT', 'NOT flat-corrected.')
for jj in range(nspec):
hdradd(hdr0, 'SUBREG%i' % jj, str(samplemask_corners[jj]))
hdradd(hdr0, 'NSUBREG', nspec)
exptime = hdr0['TRUITIME']
ncoadd = hdr0['COADDONE']
nread = hdr0['READDONE']
if linearize:
# If you edit this section, be sure to update the same
# section in spec.calibrate_stared_mosfire_spectra!
frame2, liniter = linearity_correct(frame1, nread-1, ncoadd, 1.45, exptime, fmap, linearity_mosfire, verbose=verbose, linfuncargs=dict(polycoef=polycoef), retall=True)
hdradd(hdr0, 'COMMENT', 'Linearity-corrected by MOSFIRE_speccal.')
hdradd(hdr0, 'LIN_ITER', liniter)
if isinstance(polycoef, str):
hdradd(hdr0, 'LINPOLY', os.path.split(polycoef)[1])
else:
if isinstance(polycoef, np.ndarray) and polycoef.ndim==3:
hdradd(hdr0, 'LINPOLY', 'user-input array used for linearity correction.')
else:
hdradd(hdr0, 'LINPOLY', str(polycoef))
else:
frame2 = frame1
hdradd(hdr0, 'COMMENT', 'NOT linearity-corrected.')
if skysub:
sample_indices = samplemask.nonzero()
skyval = np.median(frame2[sample_indices])
frame2 -= mastersky * skyval
hdradd(hdr0, 'COMMENT', 'Sky-subtracted by MOSFIRE_speccal')
hdr0.update('N_SKYSUB', skyval)
else:
hdr0.update('N_SKYSUB', 0)
if writefiles:
pyfits.writeto(procfile, frame2.astype(np.float32), header=hdr0, clobber=clobber)
if retmaster:
ret = mastersky, masterdomedark, masterdomeflat, mudflat, samplemask
else:
ret = None
if not writefiles:
print "'writefiles' was False, so no files were actually written to disk."
return ret
[docs]def hdradd(header, key, value):
"""Add keyword to FITS header, without overwriting existing keys.
:INPUTS:
header : FITS header object (from pyfits.getheader)
header to be updated
key : str
Keyword to use. If it exists, an integer will be appended to
form a new, non-duplicate keyword.
value : (anything)
Value to assign to keyword in header.
"""
# 2012-06-08 15:02 IJMC: Created
# 2012-11-05 13:46 IJMC: Added 'key0' to propertly keep track of initial keyword name.
iter = 0
if not key in header:
header.update(key, value)
else:
key0 = '' + key
while key in header:
iter += 1
striter = str(iter)
key = key0[0:(8-len(striter))] + striter
header.update(key, value)
return
def makebadpixelmask_dark(darkfiles, homedir=homedir):
import spec
nx, ny = pyfits.getval(darkfiles[0], 'NAXIS1'), pyfits.getval(darkfiles[0], 'NAXIS2')
bigstack = phot.subreg2(darkfiles)
masterdark = np.median(bigstack, axis=0)
pixelerr = an.stdr(bigstack, axis=0, nsigma=3)
del bigstack
pyfits.writeto(homedir+'/temp/mdark.fits', masterdark, clobber=True)
pyfits.writeto(homedir+'/temp/perr.fits', pixelerr, clobber=True)
nbin = 200
meanerr = np.median(pixelerr)
bb = np.linspace(-10*meanerr, 10*meanerr, nbin+1)
ncts, nbns0 = np.histogram(pixelerr.ravel(), bb)
xcts = .5 * (nbns0[1::] + nbns0[0:-1])
err_ncts = np.max((np.sqrt(ncts), np.ones(nbin)), axis=0)
gfit, egfit = spec.fitGaussian(ncts, err=err_ncts)
ecutoff = np.interp(gfit[2]+5*gfit[1], range(nbin), xcts)
bpm1 = pixelerr > ecutoff # Very variable pixels
bpm2 = pixelerr == 0 # Stuck pixels
badpixelmask = bpm1 + bpm2
return badpixelmask
[docs]def makebadpixelmask_science(filenamelist, readnoise, gain=1, nsigma=5, frac=0.95, medfilt=[5,1]):
"""Make bad pixel mask from science frames.
filenamelist : list of str
list of filenames, in ADU (i.e., DN).
These files should already have been flat-fielded!
readnoise : scalar
read noise, in electrons
gain : scalar
gain, in electrons / ADU
medfilt : odd int Width of 2D median filter (keep it low, or sky
lines start getting flagged!). A 2-vector (e.g., [5, 1]) can
also be set, to to avoid marking sky lines.
"""
# 2012-07-24 22:26 IJMC: Created
# 2013-10-02 21:49 IJMC: Added 'medfilt' option.
# 2014-05-22 11:22 IJMC: Added 2D medfilt option.
from scipy.signal import medfilt2d
nx, ny = pyfits.getval(filenamelist[0], 'NAXIS1'), pyfits.getval(filenamelist[0], 'NAXIS2')
nim = len(filenamelist)
outlierstack = np.zeros((nx, ny), dtype=int)
for filename in filenamelist:
frame = pyfits.getdata(filename) * gain
smooth_frame = medfilt2d(frame, kernel_size=medfilt)
err_smoothframe = np.sqrt(smooth_frame + readnoise**2)
outliers = np.abs((frame - smooth_frame) / err_smoothframe) > nsigma
outlierstack += outliers
badpixelmask = outlierstack >= (nim * frac)
return badpixelmask
[docs]def makefmap_moircs(tread=11.2e-6, treset=8e-3, dt=11.8, npix=1024*1024):
"""Make f-map for linearity correction of Subaru/MOIRCS data.
:INPUTS:
tread : scalar
time between reads of individual pixels
treset : scalar
time required to reset an individual pixel
dt : scalar
time between successive reads of the detector array
npix : int
Total number of pixels in each readout section of the detector.
:OUTPUTS:
fmap : 2D Numpy array
The variable 'f' defined in Equation 10 or 11 of the reference
given below.
:REFERENCES:
Vacca et al. 2004, PASP.
Private communications w/Dr. I. Tanaka, NAOJ/Subaru
"""
# 2012-06-05 14:32 IJMC: Created
# 2012-07-13 12:03 IJMC: Updated with new readout prescription per
# I. Tanaka.
## The following generate the old, wrong F_map -- it is kept here
## only for reference, and in case the situation as currently
## understood changes once again:
#quad = np.arange(npix).reshape(np.sqrt(npix), np.sqrt(npix))
#twoquad = np.hstack((quad, np.rot90(quad, 3)))
#m = np.vstack((twoquad, np.rot90(twoquad, 2)))
quad = np.arange(npix).reshape(np.sqrt(npix), np.sqrt(npix))[:,::-1]
twoquad = np.hstack((np.rot90(quad, 3), quad))
m = np.vstack((np.rot90(twoquad, 2), twoquad))
fglobal = 1. - (m * tread + treset) / dt
return fglobal
[docs]def makefmap_mosfire(tread=1e-5, treset=8e-3, dt=1.45, npix=64*2048):
"""Make f-map for linearity correction of Keck/MOSFIRE data.
:INPUTS:
tread : scalar
time between reads of individual pixels
treset : scalar
time required to reset an individual pixel
dt : scalar
time between successive reads of the detector array
npix : int
Total number of pixels in each readout section of the detector.
:OUTPUTS:
fmap : 2D Numpy array
The variable 'f' defined in Equation 10 or 11 of the reference
given below.
:REFERENCES:
Vacca et al. 2004, PASP.
K. Kulas, private communication, July 2012
"""
# 2012-12-28 10:29 IJMC: Created
a = np.arange(64*2048).reshape(2048, 64)
b = np.hstack((a, a[:,::-1]))
m = np.rot90(np.tile(b, (1, 16)), 1)
fglobal = 1. - (m * tread + treset) / dt
return fglobal
[docs]def linearity_mosfire(measured_counts, polycoef=None, unityminimum=True, scale=1e4):
"""Linearity correction curve for Keck/MOSFIRE detectors.
:INPUTS:
measured_counts : scalar or sequence
Measured digital counts (ADU or DN, _not_ photoelectrons)
polycoef : str, sequence, or None
sequence:
Sequence of polynomial coefficients [a_n, ... a_2, a_1, 1.],
suitable for passing to numpy.polyval, to use in the
implementation of Vacca et al.'s Eq. 20. If None, use my
defaults (which are empirically calibrated only up to
30000 ADU)
str:
Filename of 3-D array c (of size [n+1]x2048x2048) of
polynomial coefficients for each pixel, defined such that
c[:,i,j] = [a_n, ... a_2, a_1, 1.] for pixel (i,j).
polycoef='/Users/ianc/proj/transit/data/mosfire_linearity/linearity/mosfire_linearity_cnl_coefficients_bic-optimized.fits'
unityminimum : bool
If True, values less than unity (1.0) are never returned.
scale : scalar
Scale used to divide measured counts in polynomial evaluation.
This helps to avoid extremely small coefficients.
:NOTES:
This function should *not* be used to perform linearity
corrections on MOSFIRE data; for that, you should use
:func:`linearity_correct` (for which this function is only a
helper).
:SEE_ALSO:
:func:`makefmap_mosfire`, :func:`linearity_correct` :func:`linearity_moircs`
:REFERENCE:
Vacca et al. 2004, PASP.
"""
# 2012-12-28 21:12 IJMC: Created (from MOIRCS file)
# 2014-05-22 10:00 IJMC: fixed bug in call to np.array
measured_counts = np.array(measured_counts, copy=False)
if measured_counts.ndim==0:
measured_counts = measured_counts.reshape(1)
if polycoef is None:
polycoef = np.array([-0.00122734, 0.00873443, \
-0.02259769, 0.02420229, \
-0.00845038, 1.0])
if isinstance(polycoef, np.ndarray) and polycoef.ndim==1:
correction_factor = 1./np.polyval(polycoef, measured_counts/scale)
elif isinstance(polycoef, str) or (isinstance(polycoef, np.ndarray) and polycoef.ndim==3):
if isinstance(polycoef, str):
polycoef = pyfits.getdata(polycoef)
vecs = (measured_counts/scale).reshape(2048*2048)**np.arange(polycoef.shape[0]-1, -0.5, -1).reshape(polycoef.shape[0],1)
correction_factor = (vecs * polycoef.reshape(polycoef.shape[0], 2048*2048)).sum(0).reshape(2048, 2048)
# correction_factor = 0. * measured_counts
# for ii in range(measured_counts.shape[0]):
# for jj in range(measured_counts.shape[1]):
# correction_factor[ii,jj] = 1./np.polyval(polycoef[:,ii,jj], measured_counts[ii,jj]/scale)
else:
print "Unknown or invalid polycoef; exiting without linearity correction."
return -1
if unityminimum:
correction_factor[correction_factor < 1] = 1.
return correction_factor
[docs]def linearity_correct(measurement, nread, ncoadd, dt, exptime, fmap, linfunc, npix=1024*1024, tol=1e-6, maxiter=10, verbose=False, linfuncargs=None, retall=False):
"""Correct a set of IR array measurements for nonlinearity.
:INPUTS:
measurement : 1D or 2D sequence or array
Values reported back by the infrared array detector (in ADU or
DN, not photoelectrons)
nread : int
Number of reads used for measurement (for CDS, set nread = 1)
ncoadd : int
Number of coadds used for measurement
dt : scalar
Amount of time between successive reads of the detector
exptime : scalar
Integration time for measurement
fmap : sequence or array (same shape as 'measurement')
'f' values (cf. :func:`makefmap_moircs`)
linfunc : function
Function which, when passed a measured value, returns the
linearized value (the C_nl function of Vacca et al.;
cf. :func:`linearity_moircs`)
npix : int
Total number of pixels in each readout section of the detector.
tol : scalar
When all measurements have been linearized to within this
tolerance, routine will exit.
maxiter : int
If 'tol' is not reached within this many iterations,
linearization routine will exit anyway.
linfuncargs : dict
Dict of arguments, if any, to be passed to 'linfunc'.
retall : bool
If True, also return the number of iterations.
:OUTPUTS:
A Numpy array of the same shape as 'measurement,' but
appropriately linearized following the algorithm of the
reference below.
If retall is True, the linearized measurement is the first
element of a 2-tuple.
:NOTES:
Vacca et al. recommend that this algorithm be applied _after_
dark current has been subtracted and flat-fielding performed.
:SEE_ALSO:
:func:`makefmap_moircs`, :func:`linearity_moircs`
:REFERENCE:
Vacca et al. 2004, PASP.
:EXAMPLE:
::
raw_mosfire_data = 2e4 * np.ones((2048, 2048), dtype=float)
coef_fn = '/Users/ianc/proj/transit/data/mosfire_linearity/linearity/mosfire_linearity_cnl_coefficients_bic-optimized_corners-fixed.fits'
fmap = ir.makefmap_mosfire()
corrected_data = ir.linearity_correct(dat, 1, 1, 100, 101, fmap, ir.linearity_mosfire, npix=2048*2048, verbose=True, linfuncargs=dict(polycoef=coef_fn))
"""
# 2012-06-05 14:35 IJMC: Created
meas_init = np.array(measurement, copy=True)
if linfuncargs is None:
linfuncargs = dict()
corrected_measurement = meas_init.copy()
meas_old = corrected_measurement * 0.5
iter = 0
while (np.abs(corrected_measurement/meas_old - 1) > tol).any() and iter < maxiter:
p_one = (corrected_measurement / (nread * ncoadd)) * (dt / exptime) * (0.5*(nread + 1) - fmap)
s_one = meas_init / (nread * ncoadd) + p_one # _original_ value of stot
s_two = s_one * linfunc(s_one, **linfuncargs)
p_two = p_one * linfunc(p_one, **linfuncargs)
meas_old = corrected_measurement.copy()
corrected_measurement = (ncoadd * nread) * (s_two - p_two)
iter += 1
if verbose:
if iter==maxiter:
print "Reached maximum number of iterations (%i); stopping." % iter
else:
print "%i iterations to linearize within tolerance." % iter
if retall:
ret = corrected_measurement, iter
else:
ret = corrected_measurement
return ret
[docs]def moircs_spectra(fns, gain, readnoise, varfns=None, idlexec=None, dx=(200, -100), edgeboundary=25, specwidth=60, sep=100, findthresh=200, subregions='SUBREG', IDLoptions = "adjoptions={center:1,centerfit:1,centerdeg:3}, bgdeg=3, adjfunc='adjgauss'", pyrafoptions=None, tempframefn='tempframe.fits', tempbpmfn='tempbpm.fits', tempspecfn = 'tempspec.fits', tempscriptfn='temp_specextract.pro', writefns=None, findlocs='pick', badpixelmask=None, mode='pyraf', superoptions=None, trace_options=None, multiordlocs=None, verbose=False):
"""Extract spectra from a calibrated MOIRCS multi-object frame.
:INPUTS:
fns : str, or list of strs
Filenames of MOIRCS files from which to extract spectra. No,
this cannot be the MOIRCS frames as Numpy Arrays.
gain : scalar
Detector gain, in electrons / ADU
readnoise : scalar
Detector read noise, in electrons
idlexec : str
Path to the IDL executable. OPTSPECEXTR.PRO and its
associated files must be in your IDL path. If set to None,
then it will be set to: os.popen('which idl').read().strip()
:OPTIONS:
dx : 2-tuple
Extend the subregions as defined by dx[0] to the left, and by
dx[0] to the right. Something of a fudge factor (sorry!).
edgeboundary : int
Trim both sides of each subregion by this many pixels *in the
spatial direction*
varfns : str, or list of strs
Optional names of files with variances of data in 'fns'.
specwidth : int
Approximate width of spectral profile, in pixels (for IDL mode).
sep : scalar
Minimum separation of spectral traces within a particular
subregion. For finding multiple spectral traces within a given
subregion.
findthresh : scalar
Ignore all values lower than this value (cf. :func:`tools.find_peaks`)
subregions : str, or 2D numpy array
if str:
The primary header keyword from which the subregion
boundaries will be read (e.g., 'SUBREG0', 'SUBREG1', etc;
cf. :func:`moircs_speccal`)
if array:
Each row of the array corresponds to the subregion
boundaries.
In either case, the format of each subregion boundary B is
such that the subregion is:
subreg = frame[B[2]:B[3], B[0]:B[1]]
IDLoptions : str
Options to pass to OPTSPECEXTR.PRO. For example:
"adjfunc='adjgauss', adjoptions={center:1,centerfit:1,centerdeg:3}, bgdeg=3"
Note that this Python code will break if you _don't_ trace
the spectrum (adjoptions, etc.); this is an area for future
work if I ever use a spectrograph with straight traces.
tempframefn : str
If input 'frame' is an array, it will be written to this
filename in order to pass it to IDL.
tempbpmfn : str
Python will write the subregion bad pixel mask to this
filename in order to pass it to IDL
tempspecfn : str
IDL will write the spectral data to this filename in order to
pass it back to Python.
tempscriptfn : str
Filename in which a short IDL script will be written.
writefns : str or list of str
If None, write nothing.
returnspec : bool
If True, return outputs.
findlocs : str ('auto' or 'pick')
Find spectral traces automatically, or pick them manually.
badpixelmask : None or str
Filename of a 2D bad pixel mask, equal to 1 for bad pixels and
0 for good pixels.
superoptions : dict or None
options for mode 'superextract'
trace_options : dict or None
options for spectral profile tracing (see :func:`spec.traceorders`)
multiordlocs : None or 2D array (or list of lists)
My hackneyed way to automatically pass multiple 'ordlocs'
calls to :func:`spec.traceorders` (if function is invoked via
trace_options). Successively use each of the ordlocs values
in multiordlocs.
mode : str
'superextract': see :func:`spec.superExtract`
'idl' : see :func:`optspecextr_idl`
'pyraf' : see :func:`optspecextr_pyraf`
:OUTPUT:
For each frame, an object will be returned with various
interesting and useful fields. The 'specs' field will contain a
list (one per subregion), for which each element is the output
of :func:`spec.optspecextr_idl`.
:TO_DO:
Allow the routine to write files to disk, and update FITS
headers (writefns, returnspec).
Do better at finding the spectral trace.
:REQUIREMENTS:
:doc:`tools`
:SEE_ALSO:
:func:`moircs_speccal`, :func:`spec.optspecextr_idl`
:EXAMPLE:
::
import ir
obs = ir.initobs('20120429_wasp3_chip2', cal='dome')
filenames = [fn.replace('proc_', 'proc_fft-defringe_') for fn in obs.procsci]
out = ir.moircs_spectra(filenames[0], obs.gain, obs.readnoise, None, findlocs='pick', mode='superextract')
all = ir.moircs_spectra(filenames[0:2], obs.gain, obs.readnoise, None, findlocs=out.locs, subregions=out.subregions[out.spec_subregions])
::
import ir
obs = ir.initobs('20120429_wasp3_chip2', cal='dome')
filenames = [fn.replace('proc_', 'proc_fft-defringe_') for fn in obs.procsci]
speclocs = ir.moircs_spectra(filenames[0], obs.gain, obs.readnoise, findlocs='pick', mode='pickonly')
trace_options = dict(pord=3, fitwidth=60, dispaxis=0)
# Extract spectra once, just to trace spectra:
superoptions = dict(bord=2, bkg_radii=[30, 70], extract_radius=18,
polyspacing=1., pord=1, verbose=True,
qmode='fast-linear', dispaxis=0, csigma=5,
retall=True)
trace_out = ir.moircs_spectra(filenames[0], obs.gain, obs.readnoise, None, mode='superextract', superoptions=superoptions, trace_options=trace_options, findlocs=speclocs.locs, subregions=speclocs.subregions[speclocs.spec_subregions])
# Extract spectra once, just to trace spectra:
meanx = np.array(tools.flatten([range(ss.trace.size) for ss in trace_out.specs])).mean()
meany = np.array(tools.flatten([ss.trace for ss in trace_out.specs])).mean()
trace_options['ordlocs'] = [meanx, meany]
# Extract multiple spectra without any user interaction:
superoptions.update(dict(pord=2, polyspacing=0.25))
out2 = ir.moircs_spectra(filenames[0:2], obs.gain, obs.readnoise, None, mode='superextract', superoptions=superoptions, trace_options=trace_options, findlocs=speclocs.locs, subregions=speclocs.subregions[speclocs.spec_subregions])
"""
# 2012-08-18 17:48 IJMC: Created
# 2012-09-05 14:08 IJMC: Big edits to handle superExtract and tracing.
# 2012-09-10 09:21 IJMC: Now save outputs, semi-sanely.
# 2012-09-17 21:10 IJMC: Option for input of variance filenames.
# 2012-11-06 08:23 IJMC: Added verbose option
# 2012-12-13 14:09 IJMC: Now save subregion indices used for extraction.
# 2013-03-31 19:31 IJMC: Now each 'tempframe' and 'tempspec' gets
# its own unique timestamp... this seems to
# avoid random, untraceable, unrepeatable
# errors and failures.
import pylab as py
import clickinput as ci
from tools import find_peaks
import spec
from nsdata import bfixpix
import inspect
# Save inputs for output:
save_options = dict()
argnames = inspect.getargspec(moircs_spectra)[0]
for argname in argnames:
exec("save_options['%s'] = %s" % (argname, argname))
# Parse inputs:
if isinstance(fns, str):
fnwasstring = True
fns = [fns]
else:
fnwasstring = False
if isinstance(varfns, str) or varfns is None:
varfns = [varfns] * len(fns)
if idlexec is None:
idlexec = os.popen('which idl').read().strip()
if badpixelmask is None:
badpixelmask = np.zeros((2048, 2048), dtype=bool)
else:
badpixelmask = pyfits.getdata(badpixelmask).astype(bool)
# Loop through all files:
out = []
for i, fn in enumerate(fns):
if verbose: print 'Beginning file %i: %s' % (i, fn)
thisout = baseObject()
thisout.fn = fn
thisout.specs = [] # Spectral Data
frame = pyfits.getdata(fn)
hdr = pyfits.getheader(fn)
frame_specs = []
if varfns[i] is None:
variance = np.abs(frame)/gain * (readnoise/gain)**2
else:
variance = pyfits.getdata(varfns[i])
# Read in the subregion locations of the spectra.
if isinstance(subregions, str):
corners = []
moresubregions = True
nsub = 0
while moresubregions:
if ('%s%i' % (subregions,nsub)) in hdr:
corners.append(map(int, hdr['%s%i' % (subregions,nsub)].replace('[', '').replace(']', '').split()))
nsub += 1
else:
moresubregions = False
else:
corners = subregions
nsub = len(corners)
corners = np.array(corners, copy=False)
thisout.subregions = corners
thisout.locs = []
thisout.spec_subregions = []
if isinstance(trace_options, dict):
trace_options_original = dict().update(trace_options)
else:
trace_options_original = trace_options
# Loop through each subregion of the frame
if findlocs=='pick':
py.figure()
nspec = 0
for cc in range(nsub):
if verbose: print "Beginning subregion %i (of %i): %s" % (cc+1, nsub, str(corners[cc]))
data = frame[corners[cc,2]:corners[cc,3], max(0, corners[cc,0]-dx[0]):corners[cc,1]+dx[1]]
data_variance = variance[corners[cc,2]:corners[cc,3], max(0, corners[cc,0]-dx[0]):corners[cc,1]+dx[1]]
data_bpm = badpixelmask[corners[cc,2]:corners[cc,3], max(0, corners[cc,0]-dx[0]):corners[cc,1]+dx[1]]
if findlocs=='auto':
junk, locs = find_peaks((np.median(data, axis=1) - np.median(np.median(data, axis=1)))[edgeboundary:-edgeboundary], sep=sep, thresh=findthresh)
locs += edgeboundary
elif findlocs=='pick': # Pick spectral trace manually:
py.clf()
py.plot(np.median(data, axis=1) - np.median(np.median(data, axis=1)))
speclimit = data.shape[0]
axlim = py.axis()
py.xlim([0, speclimit*1.2])
py.plot([speclimit]*2, py.ylim(), '-k', linewidth=4)
locs = []
keep_picking = True
while keep_picking:
pick = ci.pickloc(resetaxes=True, titstr='Click to select a trace.\nClick to the right of the vertical black line to end.', zoom=[speclimit/2., (axlim[3]-axlim[2])/2.])
if hasattr(pick, '__iter__'):
lastloc, ypos = pick
else:
lastloc = np.inf
if lastloc < speclimit:
locs.append(int(lastloc))
py.plot([lastloc]*2, py.ylim(), '--r')
py.draw()
else:
keep_picking = False
else: # We passed in the locations, directly!
locs = [findlocs[cc]]
thisout.locs = thisout.locs + locs
thisout.spec_subregions = thisout.spec_subregions + [cc]*len(locs)
#pdb.set_trace()
# Loop through each spectrum of the subregion:
lociter = 0
for loc in locs:
lociter += 1
if verbose: print "Beginning spectrum %i of %i (in subregion)" % (lociter, len(locs))
## Replaced this (the old, confusing way of doing it):
#lowerlimit = max(loc-sep, edgeboundary)
#upperlimit = min(loc+sep,data.shape[0]-edgeboundary)
#sub = data[lowerlimit:upperlimit, :]
#subvar = data_variance[lowerlimit:upperlimit, :]
#sub_bpm = data_bpm[lowerlimit:upperlimit, :]
lh_index = max(0, corners[cc,0]-dx[0])
rh_index = corners[cc,1]+dx[1]
top_index = corners[cc,2] + max(loc - sep, edgeboundary)
bot_index = corners[cc,2] + min(loc + sep, corners[cc,3] - corners[cc,2]-edgeboundary)
sub = frame[top_index:bot_index, lh_index:rh_index]
subvar = variance[top_index:bot_index, lh_index:rh_index]
sub_bpm = badpixelmask[top_index:bot_index, lh_index:rh_index]
if mode=='superextract':
superoptions.update(goodpixelmask=True-sub_bpm.astype(bool))
trace_options.update(badpixelmask=sub_bpm.astype(bool))
if multiordlocs is not None:
trace_options.update(ordlocs=multiordlocs[nspec])
print "Automatically using multiordlocs: ", multiordlocs[nspec]
spx = spec.spextractor(sub, gain, readnoise, mode='superextract', options=superoptions, trace_options=trace_options, framevar=subvar)
junk = superoptions.pop('goodpixelmask')
junk = trace_options.pop('badpixelmask')
elif mode=='idl': # Run IDL, and transfer the data back to Python
x0 = loc - (top_index - corners[cc,2])
x1, x2 = x0 - specwidth/2., x0 + specwidth/2.
pyfits.writeto(tempframefn, sub.transpose(), clobber=True)
pyfits.writeto(tempbpmfn, 1. - sub_bpm.transpose().astype(np.int16), clobber=True)
spx = spec.optspecextr_idl(tempframefn, gain, readnoise, x1, x2, idlexec, specfn=tempspecfn, IDLoptions=IDLoptions, inmask=tempbpmfn)
elif mode=='pyraf':
sub = bfixpix(sub, sub_bpm, n=8, retdat=True)
pyfits.writeto(tempbpmfn, sub_bpm.transpose().astype(np.int16), clobber=True)
from time import time
tic = time()
tempframefn2 = tempframefn.replace('.fits', '%i.fits' % (time()*1000))
tempspecfn2 = tempspecfn.replace('.fits', '%i.fits' % (time()*1000))
pyfits.writeto(tempframefn2, sub.transpose(), clobber=True)
spx = optspecextr_pyraf(tempframefn2, tempspecfn2, gain, readnoise, badpixelmask=tempbpmfn, options=pyrafoptions, verbose=verbose)
elif mode=='pickonly':
spx = -1
else:
spx = 'No valid spectral extraction mode specified'
try:
spx.extracted_subregion = np.array([lh_index, rh_index, top_index, bot_index])
except:
pass
frame_specs.append(spx)
nspec += 1
#if len(thisout.spec_subregions)>0:
# thisout.spec_subregions = np.array([thisout.spec_subregions])
thisout.specs = frame_specs
#specs.append(frame_specs)
thisout.specs_fields = ['Spectrum [ADU]', 'Uncertainty in Spectrum [ADU]', 'Spectral Trace Location [pix]', 'Spectral Trace Width']
thisout.options = save_options
out.append(thisout)
# If we only had a single frame:
if fnwasstring:
out = out[0]
return out
[docs]def optspecextr_pyraf(filename, specname, gain, readnoise, badpixelmask=None, options=None, verbose=False, cleanup=True):
"""Python script to call Pyraf and extract spectra from a file.
:INPUTS:
filename : str
filename of frame (measured in ADU) from which to extract spectrum.
specname : str
name of Pyraf-extracted spectrum.
gain : scalar
Gain of detector (electrons per ADU). If frame is already in
units of electrons, set this to 1.0. If gain is a string, we
will try to read the eponymous FITS header keyword.
readnoise : scalar
Read noise of detector, in electrons. If gain is a string, we
will try to read the eponymous FITS header keyword.
badpixelmask : str
filename of bad pixel mask; equal to 1 at bad pixels, equal to
0 at good pixels.
options : dict
A dict of all the various options to be passed to iraf's
'apall' command.
cleanup : bool
Whether to delete working files after processing
:RETURNS:
An object with spectral bands as defined in the IRAF output file.
"""
# 2012-08-25 14:56 IJMC: Created
# 2012-09-24 12:07 IJMC: Now use iraf.fixpix instead of
# bfixpix. Also forcibly remove all files
# created during the process to ensure
# reproducability.
# 2012-11-06 09:28 IJMC: Now set IRAF directory to current Python working directory
# 2012-11-20 08:25 IJMC: Added 'cleanup' option.
# 2013-02-14 14:31 IJMC: Now explicitly load iraf packages
from pyraf import iraf
from nsdata import bfixpix
try:
from astropy.io import fits as pyfits
except:
import pyfits
import os
if badpixelmask is not None:
iraf.fixpix(filename, badpixelmask, cinterp=1)
#badpixelmask = pyfits.getdata(badpixelmask)
#bfixpix(frame, badpixelmask, n=8)
frame = pyfits.getdata(filename)
hdr = pyfits.getheader(filename)
if isinstance(gain, str):
gain = hdr[gain]
if isinstance(readnoise, str):
readnoise = hdr[readnoise]
default_options = dict()
background_options = dict(b_function='chebyshev', b_order=3, b_sample='-25:-35,25:35', bkg='yes', b_naverage=-3, b_niterate=2, background='fit')
find_options = dict(ylevel='INDEF', llimit='INDEF', ulimit=25, lower=-25, upper=25, recenter='yes',resize='no',extras='yes', nfind=1, nsubaps=1, minsep=100, find='yes', nsum=-100, line=int(frame.shape[0]*0.7))
further_options = dict(format='multispec', weights='variance', t_order=5, t_niterate=3, t_naverage=3, clean='yes', interactive=True, pfit='fit2d', review='no')
default_options.update(background_options)
default_options.update(find_options)
default_options.update(further_options)
if options is None:
options = dict()
options['output'] = specname
options['gain'] = gain
options['readnoise'] = readnoise
default_options.update(options)
if verbose:
print "Called with options:"
print default_options
iraf.unlearn('apall')
iraf.noao()
iraf.imred()
iraf.specred()
iraf.cd(os.getcwd())
if verbose: print "Changing IRAF directory to working directory, ", os.getcwd()
if os.path.isfile(specname):
iraf.delete(specname)
iraf.apall(filename, **default_options)
data = pyfits.getdata(specname)
hdrout = pyfits.getheader(specname)
# Clean up after myself:
if cleanup:
if os.path.isfile(specname):
if verbose: print "Deleting %s ..." % specname
# iraf.delete(specname)
os.remove(specname)
if os.path.isfile(filename):
if verbose: print "Deleting %s ..." % filename
# iraf.delete(filename)
os.remove(filename)
if os.path.isfile(badpixelmask):
if verbose: print "Deleting %s ..." % badpixelmask
# iraf.delete(badpixelmask)
os.remove(badpixelmask)
if os.path.isfile('database/ap%s' % filename.replace('.fits', '')):
if verbose: print "Deleting %s ..." % ('database/ap%s' % filename.replace('.fits', ''))
# iraf.delete('database/ap%s' % filename.replace('.fits', ''))
os.remove('database/ap%s' % filename.replace('.fits', ''))
if os.path.isfile('database/aplast'):
if verbose: print "Deleting database/aplast ..."
os.remove('database/aplast')
# Prepare data for delivery:
data = data.squeeze()
nband = data.shape[0]
spx = baseObject()
band_strings = [hdrout['bandid%i' % (ii+1)] for ii in range(nband)]
band_names = [bs.split('-')[0].strip() for bs in band_strings]
band_notes = [bs.split('-')[1].strip() for bs in band_strings]
spx.header = hdr
spx.header_iraf = hdrout
for ii in range(nband):
setattr(spx, band_names[ii], data[ii])
setattr(spx, band_names[ii]+'_notes', band_notes[ii])
return spx
[docs]def moircs_defringe(procfn, defringefns, subreg='SUBREG', clobber=False, plotonly=False, dx=None, dy=None, spatial_index='select', retall=False, pord=1, corners=None, defringe_options=dict()):
"""Defringe multiple subregions in MOIRCS multi-object spectroscopy frames.
subregions : str, or 2D numpy array
if str:
The primary header keyword from which the subregion
boundaries will be read (e.g., 'SUBREG0', 'SUBREG1', etc;
cf. :func:`moircs_speccal`)
defringe_options -- passed along to :func:`spec.defringe_sinusoid`
"""
# Make the filenames into iterable lists:
import tools
import pylab as py
import clickinput as ci
import spec
if isinstance(procfn, str) or not hasattr(procfn, '__iter__'):
procfn = [procfn]
if isinstance(defringefns, str) or not hasattr(defringefns, '__iter__'):
defringefns = [defringefns]
# Read in the subregion locations of the spectra.
if corners is None:
hdr0 = pyfits.getheader(procfn[0])
corners = []
nsub = 0
moresubregions = True
while moresubregions:
if ('subreg%i' % nsub) in hdr0:
corners.append(map(int, hdr0['%s%i' % (subreg,nsub)].replace('[', '').replace(']', '').split()))
nsub += 1
else:
moresubregions = False
corners = np.array(corners)
nsub = corners.shape[0]
# Select spatial indices:
if spatial_index=='select':
spatial_index = []
for cc in range(nsub):
data, junk, newcorners = tools.extractSubregion(procfn[0], corners=corners[cc], dx=dx, dy=dy, retall=True)
py.figure()
py.plot(np.median(data, axis=1) - np.median(np.median(data, axis=1)))
speclimit = data.shape[0]
axlim = py.axis()
py.xlim([-speclimit*0.2, speclimit*1.2])
py.plot([0]*2, py.ylim(), '--r', linewidth=4)
py.plot([speclimit]*2, py.ylim(), '-k', linewidth=4)
locs = []
keep_picking = True
while keep_picking:
pick = ci.pickloc(resetaxes=True, titstr='Click (twice, slowly) to select the left and right edges of the exclusion region.\nClick to the right of the vertical black line to end. Click to LHS of red dashed to reset selections.', zoom=[speclimit/2., (axlim[3]-axlim[2])/2.])
if hasattr(pick, '__iter__'):
lastloc, ypos = pick
else:
lastloc = -100
if lastloc < 0:
locs = []
py.clf()
py.plot(np.median(data, axis=1) - np.median(np.median(data, axis=1)))
speclimit = data.shape[0]
axlim = py.axis()
py.xlim([-speclimit*0.1, speclimit*1.2])
py.plot([0]*2, py.ylim(), '--r', linewidth=4)
py.plot([speclimit]*2, py.ylim(), '-k', linewidth=4)
elif lastloc < speclimit:
locs.append(int(lastloc))
py.plot([lastloc]*2, py.ylim(), '--r')
py.draw()
else:
keep_picking = False
index = np.ones(speclimit, dtype=bool)
if len(locs)>=2:
index[int(locs[0]):int(locs[1])] = False
spatial_index.append(index)
else:
pass
#spatial_index = np.array(spatial_index, copy=False)
#if spatial_index.ndim==1:
# spatial_index = spatial_index.reshape(1, spatial_index.size)
# 2012-09-19 17:28 IJMC: Created.
for fn, fn2 in zip(procfn, defringefns):
frame = pyfits.getdata(fn)
hdr = pyfits.getheader(fn)
# Loop through each subregion of the frame
for cc in range(nsub):
if plotonly:
pass
else:
data, junk, newcorners = tools.extractSubregion(fn, corners=corners[cc], dx=dx, dy=dy, retall=True)
#frame[corners[cc,2]:corners[cc,3], max(0, corners[cc,0]+dx):corners[cc,1]]
# Detrend the data to enhance the visibility of the fringes.
# Remove a quadratic from each column to account for mild spectral
# tilts. Any fringing trends that are linear will not be
# de-fringes, but in these cases the linear-fit background removal
# will account for the gradient:
trends_data = np.zeros(data.shape, dtype=float)
nx0, ny0 = data.shape
x00 = np.arange(nx0) / (nx0/2.) - 1.
for ii in range(ny0):
linfit = an.polyfitr(x00[spatial_index[cc]], data[:,ii][spatial_index[cc]], pord, 3)
trends_data[:,ii] = np.polyval(linfit, x00)
det_data0 = data - trends_data
newfringes = spec.defringe_sinusoid(det_data0, spatial_index=spatial_index[cc], **defringe_options)
frame[newcorners[2]:newcorners[3], newcorners[0]:newcorners[1]] -= newfringes
#pdb.set_trace()
if not plotonly:
hdr.update('DEFRINGE', 'Defringed w/spec.defringe_sinusoid')
pyfits.writeto(fn2, frame, hdr, clobber=clobber)
if retall:
ret = spatial_index
else:
ret = None
return ret
[docs]def mosfire_refpixvalues(fn, **kw):
"""Compute and return MOSFIRE reference pixel values, from individual readouts.
:INPUTS:
fn : str
Filename of science frame (e.g., 'm121010_0487.fits').
:OPTIONS:
mode : str
e.g., 'CDS' [ONLY CDS MODE IS VALID, FOR NOW]
retframe : bool
If true, return the tuple (refpix_values, pyfits.getdata(fn))
:OUTPUTS:
a length-32 Numpy array.
:NOTES:
For a filename 'frame.fits', reads in the file as well as the
individual CDS reads "frame_001.fits" and "frame_002.fits".
Computes the difference between (frame - (frame002-frame001)),
and uses it to return the reference pixel bias levels.
"""
# 2012-11-05 13:22 IJMC: Created
# 2013-04-08 10:33 IJMC: Removed "has_key" calls; deprecated method.
try:
from astropy.io import fits as pyfits
except:
import pyfits
if 'mode' in kw:
mode = str(kw['mode']).lower()
else:
mode = 'cds'
if 'retframe' in kw:
retframe = kw['retframe']
else:
retframe = False
data = pyfits.getdata(fn)
if mode=='cds':
fn1 = fn.replace('.fits', '_001.fits')
fn2 = fn.replace('.fits', '_002.fits')
difference_of_reads = pyfits.getdata(fn2) - pyfits.getdata(fn1)
refpix_levels = (difference_of_reads.transpose()[::-1] - data).mean(1)[::64]
else:
print "Only CDS mode is implemented!"
refpix_levels = -1
if retframe:
ret = refpix_levels, data
else:
ret = refpix_levels
return ret