1
2 """
3 Interface the spectra from the Hermes spectrograph.
4
5 The most important function is L{search}. This looks in SIMBAD for the coordinates
6 of a given object, and finds all spectra matching those within a given radius.
7 If the object's name is not recognised, it will look for correspondence between
8 the given name and the contents of the FITS header's keyword C{object}. L{search}
9 returns a record array, so that you can easily access all columns by their names.
10
11 Note that Hermes spectra retrieved in B{logscale} should B{not be corrected} for
12 the B{barycentric velocity} (the pipeline does it). The spectra retrieved in B{normal
13 wavelength scale should} be corrected. These two cases are outlined below.
14
15 Section 1. Data lookup and reading
16 ==================================
17
18 B{Example usage:} retrieve all data on HD170580
19
20 >>> mydata = search('HD170580')
21
22 Keep only those with a long enough exposure time:
23
24 >>> myselection = mydata[mydata['exptime']>500]
25
26 Now read in all the data, and plot the spectra. First, we need an extra module
27 to read the FITS file and the plotting package.
28
29 >>> from ivs.inout import fits
30 >>> import pylab as pl
31
32 Then we can easily plot the relevant data:
33
34 >>> for fname in myselection['filename']:
35 ... wave,flux = fits.read_spectrum(fname)
36 ... p = pl.plot(wave,flux)
37 >>> p = pl.ylim(0,45000)
38 >>> p = pl.xlim(4511,4513.5)
39
40 ]include figure]]ivs_catalogs_hermes_HD170580.png]
41
42 Note that you can easily shift them according to some radial velocity: as an
43 extra example, we retrieve the data in wavelength scale, shift according to the
44 barycentric velocity, convert to velocity space, and plot them:
45
46 First start a new figure and add extra modules:
47
48 >>> p = pl.figure()
49 >>> from ivs.spectra.tools import doppler_shift # to apply doppler shift
50 >>> from ivs.units import conversions # to convert to velocity space
51
52 Then get the spectra again, but now not in log scale. Make the same selection
53 on exposure time.
54
55 >>> mydata = search('HD170580',data_type='cosmicsremoved_wavelength')
56 >>> myselection = mydata[mydata['exptime']>500]
57
58 Then read them all in and shift according to the barycentric velocity. Also,
59 convert the spectra to velocity space afterwards for plotting purposes.
60
61 >>> for rv,fname in zip(myselection['bvcor'],myselection['filename']):
62 ... wave,flux = fits.read_spectrum(fname)
63 ... wave_shifted = model.doppler_shift(wave,rv)
64 ... velo_shifted = conversions.convert('angstrom','km/s',wave_shifted,wave=(4512.3,'angstrom'))
65 ... p = pl.plot(velo_shifted,flux)
66 >>> p = pl.ylim(0,45000)
67 >>> p = pl.xlim(-70,70)
68
69 ]include figure]]ivs_catalogs_hermes_HD170580_velo.png]
70
71 Section 2. Extracting information
72 =================================
73
74 If you want to list the observers of a certain target, and the amount of spectra
75 they took, you can do the following:
76
77 >>> data = search('HD50230')
78 >>> print [(i,list(data['observer']).count(i)) for i in set(data['observer'])]
79 [('Robin Lombaert', 4), ('Steven Bloemen', 6), ('Pieter Degroote', 25), ('Michel Hillen', 3)]
80
81 Section 3. Radial velocity computation with the HERMES DRS
82 ==========================================================
83
84 Make sure you have sourced the Hermes DRS (C{source ~mercator/HermesDRS.rc}) and
85 that you make a backup of your config file before proceeding.
86
87 In this example, we first make a mask file:
88
89 >>> from ivs.spectra import linelists
90 >>> teff = 12500
91 >>> logg = 4.0
92 >>> ll = linelists.get_lines(teff,logg)
93 >>> mask_file = '%.0f_%.2f.fits'%(teff,logg)
94 >>> make_mask_file(ll['wavelength'],ll['depth'],filename=mask_file)
95
96 And run the DRS:
97
98 >>> CCFList('HD170200',mask_file=mask_file)
99
100
101 Section 4. Hermes overview file
102 ===============================
103
104 To ensure a fast lookup of datafiles, an overview file C{HermesFullDataOverview.tsv}
105 is created via L{make_data_overview}. The file resides in one of the C{IVSdata}
106 directories, and there should also exist a copy at C{/STER/mercator/hermes/}.
107
108 The best strategy to keep the file up-to-date is by running this module in the
109 background in a terminal, via (probably on pleiad22)::
110
111 $:> python hermes.py update
112
113 This script will look for a C{HermesFullDataOverview.tsv} in one of the data
114 directories (see L{config} module), check until what date data files are stored,
115 and add entries for HERMES data files created since the last update. If you
116 want a complete update of all the HERMES data folders, you will need to remove
117 the file and run the update. Indeed, if C{HermesFullDataOverview.tsv} does not
118 exist, it will create a new file in the current working directory from scratch,
119 and copy it to C{/STER/mercator/hermes/}. It is the user's responsibility to
120 copy the file also to one of the IVSdata directories where the user has write
121 permission, so that next time the update is performed, the script can start from
122 the previous results.
123
124 When running the above command in the terminal, you will notice that the script
125 does not terminate until a manual C{CTRL+C} is performed. Indeed, when left running
126 the script will perform the update once a day. So once it runs, as long as the
127 filepaths do not change or computers do not shut down, you don't need to run it
128 again. If a computer does stop running, just restart again with::
129
130 $:> python hermes.py update
131
132 and everything should be fine.
133
134 @var hermesDir: Path to the directory in which the hermes folders (hermesAnalysis,
135 hermesRun, hermesDebug) are located. By default this is set to
136 '/home/user/'.
137
138 @var tempDir: Path to temporary directory where files nessessary for hermesVR to run
139 can be copied to. You need write permission in this directory. The
140 default is set to '/scratch/user/'
141 """
142 import re
143 import os
144 import sys
145 import glob
146 import time
147 import shutil
148 import logging
149 import getpass
150 import datetime
151 import subprocess
152 import numpy as np
153 import astropy.io.fits as pf
154
155 import copy
156 from lxml import etree
157 from xml.etree import ElementTree as ET
158 from collections import defaultdict
159
160 from ivs.catalogs import sesame
161 from ivs.inout import ascii, fits
162 from ivs.aux import loggers
163 from ivs.observations import airmass
164 from ivs.observations.barycentric_correction import helcorr
165 from ivs.units import conversions
166 from ivs import config
167 from ivs.spectra import tools as sptools
168
169
170 hermesDir = os.path.expanduser('~/')
171 tempDir = '/scratch/%s/'%(getpass.getuser())
172
173 logger = logging.getLogger("CAT.HERMES")
174 logger.addHandler(loggers.NullHandler())
175
176
177
178 -def search(ID=None,time_range=None,prog_ID=None,data_type='cosmicsremoved_log',
179 radius=1.,filename=None):
180 """
181 Retrieve datafiles from the Hermes catalogue.
182
183 B{If C{ID} is given}: A string search is performed to match the 'object'
184 field in the FITS headers. The coordinates are pulled from SIMBAD. If the
185 star ID is recognised by SIMBAD, an additional search is done based only on
186 the coordinates. The union of both searches is the final result.
187
188 B{If C{time_range} is given}: The search is confined within the defined
189 range. If you only give one day, the search is confined to the observations
190 made during the night starting at that day. If C{ID} is not given, all
191 observations will be returned of the given datatype.
192
193 B{If C{prog_ID} is given}: The search is performed to match the number of
194 the program. Individual stars are not queried in SIMBAD, so any information
195 that is missing in the header will not be corrected.
196
197 If you don't give either ID or time_range, the info on all data will be
198 returned. This is a huge amount of data, so it can take a while before it
199 is returned. Remember that the header of each spectrum is read in and checked.
200
201 Data type can be any of:
202 1. cosmicsremoved_log: return log merged without cosmics
203 2. cosmicsremoved_wavelength: return wavelength merged without cosmics
204 3. ext_log: return log merged with cosmics
205 4. ext_wavelength: return wavelength merged with cosmics
206 5. raw: raw files (also TECH..., i.e. any file in the raw directory)
207
208 This functions needs a C{HermesFullDataOverview.tsv} file located in one
209 of the datadirectories from C{config.py}, and subdirectory C{catalogs/hermes}.
210
211 If this file does not exist, you can create it with L{make_data_overview}.
212
213 If you want a summary file with the data you search for, you can give
214 C{filename} as an extra keyword argument. The results will be saved to that
215 file.
216
217 The columns in the returned record array are listed in L{make_data_overview},
218 but are repeated here (capital letters are directly retrieved from the
219 fits header, small letters are calculated values. The real header strings
220 are all small capitals):
221
222 1. UNSEQ
223 2. PROG_ID
224 3. OBSMODE
225 4. BVCOR
226 5. OBSERVER
227 6. OBJECT
228 7. RA
229 8. DEC
230 9. BJD
231 10. EXPTIME
232 11. PMTOTAL
233 12. DATE-AVG
234 13. OBJECT
235 14. airmass
236 15. filename
237
238 The column C{filename} contains a string with the absolute location of the
239 file. If you need any extra information from the header, you can easily
240 retrieve it.
241
242 If BVCOR or BJD are not available from the FITS header, this function will
243 attempt to calculate it. It will not succeed if the object's name is not
244 recognised by SIMBAD.
245
246 Example usage: retrieve all data on HD50230
247
248 >>> mydata = search('HD50230')
249
250 Keep only those with a long enough exposure time:
251
252 >>> myselection = mydata[mydata['exptime']>500]
253
254 Look up the 'telalt' value in the FITS headers of all these files via a fast
255 list comprehension:
256
257 >>> telalts = [pf.getheader(fname)['telalt'] for fname in myselection['filename']]
258
259 Search for all data of HD50230 taken in the night of 22 September 2009:
260
261 >>> data = hermes.search('HD50230',time_range='2009-9-22')
262
263 Or within an interval of a few days:
264
265 >>> data = hermes.search('HD50230',time_range=('2009-9-23','2009-9-30'))
266
267 Search for all data observed in a given night:
268
269 >>> data = hermes.search(time_range='2009-9-22')
270
271 B{Warning:} the heliocentric correction is not calculated when no ID is given,
272 so make sure it is present in the header if you need it, or calculate it yourself.
273
274 @param ID: ID of the star, understandable by SIMBAD
275 @type ID: str
276 @param time_range: range of dates to confine the search to
277 @type time_range: tuple strings of type '2009-09-23T04:24:35.712556' or '2009-09-23'
278 @param data_type: if None, all data will be returned. Otherwise, subset
279 'cosmicsremoved', 'merged' or 'raw'
280 @type data_type: str
281 @param radius: search radius around the coordinates (arcminutes)
282 @type radius: float
283 @param filename: write summary to outputfile if not None
284 @type filename: str
285 @return: record array with summary information on the observations, as well
286 as their location (column 'filename')
287 @rtype: numpy rec array
288 """
289
290
291 ctlFile = '/STER/mercator/hermes/HermesFullDataOverview.tsv'
292 data = ascii.read2recarray(ctlFile, splitchar='\t')
293
294 keep = np.array(np.ones(len(data)),bool)
295
296 if time_range is not None:
297 if isinstance(time_range,str):
298 time_range = _timestamp2datetime(time_range)
299 time_range = (time_range,time_range+datetime.timedelta(days=1))
300 else:
301 time_range = (_timestamp2datetime(time_range[0]),_timestamp2datetime(time_range[1]))
302 keep = keep & np.array([(time_range[0]<=_timestamp2datetime(i)<=time_range[1]) for i in data['date-avg']],bool)
303 info = None
304
305
306
307 if ID is not None:
308 info = sesame.search(ID)
309
310
311 ID = ID.replace(' ','').replace('.','').replace('+','').replace('-','').replace('*','')
312 match_names = np.array([objectn.replace(' ','').replace('.','').replace('+','').replace('-','').replace('*','') for objectn in data['object']],str)
313 keep_id = [((((ID in objectn) or (objectn in ID)) and len(objectn)) and True or False) for objectn in match_names]
314 keep_id = np.array(keep_id)
315
316 if info:
317 ra,dec = info['jradeg'],info['jdedeg']
318 keep_id = keep_id | (np.sqrt((data['ra']-ra)**2 + (data['dec']-dec)**2) < radius/60.)
319 keep = keep & keep_id
320
321 if prog_ID is not None:
322 keep = keep & (data['prog_id']==prog_ID)
323
324
325
326 if np.any(keep):
327 data = data[keep]
328
329 if data_type is not None:
330 data_type == data_type.lower()
331
332
333 if not data_type=='raw':
334 data['filename'] = [_derive_filelocation_from_raw(ff,data_type) for ff in data['filename']]
335 existing_files = np.array([ff!='naf' for ff in data['filename']],bool)
336 data = data[existing_files]
337 seqs = sorted(set(data['unseq']))
338 logger.info('ID={}/prog_ID={}: Found {:d} spectra (data type={} with unique unseqs)'.format(ID,prog_ID,len(seqs),data_type))
339 else:
340 data = data[:0]
341 logger.info('%s: Found no spectra'%(ID))
342
343
344
345
346
347 for obs in data:
348 if ID is not None and info:
349 try:
350 jd = _timestamp2jd(obs['date-avg'])
351 except ValueError:
352 logger.info('Header probably corrupted for unseq {}: no info on time or barycentric correction'.format(obs['unseq']))
353 jd = np.nan
354
355
356
357 bvcorr, hjd = helcorr(ra/360.*24, dec, jd)
358 else:
359 break
360 if np.isnan(obs['bvcor']):
361 logger.info("Corrected 'bvcor' for unseq {} (missing in header)".format(obs['unseq']))
362 obs['bvcor'] = float(bvcorr)
363 if np.isnan(obs['bjd']):
364 logger.info("Corrected 'bjd' for unseq {} (missing in header)".format(obs['unseq']))
365 obs['bjd'] = float(hjd)
366
367
368
369 if filename is not None:
370 ascii.write_array(data,filename,auto_width=True,header=True)
371 else:
372 return data
373
375 """
376 Combines HERMES spectra with sigma clipping to remove cosmics. The input spectra
377 are given as a list of filenames for the objects and for the wavelength scales.
378 The output is the wavelength scale of the first object, the merged and sigma clipped
379 flux, and the adapted header of the object. In this adapted header, the exposure
380 time is the sum of all individual exposure times, and the observing time is averaged.
381 Uses the L{ivs.spectra.tools.merge_cosmic_clipping} method to merge the spectra.
382
383 >>> data = search('KIC9540226')
384 >>> objlist = data['filename']
385 >>> wave, flux, header = merge_hermes_spectra(objlist, wscalelist=None)
386
387 @param objlist: list of OBJ or order merged filenames
388 @param wscalelist: list of wavelengthscale filenames
389 @param vrads: list of radial velocities (optional)
390 @param vrad_units: units of the radial velocities
391 @param sigma: value used for sigma clipping
392 @param window: window size used in median filter
393 @param runs: number of iterations through the spectra
394
395 @return: The combined spectrum, cosmic clipped. (wave, flux, header)
396 @rtype: [array, array, dict]
397 """
398 kwargs['full_output'] = False
399
400 if wscalelist == None:
401
402 exptime, bjd = 0, np.zeros(len(objlist))
403 waves, fluxes = [],[]
404 for i, ofile in enumerate(objlist):
405 w, f, h = fits.read_spectrum(ofile, return_header=True)
406 waves.append(w)
407 fluxes.append(f)
408 exptime += h['exptime']
409
410
411 waves, fluxes = np.array(waves), np.array(fluxes)
412 mwave, mflux = sptools.merge_cosmic_clipping(waves, fluxes, **kwargs)
413
414 else:
415
416
417 exptime, bjd = 0, np.zeros(len(objlist))
418 waves, fluxes = np.zeros((55,4608, len(objlist))), np.zeros((55,4608, len(wscalelist)))
419 for i, (ofile, wfile) in enumerate(zip(objlist, wscalelist)):
420 odata = pf.getdata(ofile, 0).T
421 oheader = pf.getheader(ofile, 0)
422 wdata = pf.getdata(wfile, 0)
423 fluxes[:,:,i] = odata
424 waves[:,:,i] = wdata
425 exptime += oheader['exptime']
426
427
428
429 mwave, mflux = waves[:,:,0], np.zeros((55, 4608))
430 for i in range(55):
431 wave, flux = sptools.merge_cosmic_clipping(waves[i,:,:].T, fluxes[i,:,:].T, **kwargs)
432 mflux[i] = flux
433
434 header = pf.getheader(objlist[0], 0)
435 header['exptime'] = exptime
436
437
438 return mwave, mflux, header
439
441 """
442 Mimics HermesTool MakeListStar.
443
444 This should work as input for HermesTool CCFList.py
445
446 The result is a file in the current working directory with name C{ID.list}.
447 If you have specified the C{direc} keyword, the file will be written inside
448 that directory. Make sure you have write permission.
449
450 The contents of the file is:
451
452 unseq, date-avg, ID, bjd, bvcor, prog_id, exptime, airmass, pmtotal
453
454 @param ID: name of the star, understandable by SIMBAD.
455 @type ID: string
456 @param direc: directory to write the file into (defaults to current working
457 directory)
458 @type direc: string
459 """
460 if direc is None:
461 direc = os.getcwd()
462 data = search(ID)
463 fname = os.path.join(direc,'%s.list'%(ID))
464 ascii.write_array([data['unseq'],data['date-avg'],[ID for i in data['unseq']],
465 data['bjd'],data['bvcor'],data['prog_id'],data['exptime'],
466 data['airmass'],data['pmtotal']],fname,axis0='cols')
467
469 """
470 Make a mask file for RV calculations with the Hermes pipeline.
471
472 See L{ivs.units.linelists} to select appropriate masks. You readily use the
473 columns C{wavelength} and C{depth} as input for this function.
474
475 @param wavelength: wavelength in angstrom
476 @type wavelength: 1D numpy array
477 @param depth: strenght of the line (1-normalised flux minimum)
478 @type depth: 1D numpy array
479 """
480 data = np.array([wavelength,depth]).T
481 hdulist = pf.PrimaryHDU(data)
482 hdulist.writeto(filename)
483
484
486 """
487 Rough calibration of Hermes spectrum.
488
489 Thanks to Roy Ostensen for computing the response function via detailed
490 modeling of Feige 66.
491 """
492
493 instrument = config.get_datafile('catalogs/hermes','hermes_20110319.resp')
494 iwave,iflux = ascii.read2array(instrument).T
495 keep = (iwave.min()<=wave) & (wave<=iwave.max())
496
497 iflux = np.interp(wave[keep],iwave,iflux)
498 return wave[keep],flux[keep]/iflux,iflux
499
500
501
502
504 """
505 Read the velocities.data file with the resuls from hermesVR. The results are
506 returned as a numpy record array. By setting unique to True, the method will
507 only return a unique set of sequence numbers. When return_latest is True it will
508 pick the latest added version, otherwise the first added is picked.
509
510 You can supply a list of unseq numbers, and then only those will be returned,
511 same goes for object. In both cases the output is ordered in the same way as
512 the list of unseq or object that you provided. These lists can be used together
513 with the unique option.
514
515 This method will search for 'velocities.data' in 'hermesDir/hermesAnalyses/'.
516 You can provide another file location by using the keyword filename.
517
518 The field in the recored array are:
519
520 - unseq: unique number
521 - object: object name
522 - hjd: HJD
523 - exptime: exposure time
524 - bvcorr: bary_corr
525 - rvdrift: RV drift (2F frames only)
526 - telloff: telluric offset
527 - tellofferr: error on telluric offset
528 - telloffwidth: width telluric fit
529 - vrad: Vrad(55-74) (bvcorr applied)
530 - vraderr: err on Vrad
531 - nlines: number of lines used
532 - ccfdepth: depth_CCF
533 - ccfdeptherr: error on depth_CCF
534 - ccfsigma: sigma_CCF
535 - ccfsigmaerr: err on sigma_CCF
536 - sn: signal to noise in spectrum (orders 55-74)
537 - gaussoc: O-C of Gaussian fit
538 - wvlfile: wavelength calibration file
539 - maskfile: mask file
540 - fffile: Flatfield file
541
542 @param unique: True is only return unique sequence numbers
543 @type unique: bool
544 @param return_latest: True to pick the last added file
545 @type return_latest: bool
546 @param unseq: list of sequence numbers to return
547 @type unseq: list
548 @param object: list of object names to return
549 @type object: list
550 @param filename: Path to the velocities.data file
551 @type filename: str
552
553 @return: array with the requested result
554 @rtype: record array
555 """
556
557
558 velocity_file = kwargs.pop('filename', hermesDir + 'hermesAnalyses/velocities.data')
559 if not os.path.isfile(velocity_file):
560 logger.warning('velocities.data file not found!')
561 return []
562 rawdata = ascii.read2list(velocity_file)[0]
563
564
565 for i, line in enumerate(rawdata):
566
567 if len(line) > 21:
568 line[20] = " ".join(line[20:])
569 rawdata[i] = tuple(line[0:21])
570
571 dtype = [('unseq','int'),('object','a20'),('hjd','f8'),('exptime','f8'),('bvcor','f8'),
572 ('rvdrift','f8'),('telloff','f8'),('tellofferr','f8'),('telloffwidth','f8'),
573 ('vrad','f8'),('vraderr','f8'),('nlines','int'),('ccfdepth','f8'),
574 ('ccfdeptherr','f8'), ('ccfsigma','f8'),('ccfsigmaerr','f8'),('sn','f8'),
575 ('gaussoc','f8'),('wvlfile','a80'),('maskfile','a80'),('fffile','a80')]
576 data = np.rec.array(rawdata, dtype=dtype)
577
578
579 if unique:
580 unique_sequence = set(data['unseq'])
581 order = -1 if return_latest else 0
582 select = []
583 for seq in unique_sequence:
584 select.append(np.where(data['unseq']==seq)[0][order])
585 data = data[(select,)]
586
587 if unseq != None:
588
589 select = np.ravel(np.hstack([np.where(data['unseq']==sq) for sq in unseq]))
590 data = data[(select,)]
591
592
593 if object != None:
594 select = np.ravel(np.hstack([np.where(data['object']==ob) for ob in object]))
595 data = data[select]
596
597 return data
598
600 """
601 Read the AllCCF.fits files written by hermesVR with the individual cross
602 correlation functions of each order. The cross-correlation functions are stored
603 in a HermesCCF object from which they can be requested.
604
605 You can read *_AllCCF.fits file by providing only the unique sequence number as
606 an integer. In this case the file will be searched in the hermesDir/hermesAnalyses/
607 folder. You can also provide the complete filename instead as a string fx.
608
609 >>> read_hermesVR_AllCCF(457640)
610 >>> read_hermesVR_AllCCF('/home/jorisv/hermesAnalysis/00457640_AllCCF.fits')
611
612 @param unseq: unique sequence number or filename to read
613 @type unseq: int or str
614
615 @return: The cross correlation functions
616 @rtype: HermesCCF object
617 """
618
619 if type(unseq) == int:
620 filename = hermesDir + 'hermesAnalyses/*{:.0f}_AllCCF.fits'.format(unseq)
621 files = glob.glob( filename )
622 if len(files) > 0:
623 filename = files[0]
624 else:
625 filename = unseq
626
627 if not os.path.isfile(filename):
628 logger.warning("Filename ''{:}'' is NOT valid!, returning empty object".format(filename))
629 return HermesCCF(filename=None)
630
631 return HermesCCF(filename=filename)
632
634 """
635 Update the hermesConfig.xml file with the specified keywords. This method will
636 return the original content of the file as a dictionary. You can use that dict
637 as an argument to this function if you want to restore the hermesConfig file to
638 the original state.
639
640 write_hermesConfig will search for hermesConfig.xml in hermesDir/hermesRun/
641 unless a path is specified in the kwarg filename.
642
643 example use:
644
645 >>> old = write_hermesConfig(CurrentNight="/STER/mercator/hermes/20101214")
646 >>> " Do some stuff "
647 >>> write_hermesConfig(**old)
648
649 @attention: This function does not check the validity of the keywords that you
650 provide, that is your own responsibility.
651
652 @param filename: complete path to the hermesConfig.xml file
653 @type filename: str
654
655 @return: The original content of hermesConfig.xml before modification.
656 @rtype: dict
657 """
658
659 config_file = kwargs.pop('filename', hermesDir + 'hermesRun/hermesConfig.xml')
660
661
662 old_config = _etree_to_dict( ET.parse(config_file).getroot() )['hermes']
663
664 if len(kwargs.keys()) == 0.0:
665 logger.debug('Nothing written to hermesConfig.xml, returning current config.')
666 return old_config
667
668
669 new_config = copy.deepcopy(old_config)
670 new_config.update(kwargs)
671
672 hermes = etree.Element("hermes")
673 for key in new_config.keys():
674 child = etree.SubElement(hermes, key)
675 child.text = new_config[key]
676
677
678 out = etree.ElementTree(element=hermes)
679 out.write(config_file, pretty_print=True)
680
681 logger.debug('Written new hermesConfig to file:\n'+etree.tostring(hermes, pretty_print=True))
682
683 return old_config
684
685
686
687
688
689
690 -def run_hermesVR(filename, mask_file=None, wvl_file=None, cosmic_clipping=True,
691 flatfield=False, vrange='default', vrot=False, version='release',
692 verbose=False, timeout=500, **kwargs):
693 """
694 Wrapper around hermesVR that will run hermesVR on any given filename. Thus not on
695 sequence number. This can be usefull if you first merge together 2 back-to-back
696 spectra and then want to calculate the radial velocity of the merged spectrum.
697
698 The original file will not be modified, it is coppied to the tempDir, and the name
699 will be changed according to the supplied flags. During the process hermesConfig.xml
700 is adapted, but after hermesVR is finished, it is restored to its old state.
701
702 When providing a wavelength scale file, provide the full path, but cut off the
703 '_HRF_TH_ext_wavelengthScale.fits' part of the filename. Thus if you want to use:
704 '/folder/2451577_HRF_TH_ext_wavelengthScale.fits', use:
705
706 >>> run_hermesVR(filename, wvl_file = '/folder/2451577')
707
708 The results of hermesVR are saved to the standard paths, as provided in the
709 hermesConfig file. This method only returns the unseq number of the filename
710 so you can find the corresponding output file, the output of hermesVR as a
711 string, and a unix return signal of running hermesVR.
712 The returncode is 0 if no errors are found. Any negative number means hermesVR
713 failed (codes can be consulted here: U{unix signals
714 <http://people.cs.pitt.edu/~alanjawi/cs449/code/shell/UnixSignals.htm>}).
715 If returncode = -35, then hermesVR completed but didn't write any output. In
716 this case consult the logs or string output to find out what went wrong.
717
718 @param filename: complete path to the input spectrum for hermesVR
719 @type filename: str
720 @param mask_file: complete path to the mask file to use or None
721 @type mask_file: str
722 @param wvl_file: path to wavelengthscale file or None
723 @type wvl_file: str
724 @param cosmic_clipping: Doesn't do anything for now
725 @type cosmic_clipping: bool
726 @param flatfield: Use flatfield devision or not
727 @type flatfield: bool
728 @param vrange: Size of velocity window ('default', 'large' (x2), 'xlarge' (x5))
729 @type vrange: str
730 @param vrot: Fit CCF with rotationaly broadend profile instead of gaussian
731 @type vrot: bool
732 @param version: which version of the pipeline to use: 'release' or 'trunk'
733 @type version: str
734 @param verbose: print output of hermesVR
735 @type verbose: bool
736 @param timeout: Maximum time hermesVR is allowed to run in seconds
737 @type timeout: int
738
739 @return: unseq used, text output of hermesVR, unix return code
740 @rtype: (int, str, int)
741
742 """
743
744 vrDir = tempDir + 'hermesvr/'
745 redDir = tempDir + 'hermesvr/reduced/'
746 delRedDir, delVrDir = False, False
747 if not os.path.isdir(vrDir):
748 delVrDir = True
749 os.mkdir(vrDir)
750 if not os.path.isdir(redDir):
751 delRedDir = True
752 os.mkdir(redDir)
753
754
755 header = pf.getheader(filename)
756 if 'unseq' in kwargs:
757 unseq = kwargs.pop('unseq')
758 logger.debug('Using provided sequence number: {:}'.format(unseq))
759 elif 'unseq' in header:
760 unseq = header['unseq']
761 logger.debug('Using sequence number from fits header: {:}'.format(unseq))
762 else:
763 unseq = 1
764 logger.debug('Using default sequence number: {:}'.format(unseq))
765
766
767 vrFile = redDir+ '{:.0f}_HRF_OBJ_ext.fits'.format(unseq)
768 shutil.copyfile(filename, vrFile )
769 logger.debug('Coppied input file to: {:}'.format(vrFile))
770
771
772 if version == 'release':
773 cmd = 'python /STER/mercator/mercator/Hermes/releases/hermes5/pipeline/run/hermesVR.py'
774 oldConfig = write_hermesConfig(Nights=tempDir, CurrentNight=vrDir, Reduced=tempDir)
775 else:
776 cmd = 'python /STER/mercator/mercator/Hermes/trunk/hermes/pipeline/run/hermesVR.py'
777 oldConfig = write_hermesConfig(Nights=vrDir, CurrentNight=vrDir, Reduced=tempDir)
778
779
780 hermesOutput = oldConfig['AnalysesResults'] + "/{:.0f}_AllCCF.data".format(unseq)
781 if os.path.isfile(hermesOutput):
782 os.remove(hermesOutput)
783 logger.debug('Removing existing hermesRV result: {:}'.format(hermesOutput))
784
785
786 runError = None
787 try:
788 cmd += " -i {:.0f}".format(unseq)
789
790 if mask_file:
791 cmd += " -m {:}".format(mask_file)
792
793 if wvl_file:
794 cmd += " -w {:}".format(wvl_file)
795
796 if flatfield and version=='release' or not flatfield and version=='trunk':
797 cmd += " -f"
798
799 if vrange == 'large':
800 cmd += " -L"
801 elif vrange == 'xlarge':
802 cmd += " -LL"
803
804 if vrot and version == 'trunk':
805 cmd += " -ROT"
806
807
808 logger.info('running hermesVR {:} version, with command:\n\t{:}'.format(version, cmd))
809 out, err, returncode = _subprocess_execute(cmd.split(), timeout)
810
811
812 if returncode == -1:
813 logger.warning('Timeout, hermesVR was terminated after' \
814 +' {:.0f}s'.format(timeout))
815 elif returncode < 0:
816 logger.warning('hermesVR terminated with signal: ' \
817 +' {:.0f}s'.format(abs(returncode)))
818 elif not os.path.isfile(hermesOutput):
819 logger.warning('hermesVR finished but did not succeed, check log for details')
820 returncode = -35
821 logger.debug(out)
822
823 if verbose: print out
824
825 except Exception, e:
826 runError = e
827
828
829 write_hermesConfig(**oldConfig)
830 os.remove(vrFile)
831 if delRedDir: os.rmdir(redDir)
832 if delVrDir: os.rmdir(vrDir)
833
834
835 if runError: raise runError
836
837 return unseq, out, returncode
838
839 -def CCFList(ID,config_dir=None,out_dir=None,mask_file=None,cosmic_clipping=True,flatfield=False):
840 """
841 Calculate radial velocities for Hermes stars.
842
843 @warning: you have to have the DRS installed locally!
844
845 @warning: this function changes your hermesConfig.xml file, you'd better
846 backup it before use
847 """
848
849 make_list_star(ID)
850 cmd = 'python /STER/mercator/mercator/HermesTools/CCFList.py -i %s.list'%(ID)
851 if mask_file is not None:
852 cmd += ' -m %s'%(os.path.abspath(mask_file))
853 if not cosmic_clipping:
854 cmd += ' -nc'
855 if flatfield:
856 cmd += ' -f'
857
858
859 if config_dir is None:
860 config_dir = os.path.expanduser('~/hermesRun')
861 if out_dir is None:
862 out_dir = os.getcwd()
863 out_dir = os.path.abspath(out_dir)
864
865 config_file = r"""<hermes>
866 <Nights>/STER/mercator/hermes</Nights>
867 <CurrentNight>/STER/mercator/hermes/20110501</CurrentNight>
868 <Reduced>/STER/mercator/hermes/20110501/reduced</Reduced>
869 <AnalysesResults>%s</AnalysesResults>
870 <DebugPath>%s</DebugPath>
871 <Raw>raw</Raw>
872 <ConsoleLogSeverity>debug</ConsoleLogSeverity>
873 <LogFileName>%s</LogFileName>
874 <LogFileFormater>%%(asctime)s :: %%(levelname)s ::
875 %%(message)s</LogFileFormater>
876 </hermes>"""%(out_dir,out_dir,os.path.join(out_dir,'hermes.log'))
877
878
879 config = open(os.path.join(config_dir,'hermesConfig.xml'),'w')
880 config.write(config_file)
881 config.close()
882
883 try:
884 retcode = subprocess.call(cmd, shell=True)
885 if retcode < 0:
886 logger.error("Child (%s) was terminated by signal %d"%(cmd,-retcode))
887 else:
888 logger.error("Child (%s) returned %d"%(cmd,retcode))
889 except OSError as e:
890 logger.error("Execution (%s) failed: %s"%(cmd,e))
891
892
893
894
895
896
897
898
899
901 """
902 Object that holds all information stored in the *_AllCCF.fits files written by hermesVR.
903 The original information is stored in the following variables that can be accessed
904 directly:
905
906 - filename: The path to the *_AllCCF.fits file
907 - header: The header of the original spectra that was the input to hermesVR
908 - groups: The summary that hermesVR prints of the 5 preset order groups stored
909 as a rec array
910 - ccfs: All cross-correlation functions stored as a array (index 0-54)
911 - vr: the velocity scale of the cross correlation functions (array)
912
913 Appart from these attributes, two functions allow you to read the ccf by order or
914 a list of orders where the order numbers are the original hermes orders (40-94).
915 """
916
918 """
919 Read the provided file and store all information.
920
921 @param filename: The complete path to the AllCCF.fits file to read.
922 @type filename: string
923 """
924 self.header = None
925 self.vr = None
926 self.ccfs = None
927 self.groups = None
928
929 if filename != None:
930 self._read_file(filename)
931 else:
932 filename = None
933
934
935
936 - def ccf(self, order):
937 """
938 Return the ccf function belonging to the requested order. The order number
939 is the real order number (40 <= order <= 94).
940
941 @param order: The hermes order of which you want the ccf
942 @type order: int
943
944 @return: The cross correlation function of the requested order
945 @rtype: array
946
947 @raise IndexError: When the order is out of bounds.
948 """
949 order = 94 - order
950
951 if order < 0 or order > 54:
952 raise IndexError('Order {:.0f} is out of bounds (40 - 94)'.format(94 - order))
953
954 return self.ccfs[order]
955
957 """
958 Sum the ccf function belonging to the requested order numbers. The order
959 numbers can be given as a list of integers or in a string representation.
960 When using the string representation, you can use the ':' sign to give a
961 range of order which will include the starting and ending order. Fx the
962 following two commands will return the same summed ccf:
963
964 >>> combine_ccf([54,55,56,57,75,76])
965 >>> combine_ccf('54:57,75,76')
966
967 @param orders: The hermes orders of which you want the summed ccf
968 @type orders: list or string
969
970 @return: The summed cross correlation function of the requested orders
971 @rtype: array
972
973 @raise IndexError: When an order is out of bounds.
974 """
975
976 if type(orders) == str:
977 orders = orders.split(',')
978 o_list = []
979 for o in orders:
980 if ':' in o:
981 o = o.split(':')
982 o_list.extend( range( int(o[0]), int(o[1])+1 ) )
983 else:
984 o_list.append( int(o) )
985 logger.debug("converted order string: ''{:}'' to list: {:}".format(orders, o_list))
986 orders = o_list
987
988 if np.any([(o<40) | (o>94) for o in orders]):
989 raise IndexError('Some/All orders {:} are out of bounds (40 - 94)'.format(orders))
990
991 ccf = np.zeros_like(self.vr)
992 for o in orders:
993 ccf += self.ccf(o)
994
995 return ccf
996
997
998
999
1000
1002 "Read the *_AllCCF.fits file "
1003 self.filename = filename
1004
1005 hdu = pf.open(filename)
1006
1007
1008 self.header = hdu[0].header
1009
1010
1011 self.ccfs = hdu[0].data.T
1012 self.vr = hdu[2].data.field('VR')
1013
1014
1015 self.groups = np.array(hdu[1].data, dtype = hdu[1].data.dtype)
1016
1017 logger.debug('Read ccf file: {:}'.format(filename))
1018
1019 return
1020
1022 "String representation of this object "
1023 if self.header == None:
1024 return "< empty Hermes CCF object >"
1025 else:
1026 obj, sq = self.header['object'], self.header['unseq']
1027 return "< CCF of {:} with unseq {:.0f} >".format(obj, sq)
1028
1029
1030
1031
1032
1034 """
1035 Summarize all Hermes data in a file for easy data retrieval.
1036
1037 The file is located in one of date data directories (see C{config.py}), in
1038 subdirectories C{catalogs/hermes/HermesFullDataOverview.tsv}. If it doesn't
1039 exist, it will be created. It contains the following columns, which are
1040 extracted from the Hermes FITS headers (except C{filename}:
1041
1042 1. UNSEQ
1043 2. PROG_ID
1044 3. OBSMODE
1045 4. BVCOR
1046 5. OBSERVER
1047 6. OBJECT
1048 7. RA
1049 8. DEC
1050 9. BJD
1051 10. EXPTIME
1052 11. PMTOTAL
1053 12. DATE-AVG
1054 13. OBJECT
1055 14. airmass
1056 15. filename
1057
1058 This file can most easily be read with the L{ivs.inout.ascii} module and the
1059 command:
1060
1061 >>> hermes_file = config.get_datafile(os.path.join('catalogs','hermes'),'HermesFullDataOverview.tsv')
1062 >>> data = ascii.read2recarray(hermes_file,splitchar='\\t')
1063
1064 """
1065 logger.info('Collecting files...')
1066
1067 dirs = sorted(glob.glob(os.path.join(config.ivs_dirs['hermes'],'20??????')))
1068 dirs = [idir for idir in dirs if os.path.isdir(idir)]
1069 obj_files = []
1070
1071 for idir in dirs:
1072 obj_files += sorted(glob.glob(os.path.join(idir,'raw','*.fits')))
1073
1074
1075
1076
1077
1078
1079 try:
1080 overview_file = config.get_datafile('catalogs/hermes','HermesFullDataOverview.tsv')
1081
1082 overview_data = ascii.read2recarray(overview_file,splitchar='\t')
1083 outfile = open(overview_file,'a')
1084 logger.info('Found %d FITS files: appending to overview file %s'%(len(obj_files),overview_file))
1085
1086 except IOError:
1087 overview_file = 'HermesFullDataOverview.tsv'
1088 outfile = open(overview_file,'w')
1089 outfile.write('#unseq prog_id obsmode bvcor observer object ra dec bjd exptime pmtotal date-avg airmass filename\n')
1090 outfile.write('#i i a20 >f8 a50 a50 >f8 >f8 >f8 >f8 >f8 a30 >f8 a200\n')
1091 overview_data = {'filename':[]}
1092 logger.info('Found %d FITS files: starting new overview file %s'%(len(obj_files),overview_file))
1093
1094
1095 existing_files = np.sort(overview_data['filename'])
1096 for i,obj_file in enumerate(obj_files):
1097 sys.stdout.write(chr(27)+'[s')
1098 sys.stdout.write(chr(27)+'[2K')
1099 sys.stdout.write('Scanning %5d / %5d FITS files'%(i+1,len(obj_files)))
1100 sys.stdout.flush()
1101
1102
1103 index = existing_files.searchsorted(obj_file)
1104 if index<len(existing_files) and existing_files[index]==obj_file:
1105 sys.stdout.write(chr(27)+'[u')
1106 continue
1107
1108
1109
1110
1111 contents = dict(unseq=-1,prog_id=-1,obsmode='nan',bvcor=np.nan,observer='nan',
1112 object='nan',ra=np.nan,dec=np.nan,
1113 bjd=np.nan,exptime=np.nan,pmtotal=np.nan,airmass=np.nan,
1114 filename=os.path.realpath(obj_file))
1115 contents['date-avg'] = 'nan'
1116 header = pf.getheader(obj_file)
1117 for key in contents:
1118 if key in header and key in ['unseq','prog_id']:
1119 try: contents[key] = int(header[key])
1120 except: pass
1121 elif key in header and key in ['obsmode','observer','object','date-avg']:
1122 contents[key] = str(header[key])
1123 elif key in header and key in ['ra','dec','exptime','pmtotal','bjd','bvcor']:
1124 contents[key] = float(header[key])
1125 elif key=='airmass' and 'telalt' in header:
1126 if float(header['telalt'])<90:
1127 try:
1128 contents[key] = airmass.airmass(90-float(header['telalt']))
1129 except ValueError:
1130 pass
1131
1132 outfile.write('%(unseq)d\t%(prog_id)d\t%(obsmode)s\t%(bvcor)f\t%(observer)s\t%(object)s\t%(ra)f\t%(dec)f\t%(bjd)f\t%(exptime)f\t%(pmtotal)f\t%(date-avg)s\t%(airmass)f\t%(filename)s\n'%contents)
1133 outfile.flush()
1134 sys.stdout.write(chr(27)+'[u')
1135 outfile.close()
1136 return overview_file
1137
1139 """
1140 Derive the location of a reduced file from the raw file.
1141 """
1142 redfile = rawfile.replace('raw','reduced')
1143 redfiledir,redfilebase = os.path.dirname(redfile),os.path.basename(redfile)
1144 base,ext = os.path.splitext(redfilebase)
1145 if data_type.lower()=='cosmicsremoved_log':
1146 redfile = os.path.join(redfiledir,'_'.join([base,'ext','CosmicsRemoved','log','merged','c']))+'.fits'
1147 elif data_type.lower()=='cosmicsremoved_wavelength':
1148 redfile = os.path.join(redfiledir,'_'.join([base,'ext','CosmicsRemoved','wavelength','merged','c']))+'.fits'
1149 elif data_type.lower()=='log':
1150 redfile = os.path.join(redfiledir,'_'.join([base,'ext','log','merged']))+'.fits'
1151 elif data_type.lower()=='wavelength':
1152 redfile = os.path.join(redfiledir,'_'.join([base,'ext','wavelength','merged']))+'.fits'
1153 else:
1154 redfile = rawfile
1155
1156 if not os.path.isfile(redfile):
1157 return 'naf'
1158 return redfile
1159
1161 """
1162 Convert the time stamp from a HERMES FITS 'date-avg' to Julian Date.
1163
1164 @param timestamp: string from 'date-avg'
1165 @type timestamp: string
1166 @return: julian date
1167 @rtype: float
1168 """
1169 date, hour = timestamp.split("T")
1170 year, month, day = date.split("-")
1171 hour, minute, second = hour.split(":")
1172 year = float(year)
1173 month = float(month)
1174 day = float(day)
1175 hour = float(hour)
1176 minute = float(minute)
1177 second = float(second)
1178 return conversions.convert("CD","JD",(year, month, day, hour, minute, second))
1179
1181 """
1182 Convert the time stamp from a HERMES FITS 'date-avg' to a datetime object.
1183
1184 @param timestamp: string from 'date-avg'
1185 @type timestamp: string
1186 @return: datetime object
1187 @rtype: datetime
1188 """
1189 if timestamp=='nan':
1190 timestamp = '1000-01-01T00:00:00'
1191
1192 timestamp = [int(i) for i in ' '.join(' '.join(' '.join(' '.join(timestamp.split('-')).split('T')).split(':')).split('.')).split()]
1193
1194 if len(timestamp)==3:
1195 timestamp += [12,0,0]
1196 return datetime.datetime(*timestamp)
1197
1199 """
1200 Convert a xml tree to a dictionary.
1201
1202 @param t: the xml tree
1203 @type t: ElementTree root
1204
1205 @return: the xml tree as a dictionary
1206 @rtype: dict
1207 """
1208 d = {t.tag: {} if t.attrib else None}
1209 children = list(t)
1210 if children:
1211 dd = defaultdict(list)
1212 for dc in map(_etree_to_dict, children):
1213 for k, v in dc.iteritems():
1214 dd[k].append(v)
1215 d = {t.tag: {k:v[0] if len(v) == 1 else v for k, v in dd.iteritems()}}
1216 if t.attrib:
1217 d[t.tag].update(('@' + k, v) for k, v in t.attrib.iteritems())
1218 if t.text:
1219 text = t.text.strip()
1220 if children or t.attrib:
1221 if text:
1222 d[t.tag]['#text'] = text
1223 else:
1224 d[t.tag] = text
1225 return d
1226
1228 """executing the command with a watchdog"""
1229
1230
1231 c = subprocess.Popen(command, stdout=subprocess.PIPE)
1232
1233
1234 t = 0
1235 while t < time_out and c.poll() is None:
1236 time.sleep(1)
1237 t += 1
1238
1239
1240 if c.poll() is None:
1241
1242 c.terminate()
1243
1244 out, err = c.communicate()
1245 returncode = -1
1246
1247 else:
1248
1249 out, err = c.communicate()
1250 returncode = c.poll()
1251
1252 return out, err, returncode
1253
1254
1255
1256 if __name__=="__main__":
1257
1258 from ivs.aux import loggers
1259 logger = loggers.get_basic_logger(clevel='debug')
1260
1261 filename = '/home/jorisv/sdB/Uli/hermesvr/reduced/2451576_HRF_OBJ_ext_CosmicsRemoved.fits'
1262 wvl_file = '/home/jorisv/sdB/Uli/hermesvr/reduced/2451577'
1263
1264 unseq,output, error = run_hermesVR(filename, wvl_file=wvl_file, verbose=False, version='trunk', timeout=500)
1265
1266 data = read_hermesVR_velocities(unseq=[unseq])
1267
1268 print data['vrad'], ' +- ', data['vraderr']
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320