# Copyright 2014 - 2015 Anna-Lea Lesage for MASCARA
#
# This file is part of the mascara package.
# mascara is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# mascara is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with mascara. If not, see <http://www.gnu.org/licenses/>.
'''
.. module:: reduction
.. moduleauthor:: Anna-Lea Lesage
Module containing the tow main reduction programs for the MASCARA pipeline.
- Reduction proceeds to the reduction of raw images following this pattern with a
cadence of 50 images:
- the fits file is read, image and header are retrieved
- the LST index is read from the header and converted in an image number, slowind,
within a set of 50
- for a given slowindex the astrometric solution is verified
- for each star within the CCD, the flux is calculated. The values are saved
in a dictionary under the star's identifier entry
- the raw image is rotated and coadded to the previous ones to create a binned image
- once theslowind reaches 50, the dictionary is saved as an hdf5 file to disk;
the binned image is averaged and save as a compressed fits file;
for each star down to magnitude 10 within the binned image, the flux is calculated;
the binned dictionary is saved as hdf5 file to disk.
- finally all values are resetted.
- SaveItAll is dedicated to the saving of the light curves. It proceeds as following:
- the last hdf5 file created is read to memory,
- the stellar identifiers are retrieved,
- the portion of light curve of each star is saved into a temporary directory
to disk in pkl format
- if the star has exited the CCD, the different pickles files are combined into
a final fits file.
- at the end of the night, all remaining pickle files per star are combined into
a fits file, and the temporary directory is deleted
'''
__version__ = '15.06.28'
### *** Latest Changes Record *** ###
### 03.12.2014 Changes in header for slow curves, removed SAPER by APER ###
### 04.12.2014 Found bug in saving binned image-> wrong header passed to Masfunct.saveCompFits
### 05.12.2014 Found bug in making slow dictionary. Need to save 27 entries, forgot lpstseq
### 07.01.2015 Removed the calculation of lstid in the fast cadence. Inserted instead
### the lpstidx value from the raw image header.
### The sequence is still calculated. lstid is recalculated in the slow
### cadence with the observation site.
### 12.01.2015 Reduced the magnitude interval for calculating the astrometry: takes
### too long when using vmag < 8.4. Trying vmag<8.
### 30.01.2015 Corrected bug line 648 with non existing header keys.
### 02.03.2015 Updated program for the new stellar catalog which now contains flag
### for blended stars and the corresponding blending values
### 02.03.2015 Implemented the new astrometric routine
### 01.04.2015 Implemented routines send_message and reduction_summary to the pipeline
### 09.04.2015 Replaced the previous stacking routine by the Stacker class.
### Cleaned up the dates in jd and calendar. Removed 3.2 seconds to the name of the
### binned images.
### 30.04.2015 Added a try?except Empty condition at the end _of_night part in Reduction to avoid
### crashed in th eprogram when the saving takes longer than planned....
### 06.05.2015 Added memory release-conditions in SaveItAll to no longer update all the slow
### databases.
### 22.05.2015 Major change in the output format of the night:
### On Linux machines, the databases are stacked together and NOT deleted at the end
### On Windows machines, the databases are not stacked and NOT deleted either!!!
### 25.06.2015 Added a condition to restart the reduction after it was interrupted
### **************************** ###
import pyfits as pf
import numpy as np
import pandas as pd
import logging
import os
from datetime import datetime, timedelta
import time
import platform
import multiprocessing
import gc
from mascara.reduction import masc, config
import mascara.funcs as mfuncs
import mascara.observer as observer
import mascara.phot as phot
from mascara.funcs.utils import QueueHandler, QueueListener, Stacker
system = platform.system().lower()
#####################################
##### SAVING FUNCTION ########
[docs]def SaveItAll(*args, **kwargs):
''' SaveItAll creates or updates an existing fits file for every star
present in the StarID dictionnary.
Inputs:
- PathfLC, PathsLC, the path root for saving the light curves (excluding \\)
- StarID, a dictionary. The key is the stellar identifier
No outputs, the files are directly written to disk
'''
saveout, event, qlog = args
qh = QueueHandler(qlog)
logger = logging.getLogger(multiprocessing.current_process().name)
logger.setLevel(logging.INFO)
logger.addHandler(qh)
tempo = []
if system == 'windows':
breakDataBase = True
else:
breakDataBase=False
savesequence = True
event.set()
sIDs = mfuncs.utils.loadStellarCatalog(filename='StarCat.fits', sIDonly=True, \
sID = 'ASCC', Vmag=[4.75, 4.8])
while savesequence:
tempo.append(saveout.get())
tmp = tempo.pop(0)
to = datetime.now()
if isinstance(tmp, str) and ((tmp.lower() == 'stop') or (tmp.lower() == 'quit')):
logger.warning('Received Quit command')
logger.warning('Closing this end of the pipe')
saveout.close()
savesequence = False
return
elif len(tmp)==4:
# Read the temporary database and save...
logger.debug('Saving the database in single stellar files...')
directory, ndb, hdr, extID = tmp
### Modified by JS
tmpPath = os.path.join(directory,'tmp', '')
PathfLC_ = os.path.join(directory, config.PathfLC, '')
PathsLC_ = os.path.join(directory, config.PathsLC, '')
if not os.path.exists(PathfLC_):
mfuncs.mio.makedir(directory, PathfLC_, PathsLC_)
##### END of the NIGHT #####
if extID == 'all':
import shutil
### Start the event to saving
event.clear()
#### Get the information of aperture used:
apr = np.array(config.apr).ravel()
naper = apr.size
sapr = np.array(config.sapr).ravel()
saper = sapr.size
### Combine the fast light curves!
night = ndb
CamID = ''.join([t for t in hdr['Station'] if t.isupper()]) + hdr['Camera'][0]
## Proceed with all the remaining high cadence files
totID = mfuncs.mio.list_directories(tmpPath, '_')
if len(totID) > 0:
logger.info('Combining the light curves to fits file ... ')
for s in totID:
sid, ndb = os.path.basename(s).split('_')
if int(ndb) < 1e3:
masc.CombinePkl(s, PathfLC_, sid)
else:
masc.CombinePkl(s, PathsLC_, sid)
logger.info('saved_all_LC')
### Save the slow light curves to database:
slowDB = mfuncs.mio.find_files(tmpPath, 'tmpdb30*', format='.h5')
slowname = 'sLC_'+night+CamID
### change NAPER parameters
for i, e in enumerate(sapr):
hdr['aper'+str(i)] = e
hdr['naper'] = saper
hdr.totextfile(os.path.join(PathsLC_, 'Header.txt'), clobber=True)
if not breakDataBase:
globdic = {}
for h,k in hdr.iteritems():
globdic[h.lower()]=k
try:
logger.info('Combining slow DB')
masc.CombineH5(PathsLC_, slowDB, globdic, name = slowname)
except Exception, e:
logger.error('Did not save slow LC', exc_info=True)
print 'Error somewhere', e
### Finally move it to its directory
for sDB in slowDB:
shutil.move(sDB, PathsLC_)
### Change NAPER parameters:
for e in range(naper, saper):
hdr.pop('aper'+str(e))
for i, e in enumerate(apr):
hdr['APER'+str(i)] = e
hdr['naper'] = naper
hdr.totextfile(os.path.join(PathfLC_, 'Header.txt'), clobber=True)
### Same for the fast light curves
fastDB = mfuncs.mio.find_files(tmpPath, 'tmpdb*', format='.h5')
fastname = 'fLC_'+night+CamID
### If single database: index the array inside
if not breakDataBase:
globdic = {}
for h,k in hdr.iteritems():
globdic[h.lower()]=k
try:
logger.info('Combining the fastDB')
masc.CombineH5(PathfLC_,fastDB, globdic, name=fastname)
except Exception:
logger.error('Fast DB', exc_info=True)
for fDB in fastDB:
shutil.move(fDB, PathfLC_)
mfuncs.mio.cleanup(tmpPath, qlog)
logger.info('moved_all_databases')
event.set()
##### DURING the NIGHT #####
else:
print ' Treating database %g' %ndb
logger.info('Connecting to database {}'.format(ndb))
print ' Connecting to database'
### Connecting to FAST database:
# extract just the infos for the selected stars.
if ndb < 1e3:
store = pd.HDFStore(os.path.join(tmpPath,'tmpdb'+str(ndb)+'.h5'), \
'r', complevel=5, complib='zlib')
cnxN = store.get('tmpdb'+str(ndb))
store.close()
### make light curves for a selection of stars for the verifications
for sid in sIDs:
if sid in cnxN.keys():
tmp = cnxN[sid]
saveloc = os.path.join(tmpPath, sid+'_100')
try:
tmp.to_pickle(os.path.join(saveloc,sid+str(ndb)+'.pkl'))
except IOError:
os.mkdir(saveloc)
tmp.to_pickle(os.path.join(saveloc,sid+str(ndb)+'.pkl'))
hdr.totextfile(os.path.join(saveloc, 'Header'+sid+'.txt'))
### Connecting to SLOW database
# Extract the infos for the selected stars
# and combines every 20 databases in one.
else:
if os.path.exists(tmpPath + 'tmpdb'+str(ndb-1)+'.h5'):
if ndb%20 ==0 or extID == 'last':
store = pd.HDFStore(os.path.join(tmpPath,'tmpdb'+str(ndb-1)+'.h5'), \
'r', complevel=5, complib='zlib')
cnxN = store.get('tmpdb'+str(ndb-1))
store.close()
### make light curves for a selection of stars for the verifications
for sid in sIDs:
if sid in cnxN.keys():
tmp = cnxN[sid]
saveloc = os.path.join(tmpPath, sid+'_10100')
try:
tmp.to_pickle(os.path.join(saveloc,sid+str(ndb)+'.pkl'))
except IOError:
os.mkdir(saveloc)
tmp.to_pickle(os.path.join(saveloc,sid+str(ndb)+'.pkl'))
hdr.totextfile(os.path.join(saveloc, 'Header'+sid+'.txt'))
else:
try:
### Load the new in memory
store = pd.HDFStore(tmpPath+'tmpdb'+str(ndb)+'.h5', 'r', \
complevel=5, complib='zlib')
cnxN = store.get('tmpdb'+str(ndb))
store.close()
cnxN = cnxN.to_dict()
### load the previous one in memory:
store = pd.HDFStore(tmpPath+'tmpdb'+str(ndb-1)+'.h5','r', \
complevel=5, complib='zlib')
cnx0 = store.get('tmpdb'+str(ndb-1))
store.close()
cnx0 = cnx0.to_dict()
logger.debug('Connection with {} and {} succeeded'.format(
str(ndb), str(ndb-1)))
for k, v in cnxN.iteritems():
if k in cnx0.keys():
d = cnx0[k]
try:
d[8] = np.vstack((d[8], v[8]))
except ValueError:
continue
cnx0[k] = d
else:
cnx0[k] = v
masc.tmpdbSave(tmpPath, cnx0, ndb, qlog)
print 'Updating the stellar database successful'
logger.info('Updated the stellar database successfully')
removeOldDatabase=True
except Exception:
logger.error('Failed in updating the database', exc_info=True)
print 'Error in updating the table'
removeOldDatabase=False
cnxN = 0
del cnx0
if removeOldDatabase:
try:
os.remove(tmpPath+'tmpdb'+str(ndb-1)+'.h5')
except Exception:
logger.error('Did not remove database', exc_info=True)
Tro = (datetime.now() - to).total_seconds()
logger.info('Finished Updating in {} seconds'.format(Tro))
else:
logger.error(' No correct input... exiting!')
savesequence=False
logger.warning('Exiting the program now')
return
##################################
##### REDUCTION FUNCTION ########
[docs]def Reduction(*args, **kwargs):
''' Reduction does the photometric reduction of the images. It either starts in
parallel to the exposure sequence. In this case, the program reduces the first
file in the Queue. Else, if no camera is connected, it reads all the files
in the directory and reduces them.
Inputs:
- directory, a string giving the path to the files
- done_queue, a Queue object which is filled by the camera process
with the name of the next image to reduce.
The done_queue is a tuple (tagstr, camera, filename). Tagstr is dark, flat...
- site, a Site object ;
- savout, a Queue object to communicate with the saving process.
- event, an multiprocessing event object
- qlog, a queue used in the logging.
Outputs: - none.
'''
import Queue
directory, done_queue, site, saveout, event, qlog = args
####### ######
### Prepare for the general loop ###
####### ######
verbose = kwargs.get('verbose', False)
quiet = kwargs.get('quiet', False)
noPhot = kwargs.get('noPhot', False)
Transmission = kwargs.get('transmission', False)
maglimit = config.maglimit
faintlimit = config.faintlimit
RFsite = masc.SiteReference()
###########################################
### Read in star catalog for astrometry:
###########################################
b = pf.open('StarCat.fits') # from Kharchenko 2009
totcat = b[1].data
b.close()
# Extract Ra, Dec from the brightcat. All those must have same length.
tmag = totcat['vmag']
# Prepare general stellar Catalog
mrange = (tmag < maglimit).nonzero()[0]
srange = (tmag < faintlimit).nonzero()[0]
Sptype = totcat['SpType'][mrange]
mag = np.float64(tmag[mrange])
bmag = np.float64(totcat['bmag'][mrange])
ra = totcat['_RAJ2000'][mrange]
dec = totcat['_DEJ2000'][mrange]
ID = totcat['ASCC'][mrange]
flagblended = totcat['flag'][mrange]
#flagblended[flagblended=='']='-'
blendedvalue = totcat['blend'][mrange]
# Removing multiple occurrence of the same ID...it represents around 200 stars
ID, gind = np.unique(ID, return_index=True)
mag, bmag, ra, dec = mag[gind], bmag[gind], ra[gind], dec[gind]
Sptype, flagblended, blendedvalue = Sptype[gind], flagblended[gind], \
blendedvalue[gind]
ID = ID.astype('S10')
## the 'very' large catalog to use on the binnec images
sSptype = totcat['SpType'][srange]
smag = np.float64(tmag[srange])
sbmag = np.float64(totcat['bmag'][srange])
sra = totcat['_RAJ2000'][srange]
sdec = totcat['_DEJ2000'][srange]
## Preparing the identifiers for the slow catalog
sID = totcat['ASCC'][srange]
sflagblended = totcat['flag'][srange]
#sflagblended[sflagblended=='']='-'
sblendedvalue = totcat['blend'][srange]
# Removing multiple occurrence of the same ID...it represents around 200 stars
sID, gind = np.unique(sID, return_index=True)
sSptype, smag, sbmag, sra, sdec = sSptype[gind], smag[gind], sbmag[gind],\
sra[gind], sdec[gind]
sflagblended, sblendedvalue = sflagblended[gind], sblendedvalue[gind]
sID = sID.astype('S10')
# delete the unused catalog from memory:
totcat = 0.
del totcat
# Remove empty strings from the Sptype chararray:
empt = np.where(Sptype == '')
Sptype[empt] = 'XX'
print 'Creating the stellar ID catalog successfull'
# Remove empty strings from the Sptype chararray:
empt = np.where(sSptype == '')
sSptype[empt] = 'XX'
print 'Creating the slow stellar ID catalog successfull'
##### ################## #####
### INITIALISATION OF THE LOOP ###
##### ################## #####
redsequence = True
localdate = datetime.now().strftime(format='%Y%m%d')
Path = os.path.join(directory, localdate + site.sname, '')
PathsI = config.PathsI
logdir = os.path.join('logs', '')
slowbinning = config.slowbinning
newnight = True
binThumbnails = config.binThumbnails
tilesize = config.tilesize
Cam = observer.Camera(name='TestCamera')
softversion = __version__
Ndarks = config.nDarks
Nflats = config.nFlats
texp = config.expTimeObs
apr = np.array([config.apr]).ravel()
naper = apr.size
sapr = np.array([config.sapr]).ravel()
nsaper = sapr.size
skyrad = config.skyrad
obstime = []
tempo = []
date = datetime.now()
##### ########### #####
### STARTING THE LOOP ###
##### ########### #####
while redsequence:
tinit = datetime.now()
## 1. Get input from camera :
try:
tempo.append(done_queue.get())
temp = tempo.pop(0)
# during daytime, the camera does not send data, the queue is empty
except multiprocessing.queues.Empty():
temp = 'wait'
### DAY TIME SEQUENCE ###
if isinstance(temp, str):
if (temp.lower()=='quit')+(temp.lower()=='stop') :
print 'Exiting program now...'
saveout.put('STOP')
if done_queue.empty():
done_queue.close()
else:
done_queue.get()
done_queue.close()
redsequence=False
saveout.close()
return
elif (temp.lower() == 'wait'):
print temp
time.sleep(6.5)
### NIGHT TIME SEQUENCE ###
elif isinstance(temp,tuple):
tagstr, Camfile, filename, Nn = temp
print '{} -- {} -- {}'.format(tagstr, Nn, filename)
### Beginning of the night:
if newnight:
## load camera instance:
Cam = observer.Camera(Camfile, filename=Camfile)
## Create the directory for the next night
localdate = tinit.strftime(format='%Y%m%d')
Path = filename.split('raw')[0]
if not quiet:
print Path
Pathlog = os.path.join(Path, logdir, '')
PathsI_ = os.path.join(Path, PathsI, '')
tmpPath = os.path.join(Path, 'tmp', '')
binThumbs = os.path.join(Path, binThumbnails)
if system == 'linux':
binThumbs = binThumbs.replace('\\', '/')
if not os.path.exists(Pathlog):
mfuncs.mio.makedir(Path, PathsI_, tmpPath, Pathlog, binThumbs)
## Start the logger for the new nigth:
qh = QueueHandler(qlog)
handler = logging.FileHandler(Pathlog+'reduc'+localdate+'.log')
handler.setLevel(logging.INFO)
logger = logging.getLogger(multiprocessing.current_process().name)
if verbose:
logger.setLevel(logging.DEBUG)
else:
logger.setLevel(logging.INFO)
logger.addHandler(qh)
ql = QueueListener(qlog, handler)
ql.start()
formatter = logging.Formatter('%(asctime)s - %(name)s - %(funcName)s - '\
+'%(levelname)s - %(message)s')
handler.setFormatter(formatter)
print 'New Night = True'
## Reset booleans
newnight = False
firstAstrometry=True
usedarks = False
highoffset = False
BadExposure = False
## and reset variables for the next night?
DBNiObs = 100
sDBNiObs = 30100
obstime = []
# Obstime is used both for checking the discontinuities in the data.
slowind = 0.
StarID = {}
timt = 0.
stacker = Stacker(site, Cam, sidex=32, sidey=32)
ndark = 0.
bdark = 0
nbias = 1
mbias=0
sbias=0
checkalt = np.array([])
checkaz = np.array([])
flagit = 0
### If the images are tagged (science images will not be tagged)
logger.debug(' Treating a %s image' %tagstr)
### If the reduction was previously interrupted and is started anew:
if tagstr == 'interrupted':
DBNiObs = filename
sDBNiObs = 30000 + DBNiObs
logger.info('Continue interrupted reduction ... ')
## And load the corresponding darks if possible
darkfilelist= masc.getDarkFilename(Cam.name)
for darkfile in darkfilelist:
if usedarks:
break
try:
dark = pf.getdata(darkfile)
usedarks = True
logger.warning('No enough valid darks, loading existing masterdark')
except Exception:
logger.error('Did not load the masterdark', exc_info=True)
ndark = 0
bdark = 0
### If the incoming image is a calibration "dark" image
elif tagstr == 'dark':
## Case 1. first dark observed
if Nn ==1:
mbias = mbias/nbias
sbias = sbias/nbias
nbias = 1
dim, HDR = pf.getdata(filename, header=True)
acqversion = HDR['AVERSION'] if 'AVERSION' in HDR.keys() \
else 'pre-14.11'
coversion = HDR['CVERSION'] if 'CVERSION' in HDR.keys() \
else 'pre-14.11'
CCDTEMP = HDR['CCDTEMP'] if 'CCDTEMP' in HDR.keys() \
else HDR['CCD-TEMP']
if not quiet:
print ' Reading in dark image for %s' %Cam.name
ndark +=1
### Here we verify that the dark is dark
if dim.mean() > mbias + 2*sbias:
dim = 0
ndark = 0
bdark = 1
## If not, we reject this dark
logger.info(' Reading in dark image for %s', Cam.name)
## Case 2. All the Nth dark observed
else:
datmp = pf.getdata(filename)
### Verification of the quality of the dark
if datmp.mean() < mbias + 1.1*sbias:
dim+=datmp
ndark+=1
else:
bdark+=1
if not quiet:
print ' Stacked %g dark images out of %g'%(ndark, Ndarks)
logger.info('Stacked %g dark images out of %g for %s' \
%(ndark, Ndarks, Cam.name))
### If there are at least 10 good darks and all expected
### darks were analyzed, save the masterdark
if (ndark+bdark) == Ndarks and ndark>10:
logger.info('Creating the master dark')
dark = dim/ndark
now = datetime.now().strftime(format='%Y%m%dT%H%M%S')
header = pf.Header()
header['STATION']= (site.name, 'MASCARA station')
header['CAMERA'] = (Cam.name, 'Camera identifier')
header['ALT0'] = (Cam.Alto, 'General Pointing altitude')
header['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
header['TH0'] = (Cam.Tho, 'Average Pointing orientation')
header['X0'] = (round(Cam.Xo), 'Average Center of the lens -- X')
header['Y0'] = (round(Cam.Yo), 'Average Center of the lens -- Y')
header['EXPTIME']= (round(HDR['EXPTIME'], 1), 'Mean exposure time')
header['NOBS'] = (Ndarks, 'Number of images averaged together')
header['CCDTEMP'] = (CCDTEMP, 'CCD temperature')
header['IMAGETYP'] = ('master'+tagstr, 'Type of image')
header['AVERSION'] = (acqversion, 'Version of the acquisition software')
header['RVERSION'] = (softversion, 'Version of the pipeline software')
header['CVERSION'] = (coversion, 'Version of the control software')
filename = now+site.sname +Cam.sname+'masterdark'
mfuncs.mio.saveCompFits(PathsI_, filename, dark, header)
dim = 0
ndark = 0
usedarks = True
masc.updateDarkFilename(Cam.name, PathsI_+filename)
# delete array to release the memory. If not sufficient we may use
gc.collect()
### If more than 10 (out of the 20) darks are bad, load an
### existing masterdark from the darkTable
elif bdark >=10:
darkfilelist= masc.getDarkFilename(Cam.name)
for darkfile in darkfilelist:
### Once we have one dark loaded, exit the for-loop
if usedarks:
break
try:
dark = pf.getdata(darkfile)
usedarks = True
logger.warning('No enough valid darks, loading existing masterdark')
except Exception:
logger.error('Did not load the masterdark', exc_info=True)
ndark = 0
bdark = 0
### BIAS files
### Use the bias files to define the quality criterion for the darks
elif (tagstr == 'bias'):
if Nn == 1:
nbias = 1
else:
nbias +=1
bias = pf.getdata(filename)
mbias += bias.mean()
sbias += bias.std()
### FLAT files
### Once all the accounted flats are loaded, they are saved
### as masterflats to disk
elif (tagstr == 'flat'):
if Nn == 1:
fim, fhdr = pf.getdata(filename, header=True)
acqversion = fhdr['AVERSION'] if 'AVERSION' in fhdr.keys() \
else 'pre-14.11'
coversion = fhdr['CVERSION'] if 'CVERSION' in fhdr.keys() \
else 'pre-14.11'
CCDTEMP = fhdr['CCDTEMP'] if 'CCDTEMP' in fhdr.keys() \
else fhdr['CCD-TEMP']
else:
fim += pf.getdata(filename)
print ' Stacked %g flat images out of %g'%(Nn, Nflats)
logger.info('Stacked %g flat images out of %g for %s' \
%(Nn, Nflats, Cam.name))
if Nn == Nflats:
logger.info('Creating master flat')
flat = fim/Nflats
now = datetime.now().strftime(format='%Y%m%dT%H%M%S')
header = pf.Header()
CCDTEMP = fhdr['CCDTEMP'] if 'CCDTEMP' in fhdr.keys() else fhdr['CCD-TEMP']
header['STATION']= (site.name, 'MASCARA station')
header['CAMERA'] = (Cam.name, 'Camera identifier')
header['ALT0'] = (Cam.Alto, 'General Pointing altitude')
header['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
header['TH0'] = (Cam.Tho, 'General Pointing orientation')
header['X0'] = (round(Cam.Xo), 'Average Center of the lens -- X')
header['Y0'] = (round(Cam.Yo), 'Average Center of the lens -- Y')
header['EXPTIME']= (fhdr['EXPTIME'], 'Exposure time')
header['NOBS'] = (Nflats, 'Number of images averaged together')
header['CCDTEMP'] = (round(CCDTEMP),'CCD temperature')
header['IMAGETYP'] = ('master'+tagstr, 'Type of image')
header['AVERSION'] = (acqversion, 'Version of the acquisition software')
header['RVERSION'] = (softversion, 'Version of the pipeline software')
header['CVERSION'] = (coversion, 'Version of the control software')
mfuncs.mio.saveCompFits(PathsI_, now+site.sname \
+Cam.sname+'masterflat', flat, header)
fim = 0.
# delete array to release the memory. If not sufficient we may use
gc.collect()
### END of the Night => Command for treating the remaining stars in catalog
### and do a combine all.
elif tagstr == 'end_of_night':
logger.info('Reaching end of the night...')
### Save the rest of the fast dictionary if not empty
if (len(StarID)>0) and DBNiObs > 100:
try:
masc.tmpdbSave(tmpPath, StarID, DBNiObs, qlog)
except Exception:
logger.error('Failed saving the dictionary. FATAL!!!', exc_info=True)
### Do photometry and Save the binned images
sdate = date + mfuncs.timing.calendar2jd(timedelta(
seconds=(slowbinning/2.-slowind-1)*texp))
slpst = site.local_sidereal_time(sdate)
slpstidx, slpstseq = RFsite.local_sidereal_time(sdate, texp*50, lst=False)
### Average the binned image
timt, Nimages = stacker.binImage(sdate)
stime = (mfuncs.timing.jd2calendar(sdate) - \
timedelta(seconds=texp/2)).strftime('%Y%m%dT%H%M%S')
name = stime + site.sname + Cam.sname
### If there are more than 1 images in the stacker instance,
### save those as binned image.
skip = False
if Nimages == 0:
logger.info('No image to save, pass ... ')
skip = True
if not skip:
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
xsalt, xsaz, xvis = site.apparent_stars(sra, sdec, sdate, nutation=True, \
precession=True, refraction=True)
tmpmag, tmpbmag, tmpspec = smag[xvis], sbmag[xvis], sSptype[xvis]
### Re-evaluate the corrections to the positions:
Cam.solve_corrections(timt, xsalt[tmpmag<7.], xsaz[tmpmag<7.])
xxpre, yypre, xyinCCD = Cam.visible_stars(xsalt, xsaz, margin=-50)
xxcor, yycor = Cam.correct_positions(xxpre, yypre)
### use indices to restrain our sample to the stars on the CCD:
xalt, xaz = xsalt[xyinCCD], xsaz[xyinCCD]
tmpra, tmpdec = sra[xvis], sdec[xvis]
sxra, sxdec = tmpra[xyinCCD], tmpdec[xyinCCD]
tmpmag, tmpbmag, tmpspec = smag[xvis], sbmag[xvis], sSptype[xvis]
vsflagblended = sflagblended[xvis]
vsblendedvalue = sblendedvalue[xvis]
sxmag, sxbmag, sxspect = tmpmag[xyinCCD], tmpbmag[xyinCCD], tmpspec[xyinCCD]
xsflagblended = vsflagblended[xyinCCD]
xsblendedvalue=vsblendedvalue[xyinCCD]
ssID = sID[xvis]
ssID = ssID[xyinCCD]
## Make one big dictionary which is saved at once.
slowStarID = {}
if not quiet:
print ' Running the photometry over %g stars in Binned mode' %xxpre.shape[0]
logger.info(' Running the photometry over {}' + \
' stars in Binned mode'.format(xxpre.shape[0]))
### Loop over the star and perform the photometry:
for w in xrange(xxpre.shape[0]):
# ----------------------------------------------------------------
# The final content of the Photometry array is:
# flag, (flux, erflux)*Naperture, sky, ersky, xc, yc,
# xpre, ypre, altitude, azimuth, CCDtemperature,
# exposure time, JDdate, lstdate, lstid, lstseq
# Total number of elements: (2*naper+6)+6
# ----------------------------------------------------------------
try:
fot = phot.OnaPhot(timt, xxcor[w], yycor[w], apr=sapr, \
skyrad=skyrad, flag=flagit)
except Exception:
logger.error('Failed running the photometry ' + \
' on {}'.format(ssID[w], exc_info=True))
continue
tmp = np.hstack((fot, xalt[w], xaz[w], CCDTEMP, \
texp, sdate, slpst, slpstidx, slpstseq))
slowStarID[ssID[w]] = [date, sxra[w], sxdec[w], \
sxmag[w], sxbmag[w], sxspect[w], \
xsflagblended[w], xsblendedvalue[w], \
np.array(tmp)]
hdrslow = pf.Header()
hdrslow['STATION']= (site.name, 'MASCARA station')
hdrslow['CAMERA'] = (Cam.name, 'Camera identifier')
hdrslow['ALT0'] = (Cam.Alto, 'General Pointing altitude')
hdrslow['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
hdrslow['TH0'] = (Cam.Tho, 'General Pointing orientation')
hdrslow['X0'] = (Cam.Xo, 'Average Center of the lens -- X')
hdrslow['Y0'] = (Cam.Yo, 'Average Center of the lens -- Y')
hdrslow['EXPTIME'] = (stacker.texp, 'Mean exposure time in seconds')
hdrslow['CCDTEMP'] = (CCDTEMP, 'CCD temperature at start of exposure')
hdrslow['AVERSION'] = (acqversion, 'Version of the acquisition software')
hdrslow['RVERSION'] = (softversion, 'Version of the pipeline software')
hdrslow['CVERSION'] = (coversion, 'Version of the control software')
hdrslow['NEXP'] = (Nimages, 'Number of images binned together')
hdrslow['NAPER'] = (nsaper, 'Number of aperture used for the photometry')
for i, r in np.ndenumerate(sapr):
hdrslow['APER'+str(i[0])] = (r, 'Aperture radii in pixel')
for j in xrange(naper):
hdrslow['STAPER'+str(j)] = ('XXX', 'Statistics of the high cadence aperture')
hdrslow['NSTAPER'] = (naper, 'Number of high cadence apertures analysed')
hdrslow['SKYRAD0'] = (skyrad[0], 'Inner sky radius')
hdrslow['SKYRAD1'] = (skyrad[1], 'Inner sky radius')
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
#### Now save the faint limit dictionary:
sDBNiObs +=1
logger.info('Save slow dictionary in database {}'.format(sDBNiObs))
try:
masc.tmpdbSave(tmpPath, slowStarID, sDBNiObs, qlog)
except Exception:
logger.error('Failed saving the slow database. ', exc_info=True)
saveout.put((Path, sDBNiObs, hdrslow, 'last'), block=False)
##### Saving the stacked image!!!
### Creating the header for the binned image:
hdrslow['LSTMID'] = (slpst, 'LST at mid-bin observation')
hdrslow['JDMID'] = (sdate, 'JD at mid-bin observation')
hdrslow.set('CRPIX1', Cam.tnx['crpix'][0], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CRPIX2', Cam.tnx['crpix'][1], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CD1_1', Cam.tnx['cd'][0,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD1_2', Cam.tnx['cd'][0,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_1', Cam.tnx['cd'][1,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_2', Cam.tnx['cd'][1,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('xterm', Cam.tnx['xterm'], 'Cross-term type', \
before='EXPTIME')
hdrslow.set('xorder', Cam.tnx['xorder'], 'Polynomial order for x', \
before='EXPTIME')
hdrslow.set('yorder', Cam.tnx['yorder'], 'Polynomial order for y', \
before='EXPTIME')
coefx = ''
coefy = ''
for j in Cam.tnx['coefx']:
coefx += '{} '.format(j)
for k in Cam.tnx['coefy']:
coefy += '{} '.format(k)
hdrslow.set('xcoefs', coefx, 'Polynomial coefficients for x', \
before='EXPTIME')
hdrslow.set('ycoefs', coefy, 'Polynomial coefficients for y', \
before='EXPTIME')
### and saving the binned image
p = multiprocessing.Process(target = mfuncs.mio.saveCompFits, \
args = (PathsI_, name, timt, hdrslow))
try:
p.start()
p.join(25)
finally:
p.terminate()
## Save a thumbnail image in png format of the binned image
mfuncs.mio.savepng(timt, binThumbs, name, binfactor=4)
timt = 0.
hdrslow= 0
skip = False
logger.debug('Generating the global header')
hdr = pf.Header()
hdr['STATION']= (site.name, 'MASCARA station')
hdr['CAMERA'] = (Cam.name, 'Camera identifier')
hdr['ALT0'] = (Cam.Alto, 'General Pointing altitude')
hdr['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
hdr['TH0'] = (Cam.Tho, 'General Pointing orientation')
hdr['X0'] = (Cam.Xo, 'Average Center of the lens -- X')
hdr['Y0'] = (Cam.Yo, 'Average Center of the lens -- Y')
hdr['AVERSION'] = (acqversion, 'Version of the acquisition software')
hdr['RVERSION'] = (softversion, 'Version of the pipeline software')
hdr['CVERSION'] = (coversion, 'Version of the control software')
### added by JS
if 'HDR' in locals() or 'HDR' in globals():
CCDTEMP = HDR['CCDTEMP'] if 'CCDTEMP' in HDR.keys() else HDR['CCD-TEMP']
hdr['EXPTIME'] = (round(HDR['EXPTIME'], 2), 'Measured exposure time in seconds')
hdr['CCDTEMP'] = (CCDTEMP, 'CCD temperature at start of exposure')
else:
hdr['EXPTIME'] = 6.4
hdr['CCDTEMP'] = -15
hdr['NAPER'] = (naper, 'Number of aperture used for the photometry')
for i, r in enumerate(apr):
hdr['APER'+str(i)] = (r, 'Aperture radii in pixel')
hdr['NSTAPER'] = (naper, 'Number of high cadence apertures analysed')
hdr['SKYRAD0'] = (skyrad[0], 'Inner sky radius')
hdr['SKYRAD1'] = (skyrad[1], 'Inner sky radius')
night = Pathlog.split(site.sname)[0][-8:]
saveout.put((Path, night, hdr, 'all'), block=False)
logger.info('end_of_night')
event.clear()
try:
time.sleep(5)
### Get from the saving process the green light:
event.wait(timeout = 3600)
logger.info('saved_all_LC')
except Queue.Empty:
logger.warning('the saving took more than one hour: timeout')
try:
Cam.readAstroCheck(make_plot=True)
masc.light_curve_overview(Path, night, \
filename='LCoverview'+night+site.sname+Cam.sname)
except:
logger.error('No astrometric plot made', exc_info=True)
message, alternate = masc.reduction_summary(Pathlog, site, Cam, night, datetime.now())
masc.send_summary(Pathlog, message, alternate, Cam, site, night)
stacker.end()
handler.close()
ql.stop()
newnight = True
else:
## Now start with the proper reduction
###### ############# #####
### Start of the reduction ###
###### ############# #####
logger.debug('Opening %s ' %filename)
im, HDR = pf.getdata(filename, header=True)
### If there is a masterdark:
if usedarks:
im = im-dark
### Get version infos from the header
acqversion = HDR['AVERSION'] if 'AVERSION' in HDR.keys() \
else 'pre-14.11'
coversion = HDR['CVERSION'] if 'CVERSION' in HDR.keys() \
else 'pre-14.11'
CCDTEMP = HDR['CCDTEMP'] if 'CCDTEMP' in HDR.keys() \
else HDR['CCD-TEMP']
texp = np.around(HDR['EXPTIME'], decimals=3)
### If one image has an exposure time larger than 6.4 sec
### mark it in the log:
if texp > 6.4:
logger.warning('Current exposure time is {}for file {}'.format(texp, filename))
### if the exposure time is really too large: flag the image.
### We don't want to use it
if texp > 6.4*1.5:
BadExposure = True
else:
BadExposure = False
else:
BadExposure = False
date = np.float64(HDR['JD']) + mfuncs.timing.calendar2jd(timedelta(seconds=texp/2))
obstime.append(date)
try:
lstid = HDR['LST-IDX'] if 'LST-IDX' in HDR.keys() else HDR['LPSTIDX']
nolst = False
logger.debug('lstid is %g ' %lstid)
except KeyError:
lstid = 0.
nolst = True
### Create a slow indexing:
if not nolst:
slowind = int(lstid-1) % slowbinning
else:
# If not lst timing, the slowindex has to reach 50 before binning...?
slowind = slowind % (slowbinning)
logger.debug('slowind is %g using LSTtiming ' %slowind)
### Check for any discontinuity in the data taking process
try:
discontinuity = mfuncs.timing.checkdiscontinuity(obstime, texp)
except Exception:
logger.error(' Failed checking for timing discontinuities ', exc_info=True)
#########################################################
##### Start discontinuity scenario : save slow data #####
if discontinuity:
if not quiet:
print 'Found an interruption in the data'
logger.warning('Found an interruption in the data, file {}'.format(filename))
# If saving time, prepare the header
tprepare = datetime.now()
### Calculate stellar position for the slow-binned data,
#### Prepare coordinates for slow cadence curves.
sdate = date + mfuncs.timing.calendar2jd(timedelta(
seconds=(slowbinning/2.-slowind-1)*texp))
slpst = site.local_sidereal_time(sdate)
slpstidx, slpstseq = RFsite.local_sidereal_time(sdate, texp*50, lst=False)
### Average the binned image
timt, Nimages = stacker.binImage(sdate)
stime = (mfuncs.timing.jd2calendar(sdate) - \
timedelta(seconds=texp/2)).strftime('%Y%m%dT%H%M%S')
name = stime + site.sname + Cam.sname
skip = False
### If there are no images in the stacker, do not save a
### binned image...
if Nimages == 0:
logger.info('No image to save, pass ... ')
skip = True
if not skip:
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
xsalt, xsaz, xvis = site.apparent_stars(sra, sdec, sdate, nutation=True, \
precession=True, refraction=True)
tmpmag, tmpbmag, tmpspec = smag[xvis], sbmag[xvis], sSptype[xvis]
### Re-evaluate the corrections to the positions:
Cam.solve_corrections(timt, xsalt[tmpmag<7.], xsaz[tmpmag<7.])
xxpre, yypre, xyinCCD = Cam.visible_stars(xsalt, xsaz, margin=-50)
xxcor, yycor = Cam.correct_positions(xxpre, yypre)
xalt, xaz = xsalt[xyinCCD], xsaz[xyinCCD]
tmpra, tmpdec = sra[xvis], sdec[xvis]
sxra, sxdec = tmpra[xyinCCD], tmpdec[xyinCCD]
tmpmag, tmpbmag, tmpspec = smag[xvis], sbmag[xvis], sSptype[xvis]
vsflagblended = sflagblended[xvis]
vsblendedvalue = sblendedvalue[xvis]
sxmag, sxbmag, sxspect = tmpmag[xyinCCD], tmpbmag[xyinCCD], tmpspec[xyinCCD]
xsflagblended = vsflagblended[xyinCCD]
xsblendedvalue=vsblendedvalue[xyinCCD]
ssID = sID[xvis]
ssID = ssID[xyinCCD]
## Make one big dictionary which is saved at once.
slowStarID = {}
if not quiet:
print ' Running the photometry over %g stars in Binned mode' %xxpre.shape[0]
logger.info(' Running the photometry over ' + \
'{} stars in Binned mode'.format(xxpre.size))
if Transmission:
## Calculating the grid points for the transmission map
lx, ly, upx, upy = mfuncs.utils.skyQuad(timt, tilesize=tilesize, aligned=True)
nulx, nuly = np.unique(lx).size, np.unique(ly).size
QuadNumber = {} ### A dictionary for the transmission map
for w in xrange(xxpre.shape[0]):
# ----------------------------------------------------------------
# The final content of the Photometry array is:
# flag, (flux, erflux)*Naperture, sky, ersky, xc, yc,
# xpre, ypre, altitude, azimuth, CCDtemperature,
# exposure time, JDdate, lstdate, lstid, lstseq
# Total number of elements: (2*naper+6)+6
# ----------------------------------------------------------------
try:
fot = phot.OnaPhot(timt, xxcor[w], yycor[w], apr=sapr, \
skyrad=skyrad, flag=flagit)
except Exception:
logger.error('Failed running the ' + \
' photometry on {}'.format(ssID[w], \
exc_info=True))
continue
if Transmission:
## Fill transmission dictionary with
## observed magnitude- apparent magnitude values
whichQuad = (xxpre[w]>=lx)*(xxpre[w]<upx)*(yypre[w]<upy)*(yypre[w]>=ly)
if whichQuad.nonzero()[0]:
try:
QuadNumber[whichQuad.nonzero()[0][0]].append([-2.5*\
np.log10(fot[1+nap*2])-sxmag[w] for nap in xrange(naper)])
except KeyError:
QuadNumber[whichQuad.nonzero()[0][0]] = \
[[-2.5*np.log10(fot[1+nap*2])-sxmag[w] for nap in xrange(naper)]]
tmp = np.hstack((fot, xalt[w], xaz[w], CCDTEMP, \
texp, sdate, slpst, slpstidx, slpstseq))
slowStarID[ssID[w]] = [date, sxra[w], sxdec[w], \
sxmag[w], sxbmag[w], sxspect[w], \
xsflagblended[w], xsblendedvalue[w], \
np.array(tmp)]
logger.debug('Preparation took %g seconds'%(datetime.now()-tprepare).total_seconds())
logger.debug('Generating the global header')
hdrslow = pf.Header()
hdrslow['STATION']= (site.name, 'MASCARA station')
hdrslow['CAMERA'] = (Cam.name, 'Camera identifier')
hdrslow['ALT0'] = (Cam.Alto, 'General Pointing altitude')
hdrslow['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
hdrslow['TH0'] = (Cam.Tho, 'General Pointing orientation')
hdrslow['X0'] = (Cam.Xo, 'Average Center of the lens -- X')
hdrslow['Y0'] = (Cam.Yo, 'Average Center of the lens -- Y')
hdrslow['EXPTIME'] = (stacker.texp, 'Mean exposure time in seconds')
hdrslow['CCDTEMP'] = (CCDTEMP, 'CCD temperature at start of exposure')
hdrslow['AVERSION'] = (acqversion, 'Version of the acquisition software')
hdrslow['RVERSION'] = (softversion, 'Version of the pipeline software')
hdrslow['CVERSION'] = (coversion, 'Version of the control software')
hdrslow['NEXP'] = (Nimages, 'Number of images binned together')
hdrslow['NAPER'] = (nsaper, 'Number of aperture used for the photometry')
for i, r in np.ndenumerate(sapr):
hdrslow['APER'+str(i[0])] = (r, 'Aperture radii in pixel')
for j in xrange(naper):
hdrslow['STAPER'+str(j)] = ('XXX', 'Statistics of the high cadence aperture')
hdrslow['NSTAPER'] = (naper, 'Number of high cadence apertures analysed')
hdrslow['SKYRAD0'] = (skyrad[0], 'Inner sky radius')
hdrslow['SKYRAD1'] = (skyrad[1], 'Inner sky radius')
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
#### Now save the faint limit dictionary:
sDBNiObs +=1
logger.info('Save slow dictionary in database {}'.format(sDBNiObs))
try:
masc.tmpdbSave(tmpPath, slowStarID, sDBNiObs, qlog)
except Exception:
logger.error('Failed saving the slow database. ', exc_info=True)
saveout.put((Path, sDBNiObs, hdrslow, ''), block=False)
if Transmission:
### Create the transmission map
logger.info('Create transmission map(s)...')
masc.buildTransmissionMap(QuadNumber, nulx, \
nuly, hdrslow, binThumbs, name+'_trans', \
writeto=True, dimension=naper)
##### Saving the stacked image!!!
hdrslow['LSTMID'] = (slpst, 'LST at mid-bin observation')
hdrslow['JDMID'] = (sdate, 'JD at mid-bin observation')
hdrslow.set('CRPIX1', Cam.tnx['crpix'][0], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CRPIX2', Cam.tnx['crpix'][1], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CD1_1', Cam.tnx['cd'][0,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD1_2', Cam.tnx['cd'][0,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_1', Cam.tnx['cd'][1,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_2', Cam.tnx['cd'][1,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('xterm', Cam.tnx['xterm'], 'Cross-term type', \
before='EXPTIME')
hdrslow.set('xorder', Cam.tnx['xorder'], 'Polynomial order for x', \
before='EXPTIME')
hdrslow.set('yorder', Cam.tnx['yorder'], 'Polynomial order for y', \
before='EXPTIME')
coefx = ''
coefy = ''
for j in Cam.tnx['coefx']:
coefx += '{} '.format(j)
for k in Cam.tnx['coefy']:
coefy += '{} '.format(k)
hdrslow.set('xcoefs', coefx, 'Polynomial coefficients for x', \
before='EXPTIME')
hdrslow.set('ycoefs', coefy, 'Polynomial coefficients for y', \
before='EXPTIME')
p = multiprocessing.Process(target = mfuncs.mio.saveCompFits, \
args = (PathsI_, name, timt, hdrslow))
try:
p.start()
p.join(25)
finally:
p.terminate()
## Save a thumbnail image in png format of the binned image
mfuncs.mio.savepng(timt, binThumbs, name, binfactor=4)
timt = 0.
skip = False
########## End of Discontinuity case #######################
############################################################
########## Continue normal procedure #######################
valt, vaz, vis = site.apparent_stars(ra, dec, date, refraction=True, \
precession=True, nutation=True)
## Discard here the stars below horizon
vmag, vbmag, vID = mag[vis], bmag[vis], ID[vis]
vra, vdec, vspect = ra[vis], dec[vis], Sptype[vis]
if discontinuity or (slowind ==0) or firstAstrometry:
valt_, vaz_, vmag_ = valt[vmag<7.], vaz[vmag<7.], vmag[vmag<7.]
if firstAstrometry or discontinuity:
highoffset = Cam.solve_corrections(im, valt_, vaz_, qlog, debug=False)
### Special case to recalculate the astrometric master solution if:
## - the offset is important
## - the image has no bright elements (Moon) or is not too bright (clouds)
## Additional safeguards are put in the solve_astrometry routine
## to exit if not enough stars are found (clouds again)
if usedarks:
brightimage = np.median(im) > np.median(dark)*0.5
else:
brightimage = True
if not brightimage and highoffset:
try:
Cam.solve_astrometry(im, valt_, vaz_, vmag_, qlog)
except Exception:
logger.error('Did not solve the master solution properly. \n'+ \
'Reload previous one', exc_info=True)
Cam.reload()
logger.info('solved_the_astrometry')
## Reset the newnight boelan
firstAstrometry=False
logger.info(' Prepare for the astrometry')
### Compute the low order distortions which are added up
### to the master solution
try:
Cam.solve_corrections(im, valt_, vaz_, qlog, debug=False)
except Exception:
logger.error('Failed running the astrometry', exc_info=True)
gc.collect()
### Save the current solution in a log
### Plus do a quick rough check of the solution
try:
ckalt, ckaz = Cam.logAstrometryCheck(date, lstid, im.shape, \
out=True, fileroot='astrocheck', logdir=Pathlog)
if checkalt.shape[0]==0:
checkalt = np.hstack((ckalt))
checkaz = np.hstack((ckaz))
else:
checkalt = np.vstack((checkalt, ckalt))
checkaz = np.vstack((checkaz, ckaz))
if (np.abs(ckalt-checkalt.mean(axis=0))>5.*checkalt.std(axis=0)).any():
flagit = 4
except Exception:
logger.error('failed ', exc_info=True)
## Do the image rotation centered on the middle image of a sequence
## of 50 => the 24.5th?
ttt = date + mfuncs.timing.calendar2jd(timedelta(
seconds=(slowbinning/2.-slowind-1)*texp))
stacker.image_rotate(im, date, ttt, exposuretime = texp)
lpsttime = site.local_sidereal_time(date)
lpstseq = RFsite.local_sidereal_time(date, texp, lstindex=False, lst=False)
### Calculate the position of the stars on CCD
xpre, ypre, inCCD = Cam.visible_stars(valt, vaz, margin=-50)
xcor, ycor = Cam.correct_positions(xpre, ypre)
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
xalt, xaz = valt[inCCD], vaz[inCCD]
xra, xdec = vra[inCCD], vdec[inCCD]
xspect, xmag, xbmag = vspect[inCCD], vmag[inCCD], vbmag[inCCD]
vflagblended = flagblended[vis]
xflagblended = vflagblended[inCCD]
vblendedvalue = blendedvalue[vis]
xblendedvalue = vblendedvalue[inCCD]
StarVis = vID[inCCD]
### The high cadence cadence part; Only photometry ###
logger.info('Running the photometry over the %g stars' %xpre.shape[0])
for j in xrange(xpre.shape[0]):
# Case 1: the star was already in the CCd and remains inside:
if StarVis[j] in StarID:
d = StarID[StarVis[j]]
Phot = d[8]
try:
fot = phot.OnaPhot(im, xcor[j], ycor[j], apr=apr, \
skyrad=skyrad, flag=flagit)
except Exception:
logger.error('Failed running the photometry on %s'%StarVis[j], \
exc_info=True)
continue
# ----------------------------------------------------------------
# The final content of the Photometry array is:
# flag, (flux, erflux)*Naperture, sky, ersky, peak, xc, yc,
# altitude, azimuth, CCD temperature,
# exposure time, JDdate, lstdate, lpstid, lpstidx
# Total number of elements: (2*naper+6)+8
# ----------------------------------------------------------------
tmp = np.hstack((fot, xalt[j], xaz[j], CCDTEMP, \
texp, date, lpsttime, lstid, lpstseq))
if len(Phot) < 1:
Phot = np.reshape(tmp, (1, (2*naper+6)+8))
else:
Phot = np.vstack((Phot, tmp))
StarID[StarVis[j]]=[d[0], d[1], d[2], d[3], d[4], \
d[5], d[6], d[7], Phot]
### Case 2: the star was not observed before, but is now on the CCD
else:
#----------------------------------------------------------------------
# The content of the dictionary is:
# start of observation, ra, dec, mag, bmag, spectral type,
# Photometric array(flag, flux, erflux, sky, ersky, xc, yc)
# xposition, yposition, altitude, azimuth,
# exposure time, jd date, lstdate, lpstid, lpstidx
#-----------------------------------------------------------------------
try:
fot = phot.OnaPhot(im, xcor[j], ycor[j], apr=apr, \
skyrad=skyrad, flag=flagit)
except Exception:
logger.error('Failed running OnaPhot on %s' %StarVis[j], exc_info=True)
continue
tmp = np.hstack((fot, xalt[j], xaz[j], CCDTEMP, \
texp, date, lpsttime, lstid, lpstseq))
StarID[StarVis[j]] = [date, xra[j], xdec[j], \
xmag[j], xbmag[j], xspect[j], \
xflagblended[j], xblendedvalue[j], tmp]
#####------------------------------------------#####
##### SAVING THE PHOTOMETRY IN TEMPORARY FILES #####
if (slowind + 1) % slowbinning == 0:
# If saving time, prepare the header
tprepare = datetime.now()
logger.debug('Generating the global header')
hdr = pf.Header()
hdr['STATION']= (site.name, 'MASCARA station')
hdr['CAMERA'] = (Cam.name, 'Camera identifier')
hdr['ALT0'] = (Cam.Alto, 'General Pointing altitude')
hdr['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
hdr['TH0'] = (Cam.Tho, 'General Pointing orientation')
hdr['X0'] = (Cam.Xo, 'Average Center of the lens -- X')
hdr['Y0'] = (Cam.Yo, 'Average Center of the lens -- Y')
hdr['EXPTIME'] = (HDR['EXPTIME'], 'Measured exposure time in seconds')
hdr['CCDTEMP'] = (CCDTEMP, 'CCD temperature at start of exposure')
hdr['AVERSION'] = (acqversion, 'Version of the acquisition software')
hdr['RVERSION'] = (softversion, 'Version of the pipeline software')
hdr['CVERSION'] = (coversion, 'Version of the control software')
hdr['NAPER'] = (naper, 'Number of apertures used for the photometry')
for i, r in np.ndenumerate(apr):
hdr['APER'+str(i[0])] = (r, 'Aperture radii in pixel')
hdr['SKYRAD0'] = (skyrad[0], 'Inner Sky radius')
hdr['SKYRAD1'] = (skyrad[1], 'Outer sky radius')
logger.debug('Preparing the dictionary for the saving sequence...')
ttt = date+ mfuncs.timing.calendar2jd( timedelta( \
seconds=(slowbinning/2.-slowind -1)*texp))
slpst = site.local_sidereal_time(ttt)
slpstidx, slpstseq = RFsite.local_sidereal_time(ttt, texp*50, lst=False)
### Average the binned image
#timt = np.float32(timt)/images_in_timt
timt, Nimages = stacker.binImage(ttt)
if not quiet:
print ' Created binned image out of %g images'%Nimages
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
### Calculate stellar position for the slow-binned data,
#### Prepare coordinates for slow cadence curves.
xsalt, xsaz, xvis = site.apparent_stars(sra, sdec, ttt, nutation=True, \
precession=True, refraction=True)
tmpmag, tmpbmag, tmpspec = smag[xvis], sbmag[xvis], sSptype[xvis]
### Re-evaluate the corrections to the positions:
Cam.solve_corrections(timt, xsalt[tmpmag<7.], xsaz[tmpmag<7.])
xxpre, yypre, xyinCCD = Cam.visible_stars(xsalt, xsaz, margin=-50)
xxcor, yycor = Cam.correct_positions(xxpre, yypre)
xalt, xaz = xsalt[xyinCCD], xsaz[xyinCCD]
tmpra, tmpdec = sra[xvis], sdec[xvis]
sxra, sxdec = tmpra[xyinCCD], tmpdec[xyinCCD]
vsflagblended = sflagblended[xvis]
vsblendedvalue = sblendedvalue[xvis]
sxmag, sxbmag, sxspect = tmpmag[xyinCCD], tmpbmag[xyinCCD], tmpspec[xyinCCD]
xsflagblended = vsflagblended[xyinCCD]
xsblendedvalue=vsblendedvalue[xyinCCD]
ssID = sID[xvis]
ssID = ssID[xyinCCD]
## Make one big dictionary which is saved at once.
slowStarID = {}
if not quiet:
print ' Running the photometry over %g stars in Binned mode' %xxpre.shape[0]
logger.info(' Running the photometry over %g stars in Binned mode' %xxpre.shape[0])
if Transmission:
## Calculating the grid points for the transmission map
lx, ly, upx, upy = mfuncs.utils.skyQuad(timt, tilesize=tilesize, aligned=True)
nulx, nuly = np.unique(lx).size, np.unique(ly).size
QuadNumber = {} ### A dictionary for the transmission map
for w in xrange(xxpre.shape[0]):
# ----------------------------------------------------------------
# The final content of the Photometry array is:
# flag, (flux, erflux)*Naperture, sky, ersky, peak, xc, yc,
# altitude, azimuth, CCD Temperature,
# exposure time, JDdate, lstdate, lstid, lstseq
# Total number of elements: (2*naper+6)+
# ----------------------------------------------------------------
try:
fot = phot.OnaPhot(timt, xxcor[w], yycor[w], apr=sapr, \
skyrad=skyrad, flag=flagit)
except Exception:
logger.error('Failed running the photometry on %s'%ssID[w], \
exc_info=True)
continue
if Transmission:
# Fill transmission dictionary with
# observed magnitude- apparent magnitude values
whichQuad = (xxpre[w]>=lx)*(xxpre[w]<upx)*(yypre[w]<upy)*(yypre[w]>=ly)
if whichQuad.nonzero()[0]:
try:
QuadNumber[whichQuad.nonzero()[0][0]].append([-2.5*\
np.log10(fot[1+nap*2])-sxmag[w] \
for nap in xrange(nsaper)])
except KeyError:
QuadNumber[whichQuad.nonzero()[0][0]] = \
[[-2.5*np.log10(fot[1+nap*2])-sxmag[w] \
for nap in xrange(nsaper)]]
tmp = np.hstack((fot, xalt[w], xaz[w], CCDTEMP, texp, \
ttt, slpst, int(slpstidx), slpstseq))
slowStarID[ssID[w]] = [date, sxra[w], sxdec[w], \
sxmag[w], sxbmag[w], sxspect[w], \
xsflagblended[w], xsblendedvalue[w],\
np.array(tmp)]
logger.debug('Preparation took %g seconds'%(datetime.now()-\
tprepare).total_seconds())
logger.info('Save high cadence dictionary in temporary database ...')
t2 = datetime.now()
try:
masc.tmpdbSave(tmpPath, StarID, DBNiObs, qlog)
except Exception:
logger.error('Failed saving the dictionary. FATAL!!!', exc_info=True)
break
#### Now save the faint limit dictionary:
logger.info('Save slow dictionary in database ...')
sDBNiObs +=1
logger.debug('Generating the global header')
hdrslow = pf.Header()
hdrslow['STATION']= (site.name, 'MASCARA station')
hdrslow['CAMERA'] = (Cam.name, 'Camera identifier')
hdrslow['ALT0'] = (Cam.Alto, 'General Pointing altitude')
hdrslow['AZ0'] = (Cam.Azo, 'General Pointing azimuth')
hdrslow['TH0'] = (Cam.Tho, 'General Pointing orientation')
hdrslow['X0'] = (Cam.Xo, 'Average Center of the lens -- X')
hdrslow['Y0'] = (Cam.Yo, 'Average Center of the lens -- Y')
hdrslow['EXPTIME'] = (stacker.texp, 'Measured exposure time in seconds')
hdrslow['CCDTEMP'] = (CCDTEMP, 'CCD temperature at start of exposure')
hdrslow['AVERSION'] = (acqversion, 'Version of the acquisition software')
hdrslow['RVERSION'] = (softversion, 'Version of the pipeline software')
hdrslow['CVERSION'] = (coversion, 'Version of the control software')
hdrslow['NEXP'] = (Nimages, 'Number of images binned together')
stime = (mfuncs.timing.jd2calendar(ttt) - \
timedelta(seconds=texp/2)).strftime('%Y%m%dT%H%M%S')
name = stime + site.sname + Cam.sname
hdrslow['NAPER'] = (nsaper, 'Number of aperture used for the photometry')
for i, r in np.ndenumerate(sapr):
hdrslow['APER'+str(i[0])] = (r, 'Aperture radii in pixel')
for j in xrange(naper):
hdrslow['STAPER'+str(j)] = ('XXX', 'Statistics of the high cadence aperture')
hdrslow['NSTAPER'] = (naper, 'Number of high cadence apertures analysed')
hdrslow['SKYRAD0'] = (skyrad[0], 'Inner sky radius')
hdrslow['SKYRAD1'] = (skyrad[1], 'Inner sky radius')
### Part to skip when not performing the photometry
if not noPhot and not BadExposure:
try:
masc.tmpdbSave(tmpPath, slowStarID, sDBNiObs, qlog)
except Exception:
logger.error('Failed saving the slow database. ', exc_info=True)
logger.debug('Complete saving done in %g' %(datetime.now()-t2).total_seconds())
## Reset the reference dictionary
logger.debug('Reset the reference dictionary...')
t4 = datetime.now()
try:
ocStarID = StarID.keys()
for k in ocStarID:
if k in StarVis:
d = StarID[k]
Phot = np.array([])
StarID.update({k:[d[0], d[1], d[2], d[3], d[4], \
d[5], d[6], d[7], Phot]})
else:
del StarID[k]
except Exception:
logger.error('Failed reseting the dictionary', exc_info=True)
## clean up unused variables:
slowStarID = None
salt, saz = None, None
xsalt, xsaz = 0, 0
xxpre, yypre, xyinCCD = 0, 0, 0
xalt, xaz = 0, 0
tmpra, tmpdec = 0, 0
sxra, sxdec = 0, 0
tmpmag, tmpbmag, tmpspec = 0, 0, 0
sxmag, sxbmag, sxspect = 0, 0, 0
ssID = None
logger.debug('Dictionary reseted in %g' %(datetime.now()-t4).total_seconds())
saveout.put((Path, DBNiObs, hdr,''), block=False)
saveout.put((Path, sDBNiObs, hdrslow,''), block=False)
DBNiObs +=1
gc.collect()
if Transmission:
##### Create transmission map(s) ####
logger.info('Create transmission map(s)...')
g = multiprocessing.Process(target=masc.buildTransmissionMap, \
args=(QuadNumber, nulx, nuly, hdrslow, \
binThumbs, name+'_trans'), \
kwargs={'writeto':True, 'dimension':nsaper})
try:
g.start()
g.join(3)
finally:
g.terminate()
##### Saving the 50stacked image!!!
hdrslow['BINNING'] = (slowbinning, 'Stacked images')
hdrslow['LSTMID'] = (slpst, 'LST at mid-bin observation')
hdrslow['JDMID'] = (ttt, 'JD at mid-bin observation')
hdrslow.set('CRPIX1', Cam.tnx['crpix'][0], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CRPIX2', Cam.tnx['crpix'][1], 'Coordinate reference pixel', \
before='EXPTIME')
hdrslow.set('CD1_1', Cam.tnx['cd'][0,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD1_2', Cam.tnx['cd'][0,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_1', Cam.tnx['cd'][1,0], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('CD2_2', Cam.tnx['cd'][1,1], 'Coordinate matrix', \
before='EXPTIME')
hdrslow.set('xterm', Cam.tnx['xterm'], 'Cross-term type', \
before='EXPTIME')
hdrslow.set('xorder', Cam.tnx['xorder'], 'Polynomial order for x', \
before='EXPTIME')
hdrslow.set('yorder', Cam.tnx['yorder'], 'Polynomial order for y', \
before='EXPTIME')
coefx = ''
coefy = ''
for j in Cam.tnx['coefx']:
coefx += '{} '.format(j)
for k in Cam.tnx['coefy']:
coefy += '{} '.format(k)
hdrslow.set('xcoefs', coefx, 'Polynomial coefficients for x', \
before='EXPTIME')
hdrslow.set('ycoefs', coefy, 'Polynomial coefficients for y', \
before='EXPTIME')
p = multiprocessing.Process(target = mfuncs.mio.saveCompFits, \
args = (PathsI_, name, timt, hdrslow))
try:
p.start()
p.join(25)
finally:
p.terminate()
p.join()
print 'Out ...................... managed saving the binned image'
## Save a thumbnail image in png format of the binned image
mfuncs.mio.savepng(timt, binThumbs, name, binfactor=4)
timt = 0.
gc.collect()
logger.info('Time needed to process one image:'+ \
' {}'.format((datetime.now() - tinit).total_seconds()))
if nolst:
slowind+=1
