Commit c0180681 authored by Monica Rainer's avatar Monica Rainer
Browse files

Faster interpolation, fixed calibration bug

parent 3dea2972
......@@ -87,6 +87,7 @@ CONFIG['RED_STR'] = os.path.join(CONFIG['RED_DIR'],'STR')
CONFIG['RES_DIR'] = os.path.join(CONFIG['APP_DIR'], 'resources')
CONFIG['BADPIX_MASK'] = os.path.join(CONFIG['RES_DIR'], 'badpix_mask.fits')
CONFIG['LOGGER'] = {'TERM':True,'FILE':True} # logger writes to terminal and/or to file
CONFIG['LOG_FILE'] = os.path.join(CONFIG['RED_DIR'], 'drs.log')
CONFIG['STRAIGHT'] = os.path.join(CONFIG['RES_DIR'],'straight_giano_2D')
CONFIG['STRAIGHT_OPT'] = 'I=1'
......@@ -148,6 +149,8 @@ CONFIG['DRS_MJD'] = (' '.join((CONFIG['KEY_DRS'],'MJD')),'MJD of the combined im
CONFIG['BERV'] = (' '.join((CONFIG['KEY_DRS'],'BERV')),'Barycentric correction [km/s]')
CONFIG['HJD'] = (' '.join((CONFIG['KEY_DRS'],'HJD')),'HJD of the combined image')
CONFIG['AIRMASS'] = (' '.join((CONFIG['KEY_DRS'],'AIRMASS')),'Airmass of the combined image')
CONFIG['MASTERFLAT'] = (' '.join((CONFIG['KEY_DRS'],'FLAT')),'Flat-field used for the reduction')
CONFIG['MASTERLAMP'] = (' '.join((CONFIG['KEY_DRS'],'LAMP')),'Calibration lamp used for the reduction')
#
# DB configuration
......@@ -189,7 +192,7 @@ CONFIG['B_POS'] = 23. # guess pixel position for B
CONFIG['C_POS'] = 15. # guess pixel position for C
CONFIG['HWTM'] = 6. # hwtm limit for contamination between A and B
CONFIG['S1D'] = True # flag for creation of s1d output (True/False)
CONFIG['S1D_STEP'] = 0.002 # constant step to use in s1d spectra
CONFIG['S1D_STEP'] = 0.001 # constant step to use in s1d spectra
#
# Wavelength calibration
......
......@@ -5,7 +5,7 @@ Creation and settings of the calibrations' database.
Connection to tre ramp processor database.
"""
import sqlite3, os
import sqlite3, os, glob
from drslib.config import CONFIG
from collections import OrderedDict
......@@ -47,6 +47,47 @@ def connect_db_ramp():
return dbramps
def remove_dbfile(night,cal=False):
"""
Removes all instances of a calibration file from the database and
from the calibrations directory.
Requires 2 args:
- night in the same format used in the filename (yyyy_mm_dd or yyyy-mm-dd)
- calibrations suffix (optional): UNE, FP, DARK, FLAT.
If missing, all the calibration files of the night will be removed
and the entire row for that night will be removed.
How to call from gofio directory:
from drslib.db import remove_dbfile
remove_dbfile('2017_06_20','UNE')
"""
if cal:
filename = ''.join(['%',night,'%',cal,'%'])
globname = ''.join([CONFIG['CALIB_DIR'],'/*',night,'*',cal,'*'])
else:
filename = ''.join(['%',night,'%'])
globname = ''.join([CONFIG['CALIB_DIR'],'/*',night,'*'])
night = night.replace('_','-')
dbcalib = sqlite3.connect(CONFIG['DB_CALIB_PATH'])
cursor = dbcalib.cursor()
#print CONFIG['DB_CALIB_COLS']
for colname in CONFIG['DB_CALIB_COLS'].values():
cursor.execute("UPDATE {table} SET {calib}=NULL WHERE {calib} LIKE ?"\
.format(table='calibrations',calib=colname),(filename,))
if not cal:
cursor.execute("DELETE FROM {table} WHERE {dcalib}=?"\
.format(table='calibrations',dcalib=CONFIG['DB_CALIB_COLS']['01_data']),(night,))
dbcalib.commit()
for f in glob.glob(globname):
os.remove(f)
return
# ---------- OPERATIONS: calibration database -----------------
......
......@@ -419,6 +419,11 @@ class GBNodding():
drs_mjd = float(heaspe[CONFIG['KEYS']['MJD']]) + (float(heaspe[CONFIG['KEYS']['EXPTIME']])/(2*86400))
heaspe[CONFIG['DRS_MJD'][0]] = (drs_mjd,CONFIG['DRS_MJD'][1])
if masterflat:
heaspe[CONFIG['MASTERFLAT'][0]] = (os.path.basename(masterflat),CONFIG['MASTERFLAT'][1])
if masterlamp:
heaspe[CONFIG['MASTERLAMP'][0]] = (os.path.basename(masterlamp),CONFIG['MASTERLAMP'][1])
snr[snr==np.inf] = 0
snr[snr==-np.inf] = 0
snr = np.nan_to_num(snr)
......@@ -462,16 +467,16 @@ class GBNodding():
if CONFIG['S1D']:
#s1d = varie.create_s1d(optSpectrum,snr,heaspe)
s1d = varie.create_s1d(optSpectrum,heaspe)
s1d, startval = varie.create_s1d(optSpectrum,heaspe)
heaspe[CONFIG['KEYS']['FILENAME']] = calname1d
heaspe['CRPIX1'] = (1.,'Reference pixel')
heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
heaspe['BUNIT'] = ('Relative Flux','Units of data values')
s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
s1dfits.writeto(calname1d,clobber=True)
......@@ -584,15 +589,15 @@ class GBNodding():
if CONFIG['S1D']:
#s1d = varie.create_s1d(abcalib,snr,heaspe)
s1d = varie.create_s1d(abcalib,heaspe)
s1d, startval = varie.create_s1d(abcalib,heaspe)
heaspe['CRPIX1'] = (1.,'Reference pixel')
heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
heaspe['BUNIT'] = ('Relative Flux','Units of data values')
s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
s1dfits.writeto(abnome1d,clobber=True)
......
......@@ -525,6 +525,11 @@ class GBStare():
heaspe = fits.Header(imstr.header)
if masterflat:
heaspe[CONFIG['MASTERFLAT'][0]] = (os.path.basename(masterflat),CONFIG['MASTERFLAT'][1])
if masterlamp:
heaspe[CONFIG['MASTERLAMP'][0]] = (os.path.basename(masterlamp),CONFIG['MASTERLAMP'][1])
snr[snr==np.inf] = 0
snr[snr==-np.inf] = 0
snr = np.nan_to_num(snr)
......@@ -563,15 +568,15 @@ class GBStare():
if CONFIG['S1D']:
#s1d = varie.create_s1d(optSpectrum,snr,heaspe)
s1d = varie.create_s1d(optSpectrum,heaspe)
s1d, startval = varie.create_s1d(optSpectrum,heaspe)
heaspe['CRPIX1'] = (1.,'Reference pixel')
heaspe['CRVAL1'] = (s1d.keys()[0],'Coordinates at reference pixel')
heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')
heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
heaspe['BUNIT'] = ('Relative Flux','Units of data values')
s1dfits = fits.PrimaryHDU(np.asarray(s1d.values()),header=heaspe)
s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
s1dfits.writeto(calname1d,clobber=True)
......
......@@ -30,7 +30,7 @@ import warnings
#import matplotlib.pyplot as plt
#from scipy import optimize, interpolate, signal
from scipy import optimize
from scipy import optimize, interpolate
from collections import OrderedDict
import time
......@@ -764,14 +764,15 @@ def wcalib(heawave,o):
#--------------------- Linear interpolation and rebinning -------------------
def rebin_linear(heawave,flux_old,heawave_old,o):
#def rebin_linear(heawave,flux_old,heawave_old,o):
def rebin_linear(wave,wave_old,flux_old):
"""
Rebin spectrum from wave_old to wave, interpolating linearly
"""
wave = wcalib(heawave,o)[::-1]
wave_old = wcalib(heawave_old,o)[::-1]
flux_old = flux_old[::-1]
#wave = wcalib(heawave,o)[::-1]
#wave_old = wcalib(heawave_old,o)[::-1]
#flux_old = flux_old[::-1]
#print wave
#print wave_old
......@@ -806,14 +807,15 @@ def rebin_linear(heawave,flux_old,heawave_old,o):
#--------------------- Parabolic interpolation and rebinning -------------------
def rebin2deg(heawave,flux_old,heawave_old,o):
#def rebin2deg(heawave,flux_old,heawave_old,o):
def rebin2deg(wave,wave_old,flux_old):
"""
Rebin spectrum from wave_old to wave, interpolating 2 degree
"""
wave = wcalib(heawave,o)[::-1]
wave_old = wcalib(heawave_old,o)[::-1]
flux_old = flux_old[::-1]
#wave = wcalib(heawave,o)[::-1]
#wave_old = wcalib(heawave_old,o)[::-1]
#flux_old = flux_old[::-1]
#print wave
#print wave_old
......@@ -840,42 +842,43 @@ def rebin2deg(heawave,flux_old,heawave_old,o):
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
flux_new.append(flux)
flux_new.append(flux)
except:
try:
try:
y0 = flux_old[iw-2]
y1 = flux_old[iw-1]
y2 = flux_old[iw]
x0 = wave_old[iw-2]
x1 = wave_old[iw-1]
x2 = wave_old[iw]
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
flux_new.append(flux)
except:
try:
except:
try:
y0 = flux_old[iw]
y1 = flux_old[iw+1]
y2 = flux_old[iw+2]
x0 = wave_old[iw]
x1 = wave_old[iw+1]
x2 = wave_old[iw+2]
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
flux_new.append(flux)
except:
except:
flux_new.append(flux_old[iw])
#plt.plot(flux_old)
#plt.plot(flux_new)
#plt.show()
......@@ -883,14 +886,53 @@ def rebin2deg(heawave,flux_old,heawave_old,o):
return np.asarray(flux_new)[::-1]
#return np.asarray(flux_new)
#--------- Switch between linear and parabolic interpolation and rebinning --------
def rebin(heawave,flux_old,heawave_old,o):
"""
Decide which rebinning to use: linear or parabolic
Decide which rebinning to use: linear, parabolic, np.interp or scipy UnivariateSpline
"""
#return rebin_linear(heawave,flux_old,heawave_old,o)
return rebin2deg(heawave,flux_old,heawave_old,o)
wave = wcalib(heawave,o)[::-1]
wave_old = wcalib(heawave_old,o)[::-1]
flux_old = flux_old[::-1]
#t1 = time.time()
flux_new = np.interp(wave,wave_old,flux_old)
return np.asarray(flux_new)[::-1]
#t2 = time.time()
#spl = interpolate.UnivariateSpline(wave_old,flux_old,k=3,s=0)
#try:
# flux_new = spl(wave)
#except:
# select_wave = wave[wave_old[0]<wave<wave_old[-1]]
# flux_new = spl(select_wave)
# try:
# wave0 = wave[wave<wave_old[0]]
# wave0.fill(flux_old[0])
# flux_new = np.append(wave0,flux_new)
# except:
# pass
# try:
# wave1 = wave[wave>wave_old[-1]]
# wave1.fill(flux_old[-1])
# flux_new = np.append(flux_new,wave1)
# except:
# pass
#t3 = time.time()
#print 'Rebin scipy spline: %s s' % str(t3-t2)
#return flux_new[::-1]
#return rebin2deg(wave,wave_old,flux_old)
#return rebin_linear(wave,wave_old,flux_old)
#--------------------- Check for keywords existence (interactive) -------------------
......@@ -910,7 +952,7 @@ def check_keyraw(header,filename):
else:
ask = False
print 'Missing keyword %s in the header of %s' % (key,filename)
#print 'Missing keyword %s in the header of %s' % (key,filename)
#ask = 'Insert %s value: ' % (key)
value = raw_input('Insert %s value: \n' % (key))
if not value:
......@@ -998,13 +1040,13 @@ def create_s1d(spectrum, header):
wave_new = {}
start = {}
end = {}
s1d = OrderedDict()
s1d = np.asarray([])
wstep = CONFIG['S1D_STEP'] # step in nm of the final s1d file
for o in reversed(xrange(CONFIG['N_ORD'])):
#t1 = time.time()
wfit = -99999
#wfit = -99999
flux_old = spectrum[o][::-1]
# use the central part of the order to normalize it
......@@ -1021,109 +1063,41 @@ def create_s1d(spectrum, header):
wave_old[o] = wcalib(header,o)[::-1]
try:
shift = int((wave_old[o][0] - end[o+1])/wstep) -1
start[o] = end[o+1] + (shift*wstep)
#shift = int((wave_old[o][0] - end[o+1])/wstep) -1
shift = max(int(round((wave_old[o][0] - end[o+1])/wstep,0)),1)
#start[o] = max(end[o+1] + (shift*wstep), end[o+1]+wstep)
start[o] = round(end[o+1] + (shift*wstep),3)
except:
start[o] = round(wave_old[o][0],3)
nstep = int(((wave_old[o][-1]-start[o])/wstep))+1
#nstep = int(round(((wave_old[o][-1]-start[o])/wstep),0))+1
nstep = int(round(((wave_old[o][-1]-start[o])/wstep),0))
end[o] = round(start[o]+(nstep*wstep),3)
wave_new[o] = np.arange(start[o],end[o]+(wstep/10),wstep)
#print o
#print start[o]
#print end[o]
#wave_new[o] = np.linspace(start[o],end[o],nstep)
wave_new[o] = np.arange(start[o],end[o],wstep)
#print wave_new[o]
#print len(wave_new[o])
spl = interpolate.UnivariateSpline(wave_old[o],flux_old,k=3,s=0)
for w in wave_new[o]:
w = round(w,3)
#print w
iw = min(np.searchsorted(wave_old[o],w),len(wave_old[o])-1)
try:
f1 = s1d[w]
continue
except:
pass
if round(wave_old[o][iw],3) == w:
flux = flux_old[iw]
else:
if wfit == iw:
try:
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
except:
flux = flux_old[iw]
else:
try:
flux_new = spl(wave_new[o])
except:
select_wave = wave_new[o][wave_old[o][0]<=wave_new<=wave_old[o][-1]]
flux_new = spl(select_wave)
start[o] = select_wave[0]
end[o] = select_wave[-1]
try:
#print 'try 1'
y0 = flux_old[iw-1]
y1 = flux_old[iw]
y2 = flux_old[iw+1]
flux_new[abs(flux_new)<1e-04] = 0
#print len(flux_new)
x0 = wave_old[o][iw-1]
x1 = wave_old[o][iw]
x2 = wave_old[o][iw+1]
if s1d.any():
gap = start[o]- (end[o+1]+wstep)
ngap = max(int(round(gap/wstep,0)),0)
s1d = np.append(s1d,np.zeros(ngap))
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
except:
try:
#print 'try 2'
y0 = flux_old[iw-2]
y1 = flux_old[iw-1]
y2 = flux_old[iw]
x0 = wave_old[o][iw-2]
x1 = wave_old[o][iw-1]
x2 = wave_old[o][iw]
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
except:
try:
#print 'try 3'
y0 = flux_old[iw]
y1 = flux_old[iw+1]
y2 = flux_old[iw+2]
x0 = wave_old[o][iw]
x1 = wave_old[o][iw+1]
x2 = wave_old[o][iw+2]
flux = y0*(w-x1)*(w-x2)/((x0-x1)*(x0-x2)) + \
y1*(w-x0)*(w-x2)/((x1-x0)*(x1-x2)) + \
y2*(w-x0)*(w-x1)/((x2-x0)*(x2-x1))
except:
#print 'no fit'
flux = flux_old[iw]
wfit = iw
if s1d:
last = next(reversed(s1d))
while w-last > wstep:
last = round(last+wstep,3)
s1d[last] = 0.0
#print 'qui'
#print w
s1d[w] = flux
s1d = np.append(s1d,flux_new)
#t2 = time.time()
#print 's1d spectrum order %s: %s s' % (str(o),str(t2-t1))
#print s1d
#print len(s1d)
return s1d
return s1d, start[CONFIG['N_ORD']-1]
......@@ -193,15 +193,15 @@ class GBWls():
self.messages.append(' *** WARNING ***')
self.messages.append('The default wavelength calibration for the order %s will be taken from the database and as such it will not be optimal for the night.' % str(x+32))
try:
wcalib = db.extract_dbfile(self.dbconn,'une_str')
wcalib = db.extract_dbfile(self.dbconn,'une_calib')
except:
wcalib = False
if not wcalib:
self.messages.append('No UNe calibration found for this night, it will be taken from the calibration database.')
db.copy_dbfile(self.dbconn,'une_str')
wcalib = db.extract_dbfile(self.dbconn,'une_str')
db.copy_dbfile(self.dbconn,'une_calib')
wcalib = db.extract_dbfile(self.dbconn,'une_calib')
try:
wlc = ccdproc.CCDData.read(wcalib, unit=u.adu)
......
......@@ -303,6 +303,15 @@ def GianoBreduce(rawfits,lists,dbcalib,working,group,do_calib,stop):
if __name__ == "__main__":
rawdir = CONFIG['RAW_DIR']
if not os.path.exists(rawdir):
drslogger.log('The RAW directory %s does not exist. GOFIO will now close down.' % rawdir)
drslogger.log('Pipeline stopped.')
sys.exit(0)
dbcalib.close()
if not CONFIG['OFFLINE']:
dbramps.close()
dbnight.close()
rawlists = {CONFIG['DARK']:[],CONFIG['FLAT']:[],CONFIG['WCAL_UNE']:[],CONFIG['WCAL_FP']:[],CONFIG['SCIENCE']:[]}
......
......@@ -17,8 +17,10 @@ class DrsLogger():
self.filelogfmt = '%(asctime)s.%(msecs)d * %(message)s'
self.conslogfmt = '%(message)s'
self.logger.setLevel(logging.INFO)
self.addfilelogger()
self.addconsolelogger()
if CFG['LOGGER']['FILE']:
self.addfilelogger()
if CFG['LOGGER']['TERM']:
self.addconsolelogger()
def addfilelogger(self):
"""The logger to the LOG_FILE is defined"""
......
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