Skip to content

GitLab

Commit 09c81308 authored by Monica Rainer's avatar Monica Rainer
Browse files
Options

Add keywords to FITS header

parent 1a11858c
Hide whitespace changes
Inline Side-by-side
Showing with 190 additions and 60 deletions
GOFIO_manual.pdf
View file @ 09c81308
No preview for this file type
gofio/__version__.py
View file @ 09c81308
__version__ = '1.6.0'
__version__ = '1.6.1'
gofio/drslib/config.py
View file @ 09c81308
......@@ -171,6 +171,15 @@ CONFIG['RED'] = (' '.join((CONFIG['KEY_DRS'],'RED')),'GOFIO reduced file')
CONFIG['RED_TYPE'] = (' '.join((CONFIG['KEY_DRS'],'RED TYPE')),'GOFIO reduced file type')
CONFIG['RED_SUBTYPE'] = (' '.join((CONFIG['KEY_DRS'],'RED SUBTYPE')),'GOFIO reduced file subtype')
CONFIG['RED_SLIT'] = (' '.join((CONFIG['KEY_DRS'],'RED SLIT')),'GOFIO reduced file slit')
CONFIG['YPOS'] = (' '.join((CONFIG['KEY_DRS'],'SIGNAL_POS')),'Pixel position of the signal maximum on the slit') # stare images
CONFIG['YPOS_LOW'] = (' '.join((CONFIG['KEY_DRS'],'SIGNAL_LOW')),'Lower pixel position of the signal on the slit') # stare images
CONFIG['YPOS_UP'] = (' '.join((CONFIG['KEY_DRS'],'SIGNAL_UP')),'Upper pixel position of the signal on the slit') # stare images
CONFIG['AYPOS'] = (' '.join((CONFIG['KEY_DRS'],'ASIGNAL_POS')),'Pixel position of the A signal maximum on the slit') # nodding images
CONFIG['AYPOS_LOW'] = (' '.join((CONFIG['KEY_DRS'],'ASIGNAL_LOW')),'Lower pixel position of the A signal on the slit') # nodding images
CONFIG['AYPOS_UP'] = (' '.join((CONFIG['KEY_DRS'],'ASIGNAL_UP')),'Upper pixel position of the A signal on the slit') # nodding images
CONFIG['BYPOS'] = (' '.join((CONFIG['KEY_DRS'],'BSIGNAL_POS')),'Pixel position of the B signal maximum on the slit') # nodding images
CONFIG['BYPOS_LOW'] = (' '.join((CONFIG['KEY_DRS'],'BSIGNAL_LOW')),'Lower pixel position of the B signal on the slit') # nodding images
CONFIG['BYPOS_UP'] = (' '.join((CONFIG['KEY_DRS'],'BSIGNAL_UP')),'Upper pixel position of the B signal on the slit') # nodding images
# Wavelength calibration
......@@ -240,7 +249,7 @@ CONFIG['MASK_C'] = os.path.join(CONFIG['RES_DIR'], 'GIANOB_MASKC.fits')
CONFIG['A_POS'] = 8. # guess pixel position for A
CONFIG['B_POS'] = 23. # guess pixel position for B
CONFIG['C_POS'] = 15. # guess pixel position for C
CONFIG['Y_POS'] = 4 # max. pixel difference between expected and real pixel position
CONFIG['Y_POS'] = 7 # max. pixel difference between expected and real pixel position
CONFIG['HWTM'] = 6. # hwtm limit for contamination between A and B
CONFIG['S1D'] = True # flag for creation of s1d output (True/False)
CONFIG['S1D_NORM'] = False # flag for normalization of s1d output (True/False)
......
gofio/drslib/gofio2ascii.py
View file @ 09c81308
"""
Written by Monica Rainer
2019-05-22
gofio2ascii.py converts the GOFIO ms1d/s1d outputs in ASCII files with or without header
The ms1d file have and 4 column: wavelength (nm), flux, S/N, and number of the echelle order
It is written in python2, but if needed it may be run with python3 by
commenting and un-commenting the appropriate rows:
"lista = raw_input/input ..."
"for o in xrange/range(50):"
The header options are:
- none: no header is written
- simple (default): only two rows with the name of the FITS files and the columns definition
- red: reduced header, with the keywords defined in the "keywords" list
If the user wants different keywords to be saved in the header,
the keywords list may be changed accordingly.
- full: the whole header is written in the ASCII files
Usage:
spa_gofio2ascii.py -h
spa_gofio2ascii.py
spa_gofio2ascii.py [--hea=<value> --split]
spa_gofio2ascii.py <list>
spa_gofio2ascii.py <list> [--hea=<value> --split]
Options:
-h,--help : show this screen
list : list of ms1d FITS data files [Default: fits.lis]
--hea=<value> : write header, options: none/simple/red/full [Default: simple] (defined by the "keywords" list)
--split : split the orders of the ms1d files in different ASCII outputs
"""
from __future__ import (absolute_import, division, print_function, unicode_literals)
from astropy.io import fits,ascii
import numpy as np
lista = raw_input('List of s1d/ms1d fits files to convert in ascii (default = fits.lis):\n')
if lista == '':
lista='fits.lis'
split = raw_input('Do you want to split the orders of the ms1d files (if any) ([y/n] default yes)? ')
if split == '':
split = 'y'
lis = open(lista,'r')
from docopt import docopt
import os
docopt_args = docopt(__doc__)
keywords = ['HIERARCH TNG OBS PROG ID' , 'OBJECT' , 'HIERARCH TNG TEL TARG ALPHA' , \
'HIERARCH TNG TEL TARG DELTA' , 'HIERARCH TNG TEL TARG EQUINOX' , 'EXPSTART' , 'MJD-OBS' , \
'AIRMASS' , 'HIERARCH TNG DRS AIRMASS' , 'DIT' , 'HIERARCH TNG TPL NPAIRS' , \
'HIERARCH TNG EXP NGROUPS' , 'NCOMBINE' , 'HIERARCH TNG DRS GAIN' , 'HIERARCH TNG DRS RON' , \
'HIERARCH TNG INS GRAY FILTER' , 'HIERARCH TNG INS DEWAR TEMPDET' , 'HIERARCH TNG METEO HUMIDITY' , \
'HIERARCH TNG DRS EXPTIME' , 'HIERARCH TNG DRS BJD' , 'HIERARCH TNG DRS BERV']
wave_fmt = '%10.3f'
flux_fmt = '%15.8f'
snr_fmt = '%10.3f'
order_fmt = '%i'
if not docopt_args['<list>']:
lista = 'fits.lis'
if not os.path.isfile(lista):
lista = raw_input('List of ms1d/s1d FITS files to convert in ASCII (default = fits.lis):\n') #python2
#lista = input('List of ms1d FITS files to convert in ASCII (default = fits.lis):\n') #python3
else:
lista = docopt_args['<list>']
#print(docopt_args['--hea'])
with open(lista) as framelist:
for frame in framelist:
frame = frame.strip()
calib_ascii = frame.replace('.fits','.txt')
calib_ascii = frame.replace('.fits','.txt')
spettro = fits.open(frame)
try:
......@@ -24,42 +72,85 @@ with open(lista) as framelist:
wcalib = sdata.field(1)
spec = sdata.field(2)
snr = sdata.field(3)
for o in xrange(50):
order = np.zeros(wcalib.shape)
hea1 = '\n'.join(( os.path.basename(frame) ,'1. wavelength (nm) - 2. flux - 3. S/N - 4. echelle order\n'))
hea2 = spettro[0].header.tostring(sep='\n')
if docopt_args['--hea'] == 'full':
hea = ''.join((hea1,hea2))
elif docopt_args['--hea'] == 'red':
heasel = []
hea3 = hea2.split('\n')
for row in hea3:
for key in keywords:
if key in row:
heasel.append(row)
break
heasel = '\n'.join(heasel)
hea = ''.join((hea1,heasel))
else:
hea = hea1
for o in xrange(50): #python2
#for o in range(50): #python3
wcalib[o] = wcalib[o][::-1]
spec[o] = spec[o][::-1]
snr[o] = snr[o][::-1]
if split == 'y':
order[o] = o+32
if docopt_args['--split']:
txt = '_' + str(o+32) + '.txt'
order_ascii = frame.replace('.fits',txt)
order = np.vstack((wcalib[o],spec[o],snr[o]))
ascii.write(np.transpose(order),order_ascii, format='no_header')
print 'Saved %s' % (order_ascii)
if split == 'n':
optCal = np.vstack((np.concatenate(np.flipud(wcalib)),np.concatenate(np.flipud(spec)),np.concatenate(np.flipud(snr))))
ascii.write(np.transpose(optCal),calib_ascii, format='no_header')
print 'Saved %s' % (calib_ascii)
order = np.vstack((wcalib[o],spec[o],snr[o],order[o]))
if docopt_args['--hea'] == 'none':
ascii.write(np.transpose(order),order_ascii, fmt=[wave_fmt, flux_fmt , snr_fmt , order_fmt])
else:
ascii.write(np.transpose(order),order_ascii, fmt=[wave_fmt, flux_fmt , snr_fmt , order_fmt], header=hea)
print('Saved %s' % (order_ascii))
if not docopt_args['--split']:
output = np.vstack( ( np.concatenate(np.flipud(wcalib)), np.concatenate(np.flipud(spec)), np.concatenate(np.flipud(snr)) , np.concatenate(np.flipud(order)) ) )
#ascii.write(np.transpose(optCal),calib_ascii, format='no_header')
if docopt_args['--hea'] == 'none':
np.savetxt(calib_ascii,np.transpose(output), fmt=[wave_fmt, flux_fmt , snr_fmt , order_fmt])
else:
np.savetxt(calib_ascii,np.transpose(output), fmt=[wave_fmt, flux_fmt , snr_fmt , order_fmt], header=hea)
print('Saved %s' % (calib_ascii))
except:
hea = spettro[0].header
head = spettro[0].header
data = spettro[0].data
start = hea['CRVAL1']
step = hea['CDELT1']
new = open(calib_ascii,'w')
k = 0
while k < len(data):
wave = start + (k*step)
new.write(str(wave) + '\t' + str(data[k]) + '\n')
k = k+1
new.close()
print 'Saved %s' % (calib_ascii)
hea1 = '\n'.join(( os.path.basename(frame) ,'1. wavelength (nm) - 2. flux\n'))
hea2 = spettro[0].header.tostring(sep='\n')
if docopt_args['--hea'] == 'full':
hea = ''.join((hea1,hea2))
elif docopt_args['--hea'] == 'red':
heasel = []
hea3 = hea2.split('\n')
for row in hea3:
for key in keywords:
if key in row:
heasel.append(row)
break
heasel = '\n'.join(heasel)
hea = ''.join((hea1,heasel))
else:
hea = hea1
start = head['CRVAL1']
step = head['CDELT1']
length = head['NAXIS1']
end = start + (step*length) - step
waves = np.linspace(start, end, length)
output = np.vstack( (waves, data ) )
if docopt_args['--hea'] == 'none':
np.savetxt(calib_ascii,np.transpose(output), fmt=[wave_fmt, flux_fmt ])
else:
np.savetxt(calib_ascii,np.transpose(output), fmt=[wave_fmt, flux_fmt ], header=hea)
print('Saved %s' % (calib_ascii))
gofio/drslib/nodding.py
View file @ 09c81308
......@@ -299,8 +299,8 @@ class GBNodding():
try:
masterflat = db.extract_dbfile(self.dbconn,'flatstr')
except:
# masterflat = False
# if not masterflat:
masterflat = False
if not masterflat:
# db.copy_dbfile(self.dbconn,'flatstr')
# masterflat = db.extract_dbfile(self.dbconn,'flatstr')
# self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
......@@ -308,6 +308,7 @@ class GBNodding():
db.copy_dbfile(self.dbconn,'flatstr')
masterflat = db.extract_dbfile(self.dbconn,'flatstr')
self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
except:
self.messages.append('No masterflat found in the calibration database, it is not possible to identify the orders. The spectra will not be reduced.')
return
......@@ -411,6 +412,10 @@ class GBNodding():
fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
x0 = np.zeros(CONFIG['N_ORD'])
x1 = np.zeros(CONFIG['N_ORD'],dtype=int)
x2 = np.zeros(CONFIG['N_ORD'],dtype=int)
all_cosmics = 0
for x in xrange(CONFIG['N_ORD']):
start = x*CONFIG['W_ORD']
......@@ -435,9 +440,9 @@ class GBNodding():
ordermasked = np.ma.MaskedArray(order,mask=omask)
goodorder = np.ma.compress_rows(ordermasked)
optSpectrum[x],varOptFlux,profile,x1,x2,cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
optSpectrum[x],varOptFlux,profile,x0[x],x1[x],x2[x],cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
all_cosmics = all_cosmics + cosmics
#print x0
#snr.append(round(max( (np.mean(optSpectrum[x][1000:1050])/np.std(optSpectrum[x][1000:1050])) ,0),2))
#optSpectrum[x] = optSpectrum[x][::-1]
#varOptFlux = varOptFlux[::-1]
......@@ -448,19 +453,21 @@ class GBNodding():
orderlamp = np.ma.MaskedArray(olamp,mask=omask)
goodlamp = np.ma.compress_rows(orderlamp)
extlamp = varie.extract(goodlamp, optSpectrum[x], x1, x2, profile, lgaineff, lroneff)
#print x1[x]
extlamp = varie.extract(goodlamp, optSpectrum[x], x1[x], x2[x], profile, lgaineff, lroneff)
#print 2
osnr = imstr.data[start:end]
if slit_pos == CONFIG['B_POS']:
osnr = -osnr
ordersnr = np.ma.MaskedArray(osnr,mask=omask)
goodsnr = np.ma.compress_rows(ordersnr)
snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1[x], x2[x], profile, gaineff, roneff)
if any(CONFIG['USE_FLAT'].values()) is True:
extflat = varie.extract(norflat, optSpectrum[x], x1, x2, profile, fgaineff, froneff)
extflat = varie.extract(norflat, optSpectrum[x], x1[x], x2[x], profile, fgaineff, froneff)
#print extflat
with np.errstate(divide='ignore', invalid='ignore'):
#fsnr[x] = (np.sqrt(froneff**2 + (fgaineff*(extflat*meanflat))))/(fgaineff*(extflat*meanflat))
......@@ -476,7 +483,7 @@ class GBNodding():
# osnr = -osnr
#ordersnr = np.ma.MaskedArray(osnr,mask=omask)
#goodsnr = np.ma.compress_rows(ordersnr)
#snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
#snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1[x], x2[x], profile, gaineff, roneff)
......@@ -521,7 +528,7 @@ class GBNodding():
optSpectrum = np.asarray(optSpectrum, dtype='float32')
self.messages.append('The nodding %s was extracted.' % str(slit),)
if all_cosmics == 1050:
if all_cosmics == (CONFIG['NCOSMIC']+1)*CONFIG['N_ORD']:
self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
else:
self.messages.append('%s cosmics were removed.' % str(all_cosmics),)
......@@ -567,7 +574,9 @@ class GBNodding():
calname1d = '_'.join((nomebase,sfx))
#calname = str(os.path.basename(fitsfile)).replace('_AB','_B')
else:
print 'Wrong slit position!'
#print 'Wrong slit position!'
self.messages.append('Wrong slit position!')
return
heaspe = fits.Header(hea_ima)
......@@ -632,6 +641,15 @@ class GBNodding():
c4 = fits.Column(name='SNR', format=''.join((str(CONFIG['YCCD']),'D')), array=snr)
heaspe[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])
if slit_pos == CONFIG['A_POS']:
heaspe[CONFIG['AYPOS'][0]] = (round((np.mean(x0)),2), CONFIG['AYPOS'][1])
heaspe[CONFIG['AYPOS_LOW'][0]] = (round((np.mean(x1)),2), CONFIG['AYPOS_LOW'][1])
heaspe[CONFIG['AYPOS_UP'][0]] = (round((np.mean(x2)),2), CONFIG['AYPOS_UP'][1])
elif slit_pos == CONFIG['B_POS']:
heaspe[CONFIG['BYPOS'][0]] = (round((np.mean(x0)),2), CONFIG['BYPOS'][1])
heaspe[CONFIG['BYPOS_LOW'][0]] = (round((np.mean(x1)),2), CONFIG['BYPOS_LOW'][1])
heaspe[CONFIG['BYPOS_UP'][0]] = (round((np.mean(x2)),2), CONFIG['BYPOS_UP'][1])
# Add metadata to header
heaspe = metadata.add_metadata(heaspe)
......@@ -842,6 +860,10 @@ class GBNodding():
heaspe[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])
heaspe[CONFIG['BYPOS'][0]] = (bwaves[CONFIG['BYPOS'][0]], CONFIG['BYPOS'][1])
heaspe[CONFIG['BYPOS_LOW'][0]] = (bwaves[CONFIG['BYPOS_LOW'][0]], CONFIG['BYPOS_LOW'][1])
heaspe[CONFIG['BYPOS_UP'][0]] = (bwaves[CONFIG['BYPOS_UP'][0]], CONFIG['BYPOS_UP'][1])
# Add metadata to header
heaspe = metadata.add_metadata(heaspe)
......
gofio/drslib/stare.py
View file @ 09c81308
......@@ -433,8 +433,8 @@ class GBStare():
try:
masterflat = db.extract_dbfile(self.dbconn,'flatstr')
except:
# masterflat = False
# if not masterflat:
masterflat = False
if not masterflat:
# db.copy_dbfile(self.dbconn,'flatstr')
# masterflat = db.extract_dbfile(self.dbconn,'flatstr')
# self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
......@@ -543,6 +543,10 @@ class GBStare():
fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
x0 = np.zeros(CONFIG['N_ORD'])
x1 = np.zeros(CONFIG['N_ORD'],dtype=int)
x2 = np.zeros(CONFIG['N_ORD'],dtype=int)
all_cosmics = 0
for x in xrange(CONFIG['N_ORD']):
......@@ -568,17 +572,17 @@ class GBStare():
goodorder = np.ma.compress_rows(ordermasked)
# call optimal extraction
optSpectrum[x],varOptFlux,profile,x1,x2,cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
optSpectrum[x],varOptFlux,profile,x0[x],x1[x],x2[x],cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
all_cosmics = all_cosmics + cosmics
olamp = mlamp.data[start:end]
orderlamp = np.ma.MaskedArray(olamp,mask=omask)
goodlamp = np.ma.compress_rows(orderlamp)
extlamp = varie.extract(goodlamp, optSpectrum[x], x1, x2, profile, lgaineff, lroneff)
extlamp = varie.extract(goodlamp, optSpectrum[x], x1[x], x2[x], profile, lgaineff, lroneff)
if any(CONFIG['USE_FLAT'].values()) is True:
extflat = varie.extract(norflat, optSpectrum[x], x1, x2, profile, fgaineff, froneff)
extflat = varie.extract(norflat, optSpectrum[x], x1[x], x2[x], profile, fgaineff, froneff)
#print extflat
with np.errstate(divide='ignore', invalid='ignore'):
#fsnr[x] = (np.sqrt(froneff**2 + (fgaineff*(extflat*meanflat))))/(fgaineff*(extflat*meanflat))
......@@ -592,7 +596,7 @@ class GBStare():
osnr = -osnr
ordersnr = np.ma.MaskedArray(osnr,mask=omask)
goodsnr = np.ma.compress_rows(ordersnr)
snrSpectrum = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
snrSpectrum = varie.extract(goodsnr, optSpectrum[x], x1[x], x2[x], profile, gaineff, roneff)
with np.errstate(divide='ignore', invalid='ignore'):
......@@ -642,7 +646,7 @@ class GBStare():
self.messages.append('The spectrum %s was extracted.' % str(os.path.basename(straight)),)
if all_cosmics == 1050:
if all_cosmics == (CONFIG['NCOSMIC']+1)*CONFIG['N_ORD']:
self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
else:
self.messages.append('%s cosmics were removed.' % str(all_cosmics),)
......@@ -716,6 +720,10 @@ class GBStare():
heaspe[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])
heaspe[CONFIG['YPOS'][0]] = (round((np.mean(x0)),2), CONFIG['YPOS'][1])
heaspe[CONFIG['YPOS_LOW'][0]] = (round((np.mean(x1)),2), CONFIG['YPOS_LOW'][1])
heaspe[CONFIG['YPOS_UP'][0]] = (round((np.mean(x2)),2), CONFIG['YPOS_UP'][1])
# Add metadata to header
heaspe = metadata.add_metadata(heaspe)
......
gofio/drslib/varie.py
View file @ 09c81308
......@@ -430,7 +430,7 @@ def optExtract(data,gain,ron,slit_pos,ordine):
#warnings.simplefilter('default',RuntimeWarning)
warnings.simplefilter("default", optimize.OptimizeWarning)
return OptFlux, varOptFlux, profile, x1, x2, stop
return OptFlux, varOptFlux, profile, x0, x1, x2, stop
#--------------------- Extraction with a pre-determined profile -------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment